チラシの裏の落書き日記

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

日々淡々と進むこと

最近思うこととしては、とにかく毎日進み続けることが何よりも大切。

ルーチンワークでやらなくてはならない雑務も多く、体力や気力も削られるかもしれないが、その中でも研究をし続け、論文をなんとしても書き続ける必要がある。

書き続けていれば、書くことの労力も減っていく。自分の中に書くということを続けることが大事。

 

スポーツをやるとこの辺りの感覚もわかる。

自分が練習しないと記録も伸びない。結局、自分との戦いだ。

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

今回は前回紹介しなかったコマンドなどについて,紹介する。しかしあまり気合が入っていないのでグダグダな感じになっている。

SCALE OF DEPENDENT VARIABLES FOR ANALYSIS

CENSORED,CATEGORICAL,NOMINAL,COUNTオプションが分析の従属変数に利用できる。指定方法は前述したgenerateと類似していいて,「CENSORED = y1(a);」などと指定する。
この指定によって,モデルと推定での変数の扱いが決まる。
独立変数としては,2値あるいは連続値が指定できる。

OPTIONS FOR DATA ANALYSIS

CLASSES, AUXILIARY, SURVIVAL のみがanalysisコマンドで利用できる。

CLASSES

CLASSESはカテゴリカル潜在変数の名前とクラス数を指定する。
これを使うためにはTYPE =MIXTUREが必要。
Betweenレベルのカテゴリカル潜在変数を使う場合には,betweenオプションを使って,宣言をする必要がある。
使い方としては,「CLASSES = c1 (2) c2 (2) c3 (3);」などと指定する。

AUXILIARY

補助変数は分析モデル内では利用されない変数である。
TYPE=MIXUREとともに,AUXILIARYオプションは自動的に3ステップアプローチ(Asparouhov and Muthén, 2014a, b 参照)を実行するために使用される。
詳細はマニュアルp.880-881を参照。

SURVIVAL

このオプションはTIMECENSOREDオプションと同時に使われないといけない。
気が向いたら追記します。
p.881-883

VARIABLES WITH SPECIAL FUNCTIONS FOR DATA GENERATION AND ANALYSIS

WITHIN

「TYPE=TWOLEVEL:」の場合は,個人レベルの変数を「WITHIN = y1 y2 x1;」と指定する。ここで指定しない変数は,within, betweenの両方でモデル化される。
「TYPE=THREELEVEL:」の場合は,「WITHIN = y1-y3 (level3) y7-y9 (level2) y4-y6;」とする。レベル1,3, 2の順に指定することに注意。y1-y3は個人レベルで,レベル1でのみモデル化される。レベル1変数は他の変数に先立ってモデル化される必要がある。y4-y6はレベル1とレベル2でモデル化される。y7-y9はレベル1とレベル3でモデル化される。
あるいは
「WITHIN = y1-y3;
WITHIN = (level2) y4-y6;
WITHIN = (level3) y7-y9;」
などと指定する。この場合には順番は気にしない。
「TYPE=CROSSCLASSIFIED:」の場合には「WITHIN = y1-y3 (level2a) y4-y6 (level2b) y7-y9;」と指定する。

BETWEEN

クラスターレベルで測定された変数の指定。
withinオプションと同様の指定方法になる。

POPULATION, COVERAGE, AND STARTING VALUES

POPULATIONと COVERAGEは前回の記事参照。
また,「STARTING = estimates.dat;」によって初期値を指定することも出来る。

SAVING DATA AND RESULTS

「REPSAVE = 1 10-15 100;」で保存するデータを指定できる。
単一のデータを保存する場合には,「SAVE = rep1.dat;」とする。
「SAVE = rep*.dat;」とすると,レプリケーションの個別のデータを指定できる。

「RESULTS = results.sav;」で各推定結果を利用できる。
しかし,すこぶる使いにくい。outファイルに以下のような出力が得られ,パラメタの順がわかる。(PARAMETER SPECIFICATIONの部分でパラメタの順も得られるので,それを参照する必要がある。)

