1. Objective

This example demonstrates:

Primary analysis under MAR using MMRM / MI.
Sensitivity analysis under MNAR using a pattern-mixture model with delta adjustment.

Clinical example:

Drug vs Placebo
HbA1c measured at Baseline, Week 4, Week 8, and Week 12
Some patients discontinue because of adverse events
Week 12 HbA1c may therefore be missing

2. Packages

library(dplyr)

## 
## Attaching package: 'dplyr'

## The following objects are masked from 'package:stats':
## 
##     filter, lag

## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union

library(tidyr)
library(nlme)

## 
## Attaching package: 'nlme'

## The following object is masked from 'package:dplyr':
## 
##     collapse

library(mice)

## Warning: package 'mice' was built under R version 4.4.3

## 
## Attaching package: 'mice'

## The following object is masked from 'package:stats':
## 
##     filter

## The following objects are masked from 'package:base':
## 
##     cbind, rbind

library(emmeans)

## Welcome to emmeans.
## Caution: You lose important information if you filter this package's results.
## See '? untidy'

set.seed(123)

3. Simulate Clinical Trial Data

n <- 300

dat <- data.frame(
  ID = 1:n,
  TRT = sample(c("Placebo", "Drug"), n, replace = TRUE),
  Baseline = rnorm(n, 8.5, 0.8)
)

dat <- dat %>%
  mutate(
    Y4 = Baseline +
      ifelse(TRT == "Drug", -0.4, -0.1) +
      rnorm(n, 0, 0.35),

    Y8 = Baseline +
      ifelse(TRT == "Drug", -0.8, -0.2) +
      rnorm(n, 0, 0.35),

    Y12_true = Baseline +
      ifelse(TRT == "Drug", -1.1, -0.3) +
      rnorm(n, 0, 0.35),

    AE = rbinom(
      n,
      1,
      ifelse(TRT == "Drug", 0.30, 0.12)
    )
  )

4. Generate MAR Dropout

Missingness depends on observed information:

Treatment
Week 8 HbA1c
Adverse event status

Therefore, this can still satisfy MAR.

# Target approximately 10% Week 12 missingness.
# Calibrate the intercept while preserving the MAR mechanism.

target_missing <- 0.10

lp_without_intercept <-
  1.2 * dat$AE +
  0.35 * (dat$Y8 - 8) +
  0.3 * (dat$TRT == "Drug")

calibrate_fun <- function(intercept) {
  mean(plogis(intercept + lp_without_intercept)) - target_missing
}

dropout_intercept <- uniroot(
  calibrate_fun,
  interval = c(-10, 2)
)$root

drop_prob <- plogis(
  dropout_intercept + lp_without_intercept
)

dat$Dropout <- rbinom(n, 1, drop_prob)

dat$Y12 <- dat$Y12_true
dat$Y12[dat$Dropout == 1] <- NA

# Overall and treatment-specific missing rates
overall_missing_rate <- mean(is.na(dat$Y12))

missing_by_trt <- dat %>%
  group_by(TRT) %>%
  summarise(
    MissingRate = mean(is.na(Y12)),
    .groups = "drop"
  )

overall_missing_rate

## [1] 0.1133333

missing_by_trt

## # A tibble: 2 × 2
##   TRT     MissingRate
##   <chr>         <dbl>
## 1 Drug         0.130 
## 2 Placebo      0.0974

5. Primary Analysis: MMRM Under MAR

longdat <- dat %>%
  select(ID, TRT, Baseline, Y4, Y8, Y12) %>%
  pivot_longer(
    cols = c(Y4, Y8, Y12),
    names_to = "Visit",
    values_to = "HbA1c"
  )

longdat$Visit <- factor(
  longdat$Visit,
  levels = c("Y4", "Y8", "Y12")
)

mmrm_fit <- lme(
  HbA1c ~ TRT * Visit + Baseline,
  random = ~1 | ID,
  data = longdat,
  na.action = na.exclude,
  method = "REML"
)

summary(mmrm_fit)

