Time-Varying G-Formula in Causal Inference

• Time-Varying G-Formula is a causal inference framework that uses sequential conditional expectations to estimate counterfactual outcomes under dynamic treatment strategies.
• It addresses treatment-confounder feedback in longitudinal studies by integrating simulation techniques like Monte Carlo and nested regression methods.
• Flexible extensions, including deep learning and Bayesian approaches, enhance its application in survival analysis, non-stationary time series, and complex treatment regimens.

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 (Lin et al., 2019, Rein et al., 2024). 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 (Naimi et al., 18 Jun 2026, Chen et al., 2024, Cai et al., 2024, Singh et al., 2021, Wang et al., 28 Apr 2026).

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 (Lin et al., 2019).

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 (Lin et al., 2019). 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 (Chen et al., 2024, Lin et al., 2019).

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 (Cai et al., 2024).

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 (Naimi et al., 18 Jun 2026).

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 (Naimi et al., 18 Jun 2026).

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}$ (Rein et al., 2024). In competing-risk settings, the formula can either treat competing events as censoring or explicitly model their hazards, leading to different causal interpretations (Lin et al., 2019).

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{A}_t$0. The coarsest balancing score is the longitudinal joint score

$\overline{A}_t$1

which suffices for static regimes, while dynamic regimes additionally require tailoring variables (Chen et al., 2024).

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 $\overline{A}_t$2, then observed and counterfactual outcomes under $\overline{A}_t$3 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 (Rein et al., 2024, Naimi et al., 18 Jun 2026).

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 (Cai et al., 2024).

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,

$\overline{A}_t$4

the $\overline{A}_t$5-lag structural direct effect, the $\overline{A}_t$6-lag effect, the $\overline{A}_t$7-step total effect,

$\overline{A}_t$8

the $\overline{A}_t$9-step general effect for arbitrary recent exposure sequences, and cumulative effects over future outcomes (Cai et al., 2024). 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 (Bartlett et al., 2023). 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 (Bartlett et al., 2023).

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 (Rein et al., 2024). 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 (Naimi et al., 18 Jun 2026).

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 (Lin et al., 2019).

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 (Bartlett et al., 2023).

Approach	Models required	Distinguishing feature
Parametric g-formula	Conditional models for covariates and outcome	Monte Carlo simulation under static or dynamic interventions (Lin et al., 2019)
NICE	Joint conditional distributions of treatment, confounders, and outcome	Plug-in two-stage estimator (Rein et al., 2024)
ICE	Outcome regressions on history	Nested expectations; no model for covariate evolution (Naimi et al., 18 Jun 2026)
MI-based g-formula	Sequential imputation models	Missing actual and counterfactual data handled in one framework (Bartlett et al., 2023)
State space g-formula	Latent time-varying coefficient models	Non-stationary time series and feedback (Cai et al., 2024)

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 (Rein et al., 2024).

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 $\overline{Z}_t$0, 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 (Li et al., 2020).

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 (Xiong et al., 2024).

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 (Chen et al., 2024).

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 (Singh et al., 2021). 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 (Wang et al., 28 Apr 2026).

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 (Cai et al., 2024). A representative outcome model is

$\overline{Z}_t$1

where $\overline{Z}_t$2 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 (Cai et al., 2024).

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 (Cai et al., 2024). 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 (Cai et al., 2024).

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 (McGrath et al., 2021). 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 (McGrath et al., 2021).

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 (McGrath et al., 2021). 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 (McGrath et al., 2021). 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 (Lin et al., 2019, Bartlett et al., 2023).

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 (Naimi et al., 18 Jun 2026, Rein et al., 2024).