causalMed and gfoRmula: A Comparative Overview

Overview

gfoRmula (McGrath et al. 2020) is the established CRAN implementation of the parametric g-formula in R, and causalMed uses it as the benchmark for its total-effect function, gformula(). This vignette sets out what the two packages share, what gfoRmula offers that causalMed does not, what causalMed adds, and how their estimates compare on the same data.

library(causalMed)
library(data.table)

Shared foundation

Both packages implement the same algorithm, that of the GFORMULA SAS macro (Logan, Young, Taubman, Hernán and colleagues), of which gfoRmula is the R port: fit parametric models to the long-format data, resample baseline records into a Monte Carlo cohort, simulate the time-varying variables forward one time step at a time under each intervention, compute survival risk as \(1 - \prod_t (1 - h_t)\) with censoring removed under intervention, and bootstrap subjects for confidence intervals.

The GFORMULA macro and gfoRmula estimate total effects only. The direct and indirect effects of causalMed go beyond both; its interventional mediator permutation follows a separate SAS macro, mGFORMULA (Lin et al. 2017).

What gfoRmula offers that causalMed does not

gfoRmula is the more mature and complete package for total-effect estimation, with a dedicated publication (McGrath et al. 2020), documentation and vignettes. causalMed is at an earlier stage of development.

More covariate types.

Covariate type	gfoRmula	causalMed
`"binary"`, `"normal"`, `"categorical"`	✓	✓
`"bounded normal"` (model fitted on a rescaled response)	✓	—
`"zero-inflated normal"`	✓	—
`"truncated normal"`	✓	—
`"absorbing"`	✓	—
`"categorical time"`	✓	—

Both packages clip simulated normal values to the observed range by default: spec_model(truncate = TRUE) in causalMed and sim_trunc = TRUE in gfoRmula, documented there as “whether to truncate simulated covariates to their range in the observed data set”. gfoRmula’s "bounded normal" is a different thing, a covariate type fitted on a response rescaled to [0, 1]. For the other types, causalMed offers only var_type = "custom" with user-written fitting and simulation functions, which leaves the choice and validation of the distribution to the user.

Competing events. gfoRmula models competing events through compevent_name; causalMed has no competing-event interface.

Intervention functions. gfoRmula interventions are R functions that receive the simulated data and the time index, and can be limited to chosen time points (int_times). causalMed uses dyn_int(), an expression evaluated in the simulated data at each step; the current time_var value is in scope (e.g. dyn_int(as.numeric(A > 0 & time >= 2))), and history-dependent rules use lagged columns maintained with the recode hooks.

What causalMed adds: mediation

gfoRmula does not do mediation analysis. mediation() decomposes the effect of a time-varying exposure into direct and indirect effects, as interventional effects (mediation_type = "I"; Lin et al. 2017), with one or more mediators, or natural effects ("N"; Zheng & van der Laan 2017), estimated by g-computation or targeted minimum loss-based estimation. See vignette("causalMed-01-overview") for a quick start and vignette("causalMed-02-mediation") for the details.

Interfaces

Aspect	causalMed	gfoRmula
Models	One `spec_model()` per variable (formula, `var_type`, `mod_type`), in a list	`covnames`, `covtypes`, `covparams`, `ymodel`, `outcome_type`
Histories	Expressions in `init_recode`, `in_recode`, `out_recode` via `recodes()`	History functions such as `lagged` and `cumavg`, via `histories` and `histvars`
Interventions	Named list of `NULL` (natural course), static values or `dyn_int()` rules	`intervention<i>.<variable>` arguments holding an intervention function and its values

Side-by-side estimates

For total effects the two packages implement the same estimator, so with the same models and truncation settings their estimates should agree up to Monte Carlo error. We compare them on a binary end-of-follow-up outcome and a survival outcome.

Data preparation

data("nonsurvivaldata", package = "causalMed")
dat <- as.data.table(nonsurvivaldata)
dat[, time := as.integer(time)]
setorder(dat, id, time)

