bayesmを使った階層ベイズモデル（線形回帰）

通常の重回帰分析を階層ベイズモデルに拡張するには、Rのbayesmを使うのが簡単です。
bayesmは、いろいろな階層モデルが用意されていますが、重回帰分析のモデルは
rhierLinearModel
となります。

MCMCの計算自体はとても簡単です。

階層モデルのパラメータDelta.meanと観測モデルのパラメータbeta.meanの行列計算が間違えるポイントかもしれません。
パラメータが2次元の行列なのと、さらに、MCMCごとのデータがあるので、3次元の構造となっています。
どこの軸を集計すると、良いのか、慣れるまで時間がかかるかもしれません。

[観測モデル]
y=bo+b1*x1+b2*x2+b3*x3+ε

[階層モデル]
b0=a1*z1+a2*z2+ε
b1=a1*z1+a2*z2+ε
b2=a1*z1+a2*z2+ε

簡単な使い方の例

library(bayesm)
R=20000

# 個人数
nreg <- length(member_id) # num of id

# Zは、階層モデル用のデータ
nz <- ncol(Z) # カラム数

# 観測モデル用のデータ
regdata=NULL
for (reg in 1:nreg) {
regdata[[reg]] <- list(y=y, X=X)
}

# MCMC
Data1 = list(regdata=regdata, Z=Z)
Mcmc1 = list(R=R, keep=1)
Mcmc_out = rhierLinearModel(Data=Data1, Mcmc=Mcmc1)

# Burn-In の設定
mcmc_start <- R/2
mcmc_end <- R

# パラメータの事後平均のデビアンス
for(j in 1:nvar_para){
for(i in 1:nreg){beta.mean[i,j] <- mean(Mcmc_out$betadraw[i, j, (mcmc_start:mcmc_end)])}
}

## 階層モデルのパラメータ
Delta.mean <- matrix(0, nz, nvar_para)

for(j in 1:nvar_para) {
for(i in 1:nz) {
Delta.mean[i,j] <- mean(Mcmc_out$Deltadraw[(mcmc_start:mcmc_end), i+(j-1)*nz])
}
}
Delta.mean