Cross-Sectional Data Imputation with MICE

Libraries

library(tidyverse)
library(mice)
library(naniar)
library(patchwork)

theme_set(
  theme_minimal(base_size = 12) +
    theme(plot.title       = element_text(face = "bold"),
          strip.text       = element_text(face = "bold"),
          panel.grid.minor = element_blank())
)

Dataset

The Salaries data contain the 2008–09 nine-month academic salaries of 397 faculty members at a college in the United States. The same data are available in R as carData::Salaries (Fox & Weisberg, 2019).

Variable	Type	Description
rank	categorical (3 levels)	Academic rank: AsstProf, AssocProf, Prof
discipline	categorical (2 levels)	Type of department: A (theoretical) or B (applied)
yrs.since.phd	numeric	Years since the PhD
yrs.service	numeric	Years of service
sex	categorical (2 levels)	Female, Male
salary	numeric	Nine-month salary (US dollars)

path <- "https://gist.githubusercontent.com/novrisuhermi/b0157aecf4e826b16698a3393fce73fe/raw/55fd29ee76f61016e75eee9b47e8d50455d2a350/salaries.csv"
df <- read.csv(path) %>%
  mutate(rank       = factor(rank, levels = c("AsstProf", "AssocProf", "Prof")),
         discipline = factor(discipline),
         sex        = factor(sex))

num_vars <- names(select(df, where(is.numeric)))   # yrs.since.phd, yrs.service, salary
cat_vars <- names(select(df, where(is.factor)))    # rank, discipline, sex

df

summary(df)

        rank     discipline yrs.since.phd    yrs.service        sex          salary      
 AsstProf : 67   A:181      Min.   : 1.00   Min.   : 0.00   Female: 39   Min.   : 57800  
 AssocProf: 64   B:216      1st Qu.:12.00   1st Qu.: 7.00   Male  :358   1st Qu.: 91000  
 Prof     :266              Median :21.00   Median :16.00                Median :107300  
                            Mean   :22.31   Mean   :17.61                Mean   :113706  
                            3rd Qu.:32.00   3rd Qu.:27.00                3rd Qu.:134185  
                            Max.   :56.00   Max.   :60.00                Max.   :231545  

The data have no missing values, so they serve as the ground truth for the experiment.

Relationships Between Variables

Imputation models can only exploit relationships that exist in the data.

Distributions

p_num <- df %>%
  pivot_longer(all_of(num_vars), names_to = "variable") %>%
  mutate(variable = factor(variable, levels = num_vars)) %>%
  ggplot(aes(value)) +
  geom_histogram(bins = 25, fill = "steelblue4", colour = "white") +
  facet_wrap(~ variable, scales = "free") +
  scale_x_continuous(labels = scales::label_number(scale_cut = scales::cut_short_scale())) +
  labs(title = "Numeric variables", x = NULL, y = "Count")

p_cat <- df %>%
  pivot_longer(all_of(cat_vars), names_to = "variable", values_to = "category") %>%
  mutate(variable = factor(variable, levels = cat_vars)) %>%
  ggplot(aes(category)) +
  geom_bar(fill = "steelblue4") +
  facet_wrap(~ variable, scales = "free_x") +
  labs(title = "Categorical variables", x = NULL, y = "Count")

p_num / p_cat

Salary, Rank and Seniority

p_scatter <- ggplot(df, aes(yrs.since.phd, salary, colour = rank)) +
  geom_point(alpha = 0.7) +
  scale_y_continuous(labels = scales::dollar) +
  scale_colour_brewer(palette = "Dark2") +
  labs(title = "Salary by seniority and rank",
       x = "Years since PhD", y = "Salary", colour = "Rank") +
  theme(legend.position = "bottom")

p_box <- df %>%
  pivot_longer(all_of(cat_vars), names_to = "variable", values_to = "category") %>%
  mutate(variable = factor(variable, levels = cat_vars)) %>%
  ggplot(aes(category, salary)) +
  geom_boxplot(fill = "steelblue4", alpha = 0.3, outlier.size = 0.8) +
  facet_wrap(~ variable, scales = "free_x") +
  scale_y_continuous(labels = scales::dollar) +
  labs(title = "Salary by category", x = NULL, y = "Salary")

