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:
- interventional direct and indirect effects
(
mediation_type = "I", the default; Lin et al. 2017, VanderWeele & Tchetgen Tchetgen 2017), with one or more mediators (Yamamuro et al. 2021); - natural direct and indirect effects (
mediation_type = "N"; Zheng & van der Laan 2017), by g-computation or targeted minimum loss-based estimation; - total effects of static or dynamic interventions with
gformula(), checked against the CRAN packagegfoRmula(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 withgfoRmula.
- Avin, C., Shpitser, I., & Pearl, J. (2005). Identifiability of path-specific effects. Proceedings of the 19th International Joint Conference on Artificial Intelligence, 357–363.
- 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. doi:10.1002/sim.7426
- 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. doi:10.1016/j.patter.2020.100008
- VanderWeele, T. J., & Tchetgen Tchetgen, E. J. (2017). Mediation analysis with time varying exposures and mediators. Journal of the Royal Statistical Society: Series B, 79(3), 917–938. doi:10.1111/rssb.12194
- Yamamuro, S., Shinozaki, T., Iimuro, S., & Matsuyama, Y. (2021). Mediational g-formula for time-varying treatment and repeated-measured multiple mediators. Statistical Methods in Medical Research, 30(8), 1782–1799. doi:10.1177/09622802211025988
- 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). doi:10.1515/jci-2016-0006