Mplusによるモンテカルロ・シミュレーションその2

Mplusを使ったモンテカルロ・シミューレーションを実行する際に使えるコマンドを，マニュアルに沿って紹介する。

Missing
- 単純なランダムな欠測
- MAR
カテゴリカル潜在変数
マルチレベル
探索的因子分析
実際のデータの利用

マニュアルの例では生存時間解析のシミュレーションのコマンドなども提示されている。

Missing

単純なランダムな欠測

まず，missingを作成する方法を示す。

TITLE: this is an example of a Monte Carlo simulation study for a CFA with covariates
(MIMIC) with continuous factor indicators
           and patterns of missing data

MONTECARLO:
NAMES ARE y1-y4 x1 x2;
NOBSERVATIONS = 500;
NREPS = 500;
SEED = 4533;
CUTPOINTS = x2(1);
PATMISS = y1(.1) y2(.2) y3(.3) y4(1) |
                     y1(1) y2(.1) y3(.2) y4(.3);
PATPROBS = .4 | .6;
save = hoge.dat;

MODEL POPULATION:
[x1-x2@0];
x1-x2@1;
f BY y1@1 y2-y4*1;
f*.5;
y1-y4*.5;
f ON x1*1 x2*.3;

MODEL:
f BY y1@1 y2-y4*1;
f*.5;
y1-y4*.5;
f ON x1*1 x2*.3;
OUTPUT: TECH9;

ランダムな欠測を発生させるためのコマンドとして，「PATPROBS」と「PATMISS」がある。
patprobsから40%がパタン1, 60%がパタン2の形でmissingする。この例では，パタン1の欠測では，y1が10%y2が20%，y3が30%，y4が100％欠測し，パタン2ではy1が100%，y2が20%，y3が20%, y4が30%欠測する．　
100%missingは計画欠損に対応する。

「CUTPOINTS = x2(1);」はx2を正規分布で1を基準として，0, 1に分ける。

