Chapter 20 showed that, with treatment-confounder feedback, traditional adjustment methods give a non-null estimate even when a time-varying treatment has no effect at all. This chapter shows that the g-methods – the g-formula, IP weighting, g-estimation, and doubly robust versions that combine them – recover the correct (null) answer in the same dataset. Throughout, we compare static strategies under the identifiability conditions of Chapter 19: sequential exchangeability, positivity, and consistency.
The running example is the sequentially randomized experiment of Chapter 20 (Table 20.1, reproduced in the book as Table 21.1), with treatments \(A_0, A_1\), a confounder \(L_1\) measured between them, and an outcome \(Y\). There is no \(L_0\) in this study.
| \(N\) | \(A_0\) | \(L_1\) | \(A_1\) | Mean \(Y\) |
|---|---|---|---|---|
| 2400 | 0 | 0 | 0 | 84 |
| 1600 | 0 | 0 | 1 | 84 |
| 2400 | 0 | 1 | 0 | 52 |
| 9600 | 0 | 1 | 1 | 52 |
| 4800 | 1 | 0 | 0 | 76 |
| 3200 | 1 | 0 | 1 | 76 |
| 1600 | 1 | 1 | 0 | 44 |
| 6400 | 1 | 1 | 1 | 44 |
For the effect of the time-fixed treatment \(A_1\) alone, the g-formula is ordinary standardization: \(\operatorname{E}\mathopen{}\left[Y^{a_1}\right]\mathclose{} = \sum_{l_1} \operatorname{E}\mathopen{}\left[Y \mid A_1 = a_1, L_1 = l_1\right]\mathclose{} \Pr[L_1 = l_1]\). For the joint treatment \((A_0, A_1)\), the weights become the distribution of \(L_1\) given past treatment.
Definition 1 (The G-Formula for Time-Varying Treatments) For two time points, the g-formula for the static strategy \((a_0, a_1)\) is
\[\sum_{l_1} \operatorname{E}\mathopen{}\left[Y \mid A_0 = a_0, A_1 = a_1, L_1 = l_1\right]\mathclose{} \, f(l_1 \mid a_0).\]
For \(K+1\) time points \(k = 0, \ldots, K\) and strategy \(\bar{a}\), it is
\[\sum_{\bar{l}} \operatorname{E}\mathopen{}\left[Y \mid \bar{A} = \bar{a}, \bar{L} = \bar{l}\right]\mathclose{} \prod_{k=0}^{K} f\mathopen{}\left(l_k \mid \bar{a}_{k-1}, \bar{l}_{k-1}\right)\mathclose{},\]
where the sum is over all covariate histories \(\bar{l}\). Under sequential exchangeability for \(Y^{\bar{a}}\) given \((\bar{L}_k, \bar{A}_{k-1})\) at each \(k\), positivity, and consistency, the g-formula equals \(\operatorname{E}\mathopen{}\left[Y^{\bar{a}}\right]\mathclose{}\).
Fine Point 21.1: Treatment and Covariate History
The relevant history at time \(k\) is the set of treatments and confounders needed for conditional exchangeability of \(A_k\). Usually this is the chronological past, but confounders can in principle lie in the temporal future of treatment (Fine Point 7.4), and adjusting for some past variables can create M-bias (Figure 7.4) (Hernán and Robins 2020, 280).
Example 1 (G-Formula in the Sequentially Randomized Experiment) Never treat, \((a_0, a_1) = (0, 0)\). Among the 16,000 individuals with \(A_0 = 0\), \(2400 + 1600 = 4000\) have \(L_1 = 0\) and \(2400 + 9600 = 12000\) have \(L_1 = 1\), so \(f(L_1 = 0 \mid A_0 = 0) = 0.25\) and \(f(L_1 = 1 \mid A_0 = 0) = 0.75\). Then
\[ \begin{aligned} \operatorname{E}\mathopen{}\left[Y^{a_0=0, a_1=0}\right]\mathclose{} &= \operatorname{E}\mathopen{}\left[Y \mid 0, 0, L_1 = 0\right]\mathclose{} \times 0.25 + \operatorname{E}\mathopen{}\left[Y \mid 0, 0, L_1 = 1\right]\mathclose{} \times 0.75 \\ &= 84 \times 0.25 + 52 \times 0.75 \\ &= 21 + 39 = 60. \end{aligned} \]
Always treat, \((a_0, a_1) = (1, 1)\). Among the 16,000 individuals with \(A_0 = 1\), \(4800 + 3200 = 8000\) have \(L_1 = 0\), so \(f(L_1 = 0 \mid A_0 = 1) = 0.5\), and
\[ \begin{aligned} \operatorname{E}\mathopen{}\left[Y^{a_0=1, a_1=1}\right]\mathclose{} &= 76 \times 0.5 + 44 \times 0.5 \\ &= 38 + 22 = 60. \end{aligned} \]
The g-formula estimate of \(\operatorname{E}\mathopen{}\left[Y^{1,1}\right]\mathclose{} - \operatorname{E}\mathopen{}\left[Y^{0,0}\right]\mathclose{}\) is \(60 - 60 = 0\), the correct null value (Hernán and Robins 2020, 278).
Remark 1 (The g-formula simulates the counterfactual joint distribution). Under sequential exchangeability for \(Y\) and \(\bar{L}\) jointly, the g-formula simulates the joint distribution of the counterfactuals \((Y^{\bar{a}}, \bar{L}^{\bar{a}})\) that would have been observed had everybody followed strategy \(\bar{a}\). On the tree graph (Figures 21.1-21.2), this means building a new tree in which everyone receives \(a_0\) and \(a_1\) with probability 1, while \(\Pr[L_1 = l_1 \mid A_0 = a_0]\) and \(\operatorname{E}\mathopen{}\left[Y \mid A_0 = a_0, A_1 = a_1, L_1 = l_1\right]\mathclose{}\) keep their observed values.
Definition 2 (Plug-in and parametric g-formula) With many confounders or time points, the components must be estimated, for example
Plugging the estimates into the formula gives the plug-in g-formula; when the estimates come from parametric models, the parametric g-formula.
Technical Point 21.1: The g-formula Density
For a static strategy \(\bar{a}\), the g-formula density of \((Y, \bar{L})\) at \((y, \bar{l})\) is \(f(y \mid \bar{a}_K, \bar{l}_K) \prod_{k=0}^{K} f(l_k \mid \bar{a}_{k-1}, \bar{l}_{k-1})\), and integrating out \(\bar{l}\) gives the g-formula density of \(Y\). For a deterministic dynamic strategy \(g\), replace each \(a_k\) by the value \(a_k^g\) that \(g\) assigns given the simulated history.
For the time-fixed treatment \(A_1\) alone (Chapter 12), the weights are \(W^{A_1} = 1/f(A_1 \mid L_1)\) or \(SW^{A_1} = f(A_1)/f(A_1 \mid L_1)\). For a time-varying treatment, the denominator becomes each individual’s probability of the treatment history they received, given their treatment and covariate history.
Definition 4 (IP Weights for Time-Varying Treatments) With \(A_{-1} \equiv 0\) by definition, the nonstabilized IP weights are
\[W^{\bar{A}} = \prod_{k=0}^{K} \frac{1}{f\mathopen{}\left(A_k \mid \bar{A}_{k-1}, \bar{L}_k\right)\mathclose{}},\]
and the stabilized IP weights are
\[SW^{\bar{A}} = \prod_{k=0}^{K} \frac{f\mathopen{}\left(A_k \mid \bar{A}_{k-1}\right)\mathclose{}}{f\mathopen{}\left(A_k \mid \bar{A}_{k-1}, \bar{L}_k\right)\mathclose{}}.\]
For \(K = 1\): \(W^{\bar{A}} = \dfrac{1}{f(A_0 \mid L_0)} \times \dfrac{1}{f(A_1 \mid A_0, L_0, L_1)}\).
Proposition 1 (IP weighting identifies the counterfactual mean) Under the identifiability conditions (sequential exchangeability, positivity, and consistency), \(\operatorname{E}\mathopen{}\left[Y^{a_0, a_1}\right]\mathclose{}\) equals the mean \(\operatorname{E}_{ps}\mathopen{}\left[Y \mid A_0 = a_0, A_1 = a_1\right]\mathclose{}\) in the pseudo-population created by either set of weights in Definition 4.
Example 2 (IP Weighting in the Sequentially Randomized Experiment) There is no \(L_0\), so the denominator is \(f(A_0) f(A_1 \mid A_0, L_1)\), with \(f(A_0 = 0) = 16000/32000 = 0.5\). For the never-treated rows of Table 1:
The 32,000 pseudo-individuals with \((A_0, A_1) = (0, 0)\) give
\[ \begin{aligned} \operatorname{E}_{ps}\mathopen{}\left[Y \mid A_0 = 0, A_1 = 0\right]\mathclose{} &= 84 \times \frac{8000}{32000} + 52 \times \frac{24000}{32000} \\ &= 21 + 39 = 60. \end{aligned} \]
The same calculation for \((1, 1)\) also gives 60, so the IP weighted effect estimate is 0. The whole pseudo-population has \(4 \times 32000 = 128000\) individuals, the 32,000 study individuals times the 4 static strategies (Hernán and Robins 2020, 282–83).
Definition 5 (Marginal structural mean model for a time-varying treatment) There are far more strategies \(\bar{a}\) than can be estimated one at a time, so we pool information with a marginal structural mean model, for example
\[\operatorname{E}\mathopen{}\left[Y^{\bar{a}}\right]\mathclose{} = \beta_0 + \beta_1 \, \mathrm{cum}(\bar{a}), \qquad \mathrm{cum}(\bar{a}) = \sum_{k=0}^{K} a_k.\]
This model is unsaturated: one unknown per strategy on the left, two parameters on the right. The average causal effect is \(\operatorname{E}\mathopen{}\left[Y^{\bar{a}}\right]\mathclose{} - \operatorname{E}\mathopen{}\left[Y^{\bar{a} = \bar{0}}\right]\mathclose{} = \beta_1 \times \mathrm{cum}(\bar{a})\).
Proposition 2 (IP weighted least squares is consistent for the marginal structural model) Suppose the marginal structural model of Definition 5 is correctly specified, and the identifiability conditions (sequential exchangeability, positivity, and consistency) hold. Let \(\hat \beta_1\) be the coefficient of \(\mathrm{cum}(\bar{A})\) in the fit of the model \(\operatorname{E}\mathopen{}\left[Y \mid \bar{A}\right]\mathclose{} = \theta_0 + \theta_1 \, \mathrm{cum}(\bar{A})\) by weighted least squares with weights \(SW^{\bar{A}}\) or \(W^{\bar{A}}\) of Definition 4 (known, or estimated from a correctly specified treatment model). Then \(\hat \beta_1\) is consistent for the causal \(\beta_1\), which in general differs from the associational \(\theta_1\) of the unweighted data.
Definition 6 (Marginal structural model with a baseline effect modifier) Effect modification: for a dichotomous baseline variable \(V\) in \(L_0\), a marginal structural model with effect modification by \(V\) is \(\operatorname{E}\mathopen{}\left[Y^{\bar{a}} \mid V\right]\mathclose{} = \beta_0 + \beta_1 \mathrm{cum}(\bar{a}) + \beta_2 V + \beta_3 \mathrm{cum}(\bar{a}) V\). Its parameters are estimated by fitting the associational model \(\operatorname{E}\mathopen{}\left[Y \mid \bar{A}, V\right]\mathclose{} = \theta_0 + \theta_1 \mathrm{cum}(\bar{A}) + \theta_2 V + \theta_3 \mathrm{cum}(\bar{A}) V\) by weighted least squares with \(W^{\bar{A}}\) or, better, \(SW^{\bar{A}}(V) = \prod_{k=0}^{K} f(A_k \mid \bar{A}_{k-1}, V) / f(A_k \mid \bar{A}_{k-1}, \bar{L}_k)\).
Proposition 3 (IP weighted least squares is consistent for the effect-modification model) Suppose the marginal structural model of Definition 6 is correctly specified, the identifiability conditions (sequential exchangeability, positivity, and consistency) hold, and \(V\) is a baseline variable in \(L_0\). Let \(\hat \beta_1\) and \(\hat \beta_3\) be the coefficients of \(\mathrm{cum}(\bar{A})\) and \(\mathrm{cum}(\bar{A}) V\) in the fit of the associational model of Definition 6 by weighted least squares with weights \(SW^{\bar{A}}(V)\) of Definition 6 or \(W^{\bar{A}}\) of Definition 4 (known, or estimated from a correctly specified model for treatment in the denominator). Then \(\hat \beta_1\) and \(\hat \beta_3\) are consistent for the causal \(\beta_1\) and \(\beta_3\).
Technical Point 21.2: IP Weighting for Dynamic Treatment Strategies
For a deterministic dynamic strategy \(g\), the g-formula equals the mean of \(Y\) among pseudo-population members who follow \(g\) in the pseudo-population created by nonstabilized weights. Stabilized weights cannot be used, because their numerator depends on \(A\).
Technical Point 21.3: The g-null Paradox
Under the sharp null with treatment-confounder feedback, the g-formula does not depend on \(\bar{a}\), but non-saturated parametric models for \(\operatorname{E}\mathopen{}\left[Y \mid \bar{A}, \bar{L}\right]\mathclose{}\) and \(f(l_k \mid \cdot)\) with variation-independent parameters cannot all be correctly specified when \(L_k\) has a discrete component (Robins and Wasserman 1997). The estimated g-formula may therefore falsely reject the null even in a sequentially randomized experiment. In practice the bias seems small relative to random variability. MSMs and structural nested mean models do not suffer from this paradox: both are correctly specified under the null whatever functional form we choose (Hernán and Robins 2020, 286).
Definition 7 (Doubly robust estimator) IP weighting needs a correct model for treatment given confounders; the g-formula needs a correct model for the outcome given treatment and confounders. A doubly robust estimator is consistent if either model is correct, without our knowing which, giving “two chances to get it right” (Hernán and Robins 2020, 286).
Algorithm 1 (Doubly robust plug-in estimator for a time-fixed treatment) For a binary \(A\), binary \(Y\), and many confounders \(L\), estimate \(\operatorname{E}\mathopen{}\left[Y^a\right]\mathclose{}\) in three steps:
The difference of the step-3 averages for \(a = 1\) and \(a = 0\) is a doubly robust estimator of \(\operatorname{E}\mathopen{}\left[Y^{a=1}\right]\mathclose{} - \operatorname{E}\mathopen{}\left[Y^{a=0}\right]\mathclose{}\).
Algorithm 2 (Doubly robust plug-in estimator for a time-varying treatment) To estimate \(\operatorname{E}\mathopen{}\left[Y^{\bar{a} = \bar{1}}\right]\mathclose{}\) (“always treated”):
Repeat with \(\bar{a} = \bar{0}\) for “never treated”; the difference estimates the average causal effect.
Proposition 4 (The time-varying doubly robust estimator is \(K + 2\) robust) Suppose each sequential outcome model \(b_m\) of alg. 2 is fit by solving the score equations of a generalized linear model with a canonical link (logistic for a binary \(Y\), linear for a continuous \(Y\)). For a binary \(Y\) and \(m < K\), the response \(\hat{B}_{m+1}\) lies in \([0, 1]\) rather than \(\{0, 1\}\), so the logistic fit is a quasi-likelihood fit. Suppose also that each \(b_m\) includes the weight \(\hat{W}^{\bar{1}_m}\) as a covariate. Suppose further that positivity holds for the strategy \(\bar{a} = \bar{1}\). Then the estimator of alg. 2 is consistent for the g-formula functional of Definition 1 at \(\bar{a} = \bar{1}\) if (i) the outcome models \(b_m\) are correct for all \(m\), or (ii) the treatment models \(\pi_k\) are correct for all \(k\). In fact it is \(K + 2\) robust: it is also consistent if, for some \(m \in \{0, \ldots, K - 1\}\), the treatment model is correct for times \(0\) to \(m\) and the outcome model is correct for times \(m + 1\) to \(K\) (Hernán and Robins 2020, 288). If, in addition, sequential exchangeability and consistency hold, that functional equals \(\operatorname{E}\mathopen{}\left[Y^{\bar{a} = \bar{1}}\right]\mathclose{}\), so the estimator is consistent for the counterfactual mean.
Technical Point 21.4: A K+2 Robust Augmented IP Weighted Estimator
For \(K = 1\), let \(\hat{b}_1(L_0, L_1) = \mathop{\hat{\operatorname{E}}}\nolimits\mathopen{}\left[Y \mid L_0, A_0 = 1, L_1, A_1 = 1\right]\mathclose{}\) and \(\hat{b}_0(L_0) = \mathop{\hat{\operatorname{E}}}\nolimits\mathopen{}\left[\hat{b}_1(L_0, L_1) \mid A_0 = 1, L_0\right]\mathclose{}\) (the iterated conditional expectation, ICE, plug-in estimator of the g-formula averages \(\hat{b}_0\)), and let \(\hat{\pi}_0, \hat{\pi}_1\) estimate \(\Pr[A_0 = 1 \mid L_0]\) and \(\Pr[A_1 = 1 \mid L_0, L_1, A_0 = 1]\). The augmented IP weighted estimator of Robins et al. (1994) is the sample average of
\[\hat{U}_{TR} = \hat{b}_0(L_0) + \frac{A_0 A_1}{\hat{\pi}_0 \hat{\pi}_1}\mathopen{}\left(Y - \hat{b}_1(L_0, L_1)\right)\mathclose{} + \frac{A_0}{\hat{\pi}_0}\mathopen{}\left(\hat{b}_1(L_0, L_1) - \hat{b}_0(L_0)\right)\mathclose{}.\]
It is consistent if (a) \(\hat{\pi}_0, \hat{\pi}_1\) are consistent (the augmentation terms then average to 0, leaving the IP weighted estimator), (b) both outcome regressions are consistent (the last two terms then average to 0, leaving the ICE estimator), or (c) \(\hat{b}_1\) and \(\hat{\pi}_0\) are consistent (Molina et al. 2017): hence “triply”, that is \(K + 2\), robust. A modification due to Tchetgen Tchetgen (2009) is \(2^{K+1}\) (quadruply) robust. Technical Point 21.5 shows that including the terms \(A_0 A_1/(\hat{\pi}_0 \hat{\pi}_1)\) and \(A_0/\hat{\pi}_0\) as covariates in the sequential logistic outcome models makes both correction terms average to exactly 0, so the plug-in estimate stays in \([0, 1]\); the estimator in the main text is an instance of this construction.
Technical Point 21.5: A Plug-In K+2 Robust Estimator
With a binary \(Y\), the sample average of \(\hat{U}_{TR}\) from Technical Point 21.4 can fall outside \([0, 1]\). The ICE estimator, the sample average of \(\hat{b}_0(L_0)\), cannot, as long as both outcome regressions are logistic: it is a plug-in estimator. To make it triply robust as well, fit the outcome models so that both correction terms of \(\hat{U}_{TR}\) have sample mean exactly 0; the average of \(\hat{U}_{TR}\) then coincides with the average of \(\hat{b}_0(L_0)\).
The doubly robust estimator of alg. 2, which uses the time-varying weight as a covariate, is an instance of this plug-in estimator (Hernán and Robins 2020, 289).
Fine Point 21.2: Representations of the g-formula
The same g-formula can be written as standardization (the joint density modeling estimator of Section 21.1), as iterated conditional expectations (the ICE estimator, the basis of this section), or as an IP weighted mean (Section 21.2). Without parametric restrictions the three coincide, but in practice each suggests a different estimator (Hernán and Robins 2020, 291).
Technical Point 21.6: A Multiply Robust Estimator
This technical point gives a multiply robust plug-in estimator of \(\operatorname{E}\mathopen{}\left[Y^g\right]\mathclose{}\) for static, dynamic, or random strategies \(g\), based on Rotnitzky et al. (2017); without the weight covariate, the same algorithm is the ICE estimator of the g-formula. The book expects such estimators, fit with machine learning and sample splitting, to become more common as software matures, especially for failure-time outcomes.
Definition 8 (Structural Nested Mean Model) A structural nested mean model (SNMM) for a time-varying treatment specifies, for each \(k = 0, \ldots, K\),
\[\operatorname{E}\mathopen{}\left[Y^{\bar{a}_{k-1}, a_k, \underline{0}_{k+1}} - Y^{\bar{a}_{k-1}, \underline{0}_k} \mid \bar{L}_k^{\bar{a}_{k-1}} = \bar{l}_k, \bar{A}_{k-1} = \bar{a}_{k-1}, A_k = a_k\right]\mathclose{} = a_k \, \gamma_k\mathopen{}\left(\bar{a}_{k-1}, \bar{l}_k, \beta\right)\mathclose{},\]
where \((\bar{a}_{k-1}, a_k, \underline{0}_{k+1})\) follows \(\bar{a}_{k-1}\) through \(k - 1\), gives \(a_k\) at \(k\), and gives no treatment from \(k + 1\) to \(K\); \(\gamma_k(\bar{a}_{k-1}, \bar{l}_k, \psi^\dagger)\) is a known function of past treatment, covariate history, and a parameter value \(\psi^\dagger\), with \(\gamma_k(\bar{a}_{k-1}, \bar{l}_k, 0) = 0\); and \(\beta\) is the unknown true value of that parameter.
Remark 3 (What an SNMM models). An SNMM models the effect on the mean of \(Y\) of a last blip of treatment \(a_k\) at time \(k\), as a function of past treatment and covariate history. Because \(\gamma_k(\cdot, 0) = 0\), the null hypothesis of no effect corresponds to \(\beta = 0\).
Example 3 (A saturated SNMM for two time points) An SNMM (Definition 8) has one equation per time point. For Table 1 (\(K = 1\)), a saturated additive SNMM is
\[ \begin{aligned} k = 0: &\quad \operatorname{E}\mathopen{}\left[Y^{a_0, a_1 = 0} - Y^{a_0 = 0, a_1 = 0} \mid A_0 = a_0\right]\mathclose{} = \beta_0 a_0, \\ k = 1: &\quad \operatorname{E}\mathopen{}\left[Y^{a_0, a_1} - Y^{a_0, a_1 = 0} \mid L_1^{a_0} = l_1, A_0 = a_0, A_1^{a_0} = a_1\right]\mathclose{} = a_1 \mathopen{}\left(\beta_{11} + \beta_{12} l_1 + \beta_{13} a_0 + \beta_{14} a_0 l_1\right)\mathclose{}. \end{aligned} \]
The first equation is the effect of treating at time 0 and never again; the second, of treating at time 1 and never again. Under sequential exchangeability, \(A_0 = a_0\) and \(A_1^{a_0} = a_1\) can be dropped from the conditioning events.
Definition 9 (Additive rank-preserving model for two time points) As in Chapter 14, g-estimation starts from the additive rank-preserving model for each individual \(i\):
\[ \begin{aligned} Y_i^{a_0, 0} &= Y_i^{0, 0} + \psi_0 a_0, \\ Y_i^{a_0, a_1} &= Y_i^{a_0, 0} + \psi_{11} a_1 + \psi_{12} a_1 L_{1,i}^{a_0} + \psi_{13} a_1 a_0 + \psi_{14} a_1 a_0 L_{1,i}^{a_0}. \end{aligned} \]
The second equation is only locally rank-preserving: ranks are preserved among individuals with the same \(a_0\) and \(l_1\).
Rank Preservation Is Implausible but Not Required
Rank preservation is biologically implausible, but it is not needed. G-estimation (the Chapter 14 method, extended here to time-varying treatment) is consistent for the SNMM parameters \(\beta\) of Definition 8 even if the rank-preserving model is wrong, provided the SNMM is correctly specified, sequential exchangeability for \(Y\) and consistency hold, and the model for treatment given past treatment and covariates is correctly specified.
Definition 10 (Candidate counterfactuals for two time points) By consistency (\(Y^{A_0, A_1} = Y\) and \(L_1^{A_0} = L_1\)), the model can be written in terms of observed data. Replacing \(\psi\) by a candidate value \(\psi^\dagger\) defines the candidate counterfactuals
\[ \begin{aligned} H_1(\psi^\dagger) &= Y - \mathopen{}\left(\psi_{11}^\dagger A_1 + \psi_{12}^\dagger A_1 L_1 + \psi_{13}^\dagger A_1 A_0 + \psi_{14}^\dagger A_1 A_0 L_1\right)\mathclose{}, \\ H_0(\psi^\dagger) &= H_1(\psi^\dagger) - \psi_0^\dagger A_0. \end{aligned} \]
Remark 4 (Candidate counterfactuals at the true parameter). When \(\psi^\dagger = \psi\), \(H_k(\psi^\dagger)\) equals the counterfactual \(Y^{\bar{A}_{k-1}, \bar{0}_k}\): observed treatment through \(k - 1\) and no treatment afterwards. Sequential exchangeability then says \(H_k(\psi)\) is independent of \(A_k\) given the past.
Fine Point 21.3: G-Estimation with a Saturated Structural Nested Model
With a saturated SNMM, g-estimation needs no search: at the true \(\psi\), sequential exchangeability forces the mean of \(H_k(\psi)\) to agree between those treated and those untreated at time \(k\) within every stratum of the past. Each such equality is one linear equation in the parameters, and solving them from the last time point backwards recovers \(\psi\) (Hernán and Robins 2020, 293). Example 4 below carries this out for Table 1.
Example 4 (G-Estimation in the Sequentially Randomized Experiment) Time 1. Within each stratum of \((A_0, L_1)\), the mean of \(H_1(\psi)\) must not depend on \(A_1\). The means of \(H_1(\psi)\) by row of Table 1 are \(84\) and \(84 - \psi_{11}\) in stratum \((0, 0)\), \(52\) and \(52 - \psi_{11} - \psi_{12}\) in \((0, 1)\), \(76\) and \(76 - \psi_{11} - \psi_{13}\) in \((1, 0)\), and \(44\) and \(44 - \psi_{11} - \psi_{12} - \psi_{13} - \psi_{14}\) in \((1, 1)\) (for \(A_1 = 0\) and \(A_1 = 1\) respectively). Equating within strata:
\[ \begin{aligned} 84 = 84 - \psi_{11} &\implies \psi_{11} = 0, \\ 52 = 52 - \psi_{11} - \psi_{12} = 52 - \psi_{12} &\implies \psi_{12} = 0, \\ 76 = 76 - \psi_{11} - \psi_{13} = 76 - \psi_{13} &\implies \psi_{13} = 0, \\ 44 = 44 - \psi_{11} - \psi_{12} - \psi_{13} - \psi_{14} = 44 - \psi_{14} &\implies \psi_{14} = 0. \end{aligned} \]
Time 0. With all \(\psi_{1j} = 0\), \(H_1(\psi) = Y\) and \(H_0(\psi) = Y - \psi_0 A_0\). Exchangeability \(Y^{0,0} \perp\!\!\!\perp A_0\) requires equal means of \(H_0(\psi)\) in the two arms of 16,000 each:
\[ \begin{aligned} A_0 = 0: &\quad 84 \times 0.25 + 52 \times 0.75 = 60, \\ A_0 = 1: &\quad (76 - \psi_0) \times 0.5 + (44 - \psi_0) \times 0.5 = 60 - \psi_0, \end{aligned} \]
so \(60 = 60 - \psi_0\) and \(\psi_0 = 0\). All parameters are 0, so \(\operatorname{E}\mathopen{}\left[Y^g\right]\mathclose{} = 60\) under every static or dynamic strategy \(g\), in agreement with the g-formula and IP weighting (Hernán and Robins 2020, 292–93).
Definition 11 (Candidate counterfactuals for a general SNMM) The candidate counterfactuals generalize to
\[H_k(\psi^\dagger) = Y - \sum_{j=k}^{K} A_j \, \gamma_j\mathopen{}\left(\bar{A}_{j-1}, \bar{L}_j, \psi^\dagger\right)\mathclose{}.\]
Technical Point 21.7: Marginal Structural Models and Structural Nested Models
An SNMM is a semiparametric MSM if and only if \(\gamma_k\) does not depend on \(\bar{l}_k\); the implied MSM is \(\operatorname{E}\mathopen{}\left[Y^{\bar{a}}\right]\mathclose{} = \alpha_0 + \sum_{k=0}^{K} a_k \gamma_k(\bar{a}_{k-1}, \beta)\) with \(\alpha_0 = \operatorname{E}\mathopen{}\left[Y^{\bar{0}}\right]\mathclose{}\). Such an SNMM additionally assumes no effect modification by past covariates. If that assumption holds, g-estimation is more efficient than IP weighting; if it fails, g-estimation is biased while IP weighting of the MSM is not: a variance-bias trade-off.
Algorithm 3 (G-estimation of a one-parameter SNMM) Take \(\gamma_k(\bar{a}_{k-1}, \bar{l}_k, \psi^\dagger) = \psi^\dagger\) (true value \(\beta_1\)), so \(H_k(\psi^\dagger) = Y - \sum_{j=k}^{K} A_j \psi^\dagger\).
Remark 5 (G-estimation as solving an estimating equation). Equivalently, \(\hat \beta\) is the \(\psi^\dagger\) at which the score test of \(\alpha_1 = 0\) has \(P = 1\), that is, the solution of
\[\sum_{i=1}^{N} \sum_{k=0}^{K} \mathopen{}\left\{A_{i,k} - \operatorname{expit}\mathopen{}\left(\hat{\alpha}_0 + \hat{\alpha}_2 W_{i,k}\right)\mathclose{}\right\}\mathclose{} H_{i,k}(\psi^\dagger) = 0,\]
with \(\hat{\alpha}_0, \hat{\alpha}_2\) from the model fit with \(\alpha_1 = 0\).
Proposition 6 (Consistency of the g-estimator) Let \(\hat \beta\) be the solution of the estimating equation of Remark 5 (equivalently, the limit of the grid-search estimate of alg. 3 as the grid step goes to 0). Then \(\hat \beta\) is consistent for the true \(\beta_1\) if (i) the one-parameter SNMM of alg. 3, \(\gamma_k(\bar{a}_{k-1}, \bar{l}_k, \beta) = \beta_1\), is correct, (ii) sequential exchangeability for \(Y\) holds, (iii) consistency holds, \(Y^{\bar{A}} = Y\), (iv) the treatment model \(\operatorname{logit}\Pr[A_k = 1 \mid \bar{L}_k, \bar{A}_{k-1}] = \alpha_0 + \alpha_2 W_k\) is correct, and (v) \(H_k(\psi^\dagger)\) enters the treatment model linearly (Technical Point 14.2).
Technical Point 21.8: A Closed Form Estimator for Linear Structural Nested Mean Models
Suppose the blip is linear in \(\beta\), \(\gamma_k(\bar{a}_{k-1}, \bar{l}_k, \beta) = {\beta}^{\top} R_k\), where \(R_k = r_k(\bar{L}_k, \bar{A}_{k-1})\) is a vector of known functions, and the treatment model is \(\operatorname{logit}\Pr[A_k = 1 \mid \bar{L}_k, \bar{A}_{k-1}] = {\alpha}^{\top} W_k\). Then \(H_k(\beta) = Y - {\beta}^{\top} S_k\) with \(S_k \stackrel{\text{def}}{=}\sum_{j=k}^{K} A_j R_j\), so the g-estimating equation
\[\sum_{i=1}^{N} \sum_{k=0}^{K} X_{i,k}(\hat{\alpha})\, Q_{i,k}\, \mathopen{}\left(Y_i - {S_{i,k}}^{\top} \beta\right)\mathclose{} = 0, \qquad X_{i,k}(\hat{\alpha}) \stackrel{\text{def}}{=}A_{i,k} - \operatorname{expit}\mathopen{}\left({\hat{\alpha}}^{\top} W_{i,k}\right)\mathclose{},\]
is linear in \(\beta\) and, when the matrix below is invertible, has the explicit solution
\[\hat \beta= \mathopen{}\left(\sum_{i=1}^{N} \sum_{k=0}^{K} X_{i,k}(\hat{\alpha})\, Q_{i,k}\, {S_{i,k}}^{\top}\right)\mathclose{}^{-1} \sum_{i=1}^{N} \sum_{k=0}^{K} Y_i\, X_{i,k}(\hat{\alpha})\, Q_{i,k}.\]
Here \(Q_{i,k} = q_k(\bar{L}_{i,k}, \bar{A}_{i,k-1})\) has the dimension of \(\beta\); its choice changes the efficiency of \(\hat \beta\) but not its consistency (Robins 1994 gives the optimal choice).
A multiply robust version adds a working model \({\varsigma}^{\top} D_k\), with \(D_k = d_k(\bar{L}_k, \bar{A}_{k-1})\), for \(\operatorname{E}\mathopen{}\left[H_k(\beta) \mid \bar{L}_k, \bar{A}_{k-1}\right]\mathclose{} = \operatorname{E}\mathopen{}\left[Y^{\bar{A}_{k-1}, \underline{0}_k} \mid \bar{L}_k, \bar{A}_{k-1}\right]\mathclose{}\), and solves jointly for \((\tilde{\beta}, \tilde{\varsigma})\)
\[\sum_{i,k} X_{i,k}(\hat{\alpha})\, Q_{i,k}\, \mathopen{}\left(H_{i,k}(\beta) - {\varsigma}^{\top} D_{i,k}\right)\mathclose{} = 0, \qquad \sum_{i,k} D_{i,k}\, \mathopen{}\left(H_{i,k}(\beta) - {\varsigma}^{\top} D_{i,k}\right)\mathclose{} = 0,\]
which are again linear, so the solution is in closed form. \(\tilde{\beta}\) is consistent and asymptotically normal if, at each \(k\), either the working outcome model or the treatment model is correct, which makes it \(2^{K+1}\) multiply robust (Hernán and Robins 2020, 296).
Remark 7 (Counterfactual means after g-estimation). Given \(\hat \beta\), how to estimate a counterfactual mean depends on the strategy and the blip function:
Technical Point 21.9: Estimating a Counterfactual Mean After g-estimation of a Structural Nested Mean Model
Estimate \(\operatorname{E}\mathopen{}\left[Y^{\bar{0}}\right]\mathclose{}\) by the average of \(H_0(\hat \beta)\); fit a model for \(f(l_k \mid \bar{a}_{k-1}, \bar{l}_{k-1})\); simulate covariate histories under \(g\); and add to \(\mathop{\hat{\operatorname{E}}}\nolimits\mathopen{}\left[Y^{\bar{0}}\right]\mathclose{}\) the average over simulations of \(\sum_{j=0}^{K} a_{v,j} \gamma_j(\bar{a}_{v,j-1}, \bar{l}_{v,j}, \hat \beta)\). The result is consistent for \(\operatorname{E}\mathopen{}\left[Y^g\right]\mathclose{}\) if the covariate model and the SNMM are correct and either the treatment model or the outcome model \(\operatorname{E}\mathopen{}\left[Y^{\bar{A}_{k-1}, \underline{0}_k} \mid \bar{L}_k, \bar{A}_{k-1}\right]\mathclose{}\) is correct. Under the null, \(\gamma_j(\cdot, \hat \beta)\) converges to 0 whenever \(\hat \beta\) is consistent for \(\beta = 0\), even when the covariate model is wrong. So, under the identifiability conditions, this procedure preserves the null, unlike the parametric g-formula, provided \(\Pr[A_k = 1 \mid \bar{L}_k, \bar{A}_{k-1}]\) is known (as in a sequentially randomized experiment) or either the treatment model or the outcome model is correct at each \(k\) (Hernán and Robins 2020, 297).
Technical Point 21.13: Formal Definition of a General Structural Nested Mean Model
This technical point generalizes SNMMs to blips followed by an arbitrary strategy \(g\) rather than “never treat” (Robins 2004), and notes identification under a time-varying parallel-trends assumption instead of sequential exchangeability (Hernán and Robins 2020, 303).
Part II treated censoring \(C\) as time-fixed.
Definition 12 (Time-varying censoring) More realistically, censoring is time-varying: \(C_1, C_2, \ldots, C_{K+1}\), where \(C_m = 0\) if the individual is still uncensored at \(m\) and 1 otherwise.
Definition 13 (The g-formula with time-varying censoring) With time-varying censoring, the g-formula of Definition 1 for strategy \(\bar{a}\) becomes
\[\sum_{\bar{l}} \operatorname{E}\mathopen{}\left[Y \mid \bar{C} = \bar{0}, \bar{A} = \bar{a}, \bar{L} = \bar{l}\right]\mathclose{} \prod_{k=0}^{K} f\mathopen{}\left(l_k \mid \bar{c}_k = \bar{0}, \bar{a}_{k-1}, \bar{l}_{k-1}\right)\mathclose{}.\]
Proposition 7 (The censored g-formula identifies the uncensored counterfactual mean) If sequential exchangeability, positivity, and consistency hold with \(A_m\) replaced by the joint treatment \((A_m, C_{m+1})\) at every \(m\), then the g-formula of Definition 13 equals \(\operatorname{E}\mathopen{}\left[Y^{\bar{a}, \bar{c} = \bar{0}}\right]\mathclose{}\): the mean outcome had everyone followed \(\bar{a}\) and nobody been lost to follow-up.
Algorithm 4 (IP weighting with time-varying censoring) Fit, for example, \(\operatorname{E}\mathopen{}\left[Y \mid \bar{A}, \bar{C} = \bar{0}\right]\mathclose{} = \theta_0 + \theta_1 \mathrm{cum}(\bar{A})\) in the pseudo-population created by \(W^{\bar{A}} \times W^{\bar{C}}\) or \(SW^{\bar{A}} \times SW^{\bar{C}}\), where
\[W^{\bar{C}} = \prod_{k=1}^{K+1} \frac{1}{\Pr\mathopen{}\left[C_k = 0 \mid C_{k-1} = 0, \bar{A}_{k-1}, \bar{L}_{k-1}\right]\mathclose{}}, \qquad SW^{\bar{C}} = \prod_{k=1}^{K+1} \frac{\Pr\mathopen{}\left[C_k = 0 \mid C_{k-1} = 0, \bar{A}_{k-1}\right]\mathclose{}}{\Pr\mathopen{}\left[C_k = 0 \mid C_{k-1} = 0, \bar{A}_{k-1}, \bar{L}_{k-1}\right]\mathclose{}}.\]
Numerator and denominator are estimated with separate logistic models.
Remark 8 (The pseudo-population created by censoring weights).
Remark 9 (G-estimation with censoring). G-estimation with censoring: first create, with nonstabilized weights \(W^{\bar{C}}\), a pseudo-population in which nobody is censored, then g-estimate in it.
Technical Point 21.10: Survival Analysis with Time-Varying Treatments
To estimate the risk \(\Pr[D_{k+1}^{\bar{a}, \bar{c} = \bar{0}} = 1]\), either use the g-formula with models for the discrete-time hazards and for the confounder densities, conditional on being event-free and uncensored (e.g., pooled logistic models, as in Chapter 17), or fit a pooled logistic hazards model in which each person-time receives the time-varying weight \(W_k^{\bar{A}} \times W_k^{\bar{C}}\) (or its stabilized version); the latter estimates a marginal structural pooled logistic model (Robins 1998) (Hernán and Robins 2020, 299).
This part of the book relies on sequential exchangeability given the measured \(L\) and identification by the g-formula, because few realistic longitudinal analyses use other identifying conditions (such as the front door criterion). But all identifying formulas are mathematically linked to the g-formula based on all variables, measured and unmeasured.
Definition 14 (The Big G-Formula) Given a causal DAG with treatment \(A\), outcome \(Y\), measured covariates \(L\), and unmeasured variables \(U\), let \(X = (L, U)\). Every parent of a treatment node is in \(A\) or \(X\), so \(X\) ensures sequential exchangeability. The big g-formula is the g-formula with \(L\) replaced by \(X\). If positivity held, it would identify \(\operatorname{E}\mathopen{}\left[Y^g\right]\mathclose{}\) under any strategy \(g\) (Hernán and Robins 2020, 300).
Definition 15 (Factuals) Factuals are variables of the actual world, in contrast to counterfactuals; some factuals, like \(U\), are unavailable for analysis because they were not measured (Hernán and Robins 2020, 300).
Remark 10 (The big g-formula cannot be computed directly). The big g-formula is a function of the distribution of the factuals \((A, L, Y, U)\) (Definition 15) alone, but it includes unmeasured \(U\), so it cannot be computed directly.
Remark 11 (When the big g-formula reduces to an observed-data formula). The key question: can the big g-formula be rewritten as a functional of the distribution of the observed \((A, L, Y)\) only? If so, the new formula reproduces the big g-formula and can be used in data analysis.
A Reduction Is Causal Only If the DAG Is
These are purely mathematical questions about distributions obeying d-separation; they acquire a causal meaning only if the DAG is truly causal, which we can never know for certain in an observational study.
Example 6 (The Front Door Formula from the Big G-Formula) Under Figure 7.14 (\(A \to M \to Y\), with \(U\) a common cause of \(A\) and \(Y\)), the big g-formula for \(\Pr[Y^a = y]\) is
\[\sum_m \sum_u \Pr[Y = y \mid M = m, A = a, U = u] \Pr[M = m \mid A = a, U = u] \Pr[U = u].\]
Step 1: by \(A \perp\!\!\!\perp Y \mid M, U\), drop \(A = a\) from the first factor; by \(U \perp\!\!\!\perp M \mid A\), replace \(\Pr[M = m \mid A = a, U = u]\) by \(\Pr[M = m \mid A = a]\); and write \(\Pr[U = u] = \sum_{a'} \Pr[U = u \mid A = a'] \Pr[A = a']\):
\[\sum_m \Pr[M = m \mid A = a] \sum_u \Pr[Y = y \mid M = m, U = u] \sum_{a'} \Pr[U = u \mid A = a'] \Pr[A = a'].\]
Step 2: by \(U \perp\!\!\!\perp M \mid A\), \(\Pr[U = u \mid A = a'] = \Pr[U = u \mid M = m, A = a']\), and by \(A \perp\!\!\!\perp Y \mid M, U\), \(\Pr[Y = y \mid M = m, U = u] = \Pr[Y = y \mid M = m, A = a', U = u]\):
\[\sum_m \Pr[M = m \mid A = a] \sum_{a'} \mathopen{}\left\{\sum_u \Pr[Y = y \mid M = m, A = a', U = u] \Pr[U = u \mid M = m, A = a']\right\}\mathclose{} \Pr[A = a'].\]
Step 3: the braces sum over \(u\) to \(\Pr[Y = y \mid M = m, A = a']\), giving the front door formula
\[\Pr[Y^a = y] = \sum_m \Pr[M = m \mid A = a] \sum_{a'} \Pr[Y = y \mid M = m, A = a'] \Pr[A = a'],\]
which involves observed variables only. This proof does not require that counterfactuals \(Y^m\) exist (Hernán and Robins 2020, 301).
Technical Point 21.11: A Big G-Formula Proof of the Front Door Formula
Technical Point 7.4 derived the front door formula using counterfactuals \(Y^m\). Starting instead from the big g-formula, the conditional independencies of Figure 7.14 alone reduce it to observed data (Example 6), so the front door formula holds even if no well-defined \(Y^m\) exists.
A second route, also free of \(Y^m\), is a coupling argument. Even if everyone agrees that \(Y^m\) is not well defined, any law of the observed data that factorizes according to Figure 7.14 can be generated by some FFRCISTG model (Technical Point 6.2) that is “as detailed as the data” in the sense of Robins and Richardson (2010), and such a model formally includes a variable \(Y^m\). Within that model, Technical Point 7.4 shows that the two formulas coincide. A factorizing law for which they differed could therefore not be generated by any such model, so no such distribution exists (Hernán and Robins 2020, 301).
Technical Point 21.12: A Front Door Formula Proof Using d-Separation of Treatment Nodes on SWIGs
A SWIG property (Shpitser et al. 2022). Let \(G(a)\) be the SWIG for strategy \(a\), assume only treatment counterfactuals are well defined, and let \(B^a\) and \(C^a\) be disjoint sets of random (not fixed) observed nodes of \(G(a)\), such as \(Y^a\), \(M^a\), or the natural value of treatment \(A\). If, in \(G(a)\), \(B^a\) and the fixed node \(a_m\) are d-separated given \(C^a\), then \(\Pr[B^a = b \mid C^a = c]\) does not depend on \(a_m\). This compares distributions from different single worlds; it is not a cross-world statement.
Front door example. In the SWIG for Figure 7.14, let \(C^a = (M^a, A)\) and \(B^a = Y^a\). The only path from the fixed node \(a\) to \(Y^a\) passes through the non-collider \(M^a\), which is conditioned on, so \(\operatorname{E}\mathopen{}\left[Y^a \mid M^a, A\right]\mathclose{} = \operatorname{E}\mathopen{}\left[Y^{a'} \mid M^{a'}, A\right]\mathclose{}\) for all \(a, a'\).
Proof of the front door formula. Following Technical Point 7.4, it remains to show \(\operatorname{E}\mathopen{}\left[Y^a \mid M^a = m\right]\mathclose{} = \sum_{a'} \operatorname{E}\mathopen{}\left[Y \mid M = m, A = a'\right]\mathclose{} \Pr[A = a']\):
\[ \begin{aligned} \operatorname{E}\mathopen{}\left[Y^a \mid M^a = m\right]\mathclose{} &= \sum_{a'} \operatorname{E}\mathopen{}\left[Y^a \mid M^a = m, A = a'\right]\mathclose{} \Pr[A = a' \mid M^a = m] \\ &= \sum_{a'} \operatorname{E}\mathopen{}\left[Y^a \mid M^a = m, A = a'\right]\mathclose{} \Pr[A = a'] \\ &= \sum_{a'} \operatorname{E}\mathopen{}\left[Y^{a'} \mid M^{a'} = m, A = a'\right]\mathclose{} \Pr[A = a'] \\ &= \sum_{a'} \operatorname{E}\mathopen{}\left[Y \mid M = m, A = a'\right]\mathclose{} \Pr[A = a']. \end{aligned} \]
The four equalities use, in turn, the law of total expectation, the d-separation of \(M^a\) from \(A\) in the SWIG, the SWIG property just shown, and consistency.
So \(\operatorname{E}\mathopen{}\left[Y^a \mid M^a\right]\mathclose{}\) is the same for every \(a\), yet it generally differs from \(\operatorname{E}\mathopen{}\left[Y \mid M\right]\mathclose{} = \sum_{a'} \operatorname{E}\mathopen{}\left[Y \mid M, A = a'\right]\mathclose{} \Pr[A = a' \mid M]\), because the factual \(M = M^A\), unlike \(M^a\), is associated with \(A\) (Hernán and Robins 2020, 302).