Measurement invariance: comparing groups with multigroup CFA

Module: Factor analysis
Tranche: core (a first course)
Before this: Factor analysis II: confirmatory models, omega, and IRT as factor analysis
Builds toward: Structural equation modeling with lavaan
Starts threads: Comparing group means needs scalar invariance (→ Structural equation modeling with lavaan)
Picks up threads: Choose the estimator from the data: categories, shape, missingness (from Factor analysis II: confirmatory models, omega, and IRT as factor analysis)
IRW tables: neurodegenerative_huizinga_2019_dass, alexandrowicz_2018_cesd
From EDUC 252: n/a
Status: draft
Authorship: Written largely by Claude (Anthropic), from Ben Domingue’s EDUC 252 materials and under Ben’s direction. Ben reviews each lesson before it is marked done.
Feedback: This page improved thanks to feedback from Lijin Zhang. All remaining errors are the responsibility of the page’s creator (BD).
Comments: open an issue on GitHub about this lesson.

What this is for

In a Dutch community sample, women report a little more anxiety than men on the DASS-21 (0.16 SD, as we’ll see). Is that a difference in anxiety, or a difference in how the items work for women and men? The scores alone can’t say. Factor analysis II fitted a factor model to one sample. This lesson fits it to two groups at once and asks, one constraint at a time, whether the same model holds in both: measurement invariance. Where it holds, the groups’ factor means can be compared. Where an item breaks it, we find the item and decide what to do about it. We start with a scale where invariance holds, the DASS-21 across women and men, and then turn to one where a single item doesn’t: crying on the CES-D.

Goals

By the end of this lesson you should be able to:

  • fit a multigroup CFA in lavaan and say how the two groups are put on one scale;
  • climb the ladder (configural, metric, scalar) with scaled difference tests and changes in CFI, and say which comparison each step supports;
  • find the constraint that fails with score tests, and fit and report a partial invariance model;
  • say what one unequal intercept does to a comparison of sum scores and of factor means;
  • say what the choice of estimator for a multigroup fit assumes, and what ordered categories change.

Core ideas

One model, two groups

Recall from Factor analysis II: a confirmatory model fixes most loadings at zero before seeing the data, and it is identified once each factor’s scale is set, by fixing one loading or the factor’s variance. Nested models are compared by a (scaled) difference in χ².

A multigroup CFA fits the model in several groups at once, with any parameter free in each group or held equal across them (Jöreskog, 1971). With intercepts, the model for item \(i\) in group \(g\) is \(x_i = \nu_{ig} + \lambda_{ig}\eta + \varepsilon_i\). Measurement invariance means that responses depend on the factor and not on the group (Meredith, 1993). It is tested in steps, each adding constraints to the one before (Vandenberg & Lance, 2000; Millsap, 2011; Putnick & Bornstein, 2016):

  • configural: the same items load on the same factors in every group, every parameter free;
  • metric: the loadings \(\lambda_i\) are equal, so a unit of the factor means the same in each group;
  • scalar: the intercepts \(\nu_i\) (for ordered items, the thresholds) are equal too, so respondents at the same factor level have the same expected response.

Meredith adds a fourth step, strict invariance, with equal unique variances as well.

Why these comparisons at these steps? A factor has no unit and no origin of its own. In the configural model each group’s factor is scaled separately, so nothing links them. Once the loadings are equal, the items carry one unit across groups, and group B’s factor variance, and its correlations and regressions, can be estimated relative to group A’s. Once the intercepts are equal too, the items carry one origin, and group B’s factor mean can be estimated relative to A’s; lavaan frees it at that step. Comparing correlations or regressions needs metric invariance; comparing means needs scalar.

In item response language, an unequal loading is non-uniform differential item functioning and an unequal intercept is uniform DIF; Differential item functioning takes up the same question item by item, matching on an observed score instead of a factor.

What an unequal intercept does

Try it: one item that works differently

Four items measure one factor in groups A and B (loadings 0.7, intercepts 0, factor SD 1). Item 4 can behave differently in group B.

With both item-4 controls at 0 the lines coincide, and the sum-score gap (0.45 SD) is the latent 0.5 SD shrunk by unreliability. Set the true gap to 0 and item 4’s intercept to 0.3: the sum scores differ by 0.10 SD between groups that don’t differ at all. A loading difference makes the lines cross, so the item favours one group at low factor levels and the other at high ones.

