Mplusで推定にWLSMV(adjusted diagonally weighted least squares)を使った測定の不変性の確認の方法について。
最尤法は正規分布と連続変数という2つの強力な仮定をすることから、順序尺度や正規分布にならないものの測定の不変性について確認する場合には、WLSMVや対角重みづけ最小二乗法などが向いている。
https://www.tandfonline.com/doi/abs/10.1080/10705511.2019.1602776
Wu and Estabrook (2016)の方法を解説している。
- Wu, H., & Estabrook, R. (2016). Identification of Confirmatory Factor Analysis Models of Different Levels of Invariance for Ordered Categorical Outcomes. Psychometrika, 81(4), 1014–1045. doi:10.1007/s11336-016-9506-0 https://link.springer.com/article/10.1007/s11336-016-9506-0
lavaanでの方法はこちらを参照のこと。内容は同じである。 ides.hatenablog.com
Mplusのコードはシンプルなものもあれば、非合理的な書き方を強いるものもあり、計算能力は優れているが、コードルールが優れているとはとても言えない。今回のコードはかなりややこしいので、もし可能であるならば、lavaanで行うことをお勧めしたい。
Wu & Estabrook(2016)の方法はよくある測定の不変性の計算と異なった考え方が異なる。
2016年に発表されたWu & Estabrookの論文では、グループ間の不変性検定の文脈で、順序付けられたカテゴリー結果のCFAにおける識別の問題を議論している。彼らのアプローチは、ME/Iを実施する現在のプラクティスとは異なる。現在のプラクティスでは、まずベースラインモデルを確立し、その後、パラメータ制限を増加させていく。Wu & Estabrookによると、現在のアプローチは最適ではないという。なぜなら、潜在的な連続反応の尺度に関してベースラインモデルがどのように特定されるかに依存しているためである。著者らは特に、閾値の古いモデルがどのように特定されるかに関心を持っている。様々な評価でリッカート型の項目が人気を博していることを考えると、データをカテゴリー的に扱うことは、検討やモデル化においてより重要になるだろう。したがって、我々の目標は、ME/Iのモデルの識別とテストに関してWu & Estabrookが提案した解決策に焦点を当てる。
サンプルのデータセットやコマンドはこちらに掲載されている。 figshare.com
データは2011年TIMSS 4thのbullyingスケールから4項目を利用(Mullis et al., 2009)。恣意的に選ばれた3つの国（31＝アゼルバイジャン、40＝オーストリア、246＝フィンランド）のデータを使用。すべての項目は0（まったくない）から3（少なくとも週に一度はある）までの4段階のリッカート尺度で測定されている。この尺度の項目では1年の間に学校で以下のようなことがどのくらいの頻度で起こったかを尋ねている。サンプルサイズはアゼルバイジャン、オーストリア、フィンランドでそれぞれ3808、4457、4520だった。
コード
ベースモデル
TITLE: Configural invariance for four items on bullying scale
DATA:
FILE IS 'BULLY.dat';
VARIABLE:
NAMES ARE IDCNTRY R09A R09B R09C R09D;
USEVARIABLES ARE IDCNTRY R09A R09B R09C R09D;
CATEGORICAL ARE R09A R09B R09C R09D;
GROUPING IS IDCNTRY (31=AZE 40=AUT 246=FIN);
ANALYSIS:
ESTIMATOR = wlsmv;
H1ITERATIONS = 3000;
MODEL:
y1 BY R09A@1; ! 負荷は常に1に固定されていなければならない
y2 BY R09B@1;
y3 BY R09C@1;
y4 BY R09D@1;
F1 BY y1-y4*;
F1@1;
[F1@0];
{R09A@1 R09B@1 R09C@1 R09D@1};
[y1-y4@0];
y1-y4@0; ! これは一種のファントム変数なので、残差分散は0である
[R09A$1-R09A$3*] ;
[R09B$1-R09B$3*] ;
[R09C$1-R09C$3*] ;
[R09D$1-R09D$3*] ;
MODEL AUT:
F1 BY y1-y4*;
F1@1;
[F1@0];
{R09A@1 R09B@1 R09C@1 R09D@1};
[y1-y4@0];
y1-y4@0;
[R09A$1-R09A$3*] ;
[R09B$1-R09B$3*] ;
[R09C$1-R09C$3*] ;
[R09D$1-R09D$3*] ;
MODEL FIN:
F1 BY y1-y4*;
F1@1;
[F1@0];
{R09A@1 R09B@1 R09C@1 R09D@1};
[y1-y4@0];
y1-y4@0;
[R09A$1-R09A$3*] ;
[R09B$1-R09B$3*] ;
[R09C$1-R09C$3*] ;
[R09D$1-R09D$3*] ;
OUTPUT:
tech1 tech4;
SAVEDATA:
DIFFTEST = prop4.dif; ! カイ二乗差検定用データの保存
データファイル(BULLY.dat)の1行目には1人目の受験者のデータが入っており、ヘッダや変数名は含まれていない。USEVARIABLESのには、解析に使用する変数が記載される。これらの変数名はMODEL部分に表示されることになっている。グループ化変数IDCNTRYには、国の識別コードに関連する値が含まれている。この場合、グループ化変数には、31（Azerbaijan）、40（Australia）、24（Finland）の3つのグループが含まれている。
GROUPING IS IDCNTRY (31=AZE 40=AUT 246=FIN);
31アゼルバイジャン、40オーストリア、246フィンランドの3つのグループが含まれている。
観測された変数がカテゴリカルであることからWLSMV(ロバスト重み付き最小二乗法)を使用する。
ESTIMATOR = wlsmv; ! 順序変数の推定
h1iterations = 3000; !Mpusのデフォルトの値。欠損データを含む無制限モデルの最大反復回数。y1 BY R09A@1; y2 BY R09B@1; y3 BY R09C@1; y4 BY R09D@1;
MODELコマンドの最初の4行は、潜在変数y*を指定するためのファントム変数を示す。負荷は常に1に固定されている。
F1@1;
[F1@0];
すべてのグループの因子平均[F1@0]と分散F1@1はそれぞれ0と1に設定されている。
{R09A@1 R09B@1 R09C@1 R09D@1};
因子尺度{R09A@1 R09B@1 R09C@1 R09D@1}はすべてのグループで1に固定されている。
[y1-y4@0];
y1-y4@0;
切片の平均[y1-y4@0]と残差分散y1-y4@0も0に固定されている。
BY y1-y4 ! すべてのグループで推定される負荷量
[R09A$1-R09A$3] ; ! すべてのグループで推定されるしきい値
[R09B$1-R09B$3] ;
[R09C$1-R09C$3] ;
[R09D$1-R09D$3*] ;
このモデルはベースライン・モデルなので、閾値と負荷量はグループ間で自由に推定されることが予想されるため、アスタリスクを使用する。
SAVEDATA: DIFFTEST = prop4.dif;
分析の次のステップでANALYSISコマンドで読み込むカイ二乗差検定に必要な情報を生する。すなわち、Proposition 4 (Wu & Estabrook, 2016)のように等しい閾値を制約する場合には、file prop4.difを入力する。
Proposition 4
thresholds invarianceを求めるコードである。
TITLE: Proposition 4 (thresholds invariance) for four items on bullying scale
DATA:
FILE IS 'BULLY.dat';
VARIABLE:
NAMES ARE IDCNTRY R09A R09B R09C R09D;
USEVARIABLES ARE IDCNTRY R09A R09B R09C R09D;
CATEGORICAL ARE R09A R09B R09C R09D;
GROUPING IS IDCNTRY (31=AZE 40=AUT 246=FIN);
ANALYSIS:
ESTIMATOR = wlsmv;
DIFFTEST = prop4.dif;
H1ITERATIONS = 3000;
MODEL:
y1 BY R09A@1; ! Loadings must always be fixed to 1 here
y2 BY R09B@1;
y3 BY R09C@1;
y4 BY R09D@1;
F1 BY y1-y4*;
F1@1; !fixed to 1 in all groups
[F1@0]; !factor means fixed to 0 in all groups
{R09A@1 R09B@1 R09C@1 R09D@1};!factor scale fixed to 1 in G1, estimated in others
[y1-y4@0]; ! intercepts means fixed to 0 in group 1, others estimated
y1-y4@0; ! Because this is a kind of phantom variable, residual variance is 0
[R09A$1-R09A$3*] (T1-T3); !constrain thresholds to be equal across groups
[R09B$1-R09B$3*] (T4-T6);
[R09C$1-R09C$3*] (T7-T9);
[R09D$1-R09D$3*] (T10-T12);
MODEL AUT:
F1 BY y1-y4*;
F1@1; !fixed to 1 in all groups
[F1@0]; !factor means fixed to 0 in all groups
{R09A* R09B* R09C* R09D*};!factor scale fixed to 1 in G1, estimated in others
[y1-y4*]; !estimated in all but G1
y1-y4@0;
[R09A$1-R09A$3*] (T1-T3);
[R09B$1-R09B$3*] (T4-T6);
[R09C$1-R09C$3*] (T7-T9);
[R09D$1-R09D$3*] (T10-T12);
MODEL FIN:
F1 BY y1-y4*;
F1@1; !fixed to 1 in all groups
[F1@0]; !factor means fixed to 0 in all groups
{R09A* R09B* R09C* R09D*}; !factor scale fixed to 1 in G1, estimated in others
[y1-y4*]; !estimated in all but G1
y1-y4@0;
[R09A$1-R09A$3*] (T1-T3);
[R09B$1-R09B$3*] (T4-T6);
[R09C$1-R09C$3*] (T7-T9);
[R09D$1-R09D$3*] (T10-T12);
OUTPUT:
tech1 tech4;
SAVEDATA:
DIFFTEST = prop7.dif;
!This produces the necessary information for the chi-square difference test that
!we will read in in the "ANALYSIS" command in the next step for equal slopes and t-holds
DIFFTEST = prop4.dif;
前のベースモデルで検討されたファイルを使用する。
モデル・アイデンティフィケーションの制約として、切片の平均は、第1グループ[y1-y4@0]では0に固定され、残りのすべてのグループ[y1-y4]では推定され、Yの分散は1に固定され、すべてのグループで次のようになる。
因子平均κg=0[F1@0]と因子分散F1@1は1に固定（diag((Φg)=I). (T1-T3)...(T10-T12)を介して閾値に制約を加える。
SAVEDATA: DIFFTEST = prop7.dif;
前回と同じように今回のステップでのカイ二乗検定を保存しておく。
Proposition 7
TITLE: Proposition 7 (thresholds and loading invariance) for four items on bullying scale
DATA:
FILE IS 'BULLY.dat';
VARIABLE:
NAMES ARE IDCNTRY R09A R09B R09C R09D;
USEVARIABLES ARE IDCNTRY R09A R09B R09C R09D;
CATEGORICAL ARE R09A R09B R09C R09D;
GROUPING IS IDCNTRY (31=AZE 40=AUT 246=FIN);
ANALYSIS:
ESTIMATOR = wlsmv;
DIFFTEST = prop7.dif;
H1ITERATIONS = 3000;
MODEL:
y1 BY R09A@1; ! Loadings must always be fixed to 1 here
y2 BY R09B@1;
y3 BY R09C@1;
y4 BY R09D@1;
F1 BY y1-y4* (L1-L4); !constrain loadings across groups;
F1@1; !fixed to 1 in G1
[F1@0]; !factor means fixed to 0 in all groups
{R09A@1 R09B@1 R09C@1 R09D@1};!factor scale fixed to 1 in G1, estimated in others
[y1-y4@1]; !estimated intercepts means fixed to 1 in group 1, others estimated
y1-y4@0; ! Because this is a kind of phantom variable, residual variance is 0
[R09A$1-R09A$3*] (T1-T3); !constrain thresholds to be equal across groups
[R09B$1-R09B$3*] (T4-T6);
[R09C$1-R09C$3*] (T7-T9);
[R09D$1-R09D$3*] (T10-T12);
MODEL AUT:
F1 BY y1-y4* (L1-L4); !constrain loadings across groups;
F1*; !estimated in all but G1
[F1@0]; !factor means fixed to 0 in all groups
{R09A* R09B* R09C* R09D*};!factor scale fixed to 1 in G1, estimated in others
[y1-y4*]; !estimated in all but G1 where it is fixed to 1 (change to 0 if not working)
y1-y4@0;
[R09A$1-R09A$3*] (T1-T3);
[R09B$1-R09B$3*] (T4-T6);
[R09C$1-R09C$3*] (T7-T9);
[R09D$1-R09D$3*] (T10-T12);
MODEL FIN:
F1 BY y1-y4* (L1-L4); !constrain loadings across groups;
F1*; !estimated in all but G1
[F1@0]; !factor means fixed to 0 in all groups
{R09A* R09B* R09C* R09D*}; !factor scale fixed to 1 in G1, estimated in others
[y1-y4*]; !estimated in all but G1
y1-y4@0;
[R09A$1-R09A$3*] (T1-T3);
[R09B$1-R09B$3*] (T4-T6);
[R09C$1-R09C$3*] (T7-T9);
[R09D$1-R09D$3*] (T10-T12);
OUTPUT:
tech1 tech4;
今回のステップでも前回のステップでのカイ二乗検定の情報を利用する。
DIFFTEST = prop7.dif;
部分開放
部分不変の例は、閾値と負荷不変のモデルの修正に基づいて取られ、修正指標によって導かれる。最も厳しいモデルから始めて、等しい閾値のモデルを保持していたので、負荷のみに焦点を当てて修正指標を要求していく。部分的な測定不変性の探索を呼び出すためのMplusのシンタックスは、OUTPUTコマンドのmodindices(3.84)である。
OUTPUT: tech1 tech4 modindices(3.84);
MODEL MODIFICATION INDICES
NOTE: Modification indices for direct effects of observed dependent variables
regressed on covariates and residual covariances among observed dependent
variables may not be included. To include these, request MODINDICES (ALL).
Minimum M.I. value for printing the modification index 3.840
M.I. E.P.C. Std E.P.C. StdYX E.P.C.
Group AZE
BY Statements
Y1 BY R09A 15.285 0.252 0.183 0.183
Y1 BY R09C 90.471 -0.629 -0.457 -0.457
Y1 BY R09D 17.423 0.262 0.190 0.190
Y2 BY R09A 15.285 0.259 0.183 0.183
Y2 BY R09C 90.471 -0.648 -0.457 -0.457
Y2 BY R09D 17.423 0.270 0.190 0.190
Y3 BY R09A 15.285 0.240 0.183 0.183
Y3 BY R09C 90.471 -0.598 -0.457 -0.457
Y3 BY R09D 17.423 0.249 0.190 0.190
Y4 BY R09A 15.285 0.250 0.183 0.183
Y4 BY R09C 90.471 -0.626 -0.457 -0.457
Y4 BY R09D 17.423 0.260 0.190 0.190
F1 BY R09A 15.285 0.183 0.183 0.183
F1 BY R09C 90.471 -0.457 -0.457 -0.457
F1 BY R09D 17.423 0.190 0.190 0.190
ここで着目するのはMI値である。MI値が最も高い部分の負荷を解放すると部分不変性が得られることが多い。
R09Cのロードを解放することで、fitが最も改善されることがわかる。 Mplusでのこれらの値は、MIAZE = 90.471、MIAUT < 3.84、MIFIN = 40.327である。
モデルの修正
F1 BY R09Cを解放するため、次のように書く。MODEL AUTの部分である。
MODEL AUT: F1 BY y1-y4* (L1-L4); !constrain loadings across groups;
F1 BY y3のみ(L3)の制約を外す。
MODEL AUT: F1 BY y1 (L1); !constrain loadings across groups; F1 BY y2 (L2); F1 BY y3; !freely estimate F1 BY y4 (L4);
MODEL FIN:も同様の表記にする。
ME/Iアプローチの最新動向に関する簡単なメモ
部分的不変性アプローチに加えて、研究者たちはME/Iの達成に失敗することに対処するいくつかの方法を提案している。以下では、ME/Iを体系的に検討しようとしたいくつかの研究(Asparouhov & Muthen, 2014; Finch & French, 2018; Pokropek, Davidov, & Schmidt, 2019; Raykov et al.,2012)のうち、ME/Iの分野で期待される手法やアプローチを利用/提案しているものを、簡潔に紹介する。 MG-CFAの比較的新しいアプローチであるアライメント法（Asparouhov & Muthen, 2014）は、測定の妥当性を確立する従来の方法に代わるである。アライメント法は、典型的なMG-CFAとは異なり、測定の不変性を前提としていない。むしろ、この方法は、グループ間のパラメーター不変性を最小化する最適なソリューションを特定する。アライメント法の特徴は、すべてのモデル・パラメータ（因子平均、因子分散、負荷量、切片/閾値、および残差分散）を推定する際に、ソリューションが既存のソリューションまたはベースライン・ソリューションと同一のフィット示すことである。これは、潜在変数の平均と分散がそれぞれ0と1の値に固定されていると仮定しているコンフィギュラルコン不変モデルとは対照的であろう。このような潜在変数の標準化は，潜在変数が同じ尺度ではないことを意味し，その結果，比較することができない。アライメント法は、比較群の数が多い場合のパラメータの等質性の検定に関連する問題を克服する実用的な方法であるが、この方法は本質的には探索的なものであり、Pokropek et al.（2019）の研究で示唆されたように、その性能は条件によって異なる。 Raykovら（2012）は、Benjamini-Hochberg（BH）偽発見率法を使用し、ファミリー全体の誤り率を事前に設定したシグニフィカンスレベルに制御する、ME/Iの検査に適した多重検定手順を概説した10。著者らが述べているように、BHアプローチは、ME/Iの全体的な評価を可能にするだけでなく、異なる、より局所的な不変性のレベルをテストすることを可能にするという点で柔軟性がある（すなわち、負荷または切片の不変性のみをテストする）。 最近、Pokropekら（2019）は、ME/Iをテストするための伝統的なアプローチと新しいアプローチを調査した。大規模なシミュレーション研究において、著者らは、MG-CFA、部分的MG-CFA、マルチグループ・ベイズ構造方程式モデリング（SEM）、部分的マルチグループ・ベイズSEM、アライメント最適化を伴うMG-CFAなど、ME/Iのテストに使用される5つの手法のパフォーマンスを調べた。Pokropekらの全体的な結論は、次のように要約できる。部分的測定不変性は、多くの項目が不変である場合に、パス係数と潜在的平均を復元するのに適した（効果的な）方法かもしれない。近似測定不変性モデルは、多くのパラメータがほぼ（正確ではないが）等しい場合に、潜在的平均を復元するのに使用するのがより適切かもしれない。一方、アライメント法は、少数の不変パラメータが存在する場合に、潜在的平均を復元するのに適切かもしれない。著者らが示唆しているように、これらの手法のいくつか（特に、より柔軟性のあるもの）について、より良い収束とより効率的なアルゴリズムの観点から、今後の研究が望まれるす。 最後に、FinchとFrench（2018）は、YuanとChan（2016）、MarcoulidesとYuan（2017）、Yuan、Chan、Marcoulides、Bentler（2016）が説明したRMSEAの同等性検定アプローチの有用性を検討しましたが、YuanとChanの研究では、単一のデータセットで有望視されていました。提案された同等性検定のアプローチは、タイプI（およびII）エラー（複数可）をコントロールする上でのカイ二乗差検定の不十分さへの対応としても生まれました。そこで、FinchとFrenchは、メトリック不変性とスカラー不変性の両方を検証するシミュレーション研究を通じて、その性能を系統的に検証しました（すなわち、不変性はそれぞれ荷重と切片に存在するとモデル化されました）。著者らは、指標が正規分布していると仮定したモデルに等価アプローチを用いることを支持しており、さらに重要な点として、不変性が存在するかどうかの程度に関する有用な情報を提供する能力を指摘しています。これは、著者らが述べているように、ME/Iに対するより伝統的なアプローチとは対照的であり、一般的には不変性の有無という結論を導き出すものである。
- Asparouhov, T., & Muthén, B. O. (2014). Multiple-group factor analysis alignment. Structural Equation Modeling: A Multidisciplinary Journal, 21(4), 495–508. doi:10.1080/10705511.2014.919210
- Finch, H. W., French, B. F., & Finch, M. E. (2018). Comparison of methods for factor invariance testing of a 1-factor model with small samples and skewed latent traits. Frontiers in Psychology, 9, 332. doi:10.3389/fpsyg.2018.00332
- Pokropek, A., Davidov, E., & Schmidt, P. (2019). A Monte Carlo Simulation Study to Assess The Appropriateness of Traditional and Newer Approaches to Test for Measurement Invariance. Advanced online publication. Structural Equation Modeling: A Multidisciplinary Journal, 1–21. doi:10.1080/10705511.2018.1561293
- Raykov, T., Marcoulides, G. A., & Millsap, R. E. (2012). Factorial invar- iance in multiple populations: A multiple testing procedure. Educational and Psychological Measurement, 73, 713–727. doi:10.1177/0013164412451978
- Yuan, K.-H., & Chan, W. (2016). Measurement invariance via multigroup SEM: Issues and solutions with chi-square difference tests. Psychological Methods, 21(3), 405–426. doi:10.1037/met0000080
- Marcoulides, K. M., & Yuan, K.-H. (2017). New ways to evaluate good- ness of fit: A note on using equivalence testing to assess structural equation models. Structural Equation Modeling. A Multidisciplinary Journal, 24, 148–153.
- Yuan, K.-H., Chan, W., Marcoulides, G. A., & Bentler, P. M. (2016). Assessing structural equation models by equivalence testing with adjusted fit indexes. Structural Equation Modeling, 23(3), 319–330. doi:10.1080/10705511.2015.1065414
