The V Lab

Scope and safety. This is a methods tutorial, not a substitute for a prespecified protocol, subject-matter knowledge, statistical review, research ethics review, or jurisdiction-specific guidance. Every example is synthetic. The code favors transparency over production-scale software. A method does not repair an incoherent question, poor measurement, or an unsupported causal assumption.

How to use this tutorial

The foundations guides explain frequency measures, elementary study designs, 2-by-2 tables, basic regression, sampling, and introductory causal diagrams. This guide starts where they stop. Its organizing sequence is:

question →\rightarrow estimand →\rightarrow identification assumptions →\rightarrow estimator →\rightarrow diagnostics →\rightarrow sensitivity analysis →\rightarrow interpretation.

Sections 2–4 develop time-to-event estimands. Sections 5–9 build causal estimators for baseline, longitudinal, incomplete, and nonrepresentative data. Sections 10–11 address self-controlled designs and residual bias. Section 12 integrates the methods into a reproducible workflow.

Most chunks require only base R, stats, and graphics. Survival examples use the recommended survival package when it is installed; each such chunk is guarded and prints an informative message otherwise. The guide deliberately does not implement specialized estimators such as Fine–Gray regression or machine- learning nuisance models by hand.

Learning objectives

After completing the tutorial, you should be able to:

  • distinguish a scientific estimand from a fitted-model coefficient;
  • estimate and interpret Kaplan–Meier survival, Nelson–Aalen cumulative hazard, restricted mean survival time (RMST), and cumulative incidence;
  • diagnose proportional-hazards problems and represent a prespecified time-varying effect;
  • construct and audit propensity, censoring, longitudinal treatment, and transportability weights;
  • implement an augmented inverse-probability weighted (AIPW) estimator and state precisely what “double robustness” does and does not mean;
  • explain treatment–confounder feedback and estimate a marginal structural model;
  • distinguish a transparent teaching demonstration of multiple imputation from a defensible production analysis;
  • analyze a case-crossover study using conditional logistic regression; and
  • conduct deterministic and probabilistic misclassification analyses and interpret negative controls cautiously.

1. Estimands and advanced assumptions: a bridge

An estimand is the quantity the study aims to learn. An estimator is the rule applied to data. An estimate is the resulting number. “Fit a Cox model” or “use propensity scores” specifies an estimator family, not a scientific target.

For treatment A∈{0,1}A\in\{0,1\}, baseline covariates WW, potential outcome YaY^a, and horizon τ\tau, common marginal estimands include:

ATERD=E(Y1)−E(Y0),ATERR=E(Y1)/E(Y0),ΔRMST(τ)=E{min⁡(T1,τ)}−E{min⁡(T0,τ)}. \begin{aligned} \text{ATE}_{RD} &= E(Y^1)-E(Y^0),\\ \text{ATE}_{RR} &= E(Y^1)/E(Y^0),\\ \Delta_{RMST}(\tau) &= E\{\min(T^1,\tau)\}-E\{\min(T^0,\tau)\}. \end{aligned}

The averaging population matters. The average treatment effect (ATE), average effect among the treated (ATT), overlap-population effect (ATO), and target- population effect (TATE) need not agree when effects vary across people.

Scientific question Example estimand Essential time/population detail
What is the 5-year event burden under each strategy? P(Ta≤5)P(T^a\le5) and a risk difference Competing events and censoring must be defined
How much event-free time is gained through year 5? ΔRMST(5)\Delta_{RMST}(5) Horizon 5; units are time
What is the effect of sustained treatment? E(Ya‾=1‾)−E(Ya‾=0‾)E(Y^{\bar a=\bar1})-E(Y^{\bar a=\bar0}) Treatment history and adherence rule
What would the effect be in another population? ET(Y1−Y0)E_T(Y^1-Y^0) Explicit target population TT

1.1 Identification is not estimation

Observed data identify a causal mean under assumptions such as:

  1. Consistency and well-defined interventions: a person’s observed outcome under the treatment received equals the corresponding potential outcome; materially different treatment versions are addressed.
  2. Exchangeability: conditional on the prespecified covariates, treatment (or censoring/selection) is independent of the relevant potential outcomes.
  3. Positivity: every treatment or observation pattern required by the estimand has positive probability within relevant covariate histories.
  4. No relevant interference, or an explicit interference structure: one person’s treatment does not alter another person’s outcome unless the estimand models that spillover.
  5. Correct temporal ordering and measurement: confounders precede the treatment decision they are intended to adjust, and variables adequately represent the causal concepts.

For longitudinal treatment, exchangeability and positivity are sequential: they must hold at every treatment time conditional on the observed past. For transportability, an additional exchangeability-over-study-selection condition is needed. None of these conditions is established by a good model fit.

The ICH E9(R1) estimand framework formalizes alignment among the clinical question, intercurrent events, analysis, and sensitivity analyses. Its discipline is valuable beyond regulatory trials.

1.2 An estimand card

Before coding, record:

  • population and eligibility;
  • treatment/exposure strategies, including versions and time zero;
  • outcome and horizon;
  • handling of death, discontinuation, rescue therapy, and other intercurrent events;
  • population-level contrast and effect scale;
  • identification assumptions and adjustment set;
  • estimator, uncertainty method, diagnostics, and sensitivity analyses.

Correctness trap: a hazard ratio is not a risk ratio. It compares instantaneous hazards among people still event-free, a selected risk set that changes over time. Even without confounding, the Cox coefficient is generally not the same estimand as a fixed-horizon marginal risk ratio or RMST difference.

1.3 Marginal is not a synonym for unadjusted

Standardization averages conditional outcome predictions over a specified covariate distribution:

μa=EW{E(Y∣A=a,W)}. \mu_a=E_W\{E(Y\mid A=a,W)\}.

Thus an adjusted regression can produce a marginal standardized contrast, while its treatment coefficient is usually conditional on modeled covariates. A crude contrast averages over the treatment groups’ observed covariate distributions; under confounding those are different distributions and the crude contrast is not the target-population causal marginal effect. Consequently, “adjusted = conditional” and “unadjusted = marginal” are not universal synonyms.

Odds ratios and hazard ratios are noncollapsible: a conditional ratio can differ from the corresponding marginal ratio even when treatment is randomized and there is no confounding. The difference alone is not evidence of bias.

# A deterministic covariate distribution removes Monte Carlo noise. Treatment is
# conceptually randomized; W is prognostic but not a confounder.
n_standard <- 100000
w_standard <- qnorm((seq_len(n_standard) - 0.5) / n_standard)
conditional_log_or <- log(2)
p0_standard <- plogis(-1.20 + 1.10 * w_standard)
p1_standard <- plogis(-1.20 + conditional_log_or + 1.10 * w_standard)
mu0 <- mean(p0_standard)
mu1 <- mean(p1_standard)
marginal_or <- (mu1 / (1 - mu1)) / (mu0 / (1 - mu0))

data.frame(
  conditional_OR = exp(conditional_log_or),
  standardized_marginal_OR = marginal_or,
  standardized_risk_A0 = mu0,
  standardized_risk_A1 = mu1,
  standardized_risk_difference = mu1 - mu0
)

2. Kaplan–Meier, Nelson–Aalen, and RMST

Let TT be event time and CC censoring time. We observe T̃=min⁡(T,C)\tilde T=\min(T,C) and Δ=I(T≤C)\Delta=I(T\le C). Standard nonparametric survival estimators require censoring to be independent of the event time, possibly after conditioning or weighting, and require event/censoring times to be measured correctly.

At ordered event times tjt_j, with djd_j events among njn_j at risk, the Kaplan–Meier estimator is

Ŝ(t)=∏tj≤t(1−djnj), \widehat S(t)=\prod_{t_j\le t}\left(1-\frac{d_j}{n_j}\right),

while the Nelson–Aalen cumulative-hazard estimator is

Ĥ(t)=∑tj≤tdjnj. \widehat H(t)=\sum_{t_j\le t}\frac{d_j}{n_j}.

exp⁡{−Ĥ(t)}\exp\{-\widehat H(t)\} and Ŝ(t)\widehat S(t) are closely related but are not identical finite-sample estimators.

set.seed(1021)
n_surv <- 1200
surv_dat <- data.frame(
  id = seq_len(n_surv),
  A = rbinom(n_surv, 1, 0.5),
  W = rnorm(n_surv)
)

# Treatment lowers the event hazard in this teaching data-generating process.
event_time <- rexp(n_surv, rate = 0.11 * exp(-0.45 * surv_dat$A +
                                              0.30 * surv_dat$W))
censor_time <- rexp(n_surv, rate = 0.035)
administrative_end <- 8
surv_dat$time <- pmin(event_time, censor_time, administrative_end)
surv_dat$status <- as.integer(event_time <= censor_time &
                                event_time <= administrative_end)

with(surv_dat, table(treatment = A, event = status))
##          event
## treatment   0   1
##         0 274 315
##         1 379 232

2.1 Kaplan–Meier survival

if (!has_survival) {
  cat("Install the recommended 'survival' package to run this example.\n")
} else {
  km_fit <- survival::survfit(
    survival::Surv(time, status) ~ A,
    data = surv_dat,
    conf.type = "log-log"
  )
  km_at <- summary(km_fit, times = c(3, 5, 8), extend = TRUE)
  km_table <- data.frame(
    group = km_at$strata,
    time = km_at$time,
    survival = km_at$surv,
    lower = km_at$lower,
    upper = km_at$upper
  )
  print(km_table, row.names = FALSE)

  plot(km_fit, col = c("#C44E52", "#2A6F97"), lwd = 2,
       xlab = "Follow-up time", ylab = "Event-free survival",
       mark.time = TRUE, conf.int = FALSE)
  legend("bottomleft", legend = c("A = 0", "A = 1"),
         col = c("#C44E52", "#2A6F97"), lwd = 2, bty = "n")
}
##  group time survival  lower  upper
##    A=0    3   0.7173 0.6781 0.7526
##    A=0    5   0.5813 0.5386 0.6216
##    A=0    8   0.4079 0.3651 0.4502
##    A=1    3   0.8279 0.7947 0.8561
##    A=1    5   0.7189 0.6798 0.7541
##    A=1    8   0.5736 0.5302 0.6145

The confidence interval describes sampling uncertainty under the estimator’s assumptions. It does not cover informative censoring, outcome misclassification, unmeasured confounding, or an inappropriate time zero.

2.2 Nelson–Aalen cumulative hazard

