MplusでMCMCした結果を使ったSEMでの情報仮説の評価

ちょっと気になったので，近年一部界隈で話題になっているベイズ情報仮説の評価をやってみる。
普通に計算すると大変なベイズファクターも，情報仮説の元では簡単に計算できる。

元ネタは以下の資料
http://www.tandfonline.com/doi/abs/10.1080/17405629.2012.732719
https://www.google.co.jp/url?sa=t&rct=j&q=&esrc=s&source=web&cd=1&cad=rja&uact=8&ved=0ahUKEwjZ9-6K-PXOAhVGpJQKHf9-B2kQFgghMAA&url=https%3A%2F%2Frensvandeschoot.files.wordpress.com%2F2012%2F07%2F02-bayes-informative-hypothesis.ppt&usg=AFQjCNHLedgBUV_FFYDOK-6AQVfL-55v3A&sig2=0kAg9LL2V3VDKPBcIF3yfg&bvm=bv.131783435,d.dGo
http://www3.psy.senshu-u.ac.jp/~ken/JJPS2014.pdf
https://www.google.co.jp/url?sa=t&rct=j&q=&esrc=s&source=web&cd=4&cad=rja&uact=8&ved=0ahUKEwjP6YGA-_XOAhUHsJQKHVFFASQQFggyMAM&url=http%3A%2F%2Fir.acc.senshu-u.ac.jp%2Findex.php%3Faction%3Dpages_view_main%26active_action%3Drepository_action_common_download%26item_id%3D9914%26item_no%3D1%26attribute_id%3D15%26file_no%3D1%26page_id%3D13%26block_id%3D52&usg=AFQjCNFFbxBCCZtGqypBtdGwYOs0BKsRmw&sig2=BKGhm-TClAlSrhKU4aKDSw&bvm=bv.131783435,d.dGo
などなど。
あまり，この手の話題に詳しくはないが，MplusでMCMCをする練習にもなるので，よいかと思う。

# 使用するパッケージ
#　install.packages("MplusAutomation")
library(MplusAutomation)
# Mplusの出力から，R上でトレースプロット等を描くためのもの。
source("http://bioconductor.org/biocLite.R")
biocLite("rhdf5")

folder <- "/Users/hoge/Dropbox/mplus/infohypo"

# 
# Mplusのinpファイルの作成
# 
cat("DATA: FILE = ex3.1.dat;
! User guideの例を利用
VARIABLE: 
NAMES ARE y1-y6 x1-x4;
USEVARIABLES ARE y1 x1-x3;
MISSING=.;

! ここの指定が大事
analysis:
! 推定方法をベイズに
estimator = bayes;
! チェイン数と連動する
processors =7;
chains = 7;
! イテレーション数。
! 何も指定しないと，半分がバーンインで捨てられる。
fbiter = 20000;
thin = 2;
! 初期値をMLで決める
stvalues = ML;
! アルゴリズムメトロポリスヘイスティングズも指定可能
algorithm = gibbs(px3);
!algorithm = MH; 

model:
y1 on x1-x3*(beta1-beta3);

! 事前分布の指定。
! なくても勝手に決めてくれる。
model priors:
beta1-beta3 ~ N(0, 1);

! 
output:
tech1 tech8 stdy;

! トレースプロットなど
plot:
type = plot2;

savedata:
! MCMCサンプルを保存
bparameters are bprm.dat;
!results are params.dat;
!ESTIMATES ARE estimate.dat;
  ",file=paste0(folder,"lm.inp"))

# モデルの実行
runModels(paste0(folder))


setwd(folder)
getwd()

# MplusAutomationの利用
btest <- getSavedata_Bparams(outfile = "lm.out", discardBurnin = TRUE)

# 変数名が独特なので，colnamesで確認
colnames(btest[[1]] )

# 仮説の表記の仕方。
# beta1 >beta2 かつ beta2 > beta3
testBParamCompoundConstraint(btest,"(Parameter.2_Y1.ON.X1..equality.label.>Parameter.3_Y1.ON.X2..equality.label.)&
  (Parameter.3_Y1.ON.X2..equality.label.>Parameter.4_Y1.ON.X3..equality.label.)")

# 制約がmetなものの割合を見ればよい？
0.2419/(1/6)

?testBParamCompoundConstraint
# beta1 < beta2 かつ　beta2 < beta3
testBParamCompoundConstraint(btest,"(Parameter.2_Y1.ON.X1..equality.label.<Parameter.3_Y1.ON.X2..equality.label.)&
  (Parameter.3_Y1.ON.X2..equality.label.<Parameter.4_Y1.ON.X3..equality.label.)")
# f_i に当たるものが算出出来る
# c_i に当たるものは，3!のなかの組合せのうち，当該の仮説の個数を計算すればよい。
# ここでは1通り。仮説が，beta1 > beta2, beta3 なら，2通りある。
0.1015571 /(1/6)


# 仮説が正しい場合と正しくない場合の比
(0.1015571 /(1/6)) /((1-0.1015571) /(1-1/6))


0.2419/(1/6) /(0.1015571 /(1/6))

# その他関数。
extractModelParameters()
extractModelSummaries()
getSavedata_Fileinfo("lm.out")

# トレースプロットをかいたりできる。
model <- readModels("lm.out")
mplus.traceplot(model, rows = 2, cols = 2, parameters_only = TRUE)
library(coda)
plot(model$bparameters$valid_draw)