Marginal Structural Models: When Today's Confounder Is Yesterday's Treatment Effect
On this page
อ่านฉบับภาษาไทย (Thai version)
Abstract
In intensive care, severe organ dysfunction makes a drug more likely, and each dose lowers the chance of organ dysfunction at the next review. Organ dysfunction is a time-varying confounder affected by earlier treatment. Leaving it out of an ordinary regression leaves every dose confounded; putting it in blocks the effect through organ dysfunction and opens a spurious path (a collider path, where two arrows meet) through unmeasured frailty. A marginal structural model (MSM) avoids both by weighting each patient by the inverse probability of the treatment history received. A hand example of 800 patients and a simulated cohort of 3,000 show both failures. In the cohort, naive regressions give +0.086 or -0.013 per treated day; the MSM gives -0.057 against a true -0.050, and -0.172 for always versus never treated against a true -0.150. This article shows how to build cumulative stabilised weights (multiplied over days and rescaled to average about 1), fit and check an MSM, and report its limits.
A regression the attending does not believe
In an intensive care unit, a fictional organ-protective drug, drug X, is reconsidered every two days: on day 0, day 2 and day 4. Severe organ dysfunction at a review makes the team more likely to give it. A dose given at one review lowers the chance of severe organ dysfunction at the next. The outcome that matters is death by day 28.
A critical care fellow assembles the unit's records into a cohort of 3,000 patients. The fellow regresses death by day 28 on the number of days each patient received drug X. Unadjusted, each treated day adds 0.068 to the risk of death, as if the drug were harmful. Adjusted for age and for severe organ dysfunction on each decision day, the estimate shrinks to 0.011 per treated day (95% confidence interval -0.006 to 0.027).
The attending believes neither number. Drug X is thought to work mainly by preventing organ dysfunction at the next review, so adjusting for organ dysfunction adjusts away the very effect the unit wants to measure. Yet leaving it out leaves the sickest patients concentrated among the treated, because organ dysfunction is why they were treated.
Both objections are right, and no choice of covariates in an ordinary regression answers both. Marginal structural models (MSMs) were built for this problem: models of the risk each treatment plan would produce, fitted with weights instead of covariate adjustment. The cohort here is simulated data, so the true answer is known. Each treated day lowers the risk of death by day 28 by 0.050.
The question: a contrast of treatment regimes
Before any model, the question needs a form that a trial could answer. A treatment regime, or treatment strategy, is a rule that sets the treatment on every decision day. The two regimes of main interest are always treated (drug X on days 0, 2 and 4) and never treated (no drug on any of the three days). A regime can also depend on the patient's state, such as "give drug X whenever severe organ dysfunction is present", but this article uses fixed regimes.
Write $A_t$ for treatment on decision day $t$ (1 = drug X given, 0 = not), where $t$ = 0, 1 and 2 stand for days 0, 2 and 4. A regime is a sequence $\bar a = (a_0, a_1, a_2)$, and the bar marks a history. Always treated is $\bar a = (1, 1, 1)$, and never treated is $\bar a = (0, 0, 0)$. Write $Y$ for death by day 28 (1 = died).
The estimand is the quantity a study aims to estimate. Here it is the risk of death by day 28 if every patient had been always treated, minus the risk if every patient had been never treated. That is the comparison a two-arm randomised trial of the two regimes would make, with every patient followed from the day-0 review.
In the usual terms it is an average treatment effect (ATE), an effect averaged over a whole population: here every patient in the ICU cohort, from day 0. The ATE here asks what would have happened had every patient followed each regime throughout. It is neither the effect of one dose nor the difference between patients who happened to follow each regime.
In the simulated cohort, 471 patients received drug X on all three days and 588 on none. The rest were treated on one day (1,009 patients) or on two days (932 patients). Comparing the 471 with the 588 gives a confounded answer, because the two groups differed in organ dysfunction on every decision day.
Drawing the problem over time
Write $L_t$ for severe organ dysfunction at the review on decision day $t$ (1 = present), recorded before that day's decision; from here on, organ dysfunction means this severe form. Write $V$ for a baseline covariate, here an age group (1 = aged 65 or older, 0 = younger), and $U$ for frailty, a patient trait that nobody in the unit records. The diagram shows how these variables cause one another over the first two decision days.
One variable, two roles
Follow $L_1$, severe organ dysfunction on day 2. The arrow $A_0 \to L_1$ says the day-0 dose changes it. The arrows $L_1 \to A_1$ and $L_1 \to Y$ say it drives the day-2 decision and the risk of death.
For the day-0 dose, $L_1$ is a mediator: a variable on a causal path from a treatment to the outcome, here $A_0 \to L_1 \to Y$. For the day-2 dose, $L_1$ is a confounder, a common cause of that dose and of death through $A_1 \leftarrow L_1 \to Y$. It is never both for the same dose. It is a result of the earlier treatment and a reason for the next one.
A covariate measured during follow-up that affects later treatment and the outcome is a time-varying confounder. When earlier treatment also changes it, as here, the pattern is called treatment-confounder feedback [1]. Patients treated at one decision also tend to be treated at the next ($A_{t-1} \to A_t$, where $A_{t-1}$ is the treatment at the previous decision). Frailty raises both the chance of organ dysfunction and the risk of death ($U \to L_t$ and $U \to Y$).
Frailty never acts on treatment except through organ dysfunction. That detail decides whether the problem can be solved with measured data, and it returns in the section on assumptions.
Why ordinary regression fails either way
Suppose the analysis is a regression of death on treatment, and the only decision is which covariates to include. Organ dysfunction can be left out or put in, and each choice breaks a different part of the diagram.
Leave organ dysfunction out: every dose stays confounded
Without $L_t$ in the model, the path $A_t \leftarrow L_t \to Y$ stays open on every decision day, including $A_0 \leftarrow L_0 \to Y$ on day 0. Patients treated on a decision day were treated because they had severe organ dysfunction, which raises the risk of death whatever the drug does. The drug inherits the risk of the patients who receive it and looks harmful. $L_0$, recorded before any treatment, is an ordinary baseline confounder that could be adjusted for safely; the trouble starts with $L_1$ and $L_2$, which earlier doses have already changed.
Adjust for organ dysfunction: the effect is blocked and a collider path opens
With $L_t$ in the model, two different things go wrong. First, $L_1$ lies on the path $A_0 \to L_1 \to Y$, so holding it fixed removes the part of the day-0 effect that works by preventing organ dysfunction. In the simulated cohort that is most of the effect. Of the true -0.050 per treated day, -0.045 runs through organ dysfunction and only -0.005 is direct.
Second, $L_1$ is a collider on another path. A collider is a variable that two arrows on a path point into, here $A_0 \to L_1 \leftarrow U$. Left alone, a collider blocks the path through it. Adjusting for it opens the path $A_0 \to L_1 \leftarrow U \to Y$, which links the day-0 dose to death through frailty.
The opened path has a direction here. Among patients with severe organ dysfunction on day 2, those treated on day 0 had it despite the drug, so more of them are frail. Frail patients die more often, so the adjusted model makes the earlier dose look worse than it is. The blocked effect and the opened path together leave the adjusted estimate far closer to zero than the truth.
The analyst is caught. Leave $L_t$ out and every dose is confounded; put it in and the earlier doses lose their real effect and gain a spurious one. No set of covariates in a single outcome regression solves both problems [1, 2].
The dilemma in one table
| Choice | What it fixes | What it breaks | Direction in the simulated cohort |
|---|---|---|---|
| Leave $L_t$ out | Keeps the effect that runs through organ dysfunction | Leaves $A_t \leftarrow L_t \to Y$ open, so each dose, including the first, is confounded by that day's organ dysfunction | Drug looks harmful |
| Adjust for $L_t$ | Removes the confounding of each dose by $L_t$ | Blocks the effect through $L_t$ and opens $A_{t-1} \to L_t \leftarrow U \to Y$ | Estimate pulled toward zero |
| Weight by the inverse probability of the treatment received (MSM) | Removes confounding by $L_t$ without conditioning on it | Needs no unmeasured confounding given the measured history, correct treatment models, and a chance above zero of each treatment for every history | Close to the truth |
Potential outcomes under a regime, and a model for them
A potential outcome is the outcome a patient would have had under a given treatment. Only the one for the treatment actually received is ever observed. Under a regime $\bar a$, write it $Y^{\bar a}$: whether the patient would have died by day 28 had they followed $\bar a$. With three decision days there are eight possible regimes, so each patient has eight potential outcomes, and here exactly one is observed, the one for the regime the patient actually followed.
A marginal structural model (MSM) is a model for the mean of these potential outcomes across regimes [2]. It is marginal because it averages over the time-varying confounders rather than conditioning on them; it may still include baseline covariates such as $V$. It is structural because it models potential outcomes, which are causal, rather than observed associations. A simple MSM for the ICU cohort is
$$E[Y^{\bar a}] = g_0 + g_1\,\mathrm{cum}(\bar a)$$Here $\mathrm{cum}(\bar a) = a_0 + a_1 + a_2$ is the number of treated days in the regime, where a treated day is a decision day on which drug X is given. In the model, $g_0$ is the risk of death if never treated, and $g_1$ is the change in that risk per treated day. The always-versus-never risk difference is then $3 g_1$. Because the model is linear in the risk, $g_1$ is itself a risk difference.
The model is an assumption: only the number of treated days matters, not their order. The simulated cohort was built so that this holds, with $g_1$ = -0.050. With real data, the form is usually chosen in advance and checked against a more flexible version, for example one term per decision day.
MSMs were introduced to epidemiology for exactly this problem [2], and an early application estimated the effect of zidovudine on the survival of HIV-positive men [3]. They are fitted after reweighting, so that in the weighted cohort treatment on each day no longer depends on the organ dysfunction history, given the earlier doses and age. The hand example below shows the reweighting with numbers small enough to check.
Hand example: a drug that works only by preventing organ dysfunction
Take 800 patients and two decision days, with no frailty and no age group. Half, chosen at random, receive drug X on the first day ($A_0 = 1$). Severe organ dysfunction on the second day ($L_1 = 1$) occurs in 50% of the untreated and 25% of the treated. On the second day the drug goes to 75% of patients with organ dysfunction and 25% of the rest ($A_1 = 1$).
The risk of death is 0.30 with organ dysfunction on the second day and 0.10 without, whatever the treatment. The drug therefore acts only by preventing organ dysfunction. Counts below are expected counts, so fractions of a death appear.
-
True risk if always treated
\[ \tfrac{1}{4} \times 0.30 + \tfrac{3}{4} \times 0.10 = 0.15 \]
If everyone were treated on the first day, a quarter would have organ dysfunction on the second.
-
True risk if never treated
\[ \tfrac{1}{2} \times 0.30 + \tfrac{1}{2} \times 0.10 = 0.20 \]
If no one were treated, half would. The true risk difference is 0.15 minus 0.20, or -0.05: the drug saves lives.
-
Who followed always treated
\[ 100 \times \tfrac{3}{4} + 300 \times \tfrac{1}{4} = 75 + 75 = 150 \]
Of the 400 patients treated on the first day, 100 have organ dysfunction on the second and 300 do not. Organ dysfunction makes the second dose likely, so the always-treated group holds 75 patients with organ dysfunction and 75 without.
-
Who followed never treated
\[ 200 \times \tfrac{1}{4} + 200 \times \tfrac{3}{4} = 50 + 150 = 200 \]
Of the untreated, 200 have organ dysfunction and 200 do not. The never-treated group holds 50 patients with organ dysfunction and 150 without.
-
Crude comparison
\[ \frac{22.5 + 7.5}{150} = 0.20, \; \frac{15 + 15}{200} = 0.15 \]
Deaths are 30 of 150 among the always treated and 30 of 200 among the never treated. The crude difference is +0.05, so the drug looks harmful: the always-treated group holds a larger share of patients with organ dysfunction.
-
Adjusting for organ dysfunction
Within each level of $L_1$, the risk is 0.30 or 0.10 whatever the treatment, so the adjusted difference is 0. The benefit runs entirely through organ dysfunction, and adjustment erases it.
-
Weight, always treated with organ dysfunction
\[ w = \frac{1}{0.5 \times 0.75} = \tfrac{8}{3} = 2.67 \]
A weight is one over the probability of the treatment history received, given the measured history. Here the first day contributes 1/0.5 = 2 and the second 1/0.75 = 1.33, a history of probability 3/8.
-
Weight, always treated without organ dysfunction
\[ w = \frac{1}{0.5 \times 0.25} = 8 \]
Without organ dysfunction, the second day contributes 1/0.25 = 4, a history of probability 1/8.
-
The weighted always-treated group
\[ \frac{60 + 60}{200 + 600} = \frac{120}{800} = 0.15 \]
The 75 always-treated patients with organ dysfunction count as 200, and the 75 without count as 600. The weighted group of 800 is a quarter with organ dysfunction, as if everyone had been always treated, and its risk is the true 0.15.
-
The weighted never-treated group
\[ \frac{120 + 40}{400 + 400} = \frac{160}{800} = 0.20 \]
Never-treated patients with organ dysfunction followed a history of probability 1/8 and get weight 8; the others get 2.67. The 50 and 150 patients count as 400 each, the group is half with organ dysfunction, and its risk is the true 0.20. The weighted difference is -0.05.
-
Stabilised weights, always treated
\[ sw = 0.1875 \times \tfrac{8}{3} = 0.5, \; sw = 0.1875 \times 8 = 1.5 \]
A stabilised weight multiplies by the probability of the same history ignoring organ dysfunction. That is 0.5 for the first dose times 0.375 for the second, since 150 of the 400 patients treated first were treated again. The numerator, 0.1875, shrinks the weights to 0.5 and 1.5.
-
Stabilised, weighted risks
\[ \frac{11.25 + 11.25}{37.5 + 112.5} = \frac{22.5}{150} = 0.15 \]
The stabilised always-treated group keeps its real size, 150, and the risk 0.15. For the never treated the numerator is 0.5 for no first dose times 0.5 for no second dose, because half of the patients untreated on the first day, 200 in all, stayed untreated. That gives 0.25, the weights become 2 and 0.67, the group keeps its size of 200, and the risk is 40 of 200, or 0.20.
Result: Weighting recovers the true risks, 0.15 and 0.20, and the true difference, -0.05, where the crude comparison gave +0.05 and adjustment gave 0. Stabilising leaves the answer unchanged and pulls the weights toward 1, from 2.67 and 8 to between 0.5 and 2.
Steps 1 and 2 are a small g-formula calculation: the risk under each regime, averaged over the organ dysfunction that regime would produce. The g-formula returns near the end of the article.
Hand example: the truth and four ways to estimate it
| Method | Risk, always treated | Risk, never treated | Difference |
|---|---|---|---|
| Truth | 0.15 | 0.20 | -0.05 |
| Crude comparison | 0.20 | 0.15 | +0.05 |
| Adjusted for $L_1$ | 0.30 or 0.10, by level of $L_1$ | 0.30 or 0.10, by level of $L_1$ | 0 |
| Inverse probability weights | 0.15 | 0.20 | -0.05 |
| Stabilised weights | 0.15 | 0.20 | -0.05 |
Cumulative stabilised weights
With three decision days, each weight is a product of one factor per day. For patient $i$, with observed treatments $a_{i0}$, $a_{i1}$ and $a_{i2}$, the cumulative stabilised weight is
$$sw_i = \prod_{t=0}^{2} \frac{P(A_t = a_{it} \mid \bar A_{t-1}, V)}{P(A_t = a_{it} \mid \bar A_{t-1}, \bar L_t, V)}$$The denominator is the probability of the treatment actually received on day $t$, given everything measured before that decision: earlier treatment $\bar A_{t-1}$ (nothing before day 0), the organ dysfunction history $\bar L_t$ and the age group $V$. The numerator is the same probability without organ dysfunction. Dropping the numerator gives the unstabilised weight $w_i$, one over the product of the denominators.
In the weighted cohort, treatment on each day no longer depends on organ dysfunction, given the earlier doses and age, as if it had been assigned without looking at it. Organ dysfunction still responds to earlier doses as before, so the effect that runs through it is kept. Confounding by $L_t$ is removed without conditioning on $L_t$, which is what regression could not do. This weighted cohort is often called a pseudo-population.
Very large weights have two common remedies besides stabilising. Truncation caps the largest weights at a chosen value, such as a high percentile. Trimming drops patients whose treatment probabilities are extreme, close to 0 or 1. All three are worked through at a single time point in Extreme Weights: What Stabilizing, Truncating and Trimming Really Change.
A point-treatment analysis has one treatment decision per patient. In a point-treatment analysis, stabilising multiplies every weight within an arm by the same constant, so an extreme patient stays extreme relative to everyone else in that arm. Stabilisation matters most in marginal structural models, where weights are multiplied over many time points. Truncation trades bias for variance; trimming changes the population the estimate describes.
In an MSM the numerator differs between treatment histories, so stabilising rescales whole histories relative to one another rather than every weight by one constant. Within one treatment history and age group the factor is the same, as in the hand example, where the always-treated weights 2.67 and 8 became 0.5 and 1.5. Unstabilised weights give every treatment history an equal share of the weighted cohort and stabilised weights keep the observed mix, so when the MSM's form is wrong the two can give different answers. Over three days the gain is already visible in the simulated cohort below, where the stabilised weights carry the information of 1,694 equally weighted patients and the unstabilised ones that of 1,456.
What may go in the numerator
The numerator may hold earlier treatment and baseline covariates, but never $L_t$ or anything else measured after treatment starts. Putting $L_t$ in the numerator would cancel it from the weight and bring the confounding back. Any baseline covariate in the numerator, here $V$, must also appear in the MSM [4]. The model then becomes
$$E[Y^{\bar a} \mid V] = g_0 + g_1\,\mathrm{cum}(\bar a) + g_2 V$$Here $g_0$ becomes the risk if never treated for patients younger than 65 ($V = 0$), and $g_2$ is the difference in risk between the age groups under any regime. In the simulated cohort a treated day has the same effect in both age groups, so $g_1$ means the same with or without $V$. To read the risk under each regime for the whole cohort, the Stata and R code further down also fits a marginal version, whose numerator drops $V$ and whose MSM holds only $\mathrm{cum}(\bar a)$.
Fitting the MSM, step by step
- Arrange the data with one row per patient per decision day, holding $A_t$, $L_t$, $V$ and the previous treatment $A_{t-1}$ (0 on day 0).
- Fit the denominator model: a logistic regression of $A_t$ on $L_t$, $A_{t-1}$, $V$ and the day, pooled over all patient-days. Pooled means one model for all days, with a term for the day.
- Fit the numerator model the same way, without $L_t$.
- On each row, take the predicted probability of treatment, $p$. The probability of the treatment actually received is $p$ if treated and $1 - p$ if not.
- Within each patient, multiply the numerator-to-denominator ratios over the three days.
- Fit the outcome model on one row per patient, weighted by the cumulative weight, with a standard error that allows for the weights.
In Stata the treatment models are logit A L A_prev V i.t and logit A A_prev V i.t. The cumulative product is a running sum of logs: by id (t): generate double sw = exp(sum(ln(fnum / fden))).
In R the denominator is den <- glm(A ~ L + A_prev + V + factor(t), data = d, family = binomial), the numerator drops L, and ave() with FUN = cumprod takes the product within each patient. Worked Stata code for this approach is published [5], and the R package ipw builds the same weights [6].
The outcome model here is a weighted linear regression of death on $\mathrm{cum}(\bar a)$ and $V$. It is a linear risk model, whose coefficients are risk differences. In Stata it is glm Y cumA V if t == 2 [pweight = sw], family(gaussian) link(identity) vce(cluster id), where t == 2 keeps each patient's day-4 row. A logistic MSM would also be valid if correctly specified and would give an odds ratio per treated day; the linear form suits this cohort because the true risk is linear by design.
A robust (sandwich) standard error estimates the variance from the spread of each patient's contribution rather than from model formulas, so it allows for unequal weights. It treats the estimated weights as if they were known, which usually makes the interval somewhat conservative [2].
With one row per patient, as here, the clustered and the ordinary robust standard errors coincide; clustering matters when the outcome model has a row per patient-day. A bootstrap, which resamples patients and repeats the whole analysis including both weight models, is the alternative when the interval matters.
In the simulated cohort, the decision on each day depends only on that day's organ dysfunction, the previous dose and age, so those terms are enough. With real data, the denominator would normally include whatever earlier history plausibly drives treatment, such as organ dysfunction at earlier reviews.
A simulated ICU cohort with a known answer
The cohort holds 3,000 simulated patients seen on days 0, 2 and 4, giving 9,000 patient-day rows. Every patient survives to day 4, so all three decisions are made, and death by day 28 is recorded for everyone. Every probability that generated the data was fixed in advance, so the true effect is known. The table below lists the design, and the next two paragraphs read it against the diagram.
Each dose lowers the chance of severe organ dysfunction at the next decision by 0.30, and each day with organ dysfunction adds 0.15 to the risk of death. A treated day therefore saves $0.30 \times 0.15 = 0.045$ through organ dysfunction; the last dose acts through organ dysfunction on day 6, which is never recorded. Adding the direct part, -0.005, gives the true effect: $g_1$ = -0.050 per treated day.
Because death risk is a straight-line function of its causes, the linear MSM is exactly right. The true risk of death by day 28 is 0.4725 if never treated and 0.3225 if always treated, a difference of -0.150. Severe organ dysfunction drives treatment strongly, with a coefficient of 2.0 on the log-odds scale, so the treated are much sicker on every decision day. Of the 3,000 patients, 1,232 died (41.1%), and 4,286 of the 9,000 patient-days were treated.
How the simulated cohort was generated
| Variable | Meaning | True model |
|---|---|---|
| $U$ | Unmeasured frailty, never recorded | 1 with probability 0.40 |
| $V$ | Aged 65 or older, measured at baseline | 1 with probability 0.45 |
| $L_t$ | Severe organ dysfunction on decision day $t$ | $P(L_t = 1) = 0.32 + 0.45\,U - 0.30\,A_{t-1}$ |
| $A_t$ | Drug X on decision day $t$ | $\operatorname{logit} P(A_t = 1) = -1.2 + 2.0\,L_t + 1.5\,A_{t-1} - 0.5\,V$ |
| $Y$ | Death by day 28 | $P(Y = 1) = 0.05 + 0.25\,U + 0.05\,V + 0.15\,(L_0 + L_1 + L_2 + L_6) - 0.005\,(A_0 + A_1 + A_2)$ |
What each analysis estimates
Six analyses were run in both Stata and R, and a seventh in R only. Four are ordinary regressions that a team might try first, and three are MSMs fitted with cumulative stabilised weights. The table gives each estimate of the effect of one treated day beside its large-sample value: where the same analysis would land in an infinitely large cohort.
| Analysis | Estimate (risk difference) | 95% CI | Large-sample value |
|---|---|---|---|
| Pooled patient-days, unadjusted | 0.086 | 0.063 to 0.109 | 0.086 |
| Pooled patient-days, adjusted for $L_t$ and $V$ | -0.013 | -0.036 to 0.011 | -0.012 |
| Number of treated days, unadjusted | 0.068 | 0.050 to 0.085 | 0.064 |
| Number of treated days, adjusted for $V$, $L_0$, $L_1$ and $L_2$ | 0.011 | -0.006 to 0.027 | 0.0005 |
| MSM, stabilised weights, with $V$ | -0.057 | -0.081 to -0.034 | -0.050 (true) |
| Marginal MSM, stabilised weights | -0.057 | -0.080 to -0.033 | -0.050 (true) |
| MSM, one treatment model per day (R, WeightIt) | -0.055 | -0.079 to -0.031 | -0.050 (true) |
Reading the table
Pooling patient-days without adjustment lands furthest from the truth: each treated day looks like +0.086 in the risk of death. Adjusting the pooled model for $L_t$ and $V$ swings the estimate to -0.013, and its large-sample value, -0.012, is a small part of the true benefit. The models that count treated days fare no better. Unadjusted, the drug looks harmful at 0.068 per day; adjusted for organ dysfunction on days 0, 2 and 4, the large-sample value is 0.0005, essentially zero.
The fellow's two regressions are the third and fourth rows. Their large-sample values show that more patients would not have helped, because the bias is built into the question each regression answers.
The MSM with stabilised weights gives -0.057 per treated day (95% CI -0.081 to -0.034), and the interval covers the true -0.050. Three treated days give an always-versus-never risk difference of -0.172 (95% CI -0.242 to -0.102), against a true -0.150. The marginal MSM gives -0.057, and weights from one treatment model per day, built with the WeightIt package in R, give -0.055.
The marginal MSM also returns the risk under each regime: 0.489 if never treated and 0.319 if always treated, against true risks of 0.4725 and 0.3225. Under the true model, one treated day gives 0.4225 and two give 0.3725, a fall of 0.050 for each day. Through its straight line, the MSM uses all 3,000 patients, not only the 471 and 588 who followed the two regimes.
Risk of death by day 28 under each regime
| Regime | True risk | Marginal MSM estimate |
|---|---|---|
| Never treated, $\bar a = (0, 0, 0)$ | 0.4725 | 0.489 |
| Always treated, $\bar a = (1, 1, 1)$ | 0.3225 | 0.319 |
What the weights look like
Weights deserve a look before the outcome model is read. The table below compares the cumulative weights at day 4, after all three decisions. The effective sample size (ESS) is the number of equally weighted patients that would carry the same information.
| Weight | Mean | Maximum | Effective sample size |
|---|---|---|---|
| Stabilised, $sw$ | 0.999 | 10.30 | 1,694 |
| Unstabilised, $w$ | 7.96 | 94.94 | 1,456 |
Reading the weights
The stabilised weights average 0.999, as they should: a correctly modelled stabilised weight has an expected value of 1 [4]. The smallest is 0.27, 99% are at or below 4.03, and the largest is 10.30. The stabilised weights of the marginal MSM behave the same way, with a mean of 0.995 and a maximum of 9.52.
The unstabilised weights average 7.96, close to eight, the number of possible treatment histories over three days, and one patient's weight reaches 94.94. Relative to its own mean, that largest weight is about as extreme as the largest stabilised weight is relative to 1, so comparing the two maxima overstates what stabilising gains.
The ESS is computed as $(\sum w)^2 / \sum w^2$; it does not change when every weight is multiplied by the same constant. It is 1,456 for the unstabilised weights against 1,694 for the stabilised ones, so with unstabilised weights the same cohort carries less information.
With daily decisions over weeks the products grow longer and unstabilised weights grow wilder, which is why stabilised weights are the usual choice for MSMs [4].
Stata: the weight models, the cumulative weights and the MSM
* Simulated ICU cohort, three treatment days: treatment-confounder feedback.
* Simulated data: not evidence about any real drug or patient.
* Question: what does each treated day of drug X do to the risk of death by day 28?
* Methods: naive regressions (with and without L_t) versus a marginal structural model (MSM)
* fitted with cumulative stabilised inverse probability of treatment weights.
* Variables: L_t = 1 when severe organ dysfunction is present on day t, A_t = 1 when the drug is given on day t,
* V = 1 for age 65 or older, Y = 1 for death by day 28.
version 18
clear all
set more off
* Check the built-in commands this file uses.
capture which logit
display "VERIFY logit " cond(_rc == 0, "available", "missing")
capture which glm
display "VERIFY glm " cond(_rc == 0, "available", "missing")
* Read the simulated cohort file (not published): one row per patient-day (days 0, 2, 4).
import delimited using ../../datasets/W2/W2.csv, clear case(preserve) varnames(1)
sort id t
count if t == 0
display "CANON w2.n " strtrim(string(r(N), "%12.0f"))
count
display "CANON w2.n_rows " strtrim(string(r(N), "%12.0f"))
count if t == 2 & Y == 1
display "CANON w2.deaths " strtrim(string(r(N), "%12.0f"))
count if A == 1
display "CANON w2.treated_days " strtrim(string(r(N), "%12.0f"))
* Naive 1: pool all patient-days and regress death on that day's treatment (linear risk model), CI clustered by patient.
glm Y A, family(gaussian) link(identity) vce(cluster id)
display "CANON w2.naive.unadj.b " strtrim(string(_b[A], "%12.4f"))
display "CANON w2.naive.unadj.se " strtrim(string(_se[A], "%12.4f"))
display "CANON w2.naive.unadj.lo " strtrim(string(_b[A] - invnormal(0.975) * _se[A], "%12.4f"))
display "CANON w2.naive.unadj.hi " strtrim(string(_b[A] + invnormal(0.975) * _se[A], "%12.4f"))
* Naive 2: the same pooled regression adjusted for that day's organ dysfunction L_t and age V.
glm Y A L V, family(gaussian) link(identity) vce(cluster id)
display "CANON w2.naive.adj.b " strtrim(string(_b[A], "%12.4f"))
display "CANON w2.naive.adj.se " strtrim(string(_se[A], "%12.4f"))
display "CANON w2.naive.adj.lo " strtrim(string(_b[A] - invnormal(0.975) * _se[A], "%12.4f"))
display "CANON w2.naive.adj.hi " strtrim(string(_b[A] + invnormal(0.975) * _se[A], "%12.4f"))
* Cumulative treatment through each day, and organ dysfunction on days 0, 2 and 4 copied to every row.
by id (t): generate cumA = sum(A)
by id (t): generate L0 = L[1]
by id (t): generate L1 = L[2]
by id (t): generate L2 = L[3]
count if t == 2 & cumA == 3
display "CANON w2.always_treated " strtrim(string(r(N), "%12.0f"))
count if t == 2 & cumA == 0
display "CANON w2.never_treated " strtrim(string(r(N), "%12.0f"))
* Naive 3: one row per patient (day 4), death on the number of treated days (no adjustment).
glm Y cumA if t == 2, family(gaussian) link(identity) vce(cluster id)
display "CANON w2.cum.crude.b " strtrim(string(_b[cumA], "%12.4f"))
display "CANON w2.cum.crude.se " strtrim(string(_se[cumA], "%12.4f"))
display "CANON w2.cum.crude.lo " strtrim(string(_b[cumA] - invnormal(0.975) * _se[cumA], "%12.4f"))
display "CANON w2.cum.crude.hi " strtrim(string(_b[cumA] + invnormal(0.975) * _se[cumA], "%12.4f"))
* Naive 4: the same, adjusted for age and organ dysfunction on days 0, 2 and 4.
glm Y cumA V L0 L1 L2 if t == 2, family(gaussian) link(identity) vce(cluster id)
display "CANON w2.cum.adj.b " strtrim(string(_b[cumA], "%12.4f"))
display "CANON w2.cum.adj.se " strtrim(string(_se[cumA], "%12.4f"))
display "CANON w2.cum.adj.lo " strtrim(string(_b[cumA] - invnormal(0.975) * _se[cumA], "%12.4f"))
display "CANON w2.cum.adj.hi " strtrim(string(_b[cumA] + invnormal(0.975) * _se[cumA], "%12.4f"))
* Treatment models, pooled over days: the denominator uses history (L_t, A_{t-1}, V).
logit A L A_prev V i.t
predict double pden, pr
* Numerator for the stabilised weight: same model without L_t.
logit A A_prev V i.t
predict double pnum, pr
* Numerator for the marginal MSM: without L_t and without V.
logit A A_prev i.t
predict double pnumm, pr
* Probability of the treatment each patient actually received on each day.
generate double fden = cond(A == 1, pden, 1 - pden)
generate double fnum = cond(A == 1, pnum, 1 - pnum)
generate double fnumm = cond(A == 1, pnumm, 1 - pnumm)
* Cumulative product over days within each patient (as a running sum of logs) gives the weight at day 4.
by id (t): generate double w = exp(sum(ln(1 / fden)))
by id (t): generate double sw = exp(sum(ln(fnum / fden)))
by id (t): generate double swm = exp(sum(ln(fnumm / fden)))
* Weight summaries at day 4; the effective sample size is (sum of weights)^2 / sum of squared weights.
summarize sw if t == 2, detail
display "CANON w2.sw.mean " strtrim(string(r(mean), "%12.4f"))
display "CANON w2.sw.sd " strtrim(string(r(sd), "%12.4f"))
display "CANON w2.sw.min " strtrim(string(r(min), "%12.4f"))
display "CANON w2.sw.p1 " strtrim(string(r(p1), "%12.4f"))
display "CANON w2.sw.p99 " strtrim(string(r(p99), "%12.4f"))
display "CANON w2.sw.max " strtrim(string(r(max), "%12.4f"))
generate double sw_sq = sw^2
summarize sw if t == 2
scalar sum_sw = r(sum)
summarize sw_sq if t == 2
display "CANON w2.sw.ess " strtrim(string(sum_sw^2 / r(sum), "%12.4f"))
generate double w_sq = w^2
summarize w if t == 2
display "CANON w2.w.mean " strtrim(string(r(mean), "%12.4f"))
display "CANON w2.w.max " strtrim(string(r(max), "%12.4f"))
scalar sum_w = r(sum)
summarize w_sq if t == 2
display "CANON w2.w.ess " strtrim(string(sum_w^2 / r(sum), "%12.4f"))
summarize swm if t == 2
display "CANON w2.swm.mean " strtrim(string(r(mean), "%12.4f"))
display "CANON w2.swm.max " strtrim(string(r(max), "%12.4f"))
* MSM: weighted regression of death on cumulative treatment and age, CI clustered by patient.
glm Y cumA V if t == 2 [pweight = sw], family(gaussian) link(identity) vce(cluster id)
display "CANON w2.msm.b " strtrim(string(_b[cumA], "%12.4f"))
display "CANON w2.msm.se " strtrim(string(_se[cumA], "%12.4f"))
display "CANON w2.msm.lo " strtrim(string(_b[cumA] - invnormal(0.975) * _se[cumA], "%12.4f"))
display "CANON w2.msm.hi " strtrim(string(_b[cumA] + invnormal(0.975) * _se[cumA], "%12.4f"))
* Always treated versus never treated is three treated days: three times the per-day effect.
display "CANON w2.msm.rd_always_never " strtrim(string(3 * _b[cumA], "%12.4f"))
display "CANON w2.msm.rd_always_never.lo " strtrim(string(3 * (_b[cumA] - invnormal(0.975) * _se[cumA]), "%12.4f"))
display "CANON w2.msm.rd_always_never.hi " strtrim(string(3 * (_b[cumA] + invnormal(0.975) * _se[cumA]), "%12.4f"))
* Marginal MSM (numerator without V, no V in the model): gives the risks under never and always treated.
glm Y cumA if t == 2 [pweight = swm], family(gaussian) link(identity) vce(cluster id)
display "CANON w2.msm_marg.b " strtrim(string(_b[cumA], "%12.4f"))
display "CANON w2.msm_marg.se " strtrim(string(_se[cumA], "%12.4f"))
display "CANON w2.msm_marg.lo " strtrim(string(_b[cumA] - invnormal(0.975) * _se[cumA], "%12.4f"))
display "CANON w2.msm_marg.hi " strtrim(string(_b[cumA] + invnormal(0.975) * _se[cumA], "%12.4f"))
display "CANON w2.msm_marg.risk_never " strtrim(string(_b[_cons], "%12.4f"))
display "CANON w2.msm_marg.risk_always " strtrim(string(_b[_cons] + 3 * _b[cumA], "%12.4f"))
display "Simulated data (ข้อมูลจำลอง): no number above is evidence about any real drug."
* ---- The true values built into the simulated cohort (exact: computed from the design, no sampling) ----
* The cohort was simulated from these probabilities. U is an unmeasured patient trait; L_6 is the day-6
* value of L (after the last treatment day, never recorded).
* P(U = 1) = pU and P(V = 1) = pV
* P(L_t = 1) = lBase + lU*U + lA*A_{t-1}, with A_{t-1} = 0 on day 0
* logit P(A_t = 1) = aBase + aL*L_t + aA*A_{t-1} + aV*V
* P(death by day 28) = yBase + yU*U + yV*V + yL*(L_0 + L_1 + L_2 + L_6) + yA*(A_0 + A_1 + A_2)
* The true values need no patient data, so start from an empty dataset.
clear
scalar pU = 0.40
scalar pV = 0.45
scalar lBase = 0.32
scalar lU = 0.45
scalar lA = -0.30
scalar aBase = -1.2
scalar aL = 2.0
scalar aA = 1.5
scalar aV = -0.5
scalar yBase = 0.05
scalar yU = 0.25
scalar yV = 0.05
scalar yL = 0.15
scalar yA = -0.005
display "CANON w2.truth.design.p_u " strtrim(string(pU, "%12.4f"))
display "CANON w2.truth.design.p_v " strtrim(string(pV, "%12.4f"))
display "CANON w2.truth.design.l_base " strtrim(string(lBase, "%12.4f"))
display "CANON w2.truth.design.l_u " strtrim(string(lU, "%12.4f"))
display "CANON w2.truth.design.l_a " strtrim(string(lA, "%12.4f"))
display "CANON w2.truth.design.a_base " strtrim(string(aBase, "%12.4f"))
display "CANON w2.truth.design.a_l " strtrim(string(aL, "%12.4f"))
display "CANON w2.truth.design.a_a " strtrim(string(aA, "%12.4f"))
display "CANON w2.truth.design.a_v " strtrim(string(aV, "%12.4f"))
display "CANON w2.truth.design.y_base " strtrim(string(yBase, "%12.4f"))
display "CANON w2.truth.design.y_u " strtrim(string(yU, "%12.4f"))
display "CANON w2.truth.design.y_v " strtrim(string(yV, "%12.4f"))
display "CANON w2.truth.design.y_l " strtrim(string(yL, "%12.4f"))
display "CANON w2.truth.design.y_a " strtrim(string(yA, "%12.4f"))
* Each treated day lowers the risk of death directly (yA) and by lowering the next day's L (yL x lA).
scalar g1 = yA + yL * lA
display "CANON w2.truth.msm.g1_per_treated_day " strtrim(string(g1, "%12.4f"))
display "CANON w2.truth.msm.g1_direct_part " strtrim(string(yA, "%12.4f"))
display "CANON w2.truth.msm.g1_through_l_part " strtrim(string(yL * lA, "%12.4f"))
* Risk of death by day 28 if every patient were treated on k of the three days (k = 0 never, k = 3 always).
* EL is P(L_t = 1) after an untreated day; L_0, L_1, L_2 and L_6 each add yL x EL, and each treated day adds g1.
scalar EL = lBase + lU * pU
scalar rNever = yBase + yU * pU + yV * pV + yL * 4 * EL
display "CANON w2.truth.msm.risk_never_treated " strtrim(string(rNever, "%12.4f"))
display "CANON w2.truth.msm.risk_cum1 " strtrim(string(rNever + g1, "%12.4f"))
display "CANON w2.truth.msm.risk_cum2 " strtrim(string(rNever + 2 * g1, "%12.4f"))
display "CANON w2.truth.msm.risk_always_treated " strtrim(string(rNever + 3 * g1, "%12.4f"))
display "CANON w2.truth.msm.rd_always_vs_never " strtrim(string(3 * g1, "%12.4f"))
* The model with age V (the one fitted above): intercept at V = 0 and the coefficient of V.
display "CANON w2.truth.msm.g0_conditional " strtrim(string(yBase + yU * pU + yL * 4 * EL, "%12.4f"))
display "CANON w2.truth.msm.g2_v " strtrim(string(yV, "%12.4f"))
* Where each naive regression lands in an infinitely large cohort: list all 512 combinations of
* U, V, L_0, A_0, L_1, A_1, L_2, A_2 and L_6 with their probabilities under the observed treatment
* process, then fit the same regressions to the expected risk, weighting each combination by its probability.
set obs 512
generate int cell = _n - 1
generate byte U = mod(cell, 2)
generate byte V = mod(floor(cell / 2), 2)
generate byte L0 = mod(floor(cell / 4), 2)
generate byte A0 = mod(floor(cell / 8), 2)
generate byte L1 = mod(floor(cell / 16), 2)
generate byte A1 = mod(floor(cell / 32), 2)
generate byte L2 = mod(floor(cell / 64), 2)
generate byte A2 = mod(floor(cell / 128), 2)
generate byte L6 = mod(floor(cell / 256), 2)
generate double pL0 = lBase + lU * U
generate double pA0 = invlogit(aBase + aL * L0 + aV * V)
generate double pL1 = lBase + lU * U + lA * A0
generate double pA1 = invlogit(aBase + aL * L1 + aA * A0 + aV * V)
generate double pL2 = lBase + lU * U + lA * A1
generate double pA2 = invlogit(aBase + aL * L2 + aA * A1 + aV * V)
generate double pL6 = lBase + lU * U + lA * A2
generate double p = cond(U == 1, pU, 1 - pU) * cond(V == 1, pV, 1 - pV)
replace p = p * cond(L0 == 1, pL0, 1 - pL0) * cond(A0 == 1, pA0, 1 - pA0)
replace p = p * cond(L1 == 1, pL1, 1 - pL1) * cond(A1 == 1, pA1, 1 - pA1)
replace p = p * cond(L2 == 1, pL2, 1 - pL2) * cond(A2 == 1, pA2, 1 - pA2)
replace p = p * cond(L6 == 1, pL6, 1 - pL6)
generate double risk = yBase + yU * U + yV * V + yL * (L0 + L1 + L2 + L6) + yA * (A0 + A1 + A2)
generate cumA = A0 + A1 + A2
* Naive 3 and 4 (one row per patient).
quietly regress risk cumA [aweight = p]
display "CANON w2.truth.naive_limit.cum_crude " strtrim(string(_b[cumA], "%12.4f"))
quietly regress risk cumA V L0 L1 L2 [aweight = p]
display "CANON w2.truth.naive_limit.cum_adj " strtrim(string(_b[cumA], "%12.4f"))
* Naive 1 and 2 (one row per patient-day): three rows per combination, one per treatment day.
expand 3
bysort cell: generate t = _n - 1
generate A = cond(t == 0, A0, cond(t == 1, A1, A2))
generate L = cond(t == 0, L0, cond(t == 1, L1, L2))
quietly regress risk A [aweight = p]
display "CANON w2.truth.naive_limit.unadj " strtrim(string(_b[A], "%12.4f"))
quietly regress risk A L V [aweight = p]
display "CANON w2.truth.naive_limit.adj " strtrim(string(_b[A], "%12.4f"))
display "These true values describe the simulation design (simulated data), not any real drug."
* ---- Follow-up length, the marginal MSM intercept, and how many days each patient was treated ----
* Death is counted up to day 28.
scalar deathDay = 28
display "CANON w2.truth.design.death_day " strtrim(string(deathDay, "%12.0f"))
* The intercept of the marginal MSM is the risk of death by day 28 if no patient were ever treated.
display "CANON w2.truth.msm.g0_marginal " strtrim(string(rNever, "%12.4f"))
* Read the simulated cohort file (not published) again: share who died, patients treated on one or two days.
import delimited using ../../datasets/W2/W2.csv, clear case(preserve) varnames(1)
sort id t
by id (t): generate cumA = sum(A)
summarize Y if t == 2
display "CANON w2.death_risk " strtrim(string(r(mean), "%12.4f"))
count if t == 2 & cumA == 1
display "CANON w2.treated_one_day " strtrim(string(r(N), "%12.0f"))
count if t == 2 & cumA == 2
display "CANON w2.treated_two_days " strtrim(string(r(N), "%12.0f"))
tabulate cumA if t == 2
display "Simulated data (ข้อมูลจำลอง): not evidence about any real drug or patient.". * MSM: weighted regression of death on cumulative treatment and age, CI clust
> ered by patient.
. glm Y cumA V if t == 2 [pweight = sw], family(gaussian) link(identity) vce(cl
> uster id)
Iteration 0: Log pseudolikelihood = -2096.0075
Generalized linear models Number of obs = 3,000
Optimization : ML Residual df = 2,997
Scale parameter = .2370902
Deviance = 710.5593289 (1/df) Deviance = .2370902
Pearson = 710.5593289 (1/df) Pearson = .2370902
Variance function: V(u) = 1 [Gaussian]
Link function : g(u) = u [Identity]
AIC = 1.399338
Log pseudolikelihood = -2096.007537 BIC = -23284.52
(Std. err. adjusted for 3,000 clusters in id)
------------------------------------------------------------------------------
| Robust
Y | Coefficient std. err. z P>|z| [95% conf. interval]
-------------+----------------------------------------------------------------
cumA | -.0572908 .0118918 -4.82 0.000 -.0805982 -.0339834
V | .0477372 .0242604 1.97 0.049 .0001876 .0952868
_cons | .4654255 .0266901 17.44 0.000 .4131138 .5177371
------------------------------------------------------------------------------R: the same analyses and the summary tables
# Simulated ICU cohort, three treatment days: treatment-confounder feedback.
# Simulated data: not evidence about any real drug or patient.
# Question: what does each treated day of drug X do to the risk of death by day 28?
# Methods: naive regressions (with and without L_t) versus a marginal structural model (MSM)
# fitted with cumulative stabilised inverse probability of treatment weights.
# Variables: L_t = 1 when severe organ dysfunction is present on day t, A_t = 1 when the drug is given on day t,
# V = 1 for age 65 or older, Y = 1 for death by day 28.
# Check the packages this script uses; install into the user library when one is missing.
need <- c("sandwich", "ipw", "WeightIt")
for (p in need) {
if (!requireNamespace(p, quietly = TRUE)) {
try(install.packages(p, repos = "https://cloud.r-project.org"), silent = TRUE)
}
cat(sprintf("VERIFY %s %s\n", p, if (requireNamespace(p, quietly = TRUE)) "available" else "missing"))
}
# Print each result on its own labelled line: 4 decimals for estimates, integers for counts.
canon <- function(key, value, digits = 4) {
cat(sprintf("CANON w2.%s %s\n", key, formatC(value, format = "f", digits = digits)))
}
# Read the simulated cohort file (not published): one row per patient-day (days 0, 2, 4).
d <- read.csv(file.path("..", "..", "datasets", "W2", "W2.csv"))
d <- d[order(d$id, d$t), ]
canon("n", length(unique(d$id)), 0)
canon("n_rows", nrow(d), 0)
canon("deaths", sum(d$Y[d$t == 2]), 0)
canon("treated_days", sum(d$A), 0)
# Cluster-robust (sandwich) 95% CI by patient for one coefficient; HC0 with the G/(G-1) factor, as Stata glm does.
ci_cl <- function(fit, term, cluster) {
V <- sandwich::vcovCL(fit, cluster = cluster, type = "HC0", cadjust = TRUE)
b <- unname(coef(fit)[term]); se <- sqrt(V[term, term])
c(b = b, se = se, lo = b - qnorm(0.975) * se, hi = b + qnorm(0.975) * se)
}
print_est <- function(key, est) {
canon(paste0(key, ".b"), est[["b"]]); canon(paste0(key, ".se"), est[["se"]])
canon(paste0(key, ".lo"), est[["lo"]]); canon(paste0(key, ".hi"), est[["hi"]])
}
# Naive 1: pool all patient-days and regress death on that day's treatment (linear risk model).
f_unadj <- glm(Y ~ A, data = d, family = gaussian)
print_est("naive.unadj", ci_cl(f_unadj, "A", ~id))
# Naive 2: the same pooled regression adjusted for that day's organ dysfunction L_t and age V.
f_adj <- glm(Y ~ A + L + V, data = d, family = gaussian)
print_est("naive.adj", ci_cl(f_adj, "A", ~id))
# Cumulative treatment through each day, and organ dysfunction on each day, in wide form.
d$cumA <- ave(d$A, d$id, FUN = cumsum)
wide <- d[d$t == 2, c("id", "V", "Y", "cumA")]
Lw <- reshape(d[, c("id", "t", "L")], idvar = "id", timevar = "t", direction = "wide")
wide <- merge(wide, Lw, by = "id")
wide <- wide[order(wide$id), ]
canon("always_treated", sum(wide$cumA == 3), 0)
canon("never_treated", sum(wide$cumA == 0), 0)
# Naive 3: one row per patient, death on the number of treated days (no adjustment).
f_cc <- glm(Y ~ cumA, data = wide, family = gaussian)
print_est("cum.crude", ci_cl(f_cc, "cumA", ~id))
# Naive 4: the same, adjusted for age and organ dysfunction on days 0, 2 and 4.
f_ca <- glm(Y ~ cumA + V + L.0 + L.1 + L.2, data = wide, family = gaussian)
print_est("cum.adj", ci_cl(f_ca, "cumA", ~id))
# Treatment models, pooled over days: denominator uses history (L_t, A_{t-1}, V); numerators drop L_t.
den <- glm(A ~ L + A_prev + V + factor(t), data = d, family = binomial)
num <- glm(A ~ A_prev + V + factor(t), data = d, family = binomial)
numm <- glm(A ~ A_prev + factor(t), data = d, family = binomial)
pa <- function(fit) { p <- fitted(fit); ifelse(d$A == 1, p, 1 - p) }
# Cumulative product over days within each patient gives the weight at day 4.
d$w <- ave(1 / pa(den), d$id, FUN = cumprod)
d$sw <- ave(pa(num) / pa(den), d$id, FUN = cumprod)
d$swm <- ave(pa(numm) / pa(den), d$id, FUN = cumprod)
last <- d[d$t == 2, ]
last <- last[order(last$id), ]
# Weight summaries at day 4 (percentiles use the same averaging rule as Stata summarize, detail).
ess <- function(x) sum(x)^2 / sum(x^2)
q <- function(x, p) unname(quantile(x, p, type = 2))
canon("sw.mean", mean(last$sw)); canon("sw.sd", sd(last$sw))
canon("sw.min", min(last$sw)); canon("sw.p1", q(last$sw, 0.01))
canon("sw.p99", q(last$sw, 0.99)); canon("sw.max", max(last$sw))
canon("sw.ess", ess(last$sw))
canon("w.mean", mean(last$w)); canon("w.max", max(last$w)); canon("w.ess", ess(last$w))
canon("swm.mean", mean(last$swm)); canon("swm.max", max(last$swm))
# MSM: weighted regression of death on cumulative treatment and age, cluster-robust CI by patient.
msm <- glm(Y ~ cumA + V, data = last, weights = sw, family = gaussian)
e_msm <- ci_cl(msm, "cumA", ~id)
print_est("msm", e_msm)
canon("msm.rd_always_never", 3 * e_msm[["b"]])
canon("msm.rd_always_never.lo", 3 * e_msm[["lo"]])
canon("msm.rd_always_never.hi", 3 * e_msm[["hi"]])
# Marginal MSM (numerator without V, no V in the model): gives the risks under never and always treated.
msmm <- glm(Y ~ cumA, data = last, weights = swm, family = gaussian)
e_m <- ci_cl(msmm, "cumA", ~id)
print_est("msm_marg", e_m)
canon("msm_marg.risk_never", unname(coef(msmm)[1]))
canon("msm_marg.risk_always", unname(coef(msmm)[1] + 3 * coef(msmm)[2]))
# Cross-check: ipw::ipwtm builds the same pooled stabilised weights from its own code.
if (requireNamespace("ipw", quietly = TRUE)) {
dd <- d
dd$tf <- factor(dd$t)
iw <- ipw::ipwtm(exposure = A, family = "binomial", link = "logit",
numerator = ~ A_prev + V + tf, denominator = ~ L + A_prev + V + tf,
id = id, timevar = t, type = "all", data = dd)
cat(sprintf("CHECK ipwtm versus hand-built weights, largest absolute difference %.2e\n",
max(abs(iw$ipw.weights - dd$sw))))
}
# Cross-check: WeightIt::weightitMSM fits one treatment model per day (not pooled), so it differs slightly.
if (requireNamespace("WeightIt", quietly = TRUE)) {
wd <- reshape(d[, c("id", "t", "V", "L", "A", "Y")], idvar = c("id", "V", "Y"),
timevar = "t", direction = "wide")
wd <- wd[order(wd$id), ]
wm <- WeightIt::weightitMSM(list(A.0 ~ L.0 + V, A.1 ~ A.0 + L.1 + V, A.2 ~ A.1 + L.2 + V),
data = wd, method = "glm", stabilize = TRUE)
wd$cumA <- wd$A.0 + wd$A.1 + wd$A.2
fw <- glm(Y ~ cumA + V, data = wd, weights = wm$weights, family = gaussian)
ew <- ci_cl(fw, "cumA", ~id)
canon("msm_weightit.b.r", ew[["b"]]); canon("msm_weightit.lo.r", ew[["lo"]]); canon("msm_weightit.hi.r", ew[["hi"]])
}
cat("Simulated data (ข้อมูลจำลอง): no number above is evidence about any real drug.\n")
# ---- The true values built into the simulated cohort (exact: computed from the design, no sampling) ----
# The cohort was simulated from these probabilities. U is an unmeasured patient trait; L_6 is the day-6
# value of L (after the last treatment day, never recorded).
# P(U = 1) = pU and P(V = 1) = pV
# P(L_t = 1) = lBase + lU*U + lA*A_{t-1}, with A_{t-1} = 0 on day 0
# logit P(A_t = 1) = aBase + aL*L_t + aA*A_{t-1} + aV*V
# P(death by day 28) = yBase + yU*U + yV*V + yL*(L_0 + L_1 + L_2 + L_6) + yA*(A_0 + A_1 + A_2)
pU <- 0.40; pV <- 0.45
lBase <- 0.32; lU <- 0.45; lA <- -0.30
aBase <- -1.2; aL <- 2.0; aA <- 1.5; aV <- -0.5
yBase <- 0.05; yU <- 0.25; yV <- 0.05; yL <- 0.15; yA <- -0.005
design <- c(p_u = pU, p_v = pV, l_base = lBase, l_u = lU, l_a = lA, a_base = aBase, a_l = aL, a_a = aA,
a_v = aV, y_base = yBase, y_u = yU, y_v = yV, y_l = yL, y_a = yA)
for (k in names(design)) canon(paste0("truth.design.", k), design[[k]])
# Each treated day lowers the risk of death directly (yA) and by lowering the next day's L (yL x lA).
g1 <- yA + yL * lA
canon("truth.msm.g1_per_treated_day", g1)
canon("truth.msm.g1_direct_part", yA)
canon("truth.msm.g1_through_l_part", yL * lA)
# Risk of death by day 28 if every patient were treated on k of the three days (k = 0 never, k = 3 always).
# EL is P(L_t = 1) after an untreated day; L_0, L_1, L_2 and L_6 each add yL x EL, and each treated day adds g1.
EL <- lBase + lU * pU
r_never <- yBase + yU * pU + yV * pV + yL * 4 * EL
canon("truth.msm.risk_never_treated", r_never)
canon("truth.msm.risk_cum1", r_never + g1)
canon("truth.msm.risk_cum2", r_never + 2 * g1)
canon("truth.msm.risk_always_treated", r_never + 3 * g1)
canon("truth.msm.rd_always_vs_never", 3 * g1)
# The model with age V (the one fitted above): intercept at V = 0 and the coefficient of V.
canon("truth.msm.g0_conditional", yBase + yU * pU + yL * 4 * EL)
canon("truth.msm.g2_v", yV)
# Where each naive regression lands in an infinitely large cohort: list all 512 combinations of
# U, V, L_0, A_0, L_1, A_1, L_2, A_2 and L_6 with their probabilities under the observed treatment
# process, then fit the same regressions to the expected risk, weighting each combination by its probability.
g <- expand.grid(U = 0:1, V = 0:1, L0 = 0:1, A0 = 0:1, L1 = 0:1, A1 = 0:1, L2 = 0:1, A2 = 0:1, L6 = 0:1)
probL <- function(U, Aprev) lBase + lU * U + lA * Aprev
probA <- function(L, Aprev, V) plogis(aBase + aL * L + aA * Aprev + aV * V)
pick <- function(x, p) ifelse(x == 1, p, 1 - p)
g$p <- with(g, pick(U, pU) * pick(V, pV) * pick(L0, probL(U, 0)) * pick(A0, probA(L0, 0, V)) *
pick(L1, probL(U, A0)) * pick(A1, probA(L1, A0, V)) *
pick(L2, probL(U, A1)) * pick(A2, probA(L2, A1, V)) * pick(L6, probL(U, A2)))
g$risk <- with(g, yBase + yU * U + yV * V + yL * (L0 + L1 + L2 + L6) + yA * (A0 + A1 + A2))
g$cumA <- g$A0 + g$A1 + g$A2
# Naive 3 and 4 (one row per patient).
lim_cum <- coef(lm(risk ~ cumA, data = g, weights = p))[["cumA"]]
lim_cum_adj <- coef(lm(risk ~ cumA + V + L0 + L1 + L2, data = g, weights = p))[["cumA"]]
canon("truth.naive_limit.cum_crude", lim_cum)
canon("truth.naive_limit.cum_adj", lim_cum_adj)
# Naive 1 and 2 (one row per patient-day): three rows per combination, one per treatment day.
gd <- rbind(transform(g, A = A0, L = L0), transform(g, A = A1, L = L1), transform(g, A = A2, L = L2))
lim_unadj <- coef(lm(risk ~ A, data = gd, weights = p))[["A"]]
lim_adj <- coef(lm(risk ~ A + L + V, data = gd, weights = p))[["A"]]
canon("truth.naive_limit.unadj", lim_unadj)
canon("truth.naive_limit.adj", lim_adj)
# ---- summary tables (simulated data) ----
# Effect of one treated day on the risk of death by day 28: each estimate with its 95% CI (clustered by
# patient) beside the value it converges to in an infinitely large cohort (for the MSM: the true effect).
est3 <- function(e) unname(e[c("b", "lo", "hi")])
t_eff <- rbind(
naive_pooled = c(est3(ci_cl(f_unadj, "A", ~id)), lim_unadj),
naive_pooled_adj_L_V = c(est3(ci_cl(f_adj, "A", ~id)), lim_adj),
naive_cumulative = c(est3(ci_cl(f_cc, "cumA", ~id)), lim_cum),
naive_cumulative_adj = c(est3(ci_cl(f_ca, "cumA", ~id)), lim_cum_adj),
msm_stabilised_weights = c(est3(e_msm), g1),
msm_marginal = c(est3(e_m), g1))
colnames(t_eff) <- c("per_day", "lo", "hi", "large_sample_value")
cat("Effect of one treated day on death by day 28, simulated data\n")
print(round(t_eff, 4))
# Cumulative weights at day 4: stabilised versus unstabilised.
t_wt <- rbind(stabilised = c(mean(last$sw), max(last$sw), ess(last$sw)),
unstabilised = c(mean(last$w), max(last$w), ess(last$w)))
colnames(t_wt) <- c("mean", "max", "ess")
cat("Cumulative weights at day 4 (3,000 patients), simulated data\n")
print(round(t_wt, 4))
# Risk of death by day 28 under each regime, from the marginal MSM, beside the true risk.
t_risk <- rbind(never_treated = c(coef(msmm)[[1]], r_never),
always_treated = c(coef(msmm)[[1]] + 3 * coef(msmm)[[2]], r_never + 3 * g1))
colnames(t_risk) <- c("msm_estimate", "true_value")
cat("Risk of death by day 28 by treatment regime, simulated data\n")
print(round(t_risk, 4))
cat("These true values describe the simulation design (simulated data), not any real drug.\n")
# ---- Follow-up length, the marginal MSM intercept, and how many days each patient was treated ----
# Death is counted up to day 28.
canon("truth.design.death_day", 28, 0)
# The intercept of the marginal MSM is the risk of death by day 28 if no patient were ever treated.
canon("truth.msm.g0_marginal", r_never)
# Share who died, and patients treated on one or two of the three days.
canon("death_risk", mean(wide$Y))
canon("treated_one_day", sum(wide$cumA == 1), 0)
canon("treated_two_days", sum(wide$cumA == 2), 0)
print(table(treated_days = wide$cumA))
cat("Simulated data (ข้อมูลจำลอง): not evidence about any real drug or patient.\n") > cat("Effect of one treated day on death by day 28, simulated data\n")
Effect of one treated day on death by day 28, simulated data
> print(round(t_eff, 4))
per_day lo hi large_sample_value
naive_pooled 0.0859 0.0631 0.1087 0.0860
naive_pooled_adj_L_V -0.0128 -0.0363 0.0107 -0.0119
naive_cumulative 0.0676 0.0498 0.0854 0.0641
naive_cumulative_adj 0.0108 -0.0057 0.0273 0.0005
msm_stabilised_weights -0.0573 -0.0806 -0.0340 -0.0500
msm_marginal -0.0566 -0.0803 -0.0329 -0.0500
> t_wt <- rbind(stabilised = c(mean(last$sw), max(last$sw),
+ ess(last$sw)), unstabilised = c(mean(last$w), max(last$w),
+ ess(last$w)))
> colnames(t_wt) <- c("mean", "max", "ess")
> cat("Cumulative weights at day 4 (3,000 patients), simulated data\n")
Cumulative weights at day 4 (3,000 patients), simulated data
> print(round(t_wt, 4))
mean max ess
stabilised 0.9994 10.3018 1694.320
unstabilised 7.9648 94.9382 1456.026
> t_risk <- rbind(never_treated = c(coef(msmm)[[1]],
+ r_never), always_treated = c(coef(msmm)[[1]] + 3 * coef(msmm)[[2]],
+ r_never + 3 * g1))
> colnames(t_risk) <- c("msm_estimate", "true_value")
> cat("Risk of death by day 28 by treatment regime, simulated data\n")
Risk of death by day 28 by treatment regime, simulated data
> print(round(t_risk, 4))
msm_estimate true_value
never_treated 0.4888 0.4725
always_treated 0.3189 0.3225 Move the arrows yourself
Each bias depends on the strength of particular arrows. If organ dysfunction barely influenced treatment, leaving it out would cost little. If the drug did not change organ dysfunction, adjusting for it would neither block an effect nor open the collider path. The widget lets you vary the three arrows that matter and compare each estimate with the truth.
When patients leave early: censoring weights in the MSM
In the simulated cohort nobody is lost before death by day 28 is known, so no censoring weight is needed, but real cohorts lose patients, for example by transfer to another hospital. When leaving depends on organ dysfunction, those who remain no longer represent those who started. Inverse probability of censoring weighting (IPCW) handles this the same way as treatment: model the probability of remaining under follow-up on each day given the measured history, weight each remaining patient by one over the cumulative probability, and multiply that weight by the treatment weight [4]. Weighting for treatment, censoring and selection is compared side by side in The IPW Family: One Idea, Three Kinds of Missing People.
Assumptions, and the checks that go with them
An MSM estimates the regime contrast only under four conditions [7]. None can be proven from the data, but each has checks that can expose a problem.
Sequential exchangeability
Sequential exchangeability means that on every decision day, among patients with the same measured history, those who received the drug and those who did not have the same distribution of potential outcomes. In plain terms, nothing unmeasured pushes treatment and death together once $\bar L_t$, $\bar A_{t-1}$ and $V$ are known. In the simulated cohort, frailty is unmeasured but acts on treatment only through organ dysfunction, so the condition holds. If the team had also judged frailty at the bedside and treated by it, no weighting on measured variables would remove the bias.
Positivity on every day
Positivity requires a chance above zero of each treatment option on every decision day, for every measured history that occurs. If a unit protocol always gave drug X whenever severe organ dysfunction was present, no patient with organ dysfunction on a decision day could represent the never-treated regime. Strictly, the risk under a regime needs positivity only along the histories that the regime passes through. Check the predicted probabilities by day: probabilities near 0 or 1 show up as very large weights.
Consistency
Consistency means a patient's observed outcome equals the potential outcome under the regime the patient actually followed. It needs well-defined regimes: "drug X on day 2" should mean one dose range and one route, not whatever was charted.
Correct models
The weights are only as good as the treatment models, and the MSM needs the right form. The stabilised weights should average close to 1 on each day (0.999 here at day 4), and their spread by day should show no extreme tail; a mean near 1 is necessary but does not show the model is right.
Balance of $L_t$ between treated and untreated patients in the weighted data, within each level of the previous dose and age group, day by day, is a further check. Without that split, correctly built stabilised weights can still leave the treated with less organ dysfunction, because they were more often treated at the previous decision and that dose lowered it; with unstabilised weights, balance is expected even without the split. A refit with truncated weights or a more flexible MSM shows whether the estimate depends on these choices [4].
MSMs and target trial emulation: a method inside a design
Target trial emulation is a design framework. Before analysing observational data, the team writes down the randomised trial it would have run: eligibility, treatment strategies, assignment, time zero (the moment follow-up starts), follow-up, outcome and causal contrast. The analysis is then built to mimic that trial [8]. An MSM is an analysis method, so the two are not alternatives.
They combine naturally. An emulated trial of "always give drug X" versus "never give it" starts every patient at the same time zero. A patient who deviates from a strategy can be artificially censored at the deviation, which removes the patient from that strategy's comparison from then on.
Deviation usually depends on the same time-varying factors, such as organ dysfunction, so this censoring is informative: leaving depends on factors that also affect death. Weights for remaining uncensored correct for it, alongside treatment weights wherever treatment is still confounded. Target trial emulation fixes the question and the start of follow-up, and the MSM and its weights handle time-varying confounding. Not every emulated trial needs an MSM, for example when the strategies are fixed once at time zero.
The parametric g-formula in brief
The parametric g-formula reaches the same target from the other side [1, 7]. Instead of modelling treatment, it models the outcome and the time-varying confounders given the history. It then simulates what would happen to each patient under a regime, drawing $L_1$ from its model given the assigned $A_0$ and so on, and averages the predicted risks.
The truth calculation in the hand example was a g-formula done by hand. The method rests on the same causal conditions as the MSM, but different models must be right: the outcome and confounder models instead of the treatment models. Agreement between the two approaches is reassuring, because they can fail in different ways.
What an MSM does not license
Three limits deserve a sentence in any report that uses an MSM.
- It does nothing about unmeasured confounding of $A_t$, because weighting balances only what the treatment models contain.
- It estimates a population contrast of regimes: the difference in risk if everyone followed one regime rather than another. It does not say how any single patient would respond.
- It needs the regimes compared to be supported by the data. Here 471 patients were always treated and 588 never treated; a regime that almost nobody followed would rest on the model's form rather than on data.
The form adds a caveat of its own. A straight line in $\mathrm{cum}(\bar a)$ treats all regimes with the same number of treated days as equal, so it would misstate the regime risks if the timing of doses mattered.
Common misreadings and their fixes
-
"An MSM solves time-varying confounding."
Weighting removes confounding only by the variables in the treatment models, and only when those models are right and every history has a chance of each treatment. Frailty does no harm in the simulated cohort because it acts on treatment through organ dysfunction alone.
Fix: An MSM removes the bias from measured time-varying confounders affected by earlier treatment, provided the weight models are right and positivity holds. It does nothing about unmeasured confounding.
-
"Adjusting for organ dysfunction on each day is enough."
Adjustment removes the confounding of each dose by that day's organ dysfunction, but it blocks the effect of the earlier doses that runs through organ dysfunction and opens the collider path through frailty. In the simulated cohort the adjusted model that counts treated days lands at 0.0005 per treated day in a very large sample, against a true -0.050.
Fix: When earlier treatment changes a confounder, handle the confounder through the weights, not as a covariate in the outcome model.
-
"Organ dysfunction is both a mediator and a confounder."
The sentence is incomplete until it names the dose. Organ dysfunction on day 2 is a mediator for the day-0 dose and a confounder for the day-2 dose.
Fix: Name the treatment time: a time-varying confounder affected by earlier treatment lies on the path from the earlier dose and confounds the later one.
-
"Stabilising changes the answer."
In the hand example, unstabilised and stabilised weights give the same risks, 0.15 and 0.20, because each regime's risk is computed on its own, so there is no model form to get wrong. Stabilising narrows the spread of the weights, which usually narrows the confidence interval, without changing the target as long as the MSM's form is right and every baseline covariate in the numerator is also in the MSM. If the MSM's form is wrong, the two weightings can give different answers.
Fix: Check that the stabilised weights average close to 1, report their spread and maximum by day, and check the MSM's form against a more flexible version.
-
"Any covariate can go in the numerator."
A covariate measured after treatment starts cancels out of the weight and brings back the confounding the weight was meant to remove. A baseline covariate in the numerator but missing from the MSM is left confounding treatment in the weighted cohort.
Fix: Put only earlier treatment and baseline covariates in the numerator, and include every baseline numerator covariate in the MSM.
-
"Regressing death on each day's treatment, pooled over days, estimates the effect of the drug."
That regression repeats each patient's single outcome on every row and compares treated with untreated days. It targets no regime contrast, and in the simulated cohort it gives +0.086 per treated day for a drug that saves lives.
Fix: State the regimes first, then fit a model for the risk under each regime, such as $E[Y^{\bar a}] = g_0 + g_1\,\mathrm{cum}(\bar a)$.
-
"Target trial emulation and MSMs are the same thing."
Target trial emulation is a design; an MSM is one analysis method that can sit inside it.
Fix: Specify the target trial first, then choose the analysis that its strategies and confounding structure require.
What to do in your own analysis
- Draw the diagram over time, and mark every covariate that an earlier dose can change.
- Write the regimes and the estimand before fitting anything, for example always versus never treated over the first three decision days.
- Build the denominator from the history that plausibly drives treatment, including recent values of each time-varying confounder, and keep $L_t$ and other post-baseline variables out of the numerator. The simulated cohort records organ dysfunction and age as yes-or-no variables, but with real data continuous confounders, such as age, are usually best kept continuous in the weight models, for example with splines (flexible curves that let the effect of age bend).
- Report the mean, spread, maximum and effective sample size of the stabilised weights by day, and consider a refit with truncated weights as a sensitivity analysis.
- Use a cluster-robust standard error or a bootstrap that refits the weight models, and add censoring weights if patients leave before the outcome is known.
- Where feasible, compare the MSM with a g-formula analysis; a disagreement points to a model worth revisiting.
Glossary
- time-varying confounder (ตัวกวนที่เปลี่ยนตามเวลา)
- A covariate measured during follow-up that affects later treatment and the outcome, and may itself be changed by earlier treatment.
- treatment-confounder feedback
- The pattern in which earlier treatment changes a time-varying confounder that then drives later treatment.
- regime
- A rule that sets the treatment at every decision point, such as drug X on every decision day.
- sequential exchangeability
- On every decision day, treated and untreated patients with the same measured history have the same distribution of potential outcomes.
- marginal structural model (แบบจำลองโครงสร้างระดับประชากร (marginal structural model, MSM))
- A model for the mean potential outcome under each treatment regime, averaged over the time-varying confounders rather than conditioned on them, usually fitted with inverse probability weights.
- stabilised weight (น้ำหนักแบบปรับเสถียร (stabilised weight))
- An inverse probability weight multiplied by the probability of the same treatment given only earlier treatment and baseline covariates; for a correctly specified MSM it keeps the target when every baseline numerator covariate is also in the MSM, and it usually narrows the spread.
- IPCW (การถ่วงน้ำหนักด้วยส่วนกลับของความน่าจะเป็นที่ไม่ถูกเซ็นเซอร์ (IPCW))
- Inverse probability of censoring weighting: weighting each patient still under follow-up by one over the probability of having remained under follow-up, given the measured history.
- target trial emulation (การจำลองการทดลองเป้าหมาย (target trial emulation))
- A design framework that specifies the randomised trial an observational analysis is meant to mimic.
- g-formula (g-formula)
- Standardisation that averages model-based predicted outcomes over the covariate history each treatment regime would produce.
- collider
- A variable that two arrows on a path point into; adjusting for it opens the path and can create a spurious association.
- positivity (ข้อสมมติ positivity)
- A chance above zero of each treatment option for every measured history that occurs.
- consistency
- The observed outcome equals the potential outcome under the regime the patient actually followed.
References
- Daniel RM, Cousens SN, De Stavola BL, Kenward MG, Sterne JAC. Methods for dealing with time-dependent confounding. Stat Med. 2013;32(9):1584-1618. doi:10.1002/sim.5686 https://doi.org/10.1002/sim.5686
- Robins JM, Hernán MA, Brumback B. Marginal structural models and causal inference in epidemiology. Epidemiology. 2000;11(5):550-560. doi:10.1097/00001648-200009000-00011 https://doi.org/10.1097/00001648-200009000-00011
- Hernán MA, Brumback B, Robins JM. Marginal structural models to estimate the causal effect of zidovudine on the survival of HIV-positive men. Epidemiology. 2000;11(5):561-570. doi:10.1097/00001648-200009000-00012 https://doi.org/10.1097/00001648-200009000-00012
- Cole SR, Hernán MA. Constructing inverse probability weights for marginal structural models. Am J Epidemiol. 2008;168(6):656-664. doi:10.1093/aje/kwn164 https://doi.org/10.1093/aje/kwn164
- Fewell Z, Hernán MA, Wolfe F, Tilling K, Choi H, Sterne JAC. Controlling for time-dependent confounding using marginal structural models. Stata J. 2004;4(4):402-420. doi:10.1177/1536867X0400400403 https://doi.org/10.1177/1536867X0400400403
- van der Wal WM, Geskus RB. ipw: an R package for inverse probability weighting. J Stat Softw. 2011;43(13):1-23. doi:10.18637/jss.v043.i13 https://doi.org/10.18637/jss.v043.i13
- Hernán MA, Robins JM. Causal inference: what if [Internet]. Boca Raton: Chapman & Hall/CRC; 2020. https://miguelhernan.org/whatifbook
- Hernán MA, Robins JM. Using big data to emulate a target trial when a randomized trial is not available. Am J Epidemiol. 2016;183(8):758-764. doi:10.1093/aje/kwv254 https://doi.org/10.1093/aje/kwv254
Key takeaways
- A covariate that earlier treatment changes, and that drives later treatment and death, makes ordinary regression fail whether or not it is adjusted for.
- A marginal structural model compares whole treatment regimes by weighting each patient by the inverse probability of the treatment history received, multiplied over decision days.
- Stabilised weights keep the target when the MSM's form is right and it holds every baseline numerator covariate; only earlier treatment and baseline covariates belong in the numerator.
- In the simulated ICU cohort, ordinary regressions made the drug look harmful or nearly useless, while the MSM gave -0.057 per treated day against a true -0.050.
- An MSM removes bias only from measured confounders, so its credibility rests on sequential exchangeability, positivity on every day and correct weight models.
Related in the wiki: [[inverse-probability-weighting-treatment-censoring-selection]] [[extreme-propensity-weights-stabilization-truncation-trimming]]