MplusAutomation:R的Mplus自動化套件(2)

我們在上一篇PO文有講過，R的lavaan套件也能做多層次的SEM模型，但是它目前的功能跟Mplus相比，還有很大一段距離。短期內想要做比較複雜的多層次結構方程式模型，Mplus仍然還是比較好的選擇。不過Mplus的操作界面不像R這麼靈活，當你需要一次建立很多個模型、在模型(參數和適配度)之間中做比較、或是對估計值做進一步的分析的時候，在R的環境裡面會方便很多。

上次我們示範過MplusAutomation這個套件，它讓我們可以在R的環境裡面操作Mplus，有效地整合這兩個軟體的優點；不過我們前次所用的那個例題，其實還不能算是多層次的結構方程式模型，為了能充分測試MplusAutomation這一個套件能力，這次我從 Heck, R. H., & Thomas, S. L. (2015). An introduction to multilevel modeling techniques: MLM and SEM approaches using Mplus. Routledge 這一本書(Chapter 6 Model 6)，挑選了一個同時含有 隨機效果 和 間接效果 的 多層次潛在變項 的 結構方程式模型，從模型的結構圖裡面(圖有點大，我放這個PO文的最下方)我們可以看到，OF2和ORGDEMOS之間，除了直接效果(e)之外，還有一個由三條路徑(f-b, d-c, f-a-c)所組成的間接效果。以下我就示範，如何透過R驅動Mplus建立模型，將模型與其分析結果讀回來，並且在R的環境裡面、使用蒙地卡羅模擬估計各別路徑和整體間接、直接效果的信心區間。(目前Mplus還不支援多層次模型的boostrap，所以它只能計算各間接、直接效果的估計值和p-value，還不能提供他們的信心區間。)

處理多層次資料的Multilevel Models，已經是社會學領域不能缺少的研究工具了；管學中心會持續推動使用R來整合各種多層次的線性迴歸和結構方程式工具，如同我們在多層次分析方法整理裡面講到的，每一種工具都有它的優點和缺點；從工具應用的角度來看，我們認為單一的操作界面有助於方法間的相互比較，也方便大家在工具軟體之間做選擇。對於多層次的線性迴歸和結構方程式有興趣的老師與同學，請在這篇PO文底下留下你的姓名、email和簡單留言；我們會把後續的相關文章發送給大家，請大家多多指教。

1. (安裝)載入套件

MplusAutomation是標準的CRAN套件，直接使用install.packages()就可以安裝。除了MplusAutomation之外，我們同時也載入幾個常會用到的套件：

# install.packages('MplusAutomation')
library(MplusAutomation)
library(MASS);library(Matrix);library(reshape2)
library(magrittr);library(ggplot2)

2. Mplus的指令檔和資料檔

把已經預先備好的輸入檔和資料檔，放在工作目錄，輸入檔(mod6a.inp)的內容是像這樣：

TITLE: Model 6: Examining an indirect effect between groups;
DATA: FILE IS Ch6SEMfull.dat;
      Format is 11f8.0,13f8.2;

VARIABLE:
  Names are orgcode deptid item1 item2 item3 item4
      item5 item6 age female deptsize dept_m orgsize empstab
      orgdemos orgqual orgprod1 orgprod2 op1 op2 op3 op4
      op5 op6;
  Usevariables are orgcode item1 item2 item3 item4
      item5 item6 female orgdemos op1 op2 op3 op4 op5 op6;
  Cluster is orgcode;
  Within = female;
  Between = orgdemos op1 op2 op3 op4 op5 op6;

ANALYSIS: TYPE = twolevel;
          Estimator is MLR;

Model:
  %Between%
      orgpro by op1 op2 op3 op4 op5 op6;
      OF1 by item1 item2 item3;
      OF2 by item4 item5 item6;
      OF1 OF2 on orgpro orgdemos;
      OF2 on OF1;
      orgpro on orgdemos;
  %Within%
      F1 by item1 item2 item3;
      F2 by item4 item5 item6;
      item6 with item3;
      F1 F2 on female;
  Model indirect:
     OF2 ind orgpro;
     OF2 ind orgdemos;

OUTPUT: sampstat stdyx tech1 tech3;

SAVEDATA:
  FILE IS fscores.txt;
  SAVE=FSCORES;

