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?

Research Question

Example 1 (NHEFS: Does Quitting Smoking Cause Weight Gain?) Population: 1566 cigarette smokers aged 25-74 years who had an NHEFS baseline visit (1971-75) and a follow-up visit about 10 years later (1982).

Treatment: \(A = 1\) if the individual reported quitting smoking before the follow-up visit, \(A = 0\) otherwise.

Outcome: \(Y\) = body weight at follow-up minus body weight at baseline (kg).

Causal estimand (average causal effect on the additive scale): \[\operatorname{E}\mathopen{}\left[Y^{a=1}\right]\mathclose{} - \operatorname{E}\mathopen{}\left[Y^{a=0}\right]\mathclose{}\]

Unadjusted Comparison

Example 2 (NHEFS: Unadjusted Comparison) Most people gained weight, but quitters gained more on average:

• \(\mathop{\hat{\operatorname{E}}}\nolimits\mathopen{}\left[Y \mid A = 1\right]\mathclose{} = 4.5\) kg among quitters
• \(\mathop{\hat{\operatorname{E}}}\nolimits\mathopen{}\left[Y \mid A = 0\right]\mathclose{} = 2.0\) kg among non-quitters
• Associational difference: \(4.5 - 2.0 = 2.5\) kg (95% CI: 1.7, 3.4)

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).

Baseline Characteristics (Table 12.1)

Table 1: Mean baseline characteristics of quitters and non-quitters in the NHEFS. Source: Hernán and Robins (2020, Table 12.1, p. 163).

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

Measured Confounders

Example 3 (NHEFS: Measured Confounders) We assume that these 9 baseline variables, collected in the vector \(L\), suffice to adjust for confounding:

• Sex (0: male, 1: female)
• Age (years)
• Race (0: white, 1: other)
• Education (5 categories)
• Intensity of smoking (cigarettes per day)
• Duration of smoking (years)
• Physical activity in daily life (3 categories)
• Recreational exercise (3 categories)
• Weight (kg)

Assumption: conditional exchangeability given \(L\): \[Y^a \perp\!\!\!\perp A \mid L\]

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).

Fine Point 12.1: Setting a Bad Example

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)

IP Weights Definition

Definition 1 (Inverse Probability Weights) The individual-specific IP weights for treatment \(A\) are

\[W^A \stackrel{\text{def}}{=}\frac{1}{f(A \mid L)}\]

i.e., the inverse of the conditional probability of receiving the treatment level the individual actually received.

For a dichotomous treatment:

• If \(A = 1\): \(W^A = \dfrac{1}{\Pr[A = 1 \mid L]}\)
• If \(A = 0\): \(W^A = \dfrac{1}{\Pr[A = 0 \mid L]} = \dfrac{1}{1 - \Pr[A = 1 \mid L]}\)

So only \(\Pr[A = 1 \mid L]\) needs to be estimated.

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.

Definition 2 (Pseudo-Population) Let each of the \(n\) individuals of a study population carry a weight \(W \geq 0\) that is a function of the observed variables, with \(0 < \operatorname{E}\mathopen{}\left[W\right]\mathclose{} < \infty\). The pseudo-population created by \(W\) is the population in which each individual is replicated \(W\) times (a non-integer number of copies is allowed). Its size is \(\sum_{i=1}^n W_i\). Its distribution is defined by its expectations: for any random variable \(h\) that is a function of the observed variables and the counterfactual outcomes, with \(\operatorname{E}\mathopen{}\left[W \mathopen{}\left|h\right|\mathclose{}\right]\mathclose{} < \infty\), \[\operatorname{E}_{ps}\mathopen{}\left[h\right]\mathclose{} \stackrel{\text{def}}{=}\frac{\operatorname{E}\mathopen{}\left[W h\right]\mathclose{}}{\operatorname{E}\mathopen{}\left[W\right]\mathclose{}}.\] Probabilities are expectations of indicators, \(\Pr_{ps}[E] \stackrel{\text{def}}{=}\operatorname{E}_{ps}\mathopen{}\left[\mathbb{1}\mathopen{}\left(E\right)\mathclose{}\right]\mathclose{}\) for any event \(E\), and conditional means are \(\operatorname{E}_{ps}\mathopen{}\left[h \mid A = a\right]\mathclose{} \stackrel{\text{def}}{=}\operatorname{E}_{ps}\mathopen{}\left[h \mathbb{1}\mathopen{}\left(A = a\right)\mathclose{}\right]\mathclose{} / \Pr_{ps}[A = a]\) whenever \(\Pr_{ps}[A = a] > 0\).

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.