if (!has_survival) {
  cat("The Nelson--Aalen example was skipped because 'survival' is unavailable.\n")
} else {
  na_fit <- survival::survfit(
    survival::Surv(time, status) ~ A,
    data = surv_dat,
    stype = 2,
    ctype = 1
  )
  na_at <- summary(na_fit, times = 5, extend = TRUE)
  data.frame(
    group = na_at$strata,
    time = na_at$time,
    nelson_aalen_H = na_at$cumhaz,
    exp_minus_H = na_at$surv
  )
}

The cumulative hazard is not a cumulative probability and can exceed 1. Its increment dj/njd_j/n_j is additive; survival is multiplicative.

2.3 Restricted mean survival time

RMST through τ\tau is the area under the survival curve:

RMST(τ)=∫0τS(t)dt. RMST(\tau)=\int_0^\tau S(t)\,dt.

It answers an absolute, horizon-specific question in units of time and does not require proportional hazards. Choose τ\tau from scientific relevance and common follow-up support, not after inspecting where curves look most different.

rmst_from_survfit <- function(fit, tau) {
  stopifnot(length(tau) == 1, is.finite(tau), tau > 0)
  keep <- fit$time < tau
  interval_edges <- c(0, fit$time[keep], tau)
  interval_survival <- c(1, fit$surv[keep])
  sum(diff(interval_edges) * interval_survival)
}

rmst_difference <- function(data, tau = 5) {
  fit0 <- survival::survfit(
    survival::Surv(time, status) ~ 1,
    data = data[data$A == 0, ]
  )
  fit1 <- survival::survfit(
    survival::Surv(time, status) ~ 1,
    data = data[data$A == 1, ]
  )
  c(A0 = rmst_from_survfit(fit0, tau),
    A1 = rmst_from_survfit(fit1, tau),
    difference = rmst_from_survfit(fit1, tau) -
      rmst_from_survfit(fit0, tau))
}

if (!has_survival) {
  cat("The RMST example was skipped because 'survival' is unavailable.\n")
} else {
  rmst_point <- rmst_difference(surv_dat, tau = 5)
  set.seed(1022)
  rmst_boot <- replicate(300, {
    index <- sample.int(nrow(surv_dat), replace = TRUE)
    rmst_difference(surv_dat[index, ], tau = 5)["difference"]
  })
  data.frame(
    estimand = c("RMST A=0", "RMST A=1", "RMST difference (A=1 minus A=0)"),
    estimate = unname(rmst_point),
    lower = c(NA, NA, unname(quantile(rmst_boot, 0.025))),
    upper = c(NA, NA, unname(quantile(rmst_boot, 0.975)))
  )
}

The percentile interval above resamples the complete person-level record. For an observational contrast, this unadjusted RMST difference is associational unless a design and adjustment strategy justify causal interpretation.

3. Cox diagnostics and time-varying effects

The Cox model specifies

h(t∣X)=h0(t)exp⁡(βTX). h(t\mid X)=h_0(t)\exp(\beta^T X).

For a binary treatment, exp⁡(βA)\exp(\beta_A) is a conditional hazard ratio under the model. A causal interpretation additionally requires a causal design and identification assumptions; proportional hazards alone is not enough.

if (!has_survival) {
  cat("The Cox example was skipped because 'survival' is unavailable.\n")
} else {
  cox_fit <- survival::coxph(
    survival::Surv(time, status) ~ A + W,
    data = surv_dat,
    x = TRUE
  )
  cox_summary <- summary(cox_fit)
  print(cbind(
    HR = exp(stats::coef(cox_fit)),
    exp(confint(cox_fit))
  ))

  ph_test <- survival::cox.zph(cox_fit, transform = "km")
  print(ph_test)
  plot(ph_test, var = "A", resid = TRUE,
       xlab = "Transformed time", ylab = "Scaled Schoenfeld residual for A")
  abline(h = 0, lty = 2, col = "grey40")
}
##      HR  2.5 % 97.5 %
## A 0.610 0.5147 0.7229
## W 1.354 1.2450 1.4724
##        chisq df    p
## A      0.103  1 0.75
## W      1.211  1 0.27
## GLOBAL 1.304  2 0.52

cox.zph() tests whether scaled Schoenfeld residuals show a systematic relation with time. Inspect the plot and the global test; do not reduce diagnostics to a single pp-value. A small study can miss important nonproportionality, while a large study can detect a clinically trivial departure. Also inspect influential observations, functional forms, event counts, and risk-set support.

3.1 A prespecified piecewise treatment effect

The next simulation has a protective early treatment effect and a harmful late effect. The cut at time 3 is part of the data-generating story; in a real protocol, such a cut should be clinically justified and prespecified.

set.seed(1031)
n_np <- 1800
np_dat <- data.frame(
  id = seq_len(n_np),
  A = rbinom(n_np, 1, 0.5),
  W = rnorm(n_np)
)
base_rate <- 0.12 * exp(0.25 * np_dat$W)
early_wait <- rexp(n_np, rate = base_rate * exp(-0.80 * np_dat$A))
late_wait <- rexp(n_np, rate = base_rate * exp(0.35 * np_dat$A))
event_np <- ifelse(early_wait <= 3, early_wait, 3 + late_wait)
censor_np <- rexp(n_np, rate = 0.025)
np_dat$time <- pmin(event_np, censor_np, 7)
np_dat$status <- as.integer(event_np <= censor_np & event_np <= 7)

if (!has_survival) {
  cat("The time-varying Cox example was skipped because 'survival' is unavailable.\n")
} else {
  # survSplit() inspects formula specials; survival 3.8.x requires Surv() to be
  # visible without a namespace qualifier on this formula's left-hand side.
  suppressPackageStartupMessages(library(survival))
  np_cox <- survival::coxph(
    survival::Surv(time, status) ~ A + W,
    data = np_dat,
    x = TRUE
  )
  print(survival::cox.zph(np_cox))

  np_long <- survival::survSplit(
    Surv(time, status) ~ .,
    data = np_dat,
    cut = 3,
    episode = "period",
    id = "split_id"
  )
  np_long$A_early <- np_long$A * (np_long$period == 1)
  np_long$A_late <- np_long$A * (np_long$period == 2)
  piecewise_fit <- survival::coxph(
    Surv(tstart, time, status) ~ A_early + A_late + W,
    data = np_long,
    cluster = id
  )
  print(cbind(
    HR = exp(stats::coef(piecewise_fit)),
    exp(confint(piecewise_fit))
  ))
}
##         chisq df       p
## A      63.330  1 1.7e-15
## W       0.273  1     0.6
## GLOBAL 63.577  2 1.6e-14
##             HR  2.5 % 97.5 %
## A_early 0.4436 0.3596 0.5473
## A_late  1.5830 1.3285 1.8864
## W       1.3276 1.2438 1.4171

The piecewise model estimates two hazard ratios; it does not turn either into a risk ratio. Flexible time interactions, stratification, or horizon-specific survival/RMST contrasts may be preferable depending on the question.

4. Competing risks and the Aalen–Johansen estimator

A competing event prevents the event of interest from subsequently occurring. For example, death from another cause competes with disease-specific death. Treating a competing event as ordinary independent censoring changes the target and generally overestimates the real-world cumulative incidence.

For cause 1, the Aalen–Johansen cumulative-incidence estimator is

F̂1(t)=∑tj≤tŜ(tj−)d1jnj, \widehat F_1(t)=\sum_{t_j\le t}\widehat S(t_j-)\frac{d_{1j}}{n_j},

where survival from all event types updates as Ŝ(tj)=Ŝ(tj−)(1−dj/nj)\widehat S(t_j)=\widehat S(t_j-)(1-d_j/n_j).

set.seed(1041)
n_cr <- 2200
comp_dat <- data.frame(
  id = seq_len(n_cr),
  A = rbinom(n_cr, 1, 0.5),
  W = rnorm(n_cr)
)
t_target <- rexp(n_cr, 0.085 * exp(-0.35 * comp_dat$A + 0.25 * comp_dat$W))
t_compete <- rexp(n_cr, 0.070 * exp(0.20 * comp_dat$A + 0.20 * comp_dat$W))
t_censor <- rexp(n_cr, 0.025)
comp_dat$time <- pmin(t_target, t_compete, t_censor, 7)
comp_dat$status <- ifelse(
  t_target <= pmin(t_compete, t_censor, 7), 1L,
  ifelse(t_compete <= pmin(t_target, t_censor, 7), 2L, 0L)
)

aj_two_cause <- function(time, status, tau = Inf) {
  event_times <- sort(unique(time[status %in% c(1L, 2L) & time <= tau]))
  state <- data.frame(time = 0, survival = 1, cif_target = 0,
                      cif_competing = 0)
  S <- 1
  F1 <- 0
  F2 <- 0
  for (tt in event_times) {
    risk <- sum(time >= tt)
    d1 <- sum(time == tt & status == 1L)
    d2 <- sum(time == tt & status == 2L)
    F1 <- F1 + S * d1 / risk
    F2 <- F2 + S * d2 / risk
    S <- S * (1 - (d1 + d2) / risk)
    state <- rbind(state, data.frame(time = tt, survival = S,
                                     cif_target = F1,
                                     cif_competing = F2))
  }
  state
}

aj_by_group <- lapply(0:1, function(a) {
  out <- aj_two_cause(comp_dat$time[comp_dat$A == a],
                      comp_dat$status[comp_dat$A == a], tau = 7)
  out$A <- a
  out
})

# Cross-check the manual estimator against survival's multistate representation.
# A factor with censoring as its first level avoids treating numeric cause codes as
# an ambiguous multistate status.
if (has_survival) {
  comp_dat$status_factor <- factor(
    comp_dat$status, levels = c(0, 1, 2),
    labels = c("censor", "target", "competing")
  )
  aj_check <- survival::survfit(
    survival::Surv(time, status_factor) ~ A,
    data = comp_dat
  )
  aj_check_at_7 <- summary(aj_check, times = 7, extend = TRUE)
  package_target_cif <- aj_check_at_7$pstate[, "target"]
  manual_target_cif <- vapply(
    aj_by_group, function(x) tail(x$cif_target, 1), numeric(1)
  )
  cat("Maximum state-probability row-sum deviation from 1:",
      max(abs(rowSums(aj_check$pstate) - 1)), "\n")
  cat("Maximum package-versus-hand target-CIF difference at time 7:",
      max(abs(package_target_cif - manual_target_cif)), "\n")
}
## Maximum state-probability row-sum deviation from 1: 3.553e-15
## Maximum package-versus-hand target-CIF difference at time 7: 8.327e-16
plot(aj_by_group[[1]]$time, aj_by_group[[1]]$cif_target,
     type = "s", lwd = 2, col = "#C44E52", xlim = c(0, 7),
     ylim = c(0, max(vapply(aj_by_group, function(x) max(x$cif_target), 0))),
     xlab = "Follow-up time", ylab = "Cumulative incidence of target event")
