7.3 Invariance Over Time#
Measurement invariance across time for quantitative item responses#
# General imports
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
# Rpy2 imports
from rpy2 import robjects as ro
from rpy2.robjects import pandas2ri, numpy2ri
from rpy2.robjects.packages import importr
# Automatic conversion of arrays and dataframes
pandas2ri.activate()
numpy2ri.activate()
# Set random seed for reproducibility
ro.r('set.seed(123)')
# Ipython extenrsion for magix plotting
%load_ext rpy2.ipython
# R imports
importr('base')
importr('lavaan')
importr('semTools')
importr('MPsychoR')
importr('lordif')
importr('Hmisc')
c:\Users\maku1542\AppData\Local\miniconda3\envs\psy126\Lib\site-packages\rpy2\robjects\packages.py:367: UserWarning: The symbol 'quartz' is not in this R namespace/package.
warnings.warn(
rpy2.robjects.packages.Package as a <module 'Hmisc'>
CFA can also be extended to longitudinal settings, where indicators are presented to the same individual at multiple points in time \(t\). Across the time points, we can consider the same measurement invariance principles as in multigroup CFA (configural, weak, and strong invariance). However, in longitudinal CFA we need to account for correlated residuals (i.e. autoregressive effects) since we cannot assume that independence holds across time points.
Load, prepare and inspect the dataset#
ro.r('data("SDOwave")')
# Convert to Python
SDOwave = pandas2ri.rpy2py(ro.globalenv['SDOwave'])
print(SDOwave.head())
I1.1996 I2.1996 I3.1996 I4.1996 I1.1997 I2.1997 I3.1997 I4.1997 \
1 2.0 2.0 4.0 4.0 1.0 1.0 4.0 1.0
2 1.0 1.0 2.0 2.0 1.0 1.0 2.0 3.0
3 1.0 1.0 1.0 1.0 1.0 1.0 1.0 1.0
4 2.0 2.0 3.0 3.0 2.0 2.0 2.0 2.0
5 1.0 1.0 6.0 4.0 7.0 1.0 1.0 1.0
I1.1998 I2.1998 I3.1998 I4.1998 I1.1999 I2.1999 I3.1999 I4.1999 \
1 2.0 1.0 4.0 2.0 2.0 3.0 3.0 1.0
2 4.0 2.0 5.0 4.0 3.0 2.0 4.0 4.0
3 1.0 1.0 2.0 1.0 1.0 1.0 6.0 1.0
4 2.0 2.0 2.0 2.0 2.0 1.0 1.0 1.0
5 1.0 1.0 4.0 1.0 4.0 1.0 3.0 2.0
I1.2000 I2.2000 I3.2000 I4.2000
1 2.0 1.0 5.0 3.0
2 1.0 1.0 3.0 3.0
3 1.0 1.0 5.0 1.0
4 1.0 1.0 1.0 1.0
5 1.0 1.0 3.0 2.0
The dataset#
To illustrate longitudinal CFA, we use a dataset on social dominance orientation (SDO; Sidanius and Pratto, 2001). SDO is assessed on the same individuals from 1996 to 2000 (Sidanius et al., 2010), involving the following four items, scored on a 7-point scale:
It is probably a good thing that certain groups are at the top and other groups are at the bottom (I1).
Inferior groups should stay in their place (I2).
We should do what we can to equalize conditions for different groups (I3, reversed).
Increased social equality is beneficial to society (I4, reversed).
For the moment, let us consider a simple latent variable structure where all four items load on a general SDO dimension. We pick 2 years only: 1996 vs. 1998.
The most relevant model within this context is the strong/scale invariance model, which allows us to compare the latent means across the two time points. We restrict the factor loadings as well as the corresponding intercepts of each item to be equal in both measurement occasions (strong/scale invariance). Also, we need to allow for residual covariances (using the ~~ symbol in the syntax) for each item across time points. With this we account for autoregressive effects, which has to be done when modeling longitudinal data.
Fit the models#
We will start with the scale invariance model and work our way down to lower levels of measurement invariance. (We exclude the strict invariance model here).
We will start with the scale invariance model and work our way down to lower levels of measurement invariance. (We exclude the strict invariance model here).
Scale Invariance#
# Specify the model
ro.r("""model_scale <- '
SDO1996 =~ 1*I1.1996 + a2*I2.1996 + a3*I3.1996 + a4*I4.1996
SDO1998 =~ 1*I1.1998 + a2*I2.1998 + a3*I3.1998 + a4*I4.1998
SDO1996 ~~ SDO1998
I1.1996 ~ int1*1; I1.1998 ~ int1*1
I2.1996 ~ int2*1; I2.1998 ~ int2*1
I3.1996 ~ int3*1; I3.1998 ~ int3*1
I4.1996 ~ int4*1; I4.1998 ~ int4*1
I1.1996 ~~ I1.1998
I2.1996 ~~ I2.1998
I3.1996 ~~ I3.1998
I4.1996 ~~ I4.1998
SDO1996 ~ 0*1
SDO1998 ~ 1'
""")
# Fit the model
ro.r('fit_scale <- lavaan::cfa(model_scale, data = SDOwave, estimator = "MLR")')
# Print the model output with fit measures and standardized estimates
summary_fit_scale = ro.r('summary(fit_scale, fit.measures=TRUE, standardized=TRUE)')
print(summary_fit_scale)
lavaan 0.6-19 ended normally after 45 iterations
Estimator ML
Optimization method NLMINB
Number of model parameters 30
Number of equality constraints 7
Number of observations 612
Model Test User Model:
Standard Scaled
Test Statistic 315.748 203.519
Degrees of freedom 21 21
P-value (Chi-square) 0.000 0.000
Scaling correction factor 1.551
Yuan-Bentler correction (Mplus variant)
Model Test Baseline Model:
Test statistic 1993.536 1098.021
Degrees of freedom 28 28
P-value 0.000 0.000
Scaling correction factor 1.816
User Model versus Baseline Model:
Comparative Fit Index (CFI) 0.850 0.829
Tucker-Lewis Index (TLI) 0.800 0.773
Robust Comparative Fit Index (CFI) 0.854
Robust Tucker-Lewis Index (TLI) 0.806
Loglikelihood and Information Criteria:
Loglikelihood user model (H0) -7384.052 -7384.052
Scaling correction factor 1.434
for the MLR correction
Loglikelihood unrestricted model (H1) -7226.178 -7226.178
Scaling correction factor 1.718
for the MLR correction
Akaike (AIC) 14814.104 14814.104
Bayesian (BIC) 14915.689 14915.689
Sample-size adjusted Bayesian (SABIC) 14842.669 14842.669
Root Mean Square Error of Approximation:
RMSEA 0.151 0.119
90 Percent confidence interval - lower 0.137 0.107
90 Percent confidence interval - upper 0.166 0.131
P-value H_0: RMSEA <= 0.050 0.000 0.000
P-value H_0: RMSEA >= 0.080 1.000 1.000
Robust RMSEA 0.148
90 Percent confidence interval - lower 0.130
90 Percent confidence interval - upper 0.167
P-value H_0: Robust RMSEA <= 0.050 0.000
P-value H_0: Robust RMSEA >= 0.080 1.000
Standardized Root Mean Square Residual:
SRMR 0.096 0.096
Parameter Estimates:
Standard errors Sandwich
Information bread Observed
Observed information based on Hessian
Latent Variables:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
SDO1996 =~
I1.1996 1.000 0.434 0.319
I2.1996 (a2) 0.790 0.101 7.842 0.000 0.343 0.302
I3.1996 (a3) 2.949 0.326 9.042 0.000 1.281 0.845
I4.1996 (a4) 2.965 0.353 8.389 0.000 1.288 0.952
SDO1998 =~
I1.1998 1.000 0.407 0.326
I2.1998 (a2) 0.790 0.101 7.842 0.000 0.322 0.308
I3.1998 (a3) 2.949 0.326 9.042 0.000 1.200 0.813
I4.1998 (a4) 2.965 0.353 8.389 0.000 1.207 0.930
Covariances:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
SDO1996 ~~
SDO1998 0.098 0.023 4.275 0.000 0.552 0.552
.I1.1996 ~~
.I1.1998 0.425 0.075 5.697 0.000 0.425 0.280
.I2.1996 ~~
.I2.1998 0.258 0.072 3.594 0.000 0.258 0.240
.I3.1996 ~~
.I3.1998 0.220 0.076 2.911 0.004 0.220 0.316
.I4.1996 ~~
.I4.1998 -0.109 0.052 -2.078 0.038 -0.109 -0.555
Intercepts:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
.I1.1996 (int1) 2.064 0.045 45.732 0.000 2.064 1.516
.I1.1998 (int1) 2.064 0.045 45.732 0.000 2.064 1.654
.I2.1996 (int2) 1.633 0.039 42.279 0.000 1.633 1.438
.I2.1998 (int2) 1.633 0.039 42.279 0.000 1.633 1.562
.I3.1996 (int3) 2.765 0.060 46.233 0.000 2.765 1.825
.I3.1998 (int3) 2.765 0.060 46.233 0.000 2.765 1.874
.I4.1996 (int4) 2.428 0.054 44.827 0.000 2.428 1.795
.I4.1998 (int4) 2.428 0.054 44.827 0.000 2.428 1.872
SDO1996 0.000 0.000 0.000
SDO1998 -0.047 0.019 -2.432 0.015 -0.115 -0.115
Variances:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
.I1.1996 1.664 0.126 13.170 0.000 1.664 0.898
.I2.1996 1.172 0.131 8.960 0.000 1.172 0.909
.I3.1996 0.656 0.103 6.377 0.000 0.656 0.286
.I4.1996 0.170 0.079 2.145 0.032 0.170 0.093
.I1.1998 1.391 0.104 13.370 0.000 1.391 0.894
.I2.1998 0.989 0.144 6.874 0.000 0.989 0.905
.I3.1998 0.739 0.102 7.275 0.000 0.739 0.339
.I4.1998 0.226 0.082 2.756 0.006 0.226 0.134
SDO1996 0.189 0.043 4.363 0.000 1.000 1.000
SDO1998 0.166 0.035 4.756 0.000 1.000 1.000
The parameters of main interest are the latent variable means (however, this has nothing to do with measurement invariance but can be used to compare latent score bewteen time points):
print(ro.r("parameterEstimates(fit_scale)[22:23,]"))
lhs op rhs label est se z pvalue ci.lower ci.upper
22 SDO1996 ~1 0.000 0.000 NA NA 0.000 0.000
23 SDO1998 ~1 -0.047 0.019 -2.432 0.015 -0.085 -0.009
We see that the mean for 1996 was fixed to 0 and is therefore used as the reference year. The second line suggests that there is a significant decrease in SDO from 1996 to 1998.
Note that a weak/metric invariance version of this model can be fitted by restricting both
latent means to 0 and freeing up the intercepts. For a configural invariance version,
the loadings need to be freed up as well (first two syntax lines). For a strict invariance version, intercepts, loadings and residuals need to be fixed.
If we do not want to do this by hand, the longInvariance function from semTools can be used
to establish such a generic sequence (similar to the measurementInvariance function used above).
Using the anova function, the models can be again tested against each other.
Let us fit a weak/metric invariance model.
Metric Invariance#
Adapt the code from the scale model to fit a metric model by freeing up the intercepts and restrict both latent variables to zero. Name your model model.metric.
# Specify the model
ro.r("""model_metric <- '
SDO1996 =~ 1*I1.1996 + a2*I2.1996 + a3*I3.1996 + a4*I4.1996
SDO1998 =~ 1*I1.1998 + a2*I2.1998 + a3*I3.1998 + a4*I4.1998
SDO1996 ~~ SDO1998
I1.1996 ~ 1; I1.1998 ~ 1
I2.1996 ~ 1; I2.1998 ~ 1
I3.1996 ~ 1; I3.1998 ~ 1
I4.1996 ~ 1; I4.1998 ~ 1
I1.1996 ~~ I1.1998
I2.1996 ~~ I2.1998
I3.1996 ~~ I3.1998
I4.1996 ~~ I4.1998
SDO1996 ~ 0*1
SDO1998 ~ 0*1'
""")
# Fit the model
ro.r('fit_metric <- lavaan::cfa(model_metric, data = SDOwave, estimator = "MLR")')
# Print the model output with fit measures and standardized estimates
summary_fit_metric = ro.r('summary(fit_metric, fit.measures=TRUE, standardized=TRUE)')
print(summary_fit_metric)
lavaan 0.6-19 ended normally after 43 iterations
Estimator ML
Optimization method NLMINB
Number of model parameters 29
Number of equality constraints 3
Number of observations 612
Model Test User Model:
Standard Scaled
Test Statistic 310.762 189.360
Degrees of freedom 18 18
P-value (Chi-square) 0.000 0.000
Scaling correction factor 1.641
Yuan-Bentler correction (Mplus variant)
Model Test Baseline Model:
Test statistic 1993.536 1098.021
Degrees of freedom 28 28
P-value 0.000 0.000
Scaling correction factor 1.816
User Model versus Baseline Model:
Comparative Fit Index (CFI) 0.851 0.840
Tucker-Lewis Index (TLI) 0.768 0.751
Robust Comparative Fit Index (CFI) 0.855
Robust Tucker-Lewis Index (TLI) 0.775
Loglikelihood and Information Criteria:
Loglikelihood user model (H0) -7381.559 -7381.559
Scaling correction factor 1.588
for the MLR correction
Loglikelihood unrestricted model (H1) -7226.178 -7226.178
Scaling correction factor 1.718
for the MLR correction
Akaike (AIC) 14815.117 14815.117
Bayesian (BIC) 14929.952 14929.952
Sample-size adjusted Bayesian (SABIC) 14847.408 14847.408
Root Mean Square Error of Approximation:
RMSEA 0.163 0.125
90 Percent confidence interval - lower 0.147 0.112
90 Percent confidence interval - upper 0.179 0.137
P-value H_0: RMSEA <= 0.050 0.000 0.000
P-value H_0: RMSEA >= 0.080 1.000 1.000
Robust RMSEA 0.160
90 Percent confidence interval - lower 0.140
90 Percent confidence interval - upper 0.181
P-value H_0: Robust RMSEA <= 0.050 0.000
P-value H_0: Robust RMSEA >= 0.080 1.000
Standardized Root Mean Square Residual:
SRMR 0.096 0.096
Parameter Estimates:
Standard errors Sandwich
Information bread Observed
Observed information based on Hessian
Latent Variables:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
SDO1996 =~
I1.1996 1.000 0.437 0.321
I2.1996 (a2) 0.781 0.101 7.726 0.000 0.341 0.301
I3.1996 (a3) 2.946 0.326 9.044 0.000 1.287 0.848
I4.1996 (a4) 2.938 0.348 8.447 0.000 1.284 0.950
SDO1998 =~
I1.1998 1.000 0.410 0.328
I2.1998 (a2) 0.781 0.101 7.726 0.000 0.320 0.306
I3.1998 (a3) 2.946 0.326 9.044 0.000 1.206 0.816
I4.1998 (a4) 2.938 0.348 8.447 0.000 1.203 0.928
Covariances:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
SDO1996 ~~
SDO1998 0.099 0.023 4.305 0.000 0.553 0.553
.I1.1996 ~~
.I1.1998 0.425 0.075 5.700 0.000 0.425 0.279
.I2.1996 ~~
.I2.1998 0.259 0.072 3.618 0.000 0.259 0.241
.I3.1996 ~~
.I3.1998 0.217 0.076 2.871 0.004 0.217 0.316
.I4.1996 ~~
.I4.1998 -0.104 0.052 -1.991 0.046 -0.104 -0.507
Intercepts:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
.I1.1996 2.042 0.056 36.593 0.000 2.042 1.500
.I1.1998 2.034 0.050 40.688 0.000 2.034 1.630
.I2.1996 1.673 0.047 35.833 0.000 1.673 1.475
.I2.1998 1.564 0.042 37.580 0.000 1.564 1.498
.I3.1996 2.742 0.061 45.150 0.000 2.742 1.807
.I3.1998 2.655 0.061 43.803 0.000 2.655 1.797
.I4.1996 2.441 0.055 44.387 0.000 2.441 1.806
.I4.1998 2.273 0.052 43.523 0.000 2.273 1.753
SDO1996 0.000 0.000 0.000
SDO1998 0.000 0.000 0.000
Variances:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
.I1.1996 1.663 0.127 13.068 0.000 1.663 0.897
.I2.1996 1.170 0.130 9.031 0.000 1.170 0.910
.I3.1996 0.646 0.105 6.166 0.000 0.646 0.280
.I4.1996 0.180 0.079 2.273 0.023 0.180 0.098
.I1.1998 1.390 0.104 13.420 0.000 1.390 0.892
.I2.1998 0.987 0.144 6.837 0.000 0.987 0.906
.I3.1998 0.729 0.100 7.310 0.000 0.729 0.334
.I4.1998 0.234 0.081 2.896 0.004 0.234 0.139
SDO1996 0.191 0.044 4.372 0.000 1.000 1.000
SDO1998 0.168 0.035 4.802 0.000 1.000 1.000
Note that the values in the intercepts section are now varying.
Again, we can extract latent variable means (however, this has again nothing to do with measurement invariance):
print(ro.r("parameterEstimates(fit_metric)[22:23,]"))
lhs op rhs label est se z pvalue ci.lower ci.upper
22 SDO1996 ~1 0 0 NA NA 0 0
23 SDO1998 ~1 0 0 NA NA 0 0
Let us continue with a configural variance model.
Configural Invariance#
Adapt the code from the metric model to fit a configural model by freeing up the loadings as well (first two syntax lines). Name your model model.configural.
# Specify the model
ro.r("""model_configural <- '
SDO1996 =~ 1*I1.1996 + I2.1996 + I3.1996 + I4.1996
SDO1998 =~ 1*I1.1998 + I2.1998 + I3.1998 + I4.1998
SDO1996 ~~ SDO1998
I1.1996 ~ 1; I1.1998 ~ 1
I2.1996 ~ 1; I2.1998 ~ 1
I3.1996 ~ 1; I3.1998 ~ 1
I4.1996 ~ 1; I4.1998 ~ 1
I1.1996 ~~ I1.1998
I2.1996 ~~ I2.1998
I3.1996 ~~ I3.1998
I4.1996 ~~ I4.1998
SDO1996 ~ 0*1
SDO1998 ~ 0*1'
""")
# Fit the model
ro.r('fit_configural <- lavaan::cfa(model_configural, data = SDOwave, estimator = "MLR")')
# Print the model output with fit measures and standardized estimates
summary_fit_configural = ro.r('summary(fit_configural, fit.measures=TRUE, standardized=TRUE)')
print(summary_fit_configural)
lavaan 0.6-19 ended normally after 61 iterations
Estimator ML
Optimization method NLMINB
Number of model parameters 29
Number of observations 612
Model Test User Model:
Standard Scaled
Test Statistic 301.687 175.637
Degrees of freedom 15 15
P-value (Chi-square) 0.000 0.000
Scaling correction factor 1.718
Yuan-Bentler correction (Mplus variant)
Model Test Baseline Model:
Test statistic 1993.536 1098.021
Degrees of freedom 28 28
P-value 0.000 0.000
Scaling correction factor 1.816
User Model versus Baseline Model:
Comparative Fit Index (CFI) 0.854 0.850
Tucker-Lewis Index (TLI) 0.728 0.720
Robust Comparative Fit Index (CFI) 0.858
Robust Tucker-Lewis Index (TLI) 0.735
Loglikelihood and Information Criteria:
Loglikelihood user model (H0) -7377.021 -7377.021
Scaling correction factor 1.719
for the MLR correction
Loglikelihood unrestricted model (H1) -7226.178 -7226.178
Scaling correction factor 1.718
for the MLR correction
Akaike (AIC) 14812.043 14812.043
Bayesian (BIC) 14940.128 14940.128
Sample-size adjusted Bayesian (SABIC) 14848.059 14848.059
Root Mean Square Error of Approximation:
RMSEA 0.177 0.132
90 Percent confidence interval - lower 0.160 0.119
90 Percent confidence interval - upper 0.194 0.146
P-value H_0: RMSEA <= 0.050 0.000 0.000
P-value H_0: RMSEA >= 0.080 1.000 1.000
Robust RMSEA 0.173
90 Percent confidence interval - lower 0.151
90 Percent confidence interval - upper 0.197
P-value H_0: Robust RMSEA <= 0.050 0.000
P-value H_0: Robust RMSEA >= 0.080 1.000
Standardized Root Mean Square Residual:
SRMR 0.092 0.092
Parameter Estimates:
Standard errors Sandwich
Information bread Observed
Observed information based on Hessian
Latent Variables:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
SDO1996 =~
I1.1996 1.000 0.493 0.358
I2.1996 0.814 0.088 9.287 0.000 0.401 0.349
I3.1996 2.562 0.314 8.164 0.000 1.264 0.842
I4.1996 2.622 0.339 7.734 0.000 1.293 0.952
SDO1998 =~
I1.1998 1.000 0.363 0.294
I2.1998 0.729 0.177 4.124 0.000 0.265 0.258
I3.1998 3.452 0.561 6.154 0.000 1.254 0.837
I4.1998 3.242 0.568 5.708 0.000 1.178 0.912
Covariances:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
SDO1996 ~~
SDO1998 0.099 0.023 4.318 0.000 0.553 0.553
.I1.1996 ~~
.I1.1998 0.423 0.074 5.745 0.000 0.423 0.279
.I2.1996 ~~
.I2.1998 0.257 0.071 3.628 0.000 0.257 0.240
.I3.1996 ~~
.I3.1998 0.210 0.075 2.813 0.005 0.210 0.316
.I4.1996 ~~
.I4.1998 -0.095 0.052 -1.842 0.065 -0.095 -0.436
Intercepts:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
.I1.1996 2.042 0.056 36.593 0.000 2.042 1.483
.I1.1998 2.034 0.050 40.688 0.000 2.034 1.649
.I2.1996 1.673 0.047 35.833 0.000 1.673 1.455
.I2.1998 1.564 0.042 37.580 0.000 1.564 1.520
.I3.1996 2.742 0.061 45.150 0.000 2.742 1.827
.I3.1998 2.655 0.061 43.803 0.000 2.655 1.771
.I4.1996 2.441 0.055 44.387 0.000 2.441 1.798
.I4.1998 2.273 0.052 43.523 0.000 2.273 1.761
SDO1996 0.000 0.000 0.000
SDO1998 0.000 0.000 0.000
Variances:
Estimate Std.Err z-value P(>|z|) Std.lv Std.all
.I1.1996 1.652 0.128 12.945 0.000 1.652 0.872
.I2.1996 1.161 0.130 8.930 0.000 1.161 0.878
.I3.1996 0.655 0.097 6.771 0.000 0.655 0.291
.I4.1996 0.171 0.078 2.210 0.027 0.171 0.093
.I1.1998 1.390 0.103 13.489 0.000 1.390 0.913
.I2.1998 0.988 0.144 6.846 0.000 0.988 0.934
.I3.1998 0.675 0.114 5.925 0.000 0.675 0.300
.I4.1998 0.279 0.097 2.866 0.004 0.279 0.167
SDO1996 0.243 0.058 4.205 0.000 1.000 1.000
SDO1998 0.132 0.041 3.257 0.001 1.000 1.000
Model Comparison#
To assess which model provides a better fit (and therefore to assess which level of measurement invariance holds) we can use the anova() function.
Let us first compare the configural model and the metric model.
print(ro.r("anova(fit_configural, 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_configural 15 14812 14940 301.69
fit_metric 18 14815 14930 310.76 7.2119 3 0.06544 .
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
We see that the metric model fits worse than the configural model. This is no suprise as we include more restrictions in the metric model. However, the fit is not significantly worse (p = .0654) and the BIC even favors the metric model. From this we may conclude that metric invariance holds but there is a huge BUT (see below).
Next, lets compare the metric model and the scale model.
print(ro.r("anova(fit_metric, fit_scale)"))
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_metric 18 14815 14930 310.76
fit_scale 21 14814 14916 315.75 4.9207 3 0.1777
From this comparison we can conclude that scale invariance holds as the scale model is not significantly worse than the metric model.
Now to the BUT: You can see in the model output of the configural model that the fit is quite bad (e.g. RMSEA = 0.177 (> .08)). This is not too surprising since research has shown that there are two subdimensions of SDO (Ho et al., 2015): anti-egalitarianism and dominance. According to the theory, the first two items are supposed to load on the dominance dimension (SDO-D), whereas the remaining two items load on anti-egalitarianism (SDO-E).
Therefore, we should be cautious with interpreting the comparisons above and first make sure that we find a good fitting model to our data.