What does the factor model make of the same item? If we impose scalar invariance anyway, the model has to explain item 4’s higher mean in group B with the one thing it has left free, group B’s factor mean. The shift is spread over every item, and the factor mean absorbs part of it. One response is partial invariance (Byrne, Shavelson & Muthén, 1989): free the offending item’s intercept and keep the others equal, so that the comparison of means rests on the items that do work alike.

Try it: what the scalar constraint does to the factor mean

Now \(k\) items load 0.7 in both groups, and the last item’s intercept is higher in group B by the amount you set. The widget computes the large-sample answer: group B’s estimated factor mean under full and under partial scalar invariance, and how often the test of the scalar step (on \(k - 1\) degrees of freedom) would catch the shifted item.

At the start, the groups don’t differ, yet the scalar model puts group B 0.107 SD above A. With equal loadings the error is the shift divided by \(k\lambda\) (problem 1), here \(0.3/(4 \times 0.7)\). The partial model recovers the truth. With 300 respondents per group, the scalar step fails 97% of the time; halve the shift and the error halves while the chance of catching it falls to 44%. Raise \(k\) to 12 and the error falls to 0.036: more invariant items dilute one shifted one. Move the true mean: the error stays the same size, whatever the groups’ real difference.

Testing the steps

Each step is nested in the one before, so a χ² difference tests it, scaled when the estimator is robust (Satorra & Bentler, 2001). As with any χ² test, a large enough sample rejects a step for differences too small to matter. So changes in fit indices are read alongside: Cheung and Rensvold (2002) proposed treating a drop in CFI of more than 0.01 as a failed step, and Chen (2007) paired it with changes in RMSEA and SRMR; Putnick and Bornstein (2016) survey how these criteria are used and questioned.

When a step fails, which constraint is responsible? A score test (lavTestScore() in lavaan) says, for each equality constraint, how much the χ² would fall if that constraint alone were released. The largest is the first suspect. But a scalar model with many items holds dozens of constraints, and a few small \(p\)-values turn up by chance.

Try it: how many flags by chance?

At 39 constraints and α = 0.01, a flag turns up by chance about a third of the time (0.32); at α = 0.05, 86% of the time. Score tests on one model aren’t independent, so these are rough guides, but they set the baseline for reading the largest test.

Choosing an estimator for multigroup fits

Recall from Factor analysis II: each estimator rests on its own assumptions about the responses. ML treats items as continuous and roughly normal; MLR keeps that model but makes standard errors and tests robust to skew; WLSMV treats items as ordered categories, through thresholds and polychoric correlations; FIML, or the graded model, uses every observed response when some are missing.

In a multigroup fit those assumptions apply group by group, and three things change. With ordered items each group gets its own thresholds, so every category of every item needs responses in each group; a category that is nearly empty in one group is a problem even when it isn’t in the pooled data (Rutkowski, Svetina & Liaw, 2019, who also show that collapsing categories changes fit). Which parameters must be fixed at each step to keep ordinal models identified and nested is subtler than with continuous items (Wu & Estabrook, 2016). And the difference tests behave differently across estimators: in simulations by Sass, Schmitt and Marsh (2014), their Type I error rates varied across ML, MLR and WLSMV, and changes in fit indices under WLSMV needed caution in misspecified models. With missing ordinal responses, WLSMV with pairwise deletion was the one method in the simulations of Chen, Wu, Garnier-Villarreal, Kite and Jia (2020) whose difference tests did not reasonably control their Type I error rate; the full-information ordinal models recovered loadings best.

With sparse categories, then, one option is MLR for the ladder, items as numbers, with an ordinal fit (on collapsed categories if need be) as a check. The two rest on different assumptions about the responses; where they agree, less rests on the choice.

Simulate

The code simulates 500 respondents in each of two groups answering six items that measure one factor, with group B’s factor mean 0.3 SD higher and item 6’s intercept 0.4 higher. It climbs the ladder with lavaan, runs the score tests on the scalar model, frees the item with the largest one, and compares group B’s estimated factor mean with the truth. The first run installs lavaan, which takes a little while; after that, each run takes a few seconds.

With the seed as written, configural to metric passes (\(p = 0.31\)) and the scalar step fails (Δχ² = 42.3 on 5 df); item 6’s intercept has by far the largest score test (40.5); and group B’s mean comes out at 0.377 under full scalar invariance and 0.319 with item 6 free (truth 0.3). Set shift <- 0 and the scalar step should pass; set n <- 150 and see how often the test still finds item 6 across a few seeds.

Download this code.

With real data

The DASS-21: women and men

