Marginal structural models and time-varying IPTW

Why this method gets its own page

There are about three ways you land here from a paper. The Methods say patients who switched treatment. They say time-updated covariates. Or they say per-protocol analysis, and mean something more than dropping the non-adherent. There is only one way you land here from your own data: patients start, stop and switch treatment during follow-up, instead of being sorted into two arms on day one and staying there.

Data of that shape has a feature no standard regression handles. The marker that decides whether a patient gets treated the second time — CD4 count, blood pressure, renal function, disease activity — is itself a consequence of the first treatment. It is a confounder of the later treatment and a consequence of the earlier one at the same time. A standard regression offers exactly two options for such a variable: put it in, or leave it out. This page is about why both are wrong, and what the third route looks like.

Three places on this site send readers here: time-varying covariates notes that a time-varying mediator-cum-confounder cannot be handled by a standard Cox model, IPTW notes that only time-varying weighting works once exposure changes over time, and target trial emulation notes that a per-protocol contrast has to deal with treatment that changes. This page is where those three lines end, so it does not hand you back to them — the weight formula, runnable code and a simulation with a known truth are all below.

Treatment and confounder feeding each other

Shrink the time axis to its minimum — two clinic visits — and the whole problem is already there:

• Visit 0: the patient is treated or not, written $A_0$. In this simulation it is randomised.
• Visit 1: a clinical marker $L_1$ is measured, and the second treatment $A_1$ is decided from it.
• Outcome $Y$: an event within, say, one year.

$L_1$ sits in an awkward place. It is a consequence of $A_0$ (treatment moved the marker) and a confounder of $A_1$ (the marker drove the second treatment decision, and it also predicts the outcome). That is treatment-confounder feedback: treatment changes the covariate, and the covariate turns around and decides the next treatment.

The graph also carries an unobserved $U$ that causes both $L_1$ and $Y$. Clinically this is everywhere: unrecorded disease severity, performance status, social circumstances, family support. They move the marker, they move the outcome, and the database has no column for them.

The generating process writes one line for each arrow above:

The cohort that comes out: 20,000 people, event rate 13.4%. Among those untreated at visit 0, 40.4% are $L_1$ positive, against 68.5% among the treated — that is the $A_0 \to L_1$ arrow showing up in the data. And 40.3% of the $L_1$ negative get the second treatment against 70.6% of the $L_1$ positive — that is $L_1 \to A_1$. Both hold at once, so feedback holds.

Putting it in and leaving it out are both wrong

Leaving it out is no better: $L_1$ is a genuine confounder of $A_1$ (it predicts $A_1$, and it predicts $Y$ through $U$), so ignoring it leaves confounding by indication untouched.

Here are the three routes, as the joint log odds ratio for treatment at both visits against treatment at neither:

The directions are predictable in advance, and the data obliges: adjusting for $L_1$ overshoots (-1.427, further from the null than the truth), no adjustment undershoots (-0.975, closer to the null), and weighting lands between them and nearest the truth (-1.173, off by 0.074). Three routes, three directions — that is what this simulation exists to show.

Weighting, the third route: what an MSM estimates

Start from the estimand. In potential-outcome notation, $Y^{a_0, a_1}$ is whether this person would have the event if they were assigned $a_0$ at visit 0 and $a_1$ at visit 1. Four static regimes, four whole-cohort risks:

A marginal structural model puts a model directly on those four risks, rather than on the observed outcome conditional on covariates:

Marginal means the left-hand side is not conditional on any covariate — it is the risk in the whole cohort, not the risk in some kind of patient. Structural means the left-hand side is a potential outcome, not the observed one. And $\beta_1 + \beta_2$ is the log odds ratio for both-treated against neither-treated in the table above, which is the number this page has been chasing.

How do you estimate that from observational data? Weight each person by the inverse of the probability of the treatment they actually received, which builds a pseudo-population in which treatment at each visit is unrelated to the covariates at that visit — as if someone had flipped a coin at every clinic appointment.

What stabilised weights look like

Unstabilised weights are the running product of inverse treatment probabilities. They explode: as soon as somebody receives, at some visit, the treatment almost nobody in their condition receives, the denominator approaches zero and that one person dominates the analysis.

Stabilised weights put back, in the numerator, a probability that uses past treatment only, not covariates:

The bars mean history up to that visit. The denominator removes the confounding carried by the covariates; the numerator puts back the part of the treatment process that depends only on past treatment, which is what keeps the weights near 1 instead of spread out. Covariates must never enter the numerator — that would add back exactly what the denominator removed.