p_scatter | p_box

Association Tests

# H0: the two variables are not associated
#   numeric - numeric         : Spearman correlation test  (effect size: |rho|)
#   numeric - categorical     : Kruskal-Wallis test        (effect size: epsilon^2 = H / (n - 1))
#   categorical - categorical : chi-squared test           (effect size: Cramer's V)
association_test <- function(x, y) {
  n <- length(x)
  if (is.numeric(x) && is.numeric(y)) {
    res <- cor.test(x, y, method = "spearman", exact = FALSE)
    tibble(test = "Spearman", effect = abs(unname(res$estimate)), p_value = res$p.value)
  } else if (is.factor(x) && is.factor(y)) {
    res <- chisq.test(table(x, y), correct = FALSE)
    k   <- min(nlevels(x), nlevels(y))
    tibble(test = "Chi-squared", effect = sqrt(unname(res$statistic) / (n * (k - 1))),
           p_value = res$p.value)
  } else {
    res <- if (is.factor(x)) kruskal.test(y, x) else kruskal.test(x, y)
    tibble(test = "Kruskal-Wallis", effect = unname(res$statistic) / (n - 1),
           p_value = res$p.value)
  }
}

associations <- combn(names(df), 2, simplify = FALSE) %>%
  map_dfr(~ bind_cols(tibble(variable_1 = .x[1], variable_2 = .x[2]),
                      association_test(df[[.x[1]]], df[[.x[2]]]))) %>%
  arrange(p_value)

associations %>%
  mutate(effect = round(effect, 3))

Predictability of Each Variable

MICE imputes each variable from all the others, so it helps to know how well each variable can be predicted from the rest.

# numeric     : R-squared of a linear regression on all other variables
# categorical : accuracy of the most likely category from a (multinomial) logistic regression,
#               next to the baseline of always predicting the most frequent category
predictability <- map_dfr(names(df), function(v) {
  f <- reformulate(setdiff(names(df), v), response = v)
  if (is.numeric(df[[v]])) {
    tibble(variable = v, measure = "R-squared",
           model = summary(lm(f, data = df))$r.squared, baseline = NA_real_)
  } else {
    fit <- nnet::multinom(f, data = df, trace = FALSE)
    tibble(variable = v, measure = "Accuracy",
           model    = mean(predict(fit, df) == df[[v]]),
           baseline = max(prop.table(table(df[[v]]))))
  }
})

predictability %>%
  mutate(across(c(model, baseline), ~ round(.x, 3)))

Simulating Missing Values

inject_missing <- function(data, prop = 0.10, mech = "MAR", type = "RIGHT") {
  vars <- names(data)
  p    <- length(vars)
  if (prop * p >= 1) stop("prop times the number of variables must be below 1.")

  patterns <- matrix(1, p, p, dimnames = list(vars, vars))
  diag(patterns) <- 0

  weights <- matrix(0, p, p, dimnames = list(vars, vars))
  weights[, "yrs.since.phd"] <- 1
  weights["yrs.since.phd", c("yrs.since.phd", "yrs.service")] <- c(0, 1)

  data_num <- mutate(data, across(where(is.factor), as.integer))
  amp <- ampute(data_num, prop = prop * p, patterns = patterns, freq = rep(1 / p, p),
                mech = mech, weights = if (mech == "MAR") weights else NULL, type = type)

  data[is.na(amp$amp)] <- NA
  data
}

set.seed(2026)
df_missing <- inject_missing(df, prop = 0.10)

df_missing

{
print(miss_var_summary(df_missing))       # missing cells per variable
print(mean(complete.cases(df_missing)))   # share of complete rows
}

# A tibble: 6 × 3
  variable      n_miss pct_miss
  <chr>          <int>    <num>
1 discipline        40    10.1 
2 yrs.since.phd     40    10.1 
3 rank              39     9.82
4 sex               39     9.82
5 yrs.service       34     8.56
6 salary            34     8.56
[1] 0.4307305

Although each variable misses only about 10% of its values, fewer than half of the rows are complete. Complete-case analysis (the default of lm()) would discard all the other rows, including the values they do contain.

# Left: location of the missing cells
p_vis <- vis_miss(df_missing) +
  labs(title = "Missing cells") +
  theme(plot.title = element_text(face = "bold"))

