Handling missing data • gllvmTMB

Missing responses and missing predictors require different models. This guide shows the two current workflows and the two extractors that go with them.

A response cell is missing, such as one unmeasured trait. Fit with miss_control(response = "include"), then inspect the masked cells with predict_missing().

One unit-level predictor is missing. Fit with mi(x), an impute model, and miss_control(predictor = "model"), then inspect it with imputed().

The examples below exercise Gaussian routes that have direct fit-equivalence and recovery tests. Binary, ordered, and unordered missing-predictor models are also implemented, but they have different output contracts and are summarized rather than expanded into more long examples here.

Before fitting: state the assumption

These workflows do not make missingness harmless. Likelihood-based inference requires an ignorable missingness process: after conditioning on variables in the model, the probability that a value is missing must not depend on its unobserved value. The response and predictor models must also be adequate. gllvmTMB does not model a missing-not-at-random process.

By default, miss_control() drops rows with missing responses and rejects missing predictors. The likelihood-based routes in this article are explicit choices:

library(gllvmTMB)

miss_control()
#> $response
#> [1] "drop"
#> 
#> $predictor
#> [1] "fail"
#> 
#> $engine
#> [1] "laplace"
miss_control(response = "include")
#> $response
#> [1] "include"
#> 
#> $predictor
#> [1] "fail"
#> 
#> $engine
#> [1] "laplace"
miss_control(predictor = "model")
#> $response
#> [1] "drop"
#> 
#> $predictor
#> [1] "model"
#> 
#> $engine
#> [1] "laplace"

Measurement error is a different problem. If an observed predictor is noisy, its true value is latent even where the recorded value is present. The mi() workflow below handles absent predictor values, not measurement error.

Missing response cells

Suppose four trait measurements are absent from a two-trait data set. With response = "include", the cells remain identifiable in the stacked data but contribute zero to the response likelihood. The fitted model can then return a plug-in fitted value for each masked cell.

set.seed(42)
n_site <- 24

dat <- data.frame(
  site = factor(seq_len(n_site)),
  z = rnorm(n_site)
)
dat$t1 <- 0.4 + 0.7 * dat$z + rnorm(n_site, sd = 0.4)
dat$t2 <- -0.2 + 0.3 * dat$z + rnorm(n_site, sd = 0.5)
dat$t1[c(3, 12)] <- NA
dat$t2[c(7, 19)] <- NA

colSums(is.na(dat[c("t1", "t2")]))
#> t1 t2 
#>  2  2

Wide and long calls

The wide call uses traits(...). The equivalent long call spells out the trait-specific intercepts and slopes. Both enter through gllvmTMB().

fit_response_wide <- gllvmTMB(
  traits(t1, t2) ~ 1 + z,
  data = dat,
  unit = "site",
  family = gaussian(),
  missing = miss_control(response = "include"),
  control = gllvmTMBcontrol(se = FALSE)
)

dat_long <- rbind(
  data.frame(site = dat$site, trait = "t1", value = dat$t1, z = dat$z),
  data.frame(site = dat$site, trait = "t2", value = dat$t2, z = dat$z)
)
dat_long$trait <- factor(dat_long$trait, levels = c("t1", "t2"))

fit_response_long <- gllvmTMB(
  value ~ 0 + trait + (0 + trait):z,
  data = dat_long,
  trait = "trait",
  unit = "site",
  family = gaussian(),
  missing = miss_control(response = "include"),
  control = gllvmTMBcontrol(se = FALSE)
)

The two forms maximize the same observed-data likelihood in this example.

data.frame(
  format = c("wide", "long"),
  log_likelihood = c(logLik(fit_response_wide), logLik(fit_response_long)),
  likelihood_rows = c(nobs(fit_response_wide), nobs(fit_response_long))
)
#>   format log_likelihood likelihood_rows
#> 1   wide      -26.20456              44
#> 2   long      -26.20456              44

Reconstruct the masked responses

