最高の初期値職人を目指して
潜在クラス分析(有限混合モデル）やら，k-meansとかいろいろと初期値に依存する分析するがある。
また，Latent variableモデリングではEMアルゴリズムを使った最尤推定を行うことが多いが，複雑なモデルになるほど解が収束しないこともしばしばある。
とくに，因子分析と潜在クラスを組合せた場合や，いわゆるGrowth mixture modelなどは収束しないだけでなく，局所解のオンパレードである。
そんな時に，いい感じの初期値を設定する方法について，まとめておく。
- モデル全体を最小二乗方による推定がうまくいくかを確認する。
上手くいけば，その解を初期値として分析を行う。だめであれば次へ。
- 一度にモデルを推定するのではなく，部分的に推定する。
このときキチンと収束することを確認する。推定方法はなんでもいい。
- すべての部分モデルを推定する。
因子分析を含んでいるのであれば，因子得点を推定して，それを顕在変数として分析する。たとえば，growth mixtureモデルの切片，傾き得点など。
ベイズ推定をした場合にも，事後分布が収束することなどをキチンと確かめる必要がある。
- その他関連する話題。
構造パラメータの段階推定について。
測定モデルを推定して，測定パラメータを固定して，潜在変数得点を推定する。
その潜在変数得点を使って，構造パラメータを推定する。この構造パラメータは一致推定量であることが示されている(光永・星野・繁桝・前川, 2005;
ci.nii.ac.jp
)。
段階推定をすることはあまり無いかもしれないが，知っておいて損ではない話題かと。
基本的には，部分的に分析してみてそれを初期値にするという戦略が重要。
初期値の指定の方法については，分析ソフトによってさまざまなので，使っているソフトのマニュアルを参考にすること。
研究について思うこと。
「研究とは何か」ということについて。
研究とはどんなゲームか。
分野にとって解く事ができる，価値のある重要な問題に効率的に解を与えるゲーム。
- 人を納得させればよい
- ゴールの見極めが大事。
- 報告はフォーマットに則って行う。
ルールをキチンと理解することがゲームを上手く進めるコツ。
list操作めも
Rのリスト操作について改めてまとめる。
特に，purrrなどを用いた操作を理解すること。
purrrの使い方についてのページ
purrr: ループ処理やapply系関数の決定版 - Heavy Watal
Hadley神のListについてのページ
Lists · R for Data Science
forループ的なlapplyやらmapの使い方が便利。
forループを使う場合には，
df <- df <- data.frame( a = rnorm(10), b = rnorm(10), c = rnorm(10), d = rnorm(10) ) temp <- numeric(length(df)) for(i in 1:length(df)){ temp[i] <- mean(df[[i]]) } temp
などと事前に結果を入れるオブジェクトを容易しないといけない。
一方，lapplyやsapplyを使うと，わりと簡潔に書ける。
lapply(df, mean) sapply(df, mean)
purrrパッケージのmapなどはlapplyを今風に処理するためのもの。
map_dfなどは型を指定して，結果をベクトルで返すもの。結果はベクトルなので注意。
library("purrr") map(df,mean) map_dbl(df, mean) map_df(df, mean)
さらに，mapやlapplyをforループ的に使うこともできる。これは言われれば納得だが，知らないと確かに書けない。今後使って行きたい。
# for loop的な使い方 map(1:3, ~df[,.]) map(1:5, runif)
また，特定の条件を満たす要素に関数を適用する，map_ifとmap_atも地味に便利。
# 0より大きい要素に対して，なんか処理をする。 -2:2 %>% map_if(.,.>0, as.factor) # 部分的に適用して，残りのリスト全体を返す。 ? map_if iris %>% map_if(is.factor, as.character) %>% str() # 全部数値 mtcars %>% str() # 1,4,5列画文字に mtcars %>% map_at(c(1, 4, 5), as.character) %>% str()
purrrパッケージにはまだ良くわかってない使い方があるので，今後もまた調べていきたいところ。
Rの関数のソースコードなど見たいときの方法
一番有名な，関数名をタイプするものでは見れないものについて。
あまり日本語の記事が無いようなきがしたので，メモを残してておく。
library("Matrix") Matrix # Matrix関数の中にある謎の関数。このままだとエラーが返ってくる。 spV2M # コロン３つなのがポイント Matrix:::spV2M # あるいは以下の関数を使うと中身がみれる。 getAnywhere("spV2M")
untar(download.packages(pkgs = "Matrix", destdir = ".", type = "source")[,2])
とすると，作業ディレクトリにパッケージがダウンロードされてuntarされて中身を見ることができるようになる。
srcフォルダにC++のソースコードなどがあるのでみてみること。
Rの組み込みのプリミティブな関数についてはpryrパッケージのshow_c_souceなんかで見れるらしい。
このページが非常に勉強になる。
stackoverflow.com
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)
出力される図ひとつめ
出力される図ふたつめ
