Why Missing Data Matters

Missing data are the rule, not the exception in longitudinal studies.

Three critical implications:

The missing data mechanism is the critical issue we must address.

Motivating Example 1: Six Cities Study

Study Design:

Main Reason for Missing Data: Moving in/out of school district

Six Cities: Two Scenarios

Question: How do we determine which scenario applies and handle each appropriately?

Motivating Example 2: Muscatine Study

Study Design:

Missing Data Patterns:

• Less than 40% had complete data at all three occasions
• Main reasons: (1) Failure to obtain consent, (2) Absence on exam day

Muscatine Study: Two Missing Data Scenarios

Key Insight: Same study, different reasons for missing data leads to different analysis strategies needed.

Central Theme of This Lecture

The missing data mechanism determines:

• Which methods yield valid inferences
• Whether we need to model the missingness process
• How to conduct sensitivity analyses

Notation: Response Vectors and Indicators

Consider a study designed to collect \(n\) measurements per subject.

Complete response vector (intended): \[Y_i = (Y_{i1}, Y_{i2}, \ldots, Y_{in})'\]

Response indicator vector: \[R_i = (R_{i1}, R_{i2}, \ldots, R_{in})'\]

where \(R_{ij} = 1\) if \(Y_{ij}\) is observed, and \(R_{ij} = 0\) if \(Y_{ij}\) is missing.

Example: For a 6-visit study (visits 0-5), if a subject drops out after visit 2: \[R_i = (1, 1, 1, 0, 0, 0)'\]

Covariates: \(X_i\) (assumed fully observed)

Partitioning the Response Vector

Given the response indicators \(R_i\), partition \(Y_i\) into:

Key observation: \(R_i\) is recorded for all individuals and stratifies the population into distinct sub-populations defined by missing data patterns.

Table 17.1: Missing Data Patterns as Stratification

Each row represents a different missing data pattern. \(R_i\) divides the target population into subpopulations.

The Three Missing Data Mechanisms

We classify missing data by considering how \(R_i\) is related to \(Y_i\):

Terminology Warning: The nomenclature is NOT intuitive!

• “Missing at Random” does NOT mean “randomly missing”
• These are technical definitions with specific mathematical meanings

MCAR: Missing Completely at Random

Formal Definition:

\[\Pr(R_i \mid Y_i^O, Y_i^M, X_i) = \Pr(R_i \mid X_i)\]

Intuition: Missingness is unrelated to both observed AND unobserved responses (given covariates).

Bivariate Example: \(Y_i = (Y_{i1}, Y_{i2})'\), \(Y_{i1}\) fully observed, \(Y_{i2}\) sometimes missing.

If \(Y_{i2}\) is MCAR: \[\Pr(R_{i2} = 1 \mid Y_{i1}, Y_{i2}, X_i) = \Pr(R_{i2} = 1 \mid X_i)\]

Probability that \(Y_{i2}\) is missing does NOT depend on \(Y_{i1}\) or \(Y_{i2}\).

MCAR: Examples

Example 1: Rotating Panel Design

• Survey design where individuals rotate in/out by design
• Timing of measurements determined a priori
• Missingness mechanism under investigator control

Example 2: Six Cities Study (Scenario A)

• Child moves due to parent’s job relocation
• Completely unrelated to child’s pulmonary function

MCAR: General vs Strict

Subtle but important distinction:

Key: MCAR requires conditional independence given all relevant covariates in \(X_i\).

MCAR: Covariate-Dependent Examples

Example 1: Clinical trial dropout by treatment arm

• Treatment group: 20% dropout, Control: 15% dropout
• Dropout unrelated to outcomes, just due to pill burden
• If treatment arm is included as covariate: MCAR given treatment

Example 2: Side effects

• Missingness related to treatment side effects
• If “side effects” is measured and in \(X_i\): data are MCAR (given \(X_i\))
• If “side effects” NOT measured: data are NOT MCAR

MCAR: Implications for Distribution

Essential Feature: Observed data are a random sample of complete data.

Consequences:

MCAR: Implications for Analysis

Under MCAR, almost any method yields valid inferences:

Caveat: Even though methods are unbiased under MCAR, there is still loss of precision (efficiency).

MAR: Missing at Random

Formal Definition:

\[\Pr(R_i \mid Y_i^O, Y_i^M, X_i) = \Pr(R_i \mid Y_i^O, X_i)\]

Intuition: Missingness depends on observed responses and covariates, but NOT on unobserved responses (given \(Y_i^O\) and \(X_i\)).

Bivariate Example: If \(Y_{i2}\) is MAR: \[\Pr(R_{i2} = 1 \mid Y_{i1}, Y_{i2}, X_i) = \Pr(R_{i2} = 1 \mid Y_{i1}, X_i)\]

Probability that \(Y_{i2}\) is missing depends on observed \(Y_{i1}\), but NOT on unobserved \(Y_{i2}\).

MAR: Examples

Example 1: Protocol-Driven Removal

• Study protocol removes subject when outcome falls outside clinical range
• Missingness in \(Y_i\) depends only on observed components of \(Y_i\)

Example 2: Six Cities Study (Scenario B - MAR version)

• Family relocates due to child’s respiratory problems
• Decision based on recorded history of pulmonary function measurements
• Missingness predictable from \(Y_i^O\) only

Non-Example (MNAR): If relocation based on unobserved variable (e.g., parent’s assessment of future health)

MAR: Stratification on Observed Values

If subjects are stratified on values of \(Y_i^O\), then within strata, missingness in \(Y_i^M\) is like a chance mechanism.

Visual analogy:

• Imagine sorting subjects by their observed trajectory \(Y_i^O\)
• Within each “trajectory group,” those missing \(Y_i^M\) are random subset
• Across groups, missingness probabilities can differ

Mathematical consequence: \[\Pr(R_i \mid Y_i^O, Y_i^M, X_i) = \Pr(R_i \mid Y_i^O, X_i)\]

MAR: Implications for Distribution

Key Difference from MCAR:

Because missingness now depends on \(Y_i^O\):

Critical Implication: Valid inferences require correctly modeling the joint distribution of \(Y_i\) (both mean AND covariance).

MAR: Predicting Missing Values

Under MAR, missing values can be validly predicted from observed data.

For multivariate normal responses:

\[E(Y_i^M \mid Y_i^O) = \mu_i^M + \Sigma_i^{MO} (\Sigma_i^O)^{-1} (Y_i^O - \mu_i^O)\]

This is the conditional expectation formula from multivariate normal theory (recall properties of partitioned multivariate normal from earlier lectures on covariance structures).

MAR: Why Covariance Matters

The prediction formula requires:

Mean vector: \[\mu_i = \begin{pmatrix} \mu_i^O \\ \mu_i^M \end{pmatrix}\]

Covariance matrix: \[\Sigma_i = \begin{pmatrix} \Sigma_i^O & \Sigma_i^{OM} \\ \Sigma_i^{MO} & \Sigma_i^M \end{pmatrix}\]

Critical insight: The off-diagonal block \(\Sigma_i^{MO}\) determines how observed values inform missing values!

Therefore: Under MAR, correct covariance specification is essential for valid inference.

MAR: Implications for Analysis

Under MAR, method validity depends on assumptions:

MAR: The “Ignorable” Terminology

MAR (and MCAR) are often called “ignorable” mechanisms.

What “ignorable” means:

• We can ignore the model for \(\Pr(R_i \mid Y_i, X_i)\) in likelihood-based inference
• We only need to model \(f(Y_i \mid X_i)\)

What “ignorable” does NOT mean:

• “We can ignore the missing data problem”
• “Any method will work”
• “Complete-case analysis is fine”

Bottom line: We still must carefully model \(f(Y_i \mid X_i)\), including covariance!

Subtle Distinction: MCAR vs MAR

Confusion is common! Let’s clarify with an example:

Scenario: In longitudinal trial, subjects with worsening symptoms more likely to drop out

MNAR: Not Missing at Random

Formal Definition:

\[\Pr(R_i \mid Y_i^O, Y_i^M, X_i)\]

depends on \(Y_i^M\) (the unobserved values).

Equivalently: MNAR means the MAR assumption is violated: \[\Pr(R_i \mid Y_i^O, Y_i^M, X_i) \neq \Pr(R_i \mid Y_i^O, X_i)\]

Intuition: Missingness related to what the values would have been had they been observed.

MNAR: Examples

Example 1: Quality-of-Life Instruments

• Patients skip QoL questionnaire when feeling poorly
• Missingness related to unobserved (low) QoL value

Example 2: Muscatine Study (Scenario B)

• Parents of obese children less likely to consent
• Missingness related to unobserved obesity status

Example 3: Six Cities Study (MNAR version)

• Family relocates based on parent’s concern about future deterioration
• Not based on observed measurements, but on unmeasured prognosis

MNAR: Implications for Distribution

Under MNAR:

Critical consequence: The specific model chosen for \(\Pr(R_i \mid Y_i, X_i)\) drives the results.

MNAR: The Identification Problem

Fundamental issue: MNAR assumptions are unverifiable from observed data.

Without external information:

• Cannot distinguish one MNAR model from another
• Observed data provide no evidence for/against MNAR mechanisms
• Different plausible MNAR models can yield very different conclusions

Implication: Under MNAR, sensitivity analysis is essential

MNAR: Implications for Analysis

Under MNAR, standard methods are all biased:

Required approach: Joint modeling of outcome AND missingness

• Selection models: \(f(Y_i, R_i \mid X_i) = f(Y_i \mid X_i) \times \Pr(R_i \mid Y_i, X_i)\)
• Pattern-mixture models: \(f(Y_i, R_i \mid X_i) = f(Y_i \mid R_i, X_i) \times \Pr(R_i \mid X_i)\)

Quick Comparison: MCAR vs MAR vs MNAR

* = conditional on correct model specification

What is Dropout?

Definition: Dropout refers to the special case where:

\[\text{If } R_{ik} = 0, \text{ then } R_{ik+1} = \cdots = R_{in} = 0\]

Characteristics:

Dropout indicator: \(D_i = k\) means subject drops out between occasions \(k-1\) and \(k\)

Figure 17.1: Monotone Pattern Visualization

library(ggplot2)
library(dplyr)

# Create monotone dropout pattern
set.seed(667)
n_subjects <- 20
n_times <- 5

dropout_data <- tibble(
  id = rep(1:n_subjects, each = n_times),
  time = rep(1:n_times, n_subjects),
  dropout_time = rep(sample(2:(n_times+1), n_subjects, replace = TRUE), each = n_times)
) %>%
  mutate(
    observed = time < dropout_time,
    dropout_group = factor(dropout_time,
                          levels = 2:(n_times+1),
                          labels = paste("Dropout after Y", 1:n_times, sep=""))
  )

ggplot(dropout_data, aes(x = time, y = reorder(id, dropout_time), fill = observed)) +
  geom_tile(aes(alpha = observed), color = "white", linewidth = 0.5) +
  geom_tile(data = dropout_data %>% filter(observed),
            aes(fill = dropout_group), color = "white", linewidth = 0.5, alpha = 1) +
  scale_fill_manual(values = c("FALSE" = "gray85",
                               "Dropout after Y1" = "#08519c",
                               "Dropout after Y2" = "#3182bd",
                               "Dropout after Y3" = "#6baed6",
                               "Dropout after Y4" = "#9ecae1",
                               "Dropout after Y5" = "#c6dbef"),
                    name = "Status") +
  scale_alpha_manual(values = c("FALSE" = 1, "TRUE" = 1), guide = "none") +
  scale_x_continuous(breaks = 1:n_times, labels = paste0("Y", 1:n_times)) +
  scale_y_discrete(breaks = seq(5, 20, 5)) +
  labs(x = "Measurement Occasion",
       y = "Individual (ordered by dropout time)",
       title = "Monotone Missing Data Pattern",
       subtitle = "Color gradient shows dropout timing (darker = earlier dropout)") +
  theme_minimal(base_size = 14) +
  theme(legend.position = "bottom",
        panel.grid.minor = element_blank())

Notice: More observations at \(Y_j\) than at \(Y_{j+1}\) for all \(j\).

Dropout Taxonomy

Applying the MCAR/MAR/MNAR hierarchy to dropout:

Dropout: The Key Question

Central issue: Do those who drop out differ from those who stay?

Simulation Study Setup

We’ll illustrate the three mechanisms via simulation.

Data Generating Process:

• \(Y_{it} = \beta_0 + \beta_1 t + \varepsilon_{it}\) (linear trend)
• \(\text{Cov}(Y_{is}, Y_{it}) = \rho^{|s-t|}\) (AR(1) with \(\rho = 0.7\))
• \(N = 1000\) subjects, \(n = 5\) time points

Dropout Model:

\[\text{logit}\{\Pr(D_i = k \mid D_i \geq k, Y_{i,k-1}, Y_{ik})\} = \theta_1 + \theta_2(Y_{i,k-1} - \mu_{k-1}) + \theta_3(Y_{ik} - \mu_k)\]

Simulation: Three Scenarios

Simulation: Data Generation Function

library(MASS)
library(dplyr)
library(ggplot2)

simulate_dropout <- function(N = 1000, n_time = 5,
                             beta0 = 5, beta1 = 0.25, rho = 0.7,
                             theta1 = -0.5, theta2 = 0, theta3 = 0) {
  # Generate AR(1) covariance
  times <- 1:n_time
  Sigma <- outer(times, times, function(s, t) rho^abs(s - t))

  # Generate complete data
  mu <- beta0 + beta1 * times
  Y_complete <- mvrnorm(N, mu = mu, Sigma = Sigma)

  # Generate dropout
  dropout_time <- rep(n_time + 1, N)  # Initialize as completers

  for (i in 1:N) {
    for (k in 2:n_time) {
      # Logistic dropout probability
      lin_pred <- theta1 +
                  theta2 * (Y_complete[i, k-1] - mu[k-1]) +
                  theta3 * (Y_complete[i, k] - mu[k])

      prob_dropout <- plogis(lin_pred)

      if (runif(1) < prob_dropout) {
        dropout_time[i] <- k
        break
      }
    }
  }

  # Create observed data (set to NA after dropout)
  Y_obs <- Y_complete
  for (i in 1:N) {
    if (dropout_time[i] <= n_time) {
      Y_obs[i, dropout_time[i]:n_time] <- NA
    }
  }

  # Return data frame
  data.frame(
    Y_complete = Y_complete,
    Y_obs = Y_obs,
    dropout_time = dropout_time,
    subject = 1:N
  )
}

Figure 17.2: Complete Data (Baseline)

set.seed(667)
beta0 <- 5
beta1 <- 0.25
n_time <- 5

# Simulate complete data (no dropout)
dat_complete <- simulate_dropout(N = 1000, theta2 = 0, theta3 = 0)

# Calculate means at each time
means_complete <- colMeans(dat_complete[, paste0("Y_complete.", 1:n_time)])

# Plot
data.frame(
  time = 1:n_time,
  mean_obs = means_complete,
  true_mean = beta0 + beta1 * (1:n_time)
) %>%
  ggplot(aes(x = time)) +
  geom_line(aes(y = true_mean), linewidth = 1, color = "black") +
  geom_point(aes(y = mean_obs), size = 3, color = "steelblue") +
  scale_x_continuous(breaks = 1:n_time) +
  labs(x = "Time", y = "Y",
       title = "Complete Data: Sample Means vs Population Line",
       subtitle = "No missing data - means coincide with population values") +
  theme_minimal(base_size = 14)

Sample means virtually coincide with population regression line.

Scenario A: MCAR Dropout

set.seed(667)
dat_mcar <- simulate_dropout(N = 1000,
                              theta1 = -0.5, theta2 = 0, theta3 = 0)

# Calculate observed means
means_mcar <- colMeans(dat_mcar[, paste0("Y_obs.", 1:n_time)], na.rm = TRUE)

# Dropout summary
dropout_summary_mcar <- table(dat_mcar$dropout_time)
cat("Dropout times:\n")

About 38% dropout at each occasion (constant hazard).

Figure 17.3(a): MCAR Dropout

data.frame(
  time = 1:n_time,
  mean_obs = means_mcar,
  true_mean = beta0 + beta1 * (1:n_time)
) %>%
  ggplot(aes(x = time)) +
  geom_line(aes(y = true_mean), linewidth = 1, color = "black") +
  geom_point(aes(y = mean_obs), size = 3, color = "steelblue") +
  scale_x_continuous(breaks = 1:n_time) +
  labs(x = "Time", y = "Y",
       title = "Scenario A: MCAR Dropout",
       subtitle = "Observed means unbiased despite heavy dropout") +
  theme_minimal(base_size = 14)

Key Observation: Even with 85% missing at time 5, sample means are unbiased!

Scenario B: MAR Dropout

set.seed(667)
dat_mar <- simulate_dropout(N = 1000,
                             theta1 = -0.5, theta2 = 1.0, theta3 = 0)

# Calculate observed means
means_mar <- colMeans(dat_mar[, paste0("Y_obs.", 1:n_time)], na.rm = TRUE)

# Dropout summary
cat("Proportion missing at each time:\n")

Dropout depends on previous response: those with high \(Y_{i,k-1}\) more likely to drop out.

Figure 17.3(b): MAR Dropout

data.frame(
  time = 1:n_time,
  mean_obs = means_mar,
  true_mean = beta0 + beta1 * (1:n_time)
) %>%
  ggplot(aes(x = time)) +
  geom_line(aes(y = true_mean), linewidth = 1, color = "black") +
  geom_point(aes(y = mean_obs), size = 3, color = "coral") +
  scale_x_continuous(breaks = 1:n_time) +
  labs(x = "Time", y = "Y",
       title = "Scenario B: MAR Dropout (depends on past Y)",
       subtitle = "Observed means BIASED downward - high responders drop out preferentially") +
  theme_minimal(base_size = 14)

Because those with large \(Y_{i,k-1}\) drop out more, observed means are biased low.

Scenario C: MNAR Dropout

set.seed(667)
dat_nmar <- simulate_dropout(N = 1000,
                              theta1 = -0.5, theta2 = 0, theta3 = 0.8)

# Calculate observed means
means_nmar <- colMeans(dat_nmar[, paste0("Y_obs.", 1:n_time)], na.rm = TRUE)

# Dropout summary
cat("Proportion missing at each time:\n")

Dropout depends on current (unobserved) response: those with high \(Y_{ik}\) more likely to drop out at time \(k\).

Figure 17.3(c): MNAR Dropout

data.frame(
  time = 1:n_time,
  mean_obs = means_nmar,
  true_mean = beta0 + beta1 * (1:n_time)
) %>%
  ggplot(aes(x = time)) +
  geom_line(aes(y = true_mean), linewidth = 1, color = "black") +
  geom_point(aes(y = mean_obs), size = 3, color = "firebrick") +
  scale_x_continuous(breaks = 1:n_time) +
  labs(x = "Time", y = "Y",
       title = "Scenario C: MNAR Dropout (depends on current unobserved Y)",
       subtitle = "Observed means STRONGLY BIASED - cannot predict from observed data alone") +
  theme_minimal(base_size = 14)

Because those with large unobserved \(Y_{ik}\) drop out, bias is substantial.

Combined Comparison: Figure 17.3 (a, b, c)

library(tidyr)

# Calculate sample sizes at each timepoint for annotations
n_obs_mcar <- colSums(!is.na(dat_mcar[, paste0("Y_obs.", 1:n_time)]))
n_obs_mar <- colSums(!is.na(dat_mar[, paste0("Y_obs.", 1:n_time)]))
n_obs_nmar <- colSums(!is.na(dat_nmar[, paste0("Y_obs.", 1:n_time)]))

combined_data <- data.frame(
  time = rep(1:n_time, 3),
  mean_obs = c(means_mcar, means_mar, means_nmar),
  mechanism = rep(c("(a) MCAR", "(b) MAR", "(c) MNAR"), each = n_time),
  true_mean = rep(beta0 + beta1 * (1:n_time), 3),
  n_obs = c(n_obs_mcar, n_obs_mar, n_obs_nmar)
)

ggplot(combined_data, aes(x = time)) +
  geom_line(aes(y = true_mean), linewidth = 0.8, color = "black") +
  geom_point(aes(y = mean_obs, color = mechanism), size = 3) +
  geom_line(aes(y = mean_obs, color = mechanism), linewidth = 0.6, linetype = "dashed") +
  geom_text(aes(label = paste0("n=", n_obs), y = mean_obs),
            vjust = -0.8, size = 2.5, color = "gray30") +
  scale_color_manual(values = c("(a) MCAR" = "steelblue",
                                 "(b) MAR" = "coral",
                                 "(c) MNAR" = "firebrick")) +
  scale_x_continuous(breaks = 1:n_time) +
  facet_wrap(~ mechanism, ncol = 3) +
  labs(x = "Time", y = "Y",
       title = "Comparison of Observed Means Under Different Dropout Mechanisms",
       subtitle = "Black line = true population mean; Points = observed sample means (n = sample size)") +
  theme_minimal(base_size = 13) +
  theme(legend.position = "none",
        strip.text = element_text(face = "bold"))

Table 17.2: ML vs GEE Parameter Estimates

Now let’s fit models and compare estimates.

library(nlme)

# Reshape MCAR data to long format
dat_mcar_long <- dat_mcar %>%
  dplyr::select(subject, starts_with("Y_obs")) %>%
  pivot_longer(cols = starts_with("Y_obs"),
               names_to = "time",
               values_to = "y") %>%
  mutate(time = as.numeric(gsub("Y_obs.", "", time))) %>%
  filter(!is.na(y))

# ML fit (with AR(1) covariance)
fit_ml_mcar <- gls(y ~ time, data = dat_mcar_long,
                   correlation = corAR1(form = ~ time | subject),
                   method = "REML")

# GEE fit (working independence)
fit_gee_mcar <- geepack::geeglm(y ~ time, data = dat_mcar_long,
                                id = subject,
                                corstr = "independence")

# Extract coefficients
ml_mcar <- summary(fit_ml_mcar)$tTable
gee_mcar <- summary(fit_gee_mcar)$coefficients

Table 17.2: Results under MCAR

Both methods provide unbiased estimates under MCAR.

Table 17.3: Results under MAR

Key result: ML remains unbiased (correct covariance), GEE substantially biased.

Table 17.4: Results under MNAR

Both methods biased under MNAR! Slope estimates far from truth (0.25).

Interpreting Model Output Under Missing Data

Table 17.5: Complete Summary

Key Takeaways from Simulation

Real-Data Case Study: Six Cities FEV1

The motivating Six Cities study is in our fev1 data: annual lung-function (log FEV1) on children, measured from about age 6 to 19. Missing data here are real, not imposed: children enter and leave the district at different ages.

# Real FLW longitudinal data (Six Cities pulmonary thread); download fallback
fev1_url  <- "https://content.sph.harvard.edu/fitzmaur/ala2e/fev1.txt"
fev1_file <- "../../data/fev1.txt"
if (!file.exists(fev1_file)) {
  fev1_file <- "data/fev1.txt"
  if (!file.exists(fev1_file)) {
    dir.create("data", showWarnings = FALSE)
    download.file(fev1_url, fev1_file)
  }
}
fev1 <- read.table(fev1_file, header = TRUE)

# Real imbalance: how many measurement occasions does each child have?
occ_per_child <- table(fev1$id)
cat("Children:", length(unique(fev1$id)), "\n")

FEV1: Complete-Case Throws Away Most Children

library(nlme)
library(dplyr)
fev1 <- fev1 %>% mutate(age_c = age - 10)   # center age at 10 years

# Available-data fit (REML): uses ALL children, including partial records
fit_avail <- gls(logfev1 ~ age_c, data = fev1,
                 correlation = corCAR1(form = ~ age | id), method = "REML")

# Complete-case: only children with all 12 occasions
cc_ids  <- as.integer(names(which(table(fev1$id) == 12)))
fev1_cc <- fev1 %>% filter(id %in% cc_ids)
fit_cc  <- gls(logfev1 ~ age_c, data = fev1_cc,
               correlation = corCAR1(form = ~ age | id), method = "REML")

avail <- summary(fit_avail)$tTable["age_c", ]
cc    <- summary(fit_cc)$tTable["age_c", ]

Available-data (REML, all 300 children): slope 0.084 (SE 0.001).

Complete-case (only 25 children): slope 0.079 (SE 0.004).

Overview of Approaches

Each method makes different assumptions and has different validity conditions.

1. Complete-Case Analysis

Approach: Exclude all subjects with any missing data

When valid: MCAR only

Recommendation: Rarely acceptable in practice

Complete-Case Example

# Using MAR data from earlier
complete_cases <- dat_mar_long %>%
  group_by(subject) %>%
  filter(n() == n_time) %>%
  ungroup()

cat("Original sample:", length(unique(dat_mar_long$subject)), "subjects\n")

Under MAR, complete-case analysis is biased even though it uses “high quality” data!

2. Available-Data Analysis

Approach: Use all observed data (not just complete cases)

Examples: GLS, Standard GEE

When valid: MCAR only

Why it fails under MAR:

Available-data methods assume sample means/covariances of \(Y_i^O\) are unbiased for population parameters. This is TRUE under MCAR but FALSE under MAR!

Available-Data: Why Bias Occurs Under MAR

Intuition:

Imagine a trial where sicker patients drop out more often.

• Remaining patients at time \(t\) are healthier than full population
• Sample mean at time \(t\) is biased upward
• GEE uses these biased sample means
• Result: Biased inference about population

Mathematical reason:

GEE estimating equations use: \[\sum_{i=1}^N D_i^T V_i^{-1} (Y_i^O - \mu_i(X_i; \beta)) = 0\]

Under MAR, \(E(Y_i^O \mid X_i) \neq \mu_i(X_i; \beta)\) in general.

3. Likelihood-Based Methods (ML/REML)

Approach: Maximize likelihood based on observed data:

\[L(\beta, \Sigma) = \prod_{i=1}^N f(Y_i^O \mid X_i; \beta, \Sigma)\]

Examples: GLS with correctly specified covariance, LMMs (lme, lmer), GLMMs

When valid: MCAR or MAR (if model is correct!)

Critical requirement under MAR: Must correctly specify both mean model and covariance model.

Why Correct Covariance Matters Under MAR

Under MAR, prediction of missing values depends critically on the covariance structure.

Recall the conditional expectation formula:

\[E(Y_i^M \mid Y_i^O) = \mu_i^M + \Sigma_i^{MO} (\Sigma_i^O)^{-1} (Y_i^O - \mu_i^O)\]

Depends on covariance!

Misspecified covariance leads to wrong conditional mean and biased \(\hat{\beta}\).

Likelihood Example: Correct vs Misspecified Covariance

# Using MAR data
# Correct model: AR(1)
fit_correct <- gls(y ~ time, data = dat_mar_long,
                   correlation = corAR1(form = ~ time | subject),
                   method = "REML")

# Misspecified model: Independence
fit_wrong <- gls(y ~ time, data = dat_mar_long,
                 method = "REML")

cat("True slope: 0.25\n\n")

Correct covariance leads to unbiased; Misspecified leads to biased (under MAR)!

SAS Reference: ML under MAR (PROC MIXED)

The likelihood-based MAR analysis (our primary approach) in SAS. PROC MIXED with REML uses all available occasions, so it is valid under MAR if the mean and covariance are correct (no imputation needed).

/* Likelihood-based analysis valid under MAR (data in long form) */
proc mixed data = long method = reml;
  class id;
  model y = time / solution;
  repeated time / subject = id type = un;   /* unstructured covariance */
run;

For the multiple-imputation thread (a Chapter 18 preview), the SAS pipeline is PROC MI then PROC MIANALYZE:

proc mi data = wide nimpute = 20 seed = 667 out = imp;
  var y1 y2 y3 y4 y5;                        /* monotone or MCMC imputation */
run;
proc mixed data = imp; by _Imputation_;
  model y = time / solution covb;
  ods output SolutionF = mixparms;
run;
proc mianalyze parms = mixparms;            /* Rubin's rules pooling */
  modeleffects intercept time;
run;

4. Imputation Methods

Basic idea: Fill in (impute) missing values, then analyze complete data.

We’ll focus heavily on LVCF because it’s widely used (and widely misunderstood).

LVCF: The Conservative Treatment Effect Myth

Statistical folklore: “LVCF gives conservative estimate of treatment effect”

LVCF: When Might It Be Appropriate?

Rare scenario where LVCF is reasonable: Dropout due to cure or recovery

• Patient improves and stops treatment
• Response likely remains at good level
• Carrying forward “cured” state makes sense

In most clinical trials: Dropout often due to:

LVCF: Additional Problems

Beyond bias, LVCF has other serious flaws:

LVCF Bias: Setup

Simple two-timepoint design:

LVCF Bias: Pattern Mixture Model (Setup)

Stratify by dropout status (completer vs dropout):

Change from baseline for controls:

\[E(Y_{i2} - Y_{i1} \mid \text{trt}_i = 0, R_{i2} = 0) = \alpha_1 \quad \text{(dropouts)}\]

\[E(Y_{i2} - Y_{i1} \mid \text{trt}_i = 0, R_{i2} = 1) = \gamma_1 \quad \text{(completers)}\]

Interpretation:

• \(\alpha_1\): How control dropouts change from baseline
• \(\gamma_1\): How control completers change from baseline

LVCF Bias: Pattern Mixture Model (Treatment Group)

Change from baseline for treatment:

\[E(Y_{i2} - Y_{i1} \mid \text{trt}_i = 1, R_{i2} = 0) = \alpha_1 + \alpha_2 \quad \text{(dropouts)}\]

\[E(Y_{i2} - Y_{i1} \mid \text{trt}_i = 1, R_{i2} = 1) = \gamma_1 + \gamma_2 \quad \text{(completers)}\]

Parameter interpretation:

• \(\alpha_2\): Treatment effect for dropouts
• \(\gamma_2\): Treatment effect for completers

LVCF Bias: True Treatment Effect

Target parameter (what we want to estimate):

\[\delta = E(Y_{i2} - Y_{i1} \mid \text{trt}_i = 1) - E(Y_{i2} - Y_{i1} \mid \text{trt}_i = 0)\]

Average over dropout status:

For treatment group: \[E(Y_{i2} - Y_{i1} \mid \text{trt}_i = 1) = (\alpha_1 + \alpha_2)(1 - \pi_1) + (\gamma_1 + \gamma_2)\pi_1\]

For control group: \[E(Y_{i2} - Y_{i1} \mid \text{trt}_i = 0) = \alpha_1(1 - \pi_0) + \gamma_1\pi_0\]

Therefore: \[\delta = \alpha_2 + (\pi_1 - \pi_0)(\gamma_1 - \alpha_1) + \pi_1(\gamma_2 - \alpha_2)\]

LVCF Bias: LVCF Estimand

LVCF assumption: \(E(Y_{i2} \mid \text{trt}_i, R_{i2} = 0) = E(Y_{i1} \mid \text{trt}_i, R_{i2} = 0)\)

Equivalently: \(E(Y_{i2} - Y_{i1} \mid \text{trt}_i, R_{i2} = 0) = 0\)

This means LVCF assumes \(\alpha_1 = \alpha_2 = 0\).

LVCF treatment effect estimand:

\[\delta_{\text{LVCF}} = (\gamma_1 + \gamma_2)\pi_1 - \gamma_1\pi_0 = (\pi_1 - \pi_0)\gamma_1 + \pi_1\gamma_2\]

LVCF Bias: General Formula

Bias = LVCF estimate - True effect:

\[\text{Bias} = \delta_{\text{LVCF}} - \delta\]

After algebra:

\[\text{Bias} = (\pi_1 - \pi_0)\alpha_1 - (1 - \pi_1)\alpha_2\]

Key insight: Bias depends on:

• Retention rate difference: \((\pi_1 - \pi_0)\), scaling the control dropout change \(\alpha_1\)
• The treatment-by-dropout interaction \(\alpha_2\), entering through \(-(1 - \pi_1)\alpha_2\) (equivalently \(\alpha_2(\pi_1 - 1)\))

LVCF Bias: Under MCAR

Special case: MCAR dropout

Under MCAR, dropouts and completers have same expected change: \[\gamma_k = \alpha_k \quad \text{for } k = 1, 2\]

True treatment effect: \(\delta = \gamma_2 = \alpha_2\)

But LVCF gives: \[\delta_{\text{LVCF}} = (\pi_1 - \pi_0)\gamma_1 + \pi_1\gamma_2\]

Therefore: \[\text{Bias}_{\text{MCAR}} = (\pi_1 - \pi_0)\gamma_1 + (\pi_1 - 1)\gamma_2\]

LVCF Bias: The Counterintuitive Result

Why does this happen?

• LVCF assumes \(\alpha_1 = 0\) (no change after dropout)
• But natural time trend \(\gamma_1\) still exists
• If retention rates differ (\(\pi_1 \neq \pi_0\)), LVCF misattributes the time trend to treatment effect

LVCF Bias: Concrete Example

Scenario: No treatment effect, but response changes over time

True treatment effect: \(\delta = 0\)

LVCF estimand: \[\delta_{\text{LVCF}} = (0.8 - 0.6) \times 2.0 + 0.8 \times 0 = 0.4\]

Bias: 0.4 units (positive bias!)

LVCF creates spurious treatment effect where none exists!

LVCF Bias: Example in Opposite Direction

Scenario: Same as before but swap dropout rates

True treatment effect: \(\delta = 0\)

LVCF estimand: \[\delta_{\text{LVCF}} = (0.6 - 0.8) \times 2.0 + 0.6 \times 0 = -0.4\]

Bias: -0.4 units (negative bias!)

LVCF: Simulation Demonstration

set.seed(667)
N <- 1000
gamma1 <- 2.0  # Both groups improve
gamma2 <- 0.0  # No treatment effect

# Scenario 1: Higher retention in treatment
pi0_scenario1 <- 0.6
pi1_scenario1 <- 0.8

# Scenario 2: Higher retention in control
pi0_scenario2 <- 0.8
pi1_scenario2 <- 0.6

simulate_lvcf_bias <- function(N, gamma1, gamma2, pi0, pi1) {
  # Generate data
  trt <- rbinom(N, 1, 0.5)
  y1 <- rnorm(N, 10, 2)  # Baseline

  # True change
  true_change <- ifelse(trt == 0, gamma1, gamma1 + gamma2) + rnorm(N, 0, 1)
  y2_true <- y1 + true_change

  # Dropout
  r2 <- rbinom(N, 1, ifelse(trt == 0, pi0, pi1))

  # LVCF: If dropout, use y1 for y2
  y2_lvcf <- ifelse(r2 == 1, y2_true, y1)

  # Estimate treatment effect
  delta_true <- mean(y2_true[trt == 1] - y1[trt == 1]) - mean(y2_true[trt == 0] - y1[trt == 0])
  delta_lvcf <- mean(y2_lvcf[trt == 1] - y1[trt == 1]) - mean(y2_lvcf[trt == 0] - y1[trt == 0])

  c(true = delta_true, lvcf = delta_lvcf, bias = delta_lvcf - delta_true)
}

result1 <- simulate_lvcf_bias(N, gamma1, gamma2, pi0_scenario1, pi1_scenario1)
result2 <- simulate_lvcf_bias(N, gamma1, gamma2, pi0_scenario2, pi1_scenario2)

cat("Scenario 1: pi0=0.6, pi1=0.8 (higher treatment retention)\n")

Bias direction depends on dropout rates! LVCF is NOT conservative.

LVCF: Summary of Problems

Three fatal flaws:

Pattern-Mixture Models: An Alternative Framework

Key idea: Model outcomes differently for different dropout patterns

Pattern-mixture factorization:

\[f(Y_i, R_i \mid X_i) = f(Y_i \mid R_i, X_i) \times \Pr(R_i \mid X_i)\]

Why useful?

• Directly models heterogeneity across dropout groups
• Can incorporate sensitivity parameters
• Natural framework for tipping point analysis

Challenge: Unobserved outcomes in dropout groups require assumptions or sensitivity analysis.

Multiple Imputation: The 3-Step Procedure

Step 1: Create \(m\) complete datasets by imputing missing values \(m\) times

• Draw from \(f(Y_i^M \mid Y_i^O, X_i)\) using proper model
• Typically \(m = 10-25\) imputations

Step 2: Analyze each of the \(m\) datasets separately

• Get \(\hat{\beta}_1, \ldots, \hat{\beta}_m\) and \(\text{SE}_1, \ldots, \text{SE}_m\)

Step 3: Pool results using Rubin’s rules

• Combine estimates accounting for imputation uncertainty

Multiple Imputation: Rubin’s Rules

Combined estimate: \[\hat{\beta}_{\text{MI}} = \frac{1}{m}\sum_{j=1}^m \hat{\beta}_j\]

Combined variance: \[\text{Var}(\hat{\beta}_{\text{MI}}) = \bar{W} + \frac{m+1}{m}B\]

Where:

Multiple Imputation: Quick Example

library(mice)

# Use MAR data from earlier (small subset for demo)
dat_mi <- dat_mar %>%
  dplyr::select(subject, Y_obs.1, Y_obs.2, Y_obs.3, Y_obs.4, Y_obs.5) %>%
  rename(y1 = Y_obs.1, y2 = Y_obs.2, y3 = Y_obs.3, y4 = Y_obs.4, y5 = Y_obs.5)

# Perform MI (m=10 imputations)
# method = "pmm": Predictive mean matching (preserves data distribution)
imp <- mice(dat_mi[, -1], m = 10, method = "pmm", printFlag = FALSE, seed = 667)

# Fit model to each imputed dataset
fit_mi <- with(imp, lm(y5 ~ y1))

# Pool results using Rubin's rules
pooled <- pool(fit_mi)
summary(pooled, conf.int = TRUE)[2, ]

MI properly accounts for imputation uncertainty. (Details in Chapter 18!)

5. Inverse Probability Weighting (IPW)

Key idea: Weight observed data to correct for under-representation

Weight for subject \(i\):

\[w_i = \frac{1}{\Pr(\text{subject } i \text{ remains in study})}\]

Subjects with low probability of remaining get high weights (represent many similar dropouts).

Calculation for complete-case:

\[w_i = \frac{1}{\pi_{i1} \times \pi_{i2} \times \cdots \times \pi_{in}}\]

where \(\pi_{ik} = \Pr(D_i > k \mid D_i \geq k, Y_{i1}, \ldots, Y_{i,k-1}, X_i)\)

IPW: Estimation Procedure

Step 1: Estimate dropout probabilities

For each time \(k = 2, \ldots, n\):

• Use subjects still at risk at time \(k-1\)
• Fit logistic regression for \(\Pr(D_i > k \mid D_i \geq k)\)
• Obtain predicted probabilities \(\hat{\pi}_{ik}\)

Step 2: Calculate weights

\[\hat{w}_i = \frac{1}{\hat{\pi}_{i1} \times \hat{\pi}_{i2} \times \cdots \times \hat{\pi}_{in}}\]

Step 3: Weighted analysis

• Fit GEE with weights \(\hat{w}_i\)

IPW: Simple Example

# Using MAR data - estimate dropout at each time
dat_mar_ipw <- dat_mar %>%
  dplyr::select(subject, dropout_time, starts_with("Y_obs")) %>%
  pivot_longer(cols = starts_with("Y_obs"),
               names_to = "time",
               values_to = "y",
               values_drop_na = FALSE) %>%
  mutate(time = as.numeric(gsub("Y_obs.", "", time)),
         observed = !is.na(y)) %>%
  arrange(subject, time) %>%
  group_by(subject) %>%
  mutate(y_lag = lag(y)) %>%
  ungroup()

# Fit dropout model: probability of observing Y at time t given Y at t-1
dropout_model <- glm(observed ~ y_lag + time,
                     data = dat_mar_ipw %>% filter(time > 1),
                     family = binomial)

# Predict probability of staying observed
dat_mar_ipw$prob_stay <- predict(dropout_model,
                                 newdata = dat_mar_ipw,
                                 type = "response")

# Calculate weights
dat_mar_ipw <- dat_mar_ipw %>%
  mutate(prob_stay = if_else(time == 1, 1, prob_stay)) %>%
  group_by(subject) %>%
  mutate(weight = 1 / cumprod(prob_stay)) %>%
  ungroup() %>%
  filter(observed)

# Weighted GEE
fit_ipw <- geepack::geeglm(y ~ time,
                           data = dat_mar_ipw,
                           id = subject,
                           weights = weight,
                           corstr = "independence")

# Compare results
cat("Comparison of slopes:\n")

IPW corrects bias! (Full implementation in Chapter 18)

Methods Summary Table

Decision Framework

                Missing Data in Longitudinal Study
                            |
                            v
        +-------------------------------------------+
        | Missingness completely unrelated to Y?    |
        +-----------------------+-------------------+
                      |
        +-------------+----------------+
        |                              |
       YES (MCAR)                     NO
        |                              |
        v                              v
  +--------------+     +-------------------------------+
  |    MCAR      |     | Predictable from observed Y?  |
  |              |     +---------------+---------------+
  | Any method:  |             |
  | - ML/REML    |    +--------+----------+
  | - GEE        |   YES (MAR)           NO (MNAR)
  | - CC valid   |    |                   |
  +--------------+    v                   v
                +--------------+    +-----------------+
                |     MAR      |    |      MNAR       |
                |              |    |                 |
                | Methods:     |    | Joint models:   |
                | - ML/REML    |    | - Selection     |
                | - MI         |    | - Pattern-mix   |
                | - IPW-GEE    |    | - Sensitivity   |
                +--------------+    +-----------------+

Practical Recommendations

Sensitivity Analysis

Always conduct sensitivity analyses!

Under MAR:

• Try different covariance structures
• Try different imputation models
• Compare ML vs MI vs IPW

Under MNAR (or suspicion thereof):

• Specify multiple plausible MNAR models
• Report range of estimates
• Use pattern-mixture models with different assumptions
• Tipping point analysis

Reporting Checklist

Describe missingness patterns:

• Table showing \(n\) observed at each time
• Reasons for missing data (if known)
• Dropout vs intermittent missingness

State assumption about missing data mechanism:

• MCAR, MAR, or MNAR
• Justification based on study design/subject matter

Describe method used:

• Complete-case, ML, MI, IPW, etc.
• Model specifications (mean + covariance)
• Software used

Reporting Checklist (continued)

Report sensitivity analyses:

• Multiple methods (if reasonable)
• Multiple models under MAR
• MNAR scenarios (if applicable)

Interpret appropriately:

• Acknowledge limitations
• Discuss plausibility of MAR
• If substantial dropout, acknowledge uncertainty

Example statement:

“We analyzed data using maximum likelihood under the assumption that data are missing at random (MAR), conditional on observed outcomes and baseline covariates. We modeled the covariance using an unstructured matrix and confirmed robustness using multiple imputation with \(m=20\) imputations.”

Common Pitfalls to Avoid

Common Pitfalls (continued)

Illustrative Workflow: Putting It All Together

Hypothetical trial: Depression treatment, 4 visits over 12 weeks

Main reasons for dropout (from exit interviews):

• Lack of improvement (40%)
• Side effects (30%)
• Moved/lost contact (20%)
• Improved and stopped (10%)

Illustrative Results Comparison (hypothetical numbers)

Conclusion: Robust effect under MAR; some sensitivity to MNAR assumptions but qualitative conclusion unchanged.

Key Takeaways

Three critical insights:

Key Takeaways (continued)

Practical wisdom:

• Assume MAR as default (not MCAR)
• Model covariance carefully under MAR
• Avoid LVCF except in rare justified cases
• Conduct sensitivity analyses always
• Report clearly what assumption you’re making
• Acknowledge uncertainty when substantial missing data

Looking Ahead: Chapter 18

Next lecture will cover:

Multiple Imputation in detail

• Proper imputation models
• Software implementation
• Combining results

Inverse Probability Weighting in detail

• Estimating weights
• Stabilized weights
• Doubly robust methods

Advanced Topics

• Joint models for outcome + missingness
• Pattern-mixture models
• Sensitivity analysis frameworks

Resources & Further Reading

Key Papers:

• Rubin (1976): MCAR/MAR/MNAR taxonomy
• Little & Rubin (2002): Statistical Analysis with Missing Data
• Laird (1988): Missing data in longitudinal studies
• Kenward & Molenberghs (1998): Methods for handling dropout

Critiques of LVCF:

• Molenberghs & Verbeke (2005, Chapter 27)
• Cook et al. (2004)
• Kenward & Molenberghs (2009)

Common Mistakes in Missing Data Analysis