---
title: Time-Varying G-Formula in Causal Inference
url: https://www.emergentmind.com/topics/time-varying-g-formula
type: topic
---

# Time-Varying G-Formula in Causal Inference

The time-varying g-formula, also called the g-formula or g-computation, is a causal inference method for estimating the effects of time-varying treatment interventions or sustained treatment strategies from longitudinal data, especially when time-varying confounders are themselves affected by past treatment and standard regression or conditioning methods fail [1908.07072] [2410.21531]. It operationalizes the law of iterated expectation to identify mean counterfactual outcomes under static, dynamic, deterministic, or random treatment regimes, and it has been developed for survival analysis, continuous or binary end-of-follow-up outcomes, non-stationary multivariate time series, dose response curves, and platform SMART designs [2606.20078] [2402.02306] [2407.17666] [2111.03950] [2604.25252].

## 1. Longitudinal causal structure and motivation

The central motivation for the time-varying g-formula is treatment-confounder feedback. In longitudinal observational studies, time-varying covariates can both affect subsequent treatment assignment and the outcome and may themselves be influenced by previous treatment. In such settings, standard regression or conditioning methods generally fail to recover causal effects, whereas estimators derived from Robins’s g-formula may recover time-varying treatment effects provided sufficient covariates are measured to control confounding by unmeasured risk factors [1908.07072].

This framework is used to estimate the causal effect of sustained treatment strategies, dynamic treatment regimes, and hypothetical time-varying interventions. It is applicable to binary or continuous time-varying treatments, to static or dynamic interventions, and to deterministic or random interventions, including interventions that depend on the natural value of treatment [1908.07072]. In survival settings, the target may be the causal survival curve or the risk by a fixed horizon under a specified intervention scenario; in other applications, the target may be a binary or continuous end-of-follow-up outcome [2402.02306] [1908.07072].

The same logic extends beyond classical cohort data. In mobile health and wearable-device settings, intensive longitudinal or multivariate time series data allow granular and dynamic evaluation of treatment effect, including persistence and accumulation over time. This has motivated new causal estimands and state space model implementations for individualized, non-stationary N-of-1 series [2407.17666].

## 2. Mathematical representations

A standard longitudinal setup indexes time by \( t = 0, 1, \ldots, T \), with treatment history \( \overline{A}_t \), covariate history \( \overline{Z}_t \) or \( \overline{L}_t \), and an outcome \( Y \). For a regime \( \overline{a}_T \), the non-iterative conditional expectation form of the g-formula can be written as
\[
E(Y^{\overline{a}_T}) =
\int \cdots \int
E(Y \mid \overline{Z}_T, \overline{A}_T = \overline{a}_T)
\prod_{t=0}^T d\mathbb{P}(Z_t \mid \overline{Z}_{t-1}, \overline{A}_{t-1} = \overline{a}_{t-1}).
\]
In this representation, conditional outcome means are averaged over all temporal histories of the time-varying covariates under the intervention of interest [2606.20078].

For two timepoints, the iterative conditional expectation form can be expressed as
\[
E(Y^{\overline{a}_1}) =
E_{Z_0}\Bigl\{
E_{Z_1 \mid Z_0, A_0 = a_0}
\left[
E(Y \mid Z_0, A_0 = a_0, Z_1, A_1 = a_1)
\right]
\Bigr\}.
\]
The NICE form and the ICE form are nonparametrically equivalent: NICE uses weighted averages with explicit modeling of the covariate distributions, whereas ICE uses a sequence of nested conditional expectations and only requires modeling of the outcome, not the covariate distribution [2606.20078].

For survival data, the deterministic g-formula averages conditional failure probabilities over all possible confounder histories, with the products of conditional survival terms and covariate densities. More general formulations additionally integrate over treatment histories under an intervention distribution \( f^{int} \) [2410.21531]. In competing-risk settings, the formula can either treat competing events as censoring or explicitly model their hazards, leading to different causal interpretations [1908.07072].

A further generalization replaces full confounder histories with longitudinal balancing scores. In discrete survival data with binary time-varying treatment, a general class of g-formulas integrates over a lower-dimensional score history \( \overline{b}_t \) rather than the full \( \overline{L}_t \). The coarsest balancing score is the longitudinal joint score
\[
e_t = P(A_t = 1, C_{t+1} = 0 \mid \overline{L}_t, \overline{A}_{t-1}, C_t = Y_t = 0),
\]
which suffices for static regimes, while dynamic regimes additionally require tailoring variables [2402.02306].

## 3. Identification assumptions and dynamic estimands

