Marginal Structural Models and G-Methods

In plain language

A marginal structural model (MSM) is a statistical method for estimating how a treatment taken over time affects an outcome, when the very lab values and health measures that predict whether someone keeps taking the treatment are also changed by that treatment. Standard regression gets this wrong because adjusting for those measures blocks part of the treatment effect and introduces a new bias at the same time. An MSM solves the problem by creating a reweighted copy of the study population where treatment choices look as if they were made by a coin flip rather than by clinical judgment, and then estimating the effect in that cleaner copy.

Marginal structural models (MSMs) and the broader g-methods

(the parametric g-formula / g-computation, IPTW-estimated MSMs, and g-estimation of structural nested models) were developed by James Robins and colleagues to solve a problem that ordinary regression cannot: estimating the effect of a sustained or time-varying treatment when the confounders that govern later treatment are themselves consequences of earlier treatment. This is treatment-confounder feedback. The canonical example is sustained statin use and myocardial infarction: LDL cholesterol confounds the statin-MI relationship (high LDL prompts treatment and raises risk), but LDL after baseline is lowered by the statin and then drives the decision to continue or stop. LDL is simultaneously a confounder for later treatment and a mediator of earlier treatment. Adjusting for it in a standard time-dependent Cox model blocks part of the true effect (it is on the causal path) and introduces collider/selection bias; not adjusting for it leaves confounding. There is no regression specification that escapes both horns — g-methods are the resolution.

Core conceptual distinction and estimand

All three g-methods target a contrast of counterfactual outcomes under interventions on a treatment strategy — e.g., "always treat" vs "never treat," or a dynamic rule such as "treat once LDL exceeds a threshold." This is fundamentally different from the conditional hazard ratio a time-dependent Cox model reports. (1) IPTW-estimated MSMs reweight each person-time record by the inverse probability of the treatment actually received, given covariate and treatment history, creating a pseudo-population in which treatment is unconfounded; the MSM (typically a weighted pooled logistic or weighted Cox) is then fit without conditioning on the time-varying confounders, so it recovers a marginal (population-averaged) effect of the strategy. (2) The parametric g-formula models the joint evolution of confounders and outcome under the natural course, then simulates the outcome distribution under each hypothetical strategy by integrating over the time-varying covariate history — uniquely able to handle dynamic regimes and to report absolute risks and risk differences, not just ratios. (3) g-estimation of structural nested models directly estimates "blip" effects (the effect of one more unit of treatment at each time, possibly modified by covariate history) via estimating equations and is the natural tool when effect modification by time-varying covariates is the question. All rest on sequential exchangeability (no unmeasured confounding at each decision point given the measured past), positivity (every treatment level has nonzero probability at every level of history), and consistency / correct model specification. Always pre-specify the estimand — an intention-to-treat-like "initiate strategy" contrast vs an as-treated/per-protocol "initiate and remain on strategy" contrast (the latter requiring censoring at deviation and inverse-probability-of-censoring weights) — because the choice changes both the model and the interpretation.

Pros, cons, and trade-offs

- vs standard time-dependent Cox / pooled logistic with covariate adjustment: g-methods are the only family that is unbiased under treatment-confounder feedback and that can express sustained/dynamic strategy effects. A time-dependent Cox model that conditions on the post-baseline confounder is biased (over-adjustment for a mediator plus collider-stratification bias); one that omits it is confounded. Cost: g-methods demand a correctly specified treatment model (IPTW) or confounder+outcome models (g-formula), strict positivity at every time, and considerably more analytic and computational effort. Prefer g-methods whenever a post-baseline variable both predicts subsequent treatment and is affected by prior treatment; otherwise a standard adjusted model is simpler and adequate. - IPTW-MSM vs parametric g-formula: MSMs are simpler to communicate and need only a treatment (and censoring) model, but are inefficient and unstable when weights are extreme, and they handle dynamic regimes awkwardly. The g-formula is far more efficient, naturally accommodates dynamic strategies and absolute-risk contrasts, and degrades gracefully under near-positivity violations, but it requires modeling the entire confounder process and is vulnerable to the g-null paradox (under the sharp null with feedback, a parametric g-formula can be guaranteed misspecified). Doubly-robust hybrids (TMLE, LTMLE; see predictive-and-causal-ml-models-rwe) combine a treatment and an outcome model so that consistency holds if either is correct, and permit machine-learning nuisance estimation. - vs g-estimation of structural nested models: g-estimation is most robust for effect modification and avoids the g-null paradox, but software support is thin and the blip-function estimand is harder to communicate to clinical and HTA audiences. Prefer IPTW-MSM or g-formula for routine RWE; reserve g-estimation for explicit effect-modification questions.