Our main example returns to neurodegenerative_huizinga_2019_dass, the Dutch DASS-21 answered by a community sample recruited by Huizinga et al. (2020), with the three-factor model of Factor analysis II (Lovibond & Lovibond, 1995). Its job is to show the ladder where it holds. Download the code to run every chunk.

Code 
library(lavaan)
# Each IRW table has a landing page (itemresponsewarehouse.org/tables/<name>/)
# with a CSV download. This link is pinned to one version of the data.
dass_url <- "https://redivis.com/api/v1/tables/datapages.item_response_warehouse:v59_0.neurodegenerative_huizinga_2019_dass/rows?format=csv"
df <- read.csv(dass_url)
df$resp <- suppressWarnings(as.numeric(df$resp))   # missing responses arrive as text
wide <- tapply(df$resp, list(df$id, df$item), function(x) x[1])
dass <- as.data.frame(wide[, paste0("DASS.", 1:21)])
names(dass) <- paste0("d", 1:21)
dass$sex <- df$cov_sex[match(rownames(dass), df$id)]
table(sex = dass$sex)                  # all 1,461 respondents
sex
  1   2 
676 785
Code 
dass <- dass[complete.cases(dass), ]   # 215 respondents skipped all 21 items; no one skipped some
table(sex = dass$sex)                  # the 1,246 who answered
sex
  1   2 
577 669
Code 
# 0 = did not apply to me at all ... 3 = applied to me very much. Every item is
# worded toward more distress, so no keying is needed. Scales from the DASS
# scoring template (Lovibond & Lovibond, 1995).
dep <- paste0("d", c(3, 5, 10, 13, 16, 17, 21))
anx <- paste0("d", c(2, 4, 7, 9, 15, 19, 20))
stress <- paste0("d", c(1, 6, 8, 11, 12, 14, 18))
m3 <- paste("dep =~", paste(dep, collapse = " + "), "\n",
            "anx =~", paste(anx, collapse = " + "), "\n",
            "str =~", paste(stress, collapse = " + "))

The table’s cov_sex codes are unlabelled in the processing script. The paper reports 785 women among its 1,461 participants (its Table 1), and code 2 holds exactly 785 of the table’s 1,461 respondents, so we read 2 as women: a match of counts, not a codebook.

Which estimator? For the single-group structure we used WLSMV: four categories, piled at zero. Group by group, the top category is thin.

Code 
# The estimator choice, group by group: with ordered items, every category of
# every item needs responses in each group, since each group gets thresholds.
top <- sapply(split(dass[1:21], dass$sex), function(g) sapply(g, function(x) sum(x == 3)))
table(responses_at_3 = top)            # items by how many responses sit in the top category
responses_at_3
 1  2  3  4  5  6  7  8 11 12 
 5  6  7  7  3  8  3  1  1  1
Code 
colSums(top == 0 | top == 1)           # items per group with 0 or 1 responses at 3
1 2 
4 1

Every item’s top category holds between 1 and 12 responses in each group, and five items (four among men, one among women) hold a single one. Thresholds resting on one response per group are fragile, so we climb the ladder with MLR, items as numbers, and fit WLSMV as the check.

Code 
# Multigroup CFA by cov_sex, three-factor model, items treated as numbers (MLR).
# Configural: same pattern in both groups. Metric: loadings equal. Scalar: loadings
# and intercepts equal; lavaan then frees group 2's factor means.
cfg <- cfa(m3, data = dass, group = "sex", estimator = "MLR")
met <- cfa(m3, data = dass, group = "sex", estimator = "MLR", group.equal = "loadings")
sca <- cfa(m3, data = dass, group = "sex", estimator = "MLR",
           group.equal = c("loadings", "intercepts"))
idx_r <- c("chisq.scaled", "df", "cfi.robust", "rmsea.robust", "srmr")
round(rbind(configural = fitMeasures(cfg, idx_r), metric = fitMeasures(met, idx_r),
            scalar = fitMeasures(sca, idx_r)), 3)
           chisq.scaled  df cfi.robust rmsea.robust  srmr
configural      874.695 372      0.931        0.062 0.044
metric          881.133 390      0.932        0.060 0.048
scalar          912.661 408      0.931        0.059 0.048
Code 
lavTestLRT(cfg, met, sca)   # scaled (Satorra-Bentler) differences
Df AIC BIC Chisq Chisq diff RMSEA Df diff Pr(>Chisq)
cfg 372 29902.65 30579.51 1548.076 NA NA NA NA
met 390 29896.78 30481.34 1578.205 13.50787 0.0000000 18 0.7605629
sca 408 29885.31 30377.57 1602.735 24.58263 0.0242018 18 0.1368512