In this simulation the two models are: numerator A1 ~ A0, denominator A1 ~ A0 + L1.

The weights themselves: mean 1.000, standard deviation 0.284, range 0.653 to 1.517, effective sample size 18,505, which is 92.5% of the nominal n. Stabilisation is why they are this well behaved. The unstabilised version of the same weights — numerator 1 instead of the marginal treatment probability — runs from 2.76 to 7.43 around a mean of 4.00, a maximum 4.90 times the stabilised one and a standard deviation 5.20 times as large. A factor of 2.00 of that is bookkeeping — the visit-0 numerator is a constant here because $A_0$ was randomised — but the effective sample size still falls to 17,597, or 88.0% of the nominal n against 92.5% for the stabilised weights.

This table repays a second look, because it turns what weighting does into something visible. Every one of the four patterns with a weight above 1 is a person whose $L_1$ and $A_1$ disagree — marker positive but untreated, or marker negative but treated — 6,881 people in total. They are a minority in the real world, and in the pseudo-population they are inflated to match everyone else, which is what decouples treatment from the marker. The four patterns with a weight below 1 are the mirror image: the majority who did follow the marker, shrunk.

Put differently: weighting prises the confounding open using the patients who did not follow clinical logic. That is also why positivity is this method’s weak point. If not a single marker-positive patient goes untreated, that cell is empty, the weight tends to infinity or cannot be estimated at all, and the mean weight will not show you — you have to look at the ESS, at the maximum, and at the count in every pattern, which is the table above.

Running it: long format, two stages

The implementation has two stages with a set of weights in between. This is where an MSM departs furthest from ordinary regression: you fit two models, and you never once look at the coefficients of the first — it exists only to produce the weights.

# Build the wide table d. U is the unmeasured confounder -- without it,
# conditioning on L1 merely over-adjusts, the naive model looks correct,
# and the page is refuted by its own numbers.
set.seed(20260823)
n  <- 20000
U  <- rnorm(n)                                   # the unmeasured confounder
A0 <- rbinom(n, 1, 0.5)
L1 <- rbinom(n, 1, plogis(-0.5 + 1.5 * A0 + 1.2 * U))
A1 <- rbinom(n, 1, plogis(-0.5 + 1.2 * L1 + 0.3 * A0))
Y  <- rbinom(n, 1, plogis(-1.5 - 0.7 * A0 - 0.7 * A1 + 1.0 * U))
d  <- data.frame(A0, L1, A1, Y)

# ── Stage 1: the weight models, on a long table (one row per person-visit) ──
# d is the wide table: one row per person, columns A0 / L1 / A1 / Y. Reshape
# first, because a weight is computed at each visit and multiplied through —
# a long table lets the loop be written once.
n <- nrow(d)
long <- rbind(
  data.frame(id = 1:n, visit = 0, A = d$A0, L = 0,     lagA = 0),
  data.frame(id = 1:n, visit = 1, A = d$A1, L = d$L1,  lagA = d$A0)
)

K  <- 1                    # index of the last visit
sw <- rep(1, n)            # one running weight per person, starting at 1

for (k in 0:K) {
  dk <- long[long$visit == k, ]
  dk <- dk[order(dk$id), ]

  # Numerator: past treatment only. Denominator: past treatment plus the
  # covariate measured now. Never put covariates in the numerator — that adds
  # back exactly the confounding the denominator removed.
  num <- glm(A ~ lagA,     data = dk, family = binomial())
  den <- glm(A ~ lagA + L, data = dk, family = binomial())

  pn <- predict(num, type = "response")   # P(A = 1 | past treatment)
  pd <- predict(den, type = "response")   # P(A = 1 | past treatment + covariate)

  # Each person contributes the probability of the treatment they actually
  # received, which is not P(A = 1) for the untreated.
  sw <- sw * ifelse(dk$A == 1, pn, 1 - pn) / ifelse(dk$A == 1, pd, 1 - pd)
}

# ⚠️ In this simulation lagA and L are constant at visit 0, so both models
# collapse to an intercept and that period's weight is exactly 1 (R reports the
# constant terms as NA coefficients, which is expected). Real data has baseline
# covariates L0 at visit 0, the denominator is A ~ L0, and the weight is not 1.

summary(sw)                       # the mean should sit near 1; far from it means a misfit
sum(sw)^2 / sum(sw^2)             # effective sample size