Mplus的輸入檔(.inp)其實只是一般的文字檔，直接用RStudio的文字編輯器就可以進行編輯。我們使用的輸入檔是從 Heck, R. H., & Thomas, S. L. (2015). Chapter 6, Model 6 改過來的。我們使用的資料檔(Ch6SEMfull.dat)可以在 Routledge.com (出版商的網站) 下載。

3. 執行Mplus的指令檔

然後執行MplusAutomation::runModels()：

runModels("mod6a.inp")

## 
## Running model: mod6a.inp 
## System command: C:\WINDOWS\system32\cmd.exe /c cd "." && "Mplus" "mod6a.inp"

它會呼叫Mplus，然後產生mod6a.out和fscores.txt這兩個輸出檔案。由於我們在指令檔裡面有要求輸出factor scores(SAVEDATA: FILE IS fscores.txt; SAVE=FSCORES;)，除了一般的.out之外，Mplus會多產生一個輸出檔案，fscores.txt。

4. 讀取Mplus的輸出檔

Mplus的輸出檔(.out)也只是一般的文字檔，所以在RStudio裡面也可以直接打開、進行檢視或編輯。重要的是，MplusAutomation:readModel()讓我們可以一次把很多個.out檔裡面的資料讀到R的環境裡面，讓後續的分析工作變得方便，好比說 …

M = readModels()

## Reading model:  C:/1/R/mplus/ch6_6/mod6a.out

會將我們工作目錄下的mod6a.out和fscores.txt這兩個輸出檔案一併讀進M這一個資料結構裡面。

5. 讀取路徑係數(Path Coefficients)和它們的共變矩陣

M是一個list，我們用R的基本語法就可以取出我們想要的路徑係數(Path Coefficients) …

cx = coef(M,'stdyx') %>% subset(Label %in% c(
  'B OF1<-ORGPRO',         # a
  'B OF2<-ORGPRO',         # b
  'B OF2<-OF1',            # c
  'B OF1<-ORGDEMOS',       # d
  'B OF2<-ORGDEMOS',       # e
  'B ORGPRO<-ORGDEMOS'     # f
  )); cx

##                 Label    est    se  pval
## 31      B OF1<-ORGPRO  0.141 0.068 0.038
## 32      B OF2<-ORGPRO  0.126 0.039 0.001
## 33         B OF2<-OF1  0.908 0.079 0.000
## 34    B OF1<-ORGDEMOS -0.784 0.040 0.000
## 35    B OF2<-ORGDEMOS -0.032 0.082 0.695
## 36 B ORGPRO<-ORGDEMOS -0.181 0.081 0.025

mu = cx$est
names(mu) = c('a','b','c','d','e','f'); mu

##      a      b      c      d      e      f 
##  0.141  0.126  0.908 -0.784 -0.032 -0.181

共變矩陣比較麻煩一點，我們先要在TECH1裡面找到模型參數(model parameters)的索引(i) …

spec = M$tech1$parameterSpecification$BETWEEN
i = c(
  spec$beta["OF1","ORGPRO"],
  spec$beta["OF2","ORGPRO"],
  spec$beta["OF2","OF1"], 
  spec$beta["OF1","ORGDEMOS"],
  spec$beta["OF2","ORGDEMOS"],
  spec$beta["ORGPRO","ORGDEMOS"]
  ); i

## [1] 51 53 54 52 55 50

然後利用這些索引，到TECH3裡面把共變矩陣抓出來 …

SG = M$tech3$paramCov[i,i]; SG

##              51         53        54           52       55           50
## 51  0.092977100         NA        NA           NA       NA  0.001185190
## 53  0.016946800  0.0752272        NA  0.003382150       NA -0.001954650
## 54 -0.000597347 -0.0217561 0.0282116 -0.000864182       NA  0.000488999
## 52  0.000421247         NA        NA  0.033095100       NA  0.000386982
## 55  0.011700700 -0.0432453 0.0636900 -0.009438120 0.181431  0.000691532
## 50           NA         NA        NA           NA       NA  0.003541740

因為我們選取參數的順序和Mplus的參數排列順序不一樣，共變數不會完全落在矩陣的左下方；在R裡面，矩陣的重整不是太大的問題，我們稍微整理一下，就可以得到完整的共變矩陣(SG)了。