Equal loadings cost almost nothing (a scaled Δχ² of 13.5 on 18 df, \(p = 0.76\)), and neither do equal intercepts (24.6 on 18 df, \(p = 0.14\)). Robust CFI is 0.931 or 0.932 at every step, far inside the 0.01 criterion. The configural fit is modest, but the question at each step is whether the constraints make it worse, and they don’t.

Code 
# Which equality constraint strains most? One score test per constraint.
score_tests <- function(fit) {
  st <- lavTestScore(fit)$uni
  pt <- parTable(fit)
  st$constraint <- sapply(st$lhs, function(l) with(pt[pt$plabel == l, ], paste(lhs, op, rhs)))
  st[order(-st$X2.scaled), c("constraint", "X2.scaled", "df", "p.value.scaled")]
}
st <- score_tests(sca)
nrow(st)
[1] 39
Code 
head(st, 3)
constraint X2.scaled df p.value.scaled
23 d16 ~1 6.921176 1 0.0085181
2 dep =~ d10 3.764150 1 0.0523625
38 d14 ~1 2.684030 1 0.1013589

The constraint under most strain is item 16’s intercept (“I was unable to become enthusiastic about anything”, in the wording of the English DASS-21; Lovibond & Lovibond, 1995), with a score test of 6.9 on 1 df. Its \(p\)-value is 0.009, and among 39 tests one that small turns up by chance about three times in ten (\(1 - 0.991^{39} \approx 0.30\), the widget’s baseline).

Code 
# Group 2's factor means under scalar invariance, in units of group 1's factor SD.
latent_means <- function(fit, factors) {
  pe <- parameterEstimates(fit)
  mu2 <- pe[pe$op == "~1" & pe$lhs %in% factors & pe$group == 2, c("lhs", "est", "se")]
  sd1 <- sqrt(diag(lavInspect(fit, "est")[[1]]$psi))
  data.frame(factor = mu2$lhs, mean_diff = round(mu2$est, 3), se = round(mu2$se, 3),
             in_group1_SDs = round(mu2$est / sd1[mu2$lhs], 2))
}
latent_means(sca, c("dep", "anx", "str"))
factor mean_diff se in_group1_SDs
dep dep 0.036 0.024 0.09
anx anx 0.039 0.016 0.16
str str 0.056 0.020 0.18

With scalar invariance the factor means can be compared: women sit 0.16 SD above men on anxiety, 0.18 SD on stress, and 0.09 SD on depression.

Code 
# The same ladder with the items as ordered categories (WLSMV): the thresholds play
# the intercepts' part. (lavaan may warn about the parameters' covariance matrix:
# some top categories hold one response in a group.)
cfg_o <- cfa(m3, data = dass, group = "sex", ordered = TRUE)
met_o <- cfa(m3, data = dass, group = "sex", ordered = TRUE, group.equal = "loadings")
sca_o <- cfa(m3, data = dass, group = "sex", ordered = TRUE,
             group.equal = c("loadings", "thresholds", "intercepts"))
idx_o <- c("df", "cfi.scaled", "rmsea.scaled", "srmr")
round(rbind(configural = fitMeasures(cfg_o, idx_o), metric = fitMeasures(met_o, idx_o),
            scalar = fitMeasures(sca_o, idx_o)), 3)
            df cfi.scaled rmsea.scaled  srmr
configural 372      0.978        0.050 0.051
metric     390      0.987        0.037 0.053
scalar     429      0.983        0.040 0.051
Code 
# Are these three models nested? Count what each one frees beyond the loadings and
# thresholds: the default scalar model frees 21 scale factors and 3 factor means
# that the default metric model fixes, so it is not the metric model plus constraints.
sapply(list(configural = cfg_o, metric = met_o, scalar = sca_o), function(f) {
  pt <- parTable(f)
  c(free_scale_factors = sum(pt$op == "~*~" & pt$free > 0),
    free_factor_means = sum(pt$op == "~1" & pt$lhs %in% c("dep", "anx", "str") & pt$free > 0))
})
                   configural metric scalar
free_scale_factors          0      0     21
free_factor_means           0      0      3

As ordered categories, every step fits well (CFI 0.978 to 0.987, RMSEA at most 0.050). We read only the indices. lavaan’s default ordinal scalar model frees 21 scale factors and 3 factor means that its metric model fixes, so the two aren’t nested and their difference isn’t a test of the thresholds alone (problem 5). The DASS-21 passes.

