Missing Data Matters (Part I): Mechanisms

Exploring Missing Data

There are lots of good resources for learning how to effectively diagnose the missing data mechanism, so in the interest of not turning this blog post into a tome, I won’t go into incredible detail about this process. Tierney and Horst (2022)’s The Missing Book is fantastic for learning about the missing data exploratory analysis workflow using the R package naniar⁴. The mice and ggmice R packages offer additional functionality for investigating missingess⁵. Instead of detailing an entire exploratory workflow, I will give a brief overview, with example plots⁶. While not exhaustive, it will hopefully serve as inspiration to get you started.

I’ve transformed the attrition dataset to include missing data spread across multiple variables. Most missing data is MCAR, but the years an employee has spent at the company is MAR. Missing values for employee tenure are associated with low job satisfaction and high monthly income (greater than $15k per month).

set.seed(123)

missing_years <-
  attrition |>
  mutate(
    across(starts_with("job") & !job_satisfaction, 
           ~ replace(.x, runif(n()) < 0.01, NA)),
    years_at_company =
      replace(
        years_at_company,
        runif(n()) < 0.75 & 
          (job_satisfaction == "Low" | monthly_income > 15000), NA
        ),
    relationship_satisfaction = 
      replace(relationship_satisfaction, runif(n()) < 0.075, NA),
    work_life_balance = replace(work_life_balance, runif(n()) < 0.05, NA)
  )

The starting point in any exploratory analysis of missing data is to get some idea of the quantity of missingness and where it is. The plot below⁷ summarises the missingness in the attrition dataset.

missing_years |> 
  summarise(across(everything(), ~ sum(is.na(.x))/n()) ) |> 
  tidyr::pivot_longer(
    cols = everything(), names_to = "variable", values_to = "value"
    ) |> 
  filter(value != 0) |>
  mutate(
    variable = 
      stringr::str_to_title(
        stringr::str_replace_all(variable, "_", " ")
        ),
    variable = forcats::fct_reorder(variable, -value, .desc = TRUE)
    ) |> 
  ggplot(aes(x = variable,  y = value)) +
  geom_col(aes(fill = as.numeric(variable) == 6), colour = "grey20") +
  geom_text(
    aes(
      label = paste0(sprintf("%2.1f", value * 100), "%"),
      colour = value > .03, hjust = if_else(value > .03, 1.2, -.2)
      ),
    size = 8, fontface = "bold", family = "Poppins"
  ) +
  geom_hline(yintercept = 0, colour = "grey20", linewidth = 0.8) +
  coord_flip() +
  scale_x_discrete(
    name = NULL, labels = scales::label_wrap(15)
    ) +
  scale_y_continuous(guide ="none", name = NULL, expand = c(1e-03, 1e-03)) +
  scale_fill_manual(values = c("grey50", "#026E99"), guide = "none") +
  scale_colour_manual(values = c("grey20", "white"), guide = "none") +
  labs(
    title = "Proportion of Values Missing for All Variables",
    subtitle = 
      stringr::str_wrap(
        glue::glue(
          "The percentage of missing values for variables that contain ",
          "missingness. Relationship Satisfaction and Work Life Balance contain ",
          "significant missing values, but Years at Company has ~3x as many as ",
          "any other variable."
        ),
        width = 93
      ),
    caption = 
      "Visualisation: Paul Johnson  |  Data: IBM Watson Analytics Lab"
  ) +
  theme(
    plot.subtitle = element_text(margin = margin(5, 0, 10, 0)),
    axis.text.y = element_text(
      colour = "grey30", size = rel(1.2), hjust = 0.5, 
      margin = margin(0, 0.5, 0, 0.5), lineheight = .4
      ),
    axis.ticks = element_blank()
    )

Having already detailed what I did when adding missingness to the data, I feel a little silly not just visualising this but now describing the results. It turns out that there are a bunch of missing values spread across six variables, but the biggest culprit is an employee’s tenure at the company. Astonishing. Whoever could have predicted this?

A rigorous analysis would look at each variable and diagnose the mechanism driving the missingness. Given I’m the underlying cause for all missing values, we will take the shortcut and focus on employee tenure.