# Right: seniority of the rows where each variable is missing versus observed
p_mar <- df_missing %>%
  bind_shadow() %>%
  select(yrs.since.phd, ends_with("_NA"), -yrs.since.phd_NA) %>%
  drop_na(yrs.since.phd) %>%
  pivot_longer(-yrs.since.phd, names_to = "variable", values_to = "status") %>%
  mutate(variable = factor(str_remove(variable, "_NA$"), levels = names(df)),
         status   = if_else(status == "NA", "Missing", "Observed")) %>%
  ggplot(aes(status, yrs.since.phd, fill = status)) +
  geom_boxplot(alpha = 0.7, outlier.size = 0.8) +
  facet_wrap(~ variable, nrow = 1) +
  scale_fill_manual(values = c(Missing = "firebrick", Observed = "steelblue4"), guide = "none") +
  labs(title = "Years since PhD by missingness status", x = NULL, y = "Years since PhD") +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))

p_vis | p_mar

# Missing-data patterns: blue = observed, red = missing;
# left = number of rows, right = number of missing variables, bottom = missing cells per variable
{md.pattern(df_missing, rotate.names = TRUE)}

    yrs.service salary rank sex discipline yrs.since.phd    
171           1      1    1   1          1             1   0
40            1      1    1   1          1             0   1
40            1      1    1   1          0             1   1
39            1      1    1   0          1             1   1
39            1      1    0   1          1             1   1
34            1      0    1   1          1             1   1
34            0      1    1   1          1             1   1
             34     34   39  39         40            40 226

# Little's test. H0: the data are missing completely at random (MCAR)
mcar_test(df_missing)

Imputation

Methods

Family	Method	Numeric variables	Categorical variables	m
Mean-based	mean	mean (mean)	mode (mode*)	1
	median	median (median*)	mode (mode*)	1
Random	random	random draw from the observed values (sample)	same	20
Regression	regression	predicted value (norm.predict)	logistic regression, most likely category (polyreg.predict*)	1
	stochastic_regression	prediction + random residual (norm.nob)	logistic regression, draw from the predicted probabilities (logreg / polyreg)	20
	bayesian_regression	Bayesian linear regression (norm)	logistic regression, draw (logreg / polyreg)	20
Predictive Mean Matching	pmm	predictive mean matching (pmm)	logistic regression, draw (logreg / polyreg)	20
	pmm_polr	predictive mean matching (pmm)	as pmm, but ordinal logistic regression for rank (polr)	20
Tree-based	cart	regression tree (cart)	classification tree (cart)	20
	rf	random forest (rf)	random forest (rf)	20

# Arguments passed by mice:
#   y : the variable to impute        ry : TRUE where y is observed
#   x : matrix of predictors          wy : TRUE where y must be imputed
mice.impute.median <- function(y, ry, x, wy = !ry, ...) {
  rep(median(y[ry]), sum(wy))
}

mice.impute.mode <- function(y, ry, x, wy = !ry, ...) {
  counts <- table(y[ry])
  rep(names(counts)[which.max(counts)], sum(wy))
}

# Most likely category from a multinomial logit (a logistic regression for two levels)
mice.impute.polyreg.predict <- function(y, ry, x, wy = !ry, ...) {
  xy  <- data.frame(y = y, x)
  fit <- nnet::multinom(y ~ ., data = xy[ry, , drop = FALSE], trace = FALSE)
  as.character(predict(fit, newdata = xy[wy, , drop = FALSE], type = "class"))
}

Running the Imputations

A dry run with maxit = 0 imputes nothing but returns the settings mice would use:

ini <- mice(df_missing, maxit = 0)

ini$method            # default method per variable

         rank    discipline yrs.since.phd   yrs.service           sex        salary 
    "polyreg"      "logreg"         "pmm"         "pmm"      "logreg"         "pmm" 

ini$predictorMatrix   # rows = imputed variables, columns = predictors (1 = used)

              rank discipline yrs.since.phd yrs.service sex salary
rank             0          1             1           1   1      1
discipline       1          0             1           1   1      1
yrs.since.phd    1          1             0           1   1      1
yrs.service      1          1             1           0   1      1
sex              1          1             1           1   0      1
salary           1          1             1           1   1      0

