mediation_analysis()

Causal mediation analysis with GAM nuisance models.

Usage

Source

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 on outcome operates through an intermediate mediator variable, versus acting directly, controlling for 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(). The mediator model always uses 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 and mediator set to m, and \hat M_i(d) is the mediator GAM’s prediction with treatment set to 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

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%
)