MICE always felt like magic to me 🪄 So I tried rebuilding it from scratch in R to understand it better. Lots of mistakes and help along the way, but my version ended up close to the mice package. Might be just luck. Still learning, but here are my messy notes! 🎃📊

Motivations

We used mulitple imputation by chained equations (MICE) quite often in imputing missing data in many occasion, but I never really understood the mechanism behind it except that it appears magical! 🧙 I want to see if I can build this from scratch as a mean to understand the foundation of MICE so that at least it will give me the illusion of understanding. 🤣 Join me in my messy journey of trying to code it from the ground up (FYI, I’m not coding OLS or MLE from the ground up…), and let’s see if we can create some “magic” here! This is going to be a bumpy one with lots of coding, and maybe less of mathematical formulae. I’ll try my best to provide a TL;DR for each section so that we can get a high level view as well. Buckle up! And oh, Happy Halloween! 🎃 What a better time to look at magical recipes of statistics! Let’s go!

Objectives:

• What Is MICE
• Simulate Our DGP
• The Procedure To Impute And Create Multiple Datasets
  • Pool them!
• Opportunities For Improvement
• Lessons Learnt

What Is MICE

TL;DR: fill in each gap several times with plausible random guesses, analyse each completed dataset, then combine the results.

When data are missing, the tempting fix is to plug in one best guess, like the mean. The problem is that a single guess pretends we know the missing value exactly, so our standard errors come out too small. MICE takes a different approach: it fills in each missing value with a plausible random draw, and does this several times to create multiple complete datasets. The differences between those datasets capture our uncertainty about what the missing values really were.

The “chained equations” part is how those draws are made. Instead of one giant model for all variables at once, MICE fits a simple regression for each incomplete variable, using all the others as predictors, and cycles through them, each one feeding its freshly imputed values into the next. After a few rounds, the imputations settle down. Finally, we analyse each completed dataset and pool the results with Rubin’s rules. ✨

Simulate Our DGP

Alright let’s simulate our data generating process (the truth) and add in missing data mechanism to create our dataset that we want to impute from (imagine this is the dataset that we have and we want to apply mice on it).

library(tidyverse)
library(broom)

set.seed(1)
n <- 1000
p1 <- rnorm(n)
p2 <- rnorm(n, 0.5*p1)
p3 <- rnorm(n, 0.4*p1 + 0.3*p2)
p4 <- rbinom(n, 1, 0.5)
p5 <- rbinom(n, 1, plogis(-0.5 + 0.3*p1 - 0.2*p2 + 0.1*p3))
outcome <- rbinom(n, 1, plogis(-0.5 + 0.3*p1 - 0.2*p2 + 0.1*p3 + 0.5*p4 + 0.4*p5))

df <- tibble(p1,p2,p3,p4,p5,outcome)

# missing mechanism 
## mcar
p1_miss <- ifelse(runif(n) < 0.1, NA, p1)
p3_miss <- ifelse(runif(n) < 0.2, NA, p3)

## mar 
p2_miss <- ifelse(runif(n) < plogis(-0.5 + 0.3*p1), NA, p2)
p4_miss <- ifelse(runif(n) < plogis(-0.5 + 0.3*p1 - 0.2*p2), NA, p4)

## mnar
p5_miss <- ifelse(runif(n) < plogis(-0.5 + 0.8*p5), NA, p5)

## recombine
df_miss <- tibble(p1 = p1_miss, p2 = p2_miss, p3 = p3_miss, p4 = p4_miss, p5 = p5_miss, outcome = outcome)

The Procedure To Impute And Create Multiple Datasets

Alright, let’s go through the procedure for MICE real quick.

1. Index all missing values and record them in a dataframe.
2. For each chain:

• Fill every missing value with a random starting value.
• For each column with missing data:
  • Bootstrap the rows where that column is observed, using the current filled-in values for the other columns.
  • Fit a model with that column as the response and all other columns, including the outcome, as predictors: linear regression for continuous columns, logistic for binary.
  • Predict the missing values, then add noise. For continuous columns, add a normal draw using the model’s residual SD. For binary columns, draw 0 or 1 from the predicted probability.
  • Write the new values back immediately, so the next column’s model uses them. This is the “chained” part.
• Repeat the column loop for 10 iterations.
• Keep the final imputed values from the 10th iteration.

3. Repeat step 2 for 10 chains.
4. Build 10 completed datasets, one per chain.
5. Fit the analysis model on each dataset.
6. Pool by averaging the coefficients, plus Rubin’s rules if you need standard errors.