When to use

Long-term or repeatedly-decided drug exposures where intermediate clinical variables (labs, vitals, adherence, disease activity, organ function) both predict future treatment and are altered by past treatment — statins and LDL, antiretrovirals and CD4/viral load, anticoagulants and renal function, immunosuppressants and disease activity. Also the analytic engine of per-protocol target-trial emulation for sustained strategies (used_with target-trial-emulation and clone-censor-weight-per-protocol), and the right tool for "when to start / when to switch" dynamic-strategy questions.

When NOT to use — and when it is actively misleading or dangerous

- Point-exposure / single-decision questions. If treatment is decided once at baseline and never revisited, there is no feedback; a new-user active-comparator design with a baseline propensity score (active-comparator-new-user, propensity-score-methods-psm-iptw) is simpler and more transparent. Reaching for an MSM here adds variance and a positivity burden for no benefit. - Structural positivity violations. If some patients can never receive a treatment level given their history (e.g., a contraindication that appears mid-follow-up and permanently rules out the drug), the "always treat" arm is undefined for them; IPT weights explode and the MSM estimates an artefact. Diagnose with the weight distribution and bounds on predicted treatment probabilities before trusting any number. Truncation hides, but does not fix, structural violations. - Sparse or extreme weights treated as if benign. Stabilized weights with a mean far from 1.0, a long right tail, or extreme maxima signal near-positivity violations or a misspecified treatment model; reporting the MSM without the weight diagnostics is the single most common way these analyses mislead. - The g-null paradox setting. Under a true sharp null with feedback, a naively specified parametric g-formula is guaranteed biased; if the null is plausible, prefer g-estimation or a doubly-robust estimator. - Unmeasured time-varying confounding. Sequential exchangeability is stronger than baseline exchangeability — it must hold at every decision point. Claims data rarely capture the labs/vitals that drive titration and switching, so an MSM built on claims alone can be confidently precise and wrong. Quantify with E-values or negative controls and consider linkage to EHR labs.

Data-source operational depth

- Claims (FFS or commercial): Excellent for constructing the long-format person-time skeleton — pharmacy fills (`fill_date`, `days_supply`, NDC) define the time-varying exposure, and diagnoses/procedures define time-varying confounders and outcomes. Failure modes: (i) Medicare Advantage / capitated person-time lacks fee-for-service claims, so a covariate or exposure that "turns off" may be unobserved rather than absent — restrict to enrollees with full A/B/D (or commercial medical+pharmacy) benefit and drop MA-only spans, or the treatment model is fit on phantom data. (ii) The lab/vital values that actually drive titration (LDL, eGFR, HbA1c, INR) are not in claims, so the most important time-varying confounders are missing — sequential exchangeability is then implausible without EHR linkage. (iii) Differential competing risks by exposure in elderly claims: death competes with the outcome and may differ by arm; handle the competing event explicitly (cause-specific vs subdistribution) rather than censoring on it. (iv) Immortal time in procedure/initiation studies: aligning the time grid so that exposure status at interval t uses only information available at the start of t prevents the look-ahead that manufactures immortal time. - EHR: Supplies the time-varying labs and vitals that claims lack, sharpening the treatment model and making sequential exchangeability defensible. Costs: irregular, visit-driven measurement (a covariate is observed only when the patient shows up), informative missingness, and incomplete capture of fills obtained out of network — last-observation-carried- forward or explicit imputation of the confounder process is usually required, and loss to follow-up must be treated as potentially informative (inverse-probability-of-censoring weights). - Registry: Often the cleanest structured longitudinal substrate (scheduled visits, adjudicated outcomes, disease severity) and ideal for validating a claims-based MSM; typically weak on complete pharmacy exposure, so link to claims for the full fill history and to a death index for the competing event. - Linked claims–EHR–vital records: The ideal substrate — EHR labs for the confounder process, claims completeness for exposure, reliable mortality for the competing event — at the price of linkage selection and reconciling order/fill/service dates before laying down the time grid.