RESULTS SAVING INFORMATION

  Order of data

    Replication number
    Parameter estimates
     (saved in order shown in Technical 1 output)
    Standard errors
     (saved in order shown in Technical 1 output)
    Chi-square : Value
    Chi-square : Degrees of Freedom
    Chi-square : P-Value
    Wald Test : Value
    Wald Test : Degrees of Freedom
    Wald Test : P-Value
    H0 Loglikelihood
    H1 Loglikelihood
    Number of Free Parameters
    Akaike (AIC)
    Bayesian (BIC)
    Sample-Size Adjusted BIC
    RMSEA : Estimate
    SRMR

推定法にbayesを指定した場合に,全てのイテレーションのパラメタの事後分布のパラメタ値が得られる。
「BPARAMETERS = bayes.dat;」

時系列解析では,ラグを指定できる。
例:「LAGGED = y (1);」
最大のラグは2である。

Mplusのモンテカルロ・シミュレーションコマンド覚書1

「MONTECARLO:」コマンドを使った場合に一緒に利用するコマンド・オプションの一部を一覧にしておく。

GENERAL SPECIFICATIONS

「NAMES = y1-y10;」は変数名を指定するオプション。

「NOBSERVATIONS = 500;」はデータ生成と分析に利用するサンプルサイズを指定するオプション。

多母集団の場合には,「NOBSERVATIONS = 500 1000;」とするとグループ1は500人,グループ2は1000人という意味になる。
さらに,「NOBSERVATIONS = 2(100) 10(500);」とすると,サンプルサイズ100人のグループが2つ,500人のグループが10つとなる。

「NGROUPS = 3;」で3つの母集団を仮定する。NGROUPSはデフォルトでは1で,1の場合にはTYPE=MIXTUREが使用できない。
なお,グループ名はg1, g2, ...と自動で名前が振られる。
TYPE=MIXTUREとともに使う場合には,既知のクラス数を指定する。ラベル「%g#1%」がはじめの既知のクラス,「%g#2%」が2番目のクラスを表す。これらのクラスを用いて,モデルを指定する。

「NREPS =100;」はシミュレーションの反復回数を表す。デフォルトは1。

「SEED = 123;」は乱数生成をするときのシード。結果の再現をしたい場合にはきちんと指定する必要がある。デフォルトは0。

DATA GENERATION

「GENERATE」コマンドは,生成する従属変数のスケールを指定する。
その他,打ち切り(censored),2値(binary),順序カテゴリ(ordered categorical; ordinal),名義(nominal),カウント(count),time-to-event変数の生成が出来る。
例えば「GENERATE = u1-u2 (1) u3 (1 p) u4 (1 l) u5 u6 (2 p) y1 (ca 1) y2 (cbi 0) u7 (n 2) u8 (ci) t1 (s 5*1);」などのように利用する。カッコ内の文字の意味を以下に説明する。

打ち切り変数は,以上・以下,インフレーションを含む・含まない,形式で生成できる。
「y1(ca 0)」は打ち切り値0より上の値を打ち切る指定。
「y1(cb 0)」は打ち切り値0より下の値を打ち切る指定。
「y1(cai 0)」はinflationを伴って,打ち切り値0より上の値を打ち切る指定。
「y1(cbi 0)」はinflationを伴って,打ち切り値0より下の値を打ち切る指定。

カテゴリカル変数の指定は以下のように行う。
「y1(1 l )」は一文字目は数字のいち, 二文字目はローマ字小文字のエルである。これで,ロジスティックモデルによる2値のカテゴリカル変数の生成が出来る。
閾値の数はカテゴリ数マイナス1である。)
「y1(1 p )」はプロビットによる2値のカテゴリカル変数の生成が出来る。(2パラメタ正規累積モデル, two-parameter normal ogive model.)
もし,lやpの指定がない場合にはロジスティックモデルが使用される。
ちなみに,「y1(1 3pl )」や「y1(1 4pl)」とすると,3パラメタ,4パラメタロジスティックモデルを指定出来る。
順序カテゴリを生成する場合には,段階反応モデル(graded response model)と同じ比例オッズモデル(proportional odds model)が使用される。
「y1(3 gpcm)」を使うと閾値3つで一般化部分採点モデル(generalized partial credit model)に従うデータを生成出来る。
また,推定方法が重み付き最小自乗法(weighted least squares estimation)を使う場合,2値と順っ女カテゴリカル変数には,プロビットモデルのみが許容される。
「y1(3 n)」は切片数3(カテゴリ数4)の名義変数を生成出来る。最後のカテゴリがレファレンスカテゴリとして切片が0と固定される。