# Lag recodes (applied inside the Monte Carlo loop)
init_rc <- recodes(lag1_A = 0, lag1_L1 = 0, lag1_L2 = 0)
in_rc   <- recodes(lag1_A = A, lag1_L1 = L1, lag1_L2 = L2)

These models use a simplified L → A ordering within each period, not the A → L process of ?nonsurvivaldata. Both packages receive the same specification, so the comparison is unaffected; in applied work the order should match the assumed data-generating process.

Example 1: Binary end-of-follow-up outcome

causalMed

models_cm <- list(
  spec_model(L1 ~ lag1_A + lag1_L1 + V + time,
             var_type = "normal", mod_type = "covariate"),
  spec_model(L2 ~ lag1_A + lag1_L2 + V + time,
             var_type = "binary", mod_type = "covariate"),
  spec_model(A  ~ lag1_A + L1 + L2 + V + time,
             var_type = "binary", mod_type = "exposure"),
  spec_model(Y_bin ~ A + L1 + L2,
             var_type = "binary", mod_type = "outcome")
)

fit_cm <- causalMed::gformula(
  data         = dat,
  id_var       = "id",
  time_var     = "time",
  base_vars    = "V",
  exposure     = "A",
  models       = models_cm,
  intervention = list(natural = NULL, always = 1, never = 0),
  ref_int      = "natural",
  init_recode  = init_rc,
  in_recode    = in_rc,
  mc_sample    = 10000,
  R            = 1,
  quiet        = TRUE,
  seed         = 20250915
)

fit_cm$effect_size
#>    Intervention       Est
#>          <fctr>     <num>
#> 1:      natural 0.2345816
#> 2:       always 0.2550366
#> 3:        never 0.1029334

gfoRmula

suppressPackageStartupMessages(requireNamespace("gfoRmula"))

fit_gf <- gfoRmula::gformula(
  obs_data     = copy(dat),
  id           = "id",
  time_name    = "time",
  time_points  = length(unique(dat$time)),
  covnames     = c("L1", "L2", "A"),
  covtypes     = c("normal", "binary", "binary"),
  basecovs     = "V",
  covparams    = list(covmodels = c(
    L1 ~ lag1_A + lag1_L1 + V + time,
    L2 ~ lag1_A + lag1_L2 + V + time,
    A  ~ lag1_A + L1 + L2 + V + time
  )),
  ymodel       = Y_bin ~ A + L1 + L2,
  outcome_name = "Y_bin",
  outcome_type = "binary_eof",
  histories    = c(gfoRmula::lagged),
  histvars     = list(c("A", "L1", "L2")),
  intervention1.A = list(gfoRmula::static, rep(1, length(unique(dat$time)))),
  intervention2.A = list(gfoRmula::static, rep(0, length(unique(dat$time)))),
  int_descript = c("Always treat", "Never treat"),
  nsimul       = 10000,
  seed         = 20250915
)

fit_gf$result[, c("k", "g-form mean")]
#>        k g-form mean
#>    <num>       <num>
#> 1:     4   0.2365039
#> 2:     4   0.2562084
#> 3:     4   0.1047391

Side-by-side at final time point

gf_res <- fit_gf$result[k == max(fit_gf$result$k),
                         c("Interv.", "g-form mean")]
gf_res[, Intervention := c("natural", "always", "never")]
setnames(gf_res, "g-form mean", "gfoRmula")
gf_res[, `:=`(gfoRmula = round(gfoRmula, 4), Interv. = NULL)]

cm_res <- fit_cm$effect_size[, .(Intervention, causalMed = round(Est, 4))]
merge(cm_res, gf_res, by = "Intervention")
#> Key: <Intervention>
#>    Intervention causalMed gfoRmula
#>          <char>     <num>    <num>
#> 1:       always    0.2550   0.2562
#> 2:      natural    0.2346   0.2365
#> 3:        never    0.1029   0.1047

Example 2: Survival (time-to-event) outcome

causalMed

data("survivaldata", package = "causalMed")
dat_s <- as.data.table(survivaldata)
setorder(dat_s, id, time)

