--- title: "Missing data with FIML" output: rmarkdown::html_vignette vignette: > %\VignetteIndexEntry{Missing data with FIML} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r, include = FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>", message = FALSE, warning = FALSE ) ``` ## 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`: ```{r data} 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")])) ``` The fit and the call to `pvalues()` look as they do for complete data: ```{r fit} fit <- cfa(model, fiml_data, missing = "fiml") pvalues(fit, c("SB", "SS", "PEBA4")) ``` ## 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: ```{r information} 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)) ) ``` 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: ```{r predictor} 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")) ``` 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`: ```{r growth-data} 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)])) ``` 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: ```{r growth} 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 coef(basis)[c("s=~t3", "s=~t4")] ``` The free loadings are estimated at `r paste(round(coef(basis)[c("s=~t3", "s=~t4")], 2), collapse = " and ")`, close to 2 and 3, and the test gives *p* = `r sprintf("%.3f", growth_test[["pall_ml"]])`. 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: ```{r exact, error = TRUE} pvalues_nested(linear, basis, A.method = "exact") ``` 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: ```{r groups} 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) ``` 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. 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. Savalei, V. (2010). Expected versus observed information in SEM with incomplete normal and nonnormal data. *Psychological Methods*, 15(4), 352–367. 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.