Mediation analysis with time-varying exposures, mediators, and confounders
causalMed uses the parametric g-formula to split the
effect of a time-varying exposure into a direct effect and an indirect
effect through one or more mediators, including when the confounders are
themselves affected by earlier exposure. It estimates:
mediation_type = "I", the default; Lin et al. 2017,
VanderWeele & Tchetgen Tchetgen 2017), with one or more mediators
(Yamamuro et al. 2021);mediation_type = "N"; Zheng & van der Laan 2017), by
g-computation or targeted minimum loss-based estimation;gformula(), checked against the CRAN package
gfoRmula (McGrath et al. 2020).Outcomes can be continuous or binary at the end of follow-up, or discrete-time survival. Confidence intervals come from a subject-level bootstrap, which can run in parallel.
# install.packages("remotes")
remotes::install_github("adayim/causalMed")nonsurvivaldata follows 3,000 subjects over five time
points: a baseline covariate V, a binary exposure
A, confounders L1 and L2 that
respond to A, a mediator M, and a binary
outcome Y_bin at the end of follow-up.
Each time-varying variable gets a model, listed in the order the
variables are generated within a time point. Lagged variables are
created inside the simulation with recodes():
init_recode sets them at the first time point and
in_recode updates them at each later one.
library(causalMed)
models <- list(
spec_model(A ~ V + lag1_A + lag1_L1 + lag1_L2 + time,
var_type = "binary", mod_type = "exposure"),
spec_model(L1 ~ V + A + lag1_L1 + time,
var_type = "normal", mod_type = "covariate"),
spec_model(L2 ~ V + A + lag1_L2 + time,
var_type = "binary", mod_type = "covariate"),
spec_model(M ~ V + A + L1 + L2 + lag1_M + time,
var_type = "normal", mod_type = "mediator"),
spec_model(Y_bin ~ V + A + M + L1 + L2,
var_type = "binary", mod_type = "outcome")
)
fit <- mediation(
data = nonsurvivaldata,
id_var = "id",
time_var = "time",
base_vars = "V",
exposure = "A",
outcome = "Y_bin",
models = models,
init_recode = recodes(lag1_A = 0, lag1_L1 = 0, lag1_L2 = 0, lag1_M = 0),
in_recode = recodes(lag1_A = A, lag1_L1 = L1, lag1_L2 = L2, lag1_M = M),
mc_sample = 10000,
R = 100, # bootstrap replicates; kept low here for speed
quiet = TRUE,
seed = 2025
)
fit
#> Call:
#> mediation(data = nonsurvivaldata, id_var = "id", base_vars = "V",
#> exposure = "A", outcome = "Y_bin", time_var = "time", models = models,
#> init_recode = recodes(lag1_A = 0, lag1_L1 = 0, lag1_L2 = 0,
#> lag1_M = 0), in_recode = recodes(lag1_A = A, lag1_L1 = L1,
#> lag1_L2 = L2, lag1_M = M), mc_sample = 10000, R = 100,
#> quiet = TRUE, seed = 2025)
#>
#> --- Analysis setup ---
#> Exposure : A
#> Mediator(s) : M
#> Outcome : Y_bin [mean outcome at t = 4, end of follow-up]
#> Time variable: time (5 time points: 0 ... 4)
#> ID variable : id
#> Baseline vars: V
#> Data : 3,000 individuals, 15,000 observations
#> Observed subjects following a (1 1 1 1 1): 1,177 of 3,000 (39.2%)
#> Observed subjects following a* (0 0 0 0 0): 5 of 3,000 (0.2%)
#> MC sample : 10000
#> Bootstrap R : 100
#> n_vw : 2 (permutation draws averaged per pool-drawing intervention)
#> Seed : 2025
#> Mediation : Interventional effects (IDE/IIE) -- Lin et al. (2017)
#>
#> --- Marginal mean outcome per intervention ---
#> Under interventional effects, each intervention draws its mediators from independently-permuted pools (G):
#> Phi11 = E[Y(a=1, G1)]: exposure=1, mediators ~ a=1 pool [reference]
#> Phi10 = E[Y(a=1, G0)]: exposure=1, mediators ~ a=0 pool [cross-regime]
#> Phi00 = E[Y(a=0, G0)]: exposure=0, mediators ~ a=0 pool [reference]
#> nat1/nat0 = E[Y(a=1)]/E[Y(a=0)]: exposure fixed, mediators natural (used for the total effect)
#> Intervention Est Sd 2.5%(pct) 97.5%(pct) 2.5%(norm) 97.5%(norm)
#> <char> <num> <num> <num> <num> <num> <num>
#> 1: nat0 0.0961 0.0151 0.0695 0.1262 0.0664 0.1257
#> 2: nat1 0.2573 0.0090 0.2400 0.2719 0.2396 0.2750
#> 3: Phi00 0.0773 0.0132 0.0547 0.1046 0.0514 0.1032
#> 4: Phi10 0.1652 0.0102 0.1501 0.1883 0.1453 0.1851
#> 5: Phi11 0.2335 0.0098 0.2157 0.2503 0.2143 0.2527
#> Observed (nonparametric) mean of Y_bin at t = 4 (end of follow-up): 0.2333
#> (informal benchmark; interventions fix the exposure, so exact agreement is not expected)
#>
#> --- Effect decomposition ---
#> Direct effect (IDE) = Phi10 - Phi00
#> Indirect effect (IIE) = Phi11 - Phi10 (sequential per mediator when N>=2)
#> IDE + IIE = Phi11 - Phi00 (interventional overall effect)
#> Total effect (TE) = nat1 - nat0 (natural plug-in g-formula)
#> TE - (Direct+Indirect)= natural TE minus interventional overall effect
#> Mediation Prop. = Indirect / (Direct + Indirect) (percentage)
#> i.e. a share of the interventional overall effect, NOT of the total
#> effect; Direct and Mediation Prop. sum to 100%
#> RD = risk difference; RR = risk ratio
#> Effect RD RR Sd(RD) RD 2.5%(pct) RD 97.5%(pct)
#> <char> <num> <num> <num> <num> <num>
#> 1: Indirect effect 0.0683 1.4132 0.0085 0.0521 0.0819
#> 2: Direct effect 0.0879 2.1378 0.0159 0.0571 0.1191
#> 3: Total effect 0.1613 2.6787 0.0179 0.1190 0.1872
#> 4: TE - (Direct + Indirect) 0.0051 NA 0.0018 -0.0001 0.0066
#> 5: Mediation Proportion 43.7076 NA 5.9414 34.0057 55.8100
#> Sd(RR) RR 2.5%(pct) RR 97.5%(pct) RD 2.5%(norm) RD 97.5%(norm) RR 2.5%(norm)
#> <num> <num> <num> <num> <num> <num>
#> 1: 0.0652 1.2857 1.5185 0.0516 0.0849 1.2854
#> 2: 0.3897 1.5371 2.9167 0.0568 0.1190 1.3739
#> 3: 0.4424 1.9739 3.6034 0.1261 0.1964 1.8117
#> 4: NA NA NA 0.0015 0.0086 NA
#> 5: NA NA NA 32.0626 55.3526 NA
#> RR 97.5%(norm)
#> <num>
#> 1: 1.5411
#> 2: 2.9016
#> 3: 3.5457
#> 4: NA
#> 5: NA
#>
#> 95% CIs: percentile (pct) and normal approximation (norm) from 100 bootstrap replicates.By default the effects compare always exposed with never exposed;
other static regimes are set with exposure_regime and
reference_regime. The rows of the decomposition are:
| Row | Meaning |
|---|---|
| Direct effect | Phi10 - Phi00: effect of the exposure with the mediator
drawn from its distribution under no exposure |
| Indirect effect | Phi11 - Phi10: with exposure, the effect of shifting
the mediator’s distribution from its unexposed to its exposed form |
| Total effect | nat1 - nat0: the ordinary g-formula total effect |
| TE - (Direct + Indirect) | the decomposition residual: direct + indirect sum to the
interventional overall effect Phi11 - Phi00, not to the
total effect |
| Mediation Proportion | indirect effect as a percentage of the interventional overall effect |
Natural effects. Set
mediation_type = "N". Zheng & van der Laan (2017, Lemma
1) identify these effects under sequential randomization and positivity;
reading them as individual-level natural effects additionally requires a
cross-world assumption that is not expected to hold when a
mediator-outcome confounder is affected by the exposure (Avin, Shpitser
& Pearl 2005), as L1 and L2 are here.
mediation() warns when a covariate model includes the
exposure. estimator = "tmle" gives the targeted estimator
of Zheng & van der Laan (2017, Section 4.3); see
?mediation for what it accepts.
Total effects. gformula() needs no
mediator model, and takes a named list of static (1,
0) or dynamic (dyn_int()) interventions:
models_te <- c(models[1:3], # exposure and confounder models from above
list(spec_model(Y_bin ~ V + A + L1 + L2,
var_type = "binary", mod_type = "outcome")))
fit_te <- gformula(
data = nonsurvivaldata, id_var = "id", time_var = "time",
base_vars = "V", exposure = "A", models = models_te,
intervention = list(natural = NULL, always = 1, never = 0,
treat_if_L1_pos = dyn_int(as.numeric(lag1_L1 > 0))),
init_recode = recodes(lag1_A = 0, lag1_L1 = 0, lag1_L2 = 0),
in_recode = recodes(lag1_A = A, lag1_L1 = L1, lag1_L2 = L2),
R = 500
)Parallel bootstrap. Set a future plan
before the call:
future::plan(future::multisession)
fit <- mediation(..., R = 500)
future::plan(future::sequential)Data must be in long format, one row per subject per time point, with a numeric time variable. For a survival outcome, remove every row after the event and after loss to follow-up. The simulation starts from the subject identifier and the baseline covariates only, so lags and other derived variables must be created with the recode hooks even if they already exist in the data.
vignette("causalMed-01-overview"): data format, model
specification and quick starts.vignette("causalMed-02-mediation"): the two estimands,
survival outcomes, multiple mediators, censoring, and natural
effects.vignette("causalMed-03-gformula"): total effects,
dynamic interventions, custom covariate distributions, and a published
replication.vignette("causalMed-04-vs-gfoRmula"): comparison with
gfoRmula.