Causal Mediation Analysis with Time-Varying Exposures, Mediators, and Confounders
causalMed implements the parametric g-formula for total effect
estimation and extends it with the survival mediational g-formula to
decompose causal effects into direct and indirect components. The
package supports:
- Interventional direct and indirect effects (IDE/IIE) — Lin et al. (2017)
- Natural direct and indirect effects (NDE/NIE) — Zheng & van der Laan (2017)
Both approaches handle time-varying exposures, mediators, and
confounders in longitudinal data, including survival outcomes. The
standard g-formula component (total effect estimation) is
cross-validated against the gfoRmula CRAN package.
# Install the development version from GitHub:
# install.packages("devtools")
devtools::install_github("adayim/causalMed")- G-formula for total effects: Monte Carlo simulation under user-defined static or dynamic (threshold) interventions for binary, continuous, and survival outcomes
- Mediational g-formula: decomposes total effects into direct and indirect pathways with time-varying mediators and exposure-induced mediator–outcome confounders
- Two mediation estimands:
mediation_type = "I": interventional IDE/IIE (VanderWeele & Tchetgen Tchetgen 2017; Lin et al. 2017; Yamamuro et al. 2021) — cross-world mediator drawn as a joint M(1:T) trajectory by row-permuting the reference (a*) cohort, matching the SASmGFORMULAmacro and the algorithm of Yamamuro et al. 2021. Does not require cross-world independence. Multiple mediator models are supported under this estimand (Yamamuro 2021 §4); per-mediator IIE(M_k) is reported alongside the overall IDE/TE.mediation_type = "N": natural NDE/NIE (Zheng & van der Laan 2017) — conditional mediator distribution with exposure swapping; requires stronger assumptions (single mediator only). Natural effects are not identifiable when a confounder of the mediator–outcome relationship is itself affected by the exposure;mediation()detects this, warns, and repeats the caveat inprint()— prefermediation_type = "I"in that setting.
- Two estimators for natural effects: the parametric g-formula
plug-in (
estimator = "gcomp", the default) or a targeted maximum likelihood estimator (estimator = "tmle", Zheng & van der Laan 2017 §4.3) — multiply robust (consistent when only certain subsets of the nuisance models are correct), with Wald confidence intervals from the efficient influence curve, so no bootstrap is needed. Available formediation_type = "N"only. - Flexible model specification: logistic regression (binary), linear
regression (normal), multinomial logistic (categorical), and custom
simulation functions. Numeric draws are clipped to the observed range
of the response by default; use
spec_model(truncate = FALSE)to draw from the untruncated fitted distribution. - Recode hooks (
init_recode,in_recode,out_recode) for lag creation, cumulative variables, and other within-loop transformations - Bootstrap confidence intervals: non-parametric bootstrap resampling individuals (preserving longitudinal structure), with percentile and normal-approximation CIs
- Parallel bootstrap: plug in
future::plan(multisession)on Windows orfuture::plan(multicore)on Unix to parallelise across bootstrap replicates
library(causalMed)
data("nonsurvivaldata", package = "causalMed")
# Manage lagged variables
init_rc <- recodes(lag1_A = 0, # At t=0, all lags initialised to 0
lag1_L1 = 0,
lag1_L2 = 0)
in_rc <- recodes(lag1_A = A, # At each subsequent step, copy current values
lag1_L1 = L1,
lag1_L2 = L2)
# Specify models in temporal order: A → L1 → L2 → Y (exposure first,
# current A in the confounder models — matching the documented DGP)
m_A <- spec_model(A ~ V + lag1_A + lag1_L1 + lag1_L2 + time,
var_type = "binary", mod_type = "exposure")
m_L1 <- spec_model(L1 ~ V + A + lag1_L1 + time,
var_type = "normal", mod_type = "covariate")
m_L2 <- spec_model(L2 ~ V + A + lag1_L2 + time,
var_type = "binary", mod_type = "covariate")
m_Y <- spec_model(Y_bin ~ V + A + L1 + L2,
var_type = "binary", mod_type = "outcome")
models_bin <- list(m_A, m_L1, m_L2, m_Y)
# Define intervention strategies
# NULL = natural course (draw exposure from its fitted model)
# 1 / 0 = always treat / never treat
ints <- list(natural = NULL, always_treat = 1, never_treat = 0)
# Run g-formula
fit_bin <- gformula(
data = nonsurvivaldata,
id_var = "id",
time_var = "time",
base_vars = "V",
exposure = "A",
models = models_bin,
intervention = ints,
ref_int = "natural",
init_recode = init_rc,
in_recode = in_rc,
mc_sample = 10000,
R = 100, # set R > 1 for bootstrap CIs; kept low here for speed
quiet = TRUE,
seed = 2025
)
#> Warning: package 'future' was built under R version 4.5.3
print(fit_bin)
#> Call:
#> gformula(data = nonsurvivaldata, id_var = "id", base_vars = "V",
#> exposure = "A", time_var = "time", models = models_bin, intervention = ints,
#> ref_int = "natural", init_recode = init_rc, in_recode = in_rc,
#> mc_sample = 10000, R = 100, quiet = TRUE, seed = 2025)
#>
#> --- Analysis setup ---
#> Exposure : A
#> 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
#> MC sample : 10000
#> Bootstrap R : 100
#> Seed : 2025
#> Reference : natural
#>
#> --- Mean outcome by intervention ---
#> Intervention Est Sd 2.5%(pct) 97.5%(pct) 2.5%(norm) 97.5%(norm)
#> <fctr> <num> <num> <num> <num> <num> <num>
#> 1: natural 0.1413 0.0066 0.1297 0.1548 0.1282 0.1543
#> 2: always_treat 0.1511 0.0075 0.1380 0.1659 0.1365 0.1657
#> 3: never_treat 0.0856 0.0136 0.0611 0.1145 0.0590 0.1123
#> Observed (nonparametric) mean of Y_bin at t = 4 (end of follow-up): 0.1407
#> (informal model check: compare with the natural-course intervention)
#>
#> --- Contrasts vs. reference intervention ---
#> Intervention Risk_type Estimate Sd 2.5%(pct) 97.5%(pct)
#> <char> <char> <num> <num> <num> <num>
#> 1: always_treat - natural Difference 0.0098 0.0024 0.0054 0.0142
#> 2: always_treat / natural Ratio 1.0694 0.0165 1.0394 1.1004
#> 3: never_treat - natural Difference -0.0556 0.0130 -0.0768 -0.0266
#> 4: never_treat / natural Ratio 0.6063 0.0910 0.4530 0.8081
#> 2.5%(norm) 97.5%(norm)
#> <num> <num>
#> 1: 0.0052 0.0145
#> 2: 1.0370 1.1018
#> 3: -0.0811 -0.0302
#> 4: 0.4279 0.7846
#>
#> 95% CIs: percentile (pct) and normal approximation (norm) from 100 bootstrap replicates.The list order must match your assumed data-generating process — here
the exposure is decided first within each period and the confounders
respond to it, as documented in ?nonsurvivaldata.
library(causalMed)
data("nonsurvivaldata", package = "causalMed")
# Model list must include a mediator model (mod_type = "mediator")
# List order must reflect the temporal DAG: A -> L -> M -> Y
init_med <- recodes(lag1_A = 0, lag1_L1 = 0, lag1_L2 = 0, lag1_M = 0)
in_med <- recodes(lag1_A = A, lag1_L1 = L1, lag1_L2 = L2, lag1_M = M)
models_med <- list(
spec_model(A ~ V + lag1_L1 + lag1_L2 + lag1_A + 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"), # <-- mediator
spec_model(Y_bin ~ V + A + M + L1 + L2,
var_type = "binary", mod_type = "outcome")
)
fit_med <- mediation(
data = nonsurvivaldata,
id_var = "id",
time_var = "time",
base_vars = "V",
exposure = "A",
outcome = "Y_bin",
models = models_med,
init_recode = init_med,
in_recode = in_med,
mediation_type = "I", # Interventional IDE/IIE
mc_sample = 10000,
R = 100,
quiet = TRUE,
seed = 2025
)
print(fit_med)
#> Call:
#> mediation(data = nonsurvivaldata, id_var = "id", base_vars = "V",
#> exposure = "A", outcome = "Y_bin", time_var = "time", models = models_med,
#> init_recode = init_med, in_recode = in_med, mc_sample = 10000,
#> mediation_type = "I", 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
#> MC sample : 10000
#> Bootstrap R : 100
#> n_vw : 2 (permutation draws averaged per cross-world 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-world]
#> Phi00 = E[Y(a=0, G0)]: exposure=0, mediators ~ a=0 pool [reference]
#> nat1/nat0 = E[Y(a=1)]/E[Y(a=0)]: natural course (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.0847 0.0135 0.0607 0.1128 0.0582 0.1112
#> 2: nat1 0.1519 0.0075 0.1383 0.1661 0.1371 0.1667
#> 3: Phi00 0.0848 0.0135 0.0608 0.1126 0.0584 0.1112
#> 4: Phi10 0.1456 0.0074 0.1329 0.1586 0.1312 0.1600
#> 5: Phi11 0.1517 0.0075 0.1389 0.1662 0.1370 0.1664
#> Observed (nonparametric) mean of Y_bin at t = 4 (end of follow-up): 0.1407
#> (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)= mediated-interaction residual (TE - overall)
#> Mediation Prop. = (Total - Direct) / Total (percentage; RR not applicable)
#> RD = risk difference; RR = risk ratio
#> Effect RD RR Sd(RD) RD 2.5%(pct)
#> <char> <num> <num> <num> <num>
#> 1: Indirect effect 0.0061 1.0421 0.0021 0.0025
#> 2: Direct effect 0.0608 1.7168 0.0150 0.0287
#> 3: Total effect 0.0672 1.7931 0.0153 0.0323
#> 4: TE - (Direct + Indirect) 0.0003 NA 0.0004 -0.0009
#> 5: Mediation Proportion 9.5299 NA 3.7089 4.0617
#> 6: Mediation Proportion (multiplicative) 9.1585 NA 3.7171 4.0309
#> RD 97.5%(pct) Sd(RR) RR 2.5%(pct) RR 97.5%(pct) RD 2.5%(norm) RD 97.5%(norm)
#> <num> <num> <num> <num> <num> <num>
#> 1: 0.0108 0.0149 1.0178 1.0732 0.0020 0.0103
#> 2: 0.0856 0.2924 1.2632 2.3606 0.0313 0.0903
#> 3: 0.0916 0.3076 1.2961 2.4283 0.0371 0.0973
#> 4: 0.0007 NA NA NA -0.0005 0.0011
#> 5: 17.6448 NA NA NA 2.2606 16.7991
#> 6: 17.7527 NA NA NA 1.8732 16.4438
#> RR 2.5%(norm) RR 97.5%(norm)
#> <num> <num>
#> 1: 1.0129 1.0713
#> 2: 1.1437 2.2899
#> 3: 1.1903 2.3959
#> 4: NA NA
#> 5: NA NA
#> 6: NA NA
#>
#> 95% CIs: percentile (pct) and normal approximation (norm) from 100 bootstrap replicates.The estimate component of fit_med contains:
| Effect | Definition |
|---|---|
| Indirect effect | Q(1,1) − Q(1,0): effect through the mediator pathway |
| Direct effect | Q(1,0) − Q(0,0): effect not through the mediator |
| Total effect | natural plug-in g-formula contrast E[Y₁] − E[Y₀] |
| TE − (Direct + Indirect) | mediated-interaction residual (interventional only; exactly 0, and omitted, for natural effects) |
| Mediation Proportion | (Total − Direct) / Total × 100% (additive) |
| Mediation Proportion (multiplicative) | RR-scale proportion mediated (Lin et al. 2017, Table 2) |
With R > 1, estimate additionally carries bootstrap SEs and
percentile/normal CIs on the RD and RR scales, and the returned object
gains boot_estimates (the per-replicate draws), data_summary, and an
observed nonparametric benchmark that print() displays next to the
simulated means.
For mediation_type = "I", Q(a₁, a₂) is the mean outcome when exposure
is set to a₁ and each mediator is a stochastic draw from its
marginal distribution under exposure a₂. The direct and indirect effects
then sum to the interventional overall effect Q(1,1) − Q(0,0), which
generally differs from the total effect; their gap is the residual row.
For mediation_type = "N" the references use each subject’s own
mediator and the decomposition sums exactly to the total effect (no
residual row).
Under mediation_type = "I", the mediators in every intervention —
the references Q(0,0) and Q(1,1) as well as the cross-world Q(1,0) — are
drawn as joint trajectories: a natural-course cohort is simulated
under each a*, and the intervention assigns each subject the full
M(1:T) of a randomly permuted reference-cohort individual (each mediator
permuted independently). This matches the SAS mGFORMULA macro (Lin et
al. 2017, eAppendix) and Yamamuro et al. 2021 (Figure 3, step 3).
Mediator values are not survival-weighted; the full reference cohort is
used. See ?mediation for details.
Change one argument:
fit_natural <- mediation(
..., # same arguments as above
mediation_type = "N" # natural NDE/NIE (Zheng & van der Laan 2017)
)Natural effects condition the mediator model on the individual’s own
covariate history but evaluate it at the alternative exposure level
(exposure swapping). They require stronger sequential
no-unmeasured-confounding assumptions than interventional effects. In
particular, natural direct and indirect effects are not identifiable
when a mediator–outcome confounder is itself affected by prior exposure
(Avin, Shpitser & Pearl 2005; VanderWeele & Tchetgen Tchetgen 2017).
mediation() detects this pattern, warns at run time, and restates the
caveat under the decomposition when you print() the result; use
mediation_type = "I" in that setting.
For natural effects only, the plug-in simulation can be replaced by the targeted maximum likelihood estimator of Zheng & van der Laan (2017, §4.3):
fit_tmle <- mediation(
..., # same arguments as above
mediation_type = "N",
estimator = "tmle" # multiply robust; EIC Wald CIs, no bootstrap
)It runs backward targeted regressions instead of forward simulation, so
it is multiply robust (consistent when only certain subsets of the
nuisance models are correct) and reports Wald CIs from the efficient
influence curve — R and mc_sample are ignored. It requires an
exposure model, a binary outcome, and lag-style recodes only
(recodes(lag_A = A); exposure lags must be first-order); derived
recodes such as splines or cumulative counts need estimator = "gcomp".
Inspect any positivity warnings it reports before trusting the affected
quantities.
library(future)
plan(multisession) # Windows; use plan(multicore) on Unix/macOS
fit_parallel <- mediation(..., R = 500)
plan(sequential) # reset after useInput data must be in long format with one row per subject per time
point. The time variable should be an ordered integer starting at 0.
Variables used in models as lags (e.g., lag1_A) must be pre-computed
and present in the data, or created via init_recode / in_recode
hooks. The nonsurvivaldata bundled with the package illustrates the
required structure.
data("nonsurvivaldata", package = "causalMed")
head(nonsurvivaldata)
#> id time V L1 L2 A M Y_cont Y_bin
#> <num> <num> <num> <num> <num> <num> <num> <num> <num>
#> 1: 1 0 0.5218287 0.5527842 0 1 -0.2196759 NA NA
#> 2: 1 1 0.5218287 -0.1758057 1 1 0.4719272 NA NA
#> 3: 1 2 0.5218287 0.5903320 0 1 -0.2237722 NA NA
#> 4: 1 3 0.5218287 1.3174098 0 1 0.3878269 NA NA
#> 5: 1 4 0.5218287 0.7974775 0 1 -0.4163290 0.1348551 0
#> 6: 2 0 0.4066671 0.5492001 1 1 1.0580998 NA NA
#> lag1_A lag1_L1 lag1_L2 lag1_M
#> <num> <num> <num> <num>
#> 1: NA NA NA NA
#> 2: 1 0.5527842 0 -0.2196759
#> 3: 1 -0.1758057 1 0.4719272
#> 4: 1 0.5903320 0 -0.2237722
#> 5: 1 1.3174098 0 0.3878269
#> 6: NA NA NA NA-
Westreich, D., Cole, S. R., Young, J. G., et al. (2012). The parametric g-formula to estimate the effect of HAART on incident AIDS or death. Statistics in Medicine, 31, 2000–2009. doi:10.1002/sim.5316
-
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
-
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
-
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
-
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: Application to atorvastatin’s effect on cardiovascular disease via cholesterol lowering and anti-inflammatory actions in elderly type 2 diabetics. Statistical Methods in Medical Research, 30(8), 1782–1799. doi:10.1177/09622802211025988