Chapter 12: IP Weighting and Marginal Structural Models
This chapter uses inverse probability (IP) weighting to estimate the average causal effect of smoking cessation on weight gain from the observational NHEFS data. Chapter 2 introduced IP weighting as a nonparametric method; here IP weights are estimated with models, which lets us handle many covariates and nondichotomous treatments, and the weighted analysis is interpreted as fitting a marginal structural model.
This chapter is based on Hernán and Robins (2020, chap. 12, pp. 163-174).
R and Stata code for this chapter’s programs (Programs 12.x) is in Tom Palmer’s cibookex-r companion (GPL-3.0), which descends from the authors’ own code.
Key concepts: IP weighting creates a pseudo-population in which the measured confounders \(L\) no longer predict treatment \(A\). Under conditional exchangeability, positivity, and consistency, association is causation in that pseudo-population, so a simple weighted regression of \(Y\) on \(A\) estimates the parameters of a marginal structural model.
NHEFS stands for the National Health and Nutrition Examination Survey Data I Epidemiologic Follow-up Study. The book’s website provides the subset of the NHEFS data used in Part II (Hernán and Robins 2020, chap. 12, p. 163).
1 12.1 The Causal Question (pp. 163-164)
Part II is organized around one question: on average, how much does quitting smoking change body weight?
1.1 Research Question
1.2 Unadjusted Comparison
The associational difference \(\operatorname{E}\mathopen{}\left[Y \mid A = 1\right]\mathclose{} - \operatorname{E}\mathopen{}\left[Y \mid A = 0\right]\mathclose{}\) need not equal the causal difference if quitters and non-quitters differ in characteristics that affect weight gain.
Age as an example: quitters were on average 4 years older than non-quitters, and older people gained less weight than younger people whether or not they quit. Age is therefore a (surrogate) confounder, and the unadjusted estimate of 2.5 kg might underestimate the causal effect (Hernán and Robins 2020, chap. 12, pp. 163-164).
1.3 Baseline Characteristics (Table 12.1)
| Mean baseline characteristic | Quitters (\(A=1\)) | Non-quitters (\(A=0\)) |
|---|---|---|
| Age, years | 46.2 | 42.8 |
| Men, % | 54.6 | 46.6 |
| White, % | 91.1 | 85.4 |
| University, % | 15.4 | 9.9 |
| Weight, kg | 72.4 | 70.3 |
| Cigarettes/day | 18.6 | 21.2 |
| Years smoking | 26.0 | 24.1 |
| Little exercise, % | 40.7 | 37.9 |
| Inactive life, % | 11.2 | 8.9 |
1.4 Measured Confounders
Chapter 18 discusses how to select confounders; here the 9-variable adjustment set is taken as given (Hernán and Robins 2020, chap. 12, p. 164).
The smoking cessation example is convenient (no deep subject-matter knowledge required, public data) but carries a potential selection bias (Hernán and Robins 2020, Fine Point 12.1, p. 164):
- Treated individuals are those who were smokers at baseline and reported having quit in the 1982 survey.
- Answering the 1982 survey requires surviving and remaining under follow-up until 1982, so inclusion in the study is conditional on an event that occurs after treatment (smoking cessation) has started.
- If treatment affects selection into the study, the analysis can suffer from the selection bias described in Chapter 8.
- Smoking cessation is really a time-varying treatment (people quit at different times); Part II ignores this, and Part III handles time-varying treatments.
A randomized trial would not have this problem, because treatment would be assigned at baseline and known even for people who did not reach the 1982 visit. Unlike the censoring of the outcome handled in Section 12.6, here the missing information concerns the treatment itself; it can be addressed through sensitivity analysis. The book uses the example deliberately to illustrate a common problem in observational analyses that emulate a target trial: misalignment of eligibility, treatment assignment, and the start of follow-up.
2 12.2 Estimating IP Weights via Modeling (pp. 164-167)
2.1 IP Weights Definition
The conditional probability of treatment \(\Pr[A = 1 \mid L]\) is the propensity score; Chapter 15 discusses propensity scores in detail (Hernán and Robins 2020, chap. 12, p. 165). Note that the propensity score is \(\Pr[A = 1 \mid L]\), not \(f(A \mid L)\): for an untreated individual, the denominator of \(W^A\) is one minus the propensity score.
IP weighting reweights individuals so that \(L\) no longer predicts \(A\). Chapter 2 did this by hand (IP weights in Chapter 2); the next definition makes the reweighted population precise for a general weight.
If the \(n\) individuals are drawn independently from the population, so that each \(W_i\) has the distribution of \(W\), the expected size of the pseudo-population is \(\operatorname{E}\mathopen{}\left[\sum_{i=1}^n W_i\right]\mathclose{} = \operatorname{E}\mathopen{}\left[W\right]\mathclose{} \cdot n\).
Different weights create different pseudo-populations; Proposition 1 is about the IP weights \(W = W^A\) of Definition 1.
Proof. As in the statement, \(k\) is the number of treatment levels.
Step 0 (property 0: the weights cancel the treatment probability). For every level \(a\) and every \(l\) with \(\Pr[L = l] > 0\), \[ \begin{aligned} \operatorname{E}\mathopen{}\left[W^A \mathbb{1}\mathopen{}\left(A = a\right)\mathclose{} \mathbb{1}\mathopen{}\left(L = l\right)\mathclose{}\right]\mathclose{} &= \operatorname{E}\mathopen{}\left[\frac{\mathbb{1}\mathopen{}\left(A = a\right)\mathclose{} \mathbb{1}\mathopen{}\left(L = l\right)\mathclose{}}{f(a \mid l)}\right]\mathclose{} && \text{(} W^A = 1/f(a \mid l) \text{ whenever } A = a \text{ and } L = l\text{)}\\ &= \frac{\operatorname{E}\mathopen{}\left[\mathbb{1}\mathopen{}\left(A = a\right)\mathclose{} \mathbb{1}\mathopen{}\left(L = l\right)\mathclose{}\right]\mathclose{}}{f(a \mid l)} && \text{(the constant } 1/f(a \mid l) > 0 \text{ factors out, by positivity)}\\ &= \frac{\Pr[A = a, L = l]}{f(a \mid l)} && \text{(expectation of an indicator)}\\ &= \frac{f(a \mid l) \Pr[L = l]}{f(a \mid l)} && \text{(definition of } f(a \mid l)\text{)}\\ &= \Pr[L = l] && \text{(cancel } f(a \mid l)\text{)}. \end{aligned} \] Summing over all \(k\) levels \(a\) and all \(l\) gives \(\operatorname{E}\mathopen{}\left[W^A\right]\mathclose{} = k \sum_l \Pr[L = l] = k\), which is property 0.
Property 1. \[ \begin{aligned} \Pr_{ps}[A = a, L = l] &= \frac{\operatorname{E}\mathopen{}\left[W^A \mathbb{1}\mathopen{}\left(A = a\right)\mathclose{} \mathbb{1}\mathopen{}\left(L = l\right)\mathclose{}\right]\mathclose{}}{\operatorname{E}\mathopen{}\left[W^A\right]\mathclose{}} && \text{(definition of the pseudo-population)}\\ &= \frac{\Pr[L = l]}{k} && \text{(Step 0, twice)}. \end{aligned} \] Summing over \(l\) gives \(\Pr_{ps}[A = a] = 1/k\); summing over \(a\) gives \(\Pr_{ps}[L = l] = \Pr[L = l]\). So \(\Pr_{ps}[A = a, L = l] = \Pr_{ps}[A = a] \Pr_{ps}[L = l]\) for all \(a\) and \(l\), which is independence.
Property 2. \[ \begin{aligned} \operatorname{E}_{ps}\mathopen{}\left[Y \mid A = a\right]\mathclose{} &= \frac{\operatorname{E}_{ps}\mathopen{}\left[Y \mathbb{1}\mathopen{}\left(A = a\right)\mathclose{}\right]\mathclose{}}{\Pr_{ps}[A = a]} && \text{(definition of the pseudo-population)}\\ &= \frac{\operatorname{E}\mathopen{}\left[W^A Y \mathbb{1}\mathopen{}\left(A = a\right)\mathclose{}\right]\mathclose{} / \operatorname{E}\mathopen{}\left[W^A\right]\mathclose{}}{1/k} && \text{(definition of the pseudo-population; Property 1)}\\ &= \frac{\operatorname{E}\mathopen{}\left[W^A Y \mathbb{1}\mathopen{}\left(A = a\right)\mathclose{}\right]\mathclose{} / k}{1/k} && \text{(} \operatorname{E}\mathopen{}\left[W^A\right]\mathclose{} = k \text{, Step 0)}\\ &= \operatorname{E}\mathopen{}\left[W^A Y \mathbb{1}\mathopen{}\left(A = a\right)\mathclose{}\right]\mathclose{} && \text{(cancel } k\text{)}\\ &= \sum_l \operatorname{E}\mathopen{}\left[W^A Y \mathbb{1}\mathopen{}\left(A = a\right)\mathclose{} \mid L = l\right]\mathclose{} \Pr[L = l] && \text{(law of total expectation over } L\text{)}\\ &= \sum_l \frac{\operatorname{E}\mathopen{}\left[Y \mathbb{1}\mathopen{}\left(A = a\right)\mathclose{} \mid L = l\right]\mathclose{}}{f(a \mid l)} \Pr[L = l] && \text{(} W^A = 1/f(a \mid l) \text{ when } A = a \text{ and } L = l\text{; pull the constant out)}\\ &= \sum_l \frac{\operatorname{E}\mathopen{}\left[Y \mid A = a, L = l\right]\mathclose{} f(a \mid l)}{f(a \mid l)} \Pr[L = l] && \text{(} \operatorname{E}\mathopen{}\left[Y \mathbb{1}\mathopen{}\left(A = a\right)\mathclose{} \mid L = l\right]\mathclose{} = \operatorname{E}\mathopen{}\left[Y \mid A = a, L = l\right]\mathclose{} \Pr[A = a \mid L = l]\text{)}\\ &= \sum_l \operatorname{E}\mathopen{}\left[Y \mid A = a, L = l\right]\mathclose{} \Pr[L = l] && \text{(cancel } f(a \mid l)\text{)}. \end{aligned} \]
Exchangeability in the pseudo-population. Fix a level \(a\), a set \(B\) of outcome values, and a level \(a'\). \[ \begin{aligned} \Pr_{ps}[Y^a \in B, A = a'] &= \frac{\operatorname{E}\mathopen{}\left[W^A \mathbb{1}\mathopen{}\left(Y^a \in B\right)\mathclose{} \mathbb{1}\mathopen{}\left(A = a'\right)\mathclose{}\right]\mathclose{}}{\operatorname{E}\mathopen{}\left[W^A\right]\mathclose{}} && \text{(definition of the pseudo-population)}\\ &= \frac{1}{k} \operatorname{E}\mathopen{}\left[W^A \mathbb{1}\mathopen{}\left(Y^a \in B\right)\mathclose{} \mathbb{1}\mathopen{}\left(A = a'\right)\mathclose{}\right]\mathclose{} && \text{(} \operatorname{E}\mathopen{}\left[W^A\right]\mathclose{} = k \text{, Step 0)}\\ &= \frac{1}{k} \sum_l \operatorname{E}\mathopen{}\left[W^A \mathbb{1}\mathopen{}\left(Y^a \in B\right)\mathclose{} \mathbb{1}\mathopen{}\left(A = a'\right)\mathclose{} \mid L = l\right]\mathclose{} \Pr[L = l] && \text{(law of total expectation over } L\text{)}\\ &= \frac{1}{k} \sum_l \frac{\operatorname{E}\mathopen{}\left[\mathbb{1}\mathopen{}\left(Y^a \in B\right)\mathclose{} \mathbb{1}\mathopen{}\left(A = a'\right)\mathclose{} \mid L = l\right]\mathclose{}}{f(a' \mid l)} \Pr[L = l] && \text{(} W^A = 1/f(a' \mid l) \text{ when } A = a' \text{ and } L = l\text{; pull the constant out)}\\ &= \frac{1}{k} \sum_l \frac{\Pr[Y^a \in B, A = a' \mid L = l]}{f(a' \mid l)} \Pr[L = l] && \text{(expectation of an indicator)}\\ &= \frac{1}{k} \sum_l \frac{\Pr[Y^a \in B \mid L = l] f(a' \mid l)}{f(a' \mid l)} \Pr[L = l] && \text{(conditional exchangeability } Y^a \perp\!\!\!\perp A \mid L\text{)}\\ &= \frac{1}{k} \sum_l \Pr[Y^a \in B \mid L = l] \Pr[L = l] && \text{(cancel } f(a' \mid l)\text{)}\\ &= \frac{1}{k} \Pr[Y^a \in B] && \text{(law of total probability)}. \end{aligned} \] Summing over the \(k\) levels \(a'\) gives \(\Pr_{ps}[Y^a \in B] = \Pr[Y^a \in B]\). Hence \(\Pr_{ps}[Y^a \in B, A = a'] = \Pr_{ps}[Y^a \in B] \Pr_{ps}[A = a']\) (using \(\Pr_{ps}[A = a'] = 1/k\) from Property 1), which is \(Y^a \perp\!\!\!\perp A\) in the pseudo-population.
Causal mean. Under conditional exchangeability, positivity, and consistency, the standardization theorem of Chapter 7 (standardization under conditional exchangeability) gives \(\operatorname{E}\mathopen{}\left[Y^a\right]\mathclose{} = \sum_l \operatorname{E}\mathopen{}\left[Y \mid L = l, A = a\right]\mathclose{} \Pr[L = l]\), which equals \(\operatorname{E}_{ps}\mathopen{}\left[Y \mid A = a\right]\mathclose{}\) by Property 2. This holds for every \(a\), so every contrast of the pseudo-population means equals the matching causal contrast.
By Proposition 1, \(\operatorname{E}\mathopen{}\left[W^A\right]\mathclose{} = k\), so for a dichotomous treatment the expected size of this pseudo-population is \(2n\), twice the study population.
2.2 Why Modeling Is Needed
2.3 The Treatment Model
What the model assumes: on the logit scale, the relation between each continuous covariate and the probability of quitting is a parabola, and each covariate contributes to the logit independently of the others. Under these parametric restrictions, we get an estimate \(\mathop{\widehat{\Pr}}\nolimits\mathopen{}\left[A = 1 \mid L\right]\mathclose{}\) for every combination of \(L\) values, and hence for each of the 1566 individuals (Hernán and Robins 2020, chap. 12, p. 165).
2.4 Estimating the Effect in the Pseudo-Population
Fit the (saturated) linear mean model \[\operatorname{E}\mathopen{}\left[Y \mid A\right]\mathclose{} = \theta_0 + \theta_1 A\] by weighted least squares with weights \(\hat{W}^A\).
The model \(\operatorname{E}\mathopen{}\left[Y \mid A\right]\mathclose{} = \theta_0 + \theta_1 A\) is saturated: it has two parameters for the two quantities \(\operatorname{E}\mathopen{}\left[Y \mid A = 1\right]\mathclose{}\) and \(\operatorname{E}\mathopen{}\left[Y \mid A = 0\right]\mathclose{}\), and \(\theta_1 = \operatorname{E}\mathopen{}\left[Y \mid A = 1\right]\mathclose{} - \operatorname{E}\mathopen{}\left[Y \mid A = 0\right]\mathclose{}\). With \(\hat{W}_i = 1\) for everyone, weighted least squares reduces to the ordinary least squares of Chapter 11 (Hernán and Robins 2020, chap. 12, pp. 165-166).
If there is no confounding in the pseudo-population (Definition 2) and the model for \(\Pr[A = 1 \mid L]\) is correct, an unbiased estimator of \(\operatorname{E}_{ps}\mathopen{}\left[Y \mid A = 1\right]\mathclose{} - \operatorname{E}_{ps}\mathopen{}\left[Y \mid A = 0\right]\mathclose{}\) is also an unbiased estimator of \(\operatorname{E}\mathopen{}\left[Y^{a=1}\right]\mathclose{} - \operatorname{E}\mathopen{}\left[Y^{a=0}\right]\mathclose{}\).
The mean of the nonstabilized weights is about 2 because \(\operatorname{E}\mathopen{}\left[W^A\right]\mathclose{} = k = 2\) for a dichotomous treatment (property 0 of Proposition 1).
2.5 Confidence Intervals for IP Weighted Estimates
If the model for \(\Pr[A = 1 \mid L]\) is misspecified, the estimates of \(\theta_0\) and \(\theta_1\) will be biased and, as in Chapter 11, the intervals may fall short of their nominal 95% coverage (Hernán and Robins 2020, chap. 12, p. 167).
Two sample estimators of the IP weighted mean \(\operatorname{E}\mathopen{}\left[Y^a\right]\mathclose{}\) differ only in how they normalize the weights.
Why the two agree in the limit under positivity (Hernán and Robins 2020, Technical Point 12.1, p. 166). The Hajek estimator is (asymptotically) unbiased for the ratio \[\frac{\operatorname{E}\mathopen{}\left[\mathbb{1}\mathopen{}\left(A = a\right)\mathclose{} Y / f(A \mid L)\right]\mathclose{}}{\operatorname{E}\mathopen{}\left[\mathbb{1}\mathopen{}\left(A = a\right)\mathclose{} / f(A \mid L)\right]\mathclose{}}.\] Under positivity, the denominator equals 1: \[ \begin{aligned} \operatorname{E}\mathopen{}\left[\frac{\mathbb{1}\mathopen{}\left(A = a\right)\mathclose{}}{f(A \mid L)}\right]\mathclose{} &= \operatorname{E}\mathopen{}\left[\operatorname{E}\mathopen{}\left[\frac{\mathbb{1}\mathopen{}\left(A = a\right)\mathclose{}}{f(a \mid L)} \,\middle|\, L\right]\mathclose{}\right]\mathclose{} && \text{(iterated expectations; } \mathbb{1}\mathopen{}\left(A=a\right)\mathclose{} \text{ forces } A = a\text{)}\\ &= \operatorname{E}\mathopen{}\left[\frac{\Pr[A = a \mid L]}{f(a \mid L)}\right]\mathclose{} && \text{(} \operatorname{E}\mathopen{}\left[\mathbb{1}\mathopen{}\left(A = a\right)\mathclose{} \mid L\right]\mathclose{} = \Pr[A = a \mid L]\text{)}\\ &= \operatorname{E}\mathopen{}\left[1\right]\mathclose{} = 1 && \text{(} f(a \mid L) = \Pr[A = a \mid L] > 0 \text{ by positivity)} \end{aligned} \] so the ratio equals the IP weighted mean \(\operatorname{E}\mathopen{}\left[\mathbb{1}\mathopen{}\left(A = a\right)\mathclose{} Y / f(A \mid L)\right]\mathclose{}\), which equals \(\operatorname{E}\mathopen{}\left[Y^a\right]\mathclose{}\) under conditional exchangeability, positivity, and consistency (Technical Point 3.1).
Why the Hajek estimator is preferred in practice: even when \(f(A \mid L)\) is replaced by predictions from a misspecified model, the Hajek estimator stays within \([0, 1]\) for a dichotomous \(Y\); the Horvitz-Thompson estimator need not.
Without positivity: let \(Q(a) = \{l : \Pr(A = a \mid L = l) > 0\}\). The Hajek ratio then equals \(\sum_l \operatorname{E}\mathopen{}\left[Y \mid A = a, L = l, L \in Q(a)\right]\mathclose{} \Pr[L = l \mid L \in Q(a)]\), which under conditional exchangeability and consistency is \(\operatorname{E}\mathopen{}\left[Y^a \mid L \in Q(a)\right]\mathclose{}\). Its denominator converges to \(\Pr[Q(a)]\) rather than 1, so the ratio of the Horvitz-Thompson limit to the Hajek limit is \(\Pr[Q(a)]\), and the difference of Hajek estimators for \(a = 1\) versus \(a = 0\) has no causal interpretation.
3 12.3 Stabilized IP Weights (pp. 167-169)
The goal of IP weighting is a pseudo-population (Definition 2) in which \(L\) does not predict \(A\). The weights \(1/f(A \mid L)\) are only one way to get there.
3.1 Why Rescaling Does Not Change the Estimate
Proof. \[ \begin{aligned} \frac{\mathop{\hat{\operatorname{E}}}\nolimits\mathopen{}\left[p \, \mathbb{1}\mathopen{}\left(A = a\right)\mathclose{} Y / f(A \mid L)\right]\mathclose{}}{\mathop{\hat{\operatorname{E}}}\nolimits\mathopen{}\left[p \, \mathbb{1}\mathopen{}\left(A = a\right)\mathclose{} / f(A \mid L)\right]\mathclose{}} &= \frac{p \, \mathop{\hat{\operatorname{E}}}\nolimits\mathopen{}\left[\mathbb{1}\mathopen{}\left(A = a\right)\mathclose{} Y / f(A \mid L)\right]\mathclose{}}{p \, \mathop{\hat{\operatorname{E}}}\nolimits\mathopen{}\left[\mathbb{1}\mathopen{}\left(A = a\right)\mathclose{} / f(A \mid L)\right]\mathclose{}} && \text{(constants factor out of averages)}\\ &= \frac{\mathop{\hat{\operatorname{E}}}\nolimits\mathopen{}\left[\mathbb{1}\mathopen{}\left(A = a\right)\mathclose{} Y / f(A \mid L)\right]\mathclose{}}{\mathop{\hat{\operatorname{E}}}\nolimits\mathopen{}\left[\mathbb{1}\mathopen{}\left(A = a\right)\mathclose{} / f(A \mid L)\right]\mathclose{}} && \text{(cancel } p\text{)} \end{aligned} \] The same holds for every treatment level \(a\), so the difference between levels is unchanged.
3.2 Stabilized Weights
The key requirement is only that treatment probability in the pseudo-population not depend on \(L\); it may differ between people. A common choice gives the treated probability \(\Pr[A = 1]\) and the untreated \(\Pr[A = 0]\), as in the original population.
The stabilized weights \(SW^A\) create their own pseudo-population (Definition 2), different from the one created by \(W^A\).
Applied to the data of Figure 2.1, where \(\Pr[A = 1] = 13/20 = 0.65\) and \(\Pr[A = 0] = 7/20 = 0.35\), the weights \(f(A)/f(A \mid L)\) produce a pseudo-population (Figure 12.1 of the book) resembling a randomized experiment in which 65% of individuals are assigned to \(A = 1\) and 35% to \(A = 0\). Preserving the 65/35 ratio requires non-integer numbers of individuals in some branches, which is no problem mathematically (Hernán and Robins 2020, chap. 12, p. 167).
With \(A\) and \(L\) each taking finitely many values and positivity holding, as in Proposition 1, and with every sum over \(l\) running over the values with \(\Pr[L = l] > 0\), the stabilized weights have mean 1: \[ \begin{aligned} \operatorname{E}\mathopen{}\left[SW^A\right]\mathclose{} &= \sum_a \sum_l \operatorname{E}\mathopen{}\left[SW^A \mathbb{1}\mathopen{}\left(A = a\right)\mathclose{} \mathbb{1}\mathopen{}\left(L = l\right)\mathclose{}\right]\mathclose{} && \text{(linearity, since } \textstyle\sum_a \sum_l \mathbb{1}\mathopen{}\left(A = a\right)\mathclose{} \mathbb{1}\mathopen{}\left(L = l\right)\mathclose{} = 1 \text{ almost surely)}\\ &= \sum_a \sum_l f(a) \operatorname{E}\mathopen{}\left[W^A \mathbb{1}\mathopen{}\left(A = a\right)\mathclose{} \mathbb{1}\mathopen{}\left(L = l\right)\mathclose{}\right]\mathclose{} && \text{(} SW^A = f(a) W^A \text{ whenever } A = a\text{)}\\ &= \sum_a \sum_l f(a) \Pr[L = l] && \text{(Step 0 of the proof of the pseudo-population proposition)}\\ &= \sum_a f(a) \sum_l \Pr[L = l] && \text{(} f(a) \text{ does not depend on } l\text{)}\\ &= \sum_a f(a) && \text{(} \textstyle\sum_l \Pr[L = l] = 1\text{)}\\ &= 1 && \text{(} \textstyle\sum_a f(a) = 1\text{)}. \end{aligned} \] So, for individuals sampled as described after Definition 2, the pseudo-population created by \(SW^A\) has expected size \(\operatorname{E}\mathopen{}\left[SW^A\right]\mathclose{} \cdot n = n\), the same as the study population, and the sample mean \(\frac{1}{n} \sum_{i=1}^n SW^A_i\) has expectation 1. When the models for \(f(A \mid L)\) and \(f(A)\) are correct and \(n\) is large, the mean of the estimated stabilized weights should therefore be close to 1, and in data analyses one should check it: a mean far from 1 suggests model misspecification or (near) violations of positivity (Hernán and Robins 2020, chap. 12, p. 168).
To estimate the average causal effect in the treated rather than in the whole population, use \(\Pr[A = 1 \mid L]\) as the numerator instead (Technical Point 4.1).
3.3 Estimation
3.4 Why Stabilize?
In the NHEFS, there are 4 white women aged 66, and none of them quit smoking: positivity is empirically violated in that stratum.
- Structural violations: some values of \(L\) make treatment (or non-treatment) impossible. Example: people off work cannot be exposed to an occupational chemical. Zero cells appear whenever we condition on such a confounder.
- Random violations: in a finite sample, stratifying on many confounders produces zero cells by chance, even though the probability of treatment is not zero in the target population. Example: no treated white women aged 66 or 67, but some treated white women aged 65 and 69.
The two kinds of violation have different consequences (Hernán and Robins 2020, Fine Point 12.2, p. 169):
- Under structural violations, neither IP weighting nor standardization supports inference about the entire population; inference has to be limited to the strata where structural positivity holds (see Technical Point 12.1).
- Under random violations, the parametric model smooths over the zeros: the logistic model of Section 12.2 estimated the probability of quitting for white women aged 66 by interpolating from everyone else.
Whenever we use parametric IP weights in the presence of zero cells, as in the estimate \(\hat\theta_1 = 3.4\), we are effectively assuming that the nonpositivity is random.
4 12.4 Marginal Structural Models (pp. 169-171)
Proof. \[ \begin{aligned} \operatorname{E}\mathopen{}\left[Y^{a=1}\right]\mathclose{} - \operatorname{E}\mathopen{}\left[Y^{a=0}\right]\mathclose{} &= (\beta_0 + \beta_1 \cdot 1) - (\beta_0 + \beta_1 \cdot 0) && \text{(evaluate the model at } a = 1, 0\text{)}\\ &= \beta_1 && \text{(cancel } \beta_0\text{)} \end{aligned} \]
So Sections 12.2-12.3 were already estimating \(\beta_1\) of a marginal structural model.
4.1 Fitting a Marginal Structural Model by IP Weighting
# weighted least squares; use robust (sandwich) standard errors
fit <- lm(weight_change ~ quit_smoking,
weights = stabilized_weights,
data = nhefs)The model \(\operatorname{E}\mathopen{}\left[Y^a\right]\mathclose{} = \beta_0 + \beta_1 a\) is saturated because \(A\) is dichotomous: two unknowns on each side (\(\operatorname{E}\mathopen{}\left[Y^{a=1}\right]\mathclose{}\) and \(\operatorname{E}\mathopen{}\left[Y^{a=0}\right]\mathclose{}\) on the left, \(\beta_0\) and \(\beta_1\) on the right). That is why sample averages in the pseudo-population sufficed (Hernán and Robins 2020, chap. 12, p. 170).
Null preservation (Chapter 9): if treatment has no average causal effect, any model for the marginal mean \(\operatorname{E}\mathopen{}\left[Y^a\right]\mathclose{}\) that has an intercept, with treatment terms whose coefficients can all be 0, is correctly specified. For example, \(\operatorname{E}\mathopen{}\left[Y^a\right]\mathclose{} = \beta_0 + \beta_1 a + \beta_2 a^2\) is then correctly specified with \(\beta_1 = \beta_2 = 0\) and \(\beta_0 = \operatorname{E}\mathopen{}\left[Y^a\right]\mathclose{}\) for every \(a\).
4.2 Continuous Treatment: Change in Smoking Intensity
New treatment: \(A\) = cigarettes per day in 1982 minus cigarettes per day at baseline (e.g., \(-25\) or \(+40\)), restricted to the 1162 individuals whose baseline smoking was at most 25 cigarettes per day.
Goal: \(\operatorname{E}\mathopen{}\left[Y^a\right]\mathclose{} - \operatorname{E}\mathopen{}\left[Y^{a'}\right]\mathclose{}\) for any \(a, a'\). A saturated model is impractical with dozens of treatment values, so we posit a nonsaturated marginal structural model (a parabolic dose-response curve):
4.3 Effect of Increasing Intensity by 20 Cigarettes/Day
To estimate \(\beta_1, \beta_2\): estimate \(SW^A\), then fit \(\operatorname{E}\mathopen{}\left[Y \mid A\right]\mathclose{} = \theta_0 + \theta_1 A + \theta_2 A^2\) to the pseudo-population.
4.4 Weights for a Continuous Treatment
For continuous \(A\), \(f(A \mid L)\) is a probability density function, which is generally hard to estimate when \(L\) is high-dimensional.
IP weighting for continuous treatments is often dangerous: the estimates can change a great deal with the choice of model or algorithm for the conditional density \(f(A \mid L)\). The constant-variance assumption seemed reasonable after inspecting a residuals plot, and other choices (e.g., a truncated normal with heteroscedasticity) gave similar estimates. Finding more stable ways to estimate these weights is ongoing research (Hernán and Robins 2020, chap. 12, pp. 170-171).
The book reports these quantities from the unrounded coefficients. Plugging in the rounded coefficients gives a slightly different value: \[ \begin{aligned} \hat \beta_0 + 20\hat \beta_1 + 400\hat \beta_2 &= 2.005 + 20(-0.109) + 400(0.003)\\ &= 2.005 - 2.180 + 1.200\\ &= 1.025, \end{aligned} \] because \(\hat \beta_2 = 0.003\) is rounded and is multiplied by 400. Use the reported 0.9 kg, not this back-of-the-envelope value.
4.5 Dichotomous Outcome: Marginal Structural Logistic Model
For a dichotomous treatment, \(\operatorname{exp}\mathopen{}\left\{\alpha_1\right\}\mathclose{}\) in Definition 8 is the causal odds ratio comparing \(a = 1\) with \(a = 0\). \(\operatorname{exp}\mathopen{}\left\{\alpha_1\right\}\mathclose{}\) is estimated by \(\operatorname{exp}\mathopen{}\left\{\hat\theta_1\right\}\mathclose{}\) from fitting \(\operatorname{logit}\Pr[D = 1 \mid A] = \theta_0 + \theta_1 A\) to the IP weighted pseudo-population (Definition 2).
This model is saturated for a dichotomous treatment; for a continuous treatment one would specify a nonsaturated logistic model (Hernán and Robins 2020, chap. 12, p. 171).
5 12.5 Effect Modification and Marginal Structural Models (pp. 171-172)
Marginal structural models for the population average causal effect include no covariates. To assess effect modification by a covariate \(V\) (which need not be a confounder), add it to the model.
Strictly, this is no longer a marginal model because it conditions on \(V\), but the name “marginal structural model” is still used (Hernán and Robins 2020, chap. 12, p. 171).
The parameter \(\beta_3\) does not generally have a causal interpretation as the effect of \(V\): exchangeability, positivity, and consistency are assumed for treatment \(A\), not for \(V\).
5.1 Fitting the Model
Either \(SW^A\) or \(SW^A(V)\) (Definition 10) can be used to fit the model.
\(SW^A(V)\) generally gives narrower confidence intervals than \(SW^A\); some intuition is that the variance of \(SW^A(V)\) is smaller than that of \(SW^A\) (Hernán and Robins 2020, chap. 12, pp. 171-172).
5.2 Example: Effect Modification by Sex
5.3 Which Covariates to Put in the Model?
For a faux marginal structural model (Definition 11), \(SW^A(L) = f[A \mid L] / f[A \mid L] = 1\), so no weighting is done: the model is fit as the unweighted outcome regression that fully adjusts for \(L\) (Chapter 15).
Teaching point: confounding and effect modification are logically distinct, but students often conflate them because stratification (Chapter 4) and regression (Chapter 15) are used for both. With marginal structural models the two tasks use distinct tools: IP weighting for confounder adjustment, and treatment-covariate product terms in the model for effect modification (Hernán and Robins 2020, chap. 12, p. 172).
Interaction versus effect modification: to study the interaction between two treatments \(A\) and \(B\) (Chapter 5), include parameters for both in the marginal structural model, use the joint probability of both treatments in the denominator of the IP weights, and assume exchangeability, positivity, and consistency for both \(A\) and \(B\).
6 12.6 Censoring and Missing Data (pp. 172-174)
The analysis so far used the 1566 individuals with a 1982 weight measurement. Another 63 eligible individuals were excluded because their 1982 weight was unknown.
The analyses of Sections 12.2-12.4 actually fit \[\operatorname{E}\mathopen{}\left[Y \mid A, C = 0\right]\mathclose{} = \theta_0 + \theta_1 A\] among individuals with \(C = 0\).
6.1 Can Censoring Cause Selection Bias Here?
6.2 The Target: A Joint Effect
The superscript \(c = 0\) makes explicit what many investigators mean by “the causal effect of \(A\)” even when they omit it.
6.3 IP Weights for Treatment and Censoring
If some components of \(L\) are affected by treatment \(A\) (as in Figure 8.4), the conditional independence \(Y^{a,c=0} \perp\!\!\!\perp(A, C) \mid L\) will not generally hold; Part III presents alternative exchangeability conditions under which IP weighting still works (Hernán and Robins 2020, chap. 12, p. 173).
Some variables in \(L\) may have zero coefficients in the model for \(f(A \mid L)\) but not in the model for \(\Pr[C = 0 \mid L, A]\), or vice versa. Even so, in large samples, efficiency is always gained by including in both models every variable in \(L\) that independently predicts the outcome.
6.4 Nonstabilized versus Stabilized Censoring Weights
| Nonstabilized \(W^C = 1/\Pr[C = 0 \mid L, A]\) | Stabilized \(SW^C = \Pr[C = 0 \mid A] / \Pr[C = 0 \mid L, A]\) | |
|---|---|---|
| Pseudo-population size | original population before censoring (about \(1566 + 63 = 1629\)) | original population after censoring (about 1566) |
| Arrows removed | from both \(L\) and \(A\) into \(C\) | from \(L\) into \(C\) |
| Result | no selection, hence no selection bias | selection at random with respect to \(L\): selection but no selection bias |
With either, fit the weighted model \(\operatorname{E}\mathopen{}\left[Y \mid A, C = 0\right]\mathclose{} = \theta_0 + \theta_1 A\) to estimate the marginal structural model \(\operatorname{E}\mathopen{}\left[Y^{a, c=0}\right]\mathclose{} = \beta_0 + \beta_1 a\); the stabilized version uses \(SW^{A,C}\) (Definition 16).
The estimated \(SW^C\) have mean 1 when the model for \(\Pr[C = 0 \mid A]\) is correctly specified (Hernán and Robins 2020, chap. 12, p. 173).
Chapter 13 describes the alternative to IP weighting for adjusting for confounding and selection bias: standardization (Hernán and Robins 2020, chap. 12, p. 174).
The stabilized weights \(f[A]/f[A \mid L]\) belong to a larger class \(g[A]/f[A \mid L]\), where \(g[A]\) is any positive function of \(A\) that is not a function of \(L\). With unsaturated marginal structural models, such weights are preferable to \(1/f[A \mid L]\) because some \(g[A]\) (often \(f[A]\)) yield more efficient estimators (Hernán and Robins 2020, Technical Point 12.2, p. 174).
The IP weighted mean \(\operatorname{E}\mathopen{}\left[g(A) \mathbb{1}\mathopen{}\left(A = a\right)\mathclose{} Y / f(A \mid L)\right]\mathclose{}\) no longer equals \(\operatorname{E}\mathopen{}\left[Y^a\right]\mathclose{}\), but the Hajek version does. Under conditional exchangeability, positivity, and consistency: \[ \begin{aligned} \operatorname{E}\mathopen{}\left[\frac{g(A) \mathbb{1}\mathopen{}\left(A = a\right)\mathclose{} Y}{f(A \mid L)}\right]\mathclose{} &= g(a) \, \operatorname{E}\mathopen{}\left[\frac{\mathbb{1}\mathopen{}\left(A = a\right)\mathclose{} Y}{f(A \mid L)}\right]\mathclose{} && \text{(} \mathbb{1}\mathopen{}\left(A = a\right)\mathclose{} \text{ forces } g(A) = g(a)\text{)}\\ &= g(a) \, \operatorname{E}\mathopen{}\left[Y^a\right]\mathclose{} && \text{(Technical Point 3.1)}\\ \operatorname{E}\mathopen{}\left[\frac{g(A) \mathbb{1}\mathopen{}\left(A = a\right)\mathclose{}}{f(A \mid L)}\right]\mathclose{} &= g(a) \, \operatorname{E}\mathopen{}\left[\frac{\mathbb{1}\mathopen{}\left(A = a\right)\mathclose{}}{f(A \mid L)}\right]\mathclose{} && \text{(same reason)}\\ &= g(a) && \text{(Technical Point 12.1)} \end{aligned} \] so the ratio of the two equals \(g(a)\operatorname{E}\mathopen{}\left[Y^a\right]\mathclose{}/g(a) = \operatorname{E}\mathopen{}\left[Y^a\right]\mathclose{}\) (since \(g(a) > 0\)).
The Hajek mean is the solution \(u\) of \(\operatorname{E}\mathopen{}\left[\frac{g[A] \mathbb{1}\mathopen{}\left(A = a\right)\mathclose{}}{f[A \mid L]} (Y - u)\right]\mathclose{} = 0\), since solving for \(u\) gives the ratio of the two expectations. (The book prints this equation without the indicator \(\mathbb{1}\mathopen{}\left(A = a\right)\mathclose{}\) (Hernán and Robins 2020, Technical Point 12.2, p. 174); without it, the solution would average \(Y\) over all treatment levels rather than estimate \(\operatorname{E}\mathopen{}\left[Y^a\right]\mathclose{}\).) Similarly, for \(\operatorname{E}\mathopen{}\left[Y^a\right]\mathclose{} = \beta_0 + \beta_1 a\), the weighted least squares estimators with weights \(g[A]/f[A \mid L]\) solve \[\mathop{\hat{\operatorname{E}}}\nolimits\mathopen{}\left[\frac{g[A]}{f[A \mid L]} \mathopen{}\left[Y - (\beta_0 + \beta_1 A)\right]\mathclose{} \begin{pmatrix} 1 \\ A \end{pmatrix}\right]\mathclose{} = 0,\] and \(\hat \beta_0\) and \(\hat \beta_0 + \hat \beta_1\) are exactly the Hajek estimators of \(\operatorname{E}\mathopen{}\left[Y^{a=0}\right]\mathclose{}\) and \(\operatorname{E}\mathopen{}\left[Y^{a=1}\right]\mathclose{}\). For a treatment \(A\) with finitely many levels, under conditional exchangeability, positivity, and consistency, in the pseudo-population created by \(W = g[A]/f[A \mid L]\) (Definition 2), the mean of \(Y\) given \(A = a\) still equals \(\operatorname{E}\mathopen{}\left[Y^a\right]\mathclose{}\).
7 Summary
Key concepts introduced:
- IP weighting via modeling: estimate \(\Pr[A = 1 \mid L]\) with a parametric model when \(L\) is high-dimensional
- Stabilized weights \(f(A)/f(A \mid L)\): same estimate as nonstabilized weights for a saturated model, generally narrower confidence intervals for nonsaturated models
- Marginal structural models: models for the mean counterfactual outcome, fit by IP weighted regression
- Effect modification: add \(V\) (and \(V \times a\) terms) to the marginal structural model; use \(SW^A(V)\)
- Censoring weights: \(W^{A,C} = W^A \times W^C\) target the joint effect \(\operatorname{E}\mathopen{}\left[Y^{a, c=0}\right]\mathclose{}\)
NHEFS results: unadjusted 2.5 kg; IP weighted 3.4 kg (95% CI: 2.4, 4.5); also adjusting for censoring 3.5 kg (95% CI: 2.5, 4.5).
- Requires a correctly specified model for treatment (and censoring)
- Positivity violations; parametric models assume random nonpositivity
- Weights for continuous treatments can be very sensitive to the density model
- Robust-variance confidence intervals are conservative
Looking ahead: Chapter 13 estimates the same effect by standardization (the parametric g-formula), and Part III extends IP weighting and marginal structural models to time-varying treatments.