以下我借用 Prof. Cheung 的模擬資料跟大家簡單示範一下，如何使用R的套件來做 SEM；大家應該會發現， 對一般的SEM模型來說，R和MPlus的語法，其實差別不大：
在R的環境下要做SEM模型，一般都是用lavvan這一個套件
library(lavaan)data1 = read.table("Exercise1.dat", sep="", header=F)
names(data1) = paste0('X', 1:24)載入資料之後，lavvan和MPlus的語法，基本上是一樣的：
model = '
BE =~ X1 + X2+ X3
MS =~ X4 + X5+ X6
MM =~ X7 + X8+ X9
ST =~ X10 + X11+ X12
CA =~ X13 + X14+ X15
RV =~ X16 + X17 + X18
BGI =~ X19 + X20 + X21
RMI =~ X22 + X23 + X24
'
fit = cfa(model, data1)在這一個簡單的CFA模型裡面，它們的分析結果也是完全相同的。
summary(fit, standardized=T, fit.measures=T)## lavaan (0.5-23.1097) converged normally after 72 iterations
##
## Number of observations 527
##
## Estimator ML
## Minimum Function Test Statistic 342.056
## Degrees of freedom 224
## P-value (Chi-square) 0.000
##
## Model test baseline model:
##
## Minimum Function Test Statistic 6135.950
## Degrees of freedom 276
## P-value 0.000
##
## User model versus baseline model:
##
## Comparative Fit Index (CFI) 0.980
## Tucker-Lewis Index (TLI) 0.975
##
## Loglikelihood and Information Criteria:
##
## Loglikelihood user model (H0) -19830.716
## Loglikelihood unrestricted model (H1) -19659.688
##
## Number of free parameters 76
## Akaike (AIC) 39813.432
## Bayesian (BIC) 40137.740
## Sample-size adjusted Bayesian (BIC) 39896.495
##
## Root Mean Square Error of Approximation:
##
## RMSEA 0.032
## 90 Percent Confidence Interval 0.025 0.038
## P-value RMSEA <= 0.05 1.000
##
## Standardized Root Mean Square Residual:
##
## SRMR 0.040
##
## Parameter Estimates:
##
## Information Expected
## Standard Errors Standard
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## BE =~
## X1 1.000 0.841 0.685
## X2 1.068 0.072 14.744 0.000 0.898 0.765
## X3 1.253 0.083 15.058 0.000 1.054 0.861
## MS =~
## X4 1.000 1.822 0.922
## X5 0.459 0.060 7.709 0.000 0.836 0.560
## X6 0.497 0.071 7.039 0.000 0.906 0.444
## MM =~
## X7 1.000 1.335 0.710
## X8 1.172 0.105 11.215 0.000 1.564 0.876
## X9 0.559 0.057 9.776 0.000 0.746 0.482
## ST =~
## X10 1.000 1.166 0.667
## X11 1.108 0.093 11.975 0.000 1.292 0.759
## X12 1.078 0.090 11.983 0.000 1.256 0.722
## CA =~
## X13 1.000 1.027 0.831
## X14 0.904 0.041 22.287 0.000 0.928 0.836
## X15 1.126 0.047 23.887 0.000 1.156 0.905
## RV =~
## X16 1.000 1.335 0.927
## X17 0.941 0.034 27.508 0.000 1.255 0.846
## X18 0.999 0.033 30.083 0.000 1.333 0.888
## BGI =~
## X19 1.000 1.347 0.884
## X20 1.018 0.035 29.073 0.000 1.371 0.905
## X21 1.070 0.038 28.363 0.000 1.441 0.890
## RMI =~
## X22 1.000 0.667 0.736
## X23 1.502 0.092 16.312 0.000 1.002 0.915
## X24 1.531 0.100 15.305 0.000 1.021 0.700
##
## Covariances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## BE ~~
## MS 0.048 0.079 0.617 0.537 0.032 0.032
## MM -0.065 0.059 -1.092 0.275 -0.058 -0.058
## ST 0.132 0.055 2.421 0.015 0.135 0.135
## CA 0.182 0.045 4.025 0.000 0.211 0.211
## RV 0.224 0.058 3.889 0.000 0.199 0.199
## BGI 0.483 0.065 7.373 0.000 0.426 0.426
## RMI 0.026 0.029 0.895 0.371 0.046 0.046
## MS ~~
## MM -0.290 0.128 -2.265 0.024 -0.119 -0.119
## ST 0.098 0.114 0.857 0.391 0.046 0.046
## CA 0.189 0.093 2.035 0.042 0.101 0.101
## RV 0.118 0.119 0.995 0.320 0.049 0.049
## BGI 0.104 0.120 0.869 0.385 0.042 0.042
## RMI 0.280 0.064 4.413 0.000 0.231 0.231
## MM ~~
## ST 0.246 0.089 2.768 0.006 0.158 0.158
## CA 0.217 0.072 3.018 0.003 0.158 0.158
## RV 0.282 0.092 3.075 0.002 0.159 0.159
## BGI 0.168 0.091 1.841 0.066 0.093 0.093
## RMI 0.187 0.049 3.822 0.000 0.210 0.210
## ST ~~
## CA 0.226 0.065 3.471 0.001 0.189 0.189
## RV 0.244 0.082 2.956 0.003 0.157 0.157
## BGI 0.081 0.082 0.991 0.322 0.052 0.052
## RMI 0.084 0.042 2.005 0.045 0.108 0.108
## CA ~~
## RV 0.542 0.072 7.558 0.000 0.395 0.395
## BGI 0.357 0.069 5.163 0.000 0.258 0.258
## RMI 0.102 0.034 2.985 0.003 0.150 0.150
## RV ~~
## BGI 0.712 0.092 7.729 0.000 0.396 0.396
## RMI 0.082 0.043 1.890 0.059 0.092 0.092
## BGI ~~
## RMI 0.039 0.044 0.889 0.374 0.043 0.043
##
## Variances:
## Estimate Std.Err z-value P(>|z|) Std.lv Std.all
## .X1 0.802 0.060 13.255 0.000 0.802 0.531
## .X2 0.570 0.052 10.997 0.000 0.570 0.414
## .X3 0.386 0.057 6.825 0.000 0.386 0.258
## .X4 0.584 0.379 1.540 0.124 0.584 0.149
## .X5 1.531 0.124 12.353 0.000 1.531 0.686
## .X6 3.339 0.227 14.679 0.000 3.339 0.803
## .X7 1.757 0.177 9.951 0.000 1.757 0.496
## .X8 0.743 0.195 3.809 0.000 0.743 0.233
## .X9 1.843 0.123 14.938 0.000 1.843 0.768
## .X10 1.695 0.142 11.971 0.000 1.695 0.555
## .X11 1.225 0.138 8.877 0.000 1.225 0.423
## .X12 1.449 0.142 10.237 0.000 1.449 0.479
## .X13 0.471 0.040 11.784 0.000 0.471 0.309
## .X14 0.371 0.032 11.594 0.000 0.371 0.301
## .X15 0.294 0.039 7.604 0.000 0.294 0.180
## .X16 0.291 0.038 7.577 0.000 0.291 0.140
## .X17 0.626 0.049 12.661 0.000 0.626 0.284
## .X18 0.476 0.045 10.556 0.000 0.476 0.211
## .X19 0.509 0.046 11.070 0.000 0.509 0.219
## .X20 0.418 0.043 9.696 0.000 0.418 0.182
## .X21 0.546 0.051 10.690 0.000 0.546 0.208
## .X22 0.377 0.031 12.175 0.000 0.377 0.459
## .X23 0.197 0.047 4.174 0.000 0.197 0.164
## .X24 1.085 0.083 13.143 0.000 1.085 0.510
## BE 0.708 0.086 8.212 0.000 1.000 1.000
## MS 3.321 0.446 7.446 0.000 1.000 1.000
## MM 1.782 0.235 7.580 0.000 1.000 1.000
## ST 1.360 0.183 7.411 0.000 1.000 1.000
## CA 1.054 0.093 11.276 0.000 1.000 1.000
## RV 1.781 0.131 13.609 0.000 1.000 1.000
## BGI 1.814 0.144 12.630 0.000 1.000 1.000
## RMI 0.445 0.049 8.997 0.000 1.000 1.000我們可以自己定義兩個function來幫我們計算信效度：
CR = function(fit, xlst) {
p = parameterEstimates(fit, standardized=T)
a = sum(p$std.all[p$op == "=~" & p$rhs %in% xlst])^2
b = sum(p$std.all[p$op == "~~" & p$rhs %in% xlst])
a / (a+b) }
AVE = function(fit, xlst) {
p = parameterEstimates(fit, standardized=T)
sum(p$std.all[p$op == "=~" & p$rhs %in% xlst]^2) / length(xlst) }這兩個function的語法乍看之下有點複雜；不過，它們在所有的CFA模型裡面都是一樣的，大家只要照抄就可以了，不必自己寫。使用這兩個function，計算信效度就變得非常簡單了：
CR(fit, c('X1','X2','X3')) # 0.8162## [1] 0.8161975CR(fit, c('X4','X5','X6')) # 0.69368## [1] 0.6936807AVE(fit, c('X1','X2','X3')) # 0.59891## [1] 0.5989105AVE(fit, c('X4','X5','X6')) # 0.45379## [1] 0.4537857各位可以對照一下看看，你會發現我們做出來的結果跟講義裡面是完全一樣的。
為了模組化，新版的lavvan把一些比較複雜的程序分割到semTools 這個套件裡面。
library(semTools)我們使用這一個套件來示範一下例題5裡面 Group Invariance 的做法。 載入資料，定義好模型之後 …
data5 = read.table("Exercise5.dat", sep="", header=F)
names(data5) = c(
'Uniq1','Uniq2','Uniq3','Value1','Value2','Value3','Value4','Value5',
'App1','App2','Self1','Self2','Self3','Happy1','Happy2','PTYPE')
model = '
PUniq =~ Uniq1+Uniq2+Uniq3
ConV =~ Value1+Value2+Value3+Value4+Value5
SApp =~ App1+App2
CSelf =~ Self1+Self2+Self3
HAPP =~ Happy1+Happy2 '要檢定 Configual、Matric (loadings)、Scalar (intercepts) Invariance，只需要一行指令：
measurementInvariance(model, data=data5, group="PTYPE")##
## Measurement invariance models:
##
## Model 1 : fit.configural
## Model 2 : fit.loadings
## Model 3 : fit.intercepts
## Model 4 : fit.means
##
## Chi Square Difference Test
##
## Df AIC BIC Chisq Chisq diff Df diff Pr(>Chisq)
## fit.configural 160 4401.2 4734.6 148.49
## fit.loadings 170 4391.8 4694.9 159.10 10.618 10 0.388
## fit.intercepts 180 4378.8 4651.5 166.03 6.930 10 0.732
## fit.means 185 4408.2 4665.7 205.43 39.391 5 1.981e-07 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Fit measures:
##
## cfi rmsea cfi.delta rmsea.delta
## fit.configural 1.000 0.000 NA NA
## fit.loadings 1.000 0.000 0.000 0.000
## fit.intercepts 1.000 0.000 0.000 0.000
## fit.means 0.981 0.038 0.019 0.038大家會看到，在fit.configual、fit.loadings、fit.intercepts這三個模型裡面，我們的 \(\chi^2\)、\(\Delta\chi^2\)、\(\Delta{df}\)、\(p\)-value，跟講義第99頁是完全一樣的；Fit Measures這一部分，跟講義第96頁基本上也是一致的。最後，在fit.means這個模型裡面semTools又分別在 Laten Variables 的 Group Means上面多加了五個 Constrains，在 Chi Square Difference Test的最後一行，我們可以清楚看到顯著的差異，而這跟我們ANOVA的檢定結果也是相符合的。除非你要做 Multi-Level 的模型，R 和 MPlus 之間的差別其實不大，對於碩士班的教學來說，免費的R軟體應該是不錯的選擇才對。