lines(aj_by_group[[2]]$time, aj_by_group[[2]]$cif_target,
      type = "s", lwd = 2, col = "#2A6F97")
legend("topleft", c("A = 0", "A = 1"),
       col = c("#C44E52", "#2A6F97"), lwd = 2, bty = "n")

vapply(aj_by_group, function(x) tail(x$cif_target, 1), numeric(1))
## [1] 0.3781 0.2725

4.1 Why 1−1-Kaplan–Meier is wrong here

if (!has_survival) {
  cat("The naive competing-risk comparison was skipped because 'survival' is unavailable.\n")
} else {
  tau_cr <- 7
  comparison <- lapply(0:1, function(a) {
    d <- comp_dat[comp_dat$A == a, ]
    aj <- aj_two_cause(d$time, d$status, tau = tau_cr)
    target_only_km <- survival::survfit(
      survival::Surv(time, status == 1L) ~ 1,
      data = d
    )
    km_s <- summary(target_only_km, times = tau_cr, extend = TRUE)$surv
    data.frame(
      A = a,
      Aalen_Johansen_CIF = tail(aj$cif_target, 1),
      naive_one_minus_KM = 1 - km_s
    )
  })
  do.call(rbind, comparison)
}

The cause-specific hazard, subdistribution hazard, cumulative incidence, and fixed-horizon risk are distinct. A Fine–Gray coefficient is a subdistribution hazard ratio, not a risk ratio and not a generic “correction” for competing events. Causal interpretation also depends on whether the target is a total effect in a world where competing events occur, or a hypothetical/direct effect under an intervention on the competing event. See Andersen et al. and Young et al..

5. Propensity scores, IPTW, and balance

The propensity score e(W)=P(A=1∣W)e(W)=P(A=1\mid W) is a treatment-assignment probability, not a disease risk. For the ATE, inverse-probability-of-treatment weights (IPTW) are

wiATE=Aiê(Wi)+1−Ai1−ê(Wi). w_i^{ATE}=\frac{A_i}{\widehat e(W_i)}+ \frac{1-A_i}{1-\widehat e(W_i)}.

The weighted sample aims to make measured baseline confounders independent of treatment. This identifies an ATE only with consistency, conditional exchangeability, positivity, and an adequate treatment model. Different weights target different populations:

Target Treated weight Untreated weight
ATE 1/e1/e 1/(1−e)1/(1-e)
ATT 11 e/(1−e)e/(1-e)
ATO (overlap) 1−e1-e ee

Changing the weights changes the estimand. Overlap weights can improve stability, but their answer concerns the population with clinical equipoise rather than the full study population.

set.seed(1051)
n_ps <- 3000
ps_dat <- data.frame(
  W1 = rnorm(n_ps),
  W2 = rbinom(n_ps, 1, 0.45),
  W3 = rnorm(n_ps)
)
ps_dat$e_true <- plogis(-0.25 + 0.80 * ps_dat$W1 + 0.65 * ps_dat$W2 -
                          0.45 * ps_dat$W3 + 0.35 * ps_dat$W1 * ps_dat$W2)
ps_dat$A <- rbinom(n_ps, 1, ps_dat$e_true)
ps_dat$p0_true <- plogis(-1.55 + 0.80 * ps_dat$W1 + 0.50 * ps_dat$W2 +
                           0.35 * ps_dat$W3)
ps_dat$p1_true <- plogis(-1.55 + 0.65 + 0.80 * ps_dat$W1 +
                           0.50 * ps_dat$W2 + 0.35 * ps_dat$W3 -
                           0.25 * ps_dat$W1)
ps_dat$Y <- rbinom(n_ps, 1,
                   ifelse(ps_dat$A == 1, ps_dat$p1_true, ps_dat$p0_true))

finite_sample_true_ate <- mean(ps_dat$p1_true - ps_dat$p0_true)
c(observed_treatment_prevalence = mean(ps_dat$A),
  finite_sample_true_ATE = finite_sample_true_ate)
## observed_treatment_prevalence        finite_sample_true_ATE
##                        0.4947                        0.1085

5.1 Construct and audit weights

Treatment-model covariates should be chosen using temporal and causal knowledge, not automated pp-value selection. Include outcome causes needed for exchangeability, flexible functional forms, and relevant interactions. Do not include post-treatment variables.

ps_fit <- glm(A ~ W1 + W2 + W3 + W1:W2,
              family = binomial(), data = ps_dat)
ps_dat$e_hat <- validate_probability(predict(ps_fit, type = "response"),
                                     "propensity score")
ps_dat$w_ate <- with(ps_dat, A / e_hat + (1 - A) / (1 - e_hat))
treatment_prevalence <- mean(ps_dat$A)
ps_dat$w_ate_stabilized <- with(
  ps_dat,
  A * treatment_prevalence / e_hat +
    (1 - A) * (1 - treatment_prevalence) / (1 - e_hat)
)
ps_dat$w_att <- with(ps_dat, A + (1 - A) * e_hat / (1 - e_hat))
ps_dat$w_ato <- with(ps_dat, A * (1 - e_hat) + (1 - A) * e_hat)

weight_summary <- function(w) {
  c(mean = mean(w), sd = sd(w), min = min(w),
    q50 = unname(quantile(w, 0.50)),
    q95 = unname(quantile(w, 0.95)),
    q99 = unname(quantile(w, 0.99)), max = max(w),
    ESS = effective_sample_size(w))
}
rbind(ATE = weight_summary(ps_dat$w_ate),
      `ATE stabilized` = weight_summary(ps_dat$w_ate_stabilized),
      ATT = weight_summary(ps_dat$w_att),
      ATO = weight_summary(ps_dat$w_ato))
##                  mean     sd      min    q50    q95    q99     max  ESS
## ATE            2.0066 1.4397 1.003788 1.6183 4.1566 7.2082 31.5377 1981
## ATE stabilized 1.0032 0.7202 0.496540 0.8090 2.0705 3.5738 15.9370 1980
## ATT            0.9856 1.0529 0.019384 1.0000 2.0291 4.5226 30.5377 1401
## ATO            0.3966 0.2021 0.003773 0.3821 0.7594 0.8613  0.9683 2382
balance_variables <- c("W1", "W2", "W3")
balance_table <- data.frame(
  variable = balance_variables,
  unweighted = vapply(balance_variables, function(v) {
    weighted_smd(ps_dat[[v]], ps_dat$A)
  }, numeric(1)),
  ATE_weighted = vapply(balance_variables, function(v) {
    weighted_smd(ps_dat[[v]], ps_dat$A, ps_dat$w_ate)
  }, numeric(1))
)
balance_table
hist(ps_dat$e_hat[ps_dat$A == 0], breaks = 25, probability = TRUE,
     col = grDevices::adjustcolor("#C44E52", alpha.f = 0.45),
     border = "white", xlim = c(0, 1),
     xlab = "Estimated propensity score", main = "Treatment overlap")
hist(ps_dat$e_hat[ps_dat$A == 1], breaks = 25, probability = TRUE,
     col = grDevices::adjustcolor("#2A6F97", alpha.f = 0.45),
     border = "white", add = TRUE)
legend("topright", c("A = 0", "A = 1"),
       fill = c(grDevices::adjustcolor("#C44E52", alpha.f = 0.45),
                grDevices::adjustcolor("#2A6F97", alpha.f = 0.45)),
       bty = "n")

Audit overlap, the full weight distribution, effective sample size (ESS), and covariate distributions before and after weighting. Standardized mean differences (SMDs) are more useful than baseline hypothesis tests, but a threshold such as |SMD|<0.1|SMD|<0.1 is only a heuristic. Check nonlinear terms, interactions, variances, and plots as well. The helper here uses the pooled within-group weighted standard deviation as its denominator. Another common convention fixes the original unweighted pooled standard deviation for every before/after comparison; state which ruler you use. A mean stabilized weight near 1 is neither necessary for the unstabilized weights above nor proof of a valid model.

weighted_effect <- function(y, a, w) {
  r1 <- weighted_mean(y[a == 1], w[a == 1])
  r0 <- weighted_mean(y[a == 0], w[a == 0])
  c(risk_A1 = r1, risk_A0 = r0, risk_difference = r1 - r0,
    risk_ratio = r1 / r0)
}

iptw_results <- rbind(
  ATE = weighted_effect(ps_dat$Y, ps_dat$A, ps_dat$w_ate),
  `ATE stabilized` = weighted_effect(ps_dat$Y, ps_dat$A,
                                     ps_dat$w_ate_stabilized),
  ATT = weighted_effect(ps_dat$Y, ps_dat$A, ps_dat$w_att),
  ATO = weighted_effect(ps_dat$Y, ps_dat$A, ps_dat$w_ato)
)

cut_points <- quantile(ps_dat$w_ate, c(0.01, 0.99))
w_ate_truncated <- pmin(pmax(ps_dat$w_ate, cut_points[1]), cut_points[2])
iptw_results <- rbind(
  iptw_results,
  `ATE weights truncated at empirical 1st/99th percentiles` =
    weighted_effect(ps_dat$Y, ps_dat$A, w_ate_truncated)
)
iptw_results
##                                                         risk_A1 risk_A0 risk_difference
## ATE                                                      0.3362  0.2442         0.09198
## ATE stabilized                                           0.3362  0.2442         0.09198
## ATT                                                      0.3666  0.3035         0.06310
## ATO                                                      0.3304  0.2305         0.09981
## ATE weights truncated at empirical 1st/99th percentiles  0.3367  0.2319         0.10478
##                                                         risk_ratio
## ATE                                                          1.377
## ATE stabilized                                               1.377
## ATT                                                          1.208
## ATO                                                          1.433
## ATE weights truncated at empirical 1st/99th percentiles      1.452

Stabilization multiplies all treated ATE weights by P(A=1)P(A=1) and all untreated weights by P(A=0)P(A=0). Therefore the arm-normalized (Hájek) risks used here, their balance, and their contrast are unchanged; the stabilized-weight mean is near 1 and its scale can improve numerical behavior. Stabilization does not cure nonpositivity. It would not be innocuous for every unnormalized estimator, so the estimating equation must be stated.

