Skip to content
Bench to Insights
Blog

Choosing the Imputation Model

On the previous page mice chose the imputation methods and you accepted them. A default is a choice someone else made, and this page reopens two of those choices and measures what each one moves. The first choice is which variables enter the imputation model at all. The second is how an imputed value is produced once the model exists: drawn from a fitted distribution, or matched to a real observed value.

Which variables belong in the imputation model

Section titled “Which variables belong in the imputation model”

The rule of thumb is generous. Keep the analysis variables, keep anything that predicts the missingness itself, and keep anything correlated with the values being imputed. The two screens below apply that rule mechanically: quickpred computes it in R, and the Python tab prints the same evidence as an explicit correlation table.

suppressPackageStartupMessages(library(mice))


complete_data <- read.csv("/opt/data/cohort_complete.csv")
missing_data <- read.csv("/opt/data/cohort_missing.csv")


# Keep the activity order so the models preserve its rank.
activity_levels <- c("low", "medium", "high")
complete_data$activity <- ordered(
  complete_data$activity,
  levels = activity_levels
)
missing_data$activity <- ordered(
  missing_data$activity,
  levels = activity_levels
)
complete_data$smoking <- factor(complete_data$smoking, levels = c(0, 1))
missing_data$smoking <- factor(missing_data$smoking, levels = c(0, 1))


benchmark_model <- lm(sbp ~ age + bmi + chol + smoking, data = complete_data)
benchmark_term <- summary(benchmark_model)$coefficients["chol", ]
cat(sprintf(
  "Complete data: chol coefficient %.4f (SE %.4f)",
  benchmark_term["Estimate"],
  benchmark_term["Std. Error"]
), "\n", sep = "")


model_data <- missing_data[, setdiff(names(missing_data), "patient_id")]
# An entry is 1 if a pairwise-complete correlation with the target, or with its
# missingness indicator, exceeds the mincor threshold.
predictor_selection <- quickpred(model_data)
cat("quickpred selection matrix (1 = kept predictor):\n")
print(predictor_selection)
Complete data: chol coefficient 3.0370 (SE 0.8299)
quickpred selection matrix (1 = kept predictor):
         age bmi activity smoking chol sbp
age        0   0        0       0    0   0
bmi        0   0        0       0    0   0
activity   1   1        0       0    1   1
smoking    1   1        0       0    1   0
chol       1   1        1       1    0   1
sbp        1   1        1       0    1   0

Read the R matrix by row: the row names are the variables with gaps, and a 1 marks a predictor kept for that variable. Cholesterol keeps everything, which is the expected answer for a variable that correlates with age at 0.464 and with blood pressure at 0.453 in the Python table. The blood pressure row drops smoking, and the correlation table shows why: at 0.088, that pair carries no signal worth imputing from. The two screens are the same evidence in two dialects.

Does the screening move the answer

Section titled “Does the screening move the answer”

Screening is worth measuring rather than assuming. Both runs below use the same seed and the same 20 imputations; the only difference is the predictor matrix.

#' Extract one pooled estimate and its standard error.
#'
#' @param pooled_model A `mice::pool()` result.
#' @param term The coefficient name to extract.
#' @return A numeric vector with the estimate and standard error.
GetPooledTerm <- function(pooled_model, term) {
  pooled_summary <- summary(pooled_model)
  term_row <- pooled_summary[pooled_summary$term == term, ]
  c(estimate = term_row$estimate, std_error = term_row$std.error)
}


#' Print a pooled cholesterol coefficient.
#'
#' @param label The text before the coefficient.
#' @param pooled_model A `mice::pool()` result.
#' @return The function returns no value.
PrintPooledChol <- function(label, pooled_model) {
  pooled_term <- GetPooledTerm(pooled_model, "chol")
  message_text <- sprintf(
    "%s: pooled chol coefficient %.4f (SE %.4f)",
    label,
    pooled_term["estimate"],
    pooled_term["std_error"]
  )
  cat(message_text, "\n", sep = "")
  invisible(NULL)
}