# ── Stage 2: the weighted outcome model — this is the MSM ──────────────────
# Treatment on the right-hand side and nothing else. No L1. The covariate
# appears only inside the stage-1 denominator.
# ⚠️ suppressWarnings is deliberate and is not hiding a misfit: once weighted,
# w_i * y_i is no longer a whole number of events, so glm prints
# non-integer #successes in a binomial glm!. The point estimates are the
# correct weighted score solutions; only glm's binomial bookkeeping objects.
msm <- suppressWarnings(
  glm(Y ~ A0 + A1, data = d, family = binomial(), weights = sw)
)
coef(msm)["A0"] + coef(msm)["A1"]   # both visits treated vs neither

# ⚠️ Do not read the standard errors from summary(msm). See the next section.

import numpy as np
import pandas as pd
import statsmodels.api as sm
import statsmodels.formula.api as smf
# Build the wide table d. U is the unmeasured confounder -- without it,
# conditioning on L1 merely over-adjusts, the naive model looks correct,
# and the page is refuted by its own numbers.

rng = np.random.default_rng(20260823)
n   = 20000
# numpy's stream is not R's, so the same seed does not draw the same
# people and these numbers are close to the table above, not equal to it.
expit = lambda x: 1 / (1 + np.exp(-x))
U  = rng.standard_normal(n)
A0 = rng.binomial(1, 0.5, n)
L1 = rng.binomial(1, expit(-0.5 + 1.5 * A0 + 1.2 * U))
A1 = rng.binomial(1, expit(-0.5 + 1.2 * L1 + 0.3 * A0))
Y  = rng.binomial(1, expit(-1.5 - 0.7 * A0 - 0.7 * A1 + 1.0 * U))
d  = pd.DataFrame({"A0": A0, "L1": L1, "A1": A1, "Y": Y})

# d is the wide DataFrame with columns A0 / L1 / A1 / Y
n = len(d)
long = pd.concat([
    pd.DataFrame({"id": range(n), "visit": 0, "A": d.A0, "L": 0,    "lagA": 0}),
    pd.DataFrame({"id": range(n), "visit": 1, "A": d.A1, "L": d.L1, "lagA": d.A0}),
])

sw = np.ones(n)
for k in (0, 1):
    dk = long[long.visit == k].sort_values("id")
    # lagA and L are constant at visit 0 and patsy drops them; real data would
    # use the baseline covariates here instead.
    rhs_num = "lagA" if dk.lagA.nunique() > 1 else "1"
    rhs_den = " + ".join(c for c in ("lagA", "L") if dk[c].nunique() > 1) or "1"
    pn = smf.glm(f"A ~ {rhs_num}", dk, family=sm.families.Binomial()).fit().fittedvalues
    pd_ = smf.glm(f"A ~ {rhs_den}", dk, family=sm.families.Binomial()).fit().fittedvalues
    a = dk.A.to_numpy()
    sw = sw * np.where(a == 1, pn, 1 - pn) / np.where(a == 1, pd_, 1 - pd_)

# var_weights, not freq_weights: a weight is not a row repeated w times
msm = smf.glm("Y ~ A0 + A1", d, family=sm.families.Binomial(),
              var_weights=sw).fit()
print(msm.params["A0"] + msm.params["A1"])
# The bse printed here is equally unusable; bootstrap it yourself

The truth is not the structural coefficient: non-collapsibility

The line that generates $Y$ gives $A_0$ and $A_1$ a coefficient of -0.70 each, summing to -1.40. But the truth this page has been comparing against is -1.247, not -1.40. The two differ, and neither is a mistake.

They are different parameters:

• -0.70 is a conditional parameter — how much treatment moves the log odds holding $U$ fixed.
• -1.247 is a marginal parameter — how much the log odds of the whole-cohort risk differ between treating everybody and treating nobody.

The odds ratio is non-collapsible: even with no confounding at all, dropping a variable that predicts the outcome (here $U$) shrinks the coefficient towards the null. That is a property of the odds ratio scale, not a bias — risk differences and risk ratios do not behave this way.

An MSM is not magic; it is also a model

In the table above, the MSM is still 0.074 from the truth. Part of that residual gap is not noise. It is misspecification — of the working model this page itself wrote down.

That model assumes the effect of $A_0$ does not depend on the value of $A_1$: the right-hand side has two main effects and no interaction term. The truth disagrees. With $A_1$ held at untreated, the $A_0$ contrast is -0.612; with $A_1$ held at treated it is -0.635, a difference of 0.022. The working model estimates a single parameter for two quantities that are not equal, so what it reports is a compromise between them.

