Chapter 7
Missing data

Lectures 13–14 · missingness mechanisms, imputation as a model, and how far the answer depends on what we cannot see

Real data sets have holes. A survey respondent skips the income question, a sensor drops out, a hospital stops reporting, a country has no system for registering deaths. The analyst is left with a table that is not rectangular, and the standard software response is to delete every row with a hole and carry on with the rows that survive. That reflex has a name, complete-case analysis, and this chapter is about when it is harmless, when it is wasteful and when it gives the wrong answer.

The Bayesian view of missing data is short enough to state in one sentence: a missing value is an unknown quantity, and unknown quantities get posterior distributions. A missing value becomes a parameter, with a prior given by a model for how such values are distributed, and is estimated along with everything else. Nothing new is needed: no special algorithm, just a model for the missing values. What is new is a question that no algorithm can answer for us: why are the values missing? As in every chapter so far, a causal graph is how we make that question precise.

• Lecture 13 is about missing outcomes. We will estimate how many people died worldwide because of the COVID-19 pandemic in 2020 and 2021, when excess deaths can be measured in only 86 of the 158 countries with more than a million people. We classify missingness with Rubin's taxonomy, see why the countries that report are a very unrepresentative sample, and build a regression model that imputes the missing countries, with honest uncertainty.
• Lecture 14 is about missing predictors and missing mechanisms. We impute missing predictors inside the model with brms's mi(), see why the popular shortcuts (filling in the mean, or a single predicted value) mislead, meet Rubin's multiple imputation, and finish with the case that no model can fix, missingness that depends on the missing value itself, and the sensitivity analysis that makes it visible.

Lecture 13 — Missing outcomes

7.1 How many people died?

Counting the dead of a pandemic sounds simple and is not. Official COVID-19 death counts depend on testing, on what doctors write on death certificates, and on whether deaths are registered at all, and all three vary enormously between countries. Demographers therefore use a different measure: excess mortality, the number of deaths from all causes beyond what would have been expected had there been no pandemic. Expected deaths come from projecting each country's pre-2020 trend in deaths forward. Excess mortality counts deaths from COVID-19 that were never diagnosed, and also deaths caused indirectly (delayed treatment, overwhelmed hospitals) or prevented (fewer traffic accidents and less influenza during lockdowns).

Excess mortality needs one thing that many countries do not have: a timely, complete count of all deaths. The World Mortality Dataset (Karlinsky and Kobak, 2021) collected weekly or monthly death counts from every country that publishes them, and Our World in Data distributes the resulting estimates of excess deaths. We combine them with reported COVID-19 deaths and a few country characteristics, and keep the 158 countries with more than a million people.

library(tidyverse)
library(brms)
library(posterior)
options(mc.cores = parallel::detectCores())

harbour <- c(ink = "#10191d", muted = "#4a5f68", grid = "#d5e3e6",
             sky = "#1d5187", amber = "#805907", rose = "#ad3b28",
             moss = "#1d6b54", lilac = "#5b3a94")
theme_harbour <- function(base_size = 11) {
  theme_minimal(base_size = base_size) +
    theme(text = element_text(colour = harbour[["ink"]]),
          axis.text = element_text(colour = harbour[["muted"]]),
          panel.grid.major = element_line(colour = harbour[["grid"]]),
          panel.grid.minor = element_blank(),
          strip.text = element_text(colour = harbour[["ink"]], hjust = 0),
          legend.position = "top")
}
theme_set(theme_harbour())

base <- "https://raw.githubusercontent.com/owid/covid-19-data/master/public/data/"
countries <- read_csv(paste0(base, "latest/owid-covid-latest.csv"), show_col_types = FALSE) |>
  filter(!startsWith(iso_code, "OWID"), population > 1e6) |>           # countries, not regions
  select(location, continent, population, median_age, gdp_per_capita, life_expectancy)
reported <- read_csv(paste0(base, "cases_deaths/total_deaths_per_million.csv"), show_col_types = FALSE) |>
  filter(date == as.Date("2021-12-31")) |>
  pivot_longer(-date, names_to = "location", values_to = "reported") |> select(-date)
excess <- read_csv(paste0(base, "excess_mortality/excess_mortality.csv"), show_col_types = FALSE) |>
  filter(date <= as.Date("2021-12-31")) |>
  group_by(location) |> slice_max(date, n = 1) |> ungroup() |>
  filter(date >= as.Date("2021-12-01")) |>                              # series must run to the end of 2021
  transmute(location, excess = cum_excess_per_million_proj_all_ages)

