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」の節で説明がある。
その他のコマンドや詳細については気が向いたときに改めて記事にすることして，今日はこのあたりでおしまい。
