Validating multiple-imputation properness: a coverage simulation
In VIM: Visualization and Imputation of Missing Values

Multiple imputation is proper when Rubin's rules give pooled standard errors that reflect the real uncertainty — so that nominal 95% confidence intervals cover the truth about 95% of the time. This vignette validates vimpute()'s MI behaviour with a known-truth simulation and, just as importantly, shows what goes wrong when a variability source is missing. The chunks below were run when the vignette was precomputed (vignettes/precompute.R in the source repository; the code is shown unchanged and runs as is) with

NSIM <- 12   # simulation replications per configuration -- demo scale!
N    <- 100  # observations
M    <- 5    # imputations

Demo scale means large Monte-Carlo error (a single coverage estimate has a standard error of about 6 percentage points at these settings); the pooled standard errors are stable much earlier and carry the qualitative message. NSIM <- 500, M <- 20 (identical code) is the paper-scale setting.

Design

Linear truth, MAR missingness in the response driven by the observed covariates:

library(VIM)
set.seed(11)

beta_x <- 2   # the estimand

sim_data <- function(n) {
  x <- rnorm(n); z <- rnorm(n)
  data.frame(x = x, z = z, y = 1 + beta_x * x - z + rnorm(n))
}

Four configurations of vimpute(method = "robust", m = M):

default (boot + pmm): bootstrap refits plus stochastic donor draws — VIM's default for multiple imputation since 7.3.0 (boot = TRUE whenever m > 1);
pmm, no boot: donor draws with a single model fit (boot = FALSE) — the single-fit variant, kept for comparison (before 7.3.0 the m > 1 default produced identical conditional-mean imputations);
boot + normalerror: bootstrap parameter uncertainty plus residual noise — the textbook-proper combination;
boot only, uncert = "none": parameter uncertainty without residual variability — deliberately improper. vimpute() warns about exactly this configuration:

w <- tryCatch(
  vimpute(sim_data(60) |> transform(y = replace(y, 1:15, NA)),
          method = "robust", m = 2, boot = TRUE, uncert = "none",
          sequential = FALSE, verbose = FALSE, seed = 1),
  warning = function(cond) conditionMessage(cond))
cat(strwrap(w, 70), sep = "\n")
#> m > 1 with boot = TRUE but uncert = 'none' and no PMM: imputations
#> use conditional means only, so the between-imputation variance is
#> underestimated and pooled standard errors will be too small. Add a
#> residual-noise mechanism (uncert = 'normalerror' or 'resid') or PMM
#> for proper multiple imputation.

The simulation

configs <- list(
  "default (boot + pmm)" = list(boot = TRUE,  uncert = "pmm"),
  "pmm, no boot"         = list(boot = FALSE, uncert = "pmm"),
  "boot + normalerror"   = list(boot = TRUE,  uncert = "normalerror"),
  "boot, uncert none"    = list(boot = TRUE,  uncert = "none")
)

sim_once <- function(cfg, seed) {
  dat <- sim_data(N)
  amp <- makeMissing(dat, prop = 0.3, mechanism = "MAR", vars = "y",
                     seed = seed)
  mi <- suppressWarnings(
    vimpute(amp, method = "robust", m = M, boot = cfg$boot,
            uncert = cfg$uncert, sequential = FALSE, verbose = FALSE))
  pooled <- summary(mice::pool(with(mi, lm(y ~ x + z))), conf.int = TRUE)
  row <- pooled[pooled$term == "x", ]
  c(est = row$estimate, se = row$std.error,
    cover = as.numeric(row$`2.5 %` <= beta_x & beta_x <= row$`97.5 %`))
}

results <- lapply(names(configs), function(cf) {
  runs <- t(vapply(seq_len(NSIM), function(i) sim_once(configs[[cf]], i),
                   numeric(3)))
  data.frame(config = cf,
             mean_pooled_se = mean(runs[, "se"]),
             sd_estimates   = sd(runs[, "est"]),
             coverage       = mean(runs[, "cover"]))
})
results <- do.call(rbind, results)

knitr::kable(results, digits = 3, row.names = FALSE,
             caption = sprintf(
  "Pooled inference for the x-coefficient over %d replications (demo scale). A calibrated method has mean pooled SE close to the empirical SD of the estimates and coverage near 0.95.",
  NSIM))

Table: Pooled inference for the x-coefficient over 12 replications (demo scale). A calibrated method has mean pooled SE close to the empirical SD of the estimates and coverage near 0.95.

|config | mean_pooled_se| sd_estimates| coverage| |:--------------------|--------------:|------------:|--------:| |default (boot + pmm) | 0.129| 0.093| 1.00| |pmm, no boot | 0.120| 0.138| 1.00| |boot + normalerror | 0.130| 0.117| 1.00| |boot, uncert none | 0.093| 0.136| 0.75|

Reading the table

The diagnostic comparison is mean pooled SE vs the empirical SD of the estimates: for calibrated inference they agree, and coverage lands near the nominal level. The deliberately improper configuration (bootstrap only, no residual noise) shows the audit-documented failure: its pooled SE is visibly smaller than the other configurations and than its own sampling variability — anti-conservative intervals whose coverage collapses at paper-scale NSIM. The stochastic configurations — the boot + pmm default, pmm without bootstrap, and the textbook boot + normalerror combination — produce pooled SEs of the right size, with the bootstrap default adding the model-uncertainty component on top of the donor draws.

At NSIM <- 500, M <- 20, N <- 300 the same code reproduces the paper's coverage table; only the three constants at the top change.