covid <- countries |> left_join(reported, by = "location") |> left_join(excess, by = "location") |>
  mutate(reported = na_if(reported, 0),       # "zero deaths" (North Korea, Turkmenistan) is not a measurement
         y  = excess / 1000,                  # excess deaths per 1,000 people, 2020-21
         r  = log(reported / 500),            # log reported COVID deaths per million, relative to 500
         a  = (median_age - 30) / 10,         # median age, in decades from 30
         g  = log(gdp_per_capita / 1e4),      # log GDP per person, relative to $10,000
         le = (life_expectancy - 70) / 10)    # life expectancy, in decades from 70
c(countries = nrow(covid), excess_observed = sum(!is.na(covid$y)))
covid |> group_by(continent) |>
  summarise(countries = n(), observed = sum(!is.na(y)),
            `% of population observed` = round(100 * sum(population[!is.na(y)]) / sum(population)))

      countries excess_observed 
            158              86 
# A tibble: 6 x 4
  continent     countries observed `% of population observed`
  <chr>             <int>    <int>                      <dbl>
1 Africa               50        6                         16
2 Asia                 45       25                         15
3 Europe               35       34                         99
4 North America        15       10                         94
5 Oceania               3        2                         76
6 South America        10        9                         93

Excess deaths through the end of 2021 are available for 86 of the 158 countries. The holes are not scattered at random. Almost all of Europe and the Americas is observed (over 90% of the population of each continent); only 6 of the 50 African countries are, covering about a sixth of Africa's population; and in Asia, the observed countries hold less than a sixth of the population, because China, India, Indonesia, Pakistan and Bangladesh are all missing. Here is a first look at what distinguishes the countries we see from the countries we do not.

covid <- covid |> mutate(status = if_else(is.na(y), "missing", "observed"))
p1 <- ggplot(covid, aes(median_age, gdp_per_capita, size = population, shape = status, colour = status)) +
  geom_point(alpha = 0.8, stroke = 0.8) +
  scale_y_log10(labels = scales::label_dollar()) + scale_size_area(max_size = 12, guide = "none") +
  scale_shape_manual(values = c(missing = 1, observed = 16)) +
  scale_colour_manual(values = c(missing = harbour[["rose"]], observed = harbour[["sky"]])) +
  labs(x = "median age (years)", y = "GDP per person", shape = NULL, colour = NULL)
p2 <- ggplot(filter(covid, !is.na(y)), aes(reported, excess)) +
  geom_line(data = tibble(reported = 10^seq(1, log10(7000), length.out = 200), excess = reported),     # excess = reported
            linetype = "dashed", colour = harbour[["muted"]]) +
  geom_hline(yintercept = 0, colour = harbour[["grid"]]) +
  geom_point(colour = harbour[["sky"]], alpha = 0.8) +
  scale_x_log10() +
  labs(x = "reported COVID-19 deaths per million", y = "excess deaths per million")
patchwork::wrap_plots(p1, p2, widths = c(1.1, 1))

The left panel is the problem of this lecture in one picture. Countries that can measure excess mortality are older and richer than countries that cannot; the median observed country has a median age of 38 and a GDP of about $23,000 per person, the median missing country a median age of 20 and about $3,600. The right panel shows why reported deaths are not enough: in almost every country, excess deaths exceeded reported COVID-19 deaths, sometimes by a factor of three or more (Russia and Serbia, for example), while a few countries that controlled the virus early, such as New Zealand, Japan and Australia, had about as many deaths as usual or fewer.

Which brings us to the obvious estimate of the world total, and why it is wrong. Take the population-weighted average excess death rate in the observed countries and multiply it by the world's population:

obs <- filter(covid, !is.na(y))
c(reported_deaths_millions = sum(covid$reported * covid$population, na.rm = TRUE) / 1e12,
  excess_in_observed_millions = sum(obs$y * obs$population) / 1e9,
  naive_world_millions = weighted.mean(obs$y, obs$population) * sum(covid$population) / 1e9)

   reported_deaths_millions excess_in_observed_millions 
                      5.429                       7.404 
       naive_world_millions 
                     21.815 