This is worth stating because it is where MSMs are most often misunderstood. Weighting handles the confounding; it does not specify your outcome model for you. You still choose whether to include an interaction, whether to model time flexibly, whether dose enters as a continuous variable. A misfit weight model distorts the weights; a misfit outcome model moves the estimand itself. They are two independent routes to a wrong answer.

Why the variance cannot be read off the model

Two honest routes:

• Bootstrap, used here. Resample people, and refit the weight models inside every replicate — refitting only the outcome model and carrying the weights through as fixed constants reproduces error 2 above inside the bootstrap. This page uses 500 replicates, and every interval in the tables comes from them.
• Robust (sandwich) standard errors. Cheap, one option in most software, but they treat the weights as known and are usually conservative for stabilised weights. Computed here as well: the MSM $A_0$ sandwich SE is 0.045 against a bootstrap SE of 0.045, and for $A_1$ it is 0.044 against 0.046. They are close on this dataset, which is a fact about this dataset rather than a general result.

IPCW: a second weight, and the two multiply

Everything above deals with who got treated. Real follow-up data has a second selection problem: who is still under observation. Patients are lost, withdraw, move hospitals — and those who leave are rarely like those who stay. If the probability of leaving depends on covariates and treatment, then still being followed is itself a conditioned-on variable, with the same structure as confounding.

The remedy has the same shape: inverse probability of censoring weighting (IPCW). The formula runs parallel to IPTW, with the probability of receiving treatment replaced by the probability of remaining uncensored:

In practice it is one more logistic model — the outcome being whether this person was censored in this interval — and everyone still under follow-up is weighted by the inverse of their probability of having stayed. The people who remain speak for those like them who left.

This simulation has no censoring, so it demonstrates only the IPTW half. Omitting the other half on real data quietly assumes that dropout is completely random — the non-informative censoring assumption from censoring and truncation, and in a population that switches and stops treatment it is usually the least defensible assumption in the analysis.

How to read it in a paper

An MSM usually survives into a paper as three or four sentences and a table. Five things to check:

1. Whether the variables in the weight models are listed, numerator and denominator separately. Adjusted using IPTW says nothing. Covariates leaking into the numerator is a common implementation error, and it shows up on no diagnostic plot.
2. Whether weights were computed at every visit, or only at baseline. An analysis with baseline weights only is not an MSM; it is baseline IPTW, and it cannot handle feedback.
3. Whether the weight distribution is reported. At minimum the mean (stabilised weights should sit near 1), the maximum, and the effective sample size. A mean alone hides the extremes; the ESS here is 92.5% of the nominal n, and dropping below about half is a signal to go back and check overlap.
4. Whether weights were truncated, where, and with a sensitivity analysis. Truncation trades a little bias for variance and is a defensible choice, but it must be stated in advance and the untruncated result shown alongside.
5. Where the standard errors came from. You want to see the word bootstrap or robust or sandwich. A weighted regression with no statement about variance is most likely reporting the unusable one.

And one that gets skipped: the outcome model’s right-hand side should carry treatment (and any baseline covariates) and nothing time-varying. An MSM with $L_k$ in its outcome model is not an MSM; it is the analysis this page warned about three sections ago.

Common misuse

Where this page sits

• Where the weighting idea starts — see propensity score weighting (IPTW). That page covers the single-timepoint version; its diagnostics (balance tables, weight distributions, ESS, truncation) apply here unchanged, except that they have to be done at every visit.
• How to draw the arrows and decide what to condition on — see DAGs and the structure of confounding. The collider bias on this page is that page’s collider rule instantiated along a time axis.
• Time-varying covariates in a Cox model — see time-varying covariates and the long table. That page teaches how the long table is built and what tmerge() does, which is the data work that precedes stage 1 here. The two part company on whether the covariate is affected by exposure: if it is not, the standard Cox model on that page is right; if it is, weighting is required.
• Where a per-protocol contrast lands — see target trial emulation. That page defines the estimand (intention-to-treat or per-protocol) and time zero; this page supplies the estimator for the per-protocol one.
• The censoring assumption — see censoring and truncation. The IPCW section replaces how non-informative censoring is handled, not the requirement itself.
• Where it shows up clinically — see the database study chapter; dispensing records in claims data are time-varying exposure by nature.

Rerun every number on this page

/opt/homebrew/bin/Rscript figures/scripts/B6-08-msm.R

Report a problem on this page