## Linear mixed-effects model fit by REML
##   Data: longdat 
##        AIC      BIC    logLik
##   695.5062 738.3082 -338.7531
## 
## Random effects:
##  Formula: ~1 | ID
##         (Intercept)  Residual
## StdDev:  0.07478946 0.3437904
## 
## Fixed effects:  HbA1c ~ TRT * Visit + Baseline 
##                          Value  Std.Error  DF   t-value p-value
## (Intercept)         -0.5278772 0.13700433 562  -3.85300   1e-04
## TRTPlacebo           0.3280141 0.04064364 297   8.07049   0e+00
## VisitY8             -0.3711867 0.04023763 562  -9.22486   0e+00
## VisitY12            -0.6905802 0.04177966 562 -16.52910   0e+00
## Baseline             1.0140474 0.01568428 297  64.65373   0e+00
## TRTPlacebo:VisitY8   0.2681477 0.05616072 562   4.77465   0e+00
## TRTPlacebo:VisitY12  0.4720353 0.05802480 562   8.13506   0e+00
##  Correlation: 
##                     (Intr) TRTPlc VistY8 VstY12 Baseln TRTP:VY8
## TRTPlacebo          -0.165                                     
## VisitY8             -0.147  0.495                              
## VisitY12            -0.150  0.477  0.482                       
## Baseline            -0.977  0.013  0.000  0.008                
## TRTPlacebo:VisitY8   0.105 -0.691 -0.716 -0.345  0.000         
## TRTPlacebo:VisitY12  0.107 -0.669 -0.347 -0.720 -0.005  0.484  
## 
## Standardized Within-Group Residuals:
##          Min           Q1          Med           Q3          Max 
## -2.693948186 -0.664785935  0.001773989  0.645792024  3.297561182 
## 
## Number of Observations: 866
## Number of Groups: 300

# Week 12 Drug - Placebo contrast
mmrm_emm <- emmeans(
  mmrm_fit,
  ~ TRT | Visit
)

mmrm_contrasts <- contrast(
  mmrm_emm,
  method = "revpairwise"
)

# Request confidence intervals explicitly.
# Depending on the model, emmeans may name CI columns
# lower.CL/upper.CL or asymp.LCL/asymp.UCL.
mmrm_contrast_summary <- as.data.frame(
  summary(
    mmrm_contrasts,
    infer = c(TRUE, TRUE)
  )
)

mmrm_week12 <- mmrm_contrast_summary %>%
  filter(Visit == "Y12")

mmrm_week12

## Visit = Y12:
##  contrast        estimate         SE  df  lower.CL  upper.CL t.ratio p.value
##  Placebo - Drug 0.8000494 0.04317906 297 0.7150738 0.8850251  18.529  <.0001
## 
## Degrees-of-freedom method: containment 
## Confidence level used: 0.95

The coefficient TRTDrug in an MMRM with a TRT * Visit interaction is the treatment difference at the reference visit, not necessarily the Week 12 effect.

Therefore, the clinically relevant Week 12 effect should be obtained from a visit-specific contrast.

6. Primary Analysis: Multiple Imputation Under MAR

mi_dat <- dat %>%
  select(TRT, Baseline, Y4, Y8, AE, Y12)

mi_dat$TRT <- factor(
  mi_dat$TRT,
  levels = c("Placebo", "Drug")
)

imp <- mice(
  mi_dat,
  m = 20,
  method = "pmm",
  seed = 123,
  printFlag = FALSE
)

fit_mi <- with(
  imp,
  lm(Y12 ~ TRT + Baseline)
)

mi_pooled <- summary(
  pool(fit_mi),
  conf.int = TRUE
)

mi_pooled

##          term   estimate  std.error  statistic       df       p.value
## 1 (Intercept) -0.6302126 0.22022253  -2.861708 254.8366  4.563488e-03
## 2     TRTDrug -0.7989755 0.04241025 -18.839208 192.4937  1.465626e-45
## 3    Baseline  1.0385687 0.02554325  40.659228 262.3539 2.921101e-115
##        2.5 %     97.5 %   conf.low  conf.high
## 1 -1.0639005 -0.1965248 -1.0639005 -0.1965248
## 2 -0.8826240 -0.7153271 -0.8826240 -0.7153271
## 3  0.9882729  1.0888646  0.9882729  1.0888646

mi_week12 <- mi_pooled %>%
  filter(term == "TRTDrug") %>%
  transmute(
    Method = "MI under MAR",
    Estimate = estimate,
    SE = std.error,
    Lower95 = `2.5 %`,
    Upper95 = `97.5 %`
  )

mi_week12

##         Method   Estimate         SE   Lower95    Upper95
## 2 MI under MAR -0.7989755 0.04241025 -0.882624 -0.7153271

The MI analysis is a Week 12 ANCOVA. Because Placebo is the reference group, TRTDrug directly estimates the adjusted Drug - Placebo difference at Week 12.

7. Compare MMRM and MI Week 12 Treatment Effects

The two methods should be compared using the same Week 12 contrast, rather than comparing the raw TRTDrug coefficient from the MMRM output.

# Identify the confidence-interval column names returned by emmeans
lower_col <- intersect(
  c("lower.CL", "asymp.LCL"),
  names(mmrm_week12)
)[1]

upper_col <- intersect(
  c("upper.CL", "asymp.UCL"),
  names(mmrm_week12)
)[1]

if (is.na(lower_col) || is.na(upper_col)) {
  stop("Confidence interval columns were not found in the emmeans output.")
}

mmrm_compare <- data.frame(
  Method = "MMRM under MAR",
  Estimate = mmrm_week12$estimate,
  SE = mmrm_week12$SE,
  Lower95 = mmrm_week12[[lower_col]],
  Upper95 = mmrm_week12[[upper_col]]
)