# One mice run per method: m = 20 imputations and 10 iterations unless stated otherwise
impute <- function(method, m = 20, maxit = 10) {
  mice(df_missing, method = method, m = m, maxit = maxit, printFlag = FALSE)
}

# One method for the numeric variables and one for the factors
# (categorical = NULL keeps the defaults logreg / polyreg)
use_methods <- function(numeric, categorical = NULL) {
  meth <- ini$method
  meth[num_vars] <- numeric
  if (!is.null(categorical)) meth[cat_vars] <- categorical
  meth
}

set.seed(2026)   # mice draws random values
imputations <- list(
  # Mean-based (single imputation): numeric -> mean or median, categorical -> mode
  mean   = impute(use_methods("mean", "mode"),   m = 1, maxit = 1),
  median = impute(use_methods("median", "mode"), m = 1, maxit = 1),

  # Random draws from the observed values of each variable
  random = impute("sample", maxit = 1),

  # Regression: deterministic prediction, prediction + noise, Bayesian draws
  regression            = impute(use_methods("norm.predict", "polyreg.predict"), m = 1),
  stochastic_regression = impute(use_methods("norm.nob")),
  bayesian_regression   = impute(use_methods("norm")),

  # Predictive mean matching (the mice default)
  pmm = impute(use_methods("pmm")),

  # Tree-based
  cart = impute("cart"),
  rf   = impute("rf"),

  # PMM with rank as an ordinal variable: proportional-odds logistic regression
  # (added last, so that the random draws of the runs above stay the same)
  pmm_polr = impute(replace(use_methods("pmm"), "rank", "polr"))
)

imputations$pmm

Class: mids
Number of multiple imputations:  20 
Imputation methods:
         rank    discipline yrs.since.phd   yrs.service           sex        salary 
    "polyreg"      "logreg"         "pmm"         "pmm"      "logreg"         "pmm" 
PredictorMatrix:
              rank discipline yrs.since.phd yrs.service sex salary
rank             0          1             1           1   1      1
discipline       1          0             1           1   1      1
yrs.since.phd    1          1             0           1   1      1
yrs.service      1          1             1           0   1      1
sex              1          1             1           1   0      1
salary           1          1             1           1   1      0

The imputations are stored in $imp: one row per missing cell and one column per completed dataset. The spread within a row shows the uncertainty about that cell.

imputations$pmm$imp$salary[1:6, 1:8]

Convergence

Every imputation is a separate chain. Healthy chains are freely intermingled and show no trend over the iterations. mice stores the mean of the imputed values per iteration and chain in chainMean (for factors, the mean of the category codes); plot(imp) draws the same trace plots with lattice.

iterative <- c("stochastic_regression", "bayesian_regression", "pmm", "pmm_polr", "cart", "rf")

chain_means <- imputations[iterative] %>%
  map_dfr(~ as.data.frame.table(.x$chainMean, responseName = "mean"), .id = "method") %>%
  rename(variable = Var1, iteration = Var2, chain = Var3) %>%
  mutate(method    = factor(method, levels = iterative),
         variable  = factor(variable, levels = names(df)),
         iteration = as.integer(as.character(iteration)))

ggplot(chain_means, aes(iteration, mean, group = chain, colour = chain)) +
  geom_line(alpha = 0.6, linewidth = 0.4) +
  facet_grid(variable ~ method, scales = "free_y") +
  scale_colour_viridis_d(guide = "none") +
  scale_x_continuous(breaks = c(1, 5, 10)) +
  scale_y_continuous(labels = scales::label_number(scale_cut = scales::cut_short_scale())) +
  labs(title = "Trace plots: mean of the imputed values per iteration",
       subtitle = "One line per chain (imputation)", x = "Iteration", y = NULL)

The chains are well mixed and show no clear trend, so 10 iterations are sufficient.

Plausibility

Without the ground truth, the next check is whether the imputed values are plausible. mice::complete(imp, "long", include = TRUE) stacks the incomplete data (.imp = 0) and all completed datasets (.imp = 1, ..., m); densityplot(imp) draws the lattice version of the plot below.