models_surv <- list(
  spec_model(L ~ V + lag1_L + time,
             var_type = "normal", mod_type = "covariate"),
  spec_model(A ~ V + lag1_A + lag1_L + L + time,
             var_type = "binary", mod_type = "exposure"),
  spec_model(Y ~ lag1_A + A + L + lag1_L + time,
             var_type = "binary", mod_type = "survival")
)

fit_cm_s <- causalMed::gformula(
  data         = dat_s,
  id_var       = "id",
  time_var     = "time",
  base_vars    = "V",
  exposure     = "A",
  models       = models_surv,
  intervention = list(natural = NULL, never = 0, always = 1),
  ref_int      = "natural",
  init_recode  = recodes(lag1_L = 0, lag1_A = 0),
  in_recode    = recodes(lag1_L = L, lag1_A = A),
  mc_sample    = 10000,
  R            = 1,
  quiet        = TRUE,
  seed         = 20250915
)

fit_cm_s$effect_size
#>    Intervention       Est
#>          <fctr>     <num>
#> 1:      natural 0.6457796
#> 2:        never 0.4036755
#> 3:       always 0.7849632

gfoRmula

T_max <- length(unique(dat_s$time))

fit_gf_s <- gfoRmula::gformula(
  obs_data     = copy(dat_s),
  id           = "id",
  time_name    = "time",
  time_points  = T_max,
  covnames     = c("L", "A"),              # same order as the causalMed list
  covtypes     = c("normal", "binary"),
  basecovs     = "V",
  covparams    = list(covmodels = c(
    L ~ V + lag1_L + time,
    A ~ V + lag1_A + lag1_L + L + time
  )),
  ymodel       = Y ~ lag1_A + A + L + lag1_L + time,
  outcome_name = "Y",
  outcome_type = "survival",
  histories    = c(gfoRmula::lagged),
  histvars     = list(c("A", "L")),
  intervention1.A = list(gfoRmula::static, rep(0, T_max)),
  intervention2.A = list(gfoRmula::static, rep(1, T_max)),
  int_descript = c("Never treat", "Always treat"),
  nsimul       = 10000,
  seed         = 20250915,
  nsamples     = 0
)

fit_gf_s$result[k == max(fit_gf_s$result$k), c("Interv.", "g-form risk")]
#>    Interv. g-form risk
#>      <num>       <num>
#> 1:       0   0.6445117
#> 2:       1   0.4036506
#> 3:       2   0.7849414

Side-by-side at final time point

gf_s_res <- fit_gf_s$result[k == max(fit_gf_s$result$k),
                              c("Interv.", "g-form risk")]
gf_s_res[, Intervention := c("natural", "never", "always")]
setnames(gf_s_res, "g-form risk", "gfoRmula")
gf_s_res[, `:=`(gfoRmula = round(gfoRmula, 4), Interv. = NULL)]

cm_s_res <- fit_cm_s$effect_size[, .(Intervention,
                                      causalMed = round(Est, 4))]
merge(cm_s_res, gf_s_res, by = "Intervention")
#> Key: <Intervention>
#>    Intervention causalMed gfoRmula
#>          <char>     <num>    <num>
#> 1:       always    0.7850   0.7849
#> 2:      natural    0.6458   0.6445
#> 3:        never    0.4037   0.4037

Example 3: Dynamic intervention

The risk under “treat only if L1 > 0”. The syntax differs; the estimand is the same.

causalMed: dyn_int() holds the rule as an expression, evaluated in the simulated data at each time step:

fit_dyn <- causalMed::gformula(
  data         = dat,
  id_var       = "id",
  time_var     = "time",
  base_vars    = "V",
  exposure     = "A",
  models       = models_cm,
  intervention = list(natural = NULL, treat_if_L1_pos = dyn_int(as.numeric(L1 > 0))),
  ref_int      = "natural",
  init_recode  = init_rc,
  in_recode    = in_rc,
  mc_sample    = 5000,
  R            = 1,
  quiet        = TRUE,
  seed         = 20250915
)