comparison_table <- bind_rows(
  mmrm_compare,
  mi_week12
)

comparison_table

##              Method   Estimate         SE    Lower95    Upper95
## ...1 MMRM under MAR  0.8000494 0.04317906  0.7150738  0.8850251
## 2      MI under MAR -0.7989755 0.04241025 -0.8826240 -0.7153271

Interpretation:

MMRM uses all available repeated outcomes and estimates the Week 12 Drug - Placebo contrast from the longitudinal model.
MI imputes missing Week 12 outcomes using observed information and then fits a Week 12 ANCOVA.
Both are MAR-based approaches, but the estimates do not have to be exactly identical because the model structures differ.
Similar estimates support consistency of the primary analysis.

8. MNAR Sensitivity Analysis: Delta Adjustment

Under MNAR, assume that patients with missing Week 12 HbA1c have outcomes worse than predicted under MAR.

For example:

\[ Y_{MNAR} = Y_{MAR} + \delta \]

where a positive delta means higher, and therefore worse, HbA1c.

run_delta <- function(delta, imp_object, original_data) {

  missing_y12 <- is.na(original_data$Y12)

  estimates <- numeric(imp_object$m)
  variances <- numeric(imp_object$m)

  for (i in 1:imp_object$m) {

    temp <- complete(imp_object, i)

    # Keep Placebo as the reference group
    temp$TRT <- factor(
      temp$TRT,
      levels = c("Placebo", "Drug")
    )

    # MNAR delta adjustment:
    # higher HbA1c is worse, so positive delta is unfavorable
    # to patients with originally missing Week 12 outcomes
    temp$Y12[missing_y12] <-
      temp$Y12[missing_y12] + delta

    fit <- lm(
      Y12 ~ TRT + Baseline,
      data = temp
    )

    estimates[i] <- unname(coef(fit)["TRTDrug"])

    vc <- vcov(fit)

    if (!"TRTDrug" %in% rownames(vc)) {
      stop(
        "TRTDrug coefficient was not found. ",
        "Check treatment factor levels."
      )
    }

    variances[i] <- vc["TRTDrug", "TRTDrug"]
  }

  # Rubin's rules
  Qbar <- mean(estimates)
  Ubar <- mean(variances)
  B <- var(estimates)

  total_var <-
    Ubar + (1 + 1 / imp_object$m) * B

  SE <- sqrt(total_var)

  data.frame(
    Delta = delta,
    Estimate = Qbar,
    SE = SE,
    Lower95 = Qbar - 1.96 * SE,
    Upper95 = Qbar + 1.96 * SE
  )
}

9. Test Several MNAR Assumptions

delta_values <- c(
  0,
  0.2,
  0.4,
  0.6,
  0.8,
  1.0
)

results <- bind_rows(
  lapply(
    delta_values,
    run_delta,
    imp_object = imp,
    original_data = mi_dat
  )
)

results

##   Delta   Estimate         SE    Lower95    Upper95
## 1   0.0 -0.7989755 0.04241025 -0.8820996 -0.7158514
## 2   0.2 -0.7925499 0.04292054 -0.8766742 -0.7084257
## 3   0.4 -0.7861243 0.04464960 -0.8736376 -0.6986111
## 4   0.6 -0.7796988 0.04746444 -0.8727290 -0.6866685
## 5   0.8 -0.7732732 0.05118622 -0.8735982 -0.6729482
## 6   1.0 -0.7668476 0.05563324 -0.8758887 -0.6578064

Interpretation:

The treatment coefficient is Drug - Placebo because Placebo is explicitly set as the reference group. For HbA1c, a negative estimate favors Drug.

Delta = 0: MAR assumption
Delta = 0.2: missing HbA1c values are assumed to be 0.2% worse than MAR predictions
Delta = 0.5: 0.5% worse
Larger delta values represent increasingly unfavorable MNAR assumptions

10. Interpretation

The overall strategy is:

\[ \text{Primary Analysis} \rightarrow \text{MMRM or MI under MAR} \]

followed by:

\[ \text{Sensitivity Analysis} \rightarrow \text{Pattern-Mixture Model + Delta Adjustment} \]

If the treatment conclusion remains similar across clinically plausible delta values, the result is considered robust to departures from MAR.

If the conclusion changes only after a sufficiently large delta, that value can be interpreted as a tipping point.

11. Key Takeaway

MAR allows missingness to depend on observed covariates and prior outcomes.

For the primary analysis, MMRM or MI can therefore be used under a reasonable MAR assumption.

Because MAR cannot be verified from the observed data alone, MNAR sensitivity analyses such as pattern-mixture models with delta adjustment are used to assess robustness.

Missing Data Analysis: MAR and MNAR Sensitivity

Dian He

2026-08-24