# TRUE where a cell was removed, in long format
miss_long <- as_tibble(is.na(df_missing)) %>%
  mutate(.id = row_number()) %>%
  pivot_longer(-.id, names_to = "variable", values_to = "was_missing")

density_data <- mice::complete(imputations$pmm, action = "long", include = TRUE) %>%
  mutate(.id = as.integer(.id)) %>%
  select(.imp, .id, all_of(num_vars)) %>%
  pivot_longer(all_of(num_vars), names_to = "variable") %>%
  left_join(miss_long, by = c(".id", "variable")) %>%
  mutate(variable = factor(variable, levels = num_vars))

ggplot() +
  geom_density(data = filter(density_data, .imp > 0, was_missing),
               aes(value, group = .imp, colour = "Imputed"), linewidth = 0.3, key_glyph = "path") +
  geom_density(data = filter(density_data, .imp == 0, !was_missing),
               aes(value, colour = "Observed"), linewidth = 1.1, key_glyph = "path") +
  facet_wrap(~ variable, scales = "free") +
  scale_colour_manual(values = c(Observed = "steelblue4", Imputed = "firebrick"),
                      breaks = c("Observed", "Imputed")) +
  scale_x_continuous(labels = scales::label_number(scale_cut = scales::cut_short_scale())) +
  labs(title = "PMM: observed values versus imputed values",
       subtitle = "One red line per completed dataset", x = NULL, y = "Density", colour = NULL) +
  theme(legend.position = "bottom", axis.text.y = element_blank())

Evaluation

The removed values are known, so every method can be compared with the truth. For the set \mathcal{M} of the n_{\mathcal{M}} removed cells of a numeric variable:

\text{RMSE} = \sqrt{\frac{1}{n_{\mathcal{M}}}\sum_{i \in \mathcal{M}} (\hat{y}_i - y_i)^2}, \qquad \text{MAE} = \frac{1}{n_{\mathcal{M}}}\sum_{i \in \mathcal{M}} \lvert \hat{y}_i - y_i \rvert, \qquad \text{Bias} = \frac{1}{n_{\mathcal{M}}}\sum_{i \in \mathcal{M}} (\hat{y}_i - y_i)

MAPE is not used because yrs.service contains zeros. The SD ratio (SD of the imputed values divided by the SD of the true values) shows whether the spread is preserved. Categorical variables are scored with the accuracy (share of correctly imputed categories) and the macro-F1, the mean of \text{F1} = 2\,\text{TP} / (2\,\text{TP} + \text{FP} + \text{FN}) over the categories. For methods with m = 20, every metric is computed for each completed dataset and then averaged.

# True and imputed values of the removed cells: one row per method, completed dataset and cell
removed_cells <- function(vars, convert = identity) {
  map_dfr(imputations, .id = "method", function(imp) {
    map_dfr(vars, function(v) {
      rows <- as.integer(rownames(imp$imp[[v]]))
      imp$imp[[v]] %>%
        mutate(across(everything(), convert),
               variable = v, row = rows, actual = convert(df[[v]][rows])) %>%
        pivot_longer(-c(variable, row, actual), names_to = ".imp", values_to = "imputed",
                     names_transform = as.integer)
    })
  })
}

df_imputed_num <- removed_cells(num_vars)
df_imputed_cat <- removed_cells(cat_vars, convert = as.character)

df_imputed_num %>%
  mutate(across(c(actual, imputed), ~ round(.x, 2)))

Numerical Variables: RMSE, MAE and Bias

method_groups <- list(
  "Mean-based"               = c("mean", "median"),
  "Random"                   = "random",
  "Regression"               = c("regression", "stochastic_regression", "bayesian_regression"),
  "Predictive Mean Matching" = c("pmm", "pmm_polr"),
  "Tree-based"               = c("cart", "rf")
)

families <- enframe(method_groups, name = "family", value = "method") %>%
  unnest(method)

evaluation_num <- df_imputed_num %>%
  group_by(method, variable, .imp) %>%
  summarise(RMSE       = sqrt(mean((imputed - actual)^2)),
            MAE        = mean(abs(imputed - actual)),
            Bias       = mean(imputed - actual),
            `SD ratio` = sd(imputed) / sd(actual),
            .groups = "drop") %>%
  group_by(method, variable) %>%
  summarise(across(c(RMSE, MAE, Bias, `SD ratio`), mean), .groups = "drop") %>%   # average over the m datasets
  left_join(families, by = "method") %>%
  mutate(variable = factor(variable, levels = num_vars)) %>%
  arrange(variable, RMSE)

