Mplusでターゲット行列を指定した回転を行う方法

Analysisコマンドで，rotation = targetを指定し，
modelコマンド内で目標とする値を~で指定する。
また，各因子の最後の(*t)などとつけて，各因子負荷が1つの因子負荷行列で表すことができることを示さないといけない。
Mplusのマニュアルは，versin 8のp681などに記載されている。正直，いろいろと回転がありすぎて何が何だかわからない。
本文をtargetで検索するのがはやい。
これを使えば，任意の値を指定したターゲット行列に向けた回転を行ってくれる。
（つまり，プロクラステス(procrustes)回転ができる。）
確認的因子分析よりもゆるく，探索的因子分析よりも仮説を反映した分析ができるので，もう少し利用されてもいい方法のように感じる。

使用したデータ
https://www.statmodel.com/usersguide/chap4/ex4.1a.dat

TITLE:
DATA:
FILE IS ex4.1a.dat;
VARIABLE:
NAMES ARE y1-y12;
ANALYSIS: 
ESTIMATOR IS ML;
 ROTATION=TARGET;
 MODEL:
f1 BY y1-y12 y1~0 y2~0 y3~0 y4~0 y5~1 y6~1(*t);
f2 BY y1-y12 y7~0 y8~0 y9~0 y10~0 y11~1 y12~1(*t);
  OUTPUT:
MODINDICES CINTERVAL;
  

参考PDF
http://www.restore.ac.uk/appliedpsychometrics/Tutorial.materials/Introduction.to.Mplus/Introduction%20to%20Mplus%20session%202.pdf

Wordの小技（大文字小文字の変換，テキストボックス内の文字列の検索）

Wordで便利な小技を見つけた。
どちらも，「高度な検索と置換」を使う。
特に，ワイルドカードと，検索対象を指定するのがポイント。
「[０-９]」で検索を行うと，大文字の数字が選択される。その上で，「A->a」で「半角」を選択すればOK。

テキストボックス内の文字列の選択では，ワイルドカードの「*」を使って，検索対象を「メイン文書」ではなく「メイン文書のテキストボックス」にして検索をする。その後，検索ウィンドを閉じて「F9」をつかって，フィールドの更新をすればよい。

Wordで数字を大文字から小文字に変換する。
office-qa.com

Wordで複数のテキストボックスの文字列を選択する。
detail.chiebukuro.yahoo.co.jp

そして，フィールドについて，何処かでまとめて勉強したほうがいい気がする。
あとで，まとめて読んでみること。

pasonyu.com

pasonyu.com

フィールドについて

フィールドについて：Word(ワード)2010基本講座

平均値の信頼区間のシミュレーションプロット

Rで，信頼区間を書いていいプロットするコードのメモ。
layout関数で，上手くプット画面を分割して，さも1枚のグラフかのように見せかけるのがポイント。
けっこう可愛い感じのグラフを描画できる。

# 
# 縦向きの，信頼区間の図。
# ポイントは，レイアウトとparの使い方。
# 
layout(matrix(c(1,2),ncol=1), heights=c(1/5,4/5))
par(mar=c(0,3,1,1))
mu <- 0
sigma <- 1
n <- 5
size <- 100

curve(dnorm(x,mu,sigma/sqrt(n)),from=-3, to = 3,axes=F, xlab = "",ylab="",lty =1, col =  rgb(1,0,1,alpha = 0.7),lwd = 2)
arrows(mu, 0, mu, dnorm(mu,mu,sigma/sqrt(n)), lwd=2, angle = 0,lty =1, col =  rgb(1,0,1,alpha = 0.7))

x <- matrix(rnorm(size*n,mean = mu,sd = sigma), ncol=n)
x_mean <- rowMeans(x)
upp <- x_mean + 1.96*sigma/sqrt(n)
low <- x_mean - 1.96*sigma/sqrt(n)
is_sig <- ((upp < 0) + (low < 0)) !=1

