Mplusによるモンテカルロ・シミュレーション入門
潜在変数モデルを統合的に扱うことができる統計分析ソフトであるMplus(https://www.statmodel.com/)はかなり広範な分析を実行可能である。
基本的な使い方については，インターネット上に多くの例が示されている。
例えば，以下の清水先生のサイトなどに日本語の解説がある。
norimune.net
また，マニュアルやMuthenらのサイトの掲示版も充実しており，実際的な応用に関しては困る場面は少ないだろう。
しかし，その一方でMplusには簡単にモンテカルロ・シミュレーションを実行できるということはあまり知られていないように感じられる。
自分でモンテカルロ・シミュレーションのプログラムを作成することもできるが，Mplus内でできるモンテカルロ・シミュレーションが実行できれば解析的に導出しなくても機械の力によって近似的な解を求めることができる。
これは，例えば複雑なサンプリングプロセスを実行する場合のサンプルサイズの決定や検定力の分析に利用できる。
(e.g., https://www.ncbi.nlm.nih.gov/pmc/articles/PMC3737006/)
今回はMplusを用いたモンテカルロ・シミュレーションの方法について，簡単な入門例を示す。
（今後もう少し複雑な例を示していきたい。）
以下，2因子を仮定した確認的因子分析モデルのモンテカルロ・シミュレーションのコードである。
それぞれのコマンドについて簡単に解説する。
title:
montecarlo:
names = y1-y6;
nobs = 200;
nreps = 100;
seed = 1234;
repsave = 100;
save = fa_monte.dat;
model population:
! Factor loadings
f1 by y1-y3*.6;
f2 by y4-y6*.6;
! Unieqness variances
y1-y6*.64;
! Factor means and variances
[f1-f2@0];
f1-f2@1;
! Covariance between factors
f1 with f2*0.5;
model:
f1 by y1-y3*.6;
f2 by y4-y6*.6;
y1-y6*.64;
[f1-f2@0];
f1-f2@1;
f1 with f2*0.5;
output:
tech9;
ちなみに，Mplusのマニュアル(ver.8の場合）のチャプター12がモンテカルロ・・シミュレーションの解説に当てられている。Mplusのモンテカルロシミュレーションでは実際に想定する母集団のモデルと分析モデルが異なっていてもよく，様々なModel misspesification状況でのパラメタ推定のバイアスやcoverageについて検証が可能となる。
「title:」には，いつもどおり，分析についてのコメントを付記する。この部分は分析の際に無視されるので，メモを書いておくとよい。
「montecarlo:」はモンテカルロ・シミュレーションを実行するオプションである。このコマンドを指定し，詳細な条件をその後に設定する
「names = y1-y6;」は分析に利用する変数名を記述する。「-」は多数の変数を指定する際に便利な省略企業なので積極的に利用したい。
「nobs = 200;」にはサンプルサイズを示す。
「nreps = 100;」には，シミュレーションの繰り返しに回数を示す。
「seed = 1234;」は乱数のシードであり，指定しなければデフォルトは0である。
「repsave = 100;」は生成したサンプルをどれだけ保存するかを示す。今回はすべての反復の結果を保存する。
「save = fa_monte.dat;」で生成したデータを保存する。
ちなみに，「repsave = all;」にして，[save = hoge*.dat」とするとhoge1.datからhoge100.datというアウトプットを作成できる。モンテカルロシミュレーションの結果を逐一再現したい場合にはこのようなデータを保存しておくとよい。また，保存したデータを使ったシミュレーションを行いたい場合，「DATA: FILE = hogereplist.dat; TYPE = MONTECARLO;」と指定することで先に保存したデータを利用できる。
これ以外にも，母集団でのクラス数や階層構造，欠測のパタンなどきめ細やかな状況の指定が可能となっている。
「model population:」は母集団でのモデルの指定を行う。
「 f1 by y1-y3*.6;」は因子f1が因子負荷量.6でy1~y3に負荷していることを示している。真値の指定は「*」か「@」によって行う。どちらでもかまわないようだ。基本的モデルの読み方は通常のMplusの分析の状況と変わらない。
「model:」コマンドでは，実際に推定に利用するモデルを指定する。ここでは，「*」と「@」が通常の分析のように，初期値と固定パラメタの意味になる。また，注意が必要なのは，ここで指定した値が，MSEやcoverageの計算に使われる母集団値としてみなされる。MSEの計算などでは，populationコマンドの中のものが利用されるわけではない点に注意が必要である（後のoutpu例を参照）。
「f1 by y1-y3*.6;」は因子f1がy1~y3に負荷していることを示していて，MSは.6からどれくらいの平均自乗誤差があるのかが計算される。
「output:」では表示するoutputの内容を指定できる。
以下，出力される結果について，簡単に説明する。
MODEL FIT INFORMATION
Number of Free Parameters 19
Loglikelihood
H0 Value
Mean -1613.161
Std Dev 22.561
Number of successful computations 100
Proportions Percentiles
Expected Observed Expected Observed
0.990 0.990 -1665.645 -1678.589
0.980 0.980 -1659.495 -1663.285
まず，この部分では対数尤度について，帰無仮説が正しい場合の100回のシミュレーション(Number of successful computations)の平均(Mean) と標準偏差（Std Dev）が記述されている。
PercentilesのExpected部分では理論的なパーセンタイル値が示されており，その値よりも極端な値が出る確率が，Proportionsに示されている。
例えば，Rで 1- pnorm(-1665.645, -1613.161, 22.561)と計算すると0.9899991という値が得られる。母平均が-1613.161，母標準偏差が 22.561の正規分布を想定していることがわかる。（Chi-Square Test of Model Fitではカイ二乗分布を想定している。）
Observedは観測されてた，結果に基づいた値とそれ以上の値を取った割合を示している。
次にパラメタ推定値の結果の一部を示す。
MODEL RESULTS
ESTIMATES S. E. M. S. E. 95% % Sig
Population Average Std. Dev. Average Cover Coeff
F1 BY
Y1 0.600 0.5963 0.0919 0.0863 0.0084 0.930 1.000
Y2 0.600 0.5931 0.0918 0.0861 0.0084 0.930 1.000
Y3 0.600 0.6177 0.0875 0.0876 0.0079 0.940 1.000「Population」の列は，modelコマンドの「*」や「@」で指定した値であることに注意が必要である。
「ESTIMATES Average」は100回の推定結果の平均値である。これを使って，100*(.6 - .5963)/.6を計算すると，バイアスの大きさの割合が計算できる(マニュアルp.472より）。
また「Std. Dev.」は100回の推定値の標準偏差である。
「S.E. Average」は100回の推定のそれぞれで得られた標準誤差の平均値である。
「M.S.E」は（Population - 各回の推定値）^2をシミュレーション回数で割ったもので，平均自乗誤差である。
「95％ Cover」は95％信頼区間が100回のシミュレーションのうち，Population列の値を含んでいた割合を示している。（ゆえに，.95に近いほうが，95%信頼区間がその役割を適切に表しているといえる。）
「% Sig Coeff」は推定された係数が0かどうかの検定を行った場合，何割が有意になるのかを示している。
これがMplusを用いたモンテカルロ・シミュレーションの簡単な実行方法と結果の読み方である。
モンテカルロ・シミュレーションの評価観点はバイアスの大きさ，MSEの大きさ，経験SEの大きさ，95％信頼区間の被覆率などがあり，Model resultsの部分でこれらの結果について情報を得ることができる。
データ生成の詳細な説明はマニュアルの「MONTE CARLO DATA GENERATION」の節に記載されている。基本的には，多次元正規分布から連続値を生成し，その後カテゴリカルな変数や打ち切りの変数があればそれに応じて値をカットしたり丸めるなどの処理を行っている。
また，データをRなどで作成し，それを用いてモンテカルロ・シミュレーションを行うこともできるようだ。
それに関しては，「MONTE CARLO DATA ANALYSIS」の節で説明がある。
その他のコマンドや詳細については気が向いたときに改めて記事にすることして，今日はこのあたりでおしまい。
考えることのヒント
適切な思考を展開するために、いくつか必要となる要素が存在する。それは、言葉の定義、前提、根拠だと思われる。
まず、言葉の定義を明確にすること。これにより、用語の境界を明確にできる。言葉を定義しないと、主題のズレや言葉の用法のズレに気がつくことができない。
さらに、議論の前提としてどのような状況を設定するのかを明示的に述べることが必要となる。言葉の定義とも重なるが、必要な状況設定、仮定を示さない限り具体的な議論を行うことは難しい。
推論を行うための根拠なく、「思われる」、「考えられる」のみで話を展開しても読者の納得は得られない。レゴブロックを組み合わせるように、１つ１つの推論に客観的な根拠や論拠をつけることで、確実な議論を展開することが出来る。
最高の初期値職人を目指して
潜在クラス分析(有限混合モデル）やら，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パッケージにはまだ良くわかってない使い方があるので，今後もまた調べていきたいところ。
