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)
Lavaanによるパラメータ推定結果を見やすくグラフにする。(Gelman's secret weaponのパクリ)
尊敬するHoxo_mさんらが翻訳された”みんなのR”にあるGelman's secret weaponを真似して自分の武器を作成してみようかと思う。
d.hatena.ne.jp
https://www.amazon.co.jp/みんなのR-データ分析と統計解析の新しい教科書-Jared-P-Lander-ebook/dp/B012Z86Q6Q
やっていることは簡単。lavaanパッケージでのパラメータ推定結果を流用して,ggplotのgeom_pointrangeを使って点推定値と95%信頼区間を描画しているだけ。
おまけで,5%水準で有意かどうか,固定パラメータかで,色の塗り分けなんかをやってみている。
これで,推定値とそのブレ,有意か否かがさくっと把握できる。
分散・共分散を除きたい場合なんかは投入するデータセットを変更すればよろしい。
色の好みなどは適宜変更すること。
# 使用するパッケージ library(ggplot2) library("lavaan") # help(package="lavaan") # sem関数の例題を利用 ## The industrialization and Political Democracy Example ## Bollen (1989), page 332 model <- ' # latent variable definitions ind60 =~ x1 + x2 + x3 dem60 =~ y1 + a*y2 + b*y3 + c*y4 dem65 =~ y5 + a*y6 + b*y7 + c*y8 # regressions dem60 ~ ind60 dem65 ~ ind60 + dem60 # residual correlations y1 ~~ y5 y2 ~~ y4 + y6 y3 ~~ y7 y4 ~~ y8 y6 ~~ y8 ' # モデルのあてはめと要約 fit <- sem(model, data=PoliticalDemocracy) summary(fit, fit.measures=TRUE) sem_params1 <- parameterEstimates(fit,standardized = T) my_weapon <- function(sem_params){ # パラメータのラベルを作成 sem_params$coef_name <- apply(sem_params[,1:3],1,function(x)paste0(x,collapse = "")) # 係数が有意か否かで色を帰るためのもの。 sem_params$is_sig <- factor(sem_params$pvalue >= 0.05, levels=c("FALSE","TRUE","Fixed")) sem_params$is_sig[is.na(sem_params$is_sig)] <- factor("Fixed") # SEと信頼区間の描画のためのもの limits <- with(sem_params, aes(ymax = est+se, ymin= est-se)) limits2 <- with(sem_params, aes(ymax = ci.upper, ymin= ci.lower)) # ggplotの描画 p <- ggplot2::ggplot(sem_params, aes(y=est,x=coef_name,colour=is_sig)) p <- p + ggplot2::geom_point() + ggplot2::geom_pointrange(limits2, alpha=0.5) p <- p + ggplot2::geom_pointrange(limits, size=0.7, alpha=0.3) # 転置させるもの。変数名が長くなりがちなので,あったほうがいいかもしれない。 ggplot2::last_plot() + ggplot2::coord_flip() } # 実際にplotを作成してみる。 my_weapon(sem_params= sem_params1) # # 他のモデルでもやってみる # HS.model <- ' visual =~ x1 + x2 + x3 textual =~ x4 + x5 + x6 speed =~ x7 + x8 + x9 ' fit <- lavaan(HS.model, data=HolzingerSwineford1939, auto.var=TRUE, auto.fix.first=TRUE, auto.cov.lv.x=TRUE) summary(fit, fit.measures=TRUE) sem_params2 <- parameterEstimates(fit,standardized = T) my_weapon(sem_params = sem_params2)
出力される図ひとつめ
出力される図ふたつめ
RからJuliaを使う
Rから今流行のスクリプト言語Juliaを使いたい。
普通にJuliaにパスを通してRのsystem関数を使えばいいのだけれど,パスを通さずにRで下のような関数を書いて作業をする。jlfileにはJuliコードを書いた拡張子.jlであるファイル名を入れる。jlfilepathは作業フォルダ的なもの。jlpathはバージョンが変わったら変わる。
runjl <- function(jlfile, jlfilepath ="hogehoge"){ jlpath <- "/Applications/Julia-0.4.5.app/Contents/Resources/julia/bin/julia" system(paste0(jlpath," ",jlfilepath,jlfile)) } # こんな感じで使う。 runjl("example.jl")
わざわざ記事を書くほどのものでも無いけど,せっかくなので。
追記:普通に.bash_profileにpathを書いて,juliaへのパスを通した方が楽。
Mplusで潜在クラス分析を行うときのめも
このノートの荒っぽいまとめ。
https://www.statmodel.com/examples/webnotes/webnote14.pdf
3ステップで潜在クラスの数を決める。
- 最適な対数尤度探し
- optseedとtech 11の挿入
- optseedとtech 14の挿入
ブートストラップ尤度比検定を行うには結構時間がかかるので,いきなりやるより,手順を踏んだほうが良いようだ。
まず,クラス数kとk-1での最適な対数尤度を探す。
この時,start = 600 120;などとすると,初期値を600個ランダムに生成して,初期の対数尤度がいけてるもののうち120個を推定を行い,その中でベストなものを最終的なアウトプットにするというもの。
普通はそこまででいいのだけれど,潜在クラス数を決める場合には,誤ったseedだと検定がうまくできない。そこで,このステップではイケてる乱数のSEEDを探すということになる。
推定が終わったら,outputの中身を"seed"などで検索して,対数尤度のランキング表を見つけること。そこにイケてるseedがある。
つぎに,このseedを使って次のステップに入る。具体的には,
starts =0;!これは乱数の初期値を使わないという意味
opted = 484687;!seedを決める
として, OUTPUT: TECH1 TECH8 TECH11;とtech 11を加える。tech 11ではLo-Mendell Rubin(LMR)検定が出来る。詳細は, Lo et al.(2001)にあるらしい。ブートストラップ尤度比検定(BLRT)より計算量が少ないので,先にやってあたりをつける感じ。
最後はBLRTの登場このseedを使って次のステップに入る。具体的には,先ほどとどうように。
starts =0;
opted = 484687;
として, OUTPUT: TECH1 TECH8 TECH14;とtech 14を加える。tech 11はもはや入れなくて良い。
これで推定しておしまい。
ところが警告が出ることがある。
そんな時は,LRTSTARTS = 0 0 100 20;などとして初期値の個数を増やす。前半2つはクラス数k-1について,後半2つはクラス数kについて。
これで,いけるはず。駄目だったらモデル自体を再考するところから始める。
Mplusメモ
Mplusはミクスチャーモデルの分析もらくらくできるし,計算も結構速い。しかし,人間業が深いもので,もっと早く分析して欲しい,なんていう欲望が出てくる。
そんな時はamalysisコマンドの中にprocessors = 4;などと記述すると,複数コアで計算を回してくれる。こんな簡単に並列計算が出来るなんていい世の中になったものだ。
このPDFにいろいろと書いてある。
https://www.statmodel.com/download/handouts/MuthenV7Part1.pdf
追記:PROCESSORS = a b; ! a = # processors, (b = # threads)らしい。
コア数とスレッド数の意味については以下のページが参考になる。
CPU選びの基本:マルチコアとマルチスレッド
Windows OS入門:第3回 プロセスとスレッド (1/2) - @IT
Windowsではsetコマンドを使うとか。
Windows機で搭載CPU数がわかるコマンド - その他([技術者向] コンピューター) 解決済 | 教えて!goo
Macでコア数などを調べる場合にはターミナルでsystem_profiler SPHardwareDataTypeとすればよい。
ちなみに,mixtureモデルを分析する場合には,outputでtech14を指定するとBootstrap Likelihood Ratio Test(BLRT)をやってくれる。詳しくは文献を参照してもらいたいが,潜在クラス数を決めるのには有効な方法である。
jagsUIなるパッケージ
MCMCサンプラーのうちの1つであるJAGSをRから利用する新しいパッケージができたようだ。基本はrjagsのラッパーのようだが,並列計算の指定などがやりやすそうだ。今後,これを中心に使っていくかもしれないのでヘルプにあった例をコピペしておく。初期値の指定が関数になっているが,いいのだろうか。詳細はのちほど検証のこと。
マニュアル
https://cran.r-project.org/web/packages/jagsUI/jagsUI.pdf
作者のgithub
GitHub - kenkellner/jagsUI: Run JAGS (Just Another Gibbs Sampler) analyses from within R.
#install.packages("jagsUI") library("jagsUI") #Analyze Longley economic data in JAGS #Number employed as a function of GNP ###################################### ## 1. Collect and Package Data ## ###################################### #Load data (built into R) data(longley) head(longley) #Separate data objects gnp <- longley$GNP employed <- longley$Employed n <- length(employed) #Input data objects must be numeric, and must be #scalars, vectors, matrices, or arrays. #Package together: several possible ways #1. A named list of the objects data <- list(gnp=gnp,employed=employed,n=n) #2. A character vector of the names of the objects data <- c( gnp , employed , n ) #3. A list of names of the objects data <- list( gnp , employed , n ) ###################################### ## 2. Write model file ## ###################################### #Write a model in the BUGS language #Generate model file directly in R #(could also read in existing model file) writeLines(" model{ #Likelihood for (i in 1:n){ employed[i] ~ dnorm(mu[i], tau) mu[i] <- alpha + beta*gnp[i] } #Priors alpha ~ dnorm(0, 0.00001) beta ~ dnorm(0, 0.00001) sigma ~ dunif(0,1000) tau <- pow(sigma,-2) } ", con="model.txt") #Identify filepath of model file; #in this case in the working directory modfile <- "model.txt" ###################################### ## 3. Initialize Parameters ## ###################################### #Best to generate initial values using function inits <- function(){ list(alpha=rnorm(1,0,1),beta=rnorm(1,0,1),sigma=runif(1,0,3)) } inits #In many cases, JAGS can pick initial values automatically; #you can leave argument inits=NULL to allow this. ###################################### ## 4. Set parameters to monitor ## ###################################### #Choose parameters you want to save output for #Only parameters in this list will appear in output object #(deviance is added automatically if DIC=TRUE) #List must be specified as a character vector params <- c('alpha' , 'beta' , 'sigma' ) ###################################### ## 5. Run Analysis ## ###################################### #Call jags function; specify number of chains, number of adaptive iterations, #the length of the burn-in period, total iterations, and the thin rate. out <- jags(data = data, inits = inits, parameters.to.save = params, model.file = modfile, n.chains = 3, n.adapt = 100, n.iter = 1000, n.burnin = 500, n.thin = 2, parallel = T, n.cores=3) #Arguments will be passed to JAGS; you will see progress bars #and other information #Examine output summary out #Look at output object elements names(out) #Plot traces and posterior densities plot(out) #Plot traces traceplot(out) #Update model another 1000 iterations out <- update(out,n.iter = 1000)
追記
JAGSの並列計算を楽にしてくれる,runjagsパッケージもあるようだ。
https://cran.r-project.org/web/packages/runjags/runjags.pdf
このへんの使い方も理解しておきたい。
MCMCをするときの計算時間が尋常じゃないので,こういうものを使って少しでも楽になればいい。並列計算のやりかたを調べたり何だりしているときの虚しさたるや。。。それはもう悲惨なので。。。
学会に行く意味?
学会に行く意味とはなんだろうか。今回はこのことに関して自分にとって何が良いことなのか考えてみる。
まず,第1に締め切りがあることによる効果が挙げられる。締め切りがあることによって,適当なところで研究をまとめなければならない。コレによって,1つの区切りをつけることができる。1つの区切りがつけば,次にやることは1つの纏まった研究にするために何が必要なのかを考えることにシフトする。この作業なしに論文を書ければいいのだが,なかなかそうはいかないのが辛いところである。
次に,学会発表によって人からのフィードバックを得ることが出来る点が挙げられる。内輪の議論ではなかなかできない外からの視点を得るためにも,学会発表は行うのが好ましいように思う。
そして,単純な業績になるということである。学会発表自体はそこまで大きなものではないが,ないよりはマシである。もちろん学会発表だけではなく,キチンと論文化まで行わなければならないのだが。
最後に,これが一番大事だと思うのだが,研究のモチベーションを保つというのがあるのではないかと思う。発表してフィードバックをもらうのが憧れの先生であるならば,これは非常に嬉しいことだ。さらに,そのような先生と懇親会で話ができるというのも楽しみなことである。先生以外にも,同じあるいは隣接した領域の知り合いが増えて,交流関係が広がるのも魅力だ。こういった,人との直接的な関わりというのは,論文を読むだけでは得られない。学会発表に足を運んでこういった交流を行うことが大事なように思う。
おまけの拾い物としては,面白い発表を聞いて,最新(?)の研究動向を知ることが出来るというのもある。
デメリットとしては,時間とお金がかかることである。しかし,ある意味では大学院生にとっての就活の一環であるので,サボるほうが将来的な損失は大きい可能性がある。最も,意味のない学会に出る必要はないのだが。