Phew. What a procedure! At least we now know the machinery behind MICE!

Index all missing values and record them in a dataframe

## procedure of simple mice
### 1. index all missing values
df_miss_idx <- df_miss |> is.na()

#### create an empty dataframe based on number of missing data
empty_row <- df_miss_idx |> as.vector() |> sum()

empty_df <- which(df_miss_idx, arr.ind = T) |> as_tibble()

data_type <- df_miss |> 
  map_chr(~class(.x)[1]) 

data_type_colnum <- 1:length(data_type)
colname <- colnames(df_miss)

data_type_df <- tibble(colnum = data_type_colnum, dtype=data_type, colname = colname)

empty_df <- empty_df |>
  left_join(data_type_df, by = c("col" = "colnum")) |>
  mutate(value = NA)

Notice that we added datatype to help us identify continuous data vs binary data?

Fill in missing values with random starting values

Let’s create a an empty dataframe then write functions to do these so we can use them on our for loop

### 2. First iteration, first chain
### 2a. randomly impute all values, use conditions per column to assess types of data before random normal, binomial distribution
fill_empty_df <- function(empty_df) { 
  df <- empty_df |>
  mutate(value = case_when(
    dtype == "numeric" ~ rnorm(n()),
    dtype == "integer" ~ rbinom(n(),1,0.5),
    TRUE ~ NA_real_
  )) |>
  mutate(model = case_when(
    dtype == "numeric" ~ "linreg",
    dtype == "integer" ~ "logreg",
    TRUE ~ NA_character_
  ))
  
  return(df)
}

grab_xvec <- function(df, outcome) {
  df |> 
    select(-all_of(outcome)) |>
    colnames()
}

gen_formula <- function(df, outcome) {
  x_vec <- grab_xvec(df, outcome)
  x <- paste(x_vec,collapse = "+")
  formula <- paste0(outcome,"~",x)
  return(formula)
}


formula_df <- empty_df |>
  distinct(colname, .keep_all = T) |>
  select(col, colname) |>
  mutate(formula = map(.x=colname, .f=~gen_formula(df_miss,.x)) |> unlist()) 


missing_fill <- function(empty_df) {
  df_miss_idx_fill <- matrix(nrow = nrow(df_miss), ncol = ncol(df_miss)) 
  
  for (i_mis in 1:empty_row) {
    rownum <- empty_df[[i_mis, "row"]]
    colnum <- empty_df[[i_mis, "col"]]
    data <- empty_df[[i_mis, "value"]]
    df_miss_idx_fill[rownum,colnum] <- data
  }
  
  df_miss_idx_fill <- df_miss_idx_fill |> as_tibble() 
  colnames(df_miss_idx_fill) <- colnames(df_miss)
  df_miss_idx_fill <- df_miss_idx_fill |>
    mutate(p4 = as.integer(p4),
           p5 = as.integer(p5))
  
  df_miss_df_miss_idx_fill <<- map2_dfc(df_miss, df_miss_idx_fill, coalesce) #learn about map2_dfc
}

empty_df contains randomly pulled number from respective datatype. grab_xvec grabs the predictor columns then gen_formula turns that into model formula for fitting. And missing_fill starts filling out the missing data 👍 Sorry if my code is not efficient. I’m sure we can have LLM make it even more efficient 🤣

Next, we iterate!

# grab the colomuns that have missing values 
col_missing <- empty_df |> distinct(col) |> pull()

# create empty df so we can fill in
empty_df_mice <- tibble()
  
