---
title: "Ordinal indicators and measurement invariance"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Ordinal indicators and measurement invariance}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r, include = FALSE}
knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>",
  message = FALSE,
  warning = FALSE
)
```

## 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.

```{r setup}
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

```{r fit}
fit <- cfa(model, ordinal_data, ordered = items, estimator = "WLSMV")
pvalues(fit, c("STD", "SB", "SS", "ALL", "PEBA4"))
```

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

```{r lavaan}
as.numeric(fitMeasures(fit, "pvalue.scaled"))
```

## Mixed indicators

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

```{r mixed}
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)
```

## 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.

```{r pairwise}
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"))
```

## 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:

```{r invariance-fits}
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:

```{r invariance-tests}
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
loadings_test
```

With `PALL`, equal thresholds give *p* =
`r sprintf("%.3f", thresholds_test[["pall_ml"]])` and adding equal loadings gives
*p* = `r sprintf("%.3f", loadings_test[["pall_ml"]])`, so neither restriction is
rejected. The configural model itself does not fit exactly:

```{r configural}
pvalues(configural)
```

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.

## 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>