The 86 observed countries alone had about 7.4 million excess deaths. Treating the other 72 as if they were like the observed ones gives almost 22 million. But they are not like the observed ones: they are much younger, and COVID-19 kills old people at hundreds of times the rate it kills the young. The naive estimate assumes the missing countries are a random sample of all countries, and the left panel of the figure shows how badly that assumption fails.

7.2 Why values go missing

Whether complete-case analysis is harmless depends on why values are missing, and that is a causal question. Add to the causal graph a missingness indicator $R_Y$: a variable that equals 1 if $Y$ was observed and 0 if not. The question is what causes $R_Y$. Donald Rubin's 1976 taxonomy distinguishes three cases, drawn below for an outcome $Y$ and a fully observed variable $X$ that affects it.

• Missing completely at random (MCAR). Missingness is unrelated to anything in the problem: a lab sample was dropped, a page of a questionnaire was lost in the post. The observed rows are a random sample of all rows. Complete-case analysis is unbiased, but it wastes the information in the incomplete rows.
• Missing at random (MAR). Missingness depends only on variables we observe. Here, whether a country can measure excess deaths depends on its income and the age of its population, which we know. The observed rows are a biased sample, so a simple average of them is biased, but the bias is repairable: comparing countries with the same income and age, the observed and missing ones are alike, so a model that conditions on income and age can predict the missing ones.
• Missing not at random (MNAR). Missingness depends on the missing value itself, or on something unobserved that also affects it. If the countries hardest hit by the pandemic were also the ones whose registration systems broke down, then even among countries with the same income and age, the missing ones would have higher excess deaths than the observed ones. No model of the observed data alone can detect or remove this bias.

The names are confusing (MAR does not mean missingness is random), but the graphs are not. Under MCAR and MAR, the missingness mechanism can be ignored when we fit a model that conditions on the variables that drive it; Rubin called such mechanisms ignorable. Under MNAR it cannot. Whether a mechanism is MAR or MNAR cannot be tested from the data, because the data say nothing about the values we do not have. It has to be argued from knowledge of how the data were collected, exactly like the arrows of any other causal graph.

Different estimates are affected differently, and the widget makes this concrete. It uses the life expectancy and GDP per person of the 154 countries with more than a million people where both are known, deletes some life expectancies according to a mechanism you choose, and compares two quantities computed from the complete cases with their true values: the average life expectancy, and the slope of life expectancy on log GDP.

Three lessons come out of a few clicks. Under MCAR, the complete cases give unbiased answers for both the mean and the slope; they are just noisier. Under MAR (missingness depending on GDP), the complete-case mean is biased upward, because poor countries with low life expectancy are missing, but the complete-case slope is still right: among countries with the same GDP, those that are observed are a random sample, so the regression of $Y$ on $X$ is unaffected. And regression imputation of the mean repairs the bias, because it uses the observed relationship to predict the missing countries from their GDP. Under MNAR, everything is biased, including the slope and the imputation. That is the situation the world excess-mortality estimate must hope to avoid.

7.3 Missing outcomes: imputing the world

If the mechanism is MAR given some observed country characteristics, then a regression of excess deaths on those characteristics, fitted to the observed countries, is unbiased (the slope lesson of the widget), and its predictions for the missing countries are the imputations we need. Which characteristics? Three are available for almost every country and plausibly affect both excess deaths and whether they are measured:

• reported COVID-19 deaths per million, which, imperfect as they are, track the size of each country's epidemic;
• median age, because the risk of death from COVID-19 rises steeply with age;
• GDP per person, which affects health care, the ability to work from home and the completeness of reporting (in poorer countries, fewer deaths are tested and reported, so the same number of reported deaths implies more real ones).

The model is a linear regression of excess deaths per 1,000 people on log reported deaths, median age and log GDP. The predictors are centred so the intercept is the excess death rate for a country with 500 reported deaths per million, a median age of 30 and GDP of $10,000 per person. We set aside the six countries with a missing predictor (Section 7.5 brings them back).

known <- filter(covid, !is.na(r), !is.na(g))     # countries with complete predictors
m_cc <- brm(y ~ r + a + g, data = filter(known, !is.na(y)),
            prior = c(prior(normal(2, 2), class = Intercept),   # excess per 1,000 for the reference country
                      prior(normal(0, 2), class = b),
                      prior(exponential(1), class = sigma)),
            seed = 7, refresh = 0, file = "fits/ch7_cc")