The CES-D: where one item breaks

Our failure case is alexandrowicz_2018_cesd: the 20-item CES-D (Radloff, 1977), answered 0 to 3 by a sample of Austrian adults (Alexandrowicz, Jahn & Wancata, 2018). Its job is to show the ladder failing and what to do next.

Code 
cesd_url <- "https://redivis.com/api/v1/tables/datapages.item_response_warehouse_4:v8_0.alexandrowicz_2018_cesd/rows?format=csv"
df <- read.csv(cesd_url)
wide <- tapply(df$resp, list(df$id, df$item), function(x) x[1])
items <- c(sprintf("ces_%02d", c(1, 2, 3, 5, 6, 7, 9, 10, 11, 13, 14, 15, 17, 18, 19, 20)),
           sprintf("cesr%02d", c(4, 8, 12, 16)))
cesd <- as.data.frame(wide[, items])
cesd$sex <- df$cov_sex[match(rownames(cesd), df$id)]
table(sex = cesd$sex)                  # all 514 respondents
sex
  1   2 
251 263
Code 
table(items_skipped = rowSums(is.na(cesd[items])))
items_skipped
  0   1   2 
508   5   1
Code 
cesd <- cesd[complete.cases(cesd), ]
table(sex = cesd$sex)                  # complete cases
sex
  1   2 
246 262
Code 
# 0 = rarely or none of the time ... 3 = most or all of the time. The four positively
# worded items (cesr04, cesr08, cesr12, cesr16) arrive reversed, so higher = more
# depressive symptoms on every item: they correlate positively with the rest.
round(range(cor(cesd[items])[lower.tri(diag(20))]), 2)
[1] -0.02  0.62
Code 
top <- sapply(split(cesd[items], cesd$sex), function(g) sapply(g, function(x) sum(x == 3)))
colSums(top <= 1)   # items per group with 0 or 1 responses in the top category
1 2 
7 5
Code 
# Two factors: the 16 symptom items, and Radloff's positive-affect items.
pos <- sprintf("cesr%02d", c(4, 8, 12, 16))
m2 <- paste("dep =~", paste(setdiff(items, pos), collapse = " + "), "\n",
            "pa =~", paste(pos, collapse = " + "))

The paper reports 264 women among 518 respondents; the table has 263 respondents coded 2 of 514, so 2 is probably women, though nothing documents it. The four positively worded items arrive reversed: no inter-item correlation is below −0.02. Six respondents skipped an item or two, so listwise deletion costs about 1%. The model has two factors: the 16 symptom items, and the four items of Radloff’s positive-affect subscale, as Alexandrowicz et al. list them (their Table 1). The estimator: the top categories are sparser than the DASS’s (seven items among men and five among women have at most one response at 3), so MLR again, with an ordinal check.

Code 
cfg_c <- cfa(m2, data = cesd, group = "sex", estimator = "MLR")
met_c <- cfa(m2, data = cesd, group = "sex", estimator = "MLR", group.equal = "loadings")
sca_c <- cfa(m2, data = cesd, group = "sex", estimator = "MLR",
             group.equal = c("loadings", "intercepts"))
round(rbind(configural = fitMeasures(cfg_c, idx_r), metric = fitMeasures(met_c, idx_r),
            scalar = fitMeasures(sca_c, idx_r)), 3)
           chisq.scaled  df cfi.robust rmsea.robust  srmr
configural      556.332 338      0.902        0.056 0.059
metric          625.794 356      0.873        0.062 0.090
scalar          663.147 374      0.865        0.063 0.092
Code 
lavTestLRT(cfg_c, met_c, sca_c)
Df AIC BIC Chisq Chisq diff RMSEA Df diff Pr(>Chisq)
cfg_c 338 15768.11 16284.23 695.6322 NA NA NA NA
met_c 356 15855.50 16295.47 819.0141 51.30654 0.1323582 18 0.0000478
sca_c 374 15858.21 16222.04 857.7328 39.10169 0.0676033 18 0.0027614

The configural model fits modestly by the conventional cutoffs (robust CFI 0.902, below 0.95; RMSEA 0.056, under 0.06). Equal loadings don’t: a scaled Δχ² of 51.3 on 18 df (\(p < .001\)), and CFI drops by 0.029, three times the 0.01 criterion. The scalar step fails as well (39.1 on 18 df).