# set up chain
for (chain in 1:10) {
  
  # internally create another df so we can use and override data
  impute_df <- fill_empty_df(empty_df)
  
# for each chain iterate up to 10 times (this is where the chain equation starts)
for (iter in 1:10) {
  
# for each iteration look at each column with missing data
for (i in col_missing) {
  
  # refresh the filled dataset with the latest imputations
  missing_fill(impute_df)
  
  # set seed for each chain, iter, and column so that we can reproduce the same for testing purposes
  seed <- paste0(chain,iter,i) |> as.numeric()
  set.seed(seed)
  
  # look for datatype of the column we're evaluating now
  dtype <- impute_df |>
    filter(col == i) |>
    distinct(model)
  
  # a long winded way of sorting what columns to use for modeling and what to readd back
  empty_df_i <- impute_df |>
    filter(col == i)
  
  empty_df_i_not <- impute_df |>
    filter(col != i)
  
  # finding out what rows are missing and remove those for modelling
  missing_vec <- empty_df_i |>
    pull(row)
  
  df_miss_i_full <- df_miss_df_miss_idx_fill |>
    slice(-missing_vec)
  
  # resample (in place of posterior sampling)
  df_miss_i <- slice_sample(df_miss_i_full, replace = T, n = nrow(df_miss_i_full))
  
  
  # find out the formula for model
  formula <- formula_df |>
    filter(col == i) |>
    pull(formula)
  
  # algo for which estimator to use
  if (dtype == "linreg") {
  model_i <- lm(formula, data = df_miss_i)
  } 
  
  if (dtype == "logreg") {
   model_i <- glm(formula, data = df_miss_i, family = "binomial")   
  }

  newdata <- df_miss_df_miss_idx_fill |>
    slice(missing_vec)
  
  if (dtype == "logreg") {
  empty_df_i <- empty_df_i |>
    mutate(value = rbinom(nrow(empty_df_i), 1, predict(model_i, newdata = newdata, type = "response")))
  } else {
  empty_df_i <- empty_df_i |>
    mutate(value = predict(model_i, newdata = newdata) + rnorm(nrow(empty_df_i),0,sigma(model_i)))
  }
  
  impute_df <- empty_df_i_not |>
    rbind(empty_df_i)
  
}
  # only save the last iter data for each chain
  if (iter == 10) {
    empty_df_mice <- empty_df_mice |>
      rbind(impute_df |> mutate(chain = chain)) 
  }
  
}
}

# write a function to create dataframe based on chains
missing_fill_full <- function(x) {
  df_miss_idx_fill <- matrix(nrow = nrow(df_miss), ncol = ncol(df_miss)) 
  
  for (i_mis in 1:empty_row) {
    rownum <- x[[i_mis, "row"]]
    colnum <- x[[i_mis, "col"]]
    data <- x[[i_mis, "value"]]
    df_miss_idx_fill[rownum,colnum] <- data
  }
  
  df_miss_idx_fill <- df_miss_idx_fill |> as_tibble() 
  colnames(df_miss_idx_fill) <- colnames(df_miss)
  df_miss_idx_fill <- df_miss_idx_fill |>
    mutate(p4 = as.integer(p4),
           p5 = as.integer(p5))
  
  # merges both df and collapse those overlap completely
  map2_dfc(df_miss, df_miss_idx_fill, coalesce) 

}

# create 10 dataframes for each chain
for (data in 1:10) {
  df_miss_mice_i <- missing_fill_full(empty_df_mice |> filter(chain == data))
  assign(paste0("df_miss_mice_",data), df_miss_mice_i)
}

Wow, OK that’s a lot of code! Hopefully the comments above will be self-explainatory. Some of my code may be redundant. Now that we have 10 imputed dataset. The next recipe of the magic is fit each dataset to a model then pool them 🏊‍♂️

Fit All!

# put all 10 dataset into a list (probably should've done this before for efficienty, oh well)
imp_list <- map(.x=1:10, ~ missing_fill_full(empty_df_mice |> filter(chain == .x)))

# fit one at a time and save
fits <- map(imp_list, ~ glm(outcome ~ p1 + p2 + p3 + p4 + p5,
                            data = .x, family = binomial))

Let’s take a peek at the first fit

fits[[1]] |> broom::tidy(conf.int = T)

Alright! Looks pretty ok.

Pool them

# assess the point estimate of coefficient
coef <- 1:6
fit_list <- 1:10
temp_coef <- c()
# temp_coef_var <- c()
model_coef <- c()
model_coef_var <- c()
for (i in coef) {
  for (j in fit_list) {
  temp_coef[j] <- fits[[j]]$coefficients[[i]]
  }
  model_coef[i] <- mean(temp_coef)
  model_coef_var[i] <- var(temp_coef)
}

# assess the pooled SE to construct 95% CI of the coefficients
se <- do.call(rbind, map(.x = fits, ~ tidy(.x) |> pull(std.error))) |> as_tibble()
within_var <- colMeans(se^2) 
between_var <- model_coef_var 
total_var <- within_var + between_var + between_var / length(fit_list)
final_se <- sqrt(total_var)

# degrees of freedom (Barnard–Rubin)
lambda <- (1 + 1/length(fit_list)) * between_var / total_var
df_com <- 1000 - 6
df_old <- (length(fit_list) - 1) / lambda^2
df_obs <- (df_com + 1) / (df_com + 3) * df_com * (1 - lambda)
degfree <- df_old * df_obs / (df_old + df_obs)

