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
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