# mediation_analysis()


Causal mediation analysis with GAM nuisance models.


Usage

``` python
mediation_analysis(
    outcome,
    treatment,
    mediator,
    confounders,
    data,
    *,
    family=None,
    fit_method="REML",
    select=False,
    n_simulations=1000,
    seed=None,
)
```


Estimates how much of the total effect of [treatment](CausalGAM.md#whittaker.CausalGAM.treatment) on [outcome](CausalGAM.md#whittaker.CausalGAM.outcome) operates through an intermediate `mediator` variable, versus acting directly, controlling for [confounders](CausalGAM.md#whittaker.CausalGAM.confounders). Uses the simulation-based approach of Imai, Keele, & Tingley (2010): a mediator GAM `E[M | D, X]` and an outcome GAM `E[Y | D, M, X]` are fit on the observed data, and then used to predict the outcome under counterfactual combinations of treatment and mediator status that isolate the direct and indirect pathways.

Use this when you have a hypothesized causal chain `treatment -> mediator -> outcome` (plus a possible direct `treatment -> outcome` path) and want to decompose the total causal effect into how much passes through the mediator versus how much does not.


## Parameters


`outcome: str`  
Name of the outcome variable.

`treatment: str`  
Name of the treatment variable (binary 0/1).

`mediator: str`  
Name of the mediator variable.

`confounders: list[str]`  
List of confounder variable names, entered as smooth terms in both the mediator and outcome models.

`data: InputData`  
Column-oriented data containing outcome, treatment, mediator, and confounder columns.

`family: Family | None = None`  
Response distribution for the outcome model. Defaults to [Gaussian()](Gaussian.md#whittaker.Gaussian). The mediator model always uses [Gaussian()](Gaussian.md#whittaker.Gaussian).

`fit_method: str = ``"REML"`  
Smoothing parameter selection method for both nuisance models.

`select: bool = ``False`  
Enable double-penalty variable selection in the nuisance models.

`n_simulations: int = ``1000`  
Number of bootstrap resamples used to estimate standard errors for the total, direct, and indirect effects.

`seed: int | None = None`  
Random seed for the bootstrap.


## Notes

Natural direct and indirect effects are computed by contrasting predicted outcomes under three counterfactual scenarios, holding treatment fixed at `d \in \{0, 1\}` and setting the mediator to its predicted value under either treatment level:

\text{indirect} = \frac{1}{n}\sum_i \left\[\hat Y_i(1, \hat M_i(1)) - \hat Y_i(1, \hat M_i(0))\right\], \qquad \text{direct} = \frac{1}{n}\sum_i \left\[\hat Y_i(1, \hat M_i(0)) - \hat Y_i(0, \hat M_i(0))\right\]

where `\hat Y_i(d, m)` is the outcome GAM's prediction with treatment set to [d](TPRS.md#whittaker.TPRS.d) and mediator set to `m`, and `\hat M_i(d)` is the mediator GAM's prediction with treatment set to [d](TPRS.md#whittaker.TPRS.d). The total effect is `indirect + direct`, and `proportion_mediated = indirect / total`. Standard errors for all three quantities come from re-running the full procedure (refitting both GAMs) on `n_simulations` bootstrap resamples of the data.


## Returns


`MediationResult`  
The total, direct, and indirect effects with bootstrap standard errors, and the proportion of the total effect that is mediated.


## Examples


``` python
import numpy as np
from whittaker.causal import mediation_analysis

rng = np.random.default_rng(0)
n = 500
x = rng.uniform(0, 1, n)
d = rng.binomial(1, 0.5, n).astype(float)
m = 0.5 * d + 0.3 * x + rng.normal(scale=0.2, size=n)
y = 0.4 * d + 0.8 * m + np.sin(2 * np.pi * x) + rng.normal(scale=0.3, size=n)

result = mediation_analysis(
    outcome="y", treatment="d", mediator="m", confounders=["x"],
    data={"x": x, "d": d, "m": m, "y": y}, n_simulations=200, seed=0,
)
print(result)
```


    MediationResult(
      total=0.8409 (SE=0.0310)
      direct=0.4257 (SE=0.0472)
      indirect=0.4152 (SE=0.0403)
      proportion_mediated=49.37%
    )