causalMed and gfoRmula: A Comparative Overview
Source:vignettes/causalMed-04-vs-gfoRmula.Rmd
causalMed-04-vs-gfoRmula.RmdOverview
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.
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
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.1029334gfoRmula
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.1047391Side-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.1047Example 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.7849632gfoRmula
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.7849414Side-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.4037Example 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.980955847gfoRmula: 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.00007The 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"