Marginal Structural Models: A Practical Guide for Clinical Researchers

Standard regression adjusts for confounders. Propensity scores balance them. But what happens when your confounder is itself affected by treatment? When a patient's CD4 count at month 6 depends on whether they received antiretrovirals at month 0, adjusting for it will block the very pathway you want to study. Marginal Structural Models solve this by weighting — not adjusting — through time-varying confounders affected by prior treatment.

The Problem: Time-Varying Confounders Affected by Treatment

In longitudinal studies, patients receive treatment over time, and their characteristics change. Those changes are often confounders for subsequent treatment decisions. But they're also affected by prior treatment. This creates a dilemma:

• If you adjust for the time-varying confounder in standard regression, you block part of the causal effect (mediation bias)
• If you don't adjust for it, residual confounding remains
• If you try stratification, you create artificial selection bias through collider stratification

This isn't a rare edge case. It's the standard situation in chronic disease management, HIV treatment, oncology, and any setting where treatment decisions are made sequentially based on evolving patient status.

The key insight of MSMs, formalized by Robins, Hernán, and Brumback in 2000, is to restructure the data using weights so that, at each time point, treatment assignment becomes independent of past covariate history — mimicking sequential randomization.

How MSMs Work: The Intuition

MSMs don't change the model. They change the data. By reweighting each patient at each time point, MSMs create a pseudo-population where:

1. At each time point, treatment assignment is independent of prior covariate history
2. The relationship between treatment and outcome can be estimated without bias from time-varying confounders
3. The "marginal" effect represents the causal effect of the treatment regime, not the effect in a particular subgroup

The weights are calculated from the probability of each patient receiving their actual treatment at each time point, conditional on their prior treatment and covariate history. These are called inverse probability of treatment weights (IPTW).

The Weight Formula

For a patient with treatment history A₀, A₁, ..., Aₖ and covariate history L₀, L₁, ..., Lₖ, the weight at time k is the product of the inverse probabilities of receiving each observed treatment given past treatment and covariates:

This weight has an important interpretation: patients who received treatment in ways that were unexpected given their history get higher weights. Patients who received treatment in expected ways get weights close to 1. In the weighted pseudo-population, the time-varying confounders are no longer associated with treatment assignment.

Stabilized Weights: Why You Need Them

The raw weights described above work, but they have a serious problem: they can be extremely large for patients with very low or very high treatment probabilities. This inflates variance and makes estimates unstable.

The solution is stabilized weights, which multiply the inverse probability by a numerator that captures the marginal probability of the observed treatment:

Stabilized weights have two important properties:

• Mean close to 1: Unlike raw weights which average to the sample size, stabilized weights should average near 1. If the mean is far from 1, something is wrong.
• Reduced variance: By multiplying by the marginal probability in the numerator, extreme weights are compressed without introducing bias.

Clinical Example: ART and Mortality in HIV

This is the canonical MSM application, and for good reason — it's where the method was born. Consider estimating the causal effect of cumulative antiretroviral therapy (cART) on mortality in HIV patients, using observational cohort data.

The Setup

At each clinic visit (quarterly), a physician decides whether to initiate or modify cART based on the patient's current CD4 count and viral load. These decisions evolve over time:

• CD4 count at time t is a confounder for mortality at time t+1
• CD4 count at time t is also affected by cART at time t-1
• Standard regression adjusting for CD4 blocks the causal effect of cART

The MSM Solution

Step-by-step:

1. Model treatment assignment: Fit a logistic regression for cART initiation at each time point, conditional on prior cART history, current CD4, viral load, age, and clinical site. These produce the denominator probabilities.
2. Model the numerator: Fit a simpler logistic regression for cART initiation conditional only on prior cART history (no covariates). These produce the numerator probabilities.
3. Compute stabilized weights: For each patient at each time point, SW = numerator / denominator. Multiply weights over time to get the cumulative stabilized weight.
4. Fit the MSM: Fit a weighted Cox model or pooled logistic regression with the stabilized weights. The model relates cumulative cART exposure to mortality — without adjusting for CD4 (which would block the effect).

Clinical Example: Chemotherapy Timing in Oncology

A second example that illustrates MSMs in a different context: does early initiation of adjuvant chemotherapy after surgery improve overall survival in stage II colon cancer?

The confounding structure: patients with postoperative complications (a confounder for survival) are less likely to receive chemotherapy on time. But postoperative complications may also be related to tumor characteristics that influenced surgical approach. Standard regression adjusting for postoperative complications would be biased.

MSM approach: at each post-operative week, model the probability of chemotherapy initiation given prior complications, performance status, and tumor stage. Weigh patients to create a pseudo-population where chemotherapy timing is independent of these factors. The MSM estimates the survival under "initiate chemo by week 8" versus "initiate chemo by week 16" — two hypothetical regimens.

Model Fitting: Weighted Regression

Once weights are computed, you fit a marginal structural model using weighted regression. The model form depends on your outcome:

Critical implementation note: when using weighted Cox models, use robust variance estimators (sandwich estimators) or robust standard errors. The weights violate the independence assumption that standard Cox standard errors rely on. In R, the survey package handles this automatically. In Python, statsmodels with cov_type='HC0' works for weighted GLMs.

