df_f <-readRDS("data/items35_N1200.rds")source('scr/rename_specific_items.R') # We decided to relabel some constructs for claritydf_f <-rename_specific_items(df_f)# F1 should be in Motivation, relabelled as M9df_f$item_label[df_f$item_label=="F1"] <-"M9"
Report CFA Model Fit for each Construct
For each construct, report (1) CFA fit, (2) measurement invariance fit indices (metric, scalar), and (3) discriminant validity indices using the nested model comparison, \(\chi^2\) (merge), see Rönkkö & Cho (2022).
fit_metric <- lavaan::cfa(model = EDR, data = tmp, estimator ="MLR",group ="VigNro",group.equal ="loadings") # includes lavaan::lavTestLRT(fit, fit_metric)
Scaled Chi-Squared Difference Test (method = "satorra.bentler.2001")
lavaan->lavTestLRT():
lavaan NOTE: The "Chisq" column contains standard test statistics, not the
robust test that should be reported per model. A robust difference test is
a function of two standard (not robust) statistics.
Df AIC BIC Chisq Chisq diff Df diff Pr(>Chisq)
fit 72 6173.2 6506.5 89.972
fit_metric 84 6160.6 6449.5 101.450 12.23 12 0.4274
fit_scalar <- lavaan::cfa(model = EDR, data = tmp, estimator ="MLR",group ="VigNro",group.equal =c("loadings","intercepts")) # includes lavaan::lavTestLRT(fit, fit_scalar)
Scaled Chi-Squared Difference Test (method = "satorra.bentler.2001")
lavaan->lavTestLRT():
lavaan NOTE: The "Chisq" column contains standard test statistics, not the
robust test that should be reported per model. A robust difference test is
a function of two standard (not robust) statistics.
Df AIC BIC Chisq Chisq diff Df diff Pr(>Chisq)
fit 72 6173.2 6506.5 89.972
fit_scalar 96 6166.3 6410.8 131.114 43.026 24 0.00988 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
fit <- lavaan::cfa(model = EDR, data = tmp, meanstructure =TRUE, estimator ="MLR") # ignores vignettesds1 <- semTools::discriminantValidity(object = fit,merge =TRUE)print(ds1)
lhs op rhs est ci.lower ci.upper Df AIC BIC Chisq
1 E ~~ D 0.3298993 0.1659311 0.4938675 26 6530.387 6634.093 251.9935
2 E ~~ R 0.4930086 0.3434819 0.6425353 26 6567.196 6670.902 288.8023
3 D ~~ R 0.5927842 0.4666596 0.7189089 26 6432.905 6536.611 154.5115
Chisq diff Df diff Pr(>Chisq)
1 110.55156 2 9.863552e-25
2 189.93508 2 5.703224e-42
3 70.95355 2 3.914113e-16
CFA - FM
set.seed(3110)FM <-"M =~ M1 + M3 + M9 # Note one Focus item (F1) has been renamed to motivation item (M9) earlierF =~ F3 + F6 + F10"items <-c("M1","M3","M9","F3","F6","F10")tmp <-CFA_items(df_f,vignette =c(4,5), items=items,vignette_name="FM",bonusquestion="HAAS")fit <- lavaan::cfa(model = FM, data = tmp, estimator ="MLR",group ="VigNro") # includes vignettes explicitlyfit.meas2 <- lavaan::fitMeasures(fit,c("tli","chisq","df","pvalue", "cfi", "rmsea", "srmr"))print(fit.meas2)
fit_metric <- lavaan::cfa(model = FM, data = tmp, estimator ="MLR",group ="VigNro",group.equal ="loadings") # includes lavaan::lavTestLRT(fit, fit_metric)
Scaled Chi-Squared Difference Test (method = "satorra.bentler.2001")
lavaan->lavTestLRT():
lavaan NOTE: The "Chisq" column contains standard test statistics, not the
robust test that should be reported per model. A robust difference test is
a function of two standard (not robust) statistics.
Df AIC BIC Chisq Chisq diff Df diff Pr(>Chisq)
fit 16 2736.5 2861.8 20.560
fit_metric 20 2736.8 2848.9 28.837 9.2554 4 0.05502 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
fit_scalar <- lavaan::cfa(model = FM, data = tmp, estimator ="MLR",group ="VigNro",group.equal =c("loadings","intercepts")) # includes lavaan::lavTestLRT(fit, fit_scalar)
Scaled Chi-Squared Difference Test (method = "satorra.bentler.2001")
lavaan->lavTestLRT():
lavaan NOTE: The "Chisq" column contains standard test statistics, not the
robust test that should be reported per model. A robust difference test is
a function of two standard (not robust) statistics.
Df AIC BIC Chisq Chisq diff Df diff Pr(>Chisq)
fit 16 2736.5 2861.8 20.560
fit_scalar 24 2749.3 2848.2 49.348 32.373 8 7.986e-05 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
fit <- lavaan::cfa(model = FM, data = tmp, meanstructure =TRUE, estimator ="MLR")ds2 <- semTools::discriminantValidity(object = fit,merge =TRUE)
Some of the latent variable variances are estimated instead of fixed to 1. The model is re-estimated by scaling the latent variables by fixing their variances and freeing all factor loadings.
print(ds2)
lhs op rhs est ci.lower ci.upper Df AIC BIC Chisq
1 M ~~ F 0.6965692 0.5431581 0.8499804 9 2860.671 2920.041 108.1906
Chisq diff Df diff Pr(>Chisq)
1 26.36594 1 2.82478e-07
fit_metric <- lavaan::cfa(model = CB, data = tmp, estimator ="MLR",group ="VigNro",group.equal ="loadings") # includes lavaan::lavTestLRT(fit, fit_metric)
Scaled Chi-Squared Difference Test (method = "satorra.bentler.2001")
lavaan->lavTestLRT():
lavaan NOTE: The "Chisq" column contains standard test statistics, not the
robust test that should be reported per model. A robust difference test is
a function of two standard (not robust) statistics.
Df AIC BIC Chisq Chisq diff Df diff Pr(>Chisq)
fit 16 3010.1 3135.5 22.136
fit_metric 20 3003.5 3115.6 23.470 0.87442 4 0.9282
fit_scalar <- lavaan::cfa(model = CB, data = tmp, estimator ="MLR",group ="VigNro",group.equal =c("loadings","intercepts")) # includes lavaan::lavTestLRT(fit, fit_scalar)
Scaled Chi-Squared Difference Test (method = "satorra.bentler.2001")
lavaan->lavTestLRT():
lavaan NOTE: The "Chisq" column contains standard test statistics, not the
robust test that should be reported per model. A robust difference test is
a function of two standard (not robust) statistics.
Df AIC BIC Chisq Chisq diff Df diff Pr(>Chisq)
fit 16 3010.1 3135.5 22.136
fit_scalar 24 3007.8 3106.8 35.803 9.6653 8 0.2893
fit <- lavaan::cfa(model = CB, data = tmp, meanstructure =TRUE,estimator ="MLR")ds3 <- semTools::discriminantValidity(object = fit,merge =TRUE)
Some of the latent variable variances are estimated instead of fixed to 1. The model is re-estimated by scaling the latent variables by fixing their variances and freeing all factor loadings.
print(ds3)
lhs op rhs est ci.lower ci.upper Df AIC BIC Chisq
1 G ~~ L 0.1375599 -0.04939183 0.3245117 9 3422.591 3481.96 285.0453
Chisq diff Df diff Pr(>Chisq)
1 95.37567 1 1.574775e-22
fit_metric <- lavaan::cfa(model = PEP, data = tmp, estimator ="MLR",group ="VigNro",group.equal ="loadings") # includes lavaan::lavTestLRT(fit, fit_metric)
Scaled Chi-Squared Difference Test (method = "satorra.bentler.2001")
lavaan->lavTestLRT():
lavaan NOTE: The "Chisq" column contains standard test statistics, not the
robust test that should be reported per model. A robust difference test is
a function of two standard (not robust) statistics.
Df AIC BIC Chisq Chisq diff Df diff Pr(>Chisq)
fit 16 2599.2 2724.6 26.872
fit_metric 20 2596.3 2708.4 31.893 5.0871 4 0.2785
fit_scalar <- lavaan::cfa(model = PEP, data = tmp, estimator ="MLR",group ="VigNro",group.equal =c("loadings","intercepts")) # includes lavaan::lavTestLRT(fit, fit_scalar)
Scaled Chi-Squared Difference Test (method = "satorra.bentler.2001")
lavaan->lavTestLRT():
lavaan NOTE: The "Chisq" column contains standard test statistics, not the
robust test that should be reported per model. A robust difference test is
a function of two standard (not robust) statistics.
Df AIC BIC Chisq Chisq diff Df diff Pr(>Chisq)
fit 16 2599.2 2724.6 26.872
fit_scalar 24 2606.4 2705.4 50.063 23.642 8 0.002631 **
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
fit <- lavaan::cfa(model = PEP, data = tmp, meanstructure =TRUE,estimator ="MLR")ds4 <- semTools::discriminantValidity(object = fit,merge =TRUE)
Some of the latent variable variances are estimated instead of fixed to 1. The model is re-estimated by scaling the latent variables by fixing their variances and freeing all factor loadings.
print(ds4)
lhs op rhs est ci.lower ci.upper Df AIC BIC Chisq
1 C ~~ X 0.7450395 0.6089644 0.8811147 9 2661.872 2721.241 61.70509
Chisq diff Df diff Pr(>Chisq)
1 21.0942 1 4.372504e-06
fit_metric <- lavaan::cfa(model = AIA, data = tmp, estimator ="MLR",group ="VigNro",group.equal ="loadings") # includes # Compare configural and metric modelslavaan::lavTestLRT(fit, fit_metric)
Scaled Chi-Squared Difference Test (method = "satorra.bentler.2001")
lavaan->lavTestLRT():
lavaan NOTE: The "Chisq" column contains standard test statistics, not the
robust test that should be reported per model. A robust difference test is
a function of two standard (not robust) statistics.
Df AIC BIC Chisq Chisq diff Df diff Pr(>Chisq)
fit 51 5474.9 5775.0 80.103
fit_metric 61 5478.3 5741.2 103.419 20.864 10 0.02206 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
fit_scalar <- lavaan::cfa(model = AIA, data = tmp, estimator ="MLR",group ="VigNro",group.equal =c("loadings","intercepts")) # includes lavaan::lavTestLRT(fit, fit_scalar)
Scaled Chi-Squared Difference Test (method = "satorra.bentler.2001")
lavaan->lavTestLRT():
lavaan NOTE: The "Chisq" column contains standard test statistics, not the
robust test that should be reported per model. A robust difference test is
a function of two standard (not robust) statistics.
Df AIC BIC Chisq Chisq diff Df diff Pr(>Chisq)
fit 51 5474.9 5775.0 80.103
fit_scalar 71 5508.9 5734.8 154.014 70.572 20 1.468e-07 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
fit <- lavaan::cfa(model = AIA, data = tmp, meanstructure =TRUE,estimator ="MLR")ds5 <- semTools::discriminantValidity(object = fit,merge =TRUE)
Some of the latent variable variances are estimated instead of fixed to 1. The model is re-estimated by scaling the latent variables by fixing their variances and freeing all factor loadings.
print(ds5)
lhs op rhs est ci.lower ci.upper Df AIC BIC Chisq
1 S ~~ I 0.3237894 0.1860453 0.4615334 19 5758.875 5851.469 260.1804
2 S ~~ B 0.5959015 0.4929085 0.6988944 19 5690.502 5783.097 191.8079
3 I ~~ B 0.5215025 0.3852497 0.6577553 19 5685.831 5778.425 187.1364
Chisq diff Df diff Pr(>Chisq)
1 92.37423 2 8.733540e-21
2 70.24699 2 5.572616e-16
3 62.68685 2 2.441870e-14
Tables for manuscript
cfa_table<-rbind(fit.meas1,fit.meas2,fit.meas3,fit.meas4,fit.meas5)rownames(cfa_table)<-c("EDR","FM","CB","PEP","AIA")print(knitr::kable(cfa_table,digits=3, caption="CFA measures of fit for constructs."))
CFA measures of fit for constructs.
tli
chisq
df
pvalue
cfi
rmsea
srmr
EDR
0.976
89.972
72
0.074
0.984
0.050
0.046
FM
0.987
20.560
16
0.196
0.993
0.053
0.033
CB
0.978
22.136
16
0.139
0.989
0.062
0.040
PEP
0.961
26.872
16
0.043
0.979
0.082
0.040
AIA
0.961
80.104
51
0.006
0.976
0.076
0.046
DS1<-as.data.frame(ds1)DS2<-as.data.frame(ds2)DS3<-as.data.frame(ds3)DS4<-as.data.frame(ds4)DS5<-as.data.frame(ds5)DS1$Episode<-"EDR"DS2$Episode<-"FM"DS3$Episode<-"CB"DS4$Episode<-"PEP"DS5$Episode<-"AIA"ds_table<-rbind(DS1,DS2,DS3,DS4,DS5)ds_table <- ds_table %>%relocate(Episode)print(knitr::kable(ds_table,digits=3, caption="Discriminant validity indices for sub-constructs."))
Discriminant validity indices for sub-constructs.
Episode
lhs
op
rhs
est
ci.lower
ci.upper
Df
AIC
BIC
Chisq
Chisq diff
Df diff
Pr(>Chisq)
EDR
E
~~
D
0.330
0.166
0.494
26
6530.387
6634.093
251.993
110.552
2
0
EDR
E
~~
R
0.493
0.343
0.643
26
6567.196
6670.902
288.802
189.935
2
0
EDR
D
~~
R
0.593
0.467
0.719
26
6432.905
6536.611
154.511
70.954
2
0
FM
M
~~
F
0.697
0.543
0.850
9
2860.671
2920.041
108.191
26.366
1
0
CB
G
~~
L
0.138
-0.049
0.325
9
3422.591
3481.960
285.045
95.376
1
0
PEP
C
~~
X
0.745
0.609
0.881
9
2661.872
2721.241
61.705
21.094
1
0
AIA
S
~~
I
0.324
0.186
0.462
19
5758.875
5851.469
260.180
92.374
2
0
AIA
S
~~
B
0.596
0.493
0.699
19
5690.502
5783.097
191.808
70.247
2
0
AIA
I
~~
B
0.522
0.385
0.658
19
5685.831
5778.425
187.136
62.687
2
0
Descriptives: means and SDs and loadings
Show summary of values (means, SDs, and scores) across the sub-constructs when run jointly (CFA model with 12 sub-constructs).
For figures in the manuscript.
lhs
op
rhs
est
se
z
pvalue
ci.lower
ci.upper
std.lv
std.all
E
~~
D
0.23
0.03
8.12
0
0.17
0.28
0.43
0.43
E
~~
R
0.26
0.03
9.98
0
0.21
0.31
0.52
0.52
D
~~
R
0.36
0.03
11.63
0
0.30
0.42
0.56
0.56
lhs
op
rhs
est
se
z
pvalue
ci.lower
ci.upper
std.lv
std.all
M
~~
F
0.58
0.03
16.6
0
0.51
0.65
0.72
0.72
lhs
op
rhs
est
se
z
pvalue
ci.lower
ci.upper
std.lv
std.all
G
~~
L
0.56
0.04
13.35
0
0.48
0.65
0.5
0.5
lhs
op
rhs
est
se
z
pvalue
ci.lower
ci.upper
std.lv
std.all
C
~~
X
0.43
0.03
13.41
0
0.37
0.49
0.67
0.67
CFA Standardized Loadings with Descriptives.
lhs
op
rhs
est
se
z
pvalue
ci.lower
ci.upper
std.lv
std.all
S
~~
I
0.58
0.04
14.56
0
0.50
0.65
0.52
0.52
S
~~
B
0.43
0.03
15.12
0
0.38
0.49
0.61
0.61
I
~~
B
0.46
0.03
16.06
0
0.41
0.52
0.75
0.75
Item
Factor
Standardized
Mean
SD
E1
E
0.75
4.16
0.85
E6
E
0.80
4.29
0.77
E7
E
0.75
3.95
0.97
D1
D
0.75
3.61
1.09
D6
D
0.74
3.30
1.24
D7
D
0.68
3.25
1.18
R1
R
0.84
3.89
0.92
R3
R
0.84
3.83
1.01
R7
R
0.86
3.94
0.93
M1
M
0.81
3.51
1.05
M3
M
0.63
3.62
1.07
M9
M
0.86
3.63
1.06
F10
F
0.80
3.64
1.02
F3
F
0.87
3.38
1.07
F6
F
0.85
3.46
1.08
G2
G
0.89
3.04
1.21
G5
G
0.88
3.17
1.22
G9
G
0.83
3.00
1.24
L5
L
0.90
3.43
1.17
L6
L
0.90
3.39
1.15
L7
L
0.69
3.55
1.09
C2
C
0.76
3.85
0.97
C3
C
0.83
4.05
0.91
C6
C
0.80
4.06
0.87
X4
X
0.81
3.38
1.07
X8
X
0.84
3.53
1.07
X9
X
0.71
3.42
1.13
S3
S
0.90
2.74
1.25
S4
S
0.89
2.66
1.23
S5
S
0.88
2.81
1.28
I1
I
0.89
3.32
1.11
I2
I
0.87
3.27
1.13
B1
B
0.74
4.15
0.85
B3
B
0.77
3.83
1.03
B5
B
0.71
3.11
1.18
Save factor scores for emotion comparison
# Extract factor scores from your CFAfactor_scores <-lavPredict(fit1, method ="regression")fs <-as.data.frame(factor_scores)fs$VigNro <- tmp$VigNrofs$ProlificID <- tmp$ProlificIDdim(tmp)
[1] 1200 37
# Save to CSV#write.csv(fs, file = here("data","factor_scores_exp2.csv"), row.names = TRUE)