カウント変数としては,Poisson, zero-inflated Poisson, negative binomial , zero-inflated negative binomial, zero-truncated negative binomial, negativi binomial hurdleの6種類が利用できる。
「x(c)」あるいは,「x(p)」でPoisson,「x(ci)」「x(pi)」でzero-inflated Poissonが指定できる。
「x(nb)」でnegative binomial, 「x(nbi)」でzero-inflated negative binomial モデル,「x(nbt)」でzero-truncated negative binomial modelが指定できる。「x(nbh)」でnegativi binomial hurdleを指定できる。

Time-to-event変数を指定するためには,文字「s」とベースライハザード関数の時間間隔の長さの数を指定する。もし,文字「s」飲みが指定された場合には,観測数と同じ長さの間隔数が指定される。
例えば,「t1 (s 5*1)」は長さ1の5つの時点間隔がデータ生成に利用されることを意味する。

MODEL POPULATIONでは,打ち切り,カウント変数のinflation partは数1に続くナンバーサイン(#)変数の打ち切り・カウント変数の名前を追加することによって参照される。

CUTPOINTS

カテゴリカル独立変数を指定するために使用される。
例えば,「CUTPOINTS = x1(0) x2(1);」とすると,x1は0, x2は1をカットポイントを持つ。カットポイントより低い値は0,高い値は1となる。
多母集団分析では,記号「|」で区切って指定を行う。例えば,「CUTPOINTS = x1 (0) x2 ( 1) | x1 (1) x2 (0);」とすると,グループ1と2で異なったカットポイントを指定出来る。

GENCLASSES

「GENCLASSES」オプションはカテゴリカル潜在変数を生成するために使用される。このオプションは,TYPE = MIXUREとともに利用される。
例えば,「GENCLASSES = c1 (3) c2 (2) c3 (3);」とすると,c1, c2, c3という名前の3つのカテゴリカル潜在変数を生成する。c1は3,c2は2,c3は3つのクラスを持つことを示す。

「GENCLASSES = cb (2 b);」のように文字「b」を着けると,between-levenのカテゴリ潜在変数を指定出来る。この変数は「BETWEEN」オプションを使用するbetween-level変数として必ず指定されなければならない。

NCSIZES

「NCSIZES」は,TYPE=TWOLEVEL,TYPE=THREELEVEL,TYPE=CROSSCLASSIFIEDとともに使用されなければならない。
例えば,TYPE=TWOLEVELに対しては,「NCSIZES = 3;」はデータ生成で使用されるユニークなクラスターサイズ3を指定している。
多母集団の場合には,「NCSIZES = 3 | 2;」とすると,グループ1では3つのユニークなクラスター,グループ2では2つのユニークなクラスターを生成する。
TYPE=THREELEVELでは,学生はクラスにネストされており,クラスは学校にネストされている。例えば,「NCSIZES=3[2];」と指定される。3(レベル3,学校)と2(レベル2,クラス)はユニークなクラスサイズとしてデータ生成に使用される。
多母集団の場合には,「NCSIZES = 3 [2] | 4 [3];」とすると,グループ1(レベル3で3つ,レベル2で2つのユニークなクラスタ)とグループ2(レベル3で4つ,レベル2で3つのユニークなクラス)で異なったクラスターサイズを指定できる。

TYPE =CROSSCLASSIFIEDを指定した場合には,生徒は学校と地域にネストされる。この場合,レベル1は生徒,レベル2aは学校,レベル2bは地域となる。
ここで「NCSIZES = 3 [2];」とすると,3はレベル2a, 2はレベル2bに対応する。分析のタイプがTHREELEVELかCROSSCLASSIFIEDかによって,NCSIZESの指定の意味が異なることに注意。

CSIZES

CSIZESはクラスターとサイズを指定する。
TYPE = TWOLEVELでは,「CSIZES = 100 (10) 30 (5) 15 (1);」と指定すると,サイズ10のクラスタを100クラスター,サイズ5のクラスターを30, サイズ1のクラスターを15がデータ生成に使用される。
多母集団では,「CSIZES = 100 (10) 30 (5) 15 (1) | 80 (10) 20 (5);」とすると,グループごとにクラスサイズを別に指定できる。
TYPE=THREELEVELの場合,「CSIZES = 40 [15(2) 10(5)] 30 [6(8)] 7 [20(2)];」と指定する。このようにすると,レベル3のクラスタ数が40, 30, 7となる。さらに,レベル3が40のものに対して,レベル1が2で15このレベル2クラスターとレベル1が5でレベル2のクラスターサイズが10こ,データ生成に利用される。この場合,3200の観測が得られる。以下同様。
多母集団の場合は「CSIZES = 30 [6(8)] 7 [20(2)] | 40 [5(6)] 20 [4(2)] ;」とする。

TYPE=CROSSCLASSIFIEDの場合,「CSIZES = 40 [15(2) 10(5)] 30 [6(8)] 7 [20(2)];」とすると,レベル2bのクラスタ数は角括弧の外の数(40, 30, 7),レベル2aのクラスタ数は角括弧の中,()の外の数(15, 10, 30, 20)を示す。レベル1のサイズは()内の数となる。

HAZARDC

HAZARDCは連続時間の生存解析の打ち切りプロセスのハザードを指定するために使用される。
0が打ち切られていない,1が右側打ち切りである打ち切りインディケータを生成するために使われる。
例えば,「HAZARDC = t1 (.5);」ここで,「t1」はtime-to-event変数名,.5は打ち切りのハザードを示す。

MISSING DATA GENERATION

Missingの生成するためには,PATMISSとPATPROBSを利用する。
例えば,「PATMISS = y1 (.2) y2 (.3) y3 (.1) |
y2 (.2) y3 (.1) y4 (.3) |
y3 (.1) y4 (.3);」とすると3種類の欠測パタンを生成できる。
「PATPROBS = .4 | .3 | .3;」とすると,PATMISSで指定した3つの欠測パタンが生じる確率を指定する。

MISSINGオプションは,どの従属変数で,欠測データが生成されるかを指定する。このオプションはMODLE MISSINGコマンドと同時に使用される。
Missing dataに観測独立変数を指定することはできない。
例えば,「MISSING = y1 y2 u1;」と指定される。
MODEL MISSINGは以前の記事で説明したように,ロジスティック回帰モデルで欠測確率を指定する。

次の記事で紹介する予定?

気が向いたら,以下のものについても紹介します。

  • SCALE OF DEPENDENT VARIABLES FOR ANALYSIS
  • OPTIONS FOR DATA ANALYSIS
  • VARIABLES WITH SPECIAL FUNCTIONS FOR DATA GENERATION AND ANALYSIS
  • POPULATION, COVERAGE, AND STARTING VALUES
  • SAVING DATA AND RESULTS

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

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

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

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)にモンテカルロコマンドの詳細がいろいろ紹介されているので,今後もう少しメモを残したい。

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」の節で説明がある。


