Multiple imputation and downstream models

Choose a route

• You want to fit a downstream model to several completed datasets. If continuous traits are missing, the posterior route may be eligible only with one row per species, no covariates, and an analysis meeting the conditions below. Its inferential validity has not been established.
• One continuous predictor is missing, while the outcome and other model variables are complete. If missingness is at random (MAR), see the separate experimental multi_impute_analysis() route documentation. Its eligibility rules differ from those for posterior draws.
• You want to explore spread around point predictions. Conformal or MC-dropout draws are for prediction diagnostics, not downstream inference; with_imputations() refuses them for that purpose.

The registered validation did not establish reliable downstream inference for either route in the tested settings. These are different workflows with different requirements.

multi_impute() offers three ways to draw the missing cells. The default, draws_method = "auto", uses "posterior" when the trait types and supplied arguments fit that method (below). If they do not, it falls back to "conformal" with a message explaining why. Inputs that cannot be imputed, such as data with no missing cells or fewer than two observed values for a trait, produce an error instead of a fallback.

The conformal and MC-dropout draws perturb each missing cell around a point prediction. They do not draw the cells jointly from their distribution given the observed data. In the conformal route, continuous, count, ordinal, and proportion traits use the reported pred$se as a scale input for fallback draws when no finite conformal score is available. For log-transformed continuous traits, log1p-transformed counts, and logit-transformed proportions, pred$se is reported on the original scale, but the draw is made on the transformed scale. Pigauto does not fully convert pred$se between these scales, so diagnostic draw spread may be too wide or too narrow. Fitted point predictions are unchanged. In a simulation of a phylogenetic regression slope (16 regimes: two traits, n = 300 or 1000, Pagel’s lambda 1 or 0.5, random or clade-structured missingness, 120 replicates each), conformal draws biased the pooled slope towards zero by 0.20 to 0.46 and the pooled 95% intervals covered the true slope in 0 to 17% of replicates. Other estimands were not measured. with_imputations() therefore refuses these draws. A list of fits built by hand is still pooled by pool_mi(), with a warning that its provenance cannot be checked; do not do this with conformal or MC-dropout draws.

This applies to the draws, not to the per-cell conformal intervals that impute() reports. Those intervals summarise uncertainty about a single imputed value and were close to 95% coverage in simulation (0.954 to 0.971, GNN off, random missingness), but coverage can fall below 95% at small n (about 100 species) and when missingness is clustered in clades.

draws_method = "posterior" fits a Bayesian multivariate phylogenetic mixed model and samples each completed dataset from the posterior predictive distribution of the missing cells. It currently requires:

• continuous traits only,
• one row per species (no species_col),
• no covariates,
• at least one missing cell and at least two observed values for each trait.

For the documented downstream route, the model must be linear in the imputed traits on the scale used for imputation, and it must use only those traits. For example, do not add an external covariate, a squared trait, or an interaction. with_imputations() runs your fitting function on each completed dataset and returns the fits. pool_mi() combines fixed-effect estimates and their uncertainty across successful fits. Neither function checks whether the model formula meets these conditions. Passing the software checks does not establish that the analysis meets them or that its inference is validated.

Evidence and limits

These results come from separate simulation campaigns and answer related, distinct questions. The earlier campaign estimates coverage under its listed conditions. The later registered campaign applies prespecified support criteria across a broader, frozen set of scenarios. Read each result within its own scope; for the current inferential evidence, see the registered results.

Later rerun of the earlier simulation campaign. With pigauto’s separation prior, continuous traits, 30% missing completely at random, 200 datasets per condition, n = 100 to 1000, Pagel’s lambda 0.3 to 1, and trait correlation 0 or 0.5, Rubin-pooled downstream PGLS slopes and phylogenetic correlations had 95% interval coverage from 0.925 to 0.975 across the tested conditions. At lambda = 1 and trait correlation = 0.5, pigauto slope coverage was 0.955, 0.940, and 0.965 at n = 100, 300, and 1000. At n = 100, pigauto’s slope bias relative to the complete-data estimate was about -0.037; the matched frequentist multiple-imputation comparator had bias about -0.021. These results describe this simulation, not general calibration.

Later registered validation. Across 44 tested inferential rows, none was supported: 12 failed and 32 were indeterminate. This includes complete-data controls and is limited to the frozen simulation settings. Within those results, five of eight posterior slope regimes failed and three were indeterminate, with empirical 95% coverage from 86.2% to 95.2%. See the registered inferential results for scope and criteria.

Sampler details. The current separation prior uses half-Cauchy priors for residual standard deviations, scaled to the observed latent traits, and an LKJ(1) prior for the full residual-correlation matrix. Before applying the numerical cutoff, the LKJ(1) prior does not give each pairwise correlation a uniform distribution when there are more than two traits. For numerical stability, the sampler gives zero prior density if the smallest eigenvalue of the residual covariance matrix, standardized by observed latent-trait standard deviations, is below 0.000001. In plain language, on this standardized scale, variance along every unit-length direction must be at least 0.000001. This cutoff narrows the prior support. Before pigauto 0.11.0.9002, the sampler used an inverse-Wishart residual prior. In the earlier campaign, that legacy prior attenuated the relationship between correlated traits at lambda near 1 (coverage 0.75 at n = 1000). It remains selectable only when reproducing that analysis: set posterior_control = list(residual_prior = "iw"). New posterior fits use the separation prior described above by default.

