Missing data with FIML

Scope

semTests supports continuous data with missing values fitted by full-information maximum likelihood (missing = "fiml"), with one or several groups. The statistic is the likelihood-ratio statistic and \(\Gamma\) is the usual biased estimate, so the UG and RLS options do not apply. Observed predictors are supported when they are modeled jointly with the other variables (fixed.x = FALSE and conditional.x = FALSE). Fixed or conditional predictors are refused.

A model with missing values

To have data that are missing at random, we delete values of x1 and x5 with probabilities that depend on the always-observed x9:

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

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

fiml_data <- HolzingerSwineford1939
set.seed(20260717)
driver <- as.numeric(scale(fiml_data$x9))
missing_x1 <- runif(nrow(fiml_data)) < plogis(qlogis(.20) + .8 * driver)
missing_x5 <- runif(nrow(fiml_data)) < plogis(qlogis(.20) - .8 * driver)
fiml_data$x1[missing_x1] <- NA
fiml_data$x5[missing_x5] <- NA
colMeans(is.na(fiml_data[c("x1", "x5")]))
#>        x1        x5 
#> 0.2624585 0.2059801

The fit and the call to pvalues() look as they do for complete data:

fit <- cfa(model, fiml_data, missing = "fiml")
pvalues(fit, c("SB", "SS", "PEBA4"))
#>        sb_ml        ss_ml     peba4_ml 
#> 2.442470e-08 1.018704e-07 6.558237e-08 
#> estimator: ML (FIML) | data: continuous | information: observed | df: 24 | FIML convention: observed

Observed information

semTests uses observed information for both the saturated and the fitted model. Expected information, as SEM software computes it, ignores the missingness mechanism and is inconsistent when data are missing at random but not completely at random (Kenward and Molenberghs, 1998). Savalei (2010) shows the consequences for SEM, and Savalei and Rosseel (2022) describe how lavaan’s options handle the choice. We also compared the conventions by simulation, for the interval coverage of a loading and a factor covariance, with normal and \(t_5\) data, complete, MCAR, and MAR data, and sample sizes from 150 to 1200. Observed information with a sandwich covariance came closest to the nominal 95%, most clearly for heavy-tailed data and misspecified models.

lavaan’s information option therefore does not change the semTests result. It does change lavaan’s own robust tests:

fit_expected <- cfa(model, fiml_data, missing = "fiml", information = "expected")
lavaan_sb <- function(fit) {
  lavTest(fit, test = "satorra.bentler")$satorra.bentler$pvalue
}
rbind(
  semTests = c(
    observed = pvalues(fit, "SB")[[1]],
    expected = pvalues(fit_expected, "SB")[[1]]
  ),
  lavaan = c(lavaan_sb(fit), lavaan_sb(fit_expected))
)
#>              observed     expected
#> semTests 2.442470e-08 2.442470e-08
#> lavaan   1.048556e-10 6.647741e-08

The printed footer reports the fit’s information option, and FIML convention: observed reports what semTests used.

Observed predictors

A predictor whose distribution is part of the model can have missing values too. Here ageyr predicts the visual factor:

predictor_data <- fiml_data
predictor_data$ageyr[seq(7, nrow(predictor_data), by = 23)] <- NA
fit_predictor <- cfa("visual =~ x1 + x2 + x3
                      visual ~ ageyr",
  predictor_data,
  missing = "fiml", fixed.x = FALSE, conditional.x = FALSE
)
pvalues(fit_predictor, c("SB", "PALL"))
#>     sb_ml   pall_ml 
#> 0.6776279 0.6775772 
#> estimator: ML (FIML) | data: continuous | information: observed | df: 2 | FIML convention: observed

The model has 2 degrees of freedom, so it cannot use the default PEBA4. Treating ageyr as random is a modeling decision. If it is a fixed design variable, a conditional analysis is the right one, and semTests does not support it yet.

Nested models: is growth linear?

Demo.growth has four repeated measurements, t1 to t4. We delete values at the last two waves with probabilities that depend on t1:

growth_data <- Demo.growth
set.seed(20260718)
driver <- as.numeric(scale(growth_data$t1))
for (wave in c("t3", "t4")) {
  dropped <- runif(nrow(growth_data)) < plogis(qlogis(.20) + .8 * driver)
  growth_data[dropped, wave] <- NA
}
colMeans(is.na(growth_data[paste0("t", 1:4)]))
#>    t1    t2    t3    t4 
#> 0.000 0.000 0.195 0.190

The linear model fixes the slope loadings at 0, 1, 2, and 3. The latent-basis model frees the last two and lets the data choose the shape of change, so the linear model is the special case and goes first:

linear <- growth("i =~ 1*t1 + 1*t2 + 1*t3 + 1*t4
                  s =~ 0*t1 + 1*t2 + 2*t3 + 3*t4",
  growth_data,
  missing = "fiml"
)
basis <- growth("i =~ 1*t1 + 1*t2 + 1*t3 + 1*t4
                 s =~ 0*t1 + 1*t2 + t3 + t4",
  growth_data,
  missing = "fiml"
)
growth_test <- pvalues_nested(linear, basis, tests = c("SB", "SS", "PALL"))
growth_test
#>     sb_ml     ss_ml   pall_ml 
#> 0.5683004 0.5672936 0.5680640 
#> estimator: ML (FIML) | data: continuous | information: observed | df: 2 | FIML convention: observed | nested (method 2000, A.method weighted)
coef(basis)[c("s=~t3", "s=~t4")]
#>    s=~t3    s=~t4 
#> 1.924029 2.869086

The free loadings are estimated at 1.92 and 2.87, close to 2 and 3, and the test gives p = 0.568. The data give no reason to prefer the latent-basis model.

The default restriction matrix, A.method = "weighted", handles this pair, where the larger model frees parameters that the smaller one fixes. A.method = "exact" needs the same free parameters in both models, as with equality constraints, and refuses it:

pvalues_nested(linear, basis, A.method = "exact")
#> Error:
#> ! The exact nested restriction requires the same underlying full free-parameter set in both models. Use `A.method = "weighted"` for models represented by different parameter sets.

When both apply, the "exact" p-value shifts slightly with the factor identification, for example std.lv = TRUE against marker loadings, while the "weighted" p-value does not. lavaan’s lavTestLRT(basis, linear, method = "satorra.2000") on MLR fits uses A.method = "delta" and follows the fits’ information options, so it need not agree with semTests.

Several groups

Grouped models use the same calls. Here the factor loadings are held equal across the two schools:

configural <- cfa(model, fiml_data, group = "school", missing = "fiml")
equal_loadings <- cfa(model, fiml_data,
  group = "school", group.equal = "loadings", missing = "fiml"
)
pvalues_nested(equal_loadings, configural)
#>   pall_ml 
#> 0.2231991 
#> estimator: ML (FIML) | data: continuous | information: observed | df: 6 | FIML convention: observed | nested (method 2000, A.method weighted)

The groups may differ in size and missingness pattern, but the two fits must use the same observations and missingness pattern within each group.

References

Kenward, M. G., & Molenberghs, G. (1998). Likelihood based frequentist inference when data are missing at random. Statistical Science, 13(3), 236–247. https://doi.org/10.1214/ss/1028905886

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

Savalei, V. (2010). Expected versus observed information in SEM with incomplete normal and nonnormal data. Psychological Methods, 15(4), 352–367. https://doi.org/10.1037/a0020143

Savalei, V., & Rosseel, Y. (2022). Computational options for standard errors and test statistics with incomplete normal and nonnormal data in SEM. Structural Equation Modeling, 29(2), 163–181. https://doi.org/10.1080/10705511.2021.1877548