Compatibility of ID.cat = "Wu.Estabrook.2016" with measEq.syntax and Multiple-Group SEM

[ahmad]

Hi,

I’m running a multiple-group SEM with both categorical and continuous indicators. Prior to this, I established measurement invariance across groups using measEq.syntax. For identification, I used:

parameterization = "theta",
estimator = "WLSMV",
ID.fac = "std.lv",
ID.cat = "Wu.Estabrook.2016",
group = "sex"

This approach is highly recommended for testing measurement invariance of categorical indicators.

After establishing measurement invariance, I now want to compare regression coefficients (structural paths) across groups, not latent means. Since “standardizing common factors by fixing their variances to 1.0 is incorrect if groups differ in their variabilities” (Kline, 2023), I’m concerned that using std.lv (which standardizes latent variances in each group) may obscure true differences in factor variances and make path comparisons less meaningful.

My question is: What is the recommended here to do? Shoudl I maintain this method while transitioning to multiple-group SEM for structural paths? Or is it more appropriate to switch to fixing the loading for the reference indicator to 1.0 (i.e., leave the factor unstandardized) in multiple-group SEM?

Thank you very much for your help in advance.

Best,
Ahmad

[Andres Perez]

Hi Ahmad,

I have found a similar problem before. There may be better solutions, but the most straightforward approach I found was to rescale the standardized factor covariance matrices resulting from the invariant model (let's refer to them as psi_std) into unstandardized ones (psi_uns). Then I used those psi_uns to run a multigroup SEM for the structural path. So, a stepwise estimation, as in the Structural-After-Measurement approach (Rosseel & Loh, 2022). In summary:

1. Test for invariance with a CFA model.

2. Extract the factor covariance matrices from the final model.

3. Rescale the standardized covariances into unstandardized covariances.

4. Use the unstandardized covariances as input for a second step where only the structural paths are estimated.

You can see Appendix A of a preprint of a paper I recently wrote (it is not the most updated version, sorry, but the rescaling did not change): https://doi.org/10.31234/osf.io/6cr74_v1 

You can rescale all the (co) variances by hand or use a function like cor2cov() to rescale the standardized factor covariances. The standard deviations are equal to the (would have been) marker loading^2 * the std factor variance.

Please see below a small example of how to do this. I did not use categorical data in this case, but the idea is the same.

library(lavaan)
# The Holzinger and Swineford (1939) example
HS.model <- ' visual  =~ x1 + x2 + x3
              textual =~ x4 + x5 + x6
              speed   =~ x7 + x8 + x9 '

fit <- cfa(model = HS.model,
           data  = HolzingerSwineford1939,
           group = "school",
           group.equal = "loadings",
           std.lv      = T)

# Extract loadings
loadings <- lapply(lavInspect(object = fit, what = "est"), "[[", "lambda")
# Extract marker loadings: first non-zero loading per factor (metric invariant model anyway)
markers  <- lapply(loadings, FUN = \(x) apply(x, 2, function(col) {col[which(col != 0)]}[1]))  

# Extract standardized factor variances (first group has a correlation matrix, but remaining groups may not have correlations depending on the constraints)
psi_std <- lapply(lavInspect(object = fit, what = "est"), "[[", "psi")

# Compute standard deviations of the unstandardized factors
sds     <- lapply(1:2, FUN = \(x) sqrt(diag(psi_std[[x]]) * markers[[x]]^2))

# Rescale the factor covariances
# To avoid computing everything by hand, use cor2cov()
# Note that the second group is not yet a correlation, so I first transform it to a correlation
psi_std[[2]] <- cov2cor(psi_std[[2]])

# Time to rescale
psi_uns <- lapply(1:2, FUN = \(x) cor2cov(psi_std[[x]], sds = sds[[x]]))


######################################
# Compare to a model where we estimate using the marker variable approach as a 'benchmark'
fit_marker <- cfa(model = HS.model,
                  data  = HolzingerSwineford1939,
                  group = "school",
                  group.equal = "loadings",
                  std.lv      = F)

psi_marker <- lapply(lavInspect(object = fit_marker, what = "est"), "[[", "psi")

# comparing
lapply(1:2, FUN = \(x) round(psi_marker[[x]], 5) == round(psi_uns[[x]], 5)) # I round the results to avoid differences in floating points

# Also possible to visually see that the covariances are the same
psi_marker
psi_uns

# You can the use psi_uns as input for the structural path model (standard errors will be biased)

References

Perez Alonso, A. F., Vermunt, J., Rosseel, Y., & De Roover, K. (2025, April 24). Mixture Multigroup Structural Equation Modeling for Ordinal Data. https://doi.org/10.31234/osf.io/6cr74_v1

Rosseel, Y., & Loh, W. W. (2022). A structural after measurement approach to structural equation modeling. Psychological Methods. https://doi.org/10.1037/met0000503

[ahmad]

Thank you for your response, Andres. I may send you a private message with some results.

Best,

A