CFA measurement invariance

graphics.off()
options(width=500)
options(digits=3)

Data simulation

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"

Data summary

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

Model spesification

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 and model fit

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

Testing for invariance

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.

Test for STRUCTURAL INVARIANCE

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

Test for CONFIGURAL INVARIANCE

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.

Test for WEAK METRIC 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

Test for STRONG SCALAR INVARIANCE

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

Test for STRICT ERROR VARIANCE INVARIANCE

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

Test for FACTOR VARIANCE COVARIANCE

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

Test for FACTOR MEAN INVARIANCE

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

Comparing models

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