Ordinal indicators and measurement invariance

Scope

semTests supports ordered and mixed indicators fitted with lavaan’s DWLS family (DWLS, WLSM, WLSMV, WLSMVS) and ULS family (ULS, ULSM, ULSMV), with one or several groups and listwise or pairwise deletion. The eigenvalues come from lavInspect(fit, "UGamma") and the statistic is lavaan’s unscaled one. The UG and RLS options apply to continuous data only.

Full WLS is refused. Its weight matrix is the inverse of \(\Gamma\), so every eigenvalue equals 1 and every correction would return the ordinary chi-square p-value.

Data

The examples cut the nine ability scores of HolzingerSwineford1939 at their quartiles, giving four ordered categories per indicator.

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

model <- "
  visual  =~ x1 + x2 + x3
  textual =~ x4 + x5 + x6
  speed   =~ x7 + x8 + x9
"
quartiles <- function(x) {
  ordered(cut(x, c(-Inf, quantile(x, c(.25, .50, .75)), Inf)))
}

items <- paste0("x", 1:9)
ordinal_data <- HolzingerSwineford1939
ordinal_data[items] <- lapply(ordinal_data[items], quartiles)

Overall fit

fit <- cfa(model, ordinal_data, ordered = items, estimator = "WLSMV")
pvalues(fit, c("STD", "SB", "SS", "ALL", "PEBA4"))
#>       std_ml        sb_ml        ss_ml       all_ml     peba4_ml 
#> 8.987813e-04 6.279085e-09 5.051364e-07 1.468198e-05 1.427912e-07 
#> estimator: DWLS (WLSMV) | data: categorical | information: expected | df: 24

SS is the mean-and-variance adjusted test that WLSMV reports, and SB is the mean-adjusted test that WLSM reports:

as.numeric(fitMeasures(fit, "pvalue.scaled"))
#> [1] 5.051364e-07

Mixed indicators

Only the variables named in ordered are treated as ordinal. Here the speed indicators stay continuous:

mixed_data <- HolzingerSwineford1939
mixed_data[items[1:6]] <- lapply(mixed_data[items[1:6]], quartiles)
fit_mixed <- cfa(model, mixed_data, ordered = items[1:6], estimator = "WLSMV")
pvalues(fit_mixed)
#>     peba4_ml 
#> 3.039647e-07 
#> estimator: DWLS (WLSMV) | data: categorical | information: expected | df: 24

Pairwise deletion

With missing = "pairwise", lavaan estimates each threshold and polychoric correlation from the cases observed on the variables involved, and semTests uses the resulting sample statistics and UGamma. Pairwise deletion is consistent when the data are missing completely at random.

pairwise_data <- ordinal_data
set.seed(314)
for (item in items) {
  pairwise_data[sample(nrow(pairwise_data), 35), item] <- NA
}
fit_pairwise <- cfa(model, pairwise_data,
  ordered = items, estimator = "WLSMV", missing = "pairwise"
)
pvalues(fit_pairwise, c("SB", "SS", "PEBA4"))
#>        sb_ml        ss_ml     peba4_ml 
#> 3.952479e-07 7.174855e-06 2.527975e-06 
#> estimator: DWLS (WLSMV) | data: categorical | information: expected | df: 24

Measurement invariance

Measurement invariance asks whether the indicators measure the factors the same way in each group. The usual approach fits a sequence of increasingly restricted models and compares each with the previous one. For ordinal indicators, Wu and Estabrook (2016) recommend restricting the thresholds before the loadings. We compare the two schools in the data:

configural <- cfa(model, ordinal_data,
  ordered = items, group = "school", estimator = "WLSMV"
)
equal_thresholds <- cfa(model, ordinal_data,
  ordered = items, group = "school", estimator = "WLSMV",
  group.equal = "thresholds"
)
equal_loadings <- cfa(model, ordinal_data,
  ordered = items, group = "school", estimator = "WLSMV",
  group.equal = c("thresholds", "loadings")
)

When the thresholds are held equal, lavaan frees the intercepts and scale factors of the latent responses in the second school, which is the identification of Wu and Estabrook. The 27 threshold equalities therefore remove 9 degrees of freedom. With three categories there would be 18 equalities and 18 freed parameters, and the threshold model would be equivalent to the configural model.

Each nested comparison takes the more restricted model first:

thresholds_test <- pvalues_nested(equal_thresholds, configural,
  tests = c("SB", "SS", "PALL")
)
loadings_test <- pvalues_nested(equal_loadings, equal_thresholds,
  tests = c("SB", "SS", "PALL")
)
thresholds_test
#>     sb_ml     ss_ml   pall_ml 
#> 0.7167006 0.7119125 0.7155130 
#> estimator: DWLS (WLSMV) | data: categorical | information: expected | df: 9 | nested (method 2000, A.method weighted)
loadings_test
#>     sb_ml     ss_ml   pall_ml 
#> 0.9971645 0.9899861 0.9970424 
#> estimator: DWLS (WLSMV) | data: categorical | information: expected | df: 6 | nested (method 2000, A.method weighted)

With PALL, equal thresholds give p = 0.716 and adding equal loadings gives p = 0.997, so neither restriction is rejected. The configural model itself does not fit exactly:

pvalues(configural)
#>     peba4_ml 
#> 1.638551e-05 
#> estimator: DWLS (WLSMV) | data: categorical | information: expected | df: 48

The nested tests assume that the larger model is correct. vignette("score-criterion", package = "semTests") describes restriction tests that allow it to be misspecified.

Nested tests for ordinal data use Satorra’s (2000) construction with the default A.method = "weighted", or lavaan’s "delta". The "exact" restriction matrix is not available here. Both fits must use the same estimator, parameterization, variables, groups, data, and missingness pattern.