predict_missing() returns one row per missing response cell. Its est column is a fitted linear predictor by default; use type = "response" for the inverse-link mean. These are plug-in fitted values based on the fitted parameters and conditional random-effect estimates. They do not currently carry reconstruction standard errors or prediction intervals.

missing_response <- predict_missing(fit_response_wide, type = "response")
missing_response
#>   original_row model_row site trait        est
#> 1            3         5    3    t1  0.6319083
#> 2            7        14    7    t2  0.2238736
#> 3           12        23   12    t1  2.1353853
#> 4           19        38   19    t2 -0.5909196

For wide input, original_row is the row of dat, while model_row is the internally stacked cell row. Join on site plus trait when identifying a cell; those columns remain stable if the stacking order changes.

nobs(fit_response_wide) is 44 because 48 cells were supplied and four were masked. Retaining a cell’s identity is not the same as adding information: the missing response itself contributes nothing to the likelihood. Estimating cross-trait dependence also requires adequate overlap among observed traits within the relevant units or covariance groups.

One missing continuous predictor

A missing predictor must be declared in both parts of the joint model:

mi(x) places the unit-level predictor in the response formula;
impute = list(x = x ~ z) specifies a model for the predictor;
miss_control(predictor = "model") activates the route.

The following wide-format example models one continuous predictor. The predictor is observed once per site and is shared by both trait responses.

set.seed(47)
n_site <- 30

pred_dat <- data.frame(
  site = factor(seq_len(n_site)),
  z = rnorm(n_site)
)
pred_dat$x <- 0.3 + 0.8 * pred_dat$z + rnorm(n_site, sd = 0.5)
pred_dat$t1 <- 0.5 + 1.1 * pred_dat$x + rnorm(n_site, sd = 0.4)
pred_dat$t2 <- -0.2 + 1.1 * pred_dat$x + rnorm(n_site, sd = 0.4)
pred_dat$x[c(4, 13, 21)] <- NA

fit_predictor <- gllvmTMB(
  traits(t1, t2) ~ 1 + mi(x),
  data = pred_dat,
  unit = "site",
  family = gaussian(),
  impute = list(x = x ~ z),
  missing = miss_control(predictor = "model"),
  control = gllvmTMBcontrol(se = TRUE)
)

imputed_wide <- imputed(fit_predictor)
imputed_wide
#>   variable level level_id original_row model_row observed   estimate std_error
#> 1        x  site        4            4         4    FALSE  0.3325135 0.2144683
#> 2        x  site       13           13        13    FALSE  0.6614312 0.2157783
#> 3        x  site       21           21        21    FALSE -0.5419708 0.2187513
#>             source uncertainty_status
#> 1 conditional_mode                 ok
#> 2 conditional_mode                 ok
#> 3 conditional_mode                 ok

For a continuous Gaussian predictor, estimate is a conditional mode and std_error is a local Gaussian/Laplace conditional standard error when the joint Hessian is available. It is not a multiple-imputation draw, a Bayesian posterior interval, or evidence of calibrated interval coverage. The level and level_id columns are the safe join key. Here they identify site; the row ordinals refer to predictor-model levels and are not general data-row keys.

The same model can start from long data. The shared mi(x) coefficient applies to both traits, while 0 + trait gives each trait its own intercept.

pred_long <- rbind(
  data.frame(site = pred_dat$site, trait = "t1", value = pred_dat$t1,
             x = pred_dat$x, z = pred_dat$z),
  data.frame(site = pred_dat$site, trait = "t2", value = pred_dat$t2,
             x = pred_dat$x, z = pred_dat$z)
)
pred_long$trait <- factor(pred_long$trait, levels = c("t1", "t2"))

fit_predictor_long <- gllvmTMB(
  value ~ 0 + trait + mi(x),
  data = pred_long,
  trait = "trait",
  unit = "site",
  family = gaussian(),
  impute = list(x = x ~ z),
  missing = miss_control(predictor = "model"),
  control = gllvmTMBcontrol(se = TRUE)
)