その他のコマンドや詳細については気が向いたときに改めて記事にすることして,今日はこのあたりでおしまい。

考えることのヒント

適切な思考を展開するために、いくつか必要となる要素が存在する。それは、言葉の定義、前提、根拠だと思われる。

 

まず、言葉の定義を明確にすること。これにより、用語の境界を明確にできる。言葉を定義しないと、主題のズレや言葉の用法のズレに気がつくことができない。

 

さらに、議論の前提としてどのような状況を設定するのかを明示的に述べることが必要となる。言葉の定義とも重なるが、必要な状況設定、仮定を示さない限り具体的な議論を行うことは難しい。

 

推論を行うための根拠なく、「思われる」、「考えられる」のみで話を展開しても読者の納得は得られない。レゴブロックを組み合わせるように、1つ1つの推論に客観的な根拠や論拠をつけることで、確実な議論を展開することが出来る。

HHKBをWindowsで使用するときにMac風のキーバインドにする

以下のページから引用

qiita.com


ポイントは次の部分:
Google日本語入力のプロパティを開き、「一般」タブの「キー設定」に「キー設定の洗濯」という項目があると思います。ここで「MS-IME」を洗濯した上で横の「編集」ボタンを押します。

この中の「モード: 直接入力, 入力キー: Henkan」のコマンドを「IMEを有効化」にします。さらに、「モード: 入力文字なし, 入力キー: Muhenkan」のコマンドを「IMEを無効化」にします。”


これで,WindowsでもMac風のキーバインドを使用できるので,快適に文字入力ができる。