The Causal Assumptions

MSMs require three key assumptions, and violating any one of them invalidates the causal interpretation:

Common Pitfalls

Implementation Checklist

Reporting Checklist

MSMs vs. Other Causal Methods

Software Implementation

R

# Load packages
library(survey)
library(survival)

# Step 1: Fit treatment models at each time point
# Denominator model (treatment ~ prior treatment + confounders)
denom_fit <- glm(
  treat_t ~ lag_treat + cd4 + viral_load + age,
  family = binomial,
  data = long_data
)
denom_prob <- predict(denom_fit, type = "response")

# Numerator model (treatment ~ prior treatment only)
num_fit <- glm(
  treat_t ~ lag_treat,
  family = binomial,
  data = long_data
)
num_prob <- predict(num_fit, type = "response")

# Step 2: Compute stabilized weights
long_data$sw <- ifelse(
  long_data$treat_t == 1,
  num_prob / denom_prob,
  (1 - num_prob) / (1 - denom_prob)
)

# Step 3: Cumulative weight (product over time)
long_data <- long_data %>%
  group_by(id) %>%
  mutate(cum_sw = cumprod(sw))

# Step 4: Truncate
long_data$sw_trunc <- pmin(long_data$cum_sw, 10)

# Step 5: Fit weighted Cox model
design <- svydesign(
  ids = ~1,
  weights = ~sw_trunc,
  data = long_data
)
svycoxph(Surv(time, death) ~ treat_cum + age, data = long_data,
         design = design)

Python

import pandas as pd
import numpy as np
from sklearn.linear_model import LogisticRegression
from lifelines import CoxPHFitter

# Step 1: Fit treatment models
def fit_treatment_model(data, covariates, treat_col):
    """Fit denominator (full) and numerator (marginal) models."""
    X_denom = data[covariates].values
    y = data[treat_col].values
    denom = LogisticRegression(max_iter=1000)
    denom.fit(X_denom, y)
    
    X_num = data[['lag_treat']].values
    num = LogisticRegression(max_iter=1000)
    num.fit(X_num, y)
    
    denom_prob = denom.predict_proba(X_denom)[:, 1]
    num_prob = num.predict_proba(X_num)[:, 1]
    return denom_prob, num_prob

# Step 2: Compute stabilized weights
denom_prob, num_prob = fit_treatment_model(
    data, ['lag_treat', 'cd4', 'viral_load', 'age'], 'treat_t'
)
sw = np.where(data['treat_t'] == 1,
              num_prob / denom_prob,
              (1 - num_prob) / (1 - denom_prob))

# Step 3: Cumulative weight
data['sw'] = sw
data['cum_sw'] = data.groupby('id')['sw'].cumprod()

# Step 4: Truncate
data['sw_trunc'] = data['cum_sw'].clip(upper=10)

# Step 5: Fit weighted Cox model
cph = CoxPHFitter()
cph.fit(data[['time', 'death', 'treat_cum', 'age', 'sw_trunc']],
        duration_col='time', event_col='death',
        weights_col='sw_trunc')
cph.print_summary()

When to Use MSMs — and When Not To

Ideal Settings

• Longitudinal data with sequential treatment decisions — HIV cohorts, oncology follow-up, chronic disease management, ICU care
• Time-varying confounders affected by prior treatment — this is the specific problem MSMs were designed to solve
• Large sample sizes — MSMs need stable weight estimates; >500 patients per treatment-time combination is ideal
• Well-defined treatment initiation points — clear temporal ordering of covariate measurement and treatment decisions

Poor Settings

• Small samples (<100) — weights will be unstable; consider parametric G-computation or stratification instead
• No time-varying confounders — if confounders are only at baseline, standard IPW or PSM is simpler and equivalent
• Positivity violations — if treatment is deterministic for some covariate values, MSMs break down
• Time-varying treatments without clear intervals — MSMs require discrete time points; continuous treatment requires specialized methods

Key References

• Robins JM, Hernán MA, Brumback B. Marginal Structural Models and Causal Inference in Epidemiology. Epidemiology. 2000;11(5):550-560.
• Hernán MA, Brumback B, Robins JM. Marginal Structural Models to Estimate the Joint Causal Effect of Nonrandomized Treatments. JASA. 2001;96(454):440-448.
• Suarez D, Borras R, Basagana X. Comparison of marginal structural cox models for the analysis of a time-varying exposure: a simulation study. Emerg Themes Epidemiol. 2007;4:14.
• Cain LE, Robins JM, Lanoy E, et al. When to initiate combination antiretroviral therapy to reduce mortality and AIDS-defining illness in HIV-infected persons in developed countries: an observational study. PLoS Med. 2011;8(8):e1001084.
• Robins JM. Marginal Structural Models versus Structural Nested Models as tools for causal inference. In: Statistical Models in Epidemiology, the Environment, and Clinical Trials. Springer; 2000:95-133.
• Hernán MA, Robins JM. Causal Inference: What If. Boca Raton: Chapman & Hall/CRC, 2020. Chapter 12. Free online