Measurement invariance with ordinal indicators

Scope

Measurement invariance asks whether a construct is measured the same way across groups. The usual strategy fits a sequence of increasingly constrained models and compares each model with the one before it.

Ordinal indicators need a little care here. Thresholds enter the measurement model alongside factor loadings, and the order of the invariance tests matters. Following Wu and Estabrook (2016), we test threshold equality first and loading equality second. The example is deliberately small, but it gives pvalues_nested() a realistic categorical job to do.

library("semTests")
library("lavaan")

model <- "
  visual  =~ x1 + x2 + x3
  textual =~ x4 + x5 + x6
  speed   =~ x7 + x8 + x9
"

ordered_names <- paste0("x", 1:9)
mi_data <- HolzingerSwineford1939
mi_data[ordered_names] <- lapply(
  mi_data[ordered_names],
  function(x) {
    ordered(cut(
      x,
      breaks = c(-Inf, quantile(x, c(.25, .50, .75)), Inf),
      include.lowest = TRUE
    ))
  }
)

The grouping variable is school, which has two levels. Cutting at the quartiles gives every indicator four ordered categories and three thresholds. This artificial split is only here to make the example reproducible. Ordinal status should come from the measurement design in a real study.

The configural model

The configural model uses the same factor pattern in both schools while leaving its measurement parameters free across schools. We can start by checking its overall fit:

configural <- cfa(
  model, mi_data,
  ordered = ordered_names,
  group = "school", estimator = "WLSMV"
)

pvalues(configural, c("SB", "SS", "ALL", "PEBA4"))
#>        sb_ml        ss_ml       all_ml     peba4_ml 
#> 2.820266e-06 8.043171e-05 3.211993e-04 1.638551e-05 
#> estimator: DWLS (WLSMV) | data: categorical | information: expected | df: 48

Thresholds first

The next model makes the thresholds equal across schools. It is more constrained than the configural model, so it goes in the first argument:

equal_thresholds <- cfa(
  model, mi_data,
  ordered = ordered_names,
  group = "school", group.equal = "thresholds", estimator = "WLSMV"
)

pvalues_nested(
  equal_thresholds, configural,
  tests = c("SB", "SS", "ALL", "PEBA4")
)
#>     sb_ml     ss_ml    all_ml  peba4_ml 
#> 0.7167006 0.7119125 0.7119828 0.7156707 
#> estimator: DWLS (WLSMV) | data: categorical | information: expected | df: 9 | nested (method 2000, A.method delta)

The test asks whether the equal-threshold restrictions create a noticeable loss of fit. A large p-value means the data give no clear reason to reject those restrictions. It does not prove that the thresholds are identical.

Then the loadings

Now retain equal thresholds and add equality constraints on the loadings:

equal_loadings <- cfa(
  model, mi_data,
  ordered = ordered_names,
  group = "school",
  group.equal = c("thresholds", "loadings"),
  estimator = "WLSMV"
)

pvalues_nested(
  equal_loadings, equal_thresholds,
  tests = c("SB", "SS", "ALL", "PEBA4")
)
#>     sb_ml     ss_ml    all_ml  peba4_ml 
#> 0.9971644 0.9899814 0.9966135 0.9970466 
#> estimator: DWLS (WLSMV) | data: categorical | information: expected | df: 6 | nested (method 2000, A.method delta)

This second comparison isolates the extra loading restrictions given equal thresholds. A small p-value would count against that step. A large one says the step did not noticeably worsen fit in this sample. In a real invariance study, the model-identification choices, effect sizes, and the purpose of the group comparison deserve attention too.

Both comparisons use Satorra’s 2000 reduction and the delta restriction map. Those are the categorical nested options supported by semTests.

attr(
  pvalues_nested(
    equal_loadings, equal_thresholds,
    tests = "PEBA4"
  ),
  "semtests"
)
#> $estimator
#> [1] "DWLS"
#> 
#> $estimator_requested
#> [1] "WLSMV"
#> 
#> $lavaan_test
#> [1] "standard"       "scaled.shifted"
#> 
#> $information
#> [1] "expected" "expected"
#> 
#> $data_type
#> [1] "categorical"
#> 
#> $missing
#> [1] "listwise"
#> 
#> $df
#> [1] 6
#> 
#> $nested
#> [1] TRUE
#> 
#> $method
#> [1] "2000"
#> 
#> $parameterization
#> [1] "delta"
#> 
#> $spectrum_source
#> [1] "lavaan UGamma"
#> 
#> $A.method
#> [1] "delta"
#> 
#> $requested_tests
#> [1] "PEBA4"
#> 
#> $statistic
#> PEBA4 
#>  "ml" 
#> 
#> $gamma
#>    PEBA4 
#> "biased"

The result footer records the categorical estimator, the reduction method, and the restriction map. That makes it easy to tell later which test you ran.

References

Satorra, A. (2000). Scaled and adjusted restricted tests in multi-sample analysis of moment structures. In R. D. H. Heijmans, D. S. G. Pollock, & A. Satorra (Eds.), Innovations in Multivariate Statistical Analysis (pp. 233–247). Kluwer Academic. https://doi.org/10.1007/978-1-4615-4603-0_17

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