--- title: "Causal Mediation Analysis with causalMed" output: rmarkdown::html_vignette: toc: true toc_depth: 3 highlight: null vignette: > %\VignetteIndexEntry{Causal Mediation Analysis with causalMed} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- ```{r setup, include=FALSE} knitr::opts_chunk$set( collapse = TRUE, comment = "#>", fig.width = 7, fig.height = 4 ) ``` ## Introduction `mediation()` splits the effect of a time-varying exposure into a **direct** effect and an **indirect** effect through one or more time-varying mediators, including when confounders are affected by earlier exposure, the setting the mediational g-formula was developed for (Lin et al. 2017; VanderWeele & Tchetgen Tchetgen 2017). This vignette works through a survival example and reads every row of the output, then covers other exposure regimes, multiple mediators, censoring, bootstrap intervals and natural effects. Data format, `spec_model()` and the `recodes()` hooks are introduced in `vignette("causalMed-01-overview")`. ```{r load, message=FALSE, warning=FALSE} library(causalMed) library(data.table) ``` ## The two estimands `mediation_type` chooses between two definitions of the direct and indirect effects. **Natural effects** (`"N"`; Zheng & van der Laan 2017) decompose the total effect exactly: $$\underbrace{E[Y_{1,M(1)}] - E[Y_{0,M(0)}]}_{\text{Total effect}} = \underbrace{E[Y_{1,M(0)}] - E[Y_{0,M(0)}]}_{\text{Direct effect}} + \underbrace{E[Y_{1,M(1)}] - E[Y_{1,M(0)}]}_{\text{Indirect effect}}$$ Here the mediator is drawn from its **conditional** distribution: at each time point it is predicted from the subject's own history, with the exposure and its first-order lags set to the other level. **Interventional effects** (`"I"`, the default; Lin et al. 2017) set the mediator to a random draw $G_{a^*}$ from its **marginal** distribution under $a^*$, drawn independently for each mediator. The direct and indirect effects then sum to the **interventional overall effect** $E[Y_{1,G_1}] - E[Y_{0,G_0}]$, which generally differs from the total effect $E[Y_1] - E[Y_0]$. `mediation()` reports the total effect separately and adds the difference as a row, `TE - (Direct + Indirect)`, the decomposition residual. | | Interventional (`"I"`) | Natural (`"N"`) | |--|----------------------|-----------------| | Mediator drawn from | Marginal distribution (permutation) | Conditional distribution (own history) | | Cross-world assumption needed | No | Only to read the effects as individual-level | | Identified with exposure-affected confounders | Yes | Yes, as the conditional-draw estimand | | Several mediators | Yes (Yamamuro et al. 2021) | No | | Reference | Lin et al. (2017) | Zheng & van der Laan (2017) | Zheng & van der Laan (2017, Lemma 1) identify the `"N"` effects under sequential randomization and positivity. Reading them as **individual-level** natural effects, contrasts of each subject's own counterfactual mediator, additionally requires a cross-world assumption that is not expected to hold when a mediator-outcome confounder is itself affected by the exposure (Avin, Shpitser & Pearl 2005; VanderWeele & Tchetgen Tchetgen 2017). For that setting VanderWeele & Tchetgen Tchetgen (2017) propose the interventional effects, which is why `"I"` is the default. With `"N"`, `mediation()` warns when a covariate model includes the exposure; the check reads formulas only and cannot confirm or rule out such a confounder. ## Requirements - **Binary exposure**, coded 0/1. - **At least one mediator model** (`mod_type = "mediator"`). Several mediators, in temporal order, are supported under `"I"` only. - **Model order** follows the order in which variables are generated within a time point, typically `A → L → M → Y`, or `A → M → L → Y` when the confounders respond to the mediator (the setting of Lin et al. 2017). A warning is given if the exposure follows the mediator or the mediator follows the outcome. - **Outcome**: binary or continuous at the end of follow-up (`mod_type = "outcome"`) or a discrete-time event (`mod_type = "survival"`). ## Survival data With a survival outcome, a `mod_type = "survival"` model estimates the discrete-time hazard among subjects still at risk. The simulation accumulates it into the cumulative incidence $1 - \prod_t (1 - h_t)$, so every reported quantity is a **risk of the event by the end of follow-up**. The data must have one row per subject per period at risk, with no rows after the event. ```{r data} data("survivaldata") dat <- as.data.table(survivaldata) head(dat, 8) ``` `survivaldata` follows 3,000 subjects over five periods (`time` 0 to 4): | Variable | Role | |----------|------| | `id`, `time` | Subject and time index | | `V` | Baseline covariate | | `A` | Binary exposure | | `L` | Continuous confounder | | `M` | Binary mediator | | `Y` | Event indicator | | `C` | Loss to follow-up | | `lag1_A`, `lag1_L`, `lag1_M` | Previous-period values | Within each period the order is **A → L → M → Y** (see `?survivaldata`). ## Single mediator: interventional effects The models follow that order, with a censoring model for loss to follow-up (see *Censoring* below) and a hazard model that may include exposure-mediator interactions: ```{r models} init_s <- recodes(lag1_A = 0, lag1_L = 0, lag1_M = 0) in_s <- recodes(lag1_A = A, lag1_L = L, lag1_M = M) models_surv <- list( spec_model(A ~ V + lag1_A + lag1_L + time, var_type = "binary", mod_type = "exposure"), spec_model(L ~ V + A + lag1_L + time, var_type = "normal", mod_type = "covariate"), spec_model(M ~ V + A + L + lag1_M + time, var_type = "binary", mod_type = "mediator"), spec_model(C ~ V + L + time, var_type = "binary", mod_type = "censor"), # loss to follow-up spec_model(Y ~ V + A + M + L + A:M + time, var_type = "binary", mod_type = "survival") # discrete-time hazard ) ``` ```{r fit-single} fit_surv <- mediation( data = dat, id_var = "id", time_var = "time", base_vars = "V", exposure = "A", outcome = "Y", models = models_surv, init_recode = init_s, in_recode = in_s, mediation_type = "I", mc_sample = 20000, R = 1, # set R > 1 for bootstrap CIs (see below) quiet = TRUE, seed = 2026 ) ``` `print(fit_surv)` shows everything below together, with a legend of the interventions, the analysis setup, and an observed (nonparametric) cumulative-incidence benchmark. ### The interventions ```{r interventions-single} fit_surv$effect_size ``` Each row is the simulated cumulative incidence under one intervention: - **`nat0` / `nat1`**: exposure fixed to the reference regime $a^*$ / the exposure regime $a$ (by default never / always exposed), mediator following its fitted model. Their contrast is the **total effect**. - **`Phi00` / `Phi11`**: exposure fixed to $a^*$ / $a$, mediator drawn from the permuted pool of mediator trajectories simulated under that same regime. - **`Phi10`**: exposure fixed to $a$, mediator drawn from the $a^*$ pool. The pools hold each simulated subject's whole mediator trajectory; each subject is given the trajectory of a randomly permuted pool member (see `?mediation`). ### The decomposition ```{r decomp-single} fit_surv$estimate ``` - **Indirect effect** = `Phi11 − Phi10`: with exposure, the change in risk from shifting the mediator's distribution from its unexposed to its exposed form. - **Direct effect** = `Phi10 − Phi00`: the effect of exposure with the mediator drawn from its unexposed distribution. - **Total effect** = `nat1 − nat0`: the ordinary g-formula total effect. - **TE − (Direct + Indirect)**: the direct and indirect effects sum to the interventional overall effect `Phi11 − Phi00`, not to the total effect; this row is the difference. It is absent under `"N"`. - **Mediation Proportion**: the indirect effect as a percentage of the interventional overall effect, the quantity Lin et al. (2017, Table 2) report. It is not a share of the total effect. `RD` is the risk difference and `RR` the risk ratio (`NA` for the last two rows). With `R > 1` both gain bootstrap standard errors and confidence limits. ### Permutation averaging: `n_vw` Each pool-drawing intervention (`Phi00`, `Phi10`, `Phi11`) averages `n_vw` independent permutations of the pool (default 2, as in the SAS `mGFORMULA` macro). `n_vw = 1` is faster, with more Monte Carlo noise. It has no effect under `"N"`, which does not permute. ```{r nvw, eval=FALSE} mediation(..., n_vw = 1) ``` ### Other exposure regimes By default the effects compare always exposed with never exposed, the example Lin et al. (2017, Section 2.2) and Zheng & van der Laan (2017, Section 2.2) give when defining the effects for arbitrary regimes. `exposure_regime` and `reference_regime` take any two static regimes as 0/1 vectors, one element per distinct time point (a scalar is recycled). Exposure at times 1 and 2 only, against never exposed: ```{r regime} fit_window <- mediation( data = dat, id_var = "id", time_var = "time", base_vars = "V", exposure = "A", outcome = "Y", models = models_surv, init_recode = init_s, in_recode = in_s, mediation_type = "I", exposure_regime = c(0, 1, 1, 0, 0), # exposed at t = 1 and 2 only reference_regime = 0, # never exposed mc_sample = 20000, R = 1, quiet = TRUE, seed = 2026 ) fit_window$estimate ``` The rows read as before, with $a$ now the windowed regime. Two points: - The pools hold whole mediator trajectories, so the indirect effect includes the window's effect on the mediator at $t$ = 3 and 4, after the window has closed. - The zeros outside the window force the exposure to 0 there; they do not let it follow its fitted model. That would be a dynamic regime, a different estimand, and is not available. Which pair of regimes answers a given question is the analyst's judgement; the package checks only that they are 0/1, of the right length, and different. `print()` also reports how many observed subjects follow each regime: ```{r regime-support} fit_window$data_summary$regime_support ``` `n_following` counts subjects whose exposure matches the regime at every time they were observed; `n_complete` counts those among subjects observed at every time point. The estimation does not use these counts: the g-formula produces an estimate whether or not anyone followed the regime. ## Multiple mediators: the Yamamuro et al. (2021) simulation Under `"I"`, several mediators are analysed by listing their models in temporal order. `yamamurodata` was simulated from the data-generating process of Yamamuro et al. (2021): a treatment `A`, a confounder `L`, two sequential mediators `M1` and `M2`, and a survival outcome over three visits, ordered **A → L → M1 → M2 → Y**. Its true effects are given in `?yamamurodata`. ```{r yam-data} data("yamamurodata") yam <- as.data.table(yamamurodata) head(yam, 6) ``` Three details of the specification: - **Visit indicators** are written `I(as.integer(time == k))` rather than `factor(time)`: each simulated time step has a single `time` value, so a factor would lose its other levels. - **`subset = time > 0`**: visit-0 values are baseline draws, so the time-varying models apply from visit 1. Visit 0 is set from the baseline columns `L0base`, `M10base`, `M20base` in `init_recode`. - The formulas include the quadratic and lag terms of the published correctly specified scenario. ```{r yam-models} models_yam <- list( spec_model(A ~ V + lag1_A + lag1_L + I(lag1_L^2) + lag1_M1 + lag1_M2 + I(as.integer(time == 2)), var_type = "binary", mod_type = "exposure", subset = time > 0), spec_model(L ~ V + A + lag1_L + I(lag1_L^2) + lag1_M1 + lag1_M2 + I(as.integer(time == 2)), var_type = "normal", mod_type = "covariate", subset = time > 0), spec_model(M1 ~ V + A + L + I(L^2) + lag1_L + I(lag1_L^2) + lag1_M1 + lag1_M2 + I(as.integer(time == 2)), var_type = "normal", mod_type = "mediator", subset = time > 0), spec_model(M2 ~ V + A + L + I(L^2) + lag1_L + I(lag1_L^2) + lag1_M1 + lag1_M2 + I(as.integer(time == 2)), var_type = "normal", mod_type = "mediator", subset = time > 0), spec_model(Y ~ V + A + L + I(L^2) + M1 + M2 + M1:M2 + I(as.integer(time == 1)) + I(as.integer(time == 2)), var_type = "binary", mod_type = "survival") ) fit_yam <- mediation( data = yam, id_var = "id", time_var = "time", base_vars = c("V", "L0base", "M10base", "M20base"), exposure = "A", outcome = "Y", models = models_yam, init_recode = recodes(L = L0base, M1 = M10base, M2 = M20base, lag1_A = 0, lag1_L = 0, lag1_M1 = 0, lag1_M2 = 0), in_recode = recodes(lag1_A = A, lag1_L = L, lag1_M1 = M1, lag1_M2 = M2), mediation_type = "I", mc_sample = 20000, R = 1, quiet = TRUE, seed = 2026 ) fit_yam$effect_size ``` With $N$ mediators there are $4 + N$ interventions: `Phi1_k` draws the first $k$ mediators from the exposed pool and the rest from the unexposed pool. The indirect effect through mediator $k$ is the **sequential contrast** $\Phi_{1,k} - \Phi_{1,k-1}$, labelled `Indirect effect ()`, and the direct effect plus all indirect effects equals `Phi11 − Phi00`. Because the decomposition is sequential, the order of the mediator models matters. Comparing with the true values (in percentage points): ```{r yam-truth} truth <- data.table( Effect = c("Total effect", "Direct effect", "Indirect effect (M1)", "Indirect effect (M2)", "TE - (Direct + Indirect)"), True = c(-6.36, -3.20, -2.29, -0.97, 0.10) # from ?yamamurodata ) est <- as.data.table(fit_yam$estimate)[, .(Effect, Estimate = RD * 100)] cmp <- merge(truth, est, by = "Effect", sort = FALSE) cmp[, Difference := Estimate - True] cmp ``` The true values are large-sample values from the data-generating process. The estimate comes from one dataset of 10,000 subjects and one Monte Carlo run, so the differences combine sampling and Monte Carlo error, neither of which this single run quantifies; a bootstrap (`R > 1`) does. ## Censoring When follow-up can end before the event, add a censoring indicator with `mod_type = "censor"`. To illustrate censoring that depends on the confounder, we add extra loss to follow-up to `survivaldata`, ending each subject's follow-up at the first censoring: ```{r censor-data} set.seed(11) dat_c <- copy(dat) dat_c[, C := rbinom(.N, 1, plogis(-3.5 + 0.8 * L))] dat_c[, after_cens := cumsum(shift(C, fill = 0)), by = id] dat_c <- dat_c[after_cens == 0][, after_cens := NULL] dat_c[C == 1, Y := 0L] # censored before the event in that period ``` The censoring model goes after the variables it depends on and before the survival model: ```{r fit-censor} models_cens <- list( spec_model(A ~ V + lag1_A + lag1_L + time, var_type = "binary", mod_type = "exposure"), spec_model(L ~ V + A + lag1_L + time, var_type = "normal", mod_type = "covariate"), spec_model(M ~ V + A + L + lag1_M + time, var_type = "binary", mod_type = "mediator"), spec_model(C ~ V + A + L + time, var_type = "binary", mod_type = "censor"), # censoring model spec_model(Y ~ V + A + M + L + A:M + time, var_type = "binary", mod_type = "survival") ) fit_cens <- mediation( data = dat_c, id_var = "id", time_var = "time", base_vars = "V", exposure = "A", outcome = "Y", models = models_cens, init_recode = init_s, in_recode = in_s, mediation_type = "I", mc_sample = 20000, R = 1, quiet = TRUE, seed = 2026 ) fit_cens$estimate ``` Under every intervention the censoring indicator is set to 0, so the reported risks are risks in the absence of censoring, as in the g-formula treatment of right-censoring (Robins 1986; Westreich et al. 2012). With the default `estimator = "gcomp"` the censoring model does not change these risks: the hazard model is fitted on the rows still at risk, which identifies them when censoring is independent of the event given the modelled history, an assumption the package cannot check. The censoring model is used by the targeted estimator (see below). ## Absorbing states with `subset` `spec_model(subset = ...)` fits and simulates a model only on the rows that meet a condition. The classic use is an **absorbing state**, such as the GvHD exposure of Keil et al. (2014), which can switch on only once: its model is fitted among the not-yet-exposed, and `out_recode` carries the value forward. ```{r subset, eval=FALSE} # the exposure can occur only while gvhdm1 == 0 ... spec_model(gvhd ~ all + cmv + male + age + ..., var_type = "binary", mod_type = "exposure", subset = gvhdm1 == 0) # ... and stays at 1 afterwards out_recode = recodes(gvhd = ifelse(gvhdm1 == 1, 1, gvhd)) ``` Rows excluded by `subset` keep their current value. The full GvHD analysis (`?gvhd`) is in `vignette("causalMed-03-gformula")`. ## Bootstrap confidence intervals Set `R > 1` for subject-level bootstrap intervals, optionally with a parallel plan: ```{r boot, eval=FALSE} future::plan(future::multisession) fit_ci <- mediation( data = dat, id_var = "id", time_var = "time", base_vars = "V", exposure = "A", outcome = "Y", models = models_surv, init_recode = init_s, in_recode = in_s, mediation_type = "I", mc_sample = 20000, R = 500, seed = 2026 ) future::plan(future::sequential) fit_ci$estimate # with Sd, percentile and normal limits fit_ci$boot_estimates$effects # the per-replicate estimates ``` ## Natural effects `mediation_type = "N"` gives the natural effects of Zheng & van der Laan (2017), for survival outcomes too. Here the confounder `L` responds to the current exposure, so `mediation()` warns: the effects are still identified (their Lemma 1), but not as individual-level natural effects. ```{r natural, warning=TRUE} fit_nat <- mediation( data = dat, id_var = "id", time_var = "time", base_vars = "V", exposure = "A", outcome = "Y", models = models_surv, init_recode = init_s, in_recode = in_s, mediation_type = "N", mc_sample = 20000, R = 1, quiet = TRUE, seed = 2026 ) ``` The choice between `"I"` and `"N"` is substantive. Miles (2023) discusses what a non-zero interventional indirect effect does and does not establish. ### Targeted estimation `estimator = "tmle"` replaces the simulation with the targeted minimum loss-based estimator of Zheng & van der Laan (2017, Section 4.3), which also handles right-censored survival outcomes. It gives Wald intervals from the efficient influence curve, so `R` and `mc_sample` are ignored. It accepts only lag-style recodes and no `subset`, `out_recode` or custom models, raising an error otherwise; `?mediation` lists the requirements. Its targeted regressions are additive main-effects working models, so the multiple robustness Zheng & van der Laan establish for correctly specified nuisance models is not claimed for it. Reusing the censored data and models from above: ```{r tmle-survival} fit_tmle_s <- mediation( data = dat_c, id_var = "id", time_var = "time", base_vars = "V", exposure = "A", outcome = "Y", models = models_cens, init_recode = init_s, in_recode = in_s, mediation_type = "N", estimator = "tmle", quiet = TRUE, seed = 2026 ) fit_tmle_s$estimate ``` - The TMLE changes the estimator, not the estimand, so the warning above still applies; `print()` repeats it and `fit$intermediate_confounders` names the covariates. - Where few subjects follow a regime (e.g. never exposed at every time), the affected targeting steps are skipped with a warning and the estimate relies on the fitted models; inspect these warnings. In a replication of the heavily censored simulation of Zheng & van der Laan (2017, Section 5), this implementation reproduced their value $J^{1,0} \approx 0.912$, but the intervals for the never-treated functionals covered about 87% rather than 95% at n = 4000, while the treated functionals covered at the nominal rate. The subject-level influence-curve values are in `fit_tmle_s$tmle_diag$eic`, for custom contrasts or diagnostics. ## References - Avin, C., Shpitser, I., & Pearl, J. (2005). Identifiability of path-specific effects. *Proceedings of the 19th International Joint Conference on Artificial Intelligence*, 357–363. - Keil, A. P., Edwards, J. K., Richardson, D. B., Naimi, A. I., & Cole, S. R. (2014). The parametric g-formula for time-to-event data: intuition and a worked example. *Epidemiology*, 25(6), 889–897. - 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, 4153–4166. - Miles, C. H. (2023). On the causal interpretation of randomised interventional indirect effects. *Journal of the Royal Statistical Society: Series B*, 85(4), 1154–1172. - Robins, J. M. (1986). A new approach to causal inference in mortality studies with a sustained exposure period. *Mathematical Modelling*, 7(9–12), 1393–1512. - 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. - Westreich, D., Cole, S. R., Young, J. G., et al. (2012). The parametric g-formula to estimate the effect of highly active antiretroviral therapy on incident AIDS or death. *Statistics in Medicine*, 31, 2000–2009. - 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. - 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).