Code 
st_met <- score_tests(met_c)   # at the metric step: which loading?
nrow(st_met)
[1] 18
Code 
head(st_met, 3)
constraint X2.scaled df p.value.scaled
12 dep =~ ces_17 11.719253 1 0.0006186
3 dep =~ ces_05 5.807588 1 0.0159572
1 dep =~ ces_02 3.420608 1 0.0643873

The crying item’s loading has the largest score test, 11.7 on 1 df among 18, twice the next. Alexandrowicz et al. (2018) found DIF by gender on the same item with a Rasch model, which is some reassurance that it isn’t an accident of this model. We free it, then look at the intercepts.

Code 
# Partial invariance: free the crying item's loading (metric), then its intercept too
# (scalar), keeping every other constraint.
met_p <- cfa(m2, data = cesd, group = "sex", estimator = "MLR", group.equal = "loadings",
             group.partial = "dep =~ ces_17")
sca_p1 <- cfa(m2, data = cesd, group = "sex", estimator = "MLR",
              group.equal = c("loadings", "intercepts"), group.partial = "dep =~ ces_17")
head(score_tests(sca_p1), 2)   # at the scalar step, with the loading already free
constraint X2.scaled df p.value.scaled
30 ces_17 ~1 7.628321 1 0.0057459
27 ces_13 ~1 6.125482 1 0.0133247
Code 
sca_p <- cfa(m2, data = cesd, group = "sex", estimator = "MLR",
             group.equal = c("loadings", "intercepts"),
             group.partial = c("dep =~ ces_17", "ces_17 ~ 1"))
lavTestLRT(cfg_c, met_p)
Df AIC BIC Chisq Chisq diff RMSEA Df diff Pr(>Chisq)
cfg_c 338 15768.11 16284.23 695.6322 NA NA NA NA
met_p 355 15801.56 16245.76 763.0753 28.05027 0.0784414 17 0.0443543
Code 
lavTestLRT(met_p, sca_p)
Df AIC BIC Chisq Chisq diff RMSEA Df diff Pr(>Chisq)
met_p 355 15801.56 16245.76 763.0753 NA NA NA NA
sca_p 372 15793.15 16165.43 788.6632 25.46969 0.0443913 17 0.0846794
Code 
round(rbind(partial_metric = fitMeasures(met_p, idx_r), partial_scalar = fitMeasures(sca_p, idx_r)), 3)
               chisq.scaled  df cfi.robust rmsea.robust  srmr
partial_metric       584.44 355      0.892        0.058 0.077
partial_scalar       610.47 372      0.889        0.057 0.078

With the crying loading free, the metric step is borderline (28.1 on 17 df, \(p = 0.04\); CFI 0.892, a drop of 0.010). At the scalar step the crying intercept is again first (7.6). Freed as well, the partial scalar model holds against the partial metric one (25.5 on 17 df, \(p = 0.08\)).

Code 
rbind(full_scalar = latent_means(sca_c, c("dep", "pa")), partial = latent_means(sca_p, c("dep", "pa")))
factor mean_diff se in_group1_SDs
full_scalar.dep dep 0.056 0.026 0.21
full_scalar.pa pa -0.045 0.047 -0.10
partial.dep dep 0.051 0.026 0.18
partial.pa pa -0.045 0.048 -0.10
Code 
# The crying item's loading and intercept in each group (partial model):
pe <- parameterEstimates(sca_p)
pe[(pe$op == "=~" & pe$rhs == "ces_17") | (pe$op == "~1" & pe$lhs == "ces_17"),
   c("lhs", "op", "rhs", "group", "est", "se")]
lhs op rhs group est se
13 dep =~ ces_17 1 0.3248080 0.1572090
56 ces_17 ~1 1 0.0853659 0.0204093
78 dep =~ ces_17 2 1.2808958 0.2350493
121 ces_17 ~1 2 0.1829560 0.0335570

Does freeing the item change the comparison? Under full scalar invariance group 2 sits 0.21 SD above group 1 on the symptom factor; with the crying item free, 0.18 SD. Positive affect is 0.10 SD lower in group 2 either way. The crying item’s loading is 0.33 in group 1 and 1.28 in group 2 (in units of the first item’s loading), and its intercept 0.085 and 0.183: in group 2, crying goes with the other symptoms far more closely, and at the factor’s zero point it is reported more often.

Code 
# As a check with an ordinal estimator: the top categories are too sparse to keep, so
# each item is scored as any symptom (1-3) against none (0), and the scalar model is
# fitted with WLSMV. Which constraints strain most?
cesd_bin <- cesd
cesd_bin[items] <- lapply(cesd_bin[items], function(x) as.integer(x > 0))
# With binary items, the theta parameterization fixes each latent response's
# residual variance at 1 in both groups, which keeps the scalar model identified
# (Wu & Estabrook, 2016); the group-2 factor means are free.
sca_b <- cfa(m2, data = cesd_bin, group = "sex", ordered = TRUE, parameterization = "theta",
             group.equal = c("loadings", "thresholds", "intercepts"))
