gcompurrA user-oriented R package for g-computation

Early gcompurr demo

Here we provide a static demo to illustrate an early version of gcompurr. The package is not yet publicly available, but a development version will be released in the upcoming weeks. A full public release is planned for early 2027.

The main function is gcomp(), which follows the logic of g-computation in this order:

Define the exposure(s) and the outcome
Specify the exposure regimes (interventions)
Specify the causal contrast
Define adjustment sets for each exposure–outcome relationship
Specify models for intermediate confounders
Specify the outcome model
Specify how to get confidence intervals and p-values
Reduce simulation noise with Monte Carlo expansion

This page shows one complete call first, then builds it up in that order. Each step has a short explanation. The expandable boxes hold extensions you can skip on a first read.

Minimal working example

A simple longitudinal setting:

l0: baseline confounder (continuous)
a0, a1: exposure at two timepoints (binary)
l1: intermediate confounder (continuous), affected by a0 and a cause of a1 and y
y: outcome (continuous)

We aim to estimate the causal mean difference comparing the mean potential outcome if everyone were exposed at both timepoints (a0 = 1, a1 = 1) with the mean potential outcome if no one were exposed at either timepoint (a0 = 0, a1 = 0).

The data are simulated as follows:

set.seed(42)
n <- 500

l0 <- rnorm(n)
a0 <- rbinom(n, 1, plogis(0.3 * l0))
l1 <- 0.5 * l0 + 0.8 * a0 + rnorm(n, sd = 0.5)
a1 <- rbinom(n, 1, plogis(0.2 * l0 + 0.5 * a0 + 0.4 * l1))
y  <- 1 + 0.5 * l0 + 1.5 * a0 + 0.8 * l1 + 1.2 * a1 + rnorm(n, sd = 1)

dat <- data.frame(l0, a0, l1, a1, y)
head(dat)

          l0 a0          l1 a1        y
1  1.3709584  0 -0.11855086  1 2.189256
2 -0.5646982  0 -1.25724098  1 0.776042
3  0.3631284  0  0.01121425  1 1.403263
4  0.6328626  1  1.20379407  1 5.811392
5  0.4042683  1 -0.13675462  0 1.797671
6 -0.1061245  0  0.09219998  0 1.361162

Because the data are simulated and the scenario is simple, the true effect can be obtained analytically: 1.5 (direct effect of a0) + 0.8 × 0.8 (through l1) + 1.2 (a1) = 3.34.

The full call

fit <- gcomp(
  data = dat,
  outcome = "y",
  exposure = c("a0", "a1"),
  reference = list(a0 = 0, a1 = 0),
  alternative = list(a0 = 1, a1 = 1),
  contrast = mean_difference(),
  adjustment_sets = list(
    a0 = "l0",
    a1 = c("l0", "a0", "l1")
  ),
  confounder_models = list(
    l1 = "linear"
  ),
  outcome_model = "linear",
  inference = "bootstrap",
  nboot = 500,
  mc_size = 100000,
  seed = 123
)

print(fit)


── G-computation results 
Outcome: "y" (linear)
Exposure: a0 → a1
Simulated confounders: "l1"
N: 500
Monte Carlo expansion: 200 copies per individual (100000 simulated rows)


── Mean potential outcomes 
reference (all set to 0): 1.0498
alternative (all set to 1): 4.3059


── Estimated effects 
Mean difference (alternative vs reference): 3.2561


── Bootstrap inference (500 replicates) 
Method: Percentile bootstrap
95% CI: [2.993, 3.4909]
Bootstrap SE: 0.1269
P-value: < 0.001

Building the call

1. Define the exposure(s) and the outcome

gcomp(
  data = dat,
  outcome = "y",
  exposure = c("a0", "a1")
  # ...
)

exposure = c("a0", "a1") makes the temporal order explicit in one place.
The exposure sequence is the backbone of the longitudinal structure. Further arguments rely on this ordering.

2. Specify the exposure regimes (interventions)

gcomp(
  data = dat,
  outcome = "y",
  exposure = c("a0", "a1"),
  reference = list(a0 = 0, a1 = 0),
  alternative = list(a0 = 1, a1 = 1)
  # ...
)

One value per exposure, keyed by variable name, so the same interface works when the exposure variable differs across timepoints.

Extensions: more intervention types

Fixed values are the simplest regime. Other options include:

Modified exposure regimes Regimes based on the natural exposure value, where the exposure is allowed to take the value it would naturally attain in the absence of an intervention.
Dynamic regimes The exact exposure value a person is assigned depends on the history of covariates. For example, only set a1 as treated if l1 > 0.5 for a given person.

# Natural course: every exposure at its natural value
reference = list(a0 = intervention_natural(), a1 = intervention_natural()),

# Shift a continuous exposure up by one SD from its natural value
alternative = list(a0 = intervention_shift(fn = sd), a1 = intervention_shift(fn = sd)),

# Dynamic rule: treat at a1 when l1 exceeds a threshold
alternative = list(a0 = 0, a1 = intervention_dynamic(~ as.numeric(l1 > 0.5))),

# Modified intervention (if a1 were numeric): reduce the natural exposure by 25%
alternative = list(a0 = intervention_dynamic(~ natural * 0.75), a1 = intervention_natural())

If modified or dynamic exposure regimes require simulating the exposure as well, this can be specified using e.g. exposure_models = list(a1 = "logistic").

3. Specify the causal contrast

gcomp(
  data = dat,
  outcome = "y",
  exposure = c("a0", "a1"),
  reference = list(a0 = 0, a1 = 0),
  alternative = list(a0 = 1, a1 = 1),
  contrast = mean_difference()
  # ...
)

contrast says how the two potential-outcome means are compared: mean_difference() or mean_ratio() for any outcome, risk_difference() or risk_ratio() as an alias for binary outcomes.

Extensions: several scenarios and contrasts at once

reference / alternative / contrast is shorthand for two scenarios and one contrast. Users can instead use the full form, which allows to specify any number of named scenarios and a list of contrasts between them:

gcomp(
  data = dat,
  outcome = "y",
  exposure = c("a0", "a1"),
  interventions = list(
    always_treated = list(a0 = 1, a1 = 1),
    never_treated = list(a0 = 0, a1 = 0),
    natural_course = list(a0 = intervention_natural(), a1 = intervention_natural())
  ),
  contrasts = list(
    mean_difference("always_treated", "never_treated"),
    mean_difference("always_treated", "natural_course")
  ),
  exposure_models = list(a1 = "logistic")
  # ...
)

4. Define adjustment sets

gcomp(
  data = dat,
  outcome = "y",
  exposure = c("a0", "a1"),
  reference = list(a0 = 0, a1 = 0),
  alternative = list(a0 = 1, a1 = 1),
  contrast = mean_difference(),
  adjustment_sets = list(
    # Adjustment set for a0 -> y
    a0 = "l0",
    # Adjustment set for a1 -> y
    a1 = c("l0", "a0", "l1")
  )
  # ...
)

One adjustment set per exposure step, for that exposure-outcome relationship. This should typically also include past exposures.
The model formulas for intermediate confounders and the outcome are derived from these sets. By default they have no interactions or non-linearities, step 5 shows how to change that.

5. Specify models for intermediate confounders

gcomp(
  data = dat,
  outcome = "y",
  exposure = c("a0", "a1"),
  reference = list(a0 = 0, a1 = 0),
  alternative = list(a0 = 1, a1 = 1),
  contrast = mean_difference(),
  adjustment_sets = list(
    a0 = "l0",
    a1 = c("l0", "a0", "l1")
  ),
  confounder_models = list(
    l1 = "linear"
  )
  # ...
)

confounder_models lists the variables that need their own model and simulation step, with a model type for each: "linear", "logistic", "poisson", "negbin" or "multinomial".
Variables in the first adjustment set (here l0) are baseline and are kept as observed, they need no entry.
The formula rule is simple: l1 first appears in the adjustment set for a1, so its model uses the adjustment set and exposure of the previous step, here l0 and a0.

Safeguard: post-baseline variables must be classified

Every variable that first appears after the first exposure step must be listed either in confounder_models (simulated) or in confounders_not_simulated (kept at its observed value, which asserts that earlier interventions cannot change it). Leaving one out is an error, not a silent default. To show this without changing the example, add a second post-baseline covariate and forget to classify it:

dat_c2 <- dat
dat_c2$c2 <- rnorm(nrow(dat_c2))

gcomp(
  data = dat_c2,
  outcome = "y",
  exposure = c("a0", "a1"),
  reference = list(a0 = 0, a1 = 0),
  alternative = list(a0 = 1, a1 = 1),
  contrast = mean_difference(),
  adjustment_sets = list(
    a0 = "l0",
    a1 = c("l0", "a0", "l1", "c2")
  ),
  confounder_models = list(
    l1 = "linear"
    # c2 is missing!
  ),
  outcome_model = "linear"
)

