REVIEW 3 major objections 4 minor 38 references
gfoRmula: An R package for estimating effects of general time-varying treatment interventions via the parametric g-formula
T0 review · 3 major / 4 minor · reviewed 2026-08-14 · deepseek-v4-flash
Pith's one-line read The gfoRmula R package implements the parametric g-formula to estimate causal effects of user-specified time-varying treatment interventions from longitudinal data.
desk verdict A useful, flexible R package for the parametric g-formula, but the competing-events estimator as printed in Eq. (9) has a real indexing error and the lack of simulation validation needs addressing before publication. read the letter →
The pith
A machine-rendered reading of the paper's core claim, the machinery that carries it, and where it could break.
The reading
What carries the argument
The load-bearing object is the parametric g-formula itself: a product of the outcome mean or hazard, the conditional confounder densities at each time, and the user's intervention density, summed over all possible treatment and confounder histories. The package's practical mechanism is Monte Carlo simulation—estimate each conditional piece from pooled generalized linear models or user-supplied fits, then draw forward covariate and treatment histories under the rule $h_{\mathrm{user}}(a_k,a^*_k,l_k)$ and average the estimated hazard or outcome contribution over many simulated histories. Expressing the intervention as a rule or density is what lets one algorithm cover static, dynamic, random, threshold, natural-course, and natural-value-dependent interventions, as well as joint interventions on multiple treatments.
What would settle it
Simulate a longitudinal dataset from a known data-generating process with a time-varying confounder affected by past treatment and a true null treatment effect, then use the package's default pooled logistic and linear models to estimate the risk under 'always treat' versus 'never treat' and compare with the known truth. A systematic deviation beyond Monte Carlo error that disappears only when the fitted models are replaced by the true models would confirm that the package's output depends on model correctness exactly where the paper places that dependence.
Extended reading notes
Core claim
The central claim is that one flexible Monte Carlo implementation of the parametric g-formula can turn an observed longitudinal dataset into estimates of the risk or mean outcome under any intervention expressible as a user-defined treatment rule. The paper formalizes the target as a g-formula sum or integral over treatment and confounder histories, multiplying the observed outcome model, the observed conditional distributions of confounders, and the intervention density $h_{\mathrm{user}}(a_k,a^*_k,l_k)$. The estimation algorithm fits pooled conditional models for covariates, outcome hazards (or end-of-follow-up means), and competing-event hazards, simulates histories forward under the rule, and averages the predicted outcome contributions. Under exchangeability, positivity, and consistency, and with correct model specification, this recovers the causal risk or mean; for rules depending on the natural value of treatment, the same machinery estimates the extended g-formula under stronger conditions.
Load-bearing premise
The package's estimates are only as good as the user's models for future covariates, the outcome, and competing events; even if all causal assumptions hold, a wrong model can bias the estimates, and under a true null effect this bias is guaranteed when past treatment affects later confounders.
Editorial extensions
If this is right
- Under the stated assumptions, a user can estimate the effects of realistic dynamic strategies such as 'start treatment within six months if CD4 count drops below 350' directly from observational data.
- For survival outcomes, treating competing events as censoring targets a controlled direct effect that eliminates competing events, while modeling them separately targets a total effect that may capture treatment effects on the competing event.
- The package can estimate joint interventions on several time-varying treatments at once and can incorporate visit processes, carry-forward of last measured values, and known deterministic restrictions to reduce model dependence.
- Bootstrap confidence intervals, risk or mean ratios and differences, and optional hazard ratios are produced from the same function call, with parallelization available for larger simulations.
Reading between the lines
- An unstated consequence is that the automatic nonparametric natural-course estimate provides a built-in gross misspecification diagnostic; a formal rule for when its divergence from the parametric estimate signals model failure would be a natural addition.
- Because the package allows custom covariate models and intervention functions, it can serve as a simulation platform for benchmarking the parametric g-formula against inverse-probability-weighted and doubly robust estimators under realistic misspecification, a comparison the paper notes has not been thoroughly studied.
- A testable extension would embed cross-validated model selection among the package's covariate types inside the simulation loop, potentially reducing the guaranteed null-state misspecification bias the Discussion describes.
Signed reviews
Editorial analysis
A structured set of objections, weighed in public.
Referee Report
Summary. This paper describes the R package gfoRmula, which implements the parametric g-formula for estimating the mean or risk of an outcome under user-specified time-varying treatment interventions. The package supports static, dynamic, deterministic, random, and natural-value-dependent interventions; survival and end-of-follow-up outcomes; and several options for handling competing events. The paper reviews the identifying assumptions and the g-formula, presents the Monte Carlo estimation algorithm, documents the package interface, and illustrates the software with four worked examples.
Significance. If the implementation is correct, gfoRmula fills a practical gap by providing a flexible, open-source R implementation of the parametric g-formula, extending the capabilities of the existing SAS macro and allowing user-supplied model-fitting functions. The manuscript is clearly written, the identification framework follows the standard Robins g-formula literature, and the Discussion honestly acknowledges the known model-misspecification limitations of the parametric g-formula, including the Robins-Wasserman result. The examples are reproducible and the package is publicly available on CRAN, which are strengths. However, the printed algorithm contains indexing and normalization errors that directly affect the claimed properties of the competing-events feature and the Monte Carlo estimator, and the paper provides no simulation study that would verify the implementation against a known estimand.
major comments (3)
- [Section 3.1, Eq. (9)] Equation (9) multiplies the event-of-interest hazard term by {1 - \hat q(a^user_{j+1,v}, l_{j+1,v})}, but the target g-formula in Eq. (8) and the definition of q_k in Step 1(c) imply that the factor should be 1 - q_j(l_j, a_j) = Pr[D_{j+1}=0 | L_j=l_j, A_j=a_j, C_{j+1}=D_j=Y_j=0]. As printed, the factor for j=k refers to q_{K+1}, which is not defined with data through time K, and for j<k it requires surviving a competing event one interval beyond the interval in which the outcome event is evaluated. If the R code follows the printed equation, the competing-events option is biased even under correctly specified models; if the code is correct, then Eq. (9) misrepresents the package. Please correct the equation and, ideally, add a simulation check that targets the competing-events estimand directly.
- [Section 3.1, Step 3] Step 2 of the algorithm generates s simulated histories, with s possibly different from n, and all worked examples use nsimul = 10000 with n = 2500. Step 3, however, computes the estimate as (1/n) \sum_{v=1}^n, which is inconsistent with the stated Monte Carlo sample size. The estimator should average over the s simulated histories, i.e., (1/s) \sum_{v=1}^s. The same issue appears in the Modified Step 3 for end-of-follow-up outcomes in Section 3.2. As written, the estimator does not use the full simulation sample when s \neq n, and the relationship between the printed algorithm and the package output is unclear.
- [Section 5, Examples] The paper includes no simulation study that verifies the package's estimates against a known g-formula value. The examples compare the parametric natural-course estimate with a nonparametric estimate (e.g., Example 1, where the g-form risk is 0.5048 versus NP risk 0.5056), but this comparison checks internal consistency of the covariate simulation and outcome model under no intervention; it does not validate the competing-events estimand, dynamic or threshold interventions, or the Monte Carlo averaging. Given the indexing issues in Eqs. (8)-(9), a simple simulation with a known data-generating process would provide decisive evidence about whether the implementation recovers the target estimand and should be added.
minor comments (4)
- [Section 2.3, Eq. (8)] Equation (8) contains the typo "C_j=D_j=0=Y_j=0" in the conditioning event; this should read "C_j=D_j=Y_j=0".
- [Section 3.2] The text says the outcome model is fitted "Using only records on line K + 1", but Section 4.2 states that the outcome value on line k = K is used for end-of-follow-up outcomes. Please clarify whether K+1 refers to the interval index or the line number.
- [Sections 6.5 and 6.7] Several code excerpts and the argument summary contain typesetting corruption, such as "/grave.ts1" in Section 6.5 and "/grave.ts1/grave.ts1categorical time" in Section 6.7. These should be repaired so that the code examples are usable.
- [Section 3.1, Steps 1(c) and 2(d)] The argument order in q_k(l_k, a_k) is inconsistent with the order in \hat q(a^user_{j+1,v}, l_{j+1,v}) in Eq. (9); while not substantive, the inconsistency makes the indexing error harder to spot.
Circularity Check
No significant circularity: the package implements a standard plug-in estimator whose identifying assumptions and models are external inputs, not outputs of the package.
full rationale
The central derivation in Section 2.3 begins with exchangeability, positivity, and consistency and defines the g-formula as a sum/integral over observed densities and the intervention density; none of these objects is defined in terms of the package's own estimates. The estimation algorithm in Section 3.1 is a Monte Carlo plug-in estimator: Step 1 fits user-specified models for covariate distributions, the outcome hazard, and the competing-event hazard; Step 2 simulates histories under the user's intervention rule; and Eq. (9) evaluates the empirical analogue of Eq. (8). No estimated parameter is relabeled as a prediction: the output risks are explicit functions of the fitted nuisance models, and the nonparametric natural-course estimate is computed separately by product-limit or Aalen-Johansen methods as an independent benchmark. The examples are illustrations on package datasets, not evidence that the package predicts its own inputs. The paper's self-citations, such as Young et al. (2014) and Young et al. (2019), support background identification formulas that are independently published, while the core g-formula is attributed to Robins (1986). The Discussion candidly acknowledges that model misspecification can bias estimates, and any indexing inconsistency in Eq. (9) would be an implementation-correctness concern rather than circularity. The derivation is therefore self-contained relative to standard causal inference assumptions and user-specified models.
Assumptions & free parameters
assumptions (5)
- domain assumption Exchangeability (no unmeasured confounding): Y^g is independent of (A_k, C_{k+1}) conditional on L_k, A_{k-1}, and C_k = 0 for all k.
- domain assumption Positivity: observed treatment and censoring probabilities are positive whenever the joint density of past treatment and covariates under the strategy is nonzero.
- domain assumption Consistency: if the observed treatment and censoring match the strategy, then the observed outcome equals the counterfactual outcome.
- domain assumption Correct specification of parametric models for conditional covariate distributions and outcome or competing-event hazards.
- domain assumption Extended g-formula identification conditions hold for interventions depending on the natural value of treatment, which are stronger than the three conditions above.
Cite this review
Pith. "Pith review of gfoRmula: An R package for estimating effects of general time-varying treatment interventions via the parametric g-formula." pith.science (2026). https://pith.science/paper/5LNNEAHS
@misc{pith2026190807072,
author = {Pith},
title = {Pith review of: gfoRmula: An R package for estimating effects of general time-varying treatment interventions via the parametric g-formula},
year = {2026},
howpublished = {\url{https://pith.science/paper/5LNNEAHS}},
note = {Machine review of arXiv:1908.07072}
}
read the original abstract
Researchers are often interested in using longitudinal data to estimate the causal effects of hypothetical time-varying treatment interventions on the mean or risk of a future outcome. Standard regression/conditioning methods for confounding control generally fail to recover causal effects when time-varying confounders are themselves affected by past treatment. In such settings, 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. The package gfoRmula implements in R one such estimator: the parametric g-formula. This estimator easily adapts to binary or continuous time-varying treatments as well as contrasts defined by static or dynamic, deterministic or random treatment interventions, as well as interventions that depend on the natural value of treatment. The package accommodates survival outcomes as well as binary or continuous end of follow-up outcomes. For survival outcomes, the package has different options for handling competing events. This paper describes the gfoRmula package, along with motivating background, features, and examples.
Reference graph
Works this paper leans on
-
[1]
An Empirical Transition Matrix for Non-Homogeneous Markov Chains Based on Censored Observations
Aalen O, Johansen S (1978). An Empirical Transition Matrix for Non-Homogeneous Markov Chains Based on Censored Observations . Scandinavian Journal of Statistics, 5, 141--150
work page 1978
-
[2]
Doubly robust estimation in missing data and causal inference models
Bang H, Robins JM (2005). Doubly robust estimation in missing data and causal inference models. Biometrics, 61, 692--972
work page 2005
-
[3]
Cain LE, Robins JM, Lanoy E, Logan R, Costagliola D, Hern\' a n MA (2010). When to start treatment? A systematic approach to the comparison of dynamic regimes using observational data. International Journal of Biostatistics, 6. Article 18
work page 2010
-
[4]
truncreg : Truncated Gaussian Regression Models
Croissant Y, Zeileis A (2018). truncreg : Truncated Gaussian Regression Models . R package version 0.2-5, ://CRAN.R-project.org/package=truncreg
work page 2018
-
[5]
Population Intervention Causal Effects Based on Stochastic Interventions
D\' i az Mu\ n oz I, van der Laan MJ (2012). Population Intervention Causal Effects Based on Stochastic Interventions. Biometrics, 68, 541--549
work page 2012
-
[6]
A proportional hazards model for the subdistribution of a competing risk
Fine J, Gray R (1999). A proportional hazards model for the subdistribution of a competing risk . Journal of the American Statistical Association, 94(446), 496--509
work page 1999
-
[7]
Estimation of the effect of interventions that modify the received treatment
Haneuse S, Rotnitzky A (2013). Estimation of the effect of interventions that modify the received treatment. Statistics in Medicine, 32, 5260--5277
work page 2013
-
[8]
Hern\' a n MA (2010). The hazards of hazard ratios . Epidemiology, 21(1), 13--15
work page 2010
Show all 38 references
-
[9]
Marginal structural models to estimate the causal effect of zidovudine on the survival of H I V -positive men
Hern\' a n MA, Brumback B, Robins JM (2000). Marginal structural models to estimate the causal effect of zidovudine on the survival of H I V -positive men. Epidemiology, 11(5), 561--570
2000
-
[10]
A structural approach to selection bias
Hern\' a n MA, Hern\' a ndez-Di\' a z S, Robins JM (2004). A structural approach to selection bias. Epidemiology, 15, 615--625
2004
-
[11]
Comparison of dynamic treatment regimes via inverse probability weighting
Hern\' a n MA, Lanoy E, Costagliola D, Robins JM (2006). Comparison of dynamic treatment regimes via inverse probability weighting. Basic & Clinical Pharmacology & Toxicology, 98, 237--242
2006
-
[12]
Observation plans in longitudinal studies with time-varying treatments
Hern\' a n MA, McAdams M, McGrath N, Lanoy E, Costagliola D (2009). Observation plans in longitudinal studies with time-varying treatments. Statistical Methods in Medical Research, 18(1), 27--52
2009
-
[13]
Causal Inference
Hern\' a n MA, Robins J (2018). Causal Inference . CRC Press. forthcoming , ://www.hsph.harvard.edu/miguel-hernan/causal-inference-book/
2018
-
[14]
The Statistical Analysis of Failure Time Data
Kalbfleisch JD, Prentice RL (1980). The Statistical Analysis of Failure Time Data . New York: John Wiley
1980
-
[15]
Nonparametric estimation from incomplete observations
Kaplan E, Meier P (1958). Nonparametric estimation from incomplete observations . Journal of the American Statistical Association, 53(282), 457--81
1958
-
[16]
Nonparametric Causal Effects Based on Incremental Propensity Score Interventions
Kennedy E (2019). Nonparametric Causal Effects Based on Incremental Propensity Score Interventions . Journal of the American Statistical Association, 114, 645--656
2019
-
[17]
ltmle : An R Package Implementing Targeted Minimum Loss-Based Estimation for Longitudinal Data
Lendle SD, Schwab J, Petersen ML, van der Laan MJ (2017). ltmle : An R Package Implementing Targeted Minimum Loss-Based Estimation for Longitudinal Data. Journal of Statistical Software, 81(1), 1--21. doi:10.18637/jss.v081.i01
2017 doi
-
[18]
gfoRmula : Parametric G-Formula
Lin V, McGrath S, Zhang Z, Logan RW, Petito LC, Young JG, Hern\' a n MA (2019). gfoRmula : Parametric G-Formula . R package version 0.2.1, ://CRAN.R-project.org/package=gfoRmula
2019
-
[19]
GFORMULA SAS MACRO
Logan R, Young J, Taubman S, Lodi S, Picciotto S, Danaei G, Hern\' a n MA (2016). GFORMULA SAS MACRO . https://www.hsph.harvard.edu/causal/software
2016
-
[20]
Marginal Mean Models for Dynamic Regimes
Murphy SA, van der Laan MJ, Robins JM (2001). Marginal Mean Models for Dynamic Regimes. Journal of the American Statistical Association, 96(456), 1410--23
2001
-
[21]
Dynamic regime marginal structural mean models for estimation of optimal dynamic treatment regimes, P art I : M ain C ontent
Orellana L, Rotnitzky A, Robins JM (2010a). Dynamic regime marginal structural mean models for estimation of optimal dynamic treatment regimes, P art I : M ain C ontent. International Journal of Biostatistics, 6. Article 7
2010
-
[22]
Dynamic regime marginal structural mean models for estimation of optimal dynamic treatment regimes, P art I I : P roofs and A dditional R esults
Orellana L, Rotnitzky A, Robins JM (2010b). Dynamic regime marginal structural mean models for estimation of optimal dynamic treatment regimes, P art I I : P roofs and A dditional R esults. International Journal of Biostatistics, 6. Article 8
2010
-
[23]
Single World Intervention Graphs (SWIGs) : a unification of the counterfactual and graphical approaches to causality
Richardson TS, Robins JM (2013). Single World Intervention Graphs (SWIGs) : a unification of the counterfactual and graphical approaches to causality. Center for Statistics and the Social Sciences, University of Washington Series. Working Paper Number 128 , http://www.csss.was...
2013
-
[24]
A new approach to causal inference in mortality studies with a sustained exposure period: application to the healthy worker survivor effect
Robins JM (1986). A new approach to causal inference in mortality studies with a sustained exposure period: application to the healthy worker survivor effect. Mathematical Modelling, 7, 1393--1512. [Errata (1987) in Computers and Mathematics with Applications 14, 917-921. Add...
1986
-
[25]
Marginal structural models
Robins JM (1998). Marginal structural models. In 1997 P roceedings of the A merican S tatistical A ssociation, S ection on B ayesian S tatistical S cience , pp. 1--10. American Statistical Association
1998
-
[26]
Estimation of the causal effects of time-varying exposures
Robins JM, Hern\' a n MA (2009). Estimation of the causal effects of time-varying exposures. In G Fitzmaurice, M Davidian, G Verbeke, G Molenberghs (eds.), A dvances in L ongitudinal D ata A nalysis , pp. 553--599. Boca Raton, FL: Chapman and Hall/CRC Press
2009
-
[27]
Marginal structural models and causal inference in epidemiology
Robins JM, Hern\' a n MA, Brumback B (2000). Marginal structural models and causal inference in epidemiology. Epidemiology, 11(5), 550--560
2000
-
[28]
Effects of multiple interventions
Robins JM, Hern\' a n MA, Siebert U (2004). Effects of multiple interventions. In M Ezzati, AD Lopez, A Rodgers, CJL Murray (eds.), C omparative Q uantification of H ealth R isks: G lobal and R egional B urden of D isease A ttributable to S elected M ajor R isk F actors . Gene...
2004
-
[29]
Estimation of effects of sequential treatments by reparameterizing directed acyclic graphs
Robins JM, Wasserman L (1997). Estimation of effects of sequential treatments by reparameterizing directed acyclic graphs. In D Geiger, P Shenoy (eds.), P roceedings of the T hirteenth C onference on U ncertainty in A rtificial I ntelligence , pp. 409--420. San Francisco: Morg...
1997
-
[30]
Intervening on risk factors for coronary heart disease: an application of the parametric g-formula
Taubman SL, Robins JM, Mittleman MA, Hern\' a n MA (2009). Intervening on risk factors for coronary heart disease: an application of the parametric g-formula. International Journal of Epidemiology, 38(6), 1599--611
2009
-
[31]
History-Adjusted Marginal Structural Models and Statically-Optimal Dynamic Treatment Regimens
van der Laan MJ, Petersen ML, Joffe MM (2005). History-Adjusted Marginal Structural Models and Statically-Optimal Dynamic Treatment Regimens. International Journal of Biostatistics, 1(1). Article 4
2005
-
[32]
Modern Applied Statistics with S
Venables WN, Ripley BD (2002). Modern Applied Statistics with S. Fourth edition. Springer, New York. ISBN 0-387-95457-0, ://www.stats.ox.ac.uk/pub/MASS4
2002
-
[33]
Inverse probability weighted estimation of risk under representative interventions in observational studies
Young J, Logan R, Robins J, Hern\' a n MA (2018). Inverse probability weighted estimation of risk under representative interventions in observational studies . Journal of the American Statistical Association, 10, In press
2018
-
[34]
A causal framework for classical statistical estimands in failure time settings with competing events
Young J, Stensrud M, Tchetgen Tchetgen E, Hern\' a n M (2019). A causal framework for classical statistical estimands in failure time settings with competing events . arXiv. ArXiv:1806.06136v2 [stat.ME]
2019 arXiv
-
[35]
Comparative effectiveness of dynamic treatment regimes: an application of the parametric g-formula
Young JG, Cain LE, Robins JM, O'Reilly EJ, Hern\' a n MA (2011). Comparative effectiveness of dynamic treatment regimes: an application of the parametric g-formula. Statistics in Biosciences, 3(1), 119--143
2011
-
[36]
Identification, estimation and approximation of risk under interventions that depend on the natural value of treatment using observational data
Young JG, Hern\' a n MA, Robins JM (2014). Identification, estimation and approximation of risk under interventions that depend on the natural value of treatment using observational data. Epidemiologic Methods, 3(1), 1--19
2014
-
[37]
, " * write output.state after.block = add.period write newline
ENTRY address archive author booktitle chapter collaboration doi edition editor eid eprint howpublished institution isbn issn journal key month note number numpages organization pages publisher school series title type url volume year label extra.label sort.label short.list IN...
-
[38]
write newline
" write newline "" before.all 'output.state := FUNCTION n.dashify 't := "" t empty not t #1 #1 substring "-" = t #1 #2 substring "--" = not "--" * t #2 global.max substring 't := t #1 #1 substring "-" = "-" * t #2 global.max substring 't := while if t #1 #1 substring * t #2 gl...
Reviewed August 14, 2026 · model on record in the stance chip above.
Discussion (0). Continue with ORCID to comment.