Structural Equation Modeling | Lab Session 3
1 What we are going to cover
- Measurement equivalence
- Multi-group SEM
2 Data
The data set used throughout is the European Social Survey ESS4-2008 Edition 4.5 was released on 1 December 2018. We will restrict the analysis to the Belgian case. Each line in the data set represents a Belgian respondent. The full dataset an documentation can be found on the ESS website
Codebook:
gvslvol Standard of living for the old, governments’ responsibility (0 Not governments’ responsibility at all - 10 Entirely governments’ responsibility)
gvslvue Standard of living for the unemployed, governments’ responsibility (0 Not governments’ responsibility at all - 10 Entirely governments’ responsibility)
gvhlthc Health care for the sick, governments’ responsibility (0 Not governments’ responsibility at all - 10 Entirely governments’ responsibility)
gvcldcr Child care services for working parents, governments’ responsibility (0 Not governments’ responsibility at all - 10 Entirely governments’ responsibility)
gvjbevn Job for everyone, governments’ responsibility (0 Not governments’ responsibility at all - 10 Entirely governments’ responsibility)
gvpdlwk Paid leave from work to care for sick family, governments’ responsibility (0 Not governments’ responsibility at all - 10 Entirely governments’ responsibility)
sbstrec Social benefits/services place too great strain on economy (1 Agree strongly - 5 Disagree strongly)
sbbsntx Social benefits/services cost businesses too much in taxes/charges (1 Agree strongly - 5 Disagree strongly)
sbprvpv Social benefits/services prevent widespread poverty (1 Agree strongly - 5 Disagree strongly)
sbeqsoc Social benefits/services lead to a more equal society (1 Agree strongly - 5 Disagree strongly)
sbcwkfm Social benefits/services make it easier to combine work and family (1 Agree strongly - 5 Disagree strongly)
sblazy Social benefits/services make people lazy (1 Agree strongly - 5 Disagree strongly)
sblwcoa Social benefits/services make people less willing care for one another (1 Agree strongly - 5 Disagree strongly)
sblwlka Social benefits/services make people less willing look after themselves/family (1 Agree strongly - 5 Disagree strongly)
In addition, we will use some other variables
agea Respondent’s age
eduyrs Years of full-time education completed
gndr Gender (1 Male, 2 Female)
hinctnta Household’s total net income, all sources (Deciles of the actual household income range in Belgium)
gincdif Government should reduce differences in income levels (1 Agree strongly - 5 Disagree strongly)
dfincac Large differences in income acceptable to reward talents and efforts (1 Agree strongly - 5 Disagree strongly)
smdfslv For fair society, differences in standard of living should be small (1 Agree strongly - 5 Disagree strongly)
3 Environment preparation
First, let’s load the necessary packages to load, manipulate, visualize and analyse the data.
# Uncomment this once if you need to install the packages on your system
### DATA MANIPULATION ###
# install.packages("haven") # data import from spss
# install.packages("dplyr") # data manipulation
# install.packages("psych") # descriptives
# install.packages("stringr") # string manipulation
# install.packages("purrr") # table manipulation
### MODELING ###
# install.packages("lavaan") # SEM modelling
### VISUALIZATION ###
# install.packages("tidySEM") # plotting SEM models
# Load the packages
### DATA MANIPULATION ###
library("haven")
library("dplyr")
library("psych")
library('stringr')
### MODELING ###
library("lavaan")
### VISUALIZATION ###
library("tidySEM")
library("purrr") 4 Measurament Equivalence
In a nutshell:
- Measurement Equivalence is achieved if a measurement instrument produces equivalent results, regardless of some unrelated properties of the test subjects.
- The absence of measurement equivalence would imply some degree of distortion of the results (“bias”). For instance, an IQ test that favors males by including “gender-biased test items”
- What is the effect of such a bias? Would we draw incorrect conclusions?
- Can we detect whether such bias is present in our data?
Pros:
- Simultaneously testing a CFA model in multiple groups
- Possible to test equality of each parameter
Cons:
- Requires a larger sample size than the MIMIC approach
- Requires more programming
- Not very elegant for comparing many groups
There are different types of measurement equivalence:
- Configural Invariance
- same factor structure in each group
- free item loadings
- free intercepts
- free residuals
- No latent mean difference is estimated (fixed to 0)
- Metric Invariance (also called “weak” invariance)
- item loadings (set to equal)
- free intercepts
- free residuals
- No latent mean difference is estimated (fixed to 0)
- Scalar Invariance (also called “strong” invariance)
- item loadings (set to equal)
- intercepts (set to equal)
- free residuals
- latent mean difference is estimated (fixed to 0 in G1)
- Strict Invariance
- item loadings (set to equal)
- intercepts (set to equal)
- residuals (fixed to 1)
- latent mean difference is estimated (fixed to 0 in G1)
- Structural Invariance
- item loadings (set to equal)
- intercepts (set to equal)
- residuals (fixed to 1)
- latent mean difference is estimated (fixed to 0 in G1)
- factor variances (set to equal)
- factor covariances (set to equal if the factors are more than 1)
In this lab, we are going to test whether welfare support is measured equivalently across genders. We start with configural invariance
4.1 Configural invariance
ess_df <- haven::read_sav("https://github.com/albertostefanelli/SEM_labs/raw/master/data/ESS4_belgium.sav")
# lavaan requires the grouping variable to be a factor
# gender is coded as 1 Male, 2 Female
ess_df$gndr <- factor(ess_df$gndr,
levels = c("1", "2"), # levels
labels = c("Male", "Female")) # labels
model_ws <-'welf_supp =~ gvslvol + gvslvue + gvhlthc'
fit_configural <- cfa(model_ws,
data = ess_df,
group = "gndr")
summary(fit_configural, standardized=TRUE)lavaan 0.6-10 ended normally after 51 iterations
Estimator ML
Optimization method NLMINB
Number of model parameters 18
Number of observations per group: Used Total
Male 861 864
Female 891 896
Model Test User Model:
Test statistic 0.000
Degrees of freedom 0
Test statistic for each group:
Male 0.000
Female 0.000
Parameter Estimates:
Standard errors Standard
Information Expected
Information saturated (h1) model Structured
Group 1 [Male]:
Latent Variables:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
welf_supp =~
gvslvol 1.000 1.215 0.802
gvslvue 0.657 0.067 9.780 0.000 0.799 0.399
gvhlthc 0.974 0.087 11.183 0.000 1.183 0.805
Intercepts:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
.gvslvol 7.851 0.052 152.113 0.000 7.851 5.184
.gvslvue 6.005 0.068 88.020 0.000 6.005 3.000
.gvhlthc 8.057 0.050 160.839 0.000 8.057 5.481
welf_supp 0.000 0.000 0.000
Variances:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
.gvslvol 0.817 0.130 6.271 0.000 0.817 0.356
.gvslvue 3.369 0.171 19.702 0.000 3.369 0.841
.gvhlthc 0.760 0.123 6.161 0.000 0.760 0.352
welf_supp 1.477 0.162 9.142 0.000 1.000 1.000
Group 2 [Female]:
Latent Variables:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
welf_supp =~
gvslvol 1.000 1.259 0.869
gvslvue 0.523 0.068 7.706 0.000 0.659 0.360
gvhlthc 0.801 0.090 8.906 0.000 1.008 0.685
Intercepts:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
.gvslvol 7.892 0.049 162.542 0.000 7.892 5.445
.gvslvue 6.112 0.061 99.650 0.000 6.112 3.338
.gvhlthc 8.007 0.049 162.306 0.000 8.007 5.437
welf_supp 0.000 0.000 0.000
Variances:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
.gvslvol 0.514 0.171 3.008 0.003 0.514 0.245
.gvslvue 2.918 0.146 20.011 0.000 2.918 0.870
.gvhlthc 1.151 0.121 9.480 0.000 1.151 0.531
welf_supp 1.586 0.195 8.141 0.000 1.000 1.000All parameters are different across groups
4.2 Metric Invariance
For metric invariance we are going to constrain factor loadings equal across groups. This shows that the construct has the same meaning across groups. Metric invariance is necessary for correlations and regressions.
Metric non-invariance:
- Meaning of the items are different across groups
- Extreme response style might be present for some items
- Or some people are more likely to choose a middle point for some itmes.
# we can simply tell lavaan what should be constrained across groups
fit_metric <- cfa(model_ws,
data = ess_df,
group = "gndr",
group.equal = c("loadings")
)
summary(fit_metric, standardized=TRUE)lavaan 0.6-10 ended normally after 45 iterations
Estimator ML
Optimization method NLMINB
Number of model parameters 18
Number of equality constraints 2
Number of observations per group: Used Total
Male 861 864
Female 891 896
Model Test User Model:
Test statistic 2.555
Degrees of freedom 2
P-value (Chi-square) 0.279
Test statistic for each group:
Male 1.155
Female 1.400
Parameter Estimates:
Standard errors Standard
Information Expected
Information saturated (h1) model Structured
Group 1 [Male]:
Latent Variables:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
welf_supp =~
gvslvol 1.000 1.268 0.834
gvslvue (.p2.) 0.595 0.048 12.472 0.000 0.754 0.379
gvhlthc (.p3.) 0.897 0.062 14.381 0.000 1.138 0.776
Intercepts:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
.gvslvol 7.851 0.052 151.451 0.000 7.851 5.161
.gvslvue 6.005 0.068 88.528 0.000 6.005 3.017
.gvhlthc 8.057 0.050 161.316 0.000 8.057 5.498
welf_supp 0.000 0.000 0.000
Variances:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
.gvslvol 0.706 0.118 5.964 0.000 0.706 0.305
.gvslvue 3.392 0.170 19.924 0.000 3.392 0.856
.gvhlthc 0.853 0.100 8.522 0.000 0.853 0.397
welf_supp 1.608 0.148 10.875 0.000 1.000 1.000
Group 2 [Female]:
Latent Variables:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
welf_supp =~
gvslvol 1.000 1.188 0.823
gvslvue (.p2.) 0.595 0.048 12.472 0.000 0.707 0.384
gvhlthc (.p3.) 0.897 0.062 14.381 0.000 1.066 0.721
Intercepts:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
.gvslvol 7.892 0.048 163.261 0.000 7.892 5.469
.gvslvue 6.112 0.062 99.094 0.000 6.112 3.320
.gvhlthc 8.007 0.050 161.742 0.000 8.007 5.419
welf_supp 0.000 0.000 0.000
Variances:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
.gvslvol 0.670 0.107 6.243 0.000 0.670 0.322
.gvslvue 2.890 0.144 20.073 0.000 2.890 0.853
.gvhlthc 1.047 0.097 10.821 0.000 1.047 0.480
welf_supp 1.412 0.131 10.815 0.000 1.000 1.000Q: What (.p2.) and (.p3.) in the output refer to ?
4.3 Scalar Invariance
Scalar invariance requires item intercepts and factor loadings equal across groups. This is important for assessing mean difference of the latent variable across groups.
Scalar non-invariance:
- A group tend to systematically give higher or lower item response. For instance, female respondents might express higher support for child care services.
- It is an additive effect. It affects the means of the observed item, hence affects the mean of the scale and the latent variable.
fit_scalar <- cfa(model_ws,
data = ess_df,
group = "gndr",
group.equal = c("loadings",
"intercepts")
)
summary(fit_scalar, standardized=TRUE)lavaan 0.6-10 ended normally after 50 iterations
Estimator ML
Optimization method NLMINB
Number of model parameters 19
Number of equality constraints 5
Number of observations per group: Used Total
Male 861 864
Female 891 896
Model Test User Model:
Test statistic 6.199
Degrees of freedom 4
P-value (Chi-square) 0.185
Test statistic for each group:
Male 3.053
Female 3.145
Parameter Estimates:
Standard errors Standard
Information Expected
Information saturated (h1) model Structured
Group 1 [Male]:
Latent Variables:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
welf_supp =~
gvslvol 1.000 1.270 0.835
gvslvue (.p2.) 0.594 0.048 12.441 0.000 0.754 0.379
gvhlthc (.p3.) 0.894 0.062 14.356 0.000 1.135 0.775
Intercepts:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
.gvslvol (.p8.) 7.865 0.050 156.652 0.000 7.865 5.171
.gvslvue (.p9.) 6.059 0.050 120.178 0.000 6.059 3.043
.gvhlthc (.10.) 8.028 0.047 171.303 0.000 8.028 5.478
wlf_spp 0.000 0.000 0.000
Variances:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
.gvslvol 0.700 0.119 5.893 0.000 0.700 0.303
.gvslvue 3.396 0.170 19.924 0.000 3.396 0.856
.gvhlthc 0.859 0.100 8.594 0.000 0.859 0.400
welf_supp 1.614 0.148 10.873 0.000 1.000 1.000
Group 2 [Female]:
Latent Variables:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
welf_supp =~
gvslvol 1.000 1.190 0.825
gvslvue (.p2.) 0.594 0.048 12.441 0.000 0.707 0.384
gvhlthc (.p3.) 0.894 0.062 14.356 0.000 1.064 0.720
Intercepts:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
.gvslvol (.p8.) 7.865 0.050 156.652 0.000 7.865 5.450
.gvslvue (.p9.) 6.059 0.050 120.178 0.000 6.059 3.290
.gvhlthc (.10.) 8.028 0.047 171.303 0.000 8.028 5.432
wlf_spp 0.014 0.066 0.207 0.836 0.012 0.012
Variances:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
.gvslvol 0.666 0.108 6.173 0.000 0.666 0.320
.gvslvue 2.893 0.144 20.074 0.000 2.893 0.853
.gvhlthc 1.053 0.097 10.888 0.000 1.053 0.482
welf_supp 1.417 0.131 10.812 0.000 1.000 1.0004.4 Strict invariance
For achieving strict invariance we need to constrain item residual variance, factor loadings, and intercepts equal across groups. If we want to use sum scores of observed items, our construct should reach strict invariance. This is because observed variance is a combination of true score variance and residual variance.
fit_strict <- cfa(model_ws,
data = ess_df,
group = "gndr",
group.equal = c("loadings",
"intercepts",
"residuals")
)
summary(fit_strict, standardized=TRUE)lavaan 0.6-10 ended normally after 40 iterations
Estimator ML
Optimization method NLMINB
Number of model parameters 19
Number of equality constraints 8
Number of observations per group: Used Total
Male 861 864
Female 891 896
Model Test User Model:
Test statistic 15.622
Degrees of freedom 7
P-value (Chi-square) 0.029
Test statistic for each group:
Male 8.016
Female 7.606
Parameter Estimates:
Standard errors Standard
Information Expected
Information saturated (h1) model Structured
Group 1 [Male]:
Latent Variables:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
welf_supp =~
gvslvol 1.000 1.273 0.841
gvslvue (.p2.) 0.591 0.048 12.408 0.000 0.753 0.391
gvhlthc (.p3.) 0.886 0.062 14.391 0.000 1.127 0.753
Intercepts:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
.gvslvol (.p8.) 7.865 0.050 157.002 0.000 7.865 5.197
.gvslvue (.p9.) 6.055 0.050 120.189 0.000 6.055 3.144
.gvhlthc (.10.) 8.025 0.047 170.293 0.000 8.025 5.363
wlf_spp 0.000 0.000 0.000
Variances:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
.gvslvol (.p4.) 0.669 0.103 6.508 0.000 0.669 0.292
.gvslvue (.p5.) 3.143 0.112 28.083 0.000 3.143 0.847
.gvhlthc (.p6.) 0.968 0.085 11.361 0.000 0.968 0.432
wlf_spp 1.621 0.146 11.085 0.000 1.000 1.000
Group 2 [Female]:
Latent Variables:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
welf_supp =~
gvslvol 1.000 1.198 0.826
gvslvue (.p2.) 0.591 0.048 12.408 0.000 0.709 0.371
gvhlthc (.p3.) 0.886 0.062 14.391 0.000 1.061 0.733
Intercepts:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
.gvslvol (.p8.) 7.865 0.050 157.002 0.000 7.865 5.421
.gvslvue (.p9.) 6.055 0.050 120.189 0.000 6.055 3.172
.gvhlthc (.10.) 8.025 0.047 170.293 0.000 8.025 5.546
wlf_spp 0.015 0.067 0.220 0.825 0.012 0.012
Variances:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
.gvslvol (.p4.) 0.669 0.103 6.508 0.000 0.669 0.318
.gvslvue (.p5.) 3.143 0.112 28.083 0.000 3.143 0.862
.gvhlthc (.p6.) 0.968 0.085 11.361 0.000 0.968 0.462
wlf_spp 1.435 0.132 10.848 0.000 1.000 1.0004.5 Structural invariance
For achieving strict invariance we need to constrain factor variances, factor covariances (if more than one latent factors), item residual variance, factor loadings, and intercepts equal across groups. Additionally, in multiple group SEM analysis we should also constrain regression path coefficients.
fit_structural <- cfa(model_ws,
data = ess_df,
group = "gndr",
group.equal = c("loadings",
"intercepts",
"residuals",
"lv.variances",
"lv.covariances")
)
summary(fit_structural, standardized=TRUE)lavaan 0.6-10 ended normally after 37 iterations
Estimator ML
Optimization method NLMINB
Number of model parameters 19
Number of equality constraints 9
Number of observations per group: Used Total
Male 861 864
Female 891 896
Model Test User Model:
Test statistic 17.621
Degrees of freedom 8
P-value (Chi-square) 0.024
Test statistic for each group:
Male 8.883
Female 8.738
Parameter Estimates:
Standard errors Standard
Information Expected
Information saturated (h1) model Structured
Group 1 [Male]:
Latent Variables:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
welf_supp =~
gvslvol 1.000 1.233 0.832
gvslvue (.p2.) 0.592 0.048 12.389 0.000 0.730 0.381
gvhlthc (.p3.) 0.889 0.062 14.289 0.000 1.096 0.745
Intercepts:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
.gvslvol (.p8.) 7.865 0.049 160.654 0.000 7.865 5.307
.gvslvue (.p9.) 6.055 0.050 121.111 0.000 6.055 3.158
.gvhlthc (.10.) 8.025 0.046 173.510 0.000 8.025 5.454
wlf_spp 0.000 0.000 0.000
Variances:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
.gvslvol (.p4.) 0.675 0.103 6.530 0.000 0.675 0.307
.gvslvue (.p5.) 3.143 0.112 28.083 0.000 3.143 0.855
.gvhlthc (.p6.) 0.963 0.086 11.196 0.000 0.963 0.445
wlf_spp (.p7.) 1.521 0.123 12.361 0.000 1.000 1.000
Group 2 [Female]:
Latent Variables:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
welf_supp =~
gvslvol 1.000 1.233 0.832
gvslvue (.p2.) 0.592 0.048 12.389 0.000 0.730 0.381
gvhlthc (.p3.) 0.889 0.062 14.289 0.000 1.096 0.745
Intercepts:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
.gvslvol (.p8.) 7.865 0.049 160.654 0.000 7.865 5.307
.gvslvue (.p9.) 6.055 0.050 121.111 0.000 6.055 3.158
.gvhlthc (.10.) 8.025 0.046 173.510 0.000 8.025 5.454
wlf_spp 0.014 0.066 0.216 0.829 0.012 0.012
Variances:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
.gvslvol (.p4.) 0.675 0.103 6.530 0.000 0.675 0.307
.gvslvue (.p5.) 3.143 0.112 28.083 0.000 3.143 0.855
.gvhlthc (.p6.) 0.963 0.086 11.196 0.000 0.963 0.445
wlf_spp (.p7.) 1.521 0.123 12.361 0.000 1.000 1.0004.6 Evaluating measurement invariance
As seen before there are two instruments to evaluate SEM models. Global and local fit measures.
Global Fit:
- Substantial decrease in goodness of fit indicates non-invariance
- It is a good practise to look at several model fit indices rather than relying on a single one
- Since all of the measurement invariance models are nested within one another, we can do a chi-square difference test.
# Let's create quick function to extract the fit indices
model_fit <- function(lavobject) {
vars <- c("df", "cfi", "tli", "rmsea", "rmsea.ci.lower", "rmsea.ci.upper", "rmsea.pvalue", "srmr")
return(fitmeasures(lavobject)[vars] %>% data.frame() %>% round(2) %>% t())
}
table_fit <-
list(model_fit(fit_configural),
model_fit(fit_metric),
model_fit(fit_scalar),
model_fit(fit_strict),
model_fit(fit_structural)) %>%
reduce(rbind)
rownames(table_fit) <- c("Configural", "Metric", "Scalar","Strict","Structural")
table_fit df cfi tli rmsea rmsea.ci.lower rmsea.ci.upper rmsea.pvalue srmr
Configural 0 1.00 1.00 0.00 0.00 0.00 NA 0.00
Metric 2 1.00 1.00 0.02 0.00 0.07 0.78 0.01
Scalar 4 1.00 1.00 0.03 0.00 0.06 0.85 0.02
Strict 7 0.99 0.99 0.04 0.01 0.06 0.77 0.03
Structural 8 0.99 0.99 0.04 0.01 0.06 0.80 0.04# Let's compare the nested model using the anova function
table_anova <- list(anova(fit_configural, fit_metric),
anova(fit_metric, fit_scalar),
anova(fit_scalar, fit_strict),
anova(fit_strict, fit_structural)) %>%
reduce(rbind) %>%
.[-c(3, 5, 7),]
table_anovaChi-Squared Difference Test
Df AIC BIC Chisq Chisq diff Df diff Pr(>Chisq)
fit_configural 0 18887 18985 0.0000
fit_metric 2 18885 18973 2.5551 2.5551 2 0.27872
fit_scalar 4 18885 18962 6.1987 3.6436 2 0.16173
fit_strict 7 18888 18949 15.6221 9.4234 3 0.02416 *
fit_structural 8 18888 18943 17.6210 1.9990 1 0.15741
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1Remember:
- Invariance is achieved if there is a non-significant chi-square change
- Invariance is NOT achieved if there is a significant chi-square change
Q: Is welfare support invariant across genders? What kind of analyses we can perform?
4.7 Partial invariance
Local Fit
- The modification indices (MI) indicates the expected decrease in chi-square if a restricted parameter is to be freed in a less restrictive model
- As done before, we can look for the largest MI value in the lavaan MI output.
- We can free one parameter at a time through an iterative process
- The usual cut-off value is 3.84, but this should be be adjusted based on sample size and number of tests conducted (type I error)
We have seen from the anova table that we reject the null hypothesis and that the scal invariance model fits better than the strict invariance model. What we can further assess whether there is partial stric invariance. We can do this by using the lavTestScore() function in the lavaan function. This function allows us to see the effect of releasing equality constraints across the groups.
lavTestScore(fit_strict)$test
total score test:
test X2 df p.value
1 score 15.555 8 0.049
$uni
univariate score tests:
lhs op rhs X2 df p.value
1 .p2. == .p13. 1.339 1 0.247
2 .p3. == .p14. 0.596 1 0.440
3 .p4. == .p15. 0.066 1 0.797
4 .p5. == .p16. 5.166 1 0.023
5 .p6. == .p17. 3.980 1 0.046
6 .p8. == .p19. 1.180 1 0.277
7 .p9. == .p20. 1.429 1 0.232
8 .p10. == .p21. 2.718 1 0.099The first output is a multivariate score test (i.e. Lagrange multiplier test). It indicates whether freeing all equality constraints represents an improvement in fit over the base model. In this case, we reject the null hypothesis (barely). The second part of the output is a univariate score tests (i.e., the chi-square difference tests) to see which equality constraints should be relaxed. In our case, the largest change in chi-square difference is produced by freeing the .p5. == .p16 and .p6. == .p17.. We can use the function parTable() to see to which parameters .p5. == .p16 and .p6. == .p17. correspond.
parTable(fit_strict) id lhs op rhs user block group free ustart exo label plabel start est se
1 1 welf_supp =~ gvslvol 1 1 1 0 1 0 .p1. 1.000 1.000 0.000
2 2 welf_supp =~ gvslvue 1 1 1 1 NA 0 .p2. .p2. 0.657 0.591 0.048
3 3 welf_supp =~ gvhlthc 1 1 1 2 NA 0 .p3. .p3. 0.974 0.886 0.062
4 4 gvslvol ~~ gvslvol 0 1 1 3 NA 0 .p4. .p4. 1.147 0.669 0.103
5 5 gvslvue ~~ gvslvue 0 1 1 4 NA 0 .p5. .p5. 2.003 3.143 0.112
6 6 gvhlthc ~~ gvhlthc 0 1 1 5 NA 0 .p6. .p6. 1.080 0.968 0.085
7 7 welf_supp ~~ welf_supp 0 1 1 6 NA 0 .p7. 0.050 1.621 0.146
8 8 gvslvol ~1 0 1 1 7 NA 0 .p8. .p8. 7.851 7.865 0.050
9 9 gvslvue ~1 0 1 1 8 NA 0 .p9. .p9. 6.005 6.055 0.050
10 10 gvhlthc ~1 0 1 1 9 NA 0 .p10. .p10. 8.057 8.025 0.047
11 11 welf_supp ~1 0 1 1 0 0 0 .p11. 0.000 0.000 0.000
12 12 welf_supp =~ gvslvol 1 2 2 0 1 0 .p12. 1.000 1.000 0.000
13 13 welf_supp =~ gvslvue 1 2 2 10 NA 0 .p2. .p13. 0.523 0.591 0.048
14 14 welf_supp =~ gvhlthc 1 2 2 11 NA 0 .p3. .p14. 0.801 0.886 0.062
15 15 gvslvol ~~ gvslvol 0 2 2 12 NA 0 .p4. .p15. 1.050 0.669 0.103
16 16 gvslvue ~~ gvslvue 0 2 2 13 NA 0 .p5. .p16. 1.676 3.143 0.112
17 17 gvhlthc ~~ gvhlthc 0 2 2 14 NA 0 .p6. .p17. 1.084 0.968 0.085
18 18 welf_supp ~~ welf_supp 0 2 2 15 NA 0 .p18. 0.050 1.435 0.132
19 19 gvslvol ~1 0 2 2 16 NA 0 .p8. .p19. 7.892 7.865 0.050
20 20 gvslvue ~1 0 2 2 17 NA 0 .p9. .p20. 6.112 6.055 0.050
21 21 gvhlthc ~1 0 2 2 18 NA 0 .p10. .p21. 8.007 8.025 0.047
22 22 welf_supp ~1 0 2 2 19 NA 0 .p22. 0.000 0.015 0.067
23 23 .p2. == .p13. 2 0 0 0 NA 0 0.000 0.000 0.000
24 24 .p3. == .p14. 2 0 0 0 NA 0 0.000 0.000 0.000
25 25 .p4. == .p15. 2 0 0 0 NA 0 0.000 0.000 0.000
26 26 .p5. == .p16. 2 0 0 0 NA 0 0.000 0.000 0.000
27 27 .p6. == .p17. 2 0 0 0 NA 0 0.000 0.000 0.000
28 28 .p8. == .p19. 2 0 0 0 NA 0 0.000 0.000 0.000
29 29 .p9. == .p20. 2 0 0 0 NA 0 0.000 0.000 0.000
30 30 .p10. == .p21. 2 0 0 0 NA 0 0.000 0.000 0.000parTable indicates that .p5. == .p16 correspond to the the residual variance of gvslvue, namely gvslvue ~~ gvslvue. As such, we start freeing the residual variance of gvslvue.
fit_strict_gvslvue <- cfa(model_ws,
data = ess_df,
group = "gndr",
group.equal = c("loadings",
"intercepts",
"residuals"),
group.partial = c(gvslvue ~~ gvslvue)
)
lavTestScore(fit_strict_gvslvue)$test
total score test:
test X2 df p.value
1 score 10.301 7 0.172
$uni
univariate score tests:
lhs op rhs X2 df p.value
1 .p2. == .p13. 0.683 1 0.409
2 .p3. == .p14. 0.466 1 0.495
3 .p4. == .p15. 0.124 1 0.725
4 .p6. == .p17. 4.022 1 0.045
5 .p8. == .p19. 1.173 1 0.279
6 .p9. == .p20. 1.411 1 0.235
7 .p10. == .p21. 2.693 1 0.101Next, we run again a multivariate score test. In this case, we fail to reject multivariate score test. This means that nothing else should be freed. Let’s now see if we achieved partial strict invariance by running an anova test.
table_anova <-
list(anova(fit_configural, fit_metric),
anova(fit_metric, fit_scalar),
anova(fit_scalar, fit_strict_gvslvue)) %>%
reduce(rbind) %>%
.[-c(3, 5),]
table_anovaChi-Squared Difference Test
Df AIC BIC Chisq Chisq diff Df diff Pr(>Chisq)
fit_configural 0 18887 18985 0.0000
fit_metric 2 18885 18973 2.5551 2.5551 2 0.2787
fit_scalar 4 18885 18962 6.1987 3.6436 2 0.1617
fit_strict_gvslvue 6 18885 18951 10.4180 4.2193 2 0.1213We did!
Q: However, what does it means practically? Can we make comparisons based on the raw scores, if the residual/error variance is not the same across groups?
5 Multi-group SEM
Similarly to measurament invariance, we can test regression paths invariance. Multi-group SEM allows us to do it and compare the results from two or more groups. These groups could reflect experimental treatments, different locations of the data collection, different genders, or any other characteristics of interest. The goal of such an analysis is to assess whether the relationships among predictor and response variables vary by group. As such, it can be thought of as a “model-wide” interaction where we identify which paths (aka effects) change based on a group of interest (i.e. moderator).
In this lab, we are going to fit the same model that we fit in the previous lab but in a multi-group setting
- Y is the dependent variable (Welfare support)
- X is the predictor (Income)
- M is a mediator (Egalitarianism)
- G is a grouping variable (Gender)
First, we fit the “unconstrained” model. In this case, all the path coefficients are free to vary across groups. Since we already tested for invariance, we are going to set the loadings equal across groups.
model_mediation_mg <- '
## Welfare Support Factor ##
welf_supp =~ gvslvol + gvslvue + gvhlthc
## Egalitarianism ##
egual =~ gincdif + dfincac + smdfslv
## Direct effect ##
welf_supp ~ c("c1", "c2")*hinctnta
## Mediator ##
egual ~ c("a1", "a2")*hinctnta
welf_supp ~ c("b1", "b2")*egual
## Indirect effect (a*b) ##
a1b1 := a1*b1
a2b2 := a2*b2
## Total effect c + (a*b) ##
total1 := c1 + (a1*b1)
total2 := c2 + (a2*b2)
'
fit_mediation_mg <- cfa(model_mediation_mg, # model formula
data = ess_df, # data frame
group = "gndr", # grouping variable (G)
group.equal = c("loadings") # equal loadings
)
summary(fit_mediation_mg, standardized=TRUE)lavaan 0.6-10 ended normally after 73 iterations
Estimator ML
Optimization method NLMINB
Number of model parameters 42
Number of equality constraints 4
Number of observations per group: Used Total
Male 779 864
Female 773 896
Model Test User Model:
Test statistic 69.425
Degrees of freedom 28
P-value (Chi-square) 0.000
Test statistic for each group:
Male 55.965
Female 13.461
Parameter Estimates:
Standard errors Standard
Information Expected
Information saturated (h1) model Structured
Group 1 [Male]:
Latent Variables:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
welf_supp =~
gvslvol 1.000 1.333 0.868
gvslvue (.p2.) 0.576 0.048 11.887 0.000 0.768 0.381
gvhlthc (.p3.) 0.833 0.056 14.948 0.000 1.110 0.758
egual =~
gincdif 1.000 0.705 0.650
dfincac (.p5.) -0.600 0.060 -9.971 0.000 -0.424 -0.403
smdfslv (.p6.) 0.871 0.082 10.659 0.000 0.614 0.611
Regressions:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
welf_supp ~
hinctnta (c1) 0.017 0.023 0.706 0.480 0.012 0.029
egual ~
hinctnta (a1) 0.062 0.014 4.422 0.000 0.089 0.205
welf_supp ~
egual (b1) -0.487 0.104 -4.693 0.000 -0.258 -0.258
Intercepts:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
.gvslvol 7.989 0.182 43.942 0.000 7.989 5.204
.gvslvue 6.033 0.123 48.896 0.000 6.033 2.989
.gvhlthc 8.150 0.154 53.056 0.000 8.150 5.565
.gincdif 1.762 0.114 15.406 0.000 1.762 1.624
.dfincac 2.826 0.077 36.771 0.000 2.826 2.691
.smdfslv 2.096 0.101 20.763 0.000 2.096 2.087
.welf_supp 0.000 0.000 0.000
.egual 0.000 0.000 0.000
Variances:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
.gvslvol 0.581 0.122 4.744 0.000 0.581 0.246
.gvslvue 3.484 0.183 19.004 0.000 3.484 0.855
.gvhlthc 0.912 0.095 9.630 0.000 0.912 0.425
.gincdif 0.681 0.062 10.969 0.000 0.681 0.578
.dfincac 0.924 0.052 17.711 0.000 0.924 0.837
.smdfslv 0.632 0.051 12.428 0.000 0.632 0.626
.welf_supp 1.662 0.155 10.735 0.000 0.936 0.936
.egual 0.477 0.063 7.596 0.000 0.958 0.958
Group 2 [Female]:
Latent Variables:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
welf_supp =~
gvslvol 1.000 1.208 0.862
gvslvue (.p2.) 0.576 0.048 11.887 0.000 0.696 0.378
gvhlthc (.p3.) 0.833 0.056 14.948 0.000 1.006 0.697
egual =~
gincdif 1.000 0.674 0.654
dfincac (.p5.) -0.600 0.060 -9.971 0.000 -0.405 -0.383
smdfslv (.p6.) 0.871 0.082 10.659 0.000 0.587 0.623
Regressions:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
welf_supp ~
hinctnta (c2) 0.001 0.021 0.069 0.945 0.001 0.003
egual ~
hinctnta (a2) 0.050 0.013 3.807 0.000 0.074 0.176
welf_supp ~
egual (b2) -0.472 0.099 -4.796 0.000 -0.264 -0.264
Intercepts:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
.gvslvol 8.079 0.156 51.758 0.000 8.079 5.768
.gvslvue 6.207 0.108 57.408 0.000 6.207 3.368
.gvhlthc 8.147 0.134 60.932 0.000 8.147 5.642
.gincdif 1.815 0.102 17.781 0.000 1.815 1.759
.dfincac 2.952 0.070 42.177 0.000 2.952 2.791
.smdfslv 2.097 0.090 23.362 0.000 2.097 2.226
.welf_supp 0.000 0.000 0.000
.egual 0.000 0.000 0.000
Variances:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
.gvslvol 0.503 0.105 4.771 0.000 0.503 0.256
.gvslvue 2.912 0.155 18.828 0.000 2.912 0.857
.gvhlthc 1.073 0.090 11.912 0.000 1.073 0.515
.gincdif 0.610 0.057 10.757 0.000 0.610 0.573
.dfincac 0.955 0.053 17.940 0.000 0.955 0.854
.smdfslv 0.543 0.046 11.930 0.000 0.543 0.612
.welf_supp 1.358 0.130 10.463 0.000 0.931 0.931
.egual 0.440 0.058 7.640 0.000 0.969 0.969
Defined Parameters:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
a1b1 -0.030 0.009 -3.297 0.001 -0.023 -0.053
a2b2 -0.023 0.008 -3.044 0.002 -0.019 -0.047
total1 -0.014 0.023 -0.611 0.541 -0.010 -0.024
total2 -0.022 0.020 -1.088 0.276 -0.018 -0.044Second, we are going to fix both the loadings and the path coefficients in each group to be the same.
model_mediation_mg_cons <- '
## Welfare Support Factor ##
welf_supp =~ gvslvol + gvslvue + gvhlthc
## Egalitarianism ##
egual =~ gincdif + dfincac + smdfslv
## Direct effect ##
welf_supp ~ c("c1", "c1")*hinctnta
## Mediator ##
egual ~ c("a1", "a1")*hinctnta
welf_supp ~ c("b1", "b1")*egual
## Indirect effect (a*b) ##
a1b1 := a1*b1
## Total effect c + (a*b) ##
total1 := c1 + (a1*b1)
'
fit_mediation_mg_cons <- cfa(model_mediation_mg_cons, # model formula
data = ess_df, # data frame
group = "gndr", # grouping variable (G)
group.equal = c("loadings") # equal loadings
)
summary(fit_mediation_mg_cons, standardized=TRUE)lavaan 0.6-10 ended normally after 73 iterations
Estimator ML
Optimization method NLMINB
Number of model parameters 42
Number of equality constraints 7
Number of observations per group: Used Total
Male 779 864
Female 773 896
Model Test User Model:
Test statistic 70.047
Degrees of freedom 31
P-value (Chi-square) 0.000
Test statistic for each group:
Male 56.319
Female 13.727
Parameter Estimates:
Standard errors Standard
Information Expected
Information saturated (h1) model Structured
Group 1 [Male]:
Latent Variables:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
welf_supp =~
gvslvol 1.000 1.332 0.868
gvslvue (.p2.) 0.576 0.048 11.887 0.000 0.768 0.381
gvhlthc (.p3.) 0.833 0.056 14.942 0.000 1.110 0.758
egual =~
gincdif 1.000 0.706 0.652
dfincac (.p5.) -0.597 0.060 -9.948 0.000 -0.422 -0.402
smdfslv (.p6.) 0.863 0.081 10.647 0.000 0.609 0.607
Regressions:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
welf_supp ~
hinctnta (c1) 0.008 0.015 0.520 0.603 0.006 0.014
egual ~
hinctnta (a1) 0.056 0.010 5.710 0.000 0.079 0.183
welf_supp ~
egual (b1) -0.476 0.073 -6.480 0.000 -0.252 -0.252
Intercepts:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
.gvslvol 8.024 0.128 62.865 0.000 8.024 5.228
.gvslvue 6.053 0.098 61.546 0.000 6.053 2.999
.gvhlthc 8.179 0.109 74.709 0.000 8.179 5.585
.gincdif 1.813 0.084 21.620 0.000 1.813 1.673
.dfincac 2.795 0.060 46.251 0.000 2.795 2.662
.smdfslv 2.144 0.074 28.884 0.000 2.144 2.137
.welf_supp 0.000 0.000 0.000
.egual 0.000 0.000 0.000
Variances:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
.gvslvol 0.581 0.122 4.740 0.000 0.581 0.246
.gvslvue 3.484 0.183 19.004 0.000 3.484 0.855
.gvhlthc 0.913 0.095 9.626 0.000 0.913 0.426
.gincdif 0.676 0.062 10.857 0.000 0.676 0.575
.dfincac 0.924 0.052 17.710 0.000 0.924 0.838
.smdfslv 0.635 0.051 12.522 0.000 0.635 0.631
.welf_supp 1.664 0.154 10.776 0.000 0.938 0.938
.egual 0.482 0.063 7.649 0.000 0.966 0.966
Group 2 [Female]:
Latent Variables:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
welf_supp =~
gvslvol 1.000 1.208 0.862
gvslvue (.p2.) 0.576 0.048 11.887 0.000 0.696 0.378
gvhlthc (.p3.) 0.833 0.056 14.942 0.000 1.006 0.697
egual =~
gincdif 1.000 0.680 0.658
dfincac (.p5.) -0.597 0.060 -9.948 0.000 -0.406 -0.384
smdfslv (.p6.) 0.863 0.081 10.647 0.000 0.586 0.622
Regressions:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
welf_supp ~
hinctnta (c1) 0.008 0.015 0.520 0.603 0.007 0.016
egual ~
hinctnta (a1) 0.056 0.010 5.710 0.000 0.082 0.197
welf_supp ~
egual (b1) -0.476 0.073 -6.480 0.000 -0.268 -0.268
Intercepts:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
.gvslvol 8.054 0.121 66.367 0.000 8.054 5.749
.gvslvue 6.192 0.092 67.259 0.000 6.192 3.360
.gvhlthc 8.126 0.106 76.856 0.000 8.126 5.627
.gincdif 1.770 0.080 22.025 0.000 1.770 1.712
.dfincac 2.978 0.059 50.328 0.000 2.978 2.814
.smdfslv 2.061 0.071 29.078 0.000 2.061 2.184
.welf_supp 0.000 0.000 0.000
.egual 0.000 0.000 0.000
Variances:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
.gvslvol 0.504 0.105 4.780 0.000 0.504 0.257
.gvslvue 2.912 0.155 18.828 0.000 2.912 0.857
.gvhlthc 1.073 0.090 11.905 0.000 1.073 0.514
.gincdif 0.606 0.057 10.671 0.000 0.606 0.567
.dfincac 0.955 0.053 17.946 0.000 0.955 0.853
.smdfslv 0.546 0.045 12.054 0.000 0.546 0.614
.welf_supp 1.356 0.129 10.498 0.000 0.930 0.930
.egual 0.444 0.058 7.683 0.000 0.961 0.961
Defined Parameters:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
a1b1 -0.027 0.006 -4.474 0.000 -0.020 -0.046
total1 -0.018 0.015 -1.224 0.221 -0.014 -0.032Next, we compare the two models using the \(\chi^2\) difference test:
anova(fit_mediation_mg, fit_mediation_mg_cons)Chi-Squared Difference Test
Df AIC BIC Chisq Chisq diff Df diff Pr(>Chisq)
fit_mediation_mg 28 29664 29867 69.425
fit_mediation_mg_cons 31 29658 29846 70.047 0.62147 3 0.8915The insignificant P-value implies that the unconstrained and the constrained models are not statistically significantly different. This means that the path coefficients vary very little by group. In this case, we can analyse the pooled data in a single global model.
On the contrary, a significant P-value would imply that some paths vary while others may not. If this is the case, we can introduce constraints to identify which path varies between groups. Note that the SEM model should still fit the data well, regardless of the constraints imposed on the paths.
This exercise of relaxing and imposing constraints is potentially very exploratory. Refrain from constraining and relaxing all paths and then choosing the most parsimonious model. Instead, choosing which paths to constrain should be motivated by solid theoretical expectations and your research questions.