fixef(m_cc) |> round(2)
posterior_summary(m_cc, variable = "sigma") |> round(2)

          Estimate Est.Error  Q2.5 Q97.5
Intercept     2.61      0.22  2.19  3.05
r             0.78      0.12  0.55  1.00
a             0.78      0.29  0.21  1.34
g            -1.41      0.25 -1.92 -0.94
      Estimate Est.Error Q2.5 Q97.5
sigma     1.41      0.11  1.2  1.65

All three predictors matter, in the expected directions. Each doubling of reported deaths adds about $0.78\times \log 2\approx 0.5$ excess deaths per 1,000 people. Ten more years of median age, at the same reported deaths and income, adds about 0.8 per 1,000. And each doubling of GDP per person lowers excess deaths by about $1.4\times \log 2\approx 1$ per 1,000 at the same reported deaths and age: richer countries reported a larger share of their real deaths, and probably prevented more. The residual standard deviation is about 1.4 deaths per 1,000, which is large: countries with the same characteristics differ a lot, and any single country's prediction will be uncertain.

Now the imputation itself. For each missing country, draw its excess death rate from the posterior predictive distribution, multiply by its population and add up. Each posterior draw is one completed world; the spread of the totals across draws is our uncertainty about the world total, and it includes both our uncertainty about the regression coefficients and the country-to-country scatter.

missing <- filter(known, is.na(y))
set.seed(7)
draws <- posterior_predict(m_cc, newdata = missing)            # draws x missing countries
total_missing <- draws %*% missing$population / 1e9            # millions of deaths, per draw
total_observed <- with(filter(known, !is.na(y)), sum(y * population)) / 1e9
world <- total_observed + total_missing
c(observed = total_observed, quantile(total_missing, c(.5, .025, .975)))
quantile(world, c(.5, .025, .975))

observed      50%     2.5%    97.5% 
 7.34739  7.22178  0.04351 14.06476 
   50%   2.5%  97.5% 
14.569  7.391 21.412 

big12 <- missing |> slice_max(population, n = 12)
pp <- posterior_predict(m_cc, newdata = big12)
q <- big12 |> mutate(med = apply(pp, 2, median), lo = apply(pp, 2, quantile, .025),
                     hi = apply(pp, 2, quantile, .975), rep1000 = reported / 1000)
p1 <- ggplot(q, aes(y = fct_reorder(location, population))) +
  geom_vline(xintercept = 0, colour = harbour[["grid"]]) +
  geom_errorbarh(aes(xmin = lo, xmax = hi), height = 0, colour = harbour[["sky"]]) +
  geom_point(aes(x = med), colour = harbour[["sky"]], size = 2) +
  geom_point(aes(x = rep1000), shape = 4, colour = harbour[["rose"]], size = 2) +
  labs(x = "excess deaths per 1,000", y = NULL)
marks <- tibble(x = c(5.4, 14.9, 18.2, 21.8), lab = c("reported", "WHO", "IHME", "naive"))
p2 <- ggplot(tibble(world = as.numeric(world)), aes(world)) +
  geom_histogram(binwidth = 0.5, fill = harbour[["sky"]], alpha = 0.6) +
  geom_vline(data = marks, aes(xintercept = x), colour = harbour[["amber"]], linetype = "dashed") +
  geom_text(data = marks, aes(x = x, y = Inf, label = lab), angle = 90, hjust = 1.1, vjust = -0.4,
            colour = harbour[["amber"]], size = 3.2) +
  labs(x = "world excess deaths (millions)", y = "posterior draws")
patchwork::wrap_plots(p1, p2, widths = c(1, 1.1))

The posterior median for the world total is about 14.6 million excess deaths in 2020–21, with a 95% interval of roughly 7.4 to 21 million. That is nearly three times the reported COVID-19 deaths, and far below the naive extrapolation. It sits close to the World Health Organization's published estimate of 14.9 million (Msemburi et al., 2023) and somewhat below the estimate of 18.2 million from the Institute for Health Metrics and Evaluation (IHME; COVID-19 Excess Mortality Collaborators, 2022). Our interval is much wider than either, because those teams also used subnational data, surveys and much richer covariates. A model with three predictors cannot match that; the point is that it lands in the right place, for the right reason, and says honestly how uncertain it is.

