無料でベイジアン構造方程式モデル

年明けに吉備国際大学大学院(修士課程)で構造方程式モデルの講義&演習するため、無料でベイジアン構造方程式モデル(Bayesian SEM、BSEM)できる方法を探しました。
ありました。

RでBSEMでできるblavaanパッケージだす。

ざっと試してみました。

普段、ぼくはMplusでBSEMしています。

でも、MplusのBSEMの情報量規準はDICとBICで、WAICがでないし、有料なんですよね。

その点、blavaanはWAICを出力してくれますし、無料です。

ただ現状、カテゴリカルデータには対応してないっぽいです(たぶん。もしかしたらできるのかも)。

年明けのSEMの講義&演習までに、わかりやすいパワーポイント資料を作成しますから、このブログではとりあえずの備忘録です。






blavaanはlavaanのベイズ版で、jagsをこき使うことでCFASEMLCMができます。

なので、blavaanの基本的なモデル式はlavaanと一緒です。

以下、blavaanとlavaanでCFAやりました。

データはlavaanに収録されているものです。

まずはblavaanです。

library(blavaan)
data(HolzingerSwineford1939)
HS.model <- '
visual =~ x1 + x2 + x3
textual =~ x4 + x5 + x6
speed =~ x7 + x8 + x9
' 
fit <- blavaan(HS.model, data=HolzingerSwineford1939, auto.var=TRUE, auto.fix.first=TRUE, auto.cov.lv.x=TRUE, jagcontrol=list(method="rjparallel"))
summary(fit,standardized=TRUE)
fitMeasures(fit)
以下、主な結果の一部を示します。

std.allは標準化係数を表しています(太字)。

収束基準は青字です。
> summary(fit,standardized=TRUE)Latent Variables:
 Post.Mean Post.SD HPD.025 HPD.975 Std.lv Std.all PSRF Prior
 visual =~ 
 x1 1.000 4.999 0.986 
 x2 1.221 0.018 1.187 1.256 6.104 0.984 1.000 dnorm(0,1e-2)
 x3 0.463 0.012 0.439 0.488 2.315 0.919 1.000 dnorm(0,1e-2)
 textual =~ 
 x4 1.000 3.205 0.979 
 x5 1.402 0.020 1.364 1.441 4.495 0.992 1.000 dnorm(0,1e-2)
 x6 0.730 0.016 0.699 0.761 2.339 0.956 1.000 dnorm(0,1e-2)
 speed =~ 
 x7 1.000 4.225 0.976 
 x8 1.319 0.020 1.279 1.357 5.571 0.991 1.000 dnorm(0,1e-2)
 x9 1.284 0.019 1.248 1.323 5.426 0.992 1.000 dnorm(0,1e-2)
次にblavaanの情報量規準です。

DICやBICに加えて、WAIC、LOOがでていますね(うれぴー)。

モデルの拡張性はMplusの完勝ですが、情報量規準はblavaanのほうがええですね。

びみょーに痒いところに手が届かない感じがなんとも^^;。

それはともかく、モデル選択は複数のモデルのWAICなどの情報量規準を比較し、その値が小さいものを選んでいきます。
> fitMeasures(fit)
 npar logl ppp bic dic p_dic waic p_waic looic 
 21.000 -4398.277 0.000 8916.333 8837.552 20.499 8837.975 20.632 8837.974 
 p_loo margloglik 
 20.631 -4481.119
次にlavaanです。
library(lavaan)
data(HolzingerSwineford1939)
HS.model <- '
visual =~ x1 + x2 + x3
textual =~ x4 + x5 + x6
speed =~ x7 + x8 + x9
' 
fit <- sem(HS.model,data=HolzingerSwineford1939, estimator = "MLR")
summary(fit, fit.measures=TRUE,standardized=TRUE)
以下、主な結果の一部を抜き出します。

std.allは標準化係数を表しています(太字)。

うーん、、、同一モデル同一データのはずなんですけど、BSEMの結果に比べると標準化係数の値がだいぶ小さいです(やり方が間違ってるんだろうか)。

今度、MplusのBSEMで結果がどうなるか調べみたいと思います。
Latent Variables:
 Estimate Std.Err Z-value P(>|z|) Std.lv Std.all
 visual =~ 
 x1 1.000 0.900                   0.772
 x2 0.554 0.132 4.191 0.000 0.498 0.424
 x3 0.729 0.141 5.170 0.000 0.656 0.581
 textual =~ 
 x4 1.000 0.990                    0.852
 x5 1.113 0.066 16.946 0.000 1.102 0.855
 x6 0.926 0.061 15.089 0.000 0.917 0.838
 speed =~ 
 x7 1.000 0.619                   0.570 x8 1.180 0.130 9.046 0.000 0.731 0.723
 x9 1.082 0.266 4.060 0.000 0.670 0.665


今日は時間切れなので、とりあえずここまで。

年明けのSEM(BSEM)の講義&演習に向けて、ぼちぼち準備していきまする。