Worked claims-style example (sustained statin strategy and MI, with feedback)

Question: among new statin initiators with type 2 diabetes, does "initiate and remain on a high-intensity statin" vs "initiate and remain on a moderate-intensity statin" change 3-year MI risk? Feedback is intrinsic: the statin lowers LDL, and the observed LDL drop then drives whether the clinician up-titrates, holds, or stops — so LDL at month t is a confounder for treatment at t and a mediator of treatment before t. (1) Time grid: one row per person-month from time zero (first qualifying fill) to first MI, death, disenrollment, or 36 months. (2) Exposure at t: on high- vs moderate-intensity, derived from `days_supply` coverage in month t with a 30-day grace period for stockpiling; deviation from the assigned strategy triggers as-treated censoring. (3) Baseline confounders (measured in the 365-day, FFS-observable lookback): age, sex, prior CVD, baseline LDL proxy, comorbidity score, baseline utilization. (4) Time-varying confounders (updated monthly from EHR-linked labs where available): LDL, eGFR, statin-intolerance flags, hospitalizations, new antidiabetic/antihypertensive starts. (5) Treatment model: monthly pooled logistic of "remains on high-intensity at t" on baseline + time-varying history; the stabilized weight for each record multiplies the cumulative product over t of [probability of the observed treatment given baseline-only history] / [probability given full history], with a parallel stabilized inverse-probability-of-censoring weight for the as-treated estimand. (6) MSM: weighted pooled logistic (or weighted Cox) of MI on the strategy indicator and follow-up time, not conditioning on the time-varying confounders; convert to a 36-month risk difference and ratio. (7) Diagnostics and sensitivity: report the stabilized-weight mean (should sit near 1.0), SD, and maximum; truncate at the 1st/99th percentiles and re-fit; vary the grace period; add a negative-control outcome; and, because death competes with MI in this older cohort, contrast cause-specific and subdistribution handling of the competing event. Build the same analysis with the parametric g-formula (gfoRmula) as a cross-check that does not rely on stable weights.

Interpreting the output

Consider the ART study above. A marginal structural Cox model fit with stabilized IPTW reports HR = 0.41 (95% CI 0.22–0.76) for 12-month opportunistic infection comparing sustained ART versus no ART.

Formal interpretation: The MSM estimate of HR 0.41 is the marginal effect of sustained antiretroviral therapy on the instantaneous rate of opportunistic infection in the IPTW pseudo-population — a population where the statistical association between monthly CD4 count and ART prescribing has been eliminated by reweighting. This targets the population-averaged (marginal) treatment effect under a sustained-exposure intervention, not a conditional effect within any CD4 stratum. The result is a valid causal estimate under three assumptions: (1) sequential exchangeability — at each time point, no unmeasured common cause of ART receipt and subsequent infection remains after conditioning on the covariates in the treatment model; (2) positivity — every patient with every covariate history had a non-zero probability of receiving or not receiving ART at each interval; and (3) consistency — the ART strategies compared are well defined and replicable. Extreme stabilized weights (above the 99th percentile) signal positivity violations and must be truncated with a sensitivity analysis.

Practical interpretation: After breaking the feedback loop between CD4 count and ART prescribing, patients who sustained ART had opportunistic infections arrive at less than half the rate of patients who received no ART. This estimate is not obtainable from a standard time-varying Cox model, which would either over-adjust for CD4 (blocking the beneficial causal pathway) or leave CD4 as an uncontrolled confounder.

Worked example

Scenario

A researcher wants to know whether staying on antiretroviral therapy (ART) for 12 months reduces the risk of opportunistic infection in HIV-positive patients. She has monthly records for 3 patients, tracking whether they took ART that month and their CD4 count at the start of each month. CD4 count is the core problem: low CD4 predicts both a higher chance of infection (outcome) and a higher chance the clinician prescribes or continues ART (treatment), making it a confounder. But CD4 also rises when a patient takes ART, meaning prior treatment changes the very confounder that predicts future treatment. This is treatment-confounder feedback. Standard logistic regression that adjusts for CD4 at each month will block part of the beneficial treatment effect (because rising CD4 is on the causal path from ART to lower infection risk) and simultaneously open a collider bias. Not adjusting leaves CD4 confounding. An MSM with IPTW escapes both horns.