# Display: ranked by RMSE within each variable
evaluation_num %>%
  group_by(variable) %>%
  mutate(Rank = row_number()) %>%
  ungroup() %>%
  mutate(across(c(RMSE, MAE, Bias, `SD ratio`), ~ round(.x, 2))) %>%
  select(Variable = variable, Rank, Family = family, Method = method, RMSE, MAE, Bias, `SD ratio`)

Categorical Variables: Accuracy and Macro-F1

macro_f1 <- function(actual, imputed) {
  classes <- union(actual, imputed)
  mean(map_dbl(classes, function(cl) {
    tp <- sum(imputed == cl & actual == cl)
    fp <- sum(imputed == cl & actual != cl)
    fn <- sum(imputed != cl & actual == cl)
    if (tp == 0) 0 else 2 * tp / (2 * tp + fp + fn)
  }))
}

# Method that mice actually used for each categorical variable in each run
cat_methods <- imputations %>%
  map_dfr(~ tibble(variable = cat_vars, model = unname(.x$method[cat_vars])), .id = "method")

evaluation_cat <- df_imputed_cat %>%
  group_by(method, variable, .imp) %>%
  summarise(Accuracy   = mean(imputed == actual),
            `Macro-F1` = macro_f1(actual, imputed),
            .groups = "drop") %>%
  group_by(method, variable) %>%
  summarise(across(c(Accuracy, `Macro-F1`), mean), .groups = "drop") %>%
  left_join(families, by = "method") %>%
  left_join(cat_methods, by = c("method", "variable")) %>%
  mutate(variable = factor(variable, levels = cat_vars)) %>%
  arrange(variable, desc(Accuracy))

# Display: ranked by accuracy within each variable
evaluation_cat %>%
  group_by(variable) %>%
  mutate(Rank = row_number()) %>%
  ungroup() %>%
  mutate(across(c(Accuracy, `Macro-F1`), ~ round(.x, 3))) %>%
  select(Variable = variable, Rank, Family = family, Method = method,
         `Imputed with` = model, Accuracy, `Macro-F1`)

# Proportional-odds check: slope of yrs.since.phd in a separate logit for each cut-off of rank
# (the first split is almost perfectly separated, hence suppressWarnings)
slope <- function(event) {
  fit <- suppressWarnings(glm(event ~ yrs.since.phd, family = binomial, data = df))
  coef(fit)[["yrs.since.phd"]]
}

c(`AsstProf | AssocProf, Prof` = slope(df$rank != "AsstProf"),
  `AsstProf, AssocProf | Prof` = slope(df$rank == "Prof"))

AsstProf | AssocProf, Prof AsstProf, AssocProf | Prof 
                 0.9745747                  0.2300914 

Imputed Values Visualization

plot_imputations <- function(group) {
  methods <- method_groups[[group]]

  data <- df_imputed_num %>%
    filter(method %in% methods, .imp == 1) %>%    # first completed dataset
    left_join(select(evaluation_num, method, variable, RMSE), by = c("method", "variable")) %>%
    mutate(method   = factor(method, levels = methods),
           variable = factor(variable, levels = num_vars)) %>%
    arrange(method, variable) %>%
    mutate(rmse  = if_else(RMSE >= 100, scales::comma(RMSE, accuracy = 1), sprintf("%.2f", RMSE)),
           panel = fct_inorder(sprintf("%s\n%s (RMSE %s)", method, variable, rmse)))

  # Same range on both axes, per variable, so that the dashed line is the diagonal
  limits <- data %>%
    group_by(variable) %>%
    summarise(low = min(actual, imputed), high = max(actual, imputed)) %>%
    right_join(distinct(data, panel, variable), by = "variable") %>%
    pivot_longer(c(low, high), values_to = "value")

  ggplot(data, aes(actual, imputed)) +
    geom_blank(data = limits, aes(value, value)) +
    geom_abline(linetype = "dashed", colour = "grey50") +
    geom_point(colour = "firebrick", alpha = 0.7, size = 1.5) +
    facet_wrap(~ panel, scales = "free", ncol = 3) +
    scale_x_continuous(labels = scales::label_number(scale_cut = scales::cut_short_scale())) +
    scale_y_continuous(labels = scales::label_number(scale_cut = scales::cut_short_scale())) +
    labs(title    = group,
         subtitle = "Removed cells, first completed dataset. Dashed line = perfect imputation",
         x = "True value", y = "Imputed value") +
    theme(aspect.ratio = 1)
}