Truncation can trade variance for bias and can blur or change the effective target population. Report its rule, the observations affected, and untruncated results; do not choose the cutoff that produces the preferred answer.

5.2 Uncertainty must include weight estimation

The next bootstrap refits the propensity model in every resample. Treating estimated weights as fixed, or using the default model-based standard error from glm(..., weights=), generally misses parts of the uncertainty.

iptw_ate_rd <- function(data) {
  fit <- glm(A ~ W1 + W2 + W3 + W1:W2,
             family = binomial(), data = data)
  e <- validate_probability(predict(fit, type = "response"),
                            "bootstrap propensity score")
  w <- with(data, A / e + (1 - A) / (1 - e))
  unname(weighted_effect(data$Y, data$A, w)["risk_difference"])
}

set.seed(1052)
iptw_boot <- replicate(300, {
  index <- sample.int(nrow(ps_dat), replace = TRUE)
  iptw_ate_rd(ps_dat[index, ])
})
data.frame(
  estimate = iptw_ate_rd(ps_dat),
  bootstrap_se = sd(iptw_boot),
  lower = unname(quantile(iptw_boot, 0.025)),
  upper = unname(quantile(iptw_boot, 0.975))
)

See Cole and Hernán for weight construction, Austin for balance diagnostics, and Petersen et al. for positivity.

6. AIPW and double robustness

Outcome standardization models ma(W)=E(Y∣A=a,W)m_a(W)=E(Y\mid A=a,W). IPTW models treatment. The augmented inverse-probability weighted estimator combines both:

ψ̂AIPW=1n∑i[m̂1(Wi)−m̂0(Wi)+Ai{Yi−m̂1(Wi)}ê(Wi)−(1−Ai){Yi−m̂0(Wi)}1−ê(Wi)]. \widehat\psi_{AIPW}=\frac{1}{n}\sum_i\left[ \widehat m_1(W_i)-\widehat m_0(W_i)+ \frac{A_i\{Y_i-\widehat m_1(W_i)\}}{\widehat e(W_i)}- \frac{(1-A_i)\{Y_i-\widehat m_0(W_i)\}}{1-\widehat e(W_i)} \right]. 

outcome_fit <- glm(
  Y ~ A * (W1 + W2 + W3),
  family = binomial(), data = ps_dat
)
data_a1 <- transform(ps_dat, A = 1)
data_a0 <- transform(ps_dat, A = 0)
m1_hat <- predict(outcome_fit, newdata = data_a1, type = "response")
m0_hat <- predict(outcome_fit, newdata = data_a0, type = "response")
e_hat <- ps_dat$e_hat

augmented_A1 <- m1_hat + ps_dat$A * (ps_dat$Y - m1_hat) / e_hat
augmented_A0 <- m0_hat + (1 - ps_dat$A) * (ps_dat$Y - m0_hat) / (1 - e_hat)
aipw_contribution <- augmented_A1 - augmented_A0
aipw_ate <- mean(aipw_contribution)
aipw_se <- sd(aipw_contribution - aipw_ate) / sqrt(nrow(ps_dat))

data.frame(
  estimate = aipw_ate,
  standard_error = aipw_se,
  lower = aipw_ate - qnorm(0.975) * aipw_se,
  upper = aipw_ate + qnorm(0.975) * aipw_se,
  finite_sample_true_ATE = finite_sample_true_ate
)

The interval above uses the empirical efficient-influence-function contribution and is illustrative when both working nuisance models have adequate probability limits and regularity conditions hold. Double robustness of the point estimator’s consistency does not make this confidence interval doubly robust when one nuisance model is misspecified; use inference justified for the fitted estimating system, such as an appropriate stacked sandwich or a bootstrap that refits both models.

Under regularity and causal-identification assumptions, AIPW is consistent if either the propensity model or the outcome model is correctly specified. This is not protection when both are wrong; it does not solve unmeasured confounding, positivity violations, measurement error, or interference. It means consistency, not exact finite-sample unbiasedness. With adaptive machine learning, sample splitting/cross-fitting and appropriate inference are normally needed; this parametric teaching example does not cover those tools.

aipw_components <- function(data, ps_formula, outcome_formula) {
  e_fit <- glm(ps_formula, family = binomial(), data = data)
  e <- validate_probability(predict(e_fit, type = "response"),
                            "AIPW propensity score")
  m_fit <- glm(outcome_formula, family = binomial(), data = data)
  d1 <- transform(data, A = 1)
  d0 <- transform(data, A = 0)
  m1 <- predict(m_fit, newdata = d1, type = "response")
  m0 <- predict(m_fit, newdata = d0, type = "response")
  c(
    standardization = mean(m1 - m0),
    IPTW = mean(data$A * data$Y / e -
                  (1 - data$A) * data$Y / (1 - e)),
    AIPW = mean(m1 - m0 + data$A * (data$Y - m1) / e -
                  (1 - data$A) * (data$Y - m0) / (1 - e))
  )
}

correct_ps <- A ~ W1 + W2 + W3 + W1:W2
wrong_ps <- A ~ W2
correct_outcome <- Y ~ A * (W1 + W2 + W3)
wrong_outcome <- Y ~ A + W2

dr_table <- rbind(
  `both working models adequate` =
    aipw_components(ps_dat, correct_ps, correct_outcome),
  `propensity adequate; outcome misspecified` =
    aipw_components(ps_dat, correct_ps, wrong_outcome),
  `propensity misspecified; outcome adequate` =
    aipw_components(ps_dat, wrong_ps, correct_outcome),
  `both misspecified` =
    aipw_components(ps_dat, wrong_ps, wrong_outcome)
)
cbind(dr_table, finite_sample_true_ATE = finite_sample_true_ate)
##                                           standardization    IPTW    AIPW
## both working models adequate                      0.09274 0.09639 0.09552
## propensity adequate; outcome misspecified         0.16861 0.09639 0.09210
## propensity misspecified; outcome adequate         0.09274 0.16862 0.09274
## both misspecified                                 0.16861 0.16862 0.16862
##                                           finite_sample_true_ATE
## both working models adequate                              0.1085
## propensity adequate; outcome misspecified                 0.1085
## propensity misspecified; outcome adequate                 0.1085
## both misspecified                                         0.1085

This single simulated sample illustrates, rather than proves, the large-sample property. Diagnostics still include overlap, balance for the treatment model, calibration and functional-form checks for the outcome model, influential contributions, and sensitivity to plausible nuisance specifications. The classic reference is Bang and Robins (2005).

7. Treatment–confounder feedback and marginal structural models

Suppose baseline treatment A0A_0 affects a later health measure L1L_1; L1L_1 then affects treatment A1A_1 and outcome YY. L1L_1 is simultaneously a mediator of A0A_0 and a confounder of A1→YA_1\rightarrow Y. Ordinary regression adjustment for L1L_1 can block part of the earlier treatment effect and induce bias; omitting it leaves later treatment confounded. Merely adding L1L_1 to an ordinary time-dependent Cox or regression model does not solve this feedback because it still conditions on a variable affected by prior treatment.

For a sustained-strategy contrast E(YA0=1,A1=1)−E(YA0=0,A1=0)E(Y^{A_0=1,A_1=1})-E(Y^{A_0=0,A_1=0}), stabilized treatment weights are products of time-specific probabilities:

SWi=P(A0i)P(A0i∣Wi)×P(A1i∣A0i)P(A1i∣A0i,Wi,L1i). SW_i=\frac{P(A_{0i})}{P(A_{0i}\mid W_i)} \times \frac{P(A_{1i}\mid A_{0i})} {P(A_{1i}\mid A_{0i},W_i,L_{1i})}. 

set.seed(1071)
n_long <- 5000
long_dat <- data.frame(W = rnorm(n_long))
long_dat$pA0 <- plogis(-0.10 + 0.70 * long_dat$W)
long_dat$A0 <- rbinom(n_long, 1, long_dat$pA0)
long_dat$pL1 <- plogis(-0.20 + 0.80 * long_dat$W + 1.00 * long_dat$A0)
long_dat$L1 <- rbinom(n_long, 1, long_dat$pL1)
long_dat$pA1 <- plogis(-0.40 + 0.50 * long_dat$W + 1.10 * long_dat$L1 +
                         0.60 * long_dat$A0)
long_dat$A1 <- rbinom(n_long, 1, long_dat$pA1)
long_dat$pY <- plogis(-2.00 + 0.45 * long_dat$A0 + 0.70 * long_dat$A1 +
                        0.90 * long_dat$L1 + 0.35 * long_dat$W -
                        0.20 * long_dat$A0 * long_dat$A1)
long_dat$Y <- rbinom(n_long, 1, long_dat$pY)

# The structural equations let us calculate the finite-sample intervention truth.
true_regime_risk <- function(w, a0, a1) {
  p_l <- plogis(-0.20 + 0.80 * w + 1.00 * a0)
  p_y_l0 <- plogis(-2.00 + 0.45 * a0 + 0.70 * a1 + 0.35 * w -
                       0.20 * a0 * a1)
  p_y_l1 <- plogis(-2.00 + 0.45 * a0 + 0.70 * a1 + 0.90 + 0.35 * w -
                       0.20 * a0 * a1)
  mean((1 - p_l) * p_y_l0 + p_l * p_y_l1)
}
true_sustained_rd <- true_regime_risk(long_dat$W, 1, 1) -
  true_regime_risk(long_dat$W, 0, 0)
c(true_risk_always = true_regime_risk(long_dat$W, 1, 1),
  true_risk_never = true_regime_risk(long_dat$W, 0, 0),
  true_risk_difference = true_sustained_rd)
##     true_risk_always      true_risk_never true_risk_difference
##               0.4023               0.1899               0.2124

7.1 Fit, diagnose, and bootstrap an MSM

num0_fit <- glm(A0 ~ 1, family = binomial(), data = long_dat)
den0_fit <- glm(A0 ~ W, family = binomial(), data = long_dat)
num1_fit <- glm(A1 ~ A0, family = binomial(), data = long_dat)
den1_fit <- glm(A1 ~ A0 + W + L1, family = binomial(), data = long_dat)

p_num0 <- predict(num0_fit, type = "response")
p_den0 <- validate_probability(predict(den0_fit, type = "response"),
                               "baseline treatment probability")