imputed_long <- imputed(fit_predictor_long)
merge(
  imputed_wide[c("level", "level_id", "estimate")],
  imputed_long[c("level", "level_id", "estimate")],
  by = c("level", "level_id"),
  suffixes = c("_wide", "_long"),
  sort = FALSE
)
#>   level level_id estimate_wide estimate_long
#> 1  site        4     0.3325135     0.3325135
#> 2  site       13     0.6614312     0.6614312
#> 3  site       21    -0.5419708    -0.5419708

Other implemented predictor types

The response formula still uses mi(x). The predictor family is set inside impute_model().

Continuous: impute_model(x ~ z, family = gaussian()), or the shorthand x ~ z, returns a conditional mode and, when available, a conditional SE.
Binary: impute_model(x ~ z, family = binomial()) returns the fitted conditional probability of level 1, without a parameter-uncertainty SE.
Ordered: impute_model(x ~ z, family = cumulative_logit()) returns a fitted conditional expected numeric score, without an SE. The public summary does not return the per-level probability vector.
Unordered: impute_model(x ~ z, family = categorical()) returns the numeric code of the fitted modal category, without an SE. The public summary does not currently return its category label or per-level probabilities.

The engine marginalizes over discrete missing states under the estimated joint model. imputed() exposes a probability for the binary case and a one-number summary for the ordered and unordered cases; it does not expose parameter uncertainty for any of them.

Gaussian predictor models support fixed effects, one grouped random intercept, one coarser level declared with mi_group(), or an intercept-only phylogenetic effect. The discrete predictor models are fixed-effect models. See impute_model() for the current call shapes. Do not extrapolate those routes to spatial, animal, count, bounded, or multiple simultaneous missing predictors.

Choose another route when needed

More than one predictor is missing: use a suitable external multiple-imputation or joint-model method.
Missingness may depend on the unobserved value: perform a sensitivity analysis with a method that models the missingness mechanism.
An observed predictor has measurement error: use a measurement-error model; mi() is not one.
You need uncertainty intervals for reconstructed responses: treat predict_missing() as point output and use a method with predictive-uncertainty support.
You need completed data sets for several downstream analyses: use a multiple-imputation workflow rather than treating imputed() as completed data.

gllvmTMB estimates its supported missing-predictor component jointly with the response model. Joint fitting makes the assumptions explicit, but it does not guarantee that the predictor model is correctly specified or automatically congenial with every scientific analysis. Multiple imputation remains a different and often preferable workflow when the goal is a set of completed data sets, pooling across downstream models, or broader multivariable missingness.

What the current evidence supports

The evidence is strongest for Gaussian response masking and for selected continuous, binary, ordered, unordered, grouped, and phylogenetic predictor models under known generating processes. It does not establish nominal interval coverage across missingness mechanisms, response families, covariance structures, or missingness rates. In particular, this article is not an MNAR robustness study.

Next decision

Proceed to covariance interpretation only after the retained-row counts, missingness assumption, predictor model (when used), and the output from check_gllvmTMB() are all reported. Use Can I trust this fit? for the post-fit checks. If the missingness mechanism or imputation model cannot be defended, stop and compare a sensitivity analysis or an external multiple- imputation workflow rather than treating one completed fit as definitive.

References

Rubin, D. B. (1976). Inference and missing data. Biometrika, 63, 581–592.
Nakagawa, S. & Freckleton, R. P. (2008). Missing inaction: the dangers of ignoring missing data. Trends in Ecology & Evolution, 23, 592–596.
Schafer, J. L. & Graham, J. W. (2002). Missing data: our view of the state of the art. Psychological Methods, 7, 147–177.
Ibrahim, J. G., Chen, M.-H., Lipsitz, S. R. & Herring, A. H. (2005). Missing-data methods for generalized linear models: a comparative review. Journal of the American Statistical Association, 100, 332–346.
Kristensen, K., Nielsen, A., Berg, C. W., Skaug, H. & Bell, B. M. (2016). TMB: automatic differentiation and Laplace approximation. Journal of Statistical Software, 70(5), 1–21.