Identification of the time-varying g-formula rests on consistency, positivity, and exchangeability. In sequential form, consistency states that if an individual’s observed treatment history matches a treatment strategy \( g \), then observed and counterfactual outcomes under \( g \) are identical. Positivity requires positive probability for receiving each treatment level in all strata of the observed covariate and treatment history. Exchangeability requires that, conditional on observed past values, treatment assignment is independent of future counterfactual outcomes [2410.21531] [2606.20078].

In some settings, additional assumptions are introduced to make identification tractable or to define controlled direct effects. For non-stationary N-of-1 time series, these include Markov independence, intermediate variable exchangeability, and periodic stable effects. Markov independence restricts dependence to a relevant window of the past as determined by the DAG and application. Intermediate variable exchangeability is invoked for controlled direct effects. Periodic stable effects allow piecewise constant effects over unknown periods, permitting non-stationarity but with stability inside periods [2407.17666].

The time-varying g-formula is associated with a broad class of causal estimands. In non-stationary N-of-1 series, these include the contemporaneous effect,
\[
\mathrm{CE}_t = Y_t(a_t=1) - Y_t(a_t=0),
\]
the \( q \)-lag structural direct effect, the \( q \)-lag effect, the \( q \)-step total effect,
\[
\mathrm{TE}_t^{(q)} = Y_t(a_{(t-q):t} = \mathbf{1}_{q+1}) - Y_t(a_{(t-q):t} = \mathbf{0}_{q+1}),
\]
the \( q \)-step general effect for arbitrary recent exposure sequences, and cumulative effects over future outcomes [2407.17666]. These estimands directly encode time, lag-effects, and dynamism, rather than only contemporaneous, temporal-average, or population-average effects.

## 4. Estimation workflows and software implementations

The standard implementation of the parametric g-formula uses Monte Carlo simulation. One fits parametric models for each time-varying confounder and for the outcome, simulates synthetic trajectories under the intervention of interest, and averages the simulated outcomes to estimate the counterfactual mean [2301.12026]. For survival outcomes, nonparametric bootstrapping has commonly been used for inference, but this can be computationally costly when both large synthetic samples and many bootstrap replicates are desired [2301.12026].

The NICE estimator is a plug-in, two-stage estimator. First, it estimates the conditional distributions of time-varying confounders, treatments, and outcome given the observed history. Second, it performs Monte Carlo simulation to integrate or sum over possible histories under the hypothetical intervention [2410.21531]. By contrast, the ICE formulation proceeds through nested regression, prediction, and averaging, moving backward through time; in practice it is often described as more tractable because it does not require explicit modeling of the distribution of the time-varying covariates [2606.20078].

The R package `gfoRmula` implements the parametric g-formula for estimating effects of general time-varying treatment interventions. It supports binary or continuous time-varying treatments, static or dynamic interventions, deterministic or random interventions, interventions that depend on the natural value of treatment, and survival, binary end-of-follow-up, or continuous end-of-follow-up outcomes. For survival outcomes, it provides different options for handling competing events [1908.07072].

Multiple imputation has been proposed as an alternative implementation of the g-formula. In this approach, synthetic counterfactual trajectories are recast as an MI problem by augmenting observed data with synthetic rows whose variables are set to missing, then imputing them with standard MI software. A modified version of Rubin’s variance estimator is used because ordinary Rubin pooling does not yield unbiased variance estimates in this setting. This formulation also allows missing actual data and missing counterfactual data to be handled in a unified procedure [2301.12026].

| Approach | Models required | Distinguishing feature |
|---|---|---|
| Parametric g-formula | Conditional models for covariates and outcome | Monte Carlo simulation under static or dynamic interventions [1908.07072] |
| NICE | Joint conditional distributions of treatment, confounders, and outcome | Plug-in two-stage estimator [2410.21531] |
| ICE | Outcome regressions on history | Nested expectations; no model for covariate evolution [2606.20078] |
| MI-based g-formula | Sequential imputation models | Missing actual and counterfactual data handled in one framework [2301.12026] |
| State space g-formula | Latent time-varying coefficient models | Non-stationary time series and feedback [2407.17666] |

## 5. Flexible and Bayesian extensions

A large recent literature replaces classical regression components of the time-varying g-formula with more flexible estimators for the required conditional distributions. One deep learning extension is a unified framework for the NICE estimator that uses multitask recurrent neural networks, specifically LSTMs, to jointly estimate the conditional distributions of the time-varying confounders and treatment, together with a separate LSTM for the outcome. Using simulated data with simple and complex temporal dependencies, this deep learning NICE estimator showed lower bias than the parametric NICE estimator, especially in settings with complex temporal dependence [2410.21531].