Proposition 1 (Properties of the IP Weighted Pseudo-Population) Suppose that

• \(A\) and \(L\) take finitely many values, and the pseudo-population (Definition 2) is created by \(W = W^A = 1/f(A \mid L)\) (Definition 1);
• positivity holds: \(f(a \mid l) = \Pr[A = a \mid L = l] > 0\) for every treatment level \(a\) and every \(l\) with \(\Pr[L = l] > 0\);
• the outcome is integrable: \(\operatorname{E}\mathopen{}\left[\mathopen{}\left|Y\right|\mathclose{}\right]\mathclose{} < \infty\), and \(\operatorname{E}\mathopen{}\left[\mathopen{}\left|Y^a\right|\mathclose{}\right]\mathclose{} < \infty\) for every \(a\).

Then, with \(k\) the number of treatment levels:

0. \(\operatorname{E}\mathopen{}\left[W^A\right]\mathclose{} = k\);
1. in the pseudo-population, \(A\) and \(L\) are statistically independent; and
2. the mean \(\operatorname{E}_{ps}\mathopen{}\left[Y \mid A = a\right]\mathclose{}\) in the pseudo-population equals the standardized mean \(\sum_l \operatorname{E}\mathopen{}\left[Y \mid A = a, L = l\right]\mathclose{} \Pr[L = l]\) in the actual population.

Here and in the proof, every sum over \(l\) runs over the values \(l\) with \(\Pr[L = l] > 0\). These three properties need no exchangeability. If, in addition, conditional exchangeability \(Y^a \perp\!\!\!\perp A \mid L\) holds in the actual population for every \(a\), and consistency holds (\(Y = Y^a\) for every individual with \(A = a\)), then:

• unconditional exchangeability \(Y^a \perp\!\!\!\perp A\) holds in the pseudo-population,
• \(\operatorname{E}\mathopen{}\left[Y^a\right]\mathclose{} = \operatorname{E}_{ps}\mathopen{}\left[Y \mid A = a\right]\mathclose{}\) for every \(a\), and
• association is causation in the pseudo-population: any contrast of the means \(\operatorname{E}_{ps}\mathopen{}\left[Y \mid A = a\right]\mathclose{}\) across levels \(a\) equals the same contrast of the counterfactual means \(\operatorname{E}\mathopen{}\left[Y^a\right]\mathclose{}\).

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.

Why Modeling Is Needed

Remark 1 (Nonparametric Weights Are Infeasible Here). In Section 2.4, \(\Pr[A = 1 \mid L]\) was estimated nonparametrically, by the proportion treated within each stratum of \(L\). That is impossible here:

• Even recoding all 9 confounders except age to at most 6 categories each gives a tree with over 2 million branches.
• With the actual ranges of smoking intensity, duration, and weight, there are many millions of strata.
• 1566 individuals cannot yield meaningful stratum-specific estimates across millions of strata.

This is the curse of dimensionality introduced in Chapter 10, so we resort to a model.

The Treatment Model

Algorithm 1 (Estimating IP Weights with a Treatment Model)  

1. Fit a logistic regression model for \(\Pr[A = 1 \mid L]\) including all 9 confounders:

   • linear and quadratic terms for the (quasi-)continuous covariates age, weight, intensity of smoking, and duration of smoking;
   • no product terms between covariates.