st_b <- lavTestScore(sca_b)$uni
pt <- parTable(sca_b)
st_b$constraint <- sapply(st_b$lhs, function(l) with(pt[pt$plabel == l, ], paste(lhs, op, rhs)))
head(st_b[order(-st_b$X2), c("constraint", "X2", "df", "p.value")], 3)
constraint X2 df p.value
12 dep =~ ces_17 33.69917 1 0.0000000
14 dep =~ ces_19 18.22343 1 0.0000196
3 dep =~ ces_05 10.92681 1 0.0009478

As a check with an ordinal estimator, on items scored as any symptom against none, the crying loading again has the largest score test (33.7).

What might a report say? One view drops the crying item, so that the scale is invariant as it stands. Another keeps it and frees its parameters, so that the comparison rests on the other 19 items while the item still counts toward each respondent’s score. A third ignores the failure because it is one item in twenty. The first changes the instrument, and its scores no longer match published norms; the third lets one item push the comparison, by 0.03 SD here; the second needs enough invariant items to carry the comparison. Which cost matters depends on what the comparison is for. Reporting it both ways, with the item and the reason named, lets a reader see how much rests on the one item: here the two give 0.21 and 0.18 SD, and the conclusion doesn’t change.

Problems

  1. The size of the error. Four items with equal loadings \(\lambda\) and unique variances \(\psi\) measure one factor in two groups, with item 4’s intercept higher by \(\delta\) in group B. (a) Show that the gap between the groups’ expected sum scores is \(4\lambda\kappa + \delta\), and express it in SDs of group A’s sum score, as the first widget does. (b) With \(k\) items, show that the scalar model’s population estimate of group B’s mean is \(\kappa + \delta/(k\lambda)\). (Hint: with the covariances fitted exactly, the scalar model’s mean part is a generalized least-squares regression of the item mean differences on \(\lambda\).)
  2. A different grouping. Split the CES-D respondents at the median of cov_age and climb the ladder with the same two-factor model. Does the crying item stand out again? Which constraints do, and what would you do about them?
  3. Significant but small. A scale is compared across two groups of 10,000 respondents each. The scalar step gives a scaled Δχ² of 180 on 20 df (\(p < .001\)) and a drop in CFI of 0.003. Is the scale invariant enough to compare means? What would you look at before deciding, and what would you report?
  4. Design. You want to compare DASS-21 scores across the education levels in cov_education. Which level of invariance do you need to compare means, and which to compare the depression–anxiety correlation? Check the group sizes first: which estimator would you use, and what would you do if one item’s intercept failed in one group?
  5. Nested or not? Using the code’s ordinal DASS fits, list what each of the three models fixes and frees (thresholds, scale factors, factor means and variances). Why is the default metric-to-scalar difference not a test of equal thresholds? Build a nested sequence as Wu and Estabrook (2016) propose (the semTools function measEq.syntax() writes the syntax), test it, and compare the verdict with the MLR ladder.
  6. Challenge (open). With thirty countries, Asparouhov and Muthén (2014) note, the scalar model tends to fit poorly and to throw up many large modification indices. They propose alignment, which estimates the countries’ factor means without requiring exact invariance. What does alignment assume in place of invariance, and how would you decide whether the means it gives are trustworthy enough to rank the countries? I don’t know of a settled answer.

Stuck, or curious about something this lesson doesn’t cover? Ask Claude about this lesson opens a new conversation that already knows what the lesson is about; add your question at the end. Claude shows a caution notice on any prompt that arrives through a link; that’s expected.