Error in `normalise_confounder_models()` at gcompurr/R/spec.R:102:3:
! Variable "c2" appears after the first exposure step in
  `adjustment_sets` but is classified in neither `confounder_models` nor
  `confounders_not_simulated`.
ℹ Every post-baseline covariate must be explicitly classified.
ℹ To model and simulate the variable: `confounder_models = list(c2 = "linear")`
ℹ To keep it at its observed values: `confounders_not_simulated = "c2"`

If c2 should not be simulated, say so:

confounder_models = list(l1 = "linear"),
confounders_not_simulated = "c2"

Extensions: custom formulas with model_spec() and add_term()

model_spec() replaces the auto-generated formula for one variable, and can add truncation bounds for the simulated draws:

confounder_models = list(
  l1 = model_spec("linear", formula = l1 ~ l0 * a0)
)

# A covariate that cannot be negative: draws are truncated at 0
confounder_models = list(
  l1 = model_spec("linear", lower = 0)
)

With many timepoints, formula_updates adds a term to several auto-generated formulas at once:

formula_updates = list(
  add_term(~ l0:a0, to = c("l1", "y")) # add l0:a0 term to models for l1 and y 
)

6. Specify the outcome model

gcomp(
  data = dat,
  outcome = "y",
  exposure = c("a0", "a1"),
  reference = list(a0 = 0, a1 = 0),
  alternative = list(a0 = 1, a1 = 1),
  contrast = mean_difference(),
  adjustment_sets = list(
    a0 = "l0",
    a1 = c("l0", "a0", "l1")
  ),
  confounder_models = list(
    l1 = "linear"
  ),
  outcome_model = "linear"
)

"linear" or "logistic". The outcome formula uses every variable in the adjustment sets plus the exposures.

Extensions: custom outcome families

model_spec() also sets a GLM family for the outcome, for example a log link for a positive, right-skewed outcome. Predictions are averaged on the response scale, so contrasts stay on the outcome’s own scale whatever the link:

outcome_model = model_spec(family = gaussian(link = "log"))

outcome_model = model_spec(
  family = Gamma(link = "log"),
  formula = y ~ l0 + a0 + l1 + a1 + a0:a1
)

7. Bootstrap inference

gcomp(
  data = dat,
  outcome = "y",
  exposure = c("a0", "a1"),
  reference = list(a0 = 0, a1 = 0),
  alternative = list(a0 = 1, a1 = 1),
  contrast = mean_difference(),
  adjustment_sets = list(
    a0 = "l0",
    a1 = c("l0", "a0", "l1")
  ),
  confounder_models = list(
    l1 = "linear"
  ),
  outcome_model = "linear",
  inference = "bootstrap",
  nboot = 500,
  seed = 123
)

inference = "bootstrap" gives percentile confidence intervals, standard errors and p-values. seed makes the simulation and the bootstrap reproducible.

Extensions: bootstrap options

inference = bootstrap(nboot = 2000, ci_type = "normal", ci_level = 0.90, parallel = TRUE)

parallel = TRUE uses furrr and needs an active future::plan(). The replicates are kept in fit$boot_replicates, one column per scenario mean and per contrast, so an interval for any other combination of scenario means can be computed from them.

8. Monte Carlo expansion

gcomp(
  # ... everything above ...
  mc_size = 100000
)

The counterfactual trajectories are simulated, so the estimates carry simulation noise. mc_size clones each person before simulating (here 200 copies each, for a simulated dataset of 100,000 rows), which reduces that noise.
The models are still fitted on the observed data and the estimand is unchanged. mc_copies sets the number of clones directly.

What was fitted

The fitted models are stored in the result, so the analysis can be checked rather than trusted:

formula(fit$fitted_models$l1)

l1 ~ l0 + a0
<environment: 0x10d0c7040>

formula(fit$fitted_models$y)

y ~ l0 + a0 + l1 + a1
<environment: 0x10d0add98>

So are the simulated datasets under each regime:

head(fit$counterfactual_data$alternative, 3)

        l0 a0       l1 a1 .obs_id .copy  .y_pred
1 1.370958  1 1.160758  1       1     1 5.307203
2 1.370958  1 1.864840  1       1     2 5.952968
3 1.370958  1 1.331626  1       1     3 5.463919