Assume MAR; Test MNAR/MCAR

Generally speaking, MCAR data is rare. Missing values can occur for many reasons, and the likelihood that the cause is independent of all other variables is low.

It is more plausible that the data is either MAR or MNAR (Nahhas 2024). However, dealing with MNAR data is particularly complicated⁸. It is reasonable to assume your data is MAR until you find enough evidence to suggest otherwise (in either direction).

In practice, identifying the type of missingness you’re dealing with requires exploratory analysis, domain expertise, and good judgment. Testing the assumption that your data is not MNAR involves exploring missing values and looking for evidence that missingness may not be random. For example, observed values may be constrained at a value lower than your expertise suggests should be the maximum, indicating some censoring of the variable. Simple summary functions and plotting distributions can uncover evidence that missing values are non-random. However, concluding that data is MNAR requires domain expertise and a good understanding of your data.

It is a little easier to test the assumption data is not MCAR because any association between observed and missing values acts as evidence data is MAR. One way to look for associations is to create a binary variable that represents a variable’s missingness using naniar::bind_shadow() and plotting the distribution of other variables conditional on missingness. Spikes in the distribution of missing values, relative to the distribution of observed values, indicate a correlation between that variable and the missing values.

The plot below visualises the monthly income distribution, conditional on employee tenure’s missingness. A thorough analysis would require digging to find associations, but I’m cheating and already know where to look.

missing_years |> 
  naniar::bind_shadow() |> 
  ggplot(aes(x = monthly_income, fill = forcats::fct_rev(years_at_company_NA))) + 
  geom_histogram(position = "dodge", colour = "grey20") +
  geom_hline(yintercept = 0, colour = "grey20", linewidth = 0.8) +
  scale_fill_manual(
    labels = c("Missing", "Not Missing"), 
    values = c("#D93649", "#026E99")
    ) +
  scale_x_continuous(labels = scales::label_currency()) +
  labs(
    title = "Years at Company Missingness Distribution",
    subtitle = 
      stringr::str_wrap(
        glue::glue(
          "The distribution of Monthly Income conditional on the missingness of ",
          "Years at Company. While the distribution of not missing values has a ",
          "is much larger scale across the majority of the income range, ",
          "there are more missing values from $15-20k."
        ),
        width = 95
      ),
    x = "Monthly Income", y = NULL,
    caption = 
      "Visualisation: Paul Johnson  |  Data: IBM Watson Analytics Lab"
  )

We expect the missing data distribution to follow a similar pattern to the observed values if it is MCAR. That is more or less the case for the monthly income distribution, but there is a spike at the higher end from around $15-20k. That suggests we might be dealing with MAR data.

We can also use naniar::geom_miss_point() to plot the association between employee tenure and income, and we can facet the data into separate plots for each level of job satisfaction (since we already know that job satisfaction is associated with tenure missingness).

missing_years |>  
  ggplot(aes(x = years_at_company, y = monthly_income)) +
  naniar::geom_miss_point(size = 1, alpha = 0.8) + 
  facet_wrap(vars(job_satisfaction)) +
  scale_colour_manual(values = c("#D93649", "#026E99")) +
  scale_y_continuous(labels = scales::label_currency()) +
  guides(colour = guide_legend(override.aes = list(size = 3))) +
  labs(
    title = "Years at Company Missingness Associations",
    subtitle = 
      stringr::str_wrap(
        glue::glue(
          "The association between Years at Company and Monthly Income, ",
          "conditional on Years at Company's missingness and split by Job ",
          "Satisfaction. For low satisfaction, missing values are distributed ",
          "across all values of monthly income. For all other levels of ",
          "satisfaction, missing values are clustered between $15k and $20k."
        ),
        width = 92
      ),
    x = "Years at Company", y = "Monthly Income",
    caption = 
      "Visualisation: Paul Johnson  |  Data: IBM Watson Analytics Lab"
  )