SG[is.na(SG)] = 0
SG = SG + t(SG)
diag(SG) = diag(SG)/2
SG

##              51          53           54           52           55           50
## 51  0.092977100  0.01694680 -0.000597347  0.000421247  0.011700700  0.001185190
## 53  0.016946800  0.07522720 -0.021756100  0.003382150 -0.043245300 -0.001954650
## 54 -0.000597347 -0.02175610  0.028211600 -0.000864182  0.063690000  0.000488999
## 52  0.000421247  0.00338215 -0.000864182  0.033095100 -0.009438120  0.000386982
## 55  0.011700700 -0.04324530  0.063690000 -0.009438120  0.181431000  0.000691532
## 50  0.001185190 -0.00195465  0.000488999  0.000386982  0.000691532  0.003541740

6. 蒙地卡羅模擬

蒙地卡羅模擬基本上就是從一個多變量常態分佈(Multivariate Normal Distribution)裡面做隨機抽樣，有了估計值和共變矩陣之後，只要下一行指令：

mvrn = mvrnorm(300000, mu, SG) %>% data.frame %>% 'names<-'(names(mu))

就可以隨機產生300,000組參數組合(mvrn)了。估計路徑和間接效果之前，我們可以先看一下個模型參數的分布狀況 …

melt(mvrn, variable.name='coef') %>% 
  ggplot(aes(x=value, col=coef)) +
  geom_density(aes(fill=coef), size=1, alpha=0.4)

## No id variables; using all as measure variables

接下來我們可以估計各路徑和間接、直接效果的95%信心區間(95% Confidence Interval) …

df = with(mvrn, data.frame(fb=f*b, dc=d*c, fac=f*a*c))
df$direct = mvrn$e
df$indirect = with(df, fb + dc + fac)
df$total = with(df, direct + indirect)
apply(df, 2, quantile, probs=c(0.025,0.500,0.975))

##                fb         dc         fac      direct   indirect       total
## 2.5%  -0.14405355 -1.1786577 -0.13622104 -0.86700142 -1.2244820 -1.44893386
## 50%   -0.01996024 -0.6972445 -0.01954723 -0.03107729 -0.7451795 -0.79895421
## 97.5%  0.07328559 -0.3334069  0.08363865  0.80419148 -0.3752029 -0.08797733

也很容易可以把它們的分布狀況畫出來 …

melt(df, variable.name='path') %>%
  ggplot(aes(x=value, col=path)) +
  geom_density(aes(fill=path), size=1, alpha=0.4)

## No id variables; using all as measure variables

從以上的圖形我們可以清楚的看到ORGDEMOS對OF2的效果：

Total Effect：整體而言，ORGDEMOS對OF2有顯著的負效果；
Direct Effect: 然而ORGDEMOS並不會直接對OF2產生作用，
Indirect Effect: 它的效果主要是透過
Path[d-c]: 這條間接路徑產生的。

7. 多層次潛在變數的分布

和以混合效果線性迴歸為基礎的多層次模型(e.g., GLMM, HLM)相比，Mplus的多層次結構方程式模型可以很精緻地估計多層次的潛在變數在各層次的值。我們用SAVEDATA指令(在mod6a.inp裡面)讓Mplus把這些值輸出來，然後就可以在不同族群(MALE,FEMALE)之間比較，各潛在變數(F1,F2)在各層次(level-1,level-2)的分布狀況。

sav = M$savedata
df = rbind(
  sav[,c('F1','F2','FEMALE')] %>% melt(id.vars="FEMALE"),
  sav[,c('OF1','OF2','FEMALE')] %>% melt(id.vars="FEMALE"))
df$LEVEL = ifelse(df$variable %in% c("F1", "F2"), 1, 2)
df$variable[df$variable == "OF1"] = "F1"  
df$variable[df$variable == "OF2"] = "F2"  
df$variable = as.factor(as.character(df$variable))
df$SEX = factor(ifelse(df$FEMALE==1, "FEMALE", "MALE"))
ggplot(df, aes(x=value, col=SEX)) +
  geom_density(aes(fill=SEX), size=1, alpha=0.4) +
  xlim(-25, 25) + ggtitle("Latent Variables, LEVEL ~ variable") +
  facet_grid(LEVEL ~ variable)

## Warning: Removed 5 rows containing non-finite values (stat_density).