The country-level imputations are where the action is. India alone contributes about 4 million excess deaths in our model (about 2.8 per 1,000 people), against about 480,000 reported. The WHO's estimate for India, about 4.7 million, was rejected by the Indian government when it was published in 2022 and became one of the most contested numbers of the pandemic; independent studies using Indian survey and death registration data (for example Jha et al., 2022) also found several times the official count. China is the opposite case: it reported almost no COVID-19 deaths in 2020–21, far outside the range of any observed country, and the model extrapolates to a negative excess (fewer deaths than usual), with a wide interval. Section 7.4 asks whether we should believe that.

7.4 Checking an imputation by extrapolating

The imputations for the missing countries are predictions for countries unlike any we have observed: younger, poorer, and with fewer reported deaths. We cannot check those predictions directly, but we can check the model's ability to extrapolate in the same direction. Split the observed countries by income, fit the model to the richer half only, and predict the poorer half, whose true values we know.

obs_k <- filter(known, !is.na(y))
rich <- filter(obs_k, g >= median(g)); poor <- filter(obs_k, g < median(g))
m_rich <- update(m_cc, newdata = rich, seed = 7, refresh = 0)
set.seed(7)
pp_poor <- posterior_predict(m_rich, newdata = poor)
tibble(countries = nrow(poor),
       mean_error = mean(colMeans(pp_poor) - poor$y),                       # predicted minus actual
       rmse = sqrt(mean((colMeans(pp_poor) - poor$y)^2)),
       coverage_95 = mean(poor$y > apply(pp_poor, 2, quantile, .025) &
                          poor$y < apply(pp_poor, 2, quantile, .975)),
       total_error_millions = sum((colMeans(pp_poor) - poor$y) * poor$population) / 1e9)

# A tibble: 1 x 5
  countries mean_error  rmse coverage_95 total_error_millions
      <int>      <dbl> <dbl>       <dbl>                <dbl>
1        42       1.04  2.06       0.857                 1.22

The check is sobering. Fitted to the 42 richer observed countries, the model over-predicts the 42 poorer ones by about 1 excess death per 1,000 on average, with a root mean squared error of about 2 per 1,000, and its 95% predictive intervals contain only 86% of the true values instead of 95%. Extrapolated toward poorer countries, the linear model overshoots, and it is overconfident. That does not mean the imputations for the missing countries are too high by the same amount: the full model is fitted to the poorer observed countries too, so it does not extrapolate as far. But the missing countries are poorer and younger still, and the check is a warning that the model's form, not just its parameters, is being stretched. It points toward a world total somewhat lower than our 14.6 million; the published estimates, which used much more information, came out at or above ours. When two honest checks pull in opposite directions, the right response is the sensitivity analysis of Section 7.7, not a confident single number.

This kind of check does not prove that the imputations for Africa and South Asia are right; nothing observed can. It does tell us whether the model's form survives extrapolation in the relevant direction, and whether its intervals are honest when it does. China is the case that most needs scrutiny: its reported death rate (4 per million) is less than half the lowest observed (New Zealand's, about 10), and its political circumstances (strict elimination of the virus throughout 2020–21, and doubts about the reported numbers) are unlike any observed country's. A careful analyst would report the world total with and without the model's extrapolation for China, as the widget in Section 7.7 lets you do.

Lecture 14 — Missing predictors and missing mechanisms

7.5 Missing predictors: imputation inside the model