fit_dyn$effect_size
#>       Intervention       Est
#>             <fctr>     <num>
#> 1:         natural 0.2348336
#> 2: treat_if_L1_pos 0.2303614
fit_dyn$estimate
#>                 Intervention  Risk_type     Estimate
#>                       <char>     <char>        <num>
#> 1: treat_if_L1_pos - natural Difference -0.004472208
#> 2: treat_if_L1_pos / natural      Ratio  0.980955847

gfoRmula: a user-written intervention function:

treat_if_L1_pos <- function(newdf, pool, intvar, intvals, time_name, t) {
  newdf[, (intvar) := as.integer(L1 > 0)]
}

gfoRmula::gformula(
  ...,
  intervention1.A = list(treat_if_L1_pos),
  int_descript    = "Treat if L1 > 0"
)

Both evaluate the rule in the simulated data at each time step.

Larger Monte Carlo sample

Repeating Example 1 with 50,000 rather than 10,000 simulated subjects shrinks the Monte Carlo part of the difference. Both packages fit the same models to the same data, so what remains is Monte Carlo variability.

fit_cv_cm <- causalMed::gformula(
  data         = dat, id_var = "id", time_var = "time", base_vars = "V",
  exposure     = "A", models = models_cm,
  intervention = list(natural = NULL, always = 1),
  ref_int      = "natural",
  init_recode  = init_rc, in_recode = in_rc,
  mc_sample    = 50000, R = 1, quiet = TRUE, seed = 20250915
)

fit_cv_gf <- gfoRmula::gformula(
  obs_data = copy(dat), id = "id", time_name = "time",
  time_points = length(unique(dat$time)),
  covnames = c("L1","L2","A"), covtypes = c("normal","binary","binary"),
  basecovs = "V",
  covparams = list(covmodels = c(
    L1 ~ lag1_A + lag1_L1 + V + time,
    L2 ~ lag1_A + lag1_L2 + V + time,
    A  ~ lag1_A + L1 + L2 + V + time
  )),
  ymodel = Y_bin ~ A + L1 + L2, outcome_name = "Y_bin",
  outcome_type = "binary_eof",
  histories = c(gfoRmula::lagged), histvars = list(c("A","L1","L2")),
  intervention1.A = list(gfoRmula::static, rep(1, length(unique(dat$time)))),
  int_descript = "Always treat",
  nsimul = 50000, seed = 20250915
)

cm_vals <- round(fit_cv_cm$effect_size$Est, 5)
gf_vals <- round(fit_cv_gf$result[k == max(fit_cv_gf$result$k)][["g-form mean"]], 5)

data.frame(
  Intervention = c("natural", "always"),
  causalMed    = cm_vals,
  gfoRmula     = gf_vals,
  Difference   = round(abs(cm_vals - gf_vals), 5)
)
#>   Intervention causalMed gfoRmula Difference
#> 1      natural   0.23445  0.23507    0.00062
#> 2       always   0.25441  0.25448    0.00007

The Difference column is what this run produced. The package’s test suite runs the same comparison against gfoRmula for a binary end-of-follow-up outcome, a survival outcome and a dynamic intervention, and requires agreement within 0.005. That agreement holds for those data-generating processes and model specifications, and does not by itself extend to others.

References

Lin, S. H., Young, J. G., Logan, R., & VanderWeele, T. J. (2017). Mediation analysis for a survival outcome with time-varying exposures, mediators, and confounders. Statistics in Medicine, 36(26), 4153–4166.
McGrath, S., Lin, V., Zhang, Z., et al. (2020). gfoRmula: An R Package for Estimating the Effects of Sustained Treatment Strategies via the Parametric g-Formula. Patterns, 1, 100008.
Zheng, W., & van der Laan, M. (2017). Longitudinal mediation analysis with time-varying mediators and exposures, with application to survival outcomes. Journal of Causal Inference, 5(2).

Package versions used for this comparison:

c(causalMed = as.character(packageVersion("causalMed")),
  gfoRmula  = if (has_gfoRmula) as.character(packageVersion("gfoRmula")) else NA)
#> causalMed  gfoRmula 
#>   "0.1.3"   "1.1.1"