While the income distribution indicated that employee tenure might contain missing values that are not MCAR, the above scatter plots provide strong evidence we’re dealing with MAR data. The scatter plots show that missing values are limited to the $15-20k monthly income range, except when job satisfaction is low, where missingness exists across all incomes. The fact that missing values group around specific ranges of observed values indicates that employee tenure’s missing values are not randomly distributed and suggests a relationship between these missing values and monthly income and job satisfaction. These results provide strong evidence that the missingness in the employee tenure variable is not MCAR.

It is also possible to compute statistical tests of the association between missing values and the observed values in the data. One option is to model missing values using logistic regression, treating missingness as the binary outcome (using naniar::bind_shadow() as before) and the rest of the dataset as explanatory variables. Table 1 shows the outputs from a regression on employee tenure’s missing values. I have rescaled monthly income (and monthly rate) using arm::rescale()⁹ because it is on a significantly larger scale than other variables, and this makes it difficult to interpret the effect size.

missing_years |> 
  naniar::bind_shadow() |> 
  select(years_at_company_NA | !years_at_company & !ends_with("NA")) |> 
  mutate(
    across(where(is.factor) & !years_at_company_NA, as.numeric),
    monthly_income = arm::rescale(monthly_income),
    monthly_rate = arm::rescale(monthly_rate)) |> 
  glm(years_at_company_NA ~ ., family = "binomial", data = _) |> 
  list("Years at Company Missingness" = _) |> 
  modelsummary::modelsummary(
    coef_rename = modelsummary::coef_rename, 
    coef_omit = c(-7, -16, -18, -19, -26),
    gof_omit =  "IC|Log", exponentiate = TRUE, 
    stars = TRUE, output = "gt"
  ) |>
  tbl_theme()

I have only included statistically significant variables¹⁰. While this is generally a bad practice, regression tables with 31 variables look like a mess in the middle of a blog post. A narrow focus on p-values is a bad idea in any regression model because effect size tells us a lot more than significance, but in this case, we’re taking some shortcuts. The extremely small p-values for monthly income and job satisfaction indicate that the probability of observing an effect size at least as extreme as our model coefficients if the data was MCAR is small¹¹. These results support our assumption that the data is MAR.

Another option is to use naniar’s mcar_test() to calculate Little (1988)’s test statistic for MCAR data¹². Table 2 shows the results of running Little’s MCAR test on our dataset.

missing_years |>  
  naniar::mcar_test() |> 
  mutate(across(everything(), ~ round(.x))) |> 
  rename(
    "Test Statistic" = statistic,
    "Degrees of Freedom" = df,
    "P-Value" = p.value,
    "Missing Patterns" = missing.patterns
  ) |> 
  gt::gt() |> 
  gt::cols_align(align = "center", columns = 1:4) |>
  tbl_theme(width = 80)

A high test statistic value and a p-value of zero means we reject the null hypothesis. If the data is MCAR, the probability of observing mean differences as or more extreme than this is very low, suggesting the missing data is either MNAR or MAR. However, there are several limitations to Little’s MCAR test. First, it tests all missing data across the entire dataset instead of testing specific missing patterns. Running logistic regressions on specific missingness is a lot more granular and can be more informative. Second, Little’s MCAR test relies on a slightly backwards approach to hypothesis testing. The null hypothesis is that the data is MCAR. Failing to reject the null is not the same as being able to conclude that the null is true, so it isn’t a test for MCAR data as much as it is a test to reject the hypothesis that the data is MCAR. I’d argue this is fine because the starting assumption is that data is MAR. Little’s MCAR test provides a broad summary of the data. A high p-value suggests further investigation is needed because data might not be MAR, and a low p-value is evidence to support the assumption that the data is MAR. However, I think logistic regressions are more informative, robust, and equally easy to compute, so I would recommend this approach instead.

Finally, another approach to identifying associations between missingness and observed values, which makes fewer assumptions of the data than either modelling approach, is calculating the mean difference of variables conditional on the missingness of the variable in question.

Ultimately, there are no shortcuts in this process. Whichever approach you take, it can only provide supporting evidence, not confirmation. A detailed exploratory analysis combined with domain expertise is the only way to understand how the data may have come to be missing.