tval <- qt(0.975, degfree)
mice_coef <- tibble(term = fits[[1]] |> tidy() |> pull(term), estimate = model_coef, std.error = final_se, tval = tval,
       conf.low = estimate - tval * std.error, conf.high = estimate + tval * std.error) |>
  select(-tval) |>
  mutate(method = "our_mice")

To combine our 10 models, we use Rubin’s rules. For the estimate, we just average the 10 coefficients. For the standard error, we add two kinds of uncertainty: the usual uncertainty within each model, and how much the 10 coefficients disagree with each other because of the missing data. Since that disagreement is estimated from only 10 datasets, we use a t value with Barnard–Rubin degrees of freedom instead of z score of 1.96 for 95% CI.

library(mice)

imp  <- mice(df_miss, m = 10, maxit = 10, seed = 1, printFlag = FALSE)
fit  <- with(imp, glm(outcome ~ p1 + p2 + p3 + p4 + p5, family = binomial))
actual_mice_coef <- tidy(pool(fit), conf.int = TRUE) |>
  select(term, estimate, std.error, conf.low, conf.high) |>
  mutate(method = "actual_mice")

For completeness sake, let’s use the professional’s method and see if our method matches up.

Grab All the Coefficients From Different Models

## function to grab stats
grab_stats <- function(df, method) {
  model <- glm(outcome ~ p1 + p2 + p3 + p4 + p5, df, family = "binomial") 
  coef <- model |> tidy(conf.int = T) |> select(-statistic, -p.value) |> mutate(method = method)
  return(coef)
}

## full data
full_data_coef <- grab_stats(df, "full_data")

## complete case
cco_coef <- grab_stats(df_miss, "complete_case")

## complete case analysis
all_coef <- rbind(mice_coef, full_data_coef, cco_coef, actual_mice_coef)

## truth
truth <- tibble(term = c("(Intercept)", "p1", "p2", "p3", "p4", "p5"),
                truth = c(-0.5, 0.3, -0.2, 0.1, 0.5, 0.4))

## visualize all
all_coef |>
  ggplot(aes(x = estimate, y = term, color = method)) +
  geom_vline(xintercept = 0, linetype = "dotted", colour = "grey60") +
  geom_pointrange(aes(xmin = conf.low, xmax = conf.high),
                  position = position_dodge(width = 0.6)) +
  geom_point(data = truth, aes(x = truth, y = term),
             inherit.aes = FALSE, shape = 4, size = 4) +
  labs(x = "Estimate (log-odds)", y = NULL,
       caption = "× = true coefficient") +
  theme_bw()

Our from-scratch MICE (purple) lines up closely with the mice package (red) for every coefficient, and both stay near the full-data fit (teal). Their intervals are slightly wider than the full-data ones, reflecting the extra uncertainty from the missing data. Complete-case analysis (green) drops every row with any missing value, so its estimates wander and its intervals are much wider. The × marks show the true coefficients from our data-generating process.

For p4, both MICE versions land a bit above the full-data estimate. Since they agree, it’s likely a quirk of this dataset rather than a coding problem. p5 is MNAR, yet both versions look fine here. That may just be luck in a single run, technically it shouldn’t really work, right? 🤔

Hooray! 🙌 We did it! Wow, that was a long journey! I started this since 10/9/26 on the procedure, draft codes etc. Met lots of hurdles over the weeks. I’d be honest with you, it was bumpy! My initial go was a wrong approach where I average out the 10 datasets imputed data (which was wrong), then realized that that is not the proper pooling. Thanks to LLM being consistent tutor and back and forth conversaiton to converge to the proper technique. Now, I know how it works under the hood!

Opportunities For Improvement

• next we should learn sampling from posterior instead of resampling
• Try predictive mean matching, mice’s default. It borrows real observed values instead of adding normal noise, which keeps imputations within realistic ranges.
• Make the code faster and tidier. Vectorise missing_fill(), keep everything in lists instead of global variables, and pass the random generator explicitly
• try out sensitivity analysis with mice

Lessons learnt

• learnt the machinery of mice!
• learnt about map2_dfc
• learnt about Barnard–Rubin degrees of freedom and how to pool using Rubin’s rule.
• learnt total variance formula

If you like this article:

• please feel free to send me a comment or visit my other blogs
• please feel free to follow me on BlueSky, twitter, GitHub or Mastodon
• if you would like collaborate please feel free to contact me