Categorical and mixed-indicator models

Scope

This guide covers ordered and mixed indicators. The implementation covers lavaan’s DWLS and ULS families, one or several groups, and listwise or pairwise missingness. The supported paths have also been checked against a separate implementation.

For these models, semTests uses lavaan’s unscaled statistic and the eigenvalues stored in lavInspect(fit, "UGamma"). The Du–Bentler UG option and the continuous-data RLS statistic belong to a different setting, so test names here do not use either suffix.

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

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

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

The cuts give a reproducible teaching example. They are not a recommendation to turn continuous measurements into categories.

Overall goodness of fit

lavaan’s WLSMV shortcut fits the DWLS family and requests a mean-and-variance-adjusted test. From there, pvalues() works as usual:

fit_ordinal <- cfa(
  model, ordinal_data,
  ordered = ordered_names,
  estimator = "WLSMV"
)

pvalues(fit_ordinal, c("STD", "SB", "SS", "ALL", "PEBA4"))
#>       std_ml        sb_ml        ss_ml       all_ml     peba4_ml 
#> 5.741040e-03 1.579645e-06 2.165026e-05 1.204187e-04 8.080871e-06 
#> estimator: DWLS (WLSMV) | data: categorical | information: expected | df: 24

The footer records both the base estimator (DWLS) and the estimator requested from lavaan (WLSMV). It also records the categorical parameterization.

attr(pvalues(fit_ordinal, "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] 24
#> 
#> $nested
#> [1] FALSE
#> 
#> $method
#> [1] NA
#> 
#> $parameterization
#> [1] "delta"
#> 
#> $spectrum_source
#> [1] "lavaan UGamma"
#> 
#> $requested_tests
#> [1] "PEBA4"
#> 
#> $statistic
#> PEBA4 
#>  "ml" 
#> 
#> $gamma
#>    PEBA4 
#> "biased"

The SB and SS results are a quick check against lavaan’s public robust corrections. The eigenvalue methods use more of the reference distribution, so they need not give the same p-value.

Estimator families

The accepted shortcuts include:

DWLS, WLSM, WLSMV, and WLSMVS
ULS, ULSM, and ULSMV

Full WLS/ADF is refused. Its weight already produces an identity spectrum, so every correction reduces to the ordinary chi-square. There is no extra correction for semTests to make.

Mixed indicators

Only variables named in ordered are treated as categorical. The remaining indicators stay continuous:

mixed_names <- paste0("x", 1:6)
mixed_data <- HolzingerSwineford1939
mixed_data[mixed_names] <- lapply(
  mixed_data[mixed_names],
  function(x) {
    ordered(cut(
      x,
      breaks = c(-Inf, quantile(x, c(.33, .67)), Inf),
      include.lowest = TRUE
    ))
  }
)

fit_mixed <- cfa(
  model, mixed_data,
  ordered = mixed_names,
  estimator = "WLSMV"
)
pvalues(fit_mixed, "PEBA4")
#>     peba4_ml 
#> 2.836827e-06 
#> estimator: DWLS (WLSMV) | data: categorical | information: expected | df: 24

Mixed models combine thresholds with the usual covariance and mean information. Lavaan has already done that bookkeeping, and semTests reads the resulting spectrum directly.

Pairwise missingness

With pairwise estimation, lavaan uses the observations available for each statistic:

pairwise_data <- ordinal_data
set.seed(314)
for (variable in ordered_names) {
  pairwise_data[sample(nrow(pairwise_data), 35), variable] <- NA
}

fit_pairwise <- cfa(
  model, pairwise_data,
  ordered = ordered_names,
  estimator = "WLSMV",
  missing = "pairwise"
)
pvalues(fit_pairwise, c("SB", "SS", "PEBA4"))
#>        sb_ml        ss_ml     peba4_ml 
#> 0.0006333295 0.0020325368 0.0010972907 
#> estimator: DWLS (WLSMV) | data: categorical | information: expected | df: 24

semTests follows lavaan’s pairwise sample-statistic and UGamma calculation. Pairwise deletion brings its own assumptions about missingness, so report and defend that choice separately. It is not FIML.

Nested categorical models

For categorical models, nested comparison uses Satorra’s 2000 construction and the delta restriction map:

constrained <- "
  visual  =~ x1 + a*x2 + a*x3
  textual =~ x4 + b*x5 + b*x6
  speed   =~ x7 + x8 + x9
"

m1 <- cfa(
  model, ordinal_data,
  ordered = ordered_names, estimator = "WLSMV"
)
m0 <- cfa(
  constrained, ordinal_data,
  ordered = ordered_names, estimator = "WLSMV"
)

pvalues_nested(
  m0, m1,
  method = "2000", A.method = "delta",
  tests = c("SB", "SS", "PALL")
)
#>     sb_ml     ss_ml   pall_ml 
#> 0.7312051 0.6986248 0.7275842 
#> estimator: DWLS (WLSMV) | data: categorical | information: expected | df: 2 | nested (method 2000, A.method delta)

Both fits must use the same requested estimator, parameterization, variables, groups, raw data, and missingness mask. Method 2001 and the exact restriction map are unavailable for this family.

Practical checks

Before reading too much into a p-value, check the basics:

Verify convergence and adequate category counts.
Look for sparse categories and, when relevant, small pairwise counts.
Confirm that the two nested fits use identical data and options.
Compare SB and SS with lavaan as a quick numerical check.
Report the estimator and any use of pairwise missingness.