graphics.off()
options(width=500)
options(digits=3)
The data are simulated in lavaan. From this simulation it is possible to extract 3 factors and 2 arbitrary groups.
# specify population model
population_model<-'f1=~x1+.6*x2+.7*x3
f2=~x4+.6*x5+.7*x6
f3=~x7+.6*x8+.7*x9'
# generate data
set.seed(1234)
model_data<-lavaan::simulateData(population_model,sample.nobs=20000)
model_data$group<-NA
model_data[1:10000,]$group<-"Group1"
model_data[10001:20000,]$group<-"Group2"
psycholatefunctions::check(model_data,name_length=5)
## Loading required package: future.apply
## Loading required package: future
## $summary
## COLLUMNS ROWS TOTAL EMPTY NULL NAN NA INFINITE FINITE FACTOR
## 1 10 20000 200000 0 0 0 0 0 0 0
##
## $check
## NAMES POS EMPTY NULL NA NOT_NA NAN INFINITE FINITE RANGE MEAN MEDIAN SD MIN MAX MODE TYPE FACTOR
## 1 x1 1 0 0 0 20000 0 0 0 19958 <NA> <NA> <NA> - character character FALSE
## 2 x2 2 0 0 0 20000 0 0 0 3349 <NA> <NA> <NA> - character character FALSE
## 3 x3 3 0 0 0 20000 0 0 0 3336 <NA> <NA> <NA> - character character FALSE
## 4 x4 4 0 0 0 20000 0 0 0 19949 <NA> <NA> <NA> - character character FALSE
## 5 x5 5 0 0 0 20000 0 0 0 3304 <NA> <NA> <NA> - character character FALSE
## 6 x6 6 0 0 0 20000 0 0 0 3371 <NA> <NA> <NA> - character character FALSE
## 7 x7 7 0 0 0 20000 0 0 0 19956 <NA> <NA> <NA> - character character FALSE
## 8 x8 8 0 0 0 20000 0 0 0 3357 <NA> <NA> <NA> - character character FALSE
## 9 x9 9 0 0 0 20000 0 0 0 3341 <NA> <NA> <NA> - character character FALSE
## 10 group 10 0 0 0 20000 0 0 0 2 <NA> <NA> <NA> G G character character FALSE
The spesified model is shown in the figure bellow.
bf_model<-'f1=~x1+x2+x3
f2=~x4+x5+x6
f3=~x7+x8+x9'
cat(bf_model)
## f1=~x1+x2+x3
## f2=~x4+x5+x6
## f3=~x7+x8+x9
bfa<-lavaan::cfa(bf_model,data=model_data)
semPlot::semPaths(bfa,"model",layout="tree2",what="est")
## Registered S3 methods overwritten by 'huge':
## method from
## plot.sim BDgraph
## print.sim BDgraph
Model estimates are provided by three main matrices:
1. The lambda matrix representing factor loadings
2. The theta matrix representing standard errors for items in the diagonal, and error covariances between items in the off diagonal
3. The psi matrix representing covariances between factors
Model estimates are derived from a converged CFA model which in turn can produce a covariance matrix by performing these computations: covariance=lambda*psi+theta, which are the estimated data. In addition, a covariance matrix can be computed directly from item responses, which are the observed data. Because of that, we can compare the covariance matrix estimated from the model and the covariance matrix computed from item responses. In a good model, these two matrices should have similar values. Fit indices, shown bellow quantify model fit. There is conflicting information and several opinions on critical values of fit indices. In our case we consider good model fit to be indicated by:
1. RMSEA < .07
2. CFI > .9
3. TLI > .9
4. SRMR < .8
5. GFI > .95
lavaan::fitMeasures(bfa,c("rmsea","cfi","tli","rni","rfi","ifi","srmr","gfi"))
## rmsea cfi tli rni rfi ifi srmr gfi
## 0.000 1.000 1.000 1.000 0.999 1.000 0.005 1.000
lavaan::inspect(bfa,"est",add.labels=TRUE,add.class=TRUE,list.by.group=TRUE,drop.list.single.group=TRUE)
## $lambda
## f1 f2 f3
## x1 1.000 0.000 0.000
## x2 0.621 0.000 0.000
## x3 0.729 0.000 0.000
## x4 0.000 1.000 0.000
## x5 0.000 0.599 0.000
## x6 0.000 0.705 0.000
## x7 0.000 0.000 1.000
## x8 0.000 0.000 0.569
## x9 0.000 0.000 0.674
##
## $theta
## x1 x2 x3 x4 x5 x6 x7 x8 x9
## x1 1.038
## x2 0.000 1.004
## x3 0.000 0.000 0.980
## x4 0.000 0.000 0.000 0.995
## x5 0.000 0.000 0.000 0.000 0.984
## x6 0.000 0.000 0.000 0.000 0.000 1.000
## x7 0.000 0.000 0.000 0.000 0.000 0.000 0.958
## x8 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.984
## x9 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 1.014
##
## $psi
## f1 f2 f3
## f1 0.956
## f2 -0.009 1.006
## f3 -0.007 -0.003 1.043
We perform 2 seperate CFA parameter estimations for each group. If these models display significantly different lambda, theta, and psi estimates, that suggests non invariant models.
For structural invariance, we create 2 models one using Group1 data, the other using the Group2 data, and their respective fit indices are inspected.
# Model with the first group responses
ft<-lavaan::cfa(bf_model,data=model_data[model_data$group %in% "Group1",])
# Model with the second group responses
st<-lavaan::cfa(bf_model,data=model_data[model_data$group %in% "Group2",])
data.frame(Group1=lavaan::fitMeasures(ft,c("rmsea","cfi","tli","rni","rfi","ifi","srmr","gfi")),
Group2=lavaan::fitMeasures(st,c("rmsea","cfi","tli","rni","rfi","ifi","srmr","gfi")))
## Group1 Group2
## rmsea 0.0000 0.00507
## cfi 1.0000 0.99942
## tli 1.0012 0.99913
## rni 1.0008 0.99942
## rfi 0.9978 0.99578
## ifi 1.0008 0.99942
## srmr 0.0049 0.00796
## gfi 0.9997 0.99933
For configural invariance we spesify a model with 2 groups. This procedure estimates 2 lambda theta psi and mu (intercepts) matrices, for each group. By comparing the two lambda and theta matrices we notice that their values are not identical. The question is if the differences between these values are significantly different. In any case, the configural invariance model bellow displays similar fit indices compared with the structural invariance models above. But more importantly, according to fit indices the model fits the data and this suggests configural invariance.
# CONFIGURAL INVARIANCE
bf<-lavaan::cfa(bf_model,data=model_data,group="group")
lavaan::fitMeasures(bf,c("rmsea","cfi","tli","rni","rfi","ifi","srmr","gfi"))
## rmsea cfi tli rni rfi ifi srmr gfi
## 0.000 1.000 1.000 1.000 0.997 1.000 0.006 0.999
lavaan::inspect(bf,"est",add.labels=TRUE,add.class=TRUE,list.by.group=TRUE,drop.list.single.group=TRUE)
## $Group1
## $Group1$lambda
## f1 f2 f3
## x1 1.000 0.000 0.00
## x2 0.617 0.000 0.00
## x3 0.711 0.000 0.00
## x4 0.000 1.000 0.00
## x5 0.000 0.594 0.00
## x6 0.000 0.695 0.00
## x7 0.000 0.000 1.00
## x8 0.000 0.000 0.58
## x9 0.000 0.000 0.68
##
## $Group1$theta
## x1 x2 x3 x4 x5 x6 x7 x8 x9
## x1 1.037
## x2 0.000 0.995
## x3 0.000 0.000 0.967
## x4 0.000 0.000 0.000 0.986
## x5 0.000 0.000 0.000 0.000 0.983
## x6 0.000 0.000 0.000 0.000 0.000 0.988
## x7 0.000 0.000 0.000 0.000 0.000 0.000 0.969
## x8 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.982
## x9 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 1.027
##
## $Group1$psi
## f1 f2 f3
## f1 0.990
## f2 -0.020 1.013
## f3 -0.006 -0.022 0.997
##
## $Group1$nu
## intrcp
## x1 -0.001
## x2 -0.017
## x3 0.007
## x4 0.006
## x5 -0.005
## x6 0.010
## x7 0.004
## x8 0.007
## x9 0.008
##
## $Group1$alpha
## intrcp
## f1 0
## f2 0
## f3 0
##
##
## $Group2
## $Group2$lambda
## f1 f2 f3
## x1 1.000 0.000 0.000
## x2 0.626 0.000 0.000
## x3 0.748 0.000 0.000
## x4 0.000 1.000 0.000
## x5 0.000 0.604 0.000
## x6 0.000 0.716 0.000
## x7 0.000 0.000 1.000
## x8 0.000 0.000 0.559
## x9 0.000 0.000 0.670
##
## $Group2$theta
## x1 x2 x3 x4 x5 x6 x7 x8 x9
## x1 1.039
## x2 0.000 1.012
## x3 0.000 0.000 0.992
## x4 0.000 0.000 0.000 1.004
## x5 0.000 0.000 0.000 0.000 0.984
## x6 0.000 0.000 0.000 0.000 0.000 1.012
## x7 0.000 0.000 0.000 0.000 0.000 0.000 0.948
## x8 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.985
## x9 0.000 0.000 0.000 0.000 0.000 0.000 0.000 0.000 1.000
##
## $Group2$psi
## f1 f2 f3
## f1 0.922
## f2 0.002 0.998
## f3 -0.009 0.016 1.088
##
## $Group2$nu
## intrcp
## x1 -0.016
## x2 0.007
## x3 -0.008
## x4 0.015
## x5 -0.006
## x6 0.021
## x7 -0.007
## x8 -0.009
## x9 -0.008
##
## $Group2$alpha
## intrcp
## f1 0
## f2 0
## f3 0
par(mfrow=c(1,2))
semPlot::semPaths(bf,"model",layout="tree2",what="est")
In order to test if differences between estimated matrices among models are significant, we have to create a series of models for which impose equality constraints for parameters between groups. By imposing constraints we impose the estimator procedure to keep a matrix between groups equal. For instance by imposing an equality constraint on loadings, we force the parameter estimator to produce models that have equal lambda matrices between groups. Finally we compare the unconstrained and constrained models to see if they have any differences. If no significant differences exist, this suggest that the models are equivalent between groups and thus the models display measurement invariance.
For weak metric invariance we impose a restriction that groups have equal loadings.
# WEAK METRIC INVARIANCE
bf_metric<-lavaan::cfa(bf_model,data=model_data,group="group",group.equal=c("loadings"))
lavaan::anova(bf,bf_metric)
## Chi Square Difference Test
##
## Df AIC BIC Chisq Chisq diff Df diff Pr(>Chisq)
## bf 48 573054 573528 46.0
## bf_metric 54 573044 573471 48.3 2.31 6 0.89
For strong scalar invariance we impose the restriction that groups have equal loadings and intercepts and then compare the strong scalar invariance model with the weak metric invariance model.
# STRONG SCALAR INVARIANCE
bf_scalar<-lavaan::cfa(bf_model,data=model_data,group="group",group.equal=c("loadings","intercepts"))
lavaan::anova(bf_metric,bf_scalar)
## Chi Square Difference Test
##
## Df AIC BIC Chisq Chisq diff Df diff Pr(>Chisq)
## bf_metric 54 573044 573471 48.3
## bf_scalar 60 573038 573417 53.9 5.56 6 0.47
For strict error variance invariance, we impose the restriction that groups have equal loadings intercepts and residuals.
# STRICT ERROR VARIANCE INVARIANCE
bf_full_strict<-lavaan::cfa(bf_model,data=model_data,group="group",group.equal=c("loadings","intercepts","residuals"))
lavaan::anova(bf_scalar,bf_full_strict)
## Chi Square Difference Test
##
## Df AIC BIC Chisq Chisq diff Df diff Pr(>Chisq)
## bf_scalar 60 573038 573417 53.9
## bf_full_strict 69 573028 573336 62.0 8.13 9 0.52
For factor variance covariance invariance we impose the restriction that groups have equal loadings intercepts residuals latent variable variances and latent variable covariances.
# FACTOR VARIANCE COVARIANCE
bf_factor_vc<-lavaan::cfa(bf_model,data=model_data,group="group",group.equal=c("loadings","intercepts","residuals","lv.variances","lv.covariances"))
lavaan::anova(bf_full_strict,bf_factor_vc)
## Chi Square Difference Test
##
## Df AIC BIC Chisq Chisq diff Df diff Pr(>Chisq)
## bf_full_strict 69 573028 573336 62.0
## bf_factor_vc 75 573025 573286 70.8 8.86 6 0.18
For factor mean invariance we impose the restriction that groups have equal loadings intercepts and residuals and we impose the restriction that factors have equal variances covariances and means between groups.
# FACTOR MEAN INVARIANCE
bf_factor_mean<-lavaan::cfa(bf_model,data=model_data,group="group",group.equal=c("loadings","intercepts","residuals","lv.variances","lv.covariances","means"))
lavaan::anova(bf_factor_vc,bf_factor_mean)
## Chi Square Difference Test
##
## Df AIC BIC Chisq Chisq diff Df diff Pr(>Chisq)
## bf_factor_vc 75 573025 573286 70.8
## bf_factor_mean 78 573020 573257 72.1 1.29 3 0.73
Finally we compare models using an ANOVA procedure where values above 0.05 indicate no significant differences and thus measurement invariance.
lavaan::anova(bf,bf_metric,bf_scalar,bf_factor_vc,bf_factor_mean,SB.classic=TRUE)
## Chi Square Difference Test
##
## Df AIC BIC Chisq Chisq diff Df diff Pr(>Chisq)
## bf 48 573054 573528 46.0
## bf_metric 54 573044 573471 48.3 2.31 6 0.89
## bf_scalar 60 573038 573417 53.9 5.56 6 0.47
## bf_factor_vc 75 573025 573286 70.8 16.99 15 0.32
## bf_factor_mean 78 573020 573257 72.1 1.29 3 0.73
In addition, the fit indices for each model above are summarized bellow.
findex<-c("rmsea","cfi","tli","rni","rfi","ifi","srmr","gfi")
fit_indices<-data.frame(rbind(lavaan::fitMeasures(bf,findex),
lavaan::fitMeasures(bf_metric,findex),
lavaan::fitMeasures(bf_scalar,findex),
lavaan::fitMeasures(bf_factor_vc,findex),
lavaan::fitMeasures(bf_factor_mean,findex)))
row.names(fit_indices)<-c("CONFIGURAL INVARIANCE",
"WEAK METRIC INVARIANCE",
"STRONG SCALAR INVARIANCE",
"FACTOR VARIANCE COVARIANCE",
"FACTOR MEAN INVARIANCE")
fit_indices
## rmsea cfi tli rni rfi ifi srmr gfi
## CONFIGURAL INVARIANCE 0 1 1 1 0.997 1 0.00587 0.999
## WEAK METRIC INVARIANCE 0 1 1 1 0.997 1 0.00605 0.999
## STRONG SCALAR INVARIANCE 0 1 1 1 0.997 1 0.00635 0.999
## FACTOR VARIANCE COVARIANCE 0 1 1 1 NA 1 0.00892 0.999
## FACTOR MEAN INVARIANCE 0 1 1 1 NA 1 0.00904 0.999