2. For each individual, compute

   • \(\hat{f}(1 \mid L) = \mathop{\widehat{\Pr}}\nolimits\mathopen{}\left[A = 1 \mid L\right]\mathclose{}\) if \(A = 1\),
   • \(\hat{f}(0 \mid L) = 1 - \mathop{\widehat{\Pr}}\nolimits\mathopen{}\left[A = 1 \mid L\right]\mathclose{}\) if \(A = 0\).

3. Set \(\hat{W}^A = 1 / \hat{f}(A \mid L)\).

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).

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\).

Definition 3 (Weighted Least Squares) Given weights \(\hat{W}_i\), the weighted least squares estimates \(\hat\theta_0, \hat\theta_1\) of the model \(\operatorname{E}\mathopen{}\left[Y \mid A\right]\mathclose{} = \theta_0 + \theta_1 A\) are the values that minimize \[\sum_i \hat{W}_i \mathopen{}\left[Y_i - (\theta_0 + \theta_1 A_i)\right]\mathclose{}^2.\]

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{}\).

Example 4 (NHEFS: Nonstabilized IP Weighting)  

• Estimated weights \(\hat{W}^A\) ranged from 1.05 to 16.7, with mean 2.00.
• \(\hat\theta_1 = 3.4\) kg: under our assumptions, quitting smoking increases weight by 3.4 kg on average.
• Conservative 95% confidence interval: (2.4, 4.5).

Compare with the unadjusted difference of 2.5 kg.

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).

Confidence Intervals for IP Weighted Estimates

Remark 2 (Confidence Intervals Must Account for the Weighting). Three options (Hernán and Robins 2020, chap. 12, p. 166):

1. Derive the variance analytically from statistical theory; the analyst must program the estimator, which standard software generally does not provide.
2. Nonparametric bootstrap; requires computing resources (or patience) for large data sets.
3. Robust (“sandwich”) variance estimator, as used for GEE models with an independent working correlation; a standard option in most software.

The robust-variance intervals are valid but conservative: their actual coverage of the super-population parameter exceeds the nominal 95%. All confidence intervals for IP weighted estimates in this chapter are of this conservative kind.

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).

Technical Point 12.1: Horvitz-Thompson and Hajek Estimators

Two sample estimators of the IP weighted mean \(\operatorname{E}\mathopen{}\left[Y^a\right]\mathclose{}\) differ only in how they normalize the weights.

Definition 4 (Horvitz-Thompson and Hajek Estimators) With \(\mathop{\hat{\operatorname{E}}}\nolimits\mathopen{}\left[\cdot\right]\mathclose{}\) the sample average:

• Horvitz-Thompson estimator (f known): \[\mathop{\hat{\operatorname{E}}}\nolimits\mathopen{}\left[\frac{\mathbb{1}\mathopen{}\left(A = a\right)\mathclose{} Y}{f(A \mid L)}\right]\mathclose{}\]
• Hajek estimator (modified Horvitz-Thompson): \[\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{}}\]

For binary \(A\), the IP weighted least squares estimate \(\hat\theta_0 + \hat\theta_1 a\) is the Hajek estimator.

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.

Remark 3 (Other Weights Give the Same Pseudo-Population Contrast).

• With weights \(1/f(A \mid L)\), a dichotomous treatment yields a pseudo-population of expected size \(\operatorname{E}\mathopen{}\left[W\right]\mathclose{} \cdot n = 2n\), twice the size of the study population: one treated copy and one untreated copy.
• With weights \(0.5/f(A \mid L)\), everyone has probability 0.5 of each treatment level regardless of \(L\); this pseudo-population has expected size \(\operatorname{E}\mathopen{}\left[W\right]\mathclose{} \cdot n = n\), the same as the study population.
• More generally, any weights \(p/f(A \mid L)\) with \(0 < p \leq 1\) give the same effect estimate.

Why Rescaling Does Not Change the Estimate

