---
title: "causalMed and gfoRmula: A Comparative Overview"
output:
  rmarkdown::html_vignette:
    toc: true
    toc_depth: 3
    highlight: null
vignette: >
  %\VignetteIndexEntry{causalMed and gfoRmula: A Comparative Overview}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include = FALSE}
knitr::opts_chunk$set(
  collapse  = TRUE,
  comment   = "#>",
  message   = FALSE,
  warning   = FALSE,
  fig.width = 7,
  fig.height = 4
)
has_gfoRmula <- requireNamespace("gfoRmula", quietly = TRUE)
```

## 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.

```{r pkgs}
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

```{r data-prep}
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**

```{r cm-binary}
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
```

**gfoRmula**

```{r gf-binary, eval = has_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")]
```

**Side-by-side at final time point**

```{r compare-binary, eval = has_gfoRmula}
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")
```

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

**causalMed**

```{r cm-survival}
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
```

**gfoRmula**

```{r gf-survival, eval = has_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")]
```

**Side-by-side at final time point**

```{r compare-survival, eval = has_gfoRmula}
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")
```

### 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:

```{r dynamic-cm}
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
fit_dyn$estimate
```

**gfoRmula**: a user-written intervention function:

```r
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.

```{r cross-val, eval = has_gfoRmula}
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)
)
```

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

1. 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.

2. 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.

3. 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:

```{r versions}
c(causalMed = as.character(packageVersion("causalMed")),
  gfoRmula  = if (has_gfoRmula) as.character(packageVersion("gfoRmula")) else NA)
```