高等量化分析
lavaan套件應用於SEM模型
Zong-Xian Huang
2018/11/12
R Markdown
This is an R Markdown presentation. Markdown is a simple formatting syntax for authoring HTML, PDF, and MS Word documents. For more details on using R Markdown see http://rmarkdown.rstudio.com.
一、開始之前
- 如果沒有相關套件，請先用
install.packages()安裝lavaan等套件 library(): 載入剛剛下載的套件，每次使用前都需要載入，但不需重新下載- lavaan套件官方網站：http://lavaan.ugent.be
#install.packages("lavaan")
#install.packages("psych")
#install.packages("foreign")
#install.packages("semPlot")
library(lavaan) ##http://lavaan.ugent.be## This is lavaan 0.6-3## lavaan is BETA software! Please report any bugs.library(psych)##
## Attaching package: 'psych'## The following object is masked from 'package:lavaan':
##
## cor2covlibrary(foreign)
library(semPlot)二、載入資料
setwd(): 設定工作環境dir(): 看工作環境內的資料有哪些rm(list=ls()): 清空Environment- 利用
read.spss()讀入清理好的.sav資料檔，並命名為mydata
setwd("/Users/huangzongxian/Desktop/lavaan")
dir()## [1] "2009_NCCU.sav" "lavaan.html" "lavaan.R" "lavaan.Rmd"
## [5] "rsconnect"rm(list=ls())
mydata <- read.spss("2009_NCCU.sav",to.data.frame = T)三、先看一下資料
dim(): 100個觀察值 (100 rows)、13個變項 (13 columns)str(): 觀察不同變數的測量尺度
dim(mydata)## [1] 100 13str(mydata)## 'data.frame': 100 obs. of 13 variables:
## $ AT1 : num 5 3 3 4 5 3 3 4 4 3 ...
## $ AT2 : num 5 4 3 5 5 5 4 4 4 4 ...
## $ AT3 : num 5 3 3 4 3 4 4 4 3 4 ...
## $ SN1 : num 2 3 4 3 1 5 2 3 4 3 ...
## $ SN2 : num 2 1 1 3 3 3 4 4 3 4 ...
## $ SN3 : num 5 4 4 2 4 4 4 3 4 3 ...
## $ PBC1: num 4 4 3 4 5 4 4 3 4 4 ...
## $ PBC2: num 2 3 3 1 2 2 2 3 3 2 ...
## $ BI1 : num 2 4 3 3 3 4 3 2 3 3 ...
## $ BI2 : num 2 3 3 2 3 3 3 2 4 3 ...
## $ B1 : num 2 2 2 2 2 2 2 2 3 2 ...
## $ B2 : num 1 2 2 1 2 2 2 2 3 3 ...
## $ B3 : num 2 1 2 1 2 2 2 2 3 2 ...head()看6筆資料summary()看基本的敘述統計
head(mydata)## AT1 AT2 AT3 SN1 SN2 SN3 PBC1 PBC2 BI1 BI2 B1 B2 B3
## 1 5 5 5 2 2 5 4 2 2 2 2 1 2
## 2 3 4 3 3 1 4 4 3 4 3 2 2 1
## 3 3 3 3 4 1 4 3 3 3 3 2 2 2
## 4 4 5 4 3 3 2 4 1 3 2 2 1 1
## 5 5 5 3 1 3 4 5 2 3 3 2 2 2
## 6 3 5 4 5 3 4 4 2 4 3 2 2 2summary(mydata)## AT1 AT2 AT3 SN1 SN2
## Min. :1.00 Min. :1.00 Min. :1.00 Min. :1.00 Min. :1.0
## 1st Qu.:3.00 1st Qu.:4.00 1st Qu.:3.00 1st Qu.:2.00 1st Qu.:3.0
## Median :4.00 Median :4.00 Median :4.00 Median :3.00 Median :3.0
## Mean :3.87 Mean :4.15 Mean :3.85 Mean :2.88 Mean :3.1
## 3rd Qu.:5.00 3rd Qu.:5.00 3rd Qu.:4.25 3rd Qu.:3.00 3rd Qu.:4.0
## Max. :5.00 Max. :5.00 Max. :5.00 Max. :5.00 Max. :5.0
## SN3 PBC1 PBC2 BI1
## Min. :1.00 Min. :1.00 Min. :1.00 Min. :1.00
## 1st Qu.:3.00 1st Qu.:3.00 1st Qu.:2.00 1st Qu.:2.00
## Median :3.50 Median :4.00 Median :2.00 Median :3.00
## Mean :3.45 Mean :3.69 Mean :2.34 Mean :2.83
## 3rd Qu.:4.00 3rd Qu.:4.00 3rd Qu.:3.00 3rd Qu.:3.25
## Max. :5.00 Max. :5.00 Max. :5.00 Max. :5.00
## BI2 B1 B2 B3
## Min. :1.00 Min. :1.00 Min. :1.00 Min. :1.00
## 1st Qu.:2.00 1st Qu.:1.75 1st Qu.:1.00 1st Qu.:1.00
## Median :3.00 Median :2.00 Median :2.00 Median :2.00
## Mean :2.84 Mean :1.80 Mean :1.81 Mean :1.74
## 3rd Qu.:4.00 3rd Qu.:2.00 3rd Qu.:2.00 3rd Qu.:2.00
## Max. :5.00 Max. :3.00 Max. :3.00 Max. :3.00- 這裡將所有變項名稱改成小寫，只是懶而已
names(mydata)<- tolower(names(mydata))四、建立模型的符號說明
=~: is measured/predicted by~: is regressed on~~: is correlated with~ 1: intercept<-: 右邊永遠等於左邊建立驗證性因素分析模型
cfamodel
cfamodel <- "
at =~ at1 + at2 + at3
sn =~ sn1 + sn2 + sn3
pbc =~ pbc1 + pbc2
bi =~ bi1 + bi2
b =~ b1 + b2 + b3 "- 建立結構方程式模型
semmodel
semmodel <- "
pbc =~ pbc1 + pbc2
sn =~ sn1 + sn2 + sn3
at =~ at1 + at2 + at3
bi =~ bi1 + bi2
b =~ b1 + b2 + b3
bi~at+sn+pbc
b~bi"- 利用
cat()來檢視模型
cat(cfamodel) # 驗證性因素分析模型##
## at =~ at1 + at2 + at3
## sn =~ sn1 + sn2 + sn3
## pbc =~ pbc1 + pbc2
## bi =~ bi1 + bi2
## b =~ b1 + b2 + b3cat(semmodel) # 結構方程式模型 ##
## pbc =~ pbc1 + pbc2
## sn =~ sn1 + sn2 + sn3
## at =~ at1 + at2 + at3
## bi =~ bi1 + bi2
## b =~ b1 + b2 + b3
## bi~at+sn+pbc
## b~bi五、主要語法function
cfa(): 驗證性因素分析functionsem(): 結構方程式function
在function之中有其他的控制碼可以調整需要的參數估計與迭代次數等
- 控制iterations數量:
control = list(iter.max=5) - 參數估計法:
estimator = c("ML","ULS","GLS","WLS")，預設為ML最大概似法，WLS就是ADF orthogonal: 是否使用正交法(也就是假設因素之間沒有相關)verbose: 要不要呈現迭代過程
cfa.fit <- cfa(cfamodel,data = mydata,
control=list(iter.max=500),
estimator = "ML",
orthogonal=FALSE,
verbose = TRUE)## Quasi-Newton steps using NLMINB:
## Objective function = 3.4960381274079513
## Objective function = Inf
## Objective function = 2.6458466664016491
## Objective function = 2.1158175289596777
## Objective function = 1.7384050842292549
## Objective function = 1.4539372098221675
## Objective function = 1.5238890970739654
## Objective function = 1.2796922684073584
## Objective function = 1.0974659521482009
## Objective function = 1.0336408797624994
## Objective function = 0.9592749551301170
## Objective function = 1.2029313729128832
## Objective function = 0.8272505736862481
## Objective function = 0.7106265015156952
## Objective function = 0.9022808427349283
## Objective function = 0.6743160849045546
## Objective function = 0.6464700322868344
## Objective function = 0.6296108343154039
## Objective function = 0.5696755523224812
## Objective function = 0.5060150990959391
## Objective function = 0.4713273932444935
## Objective function = 0.5286904113356874
## Objective function = 0.4645145596907589
## Objective function = 0.4245191571385538
## Objective function = 0.4217115330316972
## Objective function = 0.4082534329670251
## Objective function = 0.4082792638124841
## Objective function = 0.4011160417032738
## Objective function = 0.3941241082215008
## Objective function = 0.3957114025196109
## Objective function = 0.3905972346440603
## Objective function = 0.3882894160700809
## Objective function = 0.3856940334158638
## Objective function = 0.3840764062053736
## Objective function = 0.3806682126548679
## Objective function = 0.3753701023394678
## Objective function = 0.3725066467496800
## Objective function = 0.3666025853247907
## Objective function = 0.3632539484435151
## Objective function = 0.3607840949778760
## Objective function = 0.3601010640313724
## Objective function = 0.3618457382511586
## Objective function = 0.3578143814944621
## Objective function = 0.3575407759075153
## Objective function = 0.3561053142453794
## Objective function = 0.3557206451233252
## Objective function = 0.3552365530467467
## Objective function = 0.3549455801668255
## Objective function = 0.3547601233988260
## Objective function = 0.3545987858690482
## Objective function = 0.3545005172812568
## Objective function = 0.3544780343004099
## Objective function = 0.3543854281778156
## Objective function = 0.3541766942445772
## Objective function = 0.3542378684380916
## Objective function = 0.3540474687616957
## Objective function = 0.3539574184084504
## Objective function = 0.3538824542514067
## Objective function = 0.3537745736145945
## Objective function = 0.3536765855135560
## Objective function = 0.3536004442364007
## Objective function = 0.3536756316417877
## Objective function = 0.3535563834365352
## Objective function = 0.3535349439476390
## Objective function = 0.3535246562030645
## Objective function = 0.3535181703016912
## Objective function = 0.3535131924574690
## Objective function = 0.3535109132698606
## Objective function = 0.3535096757235223
## Objective function = 0.3535089491751098
## Objective function = 0.3535081811840666
## Objective function = 0.3535080170072176
## Objective function = 0.3535079079913226
## Objective function = 0.3535078759970292
## Objective function = 0.3535078643798686
## Objective function = 0.3535078594975420
## Objective function = 0.3535078591131748
## Objective function = 0.3535078590904162
## Objective function = 0.3535078590904162
## convergence status (0=ok): 0
## nlminb message says: relative convergence (4)
## number of iterations: 58
## number of function evaluations [objective, gradient]: 78 59
## Computing VCOV for se = standard ... done.
## Computing TEST for test = standard ... done.## Warning in lav_object_post_check(object): lavaan WARNING: covariance matrix of latent variables
## is not positive definite;
## use lavInspect(fit, "cov.lv") to investigate.verbose = FALSE不顯示迭代過程
sem.fit <- sem(semmodel,data = mydata,
control=list(iter.max=500),
estimator = "ML",
orthogonal=FALSE,
verbose = FALSE)## Warning in lav_object_post_check(object): lavaan WARNING: some estimated lv
## variances are negativesummary()呈現基本報表
summary(cfa.fit) # 驗證性因素分析報表## lavaan 0.6-3 ended normally after 58 iterations
##
## Optimization method NLMINB
## Number of free parameters 36
##
## Number of observations 100
##
## Estimator ML
## Model Fit Test Statistic 70.702
## Degrees of freedom 55
## P-value (Chi-square) 0.075
##
## Parameter Estimates:
##
## Information Expected
## Information saturated (h1) model Structured
## Standard Errors Standard
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|)
## at =~
## at1 1.000
## at2 0.837 0.103 8.165 0.000
## at3 0.998 0.107 9.318 0.000
## sn =~
## sn1 1.000
## sn2 1.190 0.363 3.279 0.001
## sn3 1.524 0.426 3.582 0.000
## pbc =~
## pbc1 1.000
## pbc2 0.668 0.243 2.749 0.006
## bi =~
## bi1 1.000
## bi2 1.011 0.078 12.944 0.000
## b =~
## b1 1.000
## b2 1.378 0.159 8.657 0.000
## b3 1.352 0.155 8.712 0.000
##
## Covariances:
## Estimate Std.Err z-value P(>|z|)
## at ~~
## sn -0.029 0.045 -0.638 0.523
## pbc 0.310 0.084 3.693 0.000
## bi -0.215 0.087 -2.466 0.014
## b -0.150 0.040 -3.747 0.000
## sn ~~
## pbc 0.040 0.047 0.849 0.396
## bi 0.363 0.104 3.499 0.000
## b 0.048 0.025 1.945 0.052
## pbc ~~
## bi -0.054 0.085 -0.642 0.521
## b -0.102 0.037 -2.768 0.006
## bi ~~
## b 0.148 0.045 3.293 0.001
##
## Variances:
## Estimate Std.Err z-value P(>|z|)
## .at1 0.336 0.064 5.270 0.000
## .at2 0.289 0.051 5.660 0.000
## .at3 0.133 0.044 2.991 0.003
## .sn1 0.795 0.117 6.808 0.000
## .sn2 0.697 0.107 6.512 0.000
## .sn3 0.658 0.114 5.757 0.000
## .pbc1 0.526 0.142 3.714 0.000
## .pbc2 0.747 0.119 6.300 0.000
## .bi1 0.189 0.057 3.304 0.001
## .bi2 0.121 0.054 2.227 0.026
## .b1 0.120 0.019 6.210 0.000
## .b2 0.069 0.018 3.729 0.000
## .b3 0.057 0.017 3.358 0.001
## at 0.597 0.130 4.585 0.000
## sn 0.150 0.079 1.905 0.057
## pbc 0.307 0.151 2.031 0.042
## bi 0.932 0.164 5.677 0.000
## b 0.140 0.034 4.121 0.000summary(sem.fit) # 結構方程式報表## lavaan 0.6-3 ended normally after 57 iterations
##
## Optimization method NLMINB
## Number of free parameters 33
##
## Number of observations 100
##
## Estimator ML
## Model Fit Test Statistic 89.780
## Degrees of freedom 58
## P-value (Chi-square) 0.005
##
## Parameter Estimates:
##
## Information Expected
## Information saturated (h1) model Structured
## Standard Errors Standard
##
## Latent Variables:
## Estimate Std.Err z-value P(>|z|)
## pbc =~
## pbc1 1.000
## pbc2 0.485 0.226 2.148 0.032
## sn =~
## sn1 1.000
## sn2 1.195 0.365 3.272 0.001
## sn3 1.533 0.429 3.573 0.000
## at =~
## at1 1.000
## at2 0.847 0.101 8.410 0.000
## at3 0.960 0.105 9.172 0.000
## bi =~
## bi1 1.000
## bi2 1.014 0.077 13.103 0.000
## b =~
## b1 1.000
## b2 1.387 0.163 8.506 0.000
## b3 1.368 0.160 8.550 0.000
##
## Regressions:
## Estimate Std.Err z-value P(>|z|)
## bi ~
## at -0.056 0.319 -0.176 0.860
## sn 2.521 0.827 3.049 0.002
## pbc -0.373 0.493 -0.757 0.449
## b ~
## bi 0.161 0.042 3.819 0.000
##
## Covariances:
## Estimate Std.Err z-value P(>|z|)
## pbc ~~
## sn 0.046 0.049 0.920 0.357
## at 0.335 0.087 3.861 0.000
## sn ~~
## at -0.025 0.046 -0.552 0.581
##
## Variances:
## Estimate Std.Err z-value P(>|z|)
## .pbc1 0.411 0.200 2.054 0.040
## .pbc2 0.785 0.120 6.550 0.000
## .sn1 0.797 0.117 6.815 0.000
## .sn2 0.697 0.107 6.514 0.000
## .sn3 0.657 0.114 5.753 0.000
## .at1 0.320 0.064 5.015 0.000
## .at2 0.267 0.050 5.343 0.000
## .at3 0.162 0.046 3.506 0.000
## .bi1 0.193 0.056 3.446 0.001
## .bi2 0.120 0.053 2.266 0.023
## .b1 0.122 0.020 6.217 0.000
## .b2 0.069 0.020 3.524 0.000
## .b3 0.055 0.018 3.001 0.003
## pbc 0.423 0.217 1.948 0.051
## sn 0.149 0.079 1.898 0.058
## at 0.614 0.132 4.653 0.000
## .bi -0.015 0.182 -0.084 0.933
## .b 0.113 0.028 4.020 0.000六、篩選出特定的數值
- 這裡以
semmodel作為示範 - 將標準化解存成
sfit，使後續操作較便利
standardizedSolution(sem.fit)## lhs op rhs est.std se z pvalue ci.lower ci.upper
## 1 pbc =~ pbc1 0.712 0.168 4.228 0.000 0.382 1.043
## 2 pbc =~ pbc2 0.336 0.117 2.875 0.004 0.107 0.565
## 3 sn =~ sn1 0.397 0.094 4.226 0.000 0.213 0.581
## 4 sn =~ sn2 0.484 0.090 5.374 0.000 0.307 0.660
## 5 sn =~ sn3 0.590 0.086 6.863 0.000 0.421 0.758
## 6 at =~ at1 0.811 0.045 17.894 0.000 0.722 0.900
## 7 at =~ at2 0.789 0.048 16.512 0.000 0.695 0.883
## 8 at =~ at3 0.882 0.038 22.945 0.000 0.806 0.957
## 9 bi =~ bi1 0.910 0.029 30.912 0.000 0.852 0.968
## 10 bi =~ bi2 0.943 0.027 34.773 0.000 0.890 0.996
## 11 b =~ b1 0.727 0.053 13.665 0.000 0.623 0.832
## 12 b =~ b2 0.890 0.035 25.100 0.000 0.821 0.960
## 13 b =~ b3 0.908 0.034 26.595 0.000 0.841 0.975
## 14 bi ~ at -0.046 0.260 -0.176 0.860 -0.555 0.463
## 15 bi ~ sn 1.010 0.151 6.673 0.000 0.713 1.307
## 16 bi ~ pbc -0.252 0.302 -0.835 0.404 -0.843 0.340
## 17 b ~ bi 0.419 0.091 4.581 0.000 0.240 0.598
## 18 pbc1 ~~ pbc1 0.493 0.240 2.053 0.040 0.022 0.963
## 19 pbc2 ~~ pbc2 0.887 0.078 11.312 0.000 0.734 1.041
## 20 sn1 ~~ sn1 0.842 0.075 11.296 0.000 0.696 0.989
## 21 sn2 ~~ sn2 0.766 0.087 8.795 0.000 0.595 0.937
## 22 sn3 ~~ sn3 0.652 0.101 6.438 0.000 0.454 0.851
## 23 at1 ~~ at1 0.342 0.073 4.659 0.000 0.198 0.486
## 24 at2 ~~ at2 0.377 0.075 5.004 0.000 0.230 0.525
## 25 at3 ~~ at3 0.223 0.068 3.289 0.001 0.090 0.356
## 26 bi1 ~~ bi1 0.172 0.054 3.205 0.001 0.067 0.277
## 27 bi2 ~~ bi2 0.111 0.051 2.181 0.029 0.011 0.212
## 28 b1 ~~ b1 0.471 0.077 6.082 0.000 0.319 0.623
## 29 b2 ~~ b2 0.208 0.063 3.293 0.001 0.084 0.332
## 30 b3 ~~ b3 0.176 0.062 2.840 0.005 0.055 0.297
## 31 pbc ~~ pbc 1.000 0.000 NA NA 1.000 1.000
## 32 sn ~~ sn 1.000 0.000 NA NA 1.000 1.000
## 33 at ~~ at 1.000 0.000 NA NA 1.000 1.000
## 34 bi ~~ bi -0.016 0.196 -0.084 0.933 -0.401 0.368
## 35 b ~~ b 0.825 0.077 10.772 0.000 0.675 0.975
## 36 pbc ~~ sn 0.181 0.193 0.941 0.347 -0.196 0.559
## 37 pbc ~~ at 0.657 0.164 4.005 0.000 0.336 0.979
## 38 sn ~~ at -0.084 0.150 -0.560 0.576 -0.377 0.210sfit<- standardizedSolution(sem.fit)- 只看measurement model中的 factor loading ，可用
abs()消除負號 - 將factor loading 轉換成為latent variable 可解釋的總變異量
sfit[sfit$op == "=~",]## lhs op rhs est.std se z pvalue ci.lower ci.upper
## 1 pbc =~ pbc1 0.712 0.168 4.228 0.000 0.382 1.043
## 2 pbc =~ pbc2 0.336 0.117 2.875 0.004 0.107 0.565
## 3 sn =~ sn1 0.397 0.094 4.226 0.000 0.213 0.581
## 4 sn =~ sn2 0.484 0.090 5.374 0.000 0.307 0.660
## 5 sn =~ sn3 0.590 0.086 6.863 0.000 0.421 0.758
## 6 at =~ at1 0.811 0.045 17.894 0.000 0.722 0.900
## 7 at =~ at2 0.789 0.048 16.512 0.000 0.695 0.883
## 8 at =~ at3 0.882 0.038 22.945 0.000 0.806 0.957
## 9 bi =~ bi1 0.910 0.029 30.912 0.000 0.852 0.968
## 10 bi =~ bi2 0.943 0.027 34.773 0.000 0.890 0.996
## 11 b =~ b1 0.727 0.053 13.665 0.000 0.623 0.832
## 12 b =~ b2 0.890 0.035 25.100 0.000 0.821 0.960
## 13 b =~ b3 0.908 0.034 26.595 0.000 0.841 0.975abs(sfit[sfit$op == "=~","est.std"])## [1] 0.7122840 0.3357842 0.3969866 0.4837671 0.5896533 0.8109131 0.7890624
## [8] 0.8815770 0.9100829 0.9426319 0.7273745 0.8900264 0.9077605mean(abs(sfit[sfit$op == "=~","est.std"])) ##整體可解釋的總變異量## [1] 0.7213772- 只看structural model中的 correlations
sfit[sfit$op == "~" & sfit$lhs != sfit$rhs,]## lhs op rhs est.std se z pvalue ci.lower ci.upper
## 14 bi ~ at -0.046 0.260 -0.176 0.860 -0.555 0.463
## 15 bi ~ sn 1.010 0.151 6.673 0.000 0.713 1.307
## 16 bi ~ pbc -0.252 0.302 -0.835 0.404 -0.843 0.340
## 17 b ~ bi 0.419 0.091 4.581 0.000 0.240 0.598- 透過
inspect()可以來看更多模型中的資訊 - 變項的解釋能力
r2
inspect(sem.fit,"r2")## pbc1 pbc2 sn1 sn2 sn3 at1 at2 at3 bi1 bi2 b1 b2
## 0.507 0.113 0.158 0.234 0.348 0.658 0.623 0.777 0.828 0.889 0.529 0.792
## b3 bi b
## 0.824 NA 0.175sampstat用來看模型產出的共變矩陣
inspect(sem.fit,"sampstat")## $cov
## pbc1 pbc2 sn1 sn2 sn3 at1 at2 at3 bi1 bi2
## pbc1 0.834
## pbc2 0.205 0.884
## sn1 0.043 0.061 0.946
## sn2 0.061 0.046 0.212 0.910
## sn3 0.080 -0.053 0.324 0.185 1.008
## at1 0.310 0.154 -0.126 -0.027 0.089 0.933
## at2 0.336 0.029 -0.072 -0.005 0.042 0.530 0.708
## at3 0.304 0.231 -0.138 0.015 -0.082 0.590 0.492 0.728
## bi1 -0.013 -0.062 0.380 0.407 0.576 -0.202 -0.124 -0.236 1.121
## bi2 -0.070 -0.046 0.361 0.456 0.542 -0.151 -0.166 -0.254 0.943 1.074
## b1 -0.102 -0.162 0.016 -0.020 0.090 -0.136 -0.100 -0.190 0.126 0.138
## b2 -0.099 -0.145 0.087 0.039 0.126 -0.215 -0.142 -0.218 0.208 0.200
## b3 -0.111 -0.152 0.069 0.056 0.127 -0.164 -0.131 -0.219 0.176 0.228
## b1 b2 b3
## pbc1
## pbc2
## sn1
## sn2
## sn3
## at1
## at2
## at3
## bi1
## bi2
## b1 0.260
## b2 0.192 0.334
## b3 0.188 0.261 0.312七、信度與效度指標
- 多元相關平方: SMC 反應了個別測量變項受到潛在變項影響的程度，所以會有三個數字
- 信度指標: CR (建構信度)
- 效度指標: AVE (聚合效度)
SMC <- function(fit, xlst) {
p <- parameterEstimates(fit, standardized=T)
a <- p$std.all[p$op == "=~" & p$rhs %in% xlst]^2
b <- p$std.all[p$op == "~~" & p$rhs %in% xlst]
a/(a+b)
}
SMC(sem.fit,c("at1","at2","at3"))## [1] 0.6575800 0.6226195 0.7771780CR <- 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) }
CR(sem.fit,c("at1","at2","at3")) # 建構信度## [1] 0.8672498AVE <- function(fit, xlst) {
p <- parameterEstimates(fit, standardized=T)
sum(p$std.all[p$op == "=~" & p$rhs %in% xlst]^2) / length(xlst) }
AVE(sem.fit,c("at1","at2","at3")) # 聚合效度## [1] 0.6857925八、模型評估指標
fitmeasure(): 全部都給你fitmeasure(object,c("chisq", "df",...)): 依需求分別挑出幾個重要數值- 基本上課本上面的指標都找得到，僅有「NCP」、「CAIC」找不到
- 這裡就先以
semmodel作為示範
fitmeasures(sem.fit)## npar fmin chisq
## 33.000 0.449 89.780
## df pvalue baseline.chisq
## 58.000 0.005 669.776
## baseline.df baseline.pvalue cfi
## 78.000 0.000 0.946
## tli nnfi rfi
## 0.928 0.928 0.820
## nfi pnfi ifi
## 0.866 0.644 0.948
## rni logl unrestricted.logl
## 0.946 -1324.520 -1279.630
## aic bic ntotal
## 2715.039 2801.010 100.000
## bic2 rmsea rmsea.ci.lower
## 2696.788 0.074 0.042
## rmsea.ci.upper rmsea.pvalue rmr
## 0.103 0.101 0.065
## rmr_nomean srmr srmr_bentler
## 0.065 0.116 0.116
## srmr_bentler_nomean crmr crmr_nomean
## 0.116 0.125 0.125
## srmr_mplus srmr_mplus_nomean cn_05
## 0.116 0.116 86.518
## cn_01 gfi agfi
## 96.734 0.879 0.810
## pgfi mfi ecvi
## 0.560 0.853 1.558fitMeasures(sem.fit, c("chisq", "df", "pvalue"))## chisq df pvalue
## 89.780 58.000 0.005fitmeasures(sem.fit,c("gfi","agfi","pgfi","nfi","nnfi")) # 模型契合指標## gfi agfi pgfi nfi nnfi
## 0.879 0.810 0.560 0.866 0.928fitmeasures(sem.fit,c("cfi","rmsea","aic","caic","cn_01")) # 替代指標## cfi rmsea aic cn_01
## 0.946 0.074 2715.039 96.734fitmeasures(sem.fit,c("rmr","srmr")) # 殘差指數## rmr srmr
## 0.065 0.116九、畫圖囉！
- 使用
semPaths(): 結構方程式跟驗證性因素分析都可以畫出來
semPaths(sem.fit) # 很空、什麼都沒有、會被老師揍
- 加上相關控制碼，讓圖形更符合需求
- 使用
help()來看更多功能！
semPaths(sem.fit,
what = "equality",
whatLabels = "std", # std = factor loading
groups = "latents", # 依latent 分組上色
pastel = TRUE, # 柔和色
mar = c(2, 2, 2, 2), # 邊界
intercepts = TRUE, # 殘差顯示
edge.label.cex = 1, # label 的字體大小
sizeMan = 7, # measurement variable 顯示大小
sizeLat = 10, # latent variable 顯示大小
exoCov = FALSE, # 相互影響的線不要
optimizeLatRes = TRUE, # 讓latent 的殘差長整齊一點
style="lisrel", # 讓殘差格式跟LISREL一樣
layout="tree",
rotation = 2, # 橫向
structural=FALSE #set TURE if you only need structural model
) 
semPaths(cfa.fit,
what = "equality",
whatLabels = "std", # est = 係數,std = factor loading
groups = "latents", # 依latent 分組上色
pastel = TRUE, # 柔和色
mar = c(3, 1, 3, 1), # 邊界
intercepts = TRUE, # 殘差顯示
edge.label.cex = 0.8, # label 的字體大小
sizeMan = 7, # measurement variable 顯示大小
sizeLat = 10, # latent variable 顯示大小
exoCov = FALSE, # 相互影響的線不要
optimizeLatRes = TRUE, # 讓latent 的殘差長整齊一點
style="lisrel", # 讓殘差格式跟LISREL一樣
layout="tree"
) 
help("semPaths")十、模式修飾(model modification)
modificationindices()用於模式修飾
cfamod_ind<- modificationindices(cfa.fit)
semmod_ind<- modificationindices(sem.fit)- 挑選出MI值前10高的改善手段
head(cfamod_ind[order(cfamod_ind$mi,decreasing = T),],10)## lhs op rhs mi epc sepc.lv sepc.all sepc.nox
## 111 at2 ~~ pbc2 9.840 -0.167 -0.167 -0.358 -0.358
## 165 bi1 ~~ b3 8.341 -0.048 -0.048 -0.462 -0.462
## 168 bi2 ~~ b3 8.338 0.046 0.046 0.552 0.552
## 167 bi2 ~~ b2 6.030 -0.041 -0.041 -0.448 -0.448
## 164 bi1 ~~ b2 5.449 0.041 0.041 0.357 0.357
## 110 at2 ~~ pbc1 4.989 0.109 0.109 0.281 0.281
## 98 at1 ~~ sn3 4.607 0.116 0.116 0.247 0.247
## 121 at3 ~~ pbc2 4.600 0.104 0.104 0.331 0.331
## 104 at1 ~~ b2 4.586 -0.045 -0.045 -0.297 -0.297
## 86 b =~ at3 4.352 -0.436 -0.163 -0.191 -0.191head(semmod_ind[order(semmod_ind$mi,decreasing = T),],10)## lhs op rhs mi epc sepc.lv sepc.all sepc.nox
## 177 b ~ at 16.005 -0.207 -0.437 -0.437 -0.437
## 178 b ~ sn 15.160 -1.543 -1.606 -1.606 -1.606
## 179 b ~ pbc 13.618 -0.256 -0.448 -0.448 -0.448
## 185 sn ~ b 12.997 -1.484 -1.426 -1.426 -1.426
## 172 sn ~~ b 12.997 -0.168 -1.295 -1.295 -1.295
## 107 pbc2 ~~ at2 9.258 -0.163 -0.163 -0.355 -0.355
## 174 at ~~ b 9.053 -0.079 -0.299 -0.299 -0.299
## 181 at ~ b 9.053 -0.696 -0.330 -0.330 -0.330
## 162 bi1 ~~ b3 8.920 -0.050 -0.050 -0.485 -0.485
## 88 b =~ at3 8.027 -0.432 -0.160 -0.188 -0.188- 僅挑選出MI值大於10的改善手段
cfamod_ind[cfamod_ind$mi>10,] # cfa模型中沒有MI值大於10## [1] lhs op rhs mi epc sepc.lv sepc.all sepc.nox
## <0 rows> (or 0-length row.names)semmod_ind[semmod_ind$mi>10,]## lhs op rhs mi epc sepc.lv sepc.all sepc.nox
## 172 sn ~~ b 12.997 -0.168 -1.295 -1.295 -1.295
## 177 b ~ at 16.005 -0.207 -0.437 -0.437 -0.437
## 178 b ~ sn 15.160 -1.543 -1.606 -1.606 -1.606
## 179 b ~ pbc 13.618 -0.256 -0.448 -0.448 -0.448
## 185 sn ~ b 12.997 -1.484 -1.426 -1.426 -1.426