p_num1 <- predict(num1_fit, type = "response")
p_den1 <- validate_probability(predict(den1_fit, type = "response"),
                               "follow-up treatment probability")
long_dat$sw0 <- observed_probability(p_num0, long_dat$A0) /
  observed_probability(p_den0, long_dat$A0)
long_dat$sw <- long_dat$sw0 *
  observed_probability(p_num1, long_dat$A1) /
  observed_probability(p_den1, long_dat$A1)

rbind(baseline_weight = weight_summary(long_dat$sw0),
      cumulative_weight = weight_summary(long_dat$sw))
##                    mean     sd    min    q50  q95   q99    max  ESS
## baseline_weight   1.001 0.3716 0.5085 0.9048 1.69 2.362  5.024 4394
## cumulative_weight 1.002 0.6679 0.3231 0.8053 1.92 3.564 12.849 3461
data.frame(
  diagnostic = c("W balance at A0",
                 "W balance at A1 within A0=0",
                 "L1 balance at A1 within A0=0",
                 "W balance at A1 within A0=1",
                 "L1 balance at A1 within A0=1"),
  unweighted_SMD = c(
    weighted_smd(long_dat$W, long_dat$A0),
    with(subset(long_dat, A0 == 0), weighted_smd(W, A1)),
    with(subset(long_dat, A0 == 0), weighted_smd(L1, A1)),
    with(subset(long_dat, A0 == 1), weighted_smd(W, A1)),
    with(subset(long_dat, A0 == 1), weighted_smd(L1, A1))
  ),
  weighted_SMD = c(
    weighted_smd(long_dat$W, long_dat$A0, long_dat$sw0),
    with(subset(long_dat, A0 == 0), weighted_smd(W, A1, sw)),
    with(subset(long_dat, A0 == 0), weighted_smd(L1, A1, sw)),
    with(subset(long_dat, A0 == 1), weighted_smd(W, A1, sw)),
    with(subset(long_dat, A0 == 1), weighted_smd(L1, A1, sw))
  )
)
msm_fit <- glm(Y ~ factor(A0) * factor(A1),
               family = quasibinomial(), weights = sw, data = long_dat)
regimes <- data.frame(A0 = c(0, 1), A1 = c(0, 1))
msm_risks <- predict(msm_fit, newdata = regimes, type = "response")
c(risk_never = msm_risks[1], risk_always = msm_risks[2],
  risk_difference = msm_risks[2] - msm_risks[1],
  finite_sample_true_RD = true_sustained_rd)
##          risk_never.1         risk_always.2     risk_difference.2 finite_sample_true_RD
##                0.1815                0.4062                0.2246                0.2124

Balance should be assessed at each treatment time using the history relevant to that decision, not only at baseline. Also inspect treatment probabilities within clinically important histories, cumulative-weight tails, ESS over time, and the number following each regime.

msm_sustained_rd <- function(data) {
  n0 <- glm(A0 ~ 1, family = binomial(), data = data)
  d0 <- glm(A0 ~ W, family = binomial(), data = data)
  n1 <- glm(A1 ~ A0, family = binomial(), data = data)
  d1 <- glm(A1 ~ A0 + W + L1, family = binomial(), data = data)
  sw0 <- observed_probability(predict(n0, type = "response"), data$A0) /
    observed_probability(validate_probability(predict(d0, type = "response"),
                                               "bootstrap baseline treatment probability"),
                         data$A0)
  sw <- sw0 * observed_probability(predict(n1, type = "response"), data$A1) /
    observed_probability(validate_probability(predict(d1, type = "response"),
                                               "bootstrap follow-up treatment probability"),
                         data$A1)
  fit <- glm(Y ~ factor(A0) * factor(A1), family = quasibinomial(),
             weights = sw, data = data)
  pred <- predict(fit,
                  newdata = data.frame(A0 = c(0, 1), A1 = c(0, 1)),
                  type = "response")
  unname(pred[2] - pred[1])
}

set.seed(1072)
msm_boot <- replicate(200, {
  index <- sample.int(nrow(long_dat), replace = TRUE)
  msm_sustained_rd(long_dat[index, ])
})
data.frame(
  estimate = msm_sustained_rd(long_dat),
  bootstrap_se = sd(msm_boot),
  lower = unname(quantile(msm_boot, 0.025)),
  upper = unname(quantile(msm_boot, 0.975))
)

This resamples people and rebuilds every weight model. With repeated records, resample at the independent-person or independent-cluster level. Weight truncation, alternate models, and alternate treatment definitions are sensitivity analyses, not substitutes for sequential exchangeability. See Robins, Hernán, and Brumback, Cole and Hernán, and the parametric g-formula worked example.

8. IPCW and a limited base-R multiple-imputation example

Missing outcomes and loss to follow-up are observation processes, not merely software inconveniences. Complete-case analysis targets the original population only under restrictive conditions. First distinguish:

  • the scientific estimand under complete follow-up;
  • why each value is unobserved and which predictors of observation are available;
  • the assumptions made by the primary method; and
  • sensitivity analyses for plausible departures.

8.1 Inverse probability of censoring weighting

Let R=1R=1 mean that the endpoint is observed. A stabilized observation weight is

SWiC=P(Ri=1∣Ai)P(Ri=1∣Ai,Wi)among people with Ri=1. SW_i^C=\frac{P(R_i=1\mid A_i)} {P(R_i=1\mid A_i,W_i)} \quad\text{among people with }R_i=1.

IPCW identifies the complete-follow-up contrast if observation is conditionally exchangeable given the modeled history, observation positivity holds, and the models and measurements are adequate.

set.seed(1081)
n_cens <- 2800
cens_dat <- data.frame(
  A = rbinom(n_cens, 1, 0.5),
  W = rnorm(n_cens),
  Q = rbinom(n_cens, 1, 0.4)
)
cens_dat$pY <- plogis(-1.35 + 0.55 * cens_dat$A + 0.75 * cens_dat$W +
                        0.45 * cens_dat$Q)
cens_dat$Y_full <- rbinom(n_cens, 1, cens_dat$pY)
cens_dat$pR <- plogis(1.25 - 0.50 * cens_dat$A - 0.80 * cens_dat$W +
                        0.40 * cens_dat$Q)
cens_dat$R <- rbinom(n_cens, 1, cens_dat$pR)
cens_dat$Y <- ifelse(cens_dat$R == 1, cens_dat$Y_full, NA)

num_cens <- glm(R ~ A, family = binomial(), data = cens_dat)
den_cens <- glm(R ~ A + W + Q, family = binomial(), data = cens_dat)
p_num <- predict(num_cens, type = "response")
p_den <- validate_probability(predict(den_cens, type = "response"),
                              "observation probability")
cens_dat$cw <- p_num / p_den

observed <- cens_dat$R == 1
naive_rd <- with(cens_dat[observed, ], mean(Y[A == 1]) - mean(Y[A == 0]))
ipcw_rd <- with(cens_dat[observed, ],
                weighted_mean(Y[A == 1], cw[A == 1]) -
                  weighted_mean(Y[A == 0], cw[A == 0]))
full_data_rd <- with(cens_dat,
                     mean(Y_full[A == 1]) - mean(Y_full[A == 0]))

data.frame(
  observed_fraction = mean(cens_dat$R),
  complete_case_RD = naive_rd,
  IPCW_RD = ipcw_rd,
  randomized_full_data_RD = full_data_rd,
  expected_causal_RD = mean(plogis(-1.35 + 0.55 + 0.75 * cens_dat$W +
                                      0.45 * cens_dat$Q) -
                              plogis(-1.35 + 0.75 * cens_dat$W +
                                      0.45 * cens_dat$Q))
)
rbind(observation_weights = weight_summary(cens_dat$cw[observed]))
##                       mean     sd    min    q50   q95   q99   max  ESS
## observation_weights 0.9967 0.2413 0.7061 0.9305 1.421 1.978 3.268 1944

For time-varying loss to follow-up, multiply interval-specific conditional probabilities up to each time, respecting the order “history →\rightarrow censoring decision →\rightarrow next outcome.” Combine treatment and censoring weights only when both processes require weighting, then diagnose each component and the product. The classic dependent-censoring reference is Robins and Finkelstein.

ipcw_risk_difference <- function(data) {
  num <- glm(R ~ A, family = binomial(), data = data)
  den <- glm(R ~ A + W + Q, family = binomial(), data = data)
  cw <- predict(num, type = "response") /
    validate_probability(predict(den, type = "response"),
                         "bootstrap observation probability")
  keep <- data$R == 1
  d <- data[keep, ]
  w <- cw[keep]
  weighted_mean(d$Y[d$A == 1], w[d$A == 1]) -
    weighted_mean(d$Y[d$A == 0], w[d$A == 0])
}

set.seed(1082)
ipcw_boot <- replicate(300, {
  index <- sample.int(nrow(cens_dat), replace = TRUE)
  ipcw_risk_difference(cens_dat[index, ])
})
data.frame(
  estimate = ipcw_risk_difference(cens_dat),
  bootstrap_se = sd(ipcw_boot),
  lower = unname(quantile(ipcw_boot, 0.025)),
  upper = unname(quantile(ipcw_boot, 0.975))
)

8.2 Proper multiple imputation: a transparent narrow example

The following code demonstrates Bayesian linear-regression imputation for a continuous, approximately normal outcome missing under MAR conditional on fully observed predictors. Each imputation draws the residual variance, draws regression coefficients, and draws missing outcomes from their posterior predictive distribution. Rubin’s rules combine within- and between-imputation variance.

The pooling function below uses Rubin’s original large-sample degrees-of-freedom approximation. That is reasonable for this n=900n=900 teaching example; small-sample work should use a finite-complete-data adjustment and validated MI software.

It is intentionally not a general imputation package: it does not handle categorical, bounded, multilevel, longitudinal, or survival data; passive variables; complex interactions; survey designs; or perfect prediction. Do not use mean imputation, treat one completed data set as observed, or round imputed binary values.

set.seed(1083)
n_mi <- 900
mi_dat <- data.frame(
  A = rbinom(n_mi, 1, 0.5),
  X = rnorm(n_mi),
  Z = rnorm(n_mi)
)
mi_dat$Y_full <- 1.0 + 0.75 * mi_dat$A + 0.90 * mi_dat$X -
  0.45 * mi_dat$Z + rnorm(n_mi, sd = 1.1)
mi_dat$pR <- plogis(1.05 - 0.90 * mi_dat$A + 0.55 * mi_dat$X -
                      0.35 * mi_dat$Z)
