---
title: "AIC with missing values: the silent row drop"
description: "R fits every model to whatever rows are complete for it, so an AIC table over covariates with different missingness ranks the gaps, not the ecology. Shown in R."
date: "2026-08-06 09:00"
date-modified: "2026-10-14"
categories: [R, model selection, AIC, missing data, ecology tutorial]
image: thumbnail.png
image-alt: "Four labelled points joined by a line, showing delta AIC rising steadily with the number of plots each candidate model was actually fitted to."
---
*Corrected 14 October 2026: the post said R runs such a comparison without complaint; it now says that AIC called on several models at once warns when they were fitted to different numbers of observations, though a table built one model at a time does not. It also no longer says that comparing row counts catches every case, since two models can drop the same number of different plots, and it ties densities above one to a small spread in the recorded units rather than to small units.*
*Corrected 1 October 2026: a section heading still called the direction of the artefact fixed without its condition, and now says it is fixed here; the post also said the four models separate only once multiplied by their row counts, and it now says most of the separation comes from that, since their fits per observation differ a little too.*
*Corrected 30 September 2026: the post said the direction of the row-drop artefact is fixed and general; it now says the direction holds when every row adds a negative log-likelihood term, as for counts and this response, and can reverse for a continuous response with densities above one.*
*Corrected 25 September 2026: the honest limits said complete-case analysis is unbiased only when values are missing completely at random, which is too strict for a regression. A new section measures what averaging AIC across imputed data sets does, and what to pool instead.*
Field data has holes. The soil chemistry came back for most plots but not all, the canopy measurement was skipped when the weather turned, one logger failed in July. You build a candidate model set out of the covariates you have, run `AIC`, and read off the winner.
R fits every candidate without complaint, and the answer can be entirely an artefact of where the holes are. Called on several models at once, `AIC` does warn that they were not all fitted to the same number of observations, but it still returns the scores, and a table that calls it on one model at a time, as below, prints no warning at all. Each model is fitted to the rows that are complete for the variables in that model, so different candidates are fitted to different data sets. Their log-likelihoods are sums over different numbers of rows, which makes them incomparable, and whenever each row adds a negative term, as it does for counts and for the response here, the incomparability points in a fixed direction: the model with the most missing data gets the best score.
## A candidate set where the answer is known
Two hundred plots. The response depends on `x1` and nothing else. `x2` is a covariate with no effect at all that is missing for thirty per cent of plots; `x3` is another null covariate missing for ten per cent.
```{r setup}
#| message: false
#| warning: false
library(ggplot2)
te_paper <- "#f5f4ee"
te_ink <- "#16241d"
te_body <- "#2c3a31"
te_forest <- "#275139"
te_rust <- "#b5534e"
te_gold <- "#c9b458"
te_line <- "#dad9ca"
theme_datasheet <- function() {
theme_minimal(base_size = 12) +
theme(plot.background = element_rect(fill = te_paper, colour = NA),
panel.background = element_rect(fill = te_paper, colour = NA),
panel.grid.major = element_line(colour = te_line, linewidth = 0.3),
panel.grid.minor = element_blank(),
text = element_text(colour = te_body),
plot.title = element_text(colour = te_ink, face = "bold"),
axis.text = element_text(colour = te_body))
}
make_plots <- function(n = 200, miss2 = 0.30, miss3 = 0.10) {
d <- data.frame(x1 = rnorm(n), x2 = rnorm(n), x3 = rnorm(n))
d$y <- 1.5 * d$x1 + rnorm(n) # x2 and x3 do nothing
d$x2[sample(n, round(miss2 * n))] <- NA
d$x3[sample(n, round(miss3 * n))] <- NA
d
}
set.seed(20260806)
plots <- make_plots()
colSums(is.na(plots))
```
The four candidates are the true model, the true model plus each null covariate, and the full model. Nothing here is subtle: `x2` and `x3` are noise, and any honest procedure should put `y ~ x1` first.
```{r candidates}
forms <- list(true = y ~ x1,
plus_x3 = y ~ x1 + x3,
plus_x2 = y ~ x1 + x2,
full = y ~ x1 + x2 + x3)
fits <- lapply(forms, lm, data = plots)
aic_tab <- data.frame(model = names(forms),
rows = sapply(fits, nobs),
aic = sapply(fits, AIC))
aic_tab$delta <- aic_tab$aic - min(aic_tab$aic)
aic_tab[order(aic_tab$aic), ]
```
The full model wins, by `r sprintf("%.1f", aic_tab$delta[aic_tab$model == "true"])` AIC units over the model that generated the data. On the usual reading that is overwhelming support: a difference of ten is normally treated as decisive, and this is more than twenty times that.
Look at the `rows` column. The four models were fitted to `r aic_tab$rows[aic_tab$model == "true"]`, `r aic_tab$rows[aic_tab$model == "plus_x3"]`, `r aic_tab$rows[aic_tab$model == "plus_x2"]` and `r aic_tab$rows[aic_tab$model == "full"]` plots respectively. The AIC ordering is the row-count ordering read backwards. No warning was printed, since each AIC was computed on its own, and nothing in the printed table would tell a reader which models saw which plots unless they asked.
```{r fig-mirror}
#| fig-cap: "Delta AIC for the four candidates beside the number of plots each was actually fitted to."
#| fig-alt: "Two horizontal bar panels sharing the same four model labels. In the left panel the delta AIC bars grow as the row count in the right panel grows, so the best-scoring model is the one fitted to the fewest plots."
shown <- c(true = "y ~ x1", plus_x3 = "y ~ x1 + x3",
plus_x2 = "y ~ x1 + x2", full = "y ~ x1 + x2 + x3")
mirror <- rbind(
data.frame(model = aic_tab$model, panel = "delta AIC", value = aic_tab$delta),
data.frame(model = aic_tab$model, panel = "plots used", value = aic_tab$rows))
mirror$model <- factor(shown[mirror$model],
levels = shown[aic_tab$model[order(aic_tab$rows)]])
mirror$panel <- factor(mirror$panel, levels = c("delta AIC", "plots used"))
ggplot(mirror, aes(y = model, x = value, fill = panel)) +
geom_col(width = 0.6) +
geom_text(aes(label = round(value, 1)), hjust = -0.15, size = 3.4,
colour = te_body) +
facet_wrap(~ panel, scales = "free_x") +
scale_x_continuous(expand = expansion(mult = c(0, 0.22))) +
scale_fill_manual(values = c("delta AIC" = te_rust, "plots used" = te_forest)) +
labs(x = NULL, y = NULL, fill = NULL,
title = "The ranking is the missingness, upside down") +
theme_datasheet() +
theme(legend.position = "none",
strip.text = element_text(colour = te_ink, face = "bold"))
```
## Why the direction is fixed here
The log-likelihood of a fitted model is a sum with one term per observation, and each term is negative for a continuous response with a density below one. Drop rows and you drop terms, so the log-likelihood rises towards zero and the AIC falls. The penalty for parameters is two per parameter, which is nowhere near enough to compensate for losing sixty rows.
The arithmetic is worth seeing once, because it makes clear that this is not a subtle statistical effect but a bookkeeping error.
```{r arithmetic}
per_row <- sapply(fits, function(m) as.numeric(logLik(m)) / nobs(m))
round(rbind(total_loglik = sapply(fits, function(m) as.numeric(logLik(m))),
rows = sapply(fits, nobs),
per_row = per_row), 3)
```
Per observation the four models fit much alike, spanning `r sprintf("%.3f", diff(range(per_row)))` of a log-likelihood unit. Most of the separation comes from multiplying each by its own row count, and that part is in proportion to how much data each one threw away.
## Refit on the same plots and the ordering reverses
The repair is to decide the data set before the model set. Take the rows that are complete for every variable that appears anywhere in the candidate list, and fit all four to those.
```{r complete}
vars <- c("y", "x1", "x2", "x3")
common <- plots[complete.cases(plots[, vars]), ]
fits_c <- lapply(forms, lm, data = common)
tab_c <- data.frame(model = names(forms),
rows = sapply(fits_c, nobs),
aic = sapply(fits_c, AIC))
tab_c$delta <- tab_c$aic - min(tab_c$aic)
tab_c[order(tab_c$aic), ]
```
All four are now fitted to the same `r tab_c$rows[1]` plots, and the true model is first with the two single-additions within `r sprintf("%.1f", max(tab_c$delta[tab_c$model %in% c("plus_x2", "plus_x3")]))` AIC units and the full model `r sprintf("%.1f", tab_c$delta[tab_c$model == "full"])` behind. That is the pattern AIC is supposed to produce for a null covariate: a penalty of about two per useless parameter, and no support for the elaboration.
The cost is visible too. The analysis now uses `r tab_c$rows[1]` of the original `r nrow(plots)` plots, because the union of the missingness patterns is larger than either alone. That is a real loss and it is the subject of a different decision, but it is a loss you can see and argue about, rather than one that quietly rewrites the ranking.
## How far the problem scales
The size of the artefact is a function of how much data the extra covariate is missing, and it is monotone.
```{r sweep}
sweep_miss <- function(fracs, reps = 40) {
out <- data.frame()
for (f in fracs) {
gaps <- replicate(reps, {
d <- make_plots(miss2 = f, miss3 = 0)
AIC(lm(y ~ x1, data = d)) - AIC(lm(y ~ x1 + x2, data = d))
})
out <- rbind(out, data.frame(missing = f, gap = mean(gaps),
sd = sd(gaps),
wrong = mean(gaps > 0)))
}
out
}
set.seed(2)
sw <- sweep_miss(c(0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6))
round(sw, 3)
```
With no missing values at all the true model is ahead by `r sprintf("%.1f", -sw$gap[sw$missing == 0])` AIC units on average: the two-unit penalty for the extra parameter, offset by the improvement in fit a null covariate buys by chance. The elaboration still wins `r sprintf("%.0f", 100 * sw$wrong[sw$missing == 0])` per cent of the time on that margin, which is ordinary AIC behaviour and the reason a gap of one or two units settles nothing.
Then the missingness starts. At a tenth of plots missing the elaboration already wins every one of the `r sw$wrong[sw$missing == 0.1] * 40` simulated data sets, by an average of `r sprintf("%.0f", sw$gap[sw$missing == 0.1])` AIC units. At sixty per cent missing the margin is `r sprintf("%.0f", sw$gap[sw$missing == 0.6])` units in favour of a covariate that does nothing.
```{r fig-sweep}
#| fig-cap: "Mean AIC advantage of a useless covariate over the true model, as a function of how much of that covariate is missing. Forty simulated data sets per point."
#| fig-alt: "A rising curve with error bars, starting slightly below the zero line when nothing is missing and climbing steeply above it, so a covariate with no effect gains hundreds of AIC units of apparent support purely from being missing more often."
ggplot(sw, aes(x = 100 * missing, y = gap)) +
geom_hline(yintercept = 0, colour = te_ink, linewidth = 0.5) +
geom_errorbar(aes(ymin = gap - sd, ymax = gap + sd), width = 2.2,
linewidth = 0.4, colour = te_forest) +
geom_line(linewidth = 0.9, colour = te_forest) +
geom_point(size = 2.8, colour = te_forest) +
annotate("text", x = 1, y = max(sw$gap), hjust = 0, vjust = 1, size = 3.4,
colour = te_ink,
label = "everything above the dark line is\nsupport for a covariate with no effect") +
labs(x = "per cent of plots missing the useless covariate",
y = "AIC advantage of the useless covariate",
title = "Support bought with absent data") +
theme_datasheet()
```
## Making R refuse
The habit that prevents all of this is one line at the top of an analysis script. Setting the default action on missing values to `na.fail` turns the silent drop into an error, so a model that cannot use every row of the frame you handed it stops rather than guessing.
```{r nafail}
old <- options(na.action = "na.fail")
attempt <- try(lm(y ~ x1 + x2, data = plots), silent = TRUE)
cat(class(attempt), ":", conditionMessage(attr(attempt, "condition")), "\n")
options(old)
```
With that set, the workflow becomes: build the analysis frame, subset it to complete cases over every variable in the candidate set, then fit. The models cannot disagree about their data because the frame no longer contains any disagreement.
A lighter check, if you would rather not change a global option, is to compare `nobs` across the fitted models before touching the AIC table. One line, and it catches every difference in row counts, though not two models that drop the same number of different plots; comparing the row names of each model's `model.frame` catches those too.
```{r nobscheck}
same_rows <- function(models) length(unique(sapply(models, nobs))) == 1
c(as_fitted = same_rows(fits), on_common_frame = same_rows(fits_c))
```
## If you impute: pool the test, not the AIC
Multiple imputation removes the row-count problem, since every completed data set contains every plot. It does not make AIC usable as it stands. Meng and Rubin (1992) showed how likelihood ratio tests can be combined across imputed data sets, and Wood, White and Royston (2008) compared ways of selecting variables after imputation and found that only selection based on Rubin's rules kept the type 1 error where it should be. The shortcut that is easy to reach for is different: fit each candidate to each completed data set and average the AIC.
The check below has 150 plots per study, `x1` and `x2` correlated at 0.5, and `y = 1.5 x1 + e` as before, so `x2` does nothing. `x2` goes missing more often on plots with high `x1`, which is missing at random. Each study gets 10 proper normal-linear imputations of `x2` from `x1` and `y`, built the same way as in [the chained equations post](../multiple-imputation-chained-equations/). Three rules then decide whether to keep `x2`: AIC on the common complete-case frame (the repair above), the AIC difference averaged over the 10 imputations, and the pooled Wald statistic from Rubin's rules compared with 2, which is the threshold AIC sets for one extra parameter. AIC on the data before any values were deleted is the reference. A second set of studies gives `x2` a real coefficient of 0.2, to measure power.
```{r pooltest}
mi_n <- 150
mi_m <- 10
impute_x2 <- function(x2, mi_miss, predictors) {
obs_rows <- predictors[!mi_miss, ]
qr_fit <- lm.fit(obs_rows, x2[!mi_miss])
dfree <- sum(!mi_miss) - ncol(predictors)
s2star <- sum(qr_fit$residuals^2) / rchisq(1, dfree)
bstar <- qr_fit$coefficients +
drop(t(chol(s2star * chol2inv(chol(crossprod(obs_rows))))) %*% rnorm(ncol(predictors)))
x2[mi_miss] <- drop(predictors[mi_miss, ] %*% bstar) + rnorm(sum(mi_miss), 0, sqrt(s2star))
x2
}
loglik_lm <- function(design, y) {
res <- lm.fit(design, y)$residuals
-length(y) / 2 * (log(2 * pi * mean(res^2)) + 1)
}
keeps_x2 <- function(design, y) {
(-2 * loglik_lm(design, y) + 2 * 4) < (-2 * loglik_lm(design[, 1:2], y) + 2 * 3)
}
one_study <- function(b2) {
x1 <- rnorm(mi_n)
x2 <- 0.5 * x1 + sqrt(0.75) * rnorm(mi_n)
y <- 1.5 * x1 + b2 * x2 + rnorm(mi_n)
mi_miss <- runif(mi_n) < plogis(-0.55 + 1.5 * x1)
full_design <- cbind(1, x1, x2)
delta_aic <- est <- within_var <- numeric(mi_m)
for (k in seq_len(mi_m)) {
design_k <- cbind(1, x1, impute_x2(x2, mi_miss, cbind(1, x1, y)))
delta_aic[k] <- (-2 * loglik_lm(design_k, y) + 8) - (-2 * loglik_lm(design_k[, 1:2], y) + 6)
fit_k <- lm.fit(design_k, y)
est[k] <- fit_k$coefficients[3]
within_var[k] <- sum(fit_k$residuals^2) / (mi_n - 3) * chol2inv(qr.R(fit_k$qr))[3, 3]
}
between <- var(est)
total <- mean(within_var) + (1 + 1 / mi_m) * between
c(full = keeps_x2(full_design, y),
common_cc = keeps_x2(full_design[!mi_miss, ], y[!mi_miss]),
mean_aic = mean(delta_aic) < 0,
rubin = mean(est)^2 / total > 2,
r = between / mean(within_var),
missing = mean(mi_miss))
}
n_study <- 250
set.seed(1992)
mi_null <- replicate(n_study, one_study(0))
mi_power <- replicate(n_study, one_study(0.2))
mi_keep <- cbind(null = rowMeans(mi_null[1:4, ]), power = rowMeans(mi_power[1:4, ]))
mi_r <- mi_null["r", ]
mi_closed <- mean(pchisq((2 - mi_r * (mi_m - 1) / mi_m) / (1 + mi_r * (1 + 1 / mi_m)),
1, lower.tail = FALSE))
mi_keep_se <- sqrt(mi_keep * (1 - mi_keep) / n_study)
mi_ratio <- mi_keep["mean_aic", "null"] / mi_keep["full", "null"]
mi_missing <- mean(mi_null["missing", ])
round(t(mi_keep), 3)
round(c(closed_form = mi_closed, median_r = median(mi_r)), 3)
```
Across `r sprintf("%d", n_study)` studies with `r sprintf("%.0f", 100 * mi_missing)` per cent of `x2` missing on average, the averaged AIC keeps the useless `x2` in `r sprintf("%.1f", 100 * mi_keep["mean_aic", "null"])` per cent, `r sprintf("%.1f", mi_ratio)` times the `r sprintf("%.1f", 100 * mi_keep["full", "null"])` per cent that AIC gives on the complete data (the textbook rate for one extra parameter is `r sprintf("%.1f", 100 * pchisq(2, 1, lower.tail = FALSE))` per cent). The pooled Wald keeps it in `r sprintf("%.1f", 100 * mi_keep["rubin", "null"])` per cent and the common complete-case frame in `r sprintf("%.1f", 100 * mi_keep["common_cc", "null"])` per cent, with a Monte Carlo standard error of up to `r sprintf("%.1f", 100 * max(mi_keep_se[, "null"]))` points on each rate.
The inflation is arithmetic, and the simulation only checks it. Each imputation produces its own estimate of the `x2` coefficient, scattered around the pooled estimate with between-imputation variance `B`. The likelihood ratio statistic of one completed data set is close to that estimate squared over its within-imputation variance `U`, so the average over `M` imputations exceeds the pooled estimate's own statistic by about `v (M - 1) / M`, where `v = B / U`. Under the null the pooled estimate has variance `U (1 + v (1 + 1/M))`. Averaged AIC keeps `x2` when the average statistic is above 2, so it keeps it with probability `P(chi2_1 > (2 - v (M - 1) / M) / (1 + v (1 + 1/M)))`. With the `v` measured in each study (median `r sprintf("%.3f", median(mi_r))`) that probability averages `r sprintf("%.3f", mi_closed)`, against `r sprintf("%.3f", mi_keep["mean_aic", "null"])` in the simulation. The between-imputation variance enters twice, once as extra spread in the pooled estimate and once as a positive shift in every study; Rubin's rules put it in the denominator instead. Van Ginkel and Molenaar (2026) proposed pooled estimators of AIC and BIC for multiply imputed data and found, in factor, multilevel and latent class models, that they pick the correct model more often than averaging the criteria across imputed data sets.
With a real coefficient of 0.2, AIC on the complete data keeps `x2` in `r sprintf("%.1f", 100 * mi_keep["full", "power"])` per cent of studies. The pooled Wald keeps it in `r sprintf("%.1f", 100 * mi_keep["rubin", "power"])` per cent and the common complete-case frame in `r sprintf("%.1f", 100 * mi_keep["common_cc", "power"])` per cent (standard errors up to `r sprintf("%.1f", 100 * max(mi_keep_se[, "power"]))` points). The averaged AIC reaches `r sprintf("%.1f", 100 * mi_keep["mean_aic", "power"])` per cent only because it keeps `x2` too readily whether or not it has an effect. So in this design imputation buys no power for selecting the covariate that carries the gaps, and the common frame of the earlier section does as well as the imputation rule that holds its level. If you do impute, pool the test and not the AIC: Rubin's rules for one coefficient, or Meng and Rubin's combined likelihood ratio when a term has several parameters (not run here).
## Honest limits
Restricting to complete cases fixes the comparability problem and introduces a different one. For a regression it is unbiased when whether a value is missing does not depend on the response once the covariates in the model are known: missing completely at random, or missing because of `x1` as in the imputation section above, where the common frame kept its level. When the reason a value is absent is related to the response, dropping those rows biases the estimates, and the tidy AIC table is then a comparison of models fitted to a non-random subsample. Which repair is right depends on why the data are missing, which is not a question the data can answer on their own.
Multiple imputation is the usual alternative, and the likelihood of an imputed data set is not the likelihood of the observed one, so ordinary AIC does not transfer without modification. The section above measures one case: a single incomplete covariate, `lm`, one extra parameter and missingness at random given `x1`. Several incomplete covariates, `glm` and mixed models, and selection on stacked or weighted imputed data sets were not tested, and the imputation model uses `y` whether or not the candidate contains `x2`. None of this is a reason to prefer the silent drop, which has no defensible interpretation at all.
The demonstration uses `lm` on a continuous response. The mechanism is the same for `glm` and for mixed models, and in the mixed case there is a second version of the same trap: fits with different fixed-effect structures estimated by restricted maximum likelihood are also not comparable by AIC, for a related reason about what the likelihood is a likelihood of. [AIC on REML fits: when units choose the model](../aic-on-reml-fits/) measures that second trap.
Finally, the row-drop numbers here belong to one generator with two hundred plots, a moderate effect and independent missingness. The direction of the artefact holds whenever every row adds a negative log-likelihood term, as it does for counts and for this response (a continuous response whose spread, in the units it is recorded in, is small enough for its densities to exceed one would reverse it), and the size of it is not general; run the sweep on your own frame if you want to know how much of a published ranking could have come from this.
## References
Akaike H 1974 IEEE Transactions on Automatic Control 19(6):716-723 (10.1109/TAC.1974.1100705)
Nakagawa S, Freckleton RP 2008 Trends in Ecology and Evolution 23(11):592-596 (10.1016/j.tree.2008.06.014)
Burnham KP, Anderson DR 2002 Model Selection and Multimodel Inference, 2nd edition, Springer, ISBN 978-0-387-95364-9
Meng XL, Rubin DB 1992 Biometrika 79(1):103-111 (10.1093/biomet/79.1.103)
Wood AM, White IR, Royston P 2008 Statistics in Medicine 27(17):3227-3246 (10.1002/sim.3177)
van Ginkel JR, Molenaar D 2026 Behaviormetrika 53(2):459-490 (10.1007/s41237-025-00281-6)
## Related tutorials
- [Model selection and AIC](../model-selection-aic/)
- [Missing data: MCAR, MAR and MNAR](../missing-data-mechanisms-mcar-mar-mnar/)
- [AIC on REML fits: when units choose the model](../aic-on-reml-fits/)
- [Multiple imputation by chained equations](../multiple-imputation-chained-equations/)