Proposition 2 (Rescaling the Weights Does Not Change the Estimate) For a constant \(p > 0\), the Hajek (weighted least squares) mean for treatment level \(a\) computed with weights \(p/f(A \mid L)\) equals the one computed with weights \(1/f(A \mid L)\), so the difference between treatment levels is unchanged too.

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.

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.

Definition 5 (Stabilized IP Weights) \[SW^A \stackrel{\text{def}}{=}\frac{f(A)}{f(A \mid L)}\]

For dichotomous \(A\):

• If \(A = 1\): \(SW^A = \dfrac{\Pr[A = 1]}{\Pr[A = 1 \mid L]}\)
• If \(A = 0\): \(SW^A = \dfrac{\Pr[A = 0]}{\Pr[A = 0 \mid L]} = \dfrac{1 - \Pr[A = 1]}{1 - \Pr[A = 1 \mid L]}\)

The weights \(W^A = 1/f(A \mid L)\) are called nonstabilized weights, and \(f(A)\) is the stabilizing factor.

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).

Estimation

Algorithm 2 (Estimating Stabilized IP Weights)  

• Denominator: same logistic model for \(\Pr[A = 1 \mid L]\) as in Section 12.2.
• Numerator: \(\mathop{\widehat{\Pr}}\nolimits\mathopen{}\left[A = 1\right]\mathclose{} = 403/1566\), equivalently a saturated (intercept-only) logistic model.
• Fit \(\operatorname{E}\mathopen{}\left[Y \mid A\right]\mathclose{} = \theta_0 + \theta_1 A\) by weighted least squares with weights \(\mathop{\widehat{\Pr}}\nolimits\mathopen{}\left[A = 1\right]\mathclose{}/\mathop{\widehat{\Pr}}\nolimits\mathopen{}\left[A = 1 \mid L\right]\mathclose{}\) for quitters and \((1 - \mathop{\widehat{\Pr}}\nolimits\mathopen{}\left[A = 1\right]\mathclose{})/(1 - \mathop{\widehat{\Pr}}\nolimits\mathopen{}\left[A = 1 \mid L\right]\mathclose{})\) for non-quitters.

Example 5 (NHEFS: Stabilized IP Weighting)  

• Estimated \(\hat{SW}^A\) ranged from 0.33 to 4.30, with mean 1.00 (versus 1.05 to 16.7, mean 2.00, for \(\hat{W}^A\)).
• \(\hat\theta_1 = 3.4\) kg (95% CI: 2.4, 4.5), the same estimate as with the nonstabilized weights.

Why Stabilize?

Remark 4 (When Stabilization Helps). If both kinds of weights give the same estimate, why stabilize?

• Stabilized weights typically give narrower confidence intervals.
• But this advantage can occur only when the weighted model is not saturated.
• Here, \(\operatorname{E}\mathopen{}\left[Y \mid A\right]\mathclose{} = \theta_0 + \theta_1 A\) is saturated because \(A\) is dichotomous, so the two sets of weights agree.
• With time-varying or continuous treatments the weighted model cannot be saturated, and stabilized weights are used (Section 12.4).

Fine Point 12.2: Checking Positivity

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)

Definition 6 (Marginal Structural Model) A marginal structural mean model is a model for the marginal mean of a counterfactual outcome, such as

\[\operatorname{E}\mathopen{}\left[Y^a\right]\mathclose{} = \beta_0 + \beta_1 a\]

Because the outcome \(Y^a\) is counterfactual (generally unobserved), the model cannot be fit directly to the data of any real-world study.

Proposition 3 (MSM Treatment Parameters Are Average Causal Effects) In the marginal structural model \(\operatorname{E}\mathopen{}\left[Y^a\right]\mathclose{} = \beta_0 + \beta_1 a\) of Definition 6, \[\beta_1 = \operatorname{E}\mathopen{}\left[Y^{a=1}\right]\mathclose{} - \operatorname{E}\mathopen{}\left[Y^{a=0}\right]\mathclose{}.\]

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.