mi_dat$R <- rbinom(n_mi, 1, mi_dat$pR)
mi_dat$Y <- ifelse(mi_dat$R == 1, mi_dat$Y_full, NA)

c(observed_fraction = mean(mi_dat$R),
  missing_A0 = mean(is.na(mi_dat$Y[mi_dat$A == 0])),
  missing_A1 = mean(is.na(mi_dat$Y[mi_dat$A == 1])))
## observed_fraction        missing_A0        missing_A1
##            0.6422            0.2552            0.4739
pool_rubin <- function(estimates, variances) {
  m <- nrow(estimates)
  q_bar <- colMeans(estimates)
  u_bar <- colMeans(variances)
  b <- apply(estimates, 2, var)
  total <- u_bar + (1 + 1 / m) * b
  df <- ifelse(b > 0,
               (m - 1) * (1 + u_bar / ((1 + 1 / m) * b))^2,
               Inf)
  critical <- qt(0.975, df = df)
  data.frame(
    term = colnames(estimates), estimate = q_bar,
    standard_error = sqrt(total), df = df,
    lower = q_bar - critical * sqrt(total),
    upper = q_bar + critical * sqrt(total),
    row.names = NULL
  )
}

mi_normal_outcome <- function(data, m = 40, delta = 0, seed = 1) {
  stopifnot(m >= 2, length(delta) == 1, is.finite(delta))
  set.seed(seed)
  xmat <- model.matrix(~ A + X + Z, data = data)
  observed <- !is.na(data$Y)
  missing <- !observed
  x_obs <- xmat[observed, , drop = FALSE]
  y_obs <- data$Y[observed]
  fit <- lm.fit(x = x_obs, y = y_obs)
  beta_hat <- fit$coefficients
  residual_df <- length(y_obs) - ncol(x_obs)
  s2_hat <- sum(fit$residuals^2) / residual_df
  xtx_inverse <- solve(crossprod(x_obs))

  estimates <- matrix(NA_real_, m, ncol(xmat),
                      dimnames = list(NULL, colnames(xmat)))
  variances <- estimates
  for (j in seq_len(m)) {
    sigma2_draw <- residual_df * s2_hat / rchisq(1, residual_df)
    beta_cov <- sigma2_draw * xtx_inverse
    beta_draw <- beta_hat +
      drop(t(chol(beta_cov)) %*% rnorm(length(beta_hat)))
    completed_y <- data$Y
    completed_y[missing] <- rnorm(
      sum(missing),
      mean = drop(xmat[missing, , drop = FALSE] %*% beta_draw) + delta,
      sd = sqrt(sigma2_draw)
    )
    completed <- transform(data, Y = completed_y)
    analysis <- lm(Y ~ A + X + Z, data = completed)
    estimates[j, ] <- coef(analysis)
    variances[j, ] <- diag(vcov(analysis))
  }
  list(pooled = pool_rubin(estimates, variances),
       estimates = estimates, variances = variances)
}

mi_primary <- mi_normal_outcome(mi_dat, m = 40, delta = 0, seed = 1084)
complete_case_fit <- lm(Y ~ A + X + Z, data = mi_dat)
full_data_fit <- lm(Y_full ~ A + X + Z, data = mi_dat)

rbind(
  complete_case = c(estimate = coef(complete_case_fit)["A"],
                    standard_error = sqrt(vcov(complete_case_fit)["A", "A"])),
  multiple_imputation = c(
    estimate = subset(mi_primary$pooled, term == "A")$estimate,
    standard_error = subset(mi_primary$pooled, term == "A")$standard_error
  ),
  unavailable_full_data_benchmark = c(
    estimate = coef(full_data_fit)["A"],
    standard_error = sqrt(vcov(full_data_fit)["A", "A"])
  )
)
##                                 estimate.A standard_error
## complete_case                       0.9481        0.09344
## multiple_imputation                 0.9547        0.09947
## unavailable_full_data_benchmark     0.8990        0.07215

MAR is a conditional, untestable assumption: after conditioning on the imputation variables, missingness does not depend on the missing value. The imputation model should include every analysis variable, the outcome when imputing covariates, auxiliary predictors of missingness or values, and functional forms/interactions needed by the analysis. The number of imputations should reflect the fraction of missing information and Monte Carlo error, not a ritual fixed number.

8.3 Delta-adjusted sensitivity analysis

A pattern-mixture sensitivity analysis shifts each imputed missing outcome by δ\delta relative to its MAR prediction. Negative δ\delta means missing outcomes are systematically lower, after conditioning on the imputation predictors.

delta_grid <- seq(-0.75, 0.75, by = 0.25)
delta_results <- do.call(rbind, lapply(delta_grid, function(delta) {
  result <- mi_normal_outcome(mi_dat, m = 40, delta = delta, seed = 1085)
  treatment <- subset(result$pooled, term == "A")
  data.frame(delta = delta, estimate = treatment$estimate,
             lower = treatment$lower, upper = treatment$upper)
}))
delta_results
plot(delta_results$delta, delta_results$estimate, type = "b", pch = 19,
     xlab = expression(delta~"shift for missing outcomes"),
     ylab = "Adjusted mean difference for A")
segments(delta_results$delta, delta_results$lower,
         delta_results$delta, delta_results$upper, col = "grey40")
abline(h = 0, lty = 2)

Choose delta values using outcome units, validation data, clinical knowledge, or a tipping-point question. Reusing the same random seed here isolates the delta shift from Monte Carlo noise. A single common delta is simplistic; arm-specific or covariate-dependent shifts may be more credible. See White, Royston, and Wood for practical multiple-imputation guidance.

9. Transportability weights

Internal validity does not guarantee relevance to a new population. Let S=1S=1 denote trial/study participation and S=0S=0 a representative target-population sample. For a trial with randomized AA, inverse odds of sampling weights

wS(W)=1−P(S=1∣W)P(S=1∣W) w^S(W)=\frac{1-P(S=1\mid W)}{P(S=1\mid W)}

reweight participants toward the target covariate distribution. The target average treatment effect additionally requires consistency, trial internal validity, exchangeability of potential outcomes over SS conditional on effect modifiers, selection positivity, harmonized measurements, and a representative target sample (or its design weights).

Selection-odds weights alone suffice for treatment assignment in this example because the trial uses unconditional 1:1 randomization. Unequal or covariate- adaptive randomization requires the known treatment probabilities in the estimating procedure; observational treatment additionally requires confounding adjustment.

set.seed(1091)
n_source <- 1400
n_target <- 2600
source <- data.frame(
  S = 1L,
  W = rnorm(n_source, mean = -0.30),
  Z = rbinom(n_source, 1, 0.35)
)
target <- data.frame(
  S = 0L,
  W = rnorm(n_target, mean = 0.50),
  Z = rbinom(n_target, 1, 0.65)
)
source$A <- rbinom(n_source, 1, 0.5)
source$p0 <- plogis(-1.40 + 0.50 * source$W + 0.35 * source$Z)
source$p1 <- plogis(-1.40 + 0.55 + 0.50 * source$W + 0.35 * source$Z -
                     0.35 * source$W + 0.30 * source$Z)
source$Y <- rbinom(n_source, 1,
                   ifelse(source$A == 1, source$p1, source$p0))

combined <- rbind(
  source[c("S", "W", "Z")],
  target[c("S", "W", "Z")]
)
selection_fit <- glm(S ~ W + Z, family = binomial(), data = combined)
combined$pS <- validate_probability(predict(selection_fit, type = "response"),
                                    "study-participation probability")
source$selection_odds_weight <-
  (1 - combined$pS[combined$S == 1]) / combined$pS[combined$S == 1]

naive_trial_rd <- with(source, mean(Y[A == 1]) - mean(Y[A == 0]))
transported_rd <- with(source,
  weighted_mean(Y[A == 1], selection_odds_weight[A == 1]) -
    weighted_mean(Y[A == 0], selection_odds_weight[A == 0]))
target_p0 <- plogis(-1.40 + 0.50 * target$W + 0.35 * target$Z)
target_p1 <- plogis(-1.40 + 0.55 + 0.50 * target$W + 0.35 * target$Z -
                      0.35 * target$W + 0.30 * target$Z)

all_weights <- ifelse(combined$S == 1,
                      source$selection_odds_weight, 1)
transport_balance <- data.frame(
  variable = c("W", "Z"),
  source_minus_target_SMD = c(
    weighted_smd(combined$W, combined$S),
    weighted_smd(combined$Z, combined$S)
  ),
  weighted_source_minus_target_SMD = c(
    weighted_smd(combined$W, combined$S, all_weights),
    weighted_smd(combined$Z, combined$S, all_weights)
  )
)

list(
  effects = data.frame(
    naive_trial_RD = naive_trial_rd,
    transported_RD = transported_rd,
    finite_target_true_RD = mean(target_p1 - target_p0)
  ),
  selection_weight_summary = weight_summary(source$selection_odds_weight),
  balance = transport_balance
)
## $effects
##   naive_trial_RD transported_RD finite_target_true_RD
## 1         0.1253         0.1299                0.1198
##
## $selection_weight_summary
##      mean        sd       min       q50       q95       q99       max       ESS
##   1.88616   2.79521   0.03619   0.99281   6.32314  13.50519  36.67638 438.23468
##
## $balance
##   variable source_minus_target_SMD weighted_source_minus_target_SMD
## 1        W                 -0.8508                          0.01922
## 2        Z                 -0.6665                          0.03379

The arbitrary relative sample sizes alter selection-model odds by a constant, but that constant cancels in normalized arm-specific means. With nonrandom treatment, transport and confounding adjustment must both be addressed. Inspect selection score overlap, effect-modifier balance, weights, ESS, treatment versions, outcome definitions, and health-system differences. Unmeasured effect modification can defeat transport even when measured balance is excellent. See Lesko et al. and the 2026 ISPE-endorsed framework.

10. Case-crossover studies

A case-crossover design compares a person’s exposure during a hazard window just before an acute event with exposure during referent windows for the same person. Time-invariant characteristics cancel through within-person conditioning. The design is best suited to a transient exposure, an abrupt outcome, a short and credible induction period, and no exposure carryover.

set.seed(1101)
n_cases <- 800
n_periods <- 4
cc_dat <- expand.grid(id = seq_len(n_cases), period = seq_len(n_periods))
person_tendency <- rnorm(n_cases)
cc_dat$temperature <- rnorm(nrow(cc_dat)) + 0.20 * sin(cc_dat$period)
cc_dat$exposure <- rbinom(
  nrow(cc_dat), 1,
  plogis(-0.55 + person_tendency[cc_dat$id] + 0.45 * cc_dat$temperature)
)

