{"id":"2c3d2757-9890-469f-92bd-07d603608e3e","arxiv_id":"1908.07072","paper_version":2,"verdict":"CONDITIONAL","confidence":"HIGH","novelty_score":6.0,"correctness_risk":"low","formal_verification":"none","parameter_count":0,"one_line_summary":"The gfoRmula R package implements the parametric g-formula for estimating effects of general time-varying treatment interventions, with support for survival and end-of-follow-up outcomes, competing events, and flexible user-defined interventions.","lead":"This paper presents gfoRmula, an R package for estimating causal effects of time-varying treatments from longitudinal data using the parametric g-formula. It makes a powerful but demanding method usable by epidemiologists, biostatisticians, and social scientists.","discovery_kind":"new_method","skeptic_critique":{"model":"deepseek-v4-flash","headline":"Competing-risk risk estimator in Eq. (9) indexes the competing-event survival factor one interval too late; if implemented as printed, the package's competing-events feature is biased even under correct models.","rationale":"The reader's weakest assumption was model misspecification bias, a well-known property that the paper itself acknowledges in the Discussion with the Robins and Wasserman reference. That is a genuine limitation but not a defect in the software's internal logic. My stress test found a more specific, actionable issue: Eq. (9) appears inconsistent with Eq. (8) and with the definition of q_k in §3.1 Step 1(c). If the implementation follows the printed formula, then even a perfectly specified parametric g-formula for a survival outcome with competing events would estimate the wrong quantity. That would undercut the paper's central claim for one of its advertised features. The missing known-truth validation that the reader flags is exactly what would expose this bug, but the concern here is sharper than the generic misspecification caveat: it is an internal indexing error, not a fundamental statistical limit. Because the issue is localized and a code correction plus a numerical check would resolve it, the appropriate disposition remains a conditional acceptance rather than outright rejection. I partially agree with the reader: the lack of validation is real, but the specific competing-risk formula inconsistency makes the recommended condition concrete rather than a general call for more testing.","tokens_in":32715,"tokens_out":10527,"duration_ms":114176,"concrete_test":"Simulate longitudinal data with known discrete hazards for the outcome Y and competing event D, e.g. constant hazards p=0.1 and q=0.2, a binary covariate with a known transition model, no censoring, and K+1=6. Specify the true models to gformula_survival with compevent_name set and competing events not treated as censoring, under a static intervention. Compare the package's g-formula risk estimate to the analytic risk R = sum_{k=0}^K E[ p_k(L_k,a) * prod_{j=0}^{k-1} (1-p_j(L_j,a))(1-q_j(L_j,a)) ], or, if using the paper's p_k conditional on D_{k+1}=0, the analytic formula implied by Eq. (8). If the estimate differs by roughly an extra (1-q) factor, the indexing error in Eq. (9) is present in the implementation; if it matches, the code corrects the typo and only the paper's formula needs revision.","verdict_should_be":"UNCHANGED","load_bearing_attack":"The central claim is that gfoRmula recovers the parametric g-formula estimand when the user's models and identification conditions hold. A concrete weak point is the survival/competing-events algorithm. In §3.1 Step 1(c), q_k(l_k, a_k) is defined as Pr[D_{k+1}=1 | L_k=l_k, A_k=a_k, C_{k+1}=D_k=Y_k=0]. The g-formula in (8) therefore multiplies each term by Pr[D_{j+1}=0 | L_j, A_j] = 1 - q_j. But the Monte Carlo estimator in Eq. (9) uses the factor {1 - qhat(a_{j+1,v}, l_{j+1,v})}, i.e. the next interval's competing-event hazard, rather than 1 - qhat(a_{j,v}, l_{j,v}). For j=k this factor refers to q_{K+1}, which is undefined with data through time K. For j<k it would require surviving an extra interval before the outcome event can occur, which is not what Eq. (8) says. If the R code follows the printed equation, the competing-events option does not recover the estimand in (8) even under perfectly specified models. This is an internal inconsistency in the algorithm, not merely the acknowledged model-misspecification caveat; it directly affects an advertised feature of the package.","agreement_with_reader":"partial"},"referee_report":{"model":"deepseek-v4-flash","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.","tokens_in":32978,"tokens_out":7047,"duration_ms":74454,"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":[{"comment":"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":"Section 3.1, Eq. (9)"},{"comment":"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":"Section 3.1, Step 3"},{"comment":"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.","section":"Section 5, Examples"}],"minor_comments":[{"comment":"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":"Section 2.3, Eq. (8)"},{"comment":"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.","section":"Section 3.2"},{"comment":"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":"Sections 6.5 and 6.7"},{"comment":"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.","section":"Section 3.1, Steps 1(c) and 2(d)"}],"recommendation":"major_revision","confidential_remarks":"The main reason for major revision is not the underlying methodology, which is standard and sound, but the internal inconsistency between the printed estimation algorithm and the claimed estimand in the competing-events setting, together with the Monte Carlo normalization issue. I would support acceptance after the authors correct the equations, confirm that the code implements the corrected versions, and add a validating simulation."},"author_rebuttal":null,"desk_editor":{"model":"deepseek-v4-flash","letter":"This is a solid software paper that deserves a serious referee, but it has a concrete algorithmic error in the printed competing-events estimator that has to be fixed.\n\nWhat is genuinely new: the gfoRmula package brings the parametric g-formula to R with flexibility the SAS macro lacks. Users can supply custom covariate fit/predict functions, define custom interventions, handle visit processes, impose deterministic restrictions, and simulate joint interventions. The statistical background in Sections 2 and 3 is mostly correct and well referenced. The worked examples run, and the natural-course estimates track the nonparametric estimates, which is the right sanity check to include.\n\nThe soft spots. First, the indexing error the stress-test caught is real. Equation (8) correctly includes Pr[D_{j+1}=0 | L_j, A_j] in the product for each event time k. Equation (9), however, uses 1 - qhat(a_{j+1,v}, l_{j+1,v}) — the competing-event hazard at the next interval, evaluated at next-interval covariates and treatment. For j = k this refers to a time point that is not even simulated. If the R code mirrors Eq. (9), the competing-events option does not recover the estimand in Eq. (8) even under correctly specified models. Either the equation or the code (or both) must be corrected, and the example rerun, before the paper is publishable. This is not cosmetic; it affects an advertised feature.\n\nSecond, the paper never validates the package against a known-truth simulation or against the SAS macro. The natural-course comparison is reassuring but insufficient — the g-formula estimates under interventions are never checked against an independent estimator. A small simulation study with data generated under the assumptions would substantially strengthen the paper.\n\nMinor points: self-citation is heavy but the cited work is relevant and mostly independent. The Discussion's nod to Robins and Wasserman is honest about the parametric g-formula's vulnerability to misspecification.\n\nThis paper is for epidemiologists and social scientists working with longitudinal observational data who need a practical, flexible implementation of the g-formula. Send it to peer review: the package is useful and the core ideas are sound, but the competing-events equation and the validation gap have to be addressed before it is ready.","headline":"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.","tokens_in":33501,"tokens_out":3565,"would_cite":true,"duration_ms":34218,"reading_group":"yes","serious_thinker":"yes","would_accept_peer_review":true},"rs_alignment":null,"lean_confirmation":null,"pith_extraction":{"msc":["62D20","62-04"],"pacs":[],"model":"deepseek-v4-flash","headline":"The gfoRmula R package implements the parametric g-formula to estimate causal effects of user-specified time-varying treatment interventions from longitudinal data.","keywords":["g-formula","causal inference","time-varying treatment","longitudinal data","parametric g-formula","competing events","R package","dynamic treatment regimes"],"falsifier":"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.","tokens_in":32540,"feed_emoji":"📊","tokens_out":10415,"duration_ms":98339,"temperature":0.7,"pith_summary":"The paper presents gfoRmula, an R package that implements the parametric g-formula for estimating the causal effect of user-specified time-varying treatment interventions on a future outcome from longitudinal data. This matters because standard regression methods for confounding control fail when a time-varying confounder is itself affected by past treatment, a situation common in HIV therapy studies, physical-activity research, and many clinical cohorts. The package estimates the mean or risk under static, dynamic, deterministic, random, and natural-value-dependent treatment rules, and it supports survival outcomes, competing events, and fixed end-of-follow-up outcomes. If the identifying assumptions of exchangeability, positivity, and consistency hold and the user's fitted models are correct, the package recovers the counterfactual mean or risk under each intervention.","feed_headline":"R package gfoRmula estimates causal effects of time-varying treatments","feed_subtitle":"Simulates what would happen under static, dynamic, random, or threshold interventions, even when past treatment affects later confounders.","key_machinery":"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.","core_discovery":"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.","pith_inferences":["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."],"forward_implications":["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."],"supporting_citations":[{"why":"Establishes the g-formula identification result that the package estimates: counterfactual means under regimes equal observed-data sums or integrals under exchangeability, positivity, and consistency.","marker":"Robins 1986"},{"why":"Provides the extended g-formula for interventions depending on the natural value of treatment, which the package implements.","marker":"Robins et al. 2004"},{"why":"Unifies static, dynamic, and natural-value interventions under the generic user rule and supplies identification conditions for the package's estimands.","marker":"Young et al. 2014"},{"why":"Defines how competing events shift the estimand between total and direct effects and supplies the hazard-product forms used in the survival algorithm.","marker":"Young et al. 2019"},{"why":"The SAS GFORMULA macro whose capabilities the R package adapts and extends.","marker":"Logan et al. 2016"},{"why":"Applied the parametric g-formula to dynamic regimens and established the natural-course nonparametric comparison used as a model-checking diagnostic.","marker":"Young et al. 2011"},{"why":"Sets out the causal framework and defines interventions and censoring that the package operationalizes.","marker":"Hernán and Robins 2018"},{"why":"Documents the guaranteed misspecification bias under the null that bounds the package's reliance on correct conditional models.","marker":"Robins and Wasserman 1997"}],"fun_headline_variants":["gfoRmula: causal effects of any time-varying treatment rule","Estimate intervention effects from longitudinal data with gfoRmula","R package gfoRmula: simulate outcomes under static, dynamic, or random interventions","Handle treatment-dependent confounders: gfoRmula computes g-formula estimates","Parametric g-formula made accessible: gfoRmula for R"],"cache_read_input_tokens":3200,"weakest_assumption_plain":"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.","fun_headline_variants_meta":{"raw":{"variants":["gfoRmula: causal effects of any time-varying treatment rule","Estimate intervention effects from longitudinal data with gfoRmula","R package gfoRmula: simulate outcomes under static, dynamic, or random interventions","Handle treatment-dependent confounders: gfoRmula computes g-formula estimates","Parametric g-formula made accessible: gfoRmula for R"]},"model":"deepseek-v4-flash","effort":"low","cost_usd":0.000599,"raw_usage":{"total_tokens":2784,"prompt_tokens":914,"completion_tokens":1870,"prompt_tokens_details":{"cached_tokens":384},"prompt_cache_hit_tokens":384,"prompt_cache_miss_tokens":530,"completion_tokens_details":{"reasoning_tokens":1772}},"tokens_in":530,"tokens_out":1870,"duration_ms":15411,"temperature":1.0,"reasoning_tokens":1772,"cache_read_input_tokens":384,"cache_creation_input_tokens":0},"cache_creation_input_tokens":0},"created_at":"2026-08-14T12:27:02.451700+00:00","model_set":{"reader":"deepseek-v4-flash"},"falsifier":"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.","supporting_citations":[{"cited_title":"A new approach to causal inference in mortality studies with a sustained exposure period: application to the healthy worker survivor effect","cited_arxiv_id":null,"evidence_quote":"Establishes the g-formula identification result that the package estimates: counterfactual means under regimes equal observed-data sums or integrals under exchangeability, positivity, and consistency."},{"cited_title":"Identification, estimation and approximation of risk under interventions that depend on the natural value of treatment using observational data","cited_arxiv_id":null,"evidence_quote":"Unifies static, dynamic, and natural-value interventions under the generic user rule and supplies identification conditions for the package's estimands."},{"cited_title":"GFORMULA SAS MACRO","cited_arxiv_id":null,"evidence_quote":"The SAS GFORMULA macro whose capabilities the R package adapts and extends."},{"cited_title":"Comparative effectiveness of dynamic treatment regimes: an application of the parametric g-formula","cited_arxiv_id":null,"evidence_quote":"Applied the parametric g-formula to dynamic regimens and established the natural-course nonparametric comparison used as a model-checking diagnostic."},{"cited_title":"Estimation of effects of sequential treatments by reparameterizing directed acyclic graphs","cited_arxiv_id":null,"evidence_quote":"Documents the guaranteed misspecification bias under the null that bounds the package's reliance on correct conditional models."}],"review_version":1}