Dataset

Monthly person-time records (3 patients, 2 months each). CD4 measured at start of month. ART = 1 if on therapy that month.

Steps

• Step 1 — See the feedback problem: Patient 1001 starts with CD4 = 180 (low, driving treatment) and their CD4 rises to 310 after one month on ART. Now CD4 at month 2 is both a consequence of month-1 treatment AND a predictor of month-2 treatment. That is treatment-confounder feedback in one patient.

• Step 2 — Understand why standard adjustment fails: If we include CD4 in a standard month-by-month logistic regression, we partially block the ART benefit (because higher CD4 caused by ART is on the pathway to lower infection) and open a collider path through unmeasured factors. If we leave it out, CD4 confounds the ART-infection relationship. There is no safe regression specification.

• Step 3 — Build the treatment model (denominator): For each person-month, fit a logistic model predicting the probability of receiving ART given that person's full history (baseline characteristics plus their current and prior CD4 values). This tells us how likely each treatment decision was given everything observed. Person 1002 at month 1 had CD4 = 220 and did NOT take ART; suppose the model gives P(ART=1 | CD4=220, history) = 0.60, so the probability of the observed treatment (no ART) is 1 - 0.60 = 0.40.

• Step 4 — Build the stabilized weight: Divide a simpler probability (from a model using only baseline covariates, not the time-varying CD4) by the full-history probability from Step 3. For the example record above, suppose the baseline-only model gives P(ART=1 | baseline) = 0.50, so P(no ART | baseline) = 0.50. The weight for that record is 0.50 / 0.40 = 1.25. A person who made an unlikely treatment choice given their CD4 gets a weight above 1 so their experience carries more influence in the pseudo-population.

• Step 5 — Accumulate weights over time: For each person, multiply the monthly weights together across all their follow-up months. This cumulative product is the stabilized IPTW for that person-month. A well-behaved set of weights has a mean close to 1.0 and no extreme values; extreme weights signal a positivity problem (someone almost certain to get one treatment level) or a misspecified treatment model.

• Step 6 — Fit the MSM on the pseudo-population: Run the outcome model (ART predicting infection) using the cumulative IPTW as weights, and do NOT include CD4 in this outcome model. Because the weights have already removed the statistical link between CD4 and treatment choice, the pseudo-population looks as if treatment was assigned without regard to CD4 — the feedback loop is broken. The coefficient on ART in this weighted model is the marginal effect of always taking ART versus not taking it, free from the feedback bias.

Result

In the pseudo-population created by IPTW, treatment choices are no longer driven by CD4. An MSM fit on these reweighted records estimates the population-averaged effect of sustained ART on infection risk without the over-adjustment and collider bias that would corrupt a standard CD4-adjusted model. The key diagnostic to report is the weight distribution: mean near 1.0 with a short right tail confirms the pseudo-population is well-behaved; a mean far from 1.0 or extreme maximum values signals that the approach may be unreliable for this dataset.

Runnable example

import numpy as np
import pandas as pd
import statsmodels.api as sm
import statsmodels.formula.api as smf

BASELINE = ["age", "sex", "prior_cvd"]          # fixed, measured in the lookback
TIMEVAR  = ["ldl", "egfr", "intolerance_flag"]  # updated each interval; drive titration AND are affected by treatment

