RからMplusが使えるようになるMplusAutomationパッケージについてのエントリ２つ目。
前回はデータをMplus形式に変換する方法について解説した。
今回も下記の条件で書いている。
- RStudioを使っている
- ファルダパスに日本語が混じっていない
今回はR上からMplusの分析を行う方法について解説する。
分析データと分析コード
パッケージのガイドに掲載されている簡単な例を実行する。 使用するデータはmtcarsである。Rにデフォルトで入っているので比較的有名なデータである。
mtcars function | R Documentation
head(mtcars)
以下のようなデータである。
mpg cyl disp hp drat wt qsec vs am gear carb Mazda RX4 21.0 6 160 110 3.90 2.620 16.46 0 1 4 4 Mazda RX4 Wag 21.0 6 160 110 3.90 2.875 17.02 0 1 4 4 Datsun 710 22.8 4 108 93 3.85 2.320 18.61 1 1 4 1 Hornet 4 Drive 21.4 6 258 110 3.08 3.215 19.44 1 0 3 1 Hornet Sportabout 18.7 8 360 175 3.15 3.440 17.02 0 0 3 2 Valiant 18.1 6 225 105 2.76 3.460 20.22 1 0 3 1
今回使用するのは４つの変数。
- mpg...マイル/ガロン(燃費)
- hp...Gross horsepower...総馬力
- wt...重さ
- disp...排気量
Mplusの分析コードは下記のように書く。
TITLE: Path Mode;
DATA: FILE = "Path_model.dat";
MODEL:
mpg ON hp;
wt ON disp;
OUTPUT:"CINTERVAL;
これをR上から実行する。
Mplusの分析コマンドをR上から実行
library(MplusAutomation)
pathmodel <- mplusObject(
TITLE = "Path Mode;",
MODEL = "
mpg ON hp;
wt ON disp;",
OUTPUT = "CINTERVAL;",
rdata = mtcars)
fit <- mplusModeler(pathmodel, modelout = "Path_model.inp", run = 1L, hashfilename = FALSE)
このスクリプトを走らせると、Rstudioの作業フォルダ内にデータの.datファイル、分析コードファイル.inp、分析結果のファイル.outが出力される。
変更するのは下記の点である。
- コロンはイコールにする
- ダブルクォーテーションで囲う
- 文末にコンマ
hashfilename= FALSEにすると、inpファイル名と同じデータファイル(.dat)が作成される。
結果をRに読み込む
Mplusの結果はpath_model.outに出力されている。
結果を見たい場合は、直接見てもよいだろう。
Rへ結果を返したいときには、下記の関数を使う。
res<-readModels()
結果をresに格納する。
走らせた分析が１つでout形式ファイルが１つの時にはこの方法が取れる。
結果は下記のように引数を指定すると表示できる。
まずは統計量。
res$summaries
下記のようにやや混沌とした形で出力される。
Mplus.version Title AnalysisType DataType Estimator Observations NGroups NDependentVars NIndependentVars NContinuousLatentVars
1 8 Path Mode; GENERAL INDIVIDUAL ML 32 1 2 2 0
Parameters ChiSqM_Value ChiSqM_DF ChiSqM_PValue ChiSqBaseline_Value ChiSqBaseline_DF ChiSqBaseline_PValue LL UnrestrictedLL CFI
1 7 15.564 2 4e-04 106.604 5 0 -101.059 -93.277 0.867
TLI AIC BIC aBIC RMSEA_Estimate RMSEA_90CI_LB RMSEA_90CI_UB RMSEA_pLT05 SRMR AICC Filename
1 0.666 216.117 226.377 204.555 0.46 0.266 0.685 0.001 0.136 220.7837 path_model.out
BICだけ取り出したいときは更にBICをつける。
res$summaries$BIC
こちらはシンプルでよい。
226.377
推定結果は次のように書くと返される。
res$parameters$unstandardized
出力。
paramHeader param est se est_se pval 1 MPG.ON HP -0.065 0.009 -6.889 0.000 2 WT.ON DISP 0.006 0.001 8.739 0.000 3 WT.WITH MPG -1.020 0.384 -2.657 0.008 4 Intercepts MPG 29.586 1.529 19.346 0.000 5 Intercepts WT 1.816 0.180 10.112 0.000 6 Residual.Variances MPG 14.045 3.524 3.986 0.000 7 Residual.Variances WT 0.209 0.056 3.751 0.000
95%信頼区間はci.を挟むんで$parameters$ci.unstandardizedと書くと返される。
texregパッケージできれいに表示
きれいに表示するにはtexregパッケージが使用できるようだ。
library(texreg)
screenreg(fit, summaries = c("Observations", "CFI", "SRMR"), single.row=TRUE)
下記のように整えて出力される。
==================================
Model 1
----------------------------------
MPG<-HP -0.06 (0.01) ***
WT<-DISP 0.01 (0.00) ***
WT<->MPG -1.02 (0.38) **
MPG<-Intercepts 29.59 (1.53) ***
WT<-Intercepts 1.82 (0.18) ***
MPG<->MPG 14.04 (3.52) ***
WT<->WT 0.21 (0.06) ***
----------------------------------
Observations 32
CFI 0.87
SRMR 0.14
==================================
*** p < 0.001, ** p < 0.01, * p < 0.05
モデルを比較する
Rのlmやglmで使用するupdate関数が利用できる。 下記のように新しいパスを追加する。
pathmodel2 <- update(pathmodel, MODEL = ~ . + "
mpg ON disp;
wt ON hp;")
fit2 <- mplusModeler(pathmodel2, modelout = "model2.inp", run = 1L, hashfilename = FALSE)
この結果を先ほどのモデルと比較する。
screenreg(list(
extract(fit, summaries = c("Observations", "CFI", "SRMR")),
extract(fit2, summaries = c("Observations", "CFI", "SRMR"))),
single.row=TRUE)
出力。
====================================================
Model 1 Model 2
----------------------------------------------------
MPG<-HP -0.06 (0.01) *** -0.02 (0.01)
WT<-DISP 0.01 (0.00) *** 0.01 (0.00) ***
WT<->MPG -1.02 (0.38) ** -0.73 (0.26) **
MPG<-Intercepts 29.59 (1.53) *** 30.74 (1.27) ***
WT<-Intercepts 1.82 (0.18) *** 1.68 (0.19) ***
MPG<->MPG 14.04 (3.52) *** 8.86 (2.21) ***
WT<->WT 0.21 (0.06) *** 0.19 (0.05) ***
MPG<-DISP -0.03 (0.01) ***
WT<-HP -0.00 (0.00)
----------------------------------------------------
Observations 32 32
CFI 0.87 1.00
SRMR 0.14 0.00
====================================================
*** p < 0.001, ** p < 0.01, * p < 0.05
95%信頼区間を併記する場合にはcis=TRUEをオプションに追加しておく。
screenreg(list(
extract(fit, cis=TRUE, summaries = c("Observations", "CFI", "SRMR")),
extract(fit2, cis=TRUE, summaries = c("Observations", "CFI", "SRMR"))),
single.row=TRUE)
結果省略。
パッケージの限界
上記のように書くと非常に便利な機能のように感じられるが、実際の所、すべての分析に対応しているわけではなさそうである。
因子分析を走らせてみたところ、結果のMplusが分析するとこまではできるが、Rへの結果のリターンができなかった。
Mplusに分析させるまでは、おそらくシンプルではないモデルであっても、可能だがRへ結果を返せないのだ。
Rは図表を整えるパッケージが揃っているので、RにMplusがデータを渡すことができれば、かなり重宝するだろう。
僕が理解できる範囲では一部の分析結果にとどまっている。
MplusAutomationパッケージでMplusコマンドを書くとややこしくなる。
それであれば前回の変換編で書いたように、データだけ書いてMplusで実行ファイルを書いた方がよいと思う。
RでMplusの分析コードを書いてMplusに計算させるというのは、実用性はそれほど高くなく、やはり、結果をRで読み込めることの方が重要性は高い。
結果の読み込みができるとパーフェクトに近いパッケージになるように思われた。
参考までに こちらで行った因子分析をMplusAutomationパッケージの形式で書き換えたものを掲載しておこう。
library(psych)
data(bfi)
FAmodel <- mplusObject(
TITLE = "Exploratory Factor Analysis by Mplus via R, using NEO data.;",
VARIABLE = "USEVARIABLES = A1-A5 C1-C5 E1-E5 N1-N5;",
ANALYSIS = "TYPE = EFA 4 6;
ESTIMATOR = ML;
ROTATION = oblimin;
PARALLEL = 50;",
PLOT = "TYPE = plot2;",
OUTPUT = "residual;",
usevariables = c("A1", "A2", "A3", "A4", "A5", "C1", "C2", "C3", "C4", "C5",
"E1", "E2", "E3", "E4", "E5", "N1", "N2", "N3", "N4", "N5",
"O1", "O2" ,"O3", "O4", "O5"),
rdata = bfi)
FA_fit <- mplusModeler(FAmodel, modelout = "bfi.inp", run = 1L, hashfilename = FALSE)
