Dat1 <- read.table("y.csv", header=TRUE, sep=",",na.strings="NA", dec=".", strip.white=TRUE)
#y.csvのファイルレイアウト（観測モデルで使用する）
#SID：パネルID
#y:選択結果
#x0：定数項
#x1：勤務先
#x2：仕事内容
#x3：上司との人間関係
#x4：同僚との人間関係
#x5：取引先との人間関係
#x6：給料
#x7：通勤時間
#x8：会社の将来性
#x9：労働時間
#x10：自身への会社の評価
#
IndAttr <- read.table("z.csv", header=FALSE,sep=",", na.strings="NA", dec=".", strip.white=TRUE)
#z.csvのファイルレイアウト（階層モデルで使用する）
#1列目：性別
#2列目：年齢
#3列目：独身
#4列目：夫婦
#5列目：製造業
#6列目：運輸・通信
#7列目：金融・保険
#8列目：流通サ－ビス
#9列目：情報サ－ビス
#10列目：コンサルティング
#11列目：医療・医薬関係
#12列目：研究
#13列目：企画関連
#14列目：営業
#15列目：部長以上
#16列目：課長・係長
#
#MCMCの設定
library(bayesm) #使用するライブラリのロード
#
R=20000　#MCMCの繰り返し回数
s<-15001　#事後統計量算定の初期時点
t<-20000　#事後統計量算定の最終時点
keep=1　#MCMCで発生させたサンプルを何個おきにストアするかを決めるための設定.「1」であれば発生したものを全てストアする.例えば「5」ならば5個おきにストアすることを意味する.
#
#
#個人数
reg=levels(factor(Dat1$SID))
nreg=length(reg) #分析対象個人数
#選択肢数
p=2
#選択肢共通変数数（観測モデルに含まれる説明変数数）
na=11
nvar<-11
#個人属性数（階層モデルの説明変数数）
nz=16
#
#個人ごとデータの作成
lgtdata=NULL
for (j in 1:nreg){
y=Dat1$y[Dat1$SID==reg[j]]
y=as.matrix(y)
nobs=length(y)
X<-cbind(Dat1$x0[Dat1$SID==reg[j]],
Dat1$x1[Dat1$SID==reg[j]],
Dat1$x2[Dat1$SID==reg[j]],
Dat1$x3[Dat1$SID==reg[j]],
Dat1$x4[Dat1$SID==reg[j]],
Dat1$x5[Dat1$SID==reg[j]],
Dat1$x6[Dat1$SID==reg[j]],
Dat1$x7[Dat1$SID==reg[j]],
Dat1$x8[Dat1$SID==reg[j]],
Dat1$x9[Dat1$SID==reg[j]],
Dat1$x10[Dat1$SID==reg[j]]
)
lgtdata[[j]]=list(y=y,X=X)
}
#MCMCの処理
Z=NULL
Z=as.matrix(IndAttr)
Data1=list(lgtdata=lgtdata,Z=Z)
Mcmc1=list(R=R,keep=1)
set.seed(66)
out=rhierBinLogit (Data=Data1,Mcmc=Mcmc1)　#bayesmライブラリに含まれる階層ベイズ2項ロジットを推定するための関数
#
#観測モデルのパラメータの事後平均の算定
beta.mean<-matrix(0,nrow=nreg,ncol=nvar)
for(i in 1:nreg){
beta.mean[i,1]<-mean(out$betadraw[i,1,s:t])
beta.mean[i,2]<-mean(out$betadraw[i,2,s:t])
beta.mean[i,3]<-mean(out$betadraw[i,3,s:t])
beta.mean[i,4]<-mean(out$betadraw[i,4,s:t])
beta.mean[i,5]<-mean(out$betadraw[i,5,s:t])
beta.mean[i,6]<-mean(out$betadraw[i,6,s:t])
beta.mean[i,7]<-mean(out$betadraw[i,7,s:t])
beta.mean[i,8]<-mean(out$betadraw[i,8,s:t])
beta.mean[i,9]<-mean(out$betadraw[i,9,s:t])
beta.mean[i,10]<-mean(out$betadraw[i,10,s:t])
beta.mean[i,11]<-mean(out$betadraw[i,11,s:t])
}
#観測モデルのパラメータの事後標準偏差の算定
beta.sd<-matrix(0,nrow=nreg,ncol=nvar)  
for(i in 1:nreg){
beta.sd[i,1]<-sd(out$betadraw[i,1,s:t])
beta.sd[i,2]<-sd(out$betadraw[i,2,s:t])
beta.sd[i,3]<-sd(out$betadraw[i,3,s:t])
beta.sd[i,4]<-sd(out$betadraw[i,4,s:t])
beta.sd[i,5]<-sd(out$betadraw[i,5,s:t])
beta.sd[i,6]<-sd(out$betadraw[i,6,s:t])
beta.sd[i,7]<-sd(out$betadraw[i,7,s:t])
beta.sd[i,8]<-sd(out$betadraw[i,8,s:t])
beta.sd[i,9]<-sd(out$betadraw[i,9,s:t])
beta.sd[i,10]<-sd(out$betadraw[i,10,s:t])
beta.sd[i,11]<-sd(out$betadraw[i,11,s:t])
}
#観測モデルのパラメータの事後t-値の算定
beta.t<-matrix(0,nrow=nreg,ncol=nvar)
for(i in 1:nreg){
beta.t[i,1]<-beta.mean[i,1]/beta.sd[i,1]
beta.t[i,2]<-beta.mean[i,2]/beta.sd[i,2]
beta.t[i,3]<-beta.mean[i,3]/beta.sd[i,3]
beta.t[i,4]<-beta.mean[i,4]/beta.sd[i,4]
beta.t[i,5]<-beta.mean[i,5]/beta.sd[i,5]
beta.t[i,6]<-beta.mean[i,6]/beta.sd[i,6]
beta.t[i,7]<-beta.mean[i,7]/beta.sd[i,7]
beta.t[i,8]<-beta.mean[i,8]/beta.sd[i,8]
beta.t[i,9]<-beta.mean[i,9]/beta.sd[i,9]
beta.t[i,10]<-beta.mean[i,10]/beta.sd[i,10]
beta.t[i,11]<-beta.mean[i,11]/beta.sd[i,11]
}
#観測モデルのパラメータの事後統計量の書き出し
write.table(beta.mean, file="beta_mean.txt") #個人ごとパラメータの事後平均のファイル（nreg×nvar）
write.table(beta.sd, file="beta_SD.txt") #個人ごとパラメータの事後標準偏差のファイル（nreg×nvar）
write.table(beta.t, file="beta_t.txt")#個人ごとパラメータの事後t-値のファイル（nreg×nvar）
#
kflg=nz*nvar
kflg2=s:t
#
#階層モデルのパラメータの事後平均の算定
Delta.mean<-matrix(0,nrow=1,ncol=nvar*nz)
#
for (i in 1:kflg) {
Delta.mean[1,i]<-mean(out$Deltadraw[kflg2,i])
}
Delta.mat<-matrix(0,nrow=nvar,ncol=nz)
kk<-1
for (i in 1:nvar) {
for (j in 1:nz) {
Delta.mat[i,j]<-Delta.mean[1,kk]
kk<-kk+1
}
}
write.table(Delta.mat, file="Delta_mean.txt")#階層パラメータの事後平均のファイル（nvar×nz）
#
#階層モデルのパラメータの事後標準偏差の算定
Delta.SD<-matrix(0,nrow=1,ncol=nvar*nz)
for (i in 1:kflg) {
Delta.SD[1,i]<-sd(out$Deltadraw[kflg2,i])
}
Delta.mat<-matrix(0,nrow=nvar,ncol=nz)
kk<-1
for (i in 1:nvar) {
for (j in 1:nz) {
Delta.mat[i,j]<-Delta.SD[1,kk]
kk<-kk+1
}
}
write.table(Delta.mat, file="Delta_SD.txt")#階層パラメータの事後標準偏差のファイル（nvar×nz）
#
#階層モデルのパラメータの事後t-値の算定
Delta.t<-matrix(0,nrow=1,ncol=nvar*nz)
for (i in 1:kflg) {
Delta.t[1,i]<-Delta.mean[1,i]/Delta.SD[1,i]
}
Delta.mat<-matrix(0,nrow=nvar,ncol=nz)
kk<-1
for (i in 1:nvar) {
for (j in 1:nz) {
Delta.mat[i,j]<-Delta.t[1,kk]
kk<-kk+1
}
}
write.table(Delta.mat, file="Delta_t.txt")#階層パラメータの事後t‐値のファイル（nvar×nz）
#