Fitting a Marginal Structural Model by IP Weighting

Algorithm 3 (Estimating Marginal Structural Model Parameters by IP Weighting)  

1. Use IP weights to create a pseudo-population (Definition 2) in which, under the assumptions of Proposition 1, association is causation.
2. Fit the associational model \(\operatorname{E}\mathopen{}\left[Y \mid A\right]\mathclose{} = \theta_0 + \theta_1 A\) to the pseudo-population, i.e., by IP weighted least squares.
3. \(\theta_1\) then has the same causal interpretation as \(\beta_1\), so a consistent estimator \(\hat\theta_1\) is a consistent estimator of \(\beta_1 = \operatorname{E}\mathopen{}\left[Y^{a=1}\right]\mathclose{} - \operatorname{E}\mathopen{}\left[Y^{a=0}\right]\mathclose{}\).

Listing 1: Pseudo-code for fitting a marginal structural model by IP weighted least squares.

# 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\).

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):

Definition 7 (Nonsaturated Marginal Structural Model for a Continuous Treatment) \[\operatorname{E}\mathopen{}\left[Y^a\right]\mathclose{} = \beta_0 + \beta_1 a + \beta_2 a^2\]

where \(\beta_0 = \operatorname{E}\mathopen{}\left[Y^{a=0}\right]\mathclose{}\) is the mean weight gain with no change in smoking intensity.

Effect of Increasing Intensity by 20 Cigarettes/Day

Example 6 (Effect of 20 More Cigarettes per Day) Under the model of Definition 7, \[ \begin{aligned} \operatorname{E}\mathopen{}\left[Y^{a=20}\right]\mathclose{} - \operatorname{E}\mathopen{}\left[Y^{a=0}\right]\mathclose{} &= (\beta_0 + 20\beta_1 + 20^2 \beta_2) - \beta_0 && \text{(evaluate the model at } a = 20, 0\text{)}\\ &= 20\beta_1 + 400\beta_2 && \text{(cancel } \beta_0\text{; } 20^2 = 400\text{)} \end{aligned} \]

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.

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.

Example 7 (NHEFS: Weights for Change in Smoking Intensity)  

• \(f(A \mid L)\): normal with mean \(\mu_L = \operatorname{E}\mathopen{}\left[A \mid L\right]\mathclose{}\) (estimated by linear regression) and constant variance \(\sigma^2\) (the residual variance);
• \(f(A)\) in the numerator: also normal;
• estimated \(\hat{SW}^A\) ranged from 0.19 to 5.10, with mean 1.00.

IP Weights for Continuous Treatments Are Fragile

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).

Example 8 (NHEFS: Change in Smoking Intensity) IP weighted estimates: \(\hat \beta_0 = 2.005\), \(\hat \beta_1 = -0.109\), \(\hat \beta_2 = 0.003\).

Counterfactual mean	Estimate (95% CI)
\(\operatorname{E}\mathopen{}\left[Y^{a=0}\right]\mathclose{}\) (no change in intensity)	2.0 kg (1.4, 2.6)
\(\operatorname{E}\mathopen{}\left[Y^{a=20}\right]\mathclose{}\) (+20 cigarettes/day)	0.9 kg (\(-1.7\), 3.5)
\(\operatorname{E}\mathopen{}\left[Y^{a=20}\right]\mathclose{} - \operatorname{E}\mathopen{}\left[Y^{a=0}\right]\mathclose{}\)	\(0.9 - 2.0 = -1.1\) kg

Rounded Coefficients

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.

Dichotomous Outcome: Marginal Structural Logistic Model

Definition 8 (Marginal Structural Logistic Model) For a dichotomous outcome \(D\) (1: yes, 0: no), a marginal structural logistic model is a model for the counterfactual risk on the logit scale, such as \[\operatorname{logit}\Pr[D^a = 1] = \alpha_0 + \alpha_1 a.\]

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).