As a separate, experimental option, multi_impute_analysis() declares the substantive analysis model and supports one incomplete double-precision continuous covariate under missing at random (MAR); the outcome and all other analysis variables must be fully observed. Current supported models are a Gaussian lm, binomial-logit glm, or Gaussian lmer with one random intercept. Validity also depends on a correctly specified conditional imputation model. For the binomial-logit glm (SMCFCS) and Gaussian lmer (jomo) routes, inspect the engine’s convergence diagnostics yourself: pigauto captures engine warnings and stops if either route emits them, but it does not calculate or certify sampler convergence. The Gaussian lm route uses a separate Bayesian Normal-regression sampler. This route does not use pigauto’s phylogenetic baseline or GNN and does not propagate posterior-tree uncertainty. See the documented conditions for multi_impute_analysis(). This experimental route does not validate the posterior-draw route described above.

Phylogenetic regression (nlme)

By default, multi_impute() logs each continuous trait whose observed values are all positive, then returns completed data on the original scale. In this example, the log(...) formula uses the scale on which those traits were imputed. If your input columns already contain logged values, set log_transform = FALSE and use those columns directly in the model, for example Wing.Length ~ Mass; wrapping them in log() would log them again.

fits <- with_imputations(mi, function(d) {
  d$species <- rownames(d)
  nlme::gls(log(Wing.Length) ~ log(Mass), data = d,
            correlation = ape::corBrownian(phy = tree300, form = ~species),
            method = "ML")
})
pool_mi(fits)

Compatibility example: mixed models (lme4 and glmmTMB)

These examples add a diet effect after imputation. Diet is not part of the posterior imputation model, so these fits are outside its documented downstream-inference conditions.

Variables that were not imputed can be joined to each completed dataset inside the fitting function. Here diet category (Trophic.Level) becomes a random intercept.

add_diet <- function(d) {
  d$Trophic.Level <- avonet300$Trophic.Level[
    match(rownames(d), avonet300$Species_Key)]
  d
}

fits_lmer <- with_imputations(mi, function(d) {
  lme4::lmer(log(Wing.Length) ~ log(Mass) + (1 | Trophic.Level),
             data = add_diet(d))
})
pool_mi(fits_lmer)

fits_tmb <- with_imputations(mi, function(d) {
  glmmTMB::glmmTMB(log(Wing.Length) ~ log(Mass) + (1 | Trophic.Level),
                   data = add_diet(d))
})
pool_mi(fits_tmb)

For glmmTMB fits, pool_mi() pools the conditional (mean) model only; zero-inflation and dispersion coefficients are not included.

A caveat on joined variables: the posterior imputation model conditions on the imputed traits and the tree, not on Trophic.Level. When a variable that matters for the analysis is left out of the imputation model, its association with the imputed traits is weakened in the imputed cells. Because the posterior route accepts continuous traits only, a categorical variable such as diet cannot currently be added to the imputation model; treat its estimated effects with this in mind. (Diet is also a small grouping factor here, four levels, one of them with a single species, so the random-intercept variance is poorly estimated.)

Compatibility example: location-scale model (drmTMB)

This model estimates how residual spread changes with body mass, while the imputation model assumes constant residual variance. This example demonstrates software compatibility and falls outside the documented posterior-inference conditions.

drmTMB models the mean and the residual standard deviation together. pool_mi() pools both sets of fixed effects, labelled mu: and sigma:. The posterior imputation model has a constant residual variance, so the imputed cells carry no relationship between body mass and residual spread. We expect this to weaken a pooled sigma: slope when many cells are imputed; this has not been measured.

fits_drm <- with_imputations(mi, function(d) {
  drmTMB::drmTMB(drmTMB::bf(log(Wing.Length) ~ log(Mass),
                            sigma ~ log(Mass)),
                 data = d)
})
pool_mi(fits_drm)

Compatibility example: latent-variable model (gllvmTMB)

This joint latent-variable model demonstrates software compatibility and falls outside the documented posterior-inference conditions. pool_mi() below combines its fixed trait means; it does not pool latent loadings or variance components.

gllvmTMB fits all traits jointly in long format. pool_mi() pools its fixed effects (here the trait means); latent loadings and variance components are not pooled.

to_long <- function(d) {
  tr <- names(d)
  data.frame(species = factor(rep(rownames(d), times = length(tr))),
             trait   = factor(rep(tr, each = nrow(d)), levels = tr),
             value   = log(unlist(d, use.names = FALSE)))
}

fits_gllvm <- with_imputations(mi, function(d) {
  gllvmTMB::gllvmTMB(value ~ 0 + trait +
                       latent(0 + trait | species, d = 1, unique = FALSE),
                     data = to_long(d), trait = "trait", unit = "species",
                     family = gaussian(), REML = FALSE,
                     engine = "tmb", silent = TRUE)
})
pool_mi(fits_gllvm)

Other models

For other model classes, pass extractor functions to pool_mi() through coef_fun and vcov_fun, or a tidy_fun returning term, estimate and std.error. The same named, finite fixed-effect terms must be available across the retained fits; inspect reported fit failures before using pooled estimates.