Six of the 158 countries were left out of the analysis so far because a predictor was missing: Cuba, Somalia, Syria and North Korea have no GDP figure in the data; North Korea and Turkmenistan, which both claimed no COVID-19 deaths at all, have no credible reported death rate; and Hong Kong has no reported death rate in the file. Complete-case analysis drops them silently, in two ways. Cuba and Hong Kong have observed excess death rates (Cuba's one of the highest in the Americas), and they vanish from the fit; and Somalia, Syria, North Korea and Turkmenistan, with about 72 million people between them, vanish from the world total, as if nobody died there.

The Bayesian remedy is to treat each missing predictor value as a parameter. It needs a prior, and the prior should itself be a model: here, log GDP is predicted from life expectancy and median age (which are known for every country), and log reported deaths likewise. Each missing value is then informed from two directions at once: by its imputation model, which says what values are plausible for a country with that life expectancy and age, and by the outcome model, which says what value would be consistent with the country's observed outcome. In brms, a variable with missing values is marked | mi() on the left of its own formula and enters other formulas as mi(variable). The outcome is marked | mi() too, so the missing excess death rates become parameters as well, and the whole imputation happens in one fit.

f_mi <- bf(y | mi() ~ mi(r) + a + mi(g)) +       # outcome model: missing y imputed too
        bf(r | mi() ~ le + a) +                    # imputation model for log reported deaths
        bf(g | mi() ~ le + a) +                    # imputation model for log GDP
        set_rescor(FALSE)
m_mi <- brm(f_mi, data = covid,
            prior = c(prior(normal(2, 2), class = Intercept, resp = y),
                      prior(normal(0, 2), class = b, resp = y),
                      prior(exponential(1), class = sigma, resp = y),
                      prior(normal(0, 1), class = b, resp = r),
                      prior(normal(0, 1), class = b, resp = g),
                      prior(exponential(1), class = sigma, resp = r),
                      prior(exponential(1), class = sigma, resp = g)),
            seed = 7, refresh = 0, file = "fits/ch7_mi")
fixef(m_mi)[c("y_Intercept", "y_mir", "y_a", "y_mig"), ] |> round(2)

            Estimate Est.Error  Q2.5 Q97.5
y_Intercept     2.64      0.22  2.19  3.05
y_mir           0.76      0.12  0.53  1.00
y_a             0.81      0.28  0.23  1.36
y_mig          -1.46      0.26 -1.98 -0.94

The outcome coefficients barely move from the complete-case fit (as they should, with only two more observed outcomes, Cuba's and Hong Kong's, in the fit), but now every country is in the model. The posterior draws contain parameters named Ymi_g[...] for the imputed GDPs and Ymi_y[...] for the imputed excess death rates. Here are the imputed GDPs.

ix <- which(is.na(covid$g))
ymi_g <- as_draws_df(m_mi) |> select(starts_with("Ymi_g")) |> as.matrix()
tibble(country = covid$location[ix], life_expectancy = covid$life_expectancy[ix],
       excess_observed = !is.na(covid$y[ix]),
       gdp_median = 1e4 * exp(apply(ymi_g, 2, median)),
       gdp_lo = 1e4 * exp(apply(ymi_g, 2, quantile, .05)),
       gdp_hi = 1e4 * exp(apply(ymi_g, 2, quantile, .95))) |>
  mutate(across(starts_with("gdp"), ~ round(.x, -2)))

# A tibble: 4 x 6
  country     life_expectancy excess_observed gdp_median gdp_lo gdp_hi
  <chr>                 <dbl> <lgl>                <dbl>  <dbl>  <dbl>
1 Cuba                   78.8 TRUE                 19400   7600  49500
2 North Korea            72.3 FALSE                12900   4400  36300
3 Somalia                57.4 FALSE                 1600    500   4600
4 Syria                  72.7 FALSE                 6400   2300  18600

Cuba's imputed GDP is about $19,000 per person (90% interval roughly $8,000 to $50,000), Syria's about $6,000 and Somalia's about $1,600. North Korea's is about $13,000, because its life expectancy (72 years) and median age resemble those of middle-income countries. Most outside estimates put North Korea's income per person at a small fraction of that, and Syria's, after a decade of war, is probably far below its imputed value too. These are not failures of the method. The imputation model predicts GDP from life expectancy and age, and North Korea and wartime Syria are exactly the kind of countries for which that relationship breaks: an MNAR problem, in which the reason the value is missing (an isolated or collapsed state) is related to the value itself. Here the consequence is small, because neither country has observed excess deaths and their imputed GDPs affect only their own imputed excess rates; but it is a reminder to look at the imputed values, not just at the coefficients.

For Cuba, the two sources disagree in an instructive way. Its life expectancy (79) suggests a fairly rich country; its very high excess mortality, given modest reported deaths, is what the outcome model expects of a poorer one. The posterior lands in between, and is narrower than either source alone. For Somalia, Syria and North Korea, there is no observed outcome, so the outcome model contributes nothing (try them in the widget) and the imputation is simply the prediction of the imputation model. Lowering $\sigma$ in the widget shows what would happen if the outcome model were more precise: it would dominate.

iy <- which(is.na(covid$y))
ymi_y <- as_draws_df(m_mi) |> select(starts_with("Ymi_y")) |> as.matrix()
world_mi <- ymi_y %*% covid$population[iy] / 1e9 + with(covid[-iy, ], sum(y * population)) / 1e9
c(countries_imputed = length(iy), quantile(world_mi, c(.5, .025, .975)))

countries_imputed               50%              2.5%             97.5% 
           72.000            15.013             7.823            22.244 

With all 158 countries in the model, the world total is about 15 million (roughly 8 to 22 million): a little higher than before, because Somalia, Syria, North Korea and Turkmenistan now contribute, and Cuba's large excess helps estimate the GDP coefficient.

7.6 Shortcuts that fail, and multiple imputation

Before Bayesian imputation was easy, analysts filled holes in simpler ways, and many still do. Two are especially common.

• Mean imputation replaces each missing value by the mean of the observed values. It invents data that lie on a single vertical line, which distorts every relationship involving the variable.
• Single regression imputation replaces each missing value by its prediction from a regression on other variables. Better, but every imputed value lies exactly on the regression line, so the imputed data are less noisy than real data; relationships are exaggerated, and the analysis treats invented values as if they had been measured, so its standard errors are too small.

The fix for the second problem is to impute with noise, from a distribution, and to do it several times. Each completed data set gives an estimate and a standard error; the variation between the completed data sets is the extra uncertainty due to the missing values. This is multiple imputation, and its rules for combining the results are known as Rubin's rules: with $M$ imputations giving estimates ${\widehat{q}}_m$ with variances $U_m$, report the average $\overline{q}=\frac{1}{M}{\sum }_m{\widehat{q}}_m$ with variance

The widget compares the four approaches in a setting where the truth is known: life expectancy regressed on log GDP, with some GDPs deleted.

Repeat the deletion 200 times and the pattern is stark. Complete cases give a slope that is too small, because deleting countries by their life expectancy (the outcome here) distorts the relationship. Mean imputation gives a similarly biased slope with a standard error twice as large. Single regression imputation overshoots the slope and reports a standard error smaller than the full data would give: a precision it does not have. Multiple imputation recovers the full-data slope on average, with a standard error only slightly larger than the full data's, which is exactly right, since some information really was lost.

In R, the mice package (van Buuren and Groothuis-Oudshoorn, 2011) implements multiple imputation by chained equations: each incomplete variable is imputed from all the others in turn, round and round until the imputations stabilise. It is the standard tool in epidemiology and the social sciences, and brms can fit a model to each completed data set and pool the draws with brm_multiple(). Use the joint Bayesian model when the imputation and analysis models are natural to write together, as here; use mice when there are many incomplete variables of different types, or when the imputed data will be shared with other analysts.

7.7 Missing not at random: sensitivity analysis

Everything so far assumed MAR: that among countries with the same reported deaths, median age and income, those with excess mortality data and those without had the same excess mortality on average. Is that plausible? There are reasons to doubt it in both directions. Countries whose health systems were overwhelmed may also have been the ones least able to register deaths, which would make the missing countries worse than the model predicts. On the other hand, much of the missing population lives in young, rural populations where transmission and mortality may have been lower than the observed relationships suggest.

The data cannot settle this, so a careful analysis does not pretend to. Instead, it asks how large a departure from MAR would change the conclusion. The simplest version is a delta adjustment: suppose that each missing country's excess death rate is $\delta$ per 1,000 higher than the model predicts, recompute the world total for a range of values of $\delta$, and report how the conclusion depends on it. A value of $\delta =1$ per 1,000 is a substantial departure: about the difference between the predicted excess rates of a country with median age 30 and one with median age 42.

pop_missing <- sum(covid$population[iy]) / 1e9               # billions of people with no data
tibble(delta = c(-1, -0.5, 0, 0.5, 1, 2)) |>
  mutate(world = map(delta, ~ as.numeric(world_mi) + .x * pop_missing),   # delta per 1,000 x billions = millions
         median = map_dbl(world, median),
         lo = map_dbl(world, quantile, .025), hi = map_dbl(world, quantile, .975)) |>
  select(-world) |> mutate(across(c(median, lo, hi), ~ round(.x, 1)))

# A tibble: 6 x 4
  delta median    lo    hi
  <dbl>  <dbl> <dbl> <dbl>
1  -1      9.8   2.6  17  
2  -0.5   12.4   5.2  19.6
3   0     15     7.8  22.2
4   0.5   17.6  10.4  24.9
5   1     20.3  13.1  27.5
6   2     25.5  18.3  32.7

Under MAR ($\delta =0$), the median world total is 15.0 million. If the missing countries had one more excess death per 1,000 people than predicted, it would be 20.3 million; one fewer, 9.8 million. Each unit of $\delta$ moves the total by about 5.2 million, because about 5.2 billion people live in the 72 countries whose excess deaths are imputed.

Replacing the model's negative extrapolation for China with zero raises the median to about 16.4 million and, because China's huge interval is removed, narrows the 95% interval to about 11 to 22 million. Two conclusions are robust to every plausible value of $\delta$: excess deaths were several times the reported COVID-19 deaths, and the naive extrapolation was far too high. The central number is not robust: it moves by about 5 million for every death per 1,000 by which the missing countries depart from the model. Reporting "about 15 million under MAR, and 12 to 18 million if the missing countries differ from the model's predictions by half a death per 1,000 in either direction" is less satisfying than a single number, and much more honest.

7.8 Measurement error: the same idea

A missing value is an extreme case of a measurement error: a measurement so noisy that it carries no information. Most measurements lie between perfect and absent, and the Bayesian treatment is the same throughout. Treat the true value as a parameter, give it a model, and let the noisy measurement inform it through a likelihood that says how noisy it is. Excess deaths are themselves estimates, with standard errors from the projection of expected deaths; reported COVID-19 deaths are a badly biased measurement of COVID-19 deaths; survey responses about income are noisy and sometimes rounded. In brms, an outcome measured with known standard error is written y | se(sd_y, sigma = TRUE) (the meta-analysis syntax of Chapter 6, with the residual variation kept), and a predictor measured with error gets its own formula, bf(x | mi(sd_x) ~ 1), and enters the outcome model as mi(x), which treats the observed value as a noisy measurement of the true one with standard deviation sd_x. Ignoring measurement error in a predictor biases its coefficient toward zero (the classic "regression dilution"), and the measurement-error model corrects that bias by propagating the uncertainty, just as imputation does.

7.9 Missing data in practice

7.10 Exercises

Each exercise has a hint and an answer button.

Further reading and references

1. COVID-19 Excess Mortality Collaborators (2022). Estimating excess mortality due to the COVID-19 pandemic: a systematic analysis of COVID-19-related mortality, 2020–21. The Lancet, 399(10334), 1513–1536.
2. Jha, P. et al. (2022). COVID mortality in India: national survey data and health facility deaths. Science, 375(6581), 667–671.
3. Karlinsky, A. and Kobak, D. (2021). Tracking excess mortality across countries during the COVID-19 pandemic with the World Mortality Dataset. eLife, 10, e69336. Data via Our World in Data: github.com/owid/covid-19-data.
4. Little, R. J. A. and Rubin, D. B. (2019). Statistical Analysis with Missing Data (3rd ed.). Wiley.
5. McElreath, R. (2020). Statistical Rethinking (2nd ed.). CRC Press.
6. Mohan, K. and Pearl, J. (2021). Graphical models for processing missing data. Journal of the American Statistical Association, 116(534), 1023–1037.
7. Msemburi, W., Karlinsky, A., Knutson, V., Aleshin-Guendel, S., Chatterji, S. and Wakefield, J. (2023). The WHO estimates of excess mortality associated with the COVID-19 pandemic. Nature, 613, 130–137.
8. National Research Council (2010). The Prevention and Treatment of Missing Data in Clinical Trials. National Academies Press.
9. Rubin, D. B. (1976). Inference and missing data. Biometrika, 63(3), 581–592.
10. Rubin, D. B. (1987). Multiple Imputation for Nonresponse in Surveys. Wiley.
11. Serfling, R. E. (1963). Methods for current statistical analysis of excess pneumonia-influenza deaths. Public Health Reports, 78(6), 494–506.
12. Sterne, J. A. C. et al. (2009). Multiple imputation for missing data in epidemiological and clinical research: potential and pitfalls. BMJ, 338, b2393.
13. van Buuren, S. (2018). Flexible Imputation of Missing Data (2nd ed.). CRC Press.
14. van Buuren, S. and Groothuis-Oudshoorn, K. (2011). mice: multivariate imputation by chained equations in R. Journal of Statistical Software, 45(3), 1–67.
15. von Hippel, P. T. (2007). Regression with missing Ys: an improved strategy for analyzing multiply imputed data. Sociological Methodology, 37(1), 83–117.