ここでは出てこないが，「GENERATE」に変数を指定することができる。「 GENERATE」は生成する従属変数のスケールを設定するためのコマンドである。
(p.861あたりを参照。）

MAR

つぎに，missing at random(MAR)を生成するためのコードは次に示すようになる。

TITLE: this is an example of a Monte Carlo simulation study for a linear growth model for a continuous outcome with missing data
where attrition is predicted by time-
           invariant covariates (MAR)

MONTECARLO:
NAMES ARE y1-y4 x1 x2;
NOBSERVATIONS = 500;
NREPS = 50;
SEED = 4533;
CUTPOINTS = x2(1);
MISSING = y1-y4;

MODEL POPULATION:
x1-x2@1;
[x1-x2@0];
i s | y1@0 y2@1 y3@2 y4@3;
[i*1 s*2];
i*1; s*.2; i WITH s*.1;
y1-y4*.5;
i ON x1*1 x2*.5;
s ON x1*.4 x2*.25;

MODEL MISSING:
[y1-y4@-1]; 
y1 ON x1*.4 x2*.2;
y2 ON x1*.8 x2*.4;
y3 ON x1*1.6 x2*.8;
 y4 ON x1*3.2 x2*1.6;

MODEL:
i s | y1@0 y2@1 y3@2 y4@3;
[i*1 s*2];
i*1; s*.2; i WITH s*.1;
y1-y4*.5;
i ON x1*1 x2*.5;
s ON x1*.4 x2*.25;

OUTPUT:TECH9;

「MONTECARLO」コマンドの「MISSING」オプションはどの変数でmissingするかを示す。
「MODEL MISSING:」の中で，母集団での欠測メカニズムとパラメタ値を指定する。
model missingではロジスティック回帰を指定する。
0はmissingでない，1はmissingを従属変数としている。
「[y1-y4@-1];」はx1, x2が0のときの切片が-1としてみなされる。
（x1,x2=0のとき，16%くらいの確率でmissingする。）

カテゴリカル潜在変数

つぎに，カテゴリカル潜在変数の作成する方法を示す。

TITLE: this is an example of a Monte Carlo simulation study for a growth mixture
           model with two classes and a misspecified
           model

MONTECARLO:
NAMES ARE u y1-y4 x;
NOBSERVATIONS = 500;
NREPS = 50;
SEED = 53487;
GENERATE = u (1);
CATEGORICAL = u;
GENCLASSES = c (2);
CLASSES = c (1);

MODEL POPULATION:
%OVERALL%
[x@0];
x@1;
i s | y1@0 y2@1 y3@2 y4@3;
i*.25 s*.04;
i WITH s*0;
y1*.4 y2*.35 y3*.3 y4*.25;
i ON x*.5;
 s ON x*.1;
 c#1 ON x*.2;

[c#1*0];
%c#1%
[u$1*1 i*3 s*.5];
 %c#2%
[u$1*-1 i*1 s*0];

ANALYSIS: TYPE = MIXTURE;

MODEL:
%OVERALL%
i s | y1@0 y2@1 y3@2 y4@3;
i*.25 s*.04;
i WITH s*0;
y1*.4 y2*.35 y3*.3 y4*.25;
i ON x*.5;
s ON x*.1;
! c#1 ON x*.2;
! [c#1*0];
u ON x;
%c#1%
[u$1*1 i*3 s*.5];
! %c#2%
! [u$1*-1 i*1 s*0];
OUTPUT: TECH9;

「GENERATE = u (1);」コマンドは，カテゴリカル変数を作るときはカッコ内で閾値の数を指定している。つまり，ここでは2値のカテゴリカル変数を生成する。「CATEGORICAL = u;」はカテゴリカル変数を指定するコマンドである。
「GENCLASSES = c (2);」と「CLASSES = c (1);」は cというカテゴリカル潜在変数を1つつくり， cは2つのクラスであることを表す。
また，ここでは，モデルのmisspesificationがある推定をおこなっている。

マルチレベル

つぎに，段階抽出を行なった場合のモデル（マルチレベルモデル）をシミュレーションする方法を示す。

TITLE: this is an example of a Monte Carlo simulation study for a two-level growth
           model for a continuous outcome (three-
           level analysis)

MONTECARLO:
NAMES ARE y1-y4 x w;
NOBSERVATIONS = 1000;
 NREPS = 50;
SEED = 58459;
 CUTPOINTS = x (1) w (0);
MISSING = y1-y4;
NCSIZES = 3;
CSIZES = 40 (5) 50 (10) 20 (15);

WITHIN = x;
BETWEEN = w;

MODEL POPULATION:
 %WITHIN%
 x@1;
 iw sw | y1@0 y2@1 y3@2 y4@3;
 y1-y4*.5;
w ON x*1;
sw ON x*.25;
iw*1; sw*.2;

 %BETWEEN%
w@1;
ib sb | y1@0 y2@1 y3@2 y4@3;
y1-y4@0;
ib ON w*.5;
sb ON w*.25;
[ib*1 sb*.5];
ib*.2; sb*.1;

MODEL MISSING:
[y1-y4@-1];
y1 ON x*.4;
y2 ON x*.8;
y3 ON x*1.6;
y4 ON x*3.2;
ANALYSIS: TYPE IS TWOLEVEL;
MODEL:
%WITHIN%
iw sw | y1@0 y2@1 y3@2 y4@3;
y1-y4*.5;
iw ON x*1;
sw ON x*.25;
iw*1; sw*.2;
%BETWEEN%
ib sb | y1@0 y2@1 y3@2 y4@3;
y1-y4@0;
ib ON w*.5;
sb ON w*.25;
[ib*1 sb*.5];
ib*.2; sb*.1;

OUTPUT: TECH9 NOCHISQUARE;

「NCSIZES = 3;」はユニークなクラスターのサイズが3であることが示されている。
「CSIZES = 40 (5) 50 (10) 20 (15);」はクラスターの数と大きさを指定する。ここでは，サイズ5のクラスターが40, 10が50,15が20ということを示している。
また，「WITHIN = x;」はユニットレベル，「BETWEEN = w;」はクラスターレベルの変数を指定する。 within,betweenのどちらにも入れなかった変数は，個人レベルで観測されて，分析ではwithin, betweenどちらにも使える。
「OUTPUT: 」の「NOCHISQUARE;」とするとχ2乗統計量を計算しない。こうすると，計算速度が上がる。

探索的因子分析

探索的因子分析のシミュレーションは次のスクリプトで行う。

TITLE: this is an example of a Monte Carlo simulation study for an exploratory factor analysis 
with continuous factor indicators

MONTECARLO:
NAMES ARE y1-y10;
NOBSERVATIONS = 500;
NREPS = 50;

MODEL POPULATION:
f1 BY y1-y7*.5;
f2 BY y4-y5*.25
y6-y10*.8;
 f1-f2@1;
 f1 WITH f2*.5;
 y1-y10*.36;

ANALYSIS: ROTATION = GEOMIN;

MODEL:
f1 BY y1-y7*.5
 y8-y10*0(*1);

f2 BY y1-y3*.0
 y4-y5*.25
y6-y10*.8(*1);

f1 WITH f2*.5;
y1-y10*.36;

OUTPUT: TECH9;

ポイントとしては，「f1 BY y1-y7*.5 y8-y10*0(*1);」の(*1)である。
これをつけるとEFA factorになるらしい。

実際のデータの利用

実際のデータの推定値を母集団の真値として，モンテカルロ・シミュレーションを実行することもできる。

TITLE:
this is an example of a real data analysis
of a CFA with covariates (MIMIC) for
continuous factor indicators where the
parameter estimates are saved for use in a
Monte Carlo simulation study

DATA:
FILE = ex12.7real.dat;

VARIABLE:
NAMES = y1-y10 x1 x2;

MODEL:
f1 BY y1@1 y2-y5*1;
f2 BY y6@1 y7-y10*1;
f1-f2*.5;
f1 WITH f2*.25;
y1-y5*.5;
[y1-y5*1];
y6-y10*.75;
[y6-y10*2];
f1 ON x1*.3 x2*.5;
f2 ON x1*.5 x2*.3;

OUTPUT: TECH1;

SAVEDATA: ESTIMATES = hoge_estimates.dat;

まず，実際のデータを分析する。ポイントは，「SAVEDATA: ESTIMATES = hoge_estimates.dat;」で，推定値を「hoge_estimates.dat」として保存出来る。

さらに，先に推定した結果を使うためのMplusスクリプトは次のようになる。

TITLE: this is an example of a Monte Carlo simulation study where parameter estimates
           saved from a real data analysis are used
           for population parameter values for data
           generation and coverage

MONTECARLO:
NAMES ARE y1-y10 x1 x2;
NOBSERVATIONS = 500;
NREPS = 500;
SEED = 45335;
POPULATION = hoge_estimates.dat;
COVERAGE = hoge_estimates.dat;

MODEL POPULATION:
f1 BY y1-y5;
f2 BY y6-y10;
f1 ON x1 x2;
f2 ON x1 x2;

MODEL:
f1 BY y1-y5;
f2 BY y6-y10;
f1 ON x1 x2;
f2 ON x1 x2;
OUTPUT: TECH9;

「POPULATION = hoge_estimates.dat;」で，先に推定したパラメタの結果が入ったデータファイル名を指定する。
これにより，推定した結果を母集団の値としてシミュレーション出来る。
「COVERAGE = hoge_estimates.dat;」によって，これをカバレッジの基準になる値を指定できる。ここでは，母集団のパラメタと同じファイルを使っている。カバレッジを指定しているので，モデルの推定の際には*とかつけなくていい点が楽。

チャプター19(p.859~p.891)にモンテカルロコマンドの詳細がいろいろ紹介されているので，今後もう少しメモを残したい。

チラシの裏の落書き日記

統計とか，研究とか，日常についての備忘録的なもの。