Chapter 21: G-Methods for Time-Varying Treatments
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.
This chapter is based on Hernán and Robins (2020, chap. 21, pp. 277-303).
Roadmap: each method was introduced earlier for a time-fixed treatment: IP weighting of marginal structural models (Chapter 12), the g-formula (Chapter 13), and g-estimation of structural nested models (Chapter 14) (the book’s introduction to Chapter 21 says Chapter 15 (Hernán and Robins 2020, 277), but g-estimation is the subject of Chapter 14). All three g-methods rest on the same identifying assumptions but model different parts of the data distribution: the g-formula models the outcome and covariates, IP weighting models treatment, and g-estimation models treatment plus a structural model for the effect of each “blip” of treatment. Section 21.5 adds censoring, and Section 21.6 relates the g-formula to other identifying formulas.
1 21.1 The G-Formula for Time-Varying Treatments (pp. 277-281)
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.
Robins (1986, 1987) introduced the g-formula for time-varying treatments, as cited in Hernán and Robins (2020, 278). Every factor is conditional on prior treatment and covariate history; this conditioning is unnecessary in the time-fixed case, where treatment and confounders are measured at a single time.
Positivity: the g-formula is computable only if, for every \(l_1\) with \(f(l_1 \mid a_0) \neq 0\), some individuals have \((A_0 = a_0, A_1 = a_1, L_1 = l_1)\) (the definition of positivity in Technical Point 19.2).
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).
1.1 The G-Formula as a Simulation
Two cautions (Hernán and Robins 2020, 279):
- The value of the g-formula depends on what is in \(L\). If we wrongly omit \(L_1\) (believing it is not a confounder), the g-formula reduces to \(\operatorname{E}\mathopen{}\left[Y \mid A_0 = a_0, A_1 = a_1\right]\mathclose{}\), which has no causal interpretation when \(L_1\) affects \(A_1\) (Figure 20.8).
- The g-formula can be causal even when its components are not. Under Figure 20.9, where only static sequential exchangeability holds, the g-formula including \(L_1\) still identifies \(\operatorname{E}\mathopen{}\left[Y^{\bar{a}}\right]\mathclose{}\), yet neither \(\Pr[L_1 = l_1 \mid A_0 = a_0]\) nor \(\operatorname{E}\mathopen{}\left[Y \mid A_0 = a_0, A_1 = a_1, L_1 = l_1\right]\mathclose{}\) equals its counterfactual analogue. In a sequentially randomized trial (Figures 20.1 and 20.2) the components are causal as well: \(\Pr[L_1 = l_1 \mid A_0 = a_0] = \Pr[L_1^{a_0} = l_1]\) and \(\operatorname{E}\mathopen{}\left[Y \mid A_0 = a_0, A_1 = a_1, L_1 = l_1\right]\mathclose{} = \operatorname{E}\mathopen{}\left[Y^{a_0, a_1} \mid L_1^{a_0} = l_1\right]\mathclose{}\), so the g-formula is \(\sum_{l_1} \operatorname{E}\mathopen{}\left[Y^{a_0, a_1} \mid L_1^{a_0} = l_1\right]\mathclose{} \Pr[L_1^{a_0} = l_1] = \operatorname{E}\mathopen{}\left[Y^{a_0, a_1}\right]\mathclose{}\).
1.2 Estimation: the Plug-In (Parametric) G-Formula
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.
Software: the gfoRmula R package (Lin et al. 2019) is on CRAN, and a GFORMULA SAS macro is available on GitHub (Hernán and Robins 2020, 281).
2 21.2 IP Weighting for Time-Varying Treatments (pp. 282-285)
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.
As in Technical Point 12.2, the pseudo-population mean equals \(\operatorname{E}\mathopen{}\left[W^{\bar{A}} Y I(A_0 = a_0, A_1 = a_1)\right]\mathclose{}\) (nonstabilized) or the ratio \(\operatorname{E}\mathopen{}\left[SW^{\bar{A}} Y I(\cdot)\right]\mathclose{} / \operatorname{E}\mathopen{}\left[SW^{\bar{A}} I(\cdot)\right]\mathclose{}\) (stabilized), whether or not sequential exchangeability holds (Hernán and Robins 2020, 282).
In the pseudo-population, the probability of treatment at each \(k\) is the constant \(1/2\) (nonstabilized weights) or depends at most on past treatment (stabilized weights), so sequential unconditional exchangeability holds there, and \(\operatorname{E}\mathopen{}\left[Y^{\bar{a}}\right]\mathclose{} - \operatorname{E}\mathopen{}\left[Y^{\bar{a}'}\right]\mathclose{} = \operatorname{E}_{ps}\mathopen{}\left[Y \mid \bar{A} = \bar{a}\right]\mathclose{} - \operatorname{E}_{ps}\mathopen{}\left[Y \mid \bar{A} = \bar{a}'\right]\mathclose{}\). In a true sequentially randomized trial the treatment probabilities are known by design, so nonstabilized-weight estimates are unbiased; in observational studies they must be estimated, for example by a pooled logistic model for \(\Pr[A_k = 1 \mid \bar{A}_{k-1}, \bar{L}_k]\) that includes functions of time \(k\).
A misspecified denominator model biases the estimate; a misspecified numerator model for \(f(A_k \mid \bar{A}_{k-1})\) does not.
Comparing the two parametric estimators: if the parametric g-formula and the parametric IP weighted estimates differ by more than sampling variability (quantified by bootstrapping the difference), then at least one set of models is misspecified, whether or not the identifiability assumptions hold. Similar estimates do not prove the models are correct, because both may be biased in the same direction. The book recommends estimating \(\operatorname{E}\mathopen{}\left[Y^{\bar{a}}\right]\mathclose{}\) both ways and, if they differ substantially by a prespecified criterion, revisiting the models (Hernán and Robins 2020, 284).
2.1 Marginal Structural Models
Variance: use the nonparametric bootstrap, the analytic variance, or a conservative 95% interval from the robust variance estimator. For a non-saturated model, intervals are typically narrower with \(SW^{\bar{A}}\) than with \(W^{\bar{A}}\), so \(SW^{\bar{A}}\) is preferred (Hernán and Robins 2020, 285).
Checking the MSM: fitting
\[\operatorname{E}\mathopen{}\left[Y \mid \bar{A}\right]\mathclose{} = \theta_0 + \theta_1 \mathrm{cum}(\bar{A}) + \theta_2 \mathrm{cum}_{-5}(\bar{A}) + \theta_3 \mathrm{cum}(\bar{A})^2\]
with the same weights, where \(\mathrm{cum}_{-5}\) is cumulative treatment in the final 5 months, gives a 2-degree-of-freedom Wald test of \(\theta_2 = \theta_3 = 0\), a test that the original MSM is correctly specified. In practice one may use other summaries of treatment history and flexible functions such as cubic splines.
With treatment-confounder feedback, \(V\) may include only baseline variables; if \(V\) included any \(L_k\) with \(k > 0\), the weighted-fit coefficients of \(\mathrm{cum}(\bar{A})\) and \(\mathrm{cum}(\bar{A}) V\) could be non-null even though \(\beta_1 = \beta_3 = 0\) under the null (Hernán and Robins 2020, 285).
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\).
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).
3 21.3 A Doubly Robust Estimator for Time-Varying Treatments (pp. 286-288)
3.1 Review: a Doubly Robust Plug-In Estimator for a Time-Fixed Treatment
3.2 Extension to Time-Varying Treatments
For a binary \(Y\), the time-\(m\) model can be logistic, \(\operatorname{expit}(\gamma_m X_m + \varsigma_m \hat{W}^{\bar{1}_m})\) with \(X_m\) a function of \(\bar{L}_m\): although \(\hat{B}_{m+1}\) is not 0 or 1, it lies in \([0, 1]\). For a continuous \(Y\), use a linear model (Hernán and Robins 2020, 288).
This estimator is due to Bang and Robins (2005), building on Robins (2000). It is a targeted minimum loss-based estimator (TMLE) in the terminology of van der Laan and coauthors (Hernán and Robins 2020, 287).
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.
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)\).
- Fit the logistic model for \(\hat{b}_1(L_0, L_1)\) by maximum likelihood among individuals with \(A_0 = A_1 = 1\), adding the single covariate \(A_0 A_1 / (\hat{\pi}_0 \hat{\pi}_1)\) with its own coefficient. The score equation for that coefficient is \(\sum_i \frac{A_{0i} A_{1i}}{\hat{\pi}_{0i} \hat{\pi}_{1i}} \mathopen{}\left(Y_i - \hat{b}_1(L_{0i}, L_{1i})\right)\mathclose{} = 0\), which is the first correction term (individuals outside the fitting sample contribute 0 because \(A_0 A_1 = 0\)).
- Fit the logistic model for \(\hat{b}_0(L_0)\) among individuals with \(A_0 = 1\), with response \(\hat{b}_1(L_0, L_1)\) and the extra covariate \(A_0 / \hat{\pi}_0\); the same argument makes the second correction term sum to 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).
Fact-check note: the book’s first display of \(\hat{U}_{TR}\) (Hernán and Robins 2020, 288) writes the \(\hat{b}_1\) term as \(-\{A_0 A_1/(\hat{\pi}_0 \hat{\pi}_1) - 1\}\hat{b}_1\), which disagrees with the rearranged form the same Technical Point gives (the one shown here). Expanding the rearranged form term by term gives
\[\hat{b}_0 + \frac{A_0 A_1}{\hat{\pi}_0 \hat{\pi}_1}\mathopen{}\left(Y - \hat{b}_1\right)\mathclose{} + \frac{A_0}{\hat{\pi}_0}\mathopen{}\left(\hat{b}_1 - \hat{b}_0\right)\mathclose{} = \frac{A_0 A_1 Y}{\hat{\pi}_0 \hat{\pi}_1} - \mathopen{}\left\{\frac{A_0 A_1}{\hat{\pi}_0 \hat{\pi}_1} - \frac{A_0}{\hat{\pi}_0}\right\}\mathclose{}\hat{b}_1 - \mathopen{}\left\{\frac{A_0}{\hat{\pi}_0} - 1\right\}\mathclose{}\hat{b}_0,\]
so the “\(-1\)” inside the braces on \(\hat{b}_1\) should be \(-A_0/\hat{\pi}_0\). The two versions differ by \((1 - A_0/\hat{\pi}_0)\,\hat{b}_1\), which does not average to 0 in general, so only the corrected form supports the three robustness arguments. The book’s third rearrangement (used for the Molina et al. argument) agrees with the corrected form.
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).
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.
4 21.4 G-Estimation for Time-Varying Treatments (pp. 289-296)
Both equations are saturated: \(\beta_0\) covers the single possible history at time 0, and the four \(\beta_{1j}\) cover the four histories \((A_0, L_1)\) at time 1. The effect of \(a_1\) is
- \(\beta_{11}\) when \(A_0 = 0, L_1^{a_0} = 0\);
- \(\beta_{11} + \beta_{12}\) when \(A_0 = 0, L_1^{a_0} = 1\);
- \(\beta_{11} + \beta_{13}\) when \(A_0 = 1, L_1^{a_0} = 0\);
- \(\beta_{11} + \beta_{12} + \beta_{13} + \beta_{14}\) when \(A_0 = 1, L_1^{a_0} = 1\).
The two equations nested within each other are why the model is called nested (Hernán and Robins 2020, 291).
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.
The proof is in Robins (1994). To fit an unsaturated SNMM by g-estimation, positivity is not required (Hernán and Robins 2020, 292).
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.
4.1 General Structural Nested Mean 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.
4.2 G-Estimation Algorithm (One Parameter)
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).
Fact-check note: the book’s display (Hernán and Robins 2020, 296) multiplies by \(A_{i,k}\) outside the sum and defines \(S_{i,k}\) as a sum of \(R_{i,j}\) without the treatment indicators \(A_{i,j}\). Since \(H_k(\beta) = Y - \sum_{j \ge k} A_j \gamma_j\) (Definition 11), substituting the linear blip gives \({\beta}^{\top} \sum_{j \ge k} A_j R_j\), which is the form used above; the two agree when \(K = 0\) but not in general.
4.3 From \(\hat \beta\) to Counterfactual Means
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).
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).
5 21.5 Censoring Is a Time-Varying Treatment (pp. 297-299)
Part II treated censoring \(C\) as time-fixed.
The superscript \(\bar{c} = \bar{0}\) makes explicit the contrast many investigators have in mind when they speak of “the effect of treatment” (Hernán and Robins 2020, 298).
Conditioning on being uncensored induces selection bias under the null when \(C\) is a collider on a path between \(A\) and \(Y\), or a descendant of one (Chapter 8).
5.1 IP Weighting with Censoring
After fitting the censoring models, compute the sample mean of the estimated \(SW^{\bar{C}}\). Correctly specified models give a mean close to 1 (Hernán and Robins 2020, 298), so for \(SW^{\bar{C}}\) a mean well away from 1 suggests that the numerator or the denominator model is misspecified. The same check applies to \(SW^{\bar{A}}\) with a narrower reading: given a correct denominator \(f_k\), taking iterated expectations over \(a_k\) at each time \(k\) gives \(\operatorname{E}\mathopen{}\left[\prod_k n_k / f_k\right]\mathclose{} = 1\) for any proper conditional numerator \(n_k\), so for \(SW^{\bar{A}}\) a mean away from 1 points to the denominator (treatment) model, or to extreme weights from near-violations of positivity, not to the numerator.
Fact-check note: the book states that “regardless of the type of IP weights used”, the pseudo-population has no arrows from \(L_k\) and \(A_k\) into future \(C_m\) (Hernán and Robins 2020, 299). For stabilized weights this holds for \(L_k\) only: their numerator \(\Pr[C_k = 0 \mid C_{k-1} = 0, \bar{A}_{k-1}]\) keeps the dependence of censoring on past treatment, which is harmless because the analysis conditions on \(\bar{A}\).
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).
6 21.6 The Big G-Formula (pp. 300-303)
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.
6.1 Reducing the Big G-Formula to Observed Data
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.
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).
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).
7 Summary
- In Table 1, the treatment has no effect; the g-formula, IP weighting, and g-estimation all give \(\operatorname{E}\mathopen{}\left[Y^{1,1}\right]\mathclose{} - \operatorname{E}\mathopen{}\left[Y^{0,0}\right]\mathclose{} = 60 - 60 = 0\), where traditional methods failed.
- G-formula: standardize to the distribution of each \(L_k\) given past treatment and covariates; estimate with parametric models (the parametric g-formula), which is subject to the g-null paradox.
- IP weighting: \(W^{\bar{A}}\) or \(SW^{\bar{A}}\), a product over time, fit marginal structural models by weighted regression; \(SW^{\bar{A}}\) preferred for unsaturated models.
- Doubly robust: sequential outcome regressions with the time-varying IP weight as a covariate, consistent if either set of models is correct, and in fact \(K + 2\) robust.
- G-estimation: one SNMM equation per time; find \(\psi^\dagger\) making \(H_k(\psi^\dagger)\) independent of \(A_k\) given the past.
- Censoring is handled like treatment: the joint strategy \((\bar{a}, \bar{c} = \bar{0})\) and weights \(W^{\bar{C}}\) or \(SW^{\bar{C}}\).
- The big g-formula uses all variables, measured and unmeasured; other identifying formulas, such as the front door formula, are reductions of it to observed data.