# Exactly one event window per person, sampled according to a conditional-logit DGP.
true_log_or <- log(1.8)
cc_dat$case_window <- 0L
for (i in seq_len(n_cases)) {
  rows <- which(cc_dat$id == i)
  score <- exp(true_log_or * cc_dat$exposure[rows] +
                 0.35 * cc_dat$temperature[rows])
  chosen <- sample(rows, size = 1, prob = score)
  cc_dat$case_window[chosen] <- 1L
}

if (!has_survival) {
  cat("The case-crossover model was skipped because 'survival' is unavailable.\n")
} else {
  # In survival 3.8.x, attaching the package ensures strata() is recognized as a
  # conditional-likelihood special inside clogit().
  suppressPackageStartupMessages(library(survival))
  cc_unadjusted <- clogit(case_window ~ exposure + strata(id),
                          data = cc_dat, method = "exact")
  cc_adjusted <- clogit(case_window ~ exposure + temperature + strata(id),
                        data = cc_dat, method = "exact")
  data.frame(
    model = c("Exposure only", "Adjusted for time-varying temperature"),
    exposure_OR = exp(c(coef(cc_unadjusted)["exposure"],
                        coef(cc_adjusted)["exposure"])),
    true_conditional_OR = exp(true_log_or)
  )
}

Conditional logistic regression respects the matched risk sets. With valid referent sampling and the design assumptions, its exposure coefficient estimates a within-person incidence rate ratio. It is not a fixed-horizon risk ratio.

Major threats include time trends, seasonality, overlap or carryover between windows, exposure changes caused by prodromal symptoms, time-varying confounders, referent-window selection, and event-dependent exposure opportunity. A time-stratified referent strategy can prevent overlap bias in environmental applications. The foundational paper is Maclure (1991).

11. Misclassification bias analysis and negative controls

Saying that bias is “probably toward the null” is not an analysis. Direction depends on what is misclassified, whether errors differ by exposure/outcome, the measure and adjustment structure, and other biases.

11.1 Deterministic correction for a binary outcome

Within exposure stratum aa, observed outcome risk pa*p_a^* relates to true risk pap_a through sensitivity SeaSe_a and specificity SpaSp_a:

pa*=Seapa+(1−Spa)(1−pa),pa=pa*+Spa−1Sea+Spa−1. p_a^*=Se_a p_a+(1-Sp_a)(1-p_a),\qquad p_a=\frac{p_a^*+Sp_a-1}{Se_a+Sp_a-1}.

The correction requires Sea+Spa>1Se_a+Sp_a>1 and must yield a probability in [0,1][0,1].

correct_binary_outcome_risk <- function(observed_cases, total,
                                        sensitivity, specificity) {
  stopifnot(total > 0, observed_cases >= 0, observed_cases <= total,
            sensitivity > 0, sensitivity <= 1,
            specificity > 0, specificity <= 1,
            sensitivity + specificity > 1)
  p_observed <- observed_cases / total
  p_corrected <- (p_observed + specificity - 1) /
    (sensitivity + specificity - 1)
  if (p_corrected < 0 || p_corrected > 1) return(NA_real_)
  p_corrected
}

qba_inputs <- data.frame(
  A = c(1, 0), cases_observed = c(144, 120), total = c(800, 1000),
  sensitivity = c(0.85, 0.78), specificity = c(0.98, 0.99)
)
qba_inputs$risk_observed <- with(qba_inputs, cases_observed / total)
qba_inputs$risk_corrected <- mapply(
  correct_binary_outcome_risk,
  qba_inputs$cases_observed, qba_inputs$total,
  qba_inputs$sensitivity, qba_inputs$specificity
)

qba_inputs
data.frame(
  measure = c("Risk difference", "Risk ratio"),
  observed = c(diff(rev(qba_inputs$risk_observed)),
               qba_inputs$risk_observed[1] / qba_inputs$risk_observed[2]),
  corrected = c(diff(rev(qba_inputs$risk_corrected)),
                qba_inputs$risk_corrected[1] / qba_inputs$risk_corrected[2])
)

The inputs are assumptions, not facts. Justify them from an internal validation study when possible, or assess transportability from external validation data. This example allows differential outcome misclassification because sensitivity and specificity differ by exposure stratum.

11.2 Probabilistic bias analysis

Probabilistic bias analysis replaces each bias parameter with a distribution. The example also draws observed risks from Jeffreys beta distributions to reflect binomial sampling uncertainty. The output is a simulation interval under the specified distributions—not automatically a frequentist 95% confidence interval and not an objective posterior distribution.

set.seed(1111)
bias_draws <- 10000
p_obs_1 <- rbeta(bias_draws, 144 + 0.5, 800 - 144 + 0.5)
p_obs_0 <- rbeta(bias_draws, 120 + 0.5, 1000 - 120 + 0.5)

# Illustrative beta distributions; real parameters require documented evidence.
se_1 <- rbeta(bias_draws, 85, 15)
sp_1 <- rbeta(bias_draws, 98, 2)
se_0 <- rbeta(bias_draws, 78, 22)
sp_0 <- rbeta(bias_draws, 99, 1)

p_true_1 <- (p_obs_1 + sp_1 - 1) / (se_1 + sp_1 - 1)
p_true_0 <- (p_obs_0 + sp_0 - 1) / (se_0 + sp_0 - 1)
valid <- (se_1 + sp_1 > 1) & (se_0 + sp_0 > 1) &
  is.finite(p_true_1) & is.finite(p_true_0) &
  p_true_1 > 0 & p_true_1 < 1 & p_true_0 > 0 & p_true_0 < 1
rd_draw <- p_true_1[valid] - p_true_0[valid]
rr_draw <- p_true_1[valid] / p_true_0[valid]

data.frame(
  quantity = c("Corrected risk difference", "Corrected risk ratio"),
  median = c(median(rd_draw), median(rr_draw)),
  lower_2.5 = c(quantile(rd_draw, 0.025), quantile(rr_draw, 0.025)),
  upper_97.5 = c(quantile(rd_draw, 0.975), quantile(rr_draw, 0.975)),
  valid_draw_fraction = mean(valid)
)

Correlations among sensitivities/specificities, uncertainty in validation data, and multiple simultaneous biases often matter. Report parameter sources, distributions, dependencies, invalid draws, algorithms, random seed, and results across scenarios. A substantial invalid-draw fraction means the assumed bias- parameter distributions are incompatible with the observed data; discarding those draws conditions the assumed joint distribution and is not a repair. See the STRATOS measurement-error guidance, Fox, MacLehose, and Lash (2021), and the BMJ quantitative-bias-analysis overview.

11.3 Negative controls are bias probes

A negative-control outcome should not plausibly be caused by the exposure but should share relevant confounding, selection, or measurement pathways. A negative-control exposure should not plausibly cause the outcome in the specified window but should share those bias pathways.

set.seed(1112)
n_nc <- 4500
nc_dat <- data.frame(
  W = rnorm(n_nc),
  U = rnorm(n_nc)
)
nc_dat$A <- rbinom(n_nc, 1,
                   plogis(-0.25 + 0.55 * nc_dat$W + 0.85 * nc_dat$U))
nc_dat$Y <- rbinom(n_nc, 1,
                   plogis(-1.45 + 0.45 * nc_dat$A + 0.45 * nc_dat$W +
                            0.75 * nc_dat$U))
# By construction A has no causal arrow to Y_negative.
nc_dat$Y_negative <- rbinom(n_nc, 1,
                            plogis(-1.35 + 0.40 * nc_dat$W + 0.80 * nc_dat$U))
# Future/proxy exposure shares U and W but has no causal arrow to current Y.
nc_dat$A_negative <- rbinom(n_nc, 1,
                            plogis(-0.15 + 0.50 * nc_dat$W + 0.80 * nc_dat$U))

nc_models <- list(
  `Primary A -> Y` = glm(Y ~ A + W, family = binomial(), data = nc_dat),
  `A -> negative-control outcome` =
    glm(Y_negative ~ A + W, family = binomial(), data = nc_dat),
  `Negative-control exposure -> Y` =
    glm(Y ~ A_negative + W, family = binomial(), data = nc_dat)
)
data.frame(
  contrast = names(nc_models),
  odds_ratio = exp(vapply(nc_models, function(fit) coef(fit)[2], numeric(1))),
  row.names = NULL
)

Here, associations with both negative controls flag residual bias from unmeasured UU. In real data, a non-null control can also reflect chance or failure of the “cannot cause” and shared-bias assumptions. A null control does not prove absence of bias and can have low power or weak connection to the primary bias pathway. Do not mechanically subtract a negative-control association without an identified method. See Lipsitch, Tchetgen Tchetgen, and Cohen.

12. Synthesis, reporting, exercises, and reference tools

Advanced methods are most defensible when they are treated as linked parts of a design rather than a menu of sophisticated models.

12.1 Method-to-question map

Question/data problem Example target quantity Candidate estimator Nonnegotiable diagnostics
Right-censored time to one event S(τ)S(\tau), risk, or RMST Kaplan–Meier; RMST integration Follow-up support, censoring patterns, risk sets
Covariate effects on hazard Conditional hazard ratio Cox model Schoenfeld residuals, functional form, influence, events
Event with competing causes Cause-specific cumulative incidence Aalen–Johansen State definitions, all event types, follow-up support
Baseline confounding ATE/ATT/ATO risk contrast IPTW, standardization, AIPW Overlap, weights, ESS, balance, model checks
Time-varying confounding affected by prior treatment Sustained or dynamic regime contrast MSM/IPTW or g-formula Sequential balance/support, cumulative weights, regime counts
Informative loss to follow-up Complete-follow-up contrast IPCW, MI under aligned assumptions Observation model, weights, missingness patterns, sensitivity
Trial/study differs from target TATE Selection-odds weighting or target standardization Effect-modifier overlap/balance, target sampling, harmonization
Transient exposure and abrupt event Acute conditional rate ratio Case-crossover conditional logistic model Referent strategy, time trends, carryover, time-varying confounding
Outcome misclassification Bias-adjusted risk contrast Deterministic/probabilistic QBA Bias-parameter evidence, valid draws, scenario dependence