Going further

  • Meredith, W. (1993). Measurement invariance, factor analysis and factorial invariance. Psychometrika, 58(4), 525–543. https://doi.org/10.1007/BF02294825. The definitions, from weak to strict.
  • Jöreskog, K. G. (1971). Simultaneous factor analysis in several populations. Psychometrika, 36(4), 409–426. https://doi.org/10.1007/BF02291366. The multigroup model and its tests.
  • Vandenberg, R. J., & Lance, C. E. (2000). A review and synthesis of the measurement invariance literature. Organizational Research Methods, 3(1), 4–70. https://doi.org/10.1177/109442810031002. The ladder, step by step.
  • Putnick, D. L., & Bornstein, M. H. (2016). Measurement invariance conventions and reporting: The state of the art and future directions for psychological research. Developmental Review, 41, 71–90. https://doi.org/10.1016/j.dr.2016.06.004. How invariance is tested and reported in practice, with the criteria for changes in fit.
  • Byrne, B. M., Shavelson, R. J., & Muthén, B. (1989). Testing for the equivalence of factor covariance and mean structures: The issue of partial measurement invariance. Psychological Bulletin, 105(3), 456–466. https://doi.org/10.1037/0033-2909.105.3.456. Partial invariance.
  • Cheung, G. W., & Rensvold, R. B. (2002). Evaluating goodness-of-fit indexes for testing measurement invariance. Structural Equation Modeling, 9(2), 233–255. https://doi.org/10.1207/s15328007sem0902_5. Where the change-in-CFI criterion comes from.
  • Chen, F. F. (2007). Sensitivity of goodness of fit indexes to lack of measurement invariance. Structural Equation Modeling, 14(3), 464–504. https://doi.org/10.1080/10705510701301834. Which index notices which kind of non-invariance.
  • Wu, H., & Estabrook, R. (2016). Identification of confirmatory factor analysis models of different levels of invariance for ordered categorical outcomes. Psychometrika, 81(4), 1014–1045. https://doi.org/10.1007/s11336-016-9506-0. Invariance with ordinal items, identified and nested (problem 5).
  • Svetina, D., Rutkowski, L., & Rutkowski, D. (2020). Multiple-group invariance with categorical outcomes using updated guidelines: An illustration using Mplus and the lavaan/semTools packages. Structural Equation Modeling, 27(1), 111–130. https://doi.org/10.1080/10705511.2019.1602776. A worked tutorial for the ordinal ladder in R.
  • Sass, D. A., Schmitt, T. A., & Marsh, H. W. (2014). Evaluating model fit with ordered categorical data within a measurement invariance framework: A comparison of estimators. Structural Equation Modeling, 21(2), 167–180. https://doi.org/10.1080/10705511.2014.882658. ML, MLR and WLSMV at each step.
  • Chen, P.-Y., Wu, W., Garnier-Villarreal, M., Kite, B. A., & Jia, F. (2020). Testing measurement invariance with ordinal missing data: A comparison of estimators and missing data techniques. Multivariate Behavioral Research, 55(1), 87–101. https://doi.org/10.1080/00273171.2019.1608799. Missing ordinal responses in the ladder.
  • Nye, C. D., & Drasgow, F. (2011). Effect size indices for analyses of measurement equivalence: Understanding the practical importance of differences between groups. Journal of Applied Psychology, 96(5), 966–980. https://doi.org/10.1037/a0022955. An effect size for non-invariance, for when “significant” isn’t the question (problem 3).
  • Millsap, R. E. (2011). Statistical approaches to measurement invariance. Routledge. https://doi.org/10.4324/9780203821961. The full treatment, factor models and IRT together.
  • The lavaan tutorial on multiple groups: group.equal, group.partial and the tests.

Data sources

Data from the Item Response Warehouse (each table name links to its IRW page); references from IRW biblio.

For instructors

  • Source: new material. Split from Factor analysis II on 09-25 (Ben, #58), at a colleague’s suggestion that the factor-analysis lesson was too long; multigroup invariance had moved into that lesson from Differential item functioning on 09-24. The DASS example, the CES-D failure case and the first widget come from it; the estimator section is new.
  • Timing: about 45 minutes for the core ideas and widgets; the real-data section is a 30-minute lab (the whole script runs in about a minute locally).
  • Placement: after Factor analysis II, and before Fairness in testing and Differential item functioning in the first course, which return to invariance with item response models and observed-score matching. Structural equation modeling with lavaan uses it for multigroup SEM.
  • Held back: invariance item by item (Differential item functioning); occasions and treatments as groupings (Measurement invariance under treatment and life events); many groups and alignment (problem 6 only); the effect-size view of non-invariance (in Going further).
  • Data: the DASS is in the public domain (DASS FAQ, no. 3), so one item is quoted; CES-D items are named by their content, not quoted. Which cov_sex code is women rests on counts matching each paper, not on a codebook. The CES-D’s positive-affect items are identified from Alexandrowicz et al. (2018, Table 1). The ladders are fitted with MLR, items as numbers, because some top categories hold one response or none in a group; the ordinal fits are checks.
  • Solutions: held back in solutions/measurement-invariance.qmd (not published).