Mean-based

plot_imputations("Mean-based")

Random

plot_imputations("Random")

Regression

plot_imputations("Regression")

Predictive Mean Matching

plot_imputations("Predictive Mean Matching")

Tree-based

plot_imputations("Tree-based")

Distributions and Category Shares

A good imputation method should also reproduce the distribution of the missing values, not only their centre.

family_colours <- c("Mean-based" = "#E69F00", "Random" = "#CC79A7", "Regression" = "#D55E00",
                    "Predictive Mean Matching" = "#0072B2", "Tree-based" = "#009E73",
                    "True values" = "grey70")
method_levels <- unlist(method_groups, use.names = FALSE)   # methods ordered by family

true_num <- df_imputed_num %>%
  distinct(variable, row, actual) %>%
  transmute(method = "true values", family = "True values", variable, value = actual)

df_imputed_num %>%
  left_join(families, by = "method") %>%
  transmute(method, family, variable, value = imputed) %>%    # all completed datasets
  bind_rows(true_num) %>%
  mutate(method   = factor(method, levels = rev(c("true values", method_levels))),
         family   = factor(family, levels = names(family_colours)),
         variable = factor(variable, levels = num_vars)) %>%
  ggplot(aes(value, method, fill = family)) +
  geom_boxplot(outlier.size = 0.6, alpha = 0.85) +
  facet_wrap(~ variable, scales = "free_x") +
  scale_fill_manual(values = family_colours) +
  scale_x_continuous(labels = scales::label_number(scale_cut = scales::cut_short_scale())) +
  labs(title = "Imputed values versus the removed true values",
       x = NULL, y = NULL, fill = NULL) +
  theme(legend.position = "bottom")

category_colours <- c(AsstProf = "#c6dbef", AssocProf = "#6baed6", Prof = "#08519c",
                      A = "#a1d99b", B = "#238b45", Female = "#fdae6b", Male = "#d94801")

true_cat <- df_imputed_cat %>%
  distinct(variable, row, actual) %>%
  transmute(method = "true values", variable, category = actual)

# Reference: categories of the observed (not removed) cells
obs_cat <- map_dfr(cat_vars, ~ tibble(method = "observed values", variable = .x,
                                       category = as.character(na.omit(df_missing[[.x]]))))

# Labels with the method used for the categorical variables, e.g. "pmm (polyreg/logreg)"
cat_labels <- cat_methods %>%
  group_by(method) %>%
  summarise(models = paste(unique(model), collapse = "/")) %>%
  mutate(label = if_else(models == method, method, paste0(method, " (", models, ")")))

plot_levels <- c("true values", "observed values", method_levels)
plot_labels <- c("true values", "observed values",
                 cat_labels$label[match(method_levels, cat_labels$method)])

df_imputed_cat %>%
  transmute(method, variable, category = imputed) %>%
  bind_rows(true_cat, obs_cat) %>%
  count(method, variable, category) %>%
  group_by(method, variable) %>%
  mutate(share = n / sum(n)) %>%
  ungroup() %>%
  mutate(method   = factor(method, levels = rev(plot_levels), labels = rev(plot_labels)),
         variable = factor(variable, levels = cat_vars),
         category = factor(category, levels = names(category_colours))) %>%
  ggplot(aes(share, method, fill = category)) +
  geom_col(width = 0.75, position = position_stack(reverse = TRUE)) +
  facet_wrap(~ variable) +
  scale_x_continuous(labels = scales::percent) +
  scale_fill_manual(values = category_colours) +
  guides(fill = guide_legend(nrow = 1)) +
  labs(title = "Category shares among the imputed cells",
       subtitle = "In brackets: the method used for the categorical variables",
       x = "Share", y = NULL, fill = NULL) +
  theme(legend.position = "bottom")