Robust & Scaled Estimation

Overview

Ordinary maximum likelihood (ML) for continuous data assumes multivariate normality. When that assumption is violated — skewed or heavy-tailed items, excess kurtosis — the ML point estimates remain consistent, but the model chi-square is inflated and the standard errors are biased. Robust maximum-likelihood estimators keep the ML estimates and correct the inference: they replace the naive standard errors with sandwich (robust) standard errors and rescale the chi-square so that its reference distribution is approximately correct.

Added in version 0.16 and validated against lavaan, psychonetrics implements the full robust-ML family. They are available for the lvm and varcov families and are selected either through the model constructors (estimator = …) or afterwards with setestimator().

The robust estimators

All four also report robust RMSEA, CFI and TLI. "MLR" is the most general: it is the only robust estimator here that handles missing data, using FIML with the Savalei (2010) FIML-corrected robust fit indices (matching lavaan's estimator = "MLR", missing = "fiml"). The Satorra–Bentler family (MLM/MLMV/MLMVS) is complete-data only.

When to use them

• Continuous, non-normal data. If items show meaningful skew or kurtosis, prefer a robust estimator over plain "ML". "MLR" is a safe, widely-reported default (it is the Mplus MLR / lavaan estimator = "MLR").
• Missing data + non-normality. Use "MLR", which combines FIML with robust corrections.
• Reporting a single scaled statistic. "MLMV" (scaled-and-shifted) gives a chi-square whose mean and variance both match the reference distribution; it pairs with the "satorra.2000" difference test.
• Genuinely normal data. Plain "ML" is fully efficient and the robust corrections add nothing; there is no harm in checking with "MLR", though — the scaling factor will be close to 1.

For ordinal data, the categorical estimators (WLSMV / DWLS with robust standard errors) are the appropriate choice instead; see the LVM page.

Example: MLR on HolzingerSwineford1939

We fit the classic three-factor Holzinger–Swineford CFA with robust ML. The only change from a plain ML fit is estimator = "MLR":

library("psychonetrics")
library("lavaan")
library("dplyr")
data(HolzingerSwineford1939)

# Three-factor measurement model (visual / textual / speed):
Lambda <- matrix(0, 9, 3)
Lambda[1:3, 1] <- 1   # visual:  x1, x2, x3
Lambda[4:6, 2] <- 1   # textual: x4, x5, x6
Lambda[7:9, 3] <- 1   # speed:   x7, x8, x9

mod_mlr <- lvm(HolzingerSwineford1939, lambda = Lambda,
               vars = paste0("x", 1:9),
               latents = c("visual", "textual", "speed"),
               identification = "variance",
               estimator = "MLR") %>% runmodel

mod_mlr %>% fit

Every one of these quantities matches lavaan::cfa(HS.model, data = HolzingerSwineford1939, std.lv = TRUE, estimator = "MLR") to three decimals (lavaan reports scaled χ² = 87.132, scaling factor 0.979, robust CFI 0.930, robust TLI 0.895, robust RMSEA 0.092).

The robust correction also widens the standard errors, reflecting the non-normality in these items:

# Compare naive (ML) and robust (MLR) standard errors:
parameters(mod_mlr)

The estimates are identical to a plain ML fit; only the standard errors (and hence p-values and confidence intervals) change.

Switching estimators with setestimator()

You can also start from a model fitted with the default estimator and switch afterwards. Because the robust standard errors and scaled statistics are computed from the raw data, build the original model with storedata = TRUE so the data travel with the model (the robust constructors do this automatically):

mod <- lvm(HolzingerSwineford1939, lambda = Lambda,
           vars = paste0("x", 1:9),
           latents = c("visual", "textual", "speed"),
           identification = "variance",
           storedata = TRUE) %>% runmodel        # default ML fit

# Switch to Satorra-Bentler robust ML and re-run:
mod_mlm <- mod %>% setestimator("MLM") %>% runmodel

mod_mlm %>% fit    # chisq.scaled = 80.875, scaling factor = 1.055

setestimator() resets the fit so the model is re-estimated on the next runmodel(), and propagates the estimator to the internal baseline and saturated reference models.

Scaled chi-square difference tests

Naive chi-square difference tests are not valid under robust estimation, because the difference of two scaled statistics is not itself scaled correctly. For nested models fitted with a robust or scaled estimator, compare() therefore reports a Satorra–Bentler-family scaled difference test in two extra columns, Chisq_diff_scaled and p_value_scaled. The scaled.test.method argument chooses the variant:

# Two nested MLR models: configural vs equal loadings across schools
mod_config <- lvm(HolzingerSwineford1939, lambda = Lambda,
                  vars = paste0("x", 1:9),
                  latents = c("visual", "textual", "speed"),
                  identification = "variance",
                  groups = "school", estimator = "MLR") %>% runmodel

mod_weak <- mod_config %>% groupequal("lambda") %>% runmodel

compare(configural = mod_config, weak = mod_weak)

The scaled difference (6.57, p = 0.36) reproduces lavaan::lavTestLRT(fit_config, fit_weak) (6.567, p = 0.363). Output for non-robust estimators is unchanged: the extra columns simply do not appear.

Summary

• Four robust estimators — "MLM", "MLMV", "MLMVS", "MLR" — keep the ML point estimates and correct standard errors, the test statistic, and fit indices for non-normality.
• "MLR" additionally supports missing data via FIML; the Satorra–Bentler family is complete-data only.
• Select via the constructor (estimator = "MLR") or with setestimator() (build with storedata = TRUE).
• compare() reports a scaled chi-square difference test (scaled.test.method) reproducing lavaan's lavTestLRT().
• All quantities were validated against lavaan. For converting models to and from lavaan, see the interoperability page.