G-Net implements G-computation for counterfactual outcome prediction under dynamic treatment regimes using sequential deep learning, primarily RNNs and LSTMs. It models the conditional distributions \( p(L_t \mid \bar{L}_{t-1}, \bar{A}_{t-1}) \), simulates trajectories forward under a chosen strategy, and uses Monte Carlo dropout to approximate draws from a Bayesian posterior predictive distribution. In CVSim experiments, LSTM-based conditional models had notably lower mean squared error than linear models, especially for longer time horizons, and improved calibration relative to classical approaches [2003.10551].

G-Transformer uses a Transformer-based encoder architecture to estimate the conditional distribution of relevant covariates given covariate and treatment history at each time point, then produces Monte Carlo estimates of counterfactual outcomes by simulating forward patient trajectories under treatment strategies of interest. It was evaluated on mechanistic simulated datasets and a real-world sepsis ICU dataset from MIMIC-IV, and it outperformed both classical and state-of-the-art counterfactual prediction models in these settings [2406.05504].

Flexible Bayesian g-formula estimators have also been proposed for causal survival analyses with time-dependent confounding. One formulation uses Bayesian Additive Regression Trees to model the time-evolving generative components and incorporates longitudinal balancing scores, including the longitudinal joint score. In simulation studies, the BART-based methods had substantially lower relative bias and RMSE for survival probability than the parametric method, especially with nonlinear relationships [2402.02306].

Other extensions move beyond parametric or neural-network nuisance models. Sequential kernel embedding uses kernel ridge regression to embed Robins’ g-formula in a reproducing kernel Hilbert space, allowing treatments, mediators, and covariates to be continuous in general spaces while accommodating nonlinear treatment-confounder feedback. The resulting estimators achieve nonasymptotic uniform rates, and the framework is extended to time-varying treatment effects and counterfactual distributions with semiparametric efficiency and weak convergence results [2111.03950]. In adaptive trial methodology, the Bayesian integration G-formula extends g-computation to platform SMARTs by combining information across concurrent and non-concurrent cohorts through adaptive borrowing, using log-distance or commensurate priors to account for possible time effects [2604.25252].

## 6. Non-stationary time series, diagnostics, and limitations

In non-stationary multivariate time series from N-of-1 mobile device data, the g-formula has been embedded in a state space model framework to accommodate time-varying behavior of treatment effects while accounting for exposure-outcome and outcome-covariate feedback [2407.17666]. A representative outcome model is
\[
Y_t = \beta_{0,t} + \rho_t Y_{t-1} + \beta_{1,t} A_t + \beta_{2,t} A_{t-1} + \beta_{c,t} C_{t-1} + v_t,
\]
where \( \beta_{1,t} \) is a possibly time-varying contemporaneous effect. Kalman filter and smoother algorithms are used for estimation and process tracking, and counterfactual potential outcomes are imputed by simulation from the fitted state space models [2407.17666].

This line of work has introduced diagnostic and interpretive devices tailored to longitudinal interventions. An impulse impact plot summarizes how the effect of a single exposure evolves over multiple lags. A step response plot summarizes sustained exposures. A positivity validation plot charts which exposure sequences are actually observed in the data and therefore where causal estimands are empirically supportable [2407.17666]. In an application to a multi-year smartphone observational study of bipolar patients, the method estimated the dynamic effect of phone-based communication on mood and revealed substantial heterogeneity in treatment effects over time and across individuals [2407.17666].

A principal limitation of the parametric g-formula is the g-null paradox. Under the conditions that motivate g-formula use—identifiability, time-varying confounders affected by past treatment, and the sharp null hypothesis—non-saturated parametric component models can be guaranteed to be misspecified, implying biased estimates even when the treatment effect is null [2103.03857]. Analytic examples and simulations show that bias, poor standard errors, and reduced confidence interval coverage worsen as models become less flexible and as the number of time points increases [2103.03857].

This has several practical implications. Avoiding overly parsimonious models is important; sensitivity analyses using different model specifications and different factorizations of the g-formula are recommended; and natural course validation can be informative, although it is not a guarantee against bias under interventions [2103.03857]. A common misconception is that replacing parametric regressions with machine learning automatically resolves the paradox. The available evidence does not support that view: simply switching to machine learning for the conditional models in the parametric g-formula does not resolve the paradox and can even hurt performance [2103.03857]. Positivity violations, unmeasured confounding, model misspecification, missing data, censoring, and competing events therefore remain central to interpretation and implementation across all variants of the time-varying g-formula [1908.07072] [2301.12026].

The time-varying g-formula is thus best understood not as a single estimator but as a family of identification formulas and estimation strategies for longitudinal causal inference. Its unifying function is to express counterfactual means under time-varying regimes through observable conditional distributions and nested expectations, while its practical performance depends on the strength of the identifying assumptions, the adequacy of the nuisance models, and the empirical support available in the observed histories [2606.20078] [2410.21531].

Source: https://www.emergentmind.com/topics/time-varying-g-formula