Example 9 (NHEFS: Quitting Smoking and Death by 1992) With \(A\) = quitting smoking and \(D\) = death by 1992, the IP weighted estimate of the causal odds ratio was \(\operatorname{exp}\mathopen{}\left\{\hat\theta_1\right\}\mathclose{} = 1.0\) (95% CI: 0.8, 1.4).

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.

Definition 9 (Marginal Structural Model with an Effect Modifier) \[\operatorname{E}\mathopen{}\left[Y^a \mid V\right]\mathclose{} = \beta_0 + \beta_1 a + \beta_2 V a + \beta_3 V\]

Additive effect modification by \(V\) is present if \(\beta_2 \neq 0\). The causal effect within level \(V = v\) is \[ \begin{aligned} \operatorname{E}\mathopen{}\left[Y^{a=1} \mid V = v\right]\mathclose{} - \operatorname{E}\mathopen{}\left[Y^{a=0} \mid V = v\right]\mathclose{} &= (\beta_0 + \beta_1 \cdot 1 + \beta_2 v \cdot 1 + \beta_3 v) - (\beta_0 + \beta_1 \cdot 0 + \beta_2 v \cdot 0 + \beta_3 v) && \text{(model evaluated at } a = 1 \text{ and } a = 0 \text{, with } V = v\text{)}\\ &= (\beta_0 + \beta_1 + \beta_2 v + \beta_3 v) - (\beta_0 + \beta_3 v) && \text{(multiply out the factors 1 and 0)}\\ &= \beta_1 + \beta_2 v && \text{(cancel } \beta_0 \text{ and } \beta_3 v\text{)}. \end{aligned} \]

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\).

Fitting the Model

Algorithm 4 (Fitting a Marginal Structural Model with an Effect Modifier)  

• Fit \(\operatorname{E}\mathopen{}\left[Y \mid A, V\right]\mathclose{} = \theta_0 + \theta_1 A + \theta_2 V A + \theta_3 V\) by weighted least squares with weights \(W^A\) or \(SW^A\).
• In most settings \(L\) should include \(V\); even when \(V\) is not needed for exchangeability, including it generally improves efficiency.
• Because the model is within levels of \(V\), the numerator can be \(f[A \mid V]\).

Definition 10 (Stabilized IP Weights Given \(V\)) \[SW^A(V) \stackrel{\text{def}}{=}\frac{f[A \mid V]}{f[A \mid L]}\]

The stabilized weights \(SW^A(V)\) are estimated like \(SW^A\), but with \(V\) added to the numerator’s logistic 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).

Example: Effect Modification by Sex

Example 10 (NHEFS: Does the Effect of Quitting Vary by Sex?) Let \(V\) = sex (0: male, 1: female). The 95% confidence interval for \(\hat\theta_2\) (the product-term coefficient) was \((-2.2, 1.9)\): no strong evidence of additive effect modification by sex.

Which Covariates to Put in the Model?

Remark 5 (Choosing the Covariates of a Marginal Structural Model).

• The subset \(V\) of \(L\) included in the marginal structural model should reflect only substantive interest: include \(V\) if you believe it may modify the effect and you care more about effects within levels of \(V\) than in the whole population.
• If all of \(L\) is included, IP weighting is unnecessary.

Definition 11 (Faux Marginal Structural Model) A faux marginal structural model is a marginal structural model that includes all of \(L\) as covariates. The book coins the name half-jokingly.

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).

Confounding and Effect Modification Use Different Tools

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.

Definition 12 (Censoring Indicator) The censoring indicator \(C\) equals 1 if an individual’s outcome is unmeasured (censored) and 0 if it is measured.

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\).

Can Censoring Cause Selection Bias Here?

Example 11 (NHEFS: Censoring as a Possible Source of Selection Bias) Restricting to uncensored individuals is expected to induce bias, even under the null, when \(C\) is a collider on a pathway between \(A\) and \(Y\), or a descendant of one (Figures 8.3-8.6). The NHEFS data are consistent with that structure:

• treatment predicts censoring: 5.8% of quitters versus 3.2% of non-quitters were censored;
• a predictor of \(Y\) predicts censoring: mean baseline weight was 76.6 kg in the censored versus 70.8 kg in the uncensored.

The Target: A Joint Effect

Definition 13 (Causal Effect in the Absence of Censoring) The average causal effect of \(A\) had nobody been censored is \[\operatorname{E}\mathopen{}\left[Y^{a=1, c=0}\right]\mathclose{} - \operatorname{E}\mathopen{}\left[Y^{a=0, c=0}\right]\mathclose{},\] a joint effect of \(A\) and \(C\) (Chapter 8).

The superscript \(c = 0\) makes explicit what many investigators mean by “the causal effect of \(A\)” even when they omit it.

IP Weights for Treatment and Censoring

Definition 14 (IP Weights for Censoring) \[ W^C \stackrel{\text{def}}{=} \begin{cases} \dfrac{1}{\Pr[C = 0 \mid L, A]} & \text{for uncensored individuals} \\ 0 & \text{for censored individuals} \end{cases} \]

\[W^{A,C} \stackrel{\text{def}}{=}W^A \times W^C = \frac{1}{f(A, C = 0 \mid L)}\]

using the factorization \(f(A, C = 0 \mid L) = f(A \mid L) \times \Pr[C = 0 \mid L, A]\).

Definition 15 (Identifiability Conditions for the Joint Treatment \((A, C)\)) For each treatment level \(a\) compared in the marginal structural model, the identifiability conditions for the joint treatment \((A, C)\) are:

• exchangeability: \(Y^{a, c=0} \perp\!\!\!\perp(A, C) \mid L\);
• positivity: \(\Pr[A = a, C = 0 \mid L = l] > 0\) for every \(l\) with \(\Pr[L = l] > 0\);
• consistency: \(Y = Y^{a, c=0}\) for every individual with \(A = a\) and \(C = 0\).

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.

Nonstabilized versus Stabilized Censoring Weights

Definition 16 (Stabilized IP Weights for Censoring) \[ SW^C \stackrel{\text{def}}{=} \begin{cases} \dfrac{\Pr[C = 0 \mid A]}{\Pr[C = 0 \mid L, A]} & \text{for uncensored individuals} \\ 0 & \text{for censored individuals} \end{cases} \]

\[SW^{A,C} \stackrel{\text{def}}{=}SW^A \times SW^C\]

	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).

Example 12 (NHEFS: Adjusting for Censoring)  

• Fit a logistic model for \(\Pr[C = 0 \mid L, A]\) to all 1629 individuals, with the same covariates as the treatment model; obtain \(\hat{SW}^C\) for the 1566 uncensored individuals.
• \(\hat{SW}^{A,C} = \hat{SW}^A \times \hat{SW}^C\) ranged from 0.35 to 4.09, with mean 1.00.
• \(\hat\theta_1 = 3.5\) kg (95% CI: 2.5, 4.5).

This is almost the same as the estimate with \(SW^A\) alone (3.4 kg), so either censoring introduces no selection bias here, or the measured covariates are not enough to remove it.

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).

Technical Point 12.2: More on Stabilized Weights

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:

1. IP weighting via modeling: estimate \(\Pr[A = 1 \mid L]\) with a parametric model when \(L\) is high-dimensional
2. Stabilized weights \(f(A)/f(A \mid L)\): same estimate as nonstabilized weights for a saturated model, generally narrower confidence intervals for nonsaturated models
3. Marginal structural models: models for the mean counterfactual outcome, fit by IP weighted regression
4. Effect modification: add \(V\) (and \(V \times a\) terms) to the marginal structural model; use \(SW^A(V)\)
5. 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).

Cautions

• 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.

Hernán, Miguel A, and James M Robins. 2020. Causal Inference: What If. Chapman & Hall/CRC. https://miguelhernan.org/whatifbook.