par(mar=c(3,3,0,1))
plot(x_mean,1:size, xlim = c(-3,3),col =ifelse(is_sig, rgb(1,0,0,alpha=0.7),  rgb(0,0,1,alpha=0.7)), pch = 1- is_sig ,xlab="x平均")
arrows(x_mean, 1:size, upp, 1:size, lwd=2, angle = 90,length = 0.03,col =ifelse(is_sig, rgb(1,0,0,alpha=0.7),  rgb(0,0,1,alpha=0.7)) , lty=1 + is_sig)
arrows(x_mean, 1:size, low, 1:size, lwd=2, angle = 90,length = 0.03,col =ifelse(is_sig, rgb(1,0,0,alpha=0.7),  rgb(0,0,1,alpha=0.7)) , lty=1 + is_sig)
abline(v = mu , lwd = 2,lty=2,col =  rgb(1,0,1,alpha = 0.7))

出来上がる図は次のようなもの。

#
# 横向きの，信頼区間の図。
#

layout(matrix(c(2,1),nrow=1), widths =c(4/5,1/5))
par(mar=c(3,0,1,1)) # 下，左，上，右の順で，余白を設定する。上が0なので，いい感じでくっつけられる。
mu <- 0
sigma <- 1
n <- 10
size <- 100

# 右に90度，図を転置させる。
t_dnorm <- function(x, mu, sigma){
  cbind(x, dnorm(x,mu,sigma)) %*% matrix(c(0,-1,1,0),ncol=2,byrow=T)
}
plot(t_dnorm(seq(from=-3, to =3, by=0.1), mu, sigma/sqrt(n)),type="l",axes=F, xlab = "",ylab="",lty =1, col =  rgb(1,0,1,alpha = 0.7),lwd = 2)
arrows(0, mu, dnorm(mu,mu,sigma/sqrt(n)), mu, lwd=2, angle = 0,lty =3, col =  rgb(1,0,1,alpha = 0.7))


x <- matrix(rnorm(size*n,mean = mu,sd = sigma), ncol=n)
x_mean <- rowMeans(x)
upp <- x_mean + 1.96*sigma/sqrt(n)
low <- x_mean - 1.96*sigma/sqrt(n)
is_sig <- ((upp < 0) + (low < 0)) !=1

par(mar=c(3,3,1,0)) # この余白の設定で，絶妙に図をくっつけるのがポイント。
plot(x_mean, ylim = c(-3,3),col =ifelse(is_sig, rgb(1,0,0,alpha=0.7),  rgb(0,0,1,alpha=0.7)), pch = 1- is_sig ,xlab="x平均")
arrows(1:size,x_mean, 1:size, upp, lwd=2, angle = 90,length = 0.03,col =ifelse(is_sig, rgb(1,0,0,alpha=0.7),  rgb(0,0,1,alpha=0.7)) , lty=1 + is_sig)
arrows(1:size,x_mean, 1:size, low,lwd=2, angle = 90,length = 0.03,col =ifelse(is_sig, rgb(1,0,0,alpha=0.7),  rgb(0,0,1,alpha=0.7)) , lty=1 + is_sig)
abline(h = mu , lwd = 2,lty=3,col =  rgb(1,0,1,alpha = 0.7))

出来上がる図は次のようなもの。

日々淡々と進むこと

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

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

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

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

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

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

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

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

DATA GENERATION

「GENERATE」コマンドは，生成する従属変数のスケールを指定する。
その他，打ち切り（censored)，２値(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 )」は一文字目は数字のいち, 二文字目はローマ字小文字のエルである。これで，ロジスティックモデルによる２値のカテゴリカル変数の生成が出来る。
（閾値の数はカテゴリ数マイナス1である。）
「y1(1 p )」はプロビットによる２値のカテゴリカル変数の生成が出来る。（２パラメタ正規累積モデル， 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)を使う場合，２値と順っ女カテゴリカル変数には，プロビットモデルのみが許容される。
「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（レベル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のものに対して，レベル１が2で15このレベル2クラスターとレベル1が5でレベル２のクラスターサイズが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)にモンテカルロコマンドの詳細がいろいろ紹介されているので，今後もう少しメモを残したい。