No row is automatic. For example, a Cox model may describe hazard associations while RMST answers the clinical decision question; MI and IPCW may rely on different representations of the same missingness assumption; and AIPW may be less stable than a simpler estimator under severe positivity problems.

12.2 A reproducible advanced-analysis workflow

  1. Freeze the estimand card. Name the target population, strategies, time zero, horizon, outcome, intercurrent events, effect scale, and averaging population before fitting models.
  2. Draw the longitudinal data structure. Mark when every covariate, treatment, censoring event, competing event, and outcome is measured.
  3. Write the identification argument. State consistency, exchangeability, positivity, interference, censoring, transport, and measurement assumptions separately. “Adjusted analysis” is not an argument.
  4. Lock the analytic cohort. Reproduce eligibility, exclusions, time zero, follow-up, missingness, events, and treatment histories in a flow table.
  5. Run design diagnostics before interpreting effects. Examine overlap, balance, risk sets, weights, ESS, regime support, selection support, and missingness by key variables.
  6. Estimate on interpretable scales. Prefer marginal risks, risk differences, risk ratios, cumulative incidence, or RMST alongside any hazard or odds ratio.
  7. Propagate all important uncertainty. Refit nuisance, calibration, imputation, and weight models inside the bootstrap when using bootstrap inference; resample the independent unit.
  8. Stress-test assumptions. Prespecify alternate functional forms, truncation rules, horizons, censoring models, delta values, bias parameters, target definitions, and negative controls.
  9. Interpret the ensemble. Explain why estimates move, which assumptions each analysis changes, and which limitations no analysis resolves.
  10. Archive provenance. Save code, session information, seeds, data dictionaries, model specifications, diagnostics, and a decision log.

12.3 Minimum reporting table

For each primary and sensitivity analysis, report:

Item What readers need
Estimand Population, strategies, outcome, horizon, intercurrent events, contrast
Cohort construction Eligibility, time zero, follow-up, exclusions, missingness, event counts
Identification Adjustment history and explicit causal assumptions
Models Every numerator/denominator, outcome, selection, and imputation model
Diagnostics Overlap/balance, weight summaries and ESS, risk sets, proportional hazards, invalid bias draws
Estimate Marginal arm-specific quantities, contrast, uncertainty interval, units
Sensitivity What assumption changed, plausible range/rationale, resulting estimate
Reproducibility Software/package versions, seed, code/data availability, deviations from plan

The TARGET Statement (2025) is the current reporting guideline for observational studies explicitly emulating eligible target trials. Use CONSORT 2025 and SPIRIT 2025 for randomized-trial reports and protocols. STROBE remains relevant to cohort, case-control, and cross-sectional reporting, but it is a reporting checklist, not a design-quality or risk-of-bias score. Select a guideline by design and use it in addition to, not instead of, the estimand and diagnostic details above.

12.4 Reusable analysis skeleton

Primary question and decision:
Target population and eligibility:
Treatment/exposure strategies and time zero:
Outcome, competing events, and horizon:
Primary marginal estimand and effect scale:

Identification assumptions:
  - consistency/treatment versions:
  - exchangeability adjustment set/history:
  - positivity/support:
  - censoring/missingness:
  - interference:
  - transport/measurement assumptions:

Primary estimator and uncertainty method:
Nuisance models and prespecified functional forms:
Diagnostics and acceptance/escalation rules:
Sensitivity analyses with scientific ranges:
Negative controls or validation data:
Reporting guideline and reproducibility archive:

12.5 Ten exercises with answers

Exercise 1: choose a survival estimand

A treatment has a strong early benefit and a possible late harm. Investigators propose one Cox hazard ratio over five years. Specify two more informative primary summaries and one diagnostic.

Answer

Specify five-year risks (or their difference) and five-year RMSTs (or their difference). Examine the survival curves and scaled Schoenfeld residuals; if a time-varying hazard contrast remains scientifically useful, prespecify its form. The horizon and competing-event handling must be explicit.

Exercise 2: competing risk

In a study of dementia, death without dementia is coded as censored and dementia risk is reported as 1−1-Kaplan–Meier. What is wrong, and what should replace it?

Answer

Death prevents later dementia and is not ordinary independent censoring for the real-world dementia probability. 1−1-Kaplan–Meier generally overestimates cumulative incidence. Use an Aalen–Johansen cumulative-incidence estimator and state whether the target is a total-effect risk with death present or a different hypothetical estimand.

Exercise 3: proportional hazards

cox.zph() gives p=0.30p=0.30 for treatment, but the residual plot has a pronounced curve and only 70 events occurred. May proportional hazards be declared true?

Answer

No. Failure to reject is not evidence that proportional hazards is true; the test may have low power. Consider clinical expectations, the plot, flexible or prespecified time interactions, and marginal survival/RMST summaries.

Exercise 4: extreme propensity weights

ATE weights have a 99th percentile of 12, maximum of 180, and ESS of 240 from 4,000 people. Balance is good. Is the analysis secure?

Answer

No. Good measured balance does not remove practical positivity and instability. Locate unsupported histories, verify treatment/model coding, inspect influence, and reconsider whether the full-population ATE is identifiable. Report prespecified truncation sensitivity or an overlap-population estimand rather than silently deleting high weights.

Exercise 5: double robustness

An AIPW estimate is labeled “unconfounded because the estimator is doubly robust.” Correct the statement.

Answer

Under consistency, exchangeability, positivity, regularity, and adequate data, AIPW is consistent if either the treatment model or the outcome model is correct. It is not protected when both are wrong and does not address unmeasured confounding, poor measurement, interference, or positivity violations.

Exercise 6: treatment–confounder feedback

Treatment at baseline changes blood pressure at month 3; month-3 blood pressure changes treatment at month 3 and the final outcome. Why can ordinary adjustment be biased?

Answer

Month-3 blood pressure is a confounder of later treatment and outcome but also a mediator of baseline treatment. Conditioning on it can block part of the earlier effect; omitting it leaves later treatment confounded. Under sequential assumptions, use an appropriate g-method such as an MSM with treatment-history weights or the longitudinal g-formula.

Exercise 7: missing outcomes

A researcher imputes every missing continuous outcome with its regression prediction once and runs ordinary linear regression. Identify two errors.

Answer

Single deterministic imputation suppresses residual and imputation uncertainty, so standard errors are too small, and treating completed values as observed is invalid. Proper MI draws parameters and missing values repeatedly and pools within/between-imputation variance. Its MAR/model assumptions still require delta or other MNAR sensitivity analyses.

Exercise 8: transportability

After weighting trial participants to match target age and sex, all SMDs are near zero. Can the target effect be declared unbiased?

Answer

No. Balance covers only measured variables. The analysis also needs internal trial validity, all relevant effect modifiers measured compatibly, selection positivity, representative target data, consistent treatment/outcome versions, and no important engagement or setting effects beyond the modeled variables.

Exercise 9: case-crossover

Daily air pollution rises seasonally and an acute event is compared with control days sampled from the entire year. What threat arises?

Answer

Referent days can differ systematically in season and time trend, producing confounding/overlap bias. A time-stratified referent strategy (for example, same month and weekday) plus measured time-varying confounder adjustment is often more credible. Hazard and carryover windows must also be justified.

Exercise 10: bias analysis and negative controls

A negative-control outcome is null, and the corrected risk ratio remains above 1 throughout the chosen probabilistic-bias-analysis scenarios. Does this establish causality?

Answer

No. A negative control may be weak, underpowered, or fail to share the primary bias pathway. The simulated interval covers only the bias mechanisms and parameter distributions that were specified; it does not address omitted selection, confounding, measurement, model, or positivity problems. Use evidence-based parameters, multiple diagnostics, and substantive triangulation; no single sensitivity analysis proves causality.

Glossary

Term Concise meaning
Aalen–Johansen estimator Product-integral estimator of state probabilities/cumulative incidence
AIPW Outcome regression augmented by inverse treatment weighting
ATE / ATT / ATO Average effect in the full, treated, or overlap population
Censoring exchangeability Conditional independence needed to recover outcomes lost to censoring
Competing event Event that prevents the target event from subsequently occurring
Consistency Observed outcome under received treatment equals its corresponding potential outcome
Cumulative incidence function Probability of a specified event type by time in the presence of competing events
Double robustness Consistency if either of two nuisance models is correct, given other assumptions
Effective sample size (ESS) Weight-dispersion summary (∑w)2/∑w2(\sum w)^2/\sum w^2
Estimand Precisely defined target quantity
Estimator Rule mapping observed data to an estimate
Exchangeability No residual confounding/selection after the specified conditioning history
Hazard Instantaneous event rate among those still at risk
IPCW Inverse probability weighting for remaining observed/uncensored
IPTW Inverse probability weighting for treatment received
MAR Missingness independent of missing values conditional on observed data in the model
Marginal structural model Model for marginal potential-outcome means under treatment histories
Nelson–Aalen estimator Sum of event/risk-set increments estimating cumulative hazard
Negative control Exposure/outcome designed to probe shared bias without the primary causal relation
Positivity Required treatment/observation patterns occur within relevant histories
Propensity score Conditional probability of treatment given baseline covariates
Quantitative bias analysis Explicit propagation of assumptions about systematic error
RMST Expected event-free time restricted to a prespecified horizon
Sequential exchangeability Exchangeability at each decision time given observed past history
Stabilized weight Probability ratio with a numerator chosen to reduce variability/preserve a marginal structure
TATE Average treatment effect in an explicit external target population
Treatment–confounder feedback Prior treatment changes a later confounder of subsequent treatment

Authoritative references and further study

Causal estimands and g-methods

Survival and competing events

Missing data, transportability, and self-controlled designs

Measurement error, bias analysis, and reporting

  • Shaw PA et al. STRATOS guidance on measurement error and misclassification, Part 1. Statistics in Medicine 2020.
  • Fox MP, MacLehose RF, Lash TL. Applying Quantitative Bias Analysis to Epidemiologic Data, 2nd ed. Springer 2021.
  • Brown JP et al. Quantifying possible bias with quantitative bias analysis. BMJ 2024.
  • Cashin AG et al. Transparent reporting of observational studies emulating a target trial: the TARGET Statement. BMJ 2025.

Final takeaway

An advanced analysis is convincing when its target is explicit, its assumptions match the data-generating timeline, its diagnostics interrogate those assumptions, and its sensitivity analyses change scientifically meaningful inputs. Complexity is not rigor. A transparent marginal estimate with visible limitations is more useful than an opaque coefficient from a method whose target and support were never defined.