def fit_iptw_msm(panel: pd.DataFrame, trunc=(1, 99)) -> dict:
    df = panel.sort_values(["person_id", "t"]).copy()

    # --- Treatment models for the weight (predict P(observed treatment | history) at each interval) ---
    denom_rhs = " + ".join(["t"] + BASELINE + TIMEVAR)           # full history -> denominator
    numer_rhs = " + ".join(["t"] + BASELINE)                     # baseline only -> numerator (stabilization)
    m_denom = smf.glm("treat ~ " + denom_rhs, df, family=sm.families.Binomial()).fit()
    m_numer = smf.glm("treat ~ " + numer_rhs, df, family=sm.families.Binomial()).fit()

    # P(treat actually received): use the fitted prob when treat=1, its complement when treat=0.
    pd_ = m_denom.predict(df); pn_ = m_numer.predict(df)
    df["p_denom"] = np.where(df["treat"] == 1, pd_, 1 - pd_)
    df["p_numer"] = np.where(df["treat"] == 1, pn_, 1 - pn_)

    # Stabilized weight = running product over time of (numerator / denominator) within person.
    df["ratio"] = df["p_numer"] / df["p_denom"]
    df["sw"] = df.groupby("person_id")["ratio"].cumprod()

    lo, hi = np.percentile(df["sw"], trunc)                       # truncate extreme weights (report both)
    df["sw_trunc"] = df["sw"].clip(lo, hi)

    # --- MSM: weighted pooled logistic of the outcome on strategy + time; NO time-varying confounders here. ---
    # For an as-treated estimand, multiply sw by an analogous stabilized inverse-probability-of-censoring
    # weight built from a model of remaining uncensored (not on treatment) given history.
    # var_weights (not freq_weights) carries the non-integer IPTW correctly; cluster SE fixes inference.
    msm = smf.glm("event ~ treat + t + I(t**2)", df,
                  family=sm.families.Binomial(), var_weights=df["sw_trunc"]).fit(
                  cov_type="cluster", cov_kwds={"groups": df["person_id"]})  # robust SE for the pseudo-population
    return {
        "weight_mean": df["sw"].mean(), "weight_sd": df["sw"].std(), "weight_max": df["sw"].max(),
        "trunc_at": (lo, hi), "msm": msm,
        "log_or_per_interval": msm.params["treat"],
    }

res = fit_iptw_msm(panel)
print(f"stabilized weight  mean={res['weight_mean']:.3f}  sd={res['weight_sd']:.3f}  max={res['weight_max']:.2f}")
print(res["msm"].summary())

library(ipw); library(geepack); library(data.table)
setDT(panel)

## (A) Stabilized IPT weights for time-varying treatment, then the MSM ----------------------
w <- ipwtm(
  exposure   = treat,
  family     = "binomial", link = "logit",
  numerator  = ~ t + age + sex + prior_cvd,                 # baseline only -> stabilization
  denominator= ~ t + age + sex + prior_cvd + ldl + egfr + intolerance_flag,  # + time-varying history
  id         = person_id, timevar = t, type = "all",
  data       = as.data.frame(panel)
)
panel[, sw := w$ipw.weights]
# For an as-treated estimand, multiply sw by an analogous stabilized IPC weight (ipwtm on the
# censoring indicator) built from a model of remaining uncensored given history.
panel[, sw_tr := pmin(pmax(sw, quantile(sw, .01)), quantile(sw, .99))]   # truncate; report both
cat(sprintf("sw mean=%.3f sd=%.3f max=%.2f\n", mean(panel$sw), sd(panel$sw), max(panel$sw)))

## MSM = weighted pooled logistic; NO time-varying confounders on the RHS. Cluster-robust SE via GEE.
msm <- geeglm(event ~ treat + t + I(t^2), family = binomial,
              weights = sw_tr, id = person_id, corstr = "independence", data = panel)
summary(msm)

## (B) Parametric g-formula for the "always high-intensity vs always moderate" contrast -------
library(gfoRmula)
gf <- gformula_survival(
  obs_data = panel, id = "person_id", time_name = "t",
  covnames = c("ldl", "egfr", "intolerance_flag", "treat"),
  covtypes = c("normal", "normal", "binary", "binary"),
  covparams = list(covmodels = c(
    ldl ~ lag1_treat + lag1_ldl + t,
    egfr ~ lag1_treat + lag1_egfr + t,
    intolerance_flag ~ lag1_treat + t,
    treat ~ lag1_treat + ldl + egfr + intolerance_flag + age + sex + prior_cvd + t)),
  outcome_name = "event", ymodel = event ~ treat + ldl + egfr + t + age + sex + prior_cvd,
  intvars = list("treat", "treat"),
  interventions = list(list(c(static, rep(1, max(panel$t) + 1))),
                       list(c(static, rep(0, max(panel$t) + 1)))),
  int_descript = c("always high-intensity", "always moderate"),
  nsimul = 10000, seed = 1)
print(gf)   # cumulative-incidence risk difference / ratio under each strategy