---
title: "Single imputation: bias and variance"
description: "Why single imputation misleads: mean imputation attenuates the slope, regression imputation inflates correlations, and one dataset hides the real uncertainty."
date: "2026-05-28 11:00"
categories: [R, missing data, ggplot2, ecology tutorial]
image: thumbnail.png
image-alt: "Three scatter panels showing missing responses filled by mean imputation as a flat stripe, by regression imputation exactly on the fitted line, and by stochastic imputation scattered around it."
---
Filling a gap with a single plausible number and carrying on as if the value were real is the most tempting fix for missing data, and the most quietly damaging. This post works through three single-imputation methods, mean imputation, regression imputation, and its stochastic version, and shows two separate failures: some methods bias the estimate you care about, and *all* of them understate its uncertainty. The second failure is the harder one to see, and it is the reason the next post moves to multiple imputation.
The setting is a response `y` that is missing at random given an observed predictor `x`, with the slope of `y` on `x` as the quantity of interest (true value 1.5).
```{r}
#| label: setup
#| message: false
#| warning: false
#| code-fold: true
#| code-summary: "Setup: packages, colours and plot theme"
library(ggplot2); library(dplyr); library(tidyr)
ink <- "#2b2b2b"; grid_col <- "#d9d5cc"; paper <- "#f5f4ee"
imp_col <- c("Mean imputation" = "#c08a3e",
"Regression imputation" = "#a24b4b",
"Stochastic imputation" = "#4a7a8c")
theme_te <- function(base = 12) {
theme_minimal(base_size = base) +
theme(panel.grid.minor = element_blank(),
panel.grid.major = element_line(colour = grid_col, linewidth = 0.3),
plot.background = element_rect(fill = paper, colour = NA),
panel.background = element_rect(fill = paper, colour = NA),
strip.text = element_text(face = "bold", colour = ink),
plot.title = element_text(face = "bold", colour = ink),
axis.title = element_text(colour = ink),
legend.position = "top")
}
```
## The data and the three fills
```{r}
#| label: dgp
set.seed(2947)
n <- 500; c0 <- 1; c1 <- 1.5; sdx <- 1.2; sde <- 1.6
x <- rnorm(n, 0, sdx)
y <- c0 + c1 * x + rnorm(n, 0, sde)
full_slope <- unname(coef(lm(y ~ x))[2]); full_r <- cor(x, y); full_sdY <- sd(y)
## missing at random: higher x makes y more likely to be missing (about 40%)
a1 <- 1.1
a0 <- uniroot(function(a) mean(plogis(a + a1 * x)) - 0.40, c(-5, 5))$root
miss <- runif(n) < plogis(a0 + a1 * x); obs <- !miss
impute <- function(kind) {
yo <- y; yo[miss] <- NA
fit <- lm(yo[obs] ~ x[obs]); sig <- summary(fit)$sigma
yhat <- coef(fit)[1] + coef(fit)[2] * x[miss]
if (kind == "mean") yo[miss] <- mean(y[obs]) # constant fill
if (kind == "detreg") yo[miss] <- yhat # on the fitted line
if (kind == "stoch") yo[miss] <- yhat + rnorm(sum(miss), 0, sig) # line plus residual noise
yo
}
```
About `r round(100 * mean(miss))`% of the responses are blanked. Mean imputation replaces every gap with the observed mean of `y`. Regression imputation replaces each gap with its prediction from a line fitted to the complete cases. Stochastic regression imputation adds a random residual to that prediction, so the filled points scatter around the line instead of lying on it. Figure 1 shows the three fills against the observed points.
```{r}
#| label: fig-fills
#| fig-cap: "The same missing responses filled three ways. Mean imputation places them in a flat horizontal band; regression imputation places them exactly on the fitted line; stochastic imputation scatters them around it. The dashed line is the full-data fit."
#| fig-alt: "Three scatter panels of y against x. Mean imputation shows imputed points as a horizontal stripe; regression imputation shows them as a straight line of points; stochastic imputation shows them as a cloud around the line."
#| fig-width: 7.5
#| fig-height: 4.4
lev <- c("Mean imputation", "Regression imputation", "Stochastic imputation")
mkp <- function(yo, nm) data.frame(x = x, y = yo, nm = nm,
kind = ifelse(miss, "Imputed", "Observed"))
df1 <- rbind(mkp(impute("mean"), lev[1]), mkp(impute("detreg"), lev[2]),
mkp(impute("stoch"), lev[3]))
df1$nm <- factor(df1$nm, levels = lev)
fl <- coef(lm(y ~ x))
ggplot(df1, aes(x, y)) +
geom_point(data = subset(df1, kind == "Observed"), colour = "#8a8578", size = 0.9, alpha = 0.6) +
geom_point(data = subset(df1, kind == "Imputed"), aes(colour = nm), size = 1.1, alpha = 0.85) +
geom_abline(intercept = fl[1], slope = fl[2], colour = ink, linewidth = 0.5, linetype = "dashed") +
facet_wrap(~nm) + scale_colour_manual(values = imp_col, guide = "none") +
labs(title = "How each single-imputation method fills the gaps",
subtitle = "Grey: observed. Coloured: imputed",
x = "Predictor x", y = "Response y") + theme_te()
```
## What the fills do to the estimate
Estimate the slope, the correlation, and the spread of `y` on each completed dataset.
```{r}
#| label: point-estimates
estim <- function(yo) { ok <- !is.na(yo); f <- lm(yo[ok] ~ x[ok])
c(slope = unname(coef(f)[2]), r = cor(x[ok], yo[ok]), sdY = sd(yo[ok])) }
set.seed(29471)
tab <- rbind(
Full = estim(y),
`Complete-case` = { z <- y; z[miss] <- NA; estim(z) },
`Mean imputation` = estim(impute("mean")),
`Regression imputation` = estim(impute("detreg")),
`Stochastic imputation` = estim(impute("stoch")))
round(tab, 3)
```
The full data give a slope of `r round(full_slope, 3)`, a correlation of `r round(full_r, 3)`, and a standard deviation of `r round(full_sdY, 3)`. Mean imputation drags all three down: the slope falls to `r round(tab["Mean imputation","slope"], 3)`, the correlation to `r round(tab["Mean imputation","r"], 3)`, and the spread to `r round(tab["Mean imputation","sdY"], 3)`. Piling filled values onto a single horizontal line flattens the relationship and shrinks the variance, so every association is diluted.
Regression imputation does the opposite to the correlation. Its slope, `r round(tab["Regression imputation","slope"], 3)`, matches the complete-case fit, but the correlation climbs to `r round(tab["Regression imputation","r"], 3)`, above even the full-data value, because the imputed points sit perfectly on the line and manufacture agreement that is not in the data. Only stochastic imputation keeps the spread (`r round(tab["Stochastic imputation","sdY"], 3)`) and correlation (`r round(tab["Stochastic imputation","r"], 3)`) near their full-data values, because the added noise restores the scatter the other two erase.
So for an unbiased point estimate of the slope, mean imputation is out, and regression and stochastic imputation both land on the right answer. That would seem to settle it. It does not.
## The uncertainty that single imputation hides
Repeat the whole process 2000 times, recording each method's slope, its reported standard error, and whether the 95% interval covers the truth.
```{r}
#| label: monte-carlo
M <- 2000; keep <- c("cc", "mean", "detreg", "stoch")
sl <- se <- cv <- matrix(NA, M, 4, dimnames = list(NULL, keep))
set.seed(11072026)
for (i in 1:M) {
xx <- rnorm(n, 0, sdx); yy <- c0 + c1 * xx + rnorm(n, 0, sde)
a0i <- uniroot(function(a) mean(plogis(a + a1 * xx)) - 0.40, c(-6, 6))$root
mi <- runif(n) < plogis(a0i + a1 * xx); ob <- !mi
fit <- lm(yy[ob] ~ xx[ob]); sig <- summary(fit)$sigma
yh <- coef(fit)[1] + coef(fit)[2] * xx[mi]
comp <- list(cc = { z <- yy; z[mi] <- NA; z },
mean = { z <- yy; z[mi] <- mean(yy[ob]); z },
detreg = { z <- yy; z[mi] <- yh; z },
stoch = { z <- yy; z[mi] <- yh + rnorm(sum(mi), 0, sig); z })
for (k in keep) { z <- comp[[k]]; ok <- !is.na(z); f <- lm(z[ok] ~ xx[ok])
b <- unname(coef(f)[2]); s <- summary(f)$coef[2, 2]
sl[i, k] <- b; se[i, k] <- s
cv[i, k] <- (b - 1.96 * s <= c1) & (c1 <= b + 1.96 * s) }
}
summ <- rbind(bias = colMeans(sl) - c1, emp_SD = apply(sl, 2, sd),
mean_SE = colMeans(se), coverage = colMeans(cv))
round(summ, 3)
```
Complete-case analysis is honest: bias `r sprintf("%.3f", summ["bias","cc"])`, and its reported error (`r round(summ["mean_SE","cc"], 3)`) matches the true spread of estimates (`r round(summ["emp_SD","cc"], 3)`), so coverage is `r round(100 * summ["coverage","cc"])`%. Mean imputation fails on bias, `r sprintf("%.3f", summ["bias","mean"])`, and never covers.
The two unbiased imputations fail on a subtler count. Regression imputation reports an average error of `r round(summ["mean_SE","detreg"], 3)`, roughly half the true spread of `r round(summ["emp_SD","detreg"], 3)`, so its intervals cover only `r round(100 * summ["coverage","detreg"])`% of the time. Stochastic imputation is better but not fixed: its reported error `r round(summ["mean_SE","stoch"], 3)` still falls short of the true `r round(summ["emp_SD","stoch"], 3)`, and coverage sits at `r round(100 * summ["coverage","stoch"])`%. Figure 2 puts the reported error against the true spread for each method.
```{r}
#| label: fig-se
#| fig-cap: "Empirical spread of the slope over 2000 datasets against the average standard error each method reports, with interval coverage. Where the reported error falls short of the true spread, coverage drops below 95%."
#| fig-alt: "A dumbbell plot by method. Complete-case has its reported SE and empirical SD nearly equal at about 95% coverage; regression and stochastic imputation report SEs well below their empirical SDs, with coverage well short of 95%."
#| fig-width: 7.5
#| fig-height: 4.4
labs_m <- c(cc = "Complete-case", mean = "Mean imputation",
detreg = "Regression imputation", stoch = "Stochastic imputation")
df2 <- data.frame(method = factor(labs_m, levels = rev(labs_m)),
emp_SD = summ["emp_SD", ], mean_SE = summ["mean_SE", ],
coverage = summ["coverage", ])
dfl <- pivot_longer(df2, c(emp_SD, mean_SE), names_to = "q", values_to = "v")
dfl$q <- factor(dfl$q, levels = c("emp_SD", "mean_SE"),
labels = c("Empirical SD of slope", "Mean reported SE"))
ggplot(df2, aes(y = method)) +
geom_segment(aes(x = mean_SE, xend = emp_SD, yend = method), colour = "#b9b4a8", linewidth = 1.1) +
geom_point(data = dfl, aes(x = v, colour = q), size = 3) +
geom_text(aes(x = pmax(emp_SD, mean_SE) + 0.006,
label = sprintf("%.0f%% cover", 100 * coverage)),
hjust = 0, size = 3.3, colour = ink) +
scale_colour_manual(values = c("Empirical SD of slope" = "#2b2b2b", "Mean reported SE" = "#a24b4b")) +
coord_cartesian(xlim = c(0.03, 0.135)) +
labs(title = "A right point estimate can still report the wrong uncertainty",
subtitle = "Reported error below the true spread means intervals undercover",
x = "Standard error of the slope", y = NULL, colour = NULL) + theme_te()
```
## The honest limit
Mean imputation biases the estimate; regression imputation manufactures correlation; both understate error badly. Stochastic imputation clears the bias and restores the variance, yet still reports too little uncertainty. The reason is structural, not a matter of tuning. A single completed dataset treats the filled values as if they were the true ones, so the standard error reflects only the sampling variability of that one dataset, and not the extra uncertainty from having guessed the missing values in the first place. No single imputation, however clever the fill, can express that second layer of doubt, because it commits to one set of guesses.
The fix is to stop pretending there is one right fill. Multiple imputation draws several completed datasets, analyses each, and then combines the results so that the disagreement between them feeds back into the standard error. That is the subject of the next post.
## References
Little RJA, Rubin DB 2019. Statistical Analysis with Missing Data, 3rd edn. Wiley. ISBN 978-0-470-52679-8.
Enders CK 2010. Applied Missing Data Analysis. Guilford Press. ISBN 978-1-60623-639-0.
Nakagawa S, Freckleton RP 2008. Trends in Ecology and Evolution 23(11):592-596 (10.1016/j.tree.2008.06.014).
Rubin DB 1987. Multiple Imputation for Nonresponse in Surveys. Wiley. ISBN 978-0-471-08705-2.
Schafer JL 1997. Analysis of Incomplete Multivariate Data. Chapman and Hall/CRC. ISBN 978-0-412-04061-0.
## Related tutorials
- [Missing data: MCAR, MAR and MNAR](../missing-data-mechanisms-mcar-mar-mnar/)
- [Multiple imputation by chained equations](../multiple-imputation-chained-equations/)
- [Checking missing-data assumptions](../checking-missing-data-assumptions/)
- [Gap filling a flux time series](../gap-filling-a-flux-time-series/)