all_predictors <- mice(
  model_data,
  m = 20,
  seed = 20260820,
  printFlag = FALSE
)
selected_predictors <- mice(
  model_data,
  m = 20,
  seed = 20260820,
  predictorMatrix = predictor_selection,
  printFlag = FALSE
)
all_pooled <- pool(with(all_predictors, lm(sbp ~ age + bmi + chol + smoking)))
selected_pooled <- pool(
  with(selected_predictors, lm(sbp ~ age + bmi + chol + smoking))
)
PrintPooledChol("All predictors", all_pooled)
PrintPooledChol("quickpred selection", selected_pooled)
All predictors: pooled chol coefficient 2.8646 (SE 0.9261)
quickpred selection: pooled chol coefficient 2.8596 (SE 0.9568)

Here it barely moves: 2.8646 with the full matrix against 2.8596 with the screened one. That is the honest result on a six-variable cohort where nearly everything correlates with something. Screening pays its rent on wide, noisy tables, where dozens of junk predictors dilute the imputation model; with six variables there is little to dilute. The habit to keep is running the screen and reading it, not expecting it to rescue every data set.

How the imputed value is drawn

Section titled “How the imputed value is drawn”

The second choice is subtler. A parametric method draws from a fitted distribution, and the draw can land anywhere the distribution has mass, including outside any range the observed data ever touched. Predictive mean matching instead takes the patients whose predicted values are closest and returns one of their observed measurements, so an imputed value is always a real patient’s real number. mice’s default for numeric variables is pmm; the norm method is the pure draw. The same fork exists in miceforest: mean_match_candidates set to 5 matches, and set to 0 draws.

#' Print the share outside an observed range.
#'
#' @param label The text before the share.
#' @param values The imputed values to assess.
#' @param observed_range A numeric vector with the observed minimum and maximum.
#' @return The function returns no value.
PrintOutsideRangeShare <- function(label, values, observed_range) {
  outside_share <- mean(
    values < observed_range[1] | values > observed_range[2]
  ) * 100
  cat(sprintf(
    "%s: imputed sbp outside observed range %.2f%%",
    label,
    outside_share
  ), "\n", sep = "")
  invisible(NULL)
}


sbp_missing <- is.na(model_data$sbp)
observed_sbp_range <- range(model_data$sbp, na.rm = TRUE)
pmm_imputation <- mice(
  model_data,
  m = 20,
  seed = 20260823,
  printFlag = FALSE
)
# Keep the default method for every variable except sbp.
default_methods <- mice(
  model_data,
  m = 1,
  maxit = 0,
  seed = 20260823,
  printFlag = FALSE
)$method
normal_methods <- default_methods
normal_methods["sbp"] <- "norm"
normal_imputation <- mice(
  model_data,
  m = 20,
  seed = 20260823,
  method = normal_methods,
  printFlag = FALSE
)


pmm_long <- complete(pmm_imputation, "long")
normal_long <- complete(normal_imputation, "long")
pmm_sbp <- pmm_long$sbp[sbp_missing[as.integer(pmm_long$.id)]]
normal_sbp <- normal_long$sbp[sbp_missing[as.integer(normal_long$.id)]]
PrintOutsideRangeShare("pmm on sbp", pmm_sbp, observed_sbp_range)
PrintOutsideRangeShare("norm on sbp", normal_sbp, observed_sbp_range)


pmm_pooled <- pool(with(pmm_imputation, lm(sbp ~ age + bmi + chol + smoking)))
normal_pooled <- pool(
  with(normal_imputation, lm(sbp ~ age + bmi + chol + smoking))
)
PrintPooledChol("pmm on sbp", pmm_pooled)
PrintPooledChol("norm on sbp", normal_pooled)
pmm on sbp: imputed sbp outside observed range 0.00%
norm on sbp: imputed sbp outside observed range 0.10%
pmm on sbp: pooled chol coefficient 2.9565 (SE 0.9658)
norm on sbp: pooled chol coefficient 2.8114 (SE 1.0135)

The structural difference shows in the measurements. Under norm, 0.10 percent of the imputed blood pressures fall outside the observed range; under pmm, the share is 0.00 percent, and it cannot be anything else, because the method only returns observed values. In the Python tab the same fork is starker: a posterior draw for a binary variable is fractional in 100.00 percent of the imputed cells, patients who smoke 0.3 of a cigarette, while mean matching never fabricates one. The pooled answers move only modestly in return, 2.9565 against 2.8114 in R and 3.7856 against 3.8964 in Python. You do not choose matching because it wins on this cohort. You choose it because of what it is incapable of inventing.

Both engines then agree with the last page in the only way that matters: the model choice moves the estimate less than the decision to impute multiply did, and more than nothing. What remains is checking the work itself, which is the final page.