19-mediation.Rmd

# Mediation

## Traditional Approach

The classical mediation analysis follows the approach introduced by @baron1986moderator, though it has limitations, particularly in requiring the first step ($X \to Y$) to be significant. Despite its shortcomings, this framework provides a useful foundation.

### Steps in the Traditional Mediation Model

Mediation is typically assessed through three regression models:

1.  **Total Effect**: $X \to Y$
2.  **Path** $a$: $X \to M$
3.  **Path** $b$ and Direct Effect ($c'$): $X + M \to Y$

where:

-   $X$ = independent (causal) variable
-   $Y$ = dependent (outcome) variable
-   $M$ = mediating variable

Originally, @baron1986moderator required the direct path $X \to Y$ to be significant. However, mediation can still occur even if this direct effect is not significant. For example:

-   The effect of $X$ on $Y$ might be **fully absorbed** by $M$.

-   Multiple mediators ($M_1, M_2$) with **opposing effects** could cancel each other out, leading to a non-significant direct effect.

### Graphical Representation of Mediation

#### Unmediated Model

![Unmediated model](images/mediation/direct_mod.png){width="90%"}

Here, $c$ represents the **total effect** of $X$ on $Y$.

#### Mediated Model

![](images/mediation/full_mod.png){width="90%"}

Here:

-   $c'$ = **direct effect** (effect of $X$ on $Y$ after accounting for mediation)
-   $ab$ = **indirect effect** (mediation pathway)

Thus, we can express:

$$
\text{total effect} = \text{direct effect} + \text{indirect effect}
$$

or,

$$
c = c' + ab
$$

This equation holds under standard linear models but not necessarily in cases such as:

1.  Latent variable models
2.  Logistic regression (only an approximation)
3.  Multilevel models [@bauer2006conceptualizing]

------------------------------------------------------------------------

### Measuring Mediation

Several approaches exist for quantifying the **indirect effect** ($ab$):

1.  **Proportional Reduction Approach**:\
    $$1 - \frac{c'}{c}$$\
    *Not recommended* due to high instability, especially when $c$ is small [@mackinnon1995simulation].

2.  **Product Method**:\
    $$a \times b$$\
    The most common approach.

3.  **Difference Method**:\
    $$c - c'$$\
    Conceptually similar to the product method but less precise in small samples.

------------------------------------------------------------------------

### Assumptions in Linear Mediation Models

For valid mediation analysis, the following assumptions should hold:

1.  **No unmeasured confounders** between $X-Y$, $X-M$, and $M-Y$.
2.  **No reverse causality**: $X$ should not be influenced by a confounder ($C$) that also affects $M-Y$.
3.  **Measurement reliability**: $M$ should be measured without error (if not, consider errors-in-variables models).

------------------------------------------------------------------------

**Regression Equations for Mediation Steps**

Step 1: Total Effect of $X$ on $Y$

$$
Y = \beta_0 + cX + \epsilon
$$

-   The significance of $c$ is **not required** for mediation to occur.

Step 2: Effect of $X$ on $M$

$$
M = \alpha_0 + aX + \epsilon
$$

-   The coefficient $a$ must be significant for mediation analysis to proceed.

Step 3: Effect of $M$ on $Y$ (Including $X$)

$$
Y = \gamma_0 + c'X + bM + \epsilon
$$

-   If $c'$ becomes non-significant after including $M$, full mediation occurs.
-   If $c'$ is reduced but remains significant, partial mediation is present.

| Effect of $X$ on $Y$                            | Mediation Type    |
|-------------------------------------------------|-------------------|
| $b_4$ (from Step 3) is **insignificant**        | Full mediation    |
| $b_4 < b_1$ (from Step 1) but still significant | Partial mediation |

: Interpretation of Mediation Outcomes

------------------------------------------------------------------------

### Testing for Mediation

Several statistical tests exist to assess whether the indirect effect ($ab$) is significant:

1.  **Sobel Test** [@sobel1982asymptotic]
    -   Based on the standard error of $ab$.
    -   **Limitation**: Assumes normality of $ab$, which may not hold in small samples.
2.  **Joint Significance Test**
    -   If both $a$ and $b$ are significant, mediation is likely.
3.  **Bootstrapping (Preferred)** [@preacher2004spss, @shrout2002mediation]
    -   Estimates the confidence interval for $ab$.
    -   Does not assume normality.
    -   Recommended for small-to-moderate sample sizes.

------------------------------------------------------------------------

### Additional Considerations

-   Proximal mediation (where path $a$ exceeds path $b$) can lead to multicollinearity and reduced statistical power. In contrast, distal mediation (where path $b$ exceeds path $a$) tends to maximize power. In fact, slightly distal mediators---where $b$ is somewhat larger than $a$---often strike an ideal balance for power in mediation analyses [@hoyle1999statistical].
-   Tests of direct effects ($c$ and $c'$) generally have lower power than tests of the indirect effect ($ab$). As a result, it is possible for the indirect effect ($ab$) to be statistically significant even when the direct effect ($c$) is not. This situation can appear to indicate "complete mediation," yet the lack of a statistically significant direct effect between $X$ and $Y$ (i.e., $c'$) does not definitively rule out other possibilities [@kenny2014power].
-   Because testing $ab$ essentially combines two tests, it often provides a power advantage over testing $c'$ alone. However, using a non-significant $c'$ as the sole criterion for claiming complete mediation should be done cautiously---if at all---given the importance of adequate sample size and power. Indeed, @hayes2013relative recommend avoiding claims of complete mediation based solely on a non-significant $c'$, particularly when partial mediation may still be present.

------------------------------------------------------------------------

### Assumptions in Mediation Analysis

Valid mediation analysis requires several key assumptions, which can be categorized into **causal direction**, **interaction effects**, **measurement reliability**, and **confounding control**.

------------------------------------------------------------------------

#### Direction

-   Causal Order of Variables

-   A simple but weak solution is to measure $X$ before $M$ and $Y$ to prevent reverse causality (i.e., $M$ or $Y$ causing $X$). Similarly, measuring $M$ before $Y$ avoids feedback effects of $Y$ on $M$.

-   However, causal feedback loops between $M$ and $Y$ may still exist.

    -   If we assume full mediation ($c' = 0$), models with reciprocal causal effects between $M$ and $Y$ can be estimated using instrumental variables (IV).

    -   @smith1982beliefs suggests treating both $M$ and $Y$ as potential mediators of each other, requiring distinct instrumental variables for each to avoid cross-contamination of causal effects.

------------------------------------------------------------------------

#### Interaction Effects in Mediation

-   If $M$ interacts with $X$ in predicting $Y$, then $M$ is both a mediator and a moderator [@baron1986moderator].

-   **The interaction term** $X \times M$ should always be included in the model to account for possible moderation effects.

-   For interpreting such interactions in mediation models, see [@vanderweele2015explanation], who provides a framework for **moderated mediation** analysis.

------------------------------------------------------------------------

#### Reliability

Measurement error in any of the three key variables ($X, M, Y$) can bias estimates of mediation effects.

1.  Measurement Error in the Mediator ($M$):
    -   Biases both $b$ and $c'$.
    -   Potential solution: Model $M$ as a latent variable (reduces bias but may decrease statistical power) [@ledgerwood2011trade].
    -   Specific effects:
        -   $b$ is attenuated (biased toward 0).
        -   $c'$ is:
            -   Overestimated if $ab > 0$.
            -   Underestimated if $ab < 0$.
2.  Measurement Error in the Treatment ($X$):
    -   Biases both $a$ and $b$.
    -   Specific effects:
        -   $a$ is attenuated.
        -   $b$ is:
            -   Overestimated if $ac' > 0$.
            -   Underestimated if $ac' < 0$.
3.  Measurement Error in the Outcome ($Y$):
    -   If unstandardized, there is no bias.
    -   If standardized, there is attenuation bias (reduced effect sizes due to error variance).

------------------------------------------------------------------------

#### Confounding in Mediation Analysis

Omitted variable bias can distort any of the three core relationships ($X \to Y$, $X \to M$, $M \to Y$). Addressing confounding requires either design-based or statistical solutions.

**Design-Based Strategies (Preferred if Feasible)**

-   Randomization of the independent variable ($X$) reduces confounding bias.
-   Randomization of the mediator ($M$), if possible, further strengthens causal claims.
-   Controlling for measured confounders, though this only addresses observable confounding.

**Statistical Strategies (When Randomization is Not Possible)**

1.  Instrumental Variables Approach:
    -   Used when a confounder affects both $M$ and $Y$.
    -   Front-door adjustment can be applied if there exists a third variable that fully mediates the effect of $M$ on $Y$ while being independent of the confounder.
2.  Weighting Methods (e.g., Inverse Probability Weighting - IPW):
    -   Corrects for confounding by reweighting observations to balance confounders across treatment groups.
    -   Requires the strong ignorability assumption: All confounders must be measured and correctly specified [@westfall2016statistically].
    -   While this assumption cannot be formally tested, sensitivity analyses can help assess robustness.
    -   See [Heiss](https://www.andrewheiss.com/blog/2020/12/01/ipw-binary-continuous/) for R code on implementing IPW in mediation models.

------------------------------------------------------------------------

### Indirect Effect Tests

Testing the indirect effect ($ab$) is crucial in mediation analysis. Several methods exist, each with its advantages and limitations.

------------------------------------------------------------------------

#### Sobel Test (Delta Method)

-   Developed by @sobel1982asymptotic.
-   Also known as the delta method.
-   Not recommended because it assumes the sampling distribution of $ab$ is normal, which often does not hold [@mackinnon1995simulation].

The standard error (SE) of the indirect effect is:

$$
SE_{ab} = \sqrt{\hat{b}^2 s_{\hat{a}}^2 + \hat{a}^2 s_{\hat{b}}^2}
$$

The **Z-statistic** for testing whether $ab$ is significantly different from 0 is:

$$
z = \frac{\hat{ab}}{\sqrt{\hat{b}^2 s_{\hat{a}}^2 + \hat{a}^2 s_{\hat{b}}^2}}
$$

**Disadvantages**

-   Assumes $a$ and $b$ are independent.
-   Assumes $ab$ follows a normal distribution.
-   Poor performance in small samples.
-   Lower power and more conservative than bootstrapping.

**Special Case: Inconsistent Mediation**

-   Mediation can occur even when **direct and indirect effects** have **opposite signs**, known as **inconsistent mediation** [@mackinnon2007mediation].
-   This happens when the **mediator acts as a suppressor variable**, leading to counteracting paths.

```{r}
library(bda)
library(mediation)
data("boundsdata")

# Sobel Test for Mediation
bda::mediation.test(boundsdata$med, boundsdata$ttt, boundsdata$out) |>
    tibble::rownames_to_column() |>
    causalverse::nice_tab(2)
```

------------------------------------------------------------------------

#### Joint Significance Test

-   Tests if the indirect effect is nonzero by checking whether both $a$ and $b$ are statistically significant.
-   Assumes $a \perp b$ (independence of paths).
-   Performs similarly to bootstrapping [@hayes2013relative].
-   More robust to non-normality but can be sensitive to heteroscedasticity [@fossum2023use].
-   Does not provide confidence intervals, making effect size interpretation harder.

------------------------------------------------------------------------

#### Bootstrapping (Preferred Method)

-   First applied to mediation by @bollen1990direct.
-   Uses resampling to empirically estimate the sampling distribution of the indirect effect.
-   Does not require normality assumptions or $a \perp b$ independence.
-   Works well with small samples.
-   Can handle complex models.

**Which Bootstrapping Method?**

-   **Percentile bootstrap** is preferred due to better **Type I error rates** [@tibbe2022correcting].
-   **Bias-corrected bootstrapping** can be too **liberal** (inflates Type I errors) [@fritz2012explanation].

**Special Case: Meta-Analytic Bootstrapping**

-   Bootstrapping can be applied **without raw data**, using only $a, b, var(a), var(b), cov(a,b)$ from multiple studies.

```{r}
# Meta-Analytic Bootstrapping for Mediation
library(causalverse)

result <- causalverse::med_ind(
    a = 0.5, 
    b = 0.7, 
    var_a = 0.04, 
    var_b = 0.05, 
    cov_ab = 0.01
)
result$plot
```

When an **instrumental variable (IV)** is available, the causal effect can be estimated more reliably. Below are visual representations.

```{r, eval = FALSE}
library(DiagrammeR)

# Simple Treatment-Outcome Model
grViz("
digraph {
  graph []
  node [shape = plaintext]
    X [label = 'Treatment']
    Y [label = 'Outcome']
  edge [minlen = 2]
    X->Y
  { rank = same; X; Y }
}")

# Mediation Model with an Instrument
grViz("
digraph {
  graph []
  node [shape = plaintext]
    X [label ='Treatment', shape = box]
    Y [label ='Outcome', shape = box]
    M [label ='Mediator', shape = box]
    IV [label ='Instrument', shape = box]
  edge [minlen = 2]
    IV->X
    X->M  
    M->Y 
    X->Y 
  { rank = same; X; Y; M }
}")

```

Mediation Analysis with Fixed Effects Models

```{r}
library(mediation)
library(fixest)

data("boundsdata")

# Step 1: Total Effect (c)
out1 <- feols(out ~ ttt, data = boundsdata)

# Step 2: Indirect Effect (a)
out2 <- feols(med ~ ttt, data = boundsdata)

# Step 3: Direct & Indirect Effect (c' & b)
out3 <- feols(out ~ med + ttt, data = boundsdata)

# Proportion of Mediation
coef(out2)['ttt'] * coef(out3)['med'] / coef(out1)['ttt'] * 100
```

Bootstrapped Mediation Analysis

```{r}
library(boot)
set.seed(1)

# Define the bootstrapping function
mediation_fn <- function(data, i) {
    df <- data[i,]
    
    a_path <- feols(med ~ ttt, data = df)
    a <- coef(a_path)['ttt']
    
    b_path <- feols(out ~ med + ttt, data = df)
    b <- coef(b_path)['med']
    
    cp <- coef(b_path)['ttt']
    
    # Indirect Effect (a * b)
    ind_ef <- a * b
    total_ef <- a * b + cp
    return(c(ind_ef, total_ef))
}

# Perform Bootstrapping
boot_med <- boot(boundsdata, mediation_fn, R = 100, parallel = "multicore", ncpus = 2)
boot_med 

# Summary and Confidence Intervals
summary(boot_med) |> causalverse::nice_tab()

# Confidence Intervals (percentile bootstrap preferred)
boot.ci(boot_med, type = c("norm", "perc"))

# Point Estimates (Indirect and Total Effects)
colMeans(boot_med$t)

```

Alternatively, use the `robmed` package for **robust** mediation analysis:

```{r}
library(robmed)

```

### Power Analysis for Mediation

To assess whether the study has sufficient **power** to detect mediation effects, use:

```{r}
library(pwr2ppl)

# Power analysis for the indirect effect (ab path)
medjs(
    rx1m1 = .3,  # Correlation: X → M (path a)
    rx1y  = .1,  # Correlation: X → Y (path c')
    rym1  = .3,  # Correlation: M → Y (path b)
    n     = 100, # Sample size
    alpha = 0.05,
    mvars = 1,   # Number of mediators
    rep   = 1000 # Replications (use 10,000 for accuracy)
)

```

For **interactive** power analysis, see [Kenny's Mediation Power App](https://davidakenny.shinyapps.io/MedPower/).

**Summary of Indirect Effect Tests**

| **Test**                        | **Pros**                                       | **Cons**                         |
|---------------------|-----------------------------|----------------------|
| **Sobel Test**                  | Simple, fast                                   | Assumes normality, low power     |
| **Joint Significance Test**     | Robust to non-normality                        | No confidence interval           |
| **Bootstrapping** (Recommended) | No normality assumption, handles small samples | May be liberal if bias-corrected |

### Multiple Mediation Analysis

In some cases, a single mediator ($M$) does not fully capture the indirect effect of $X$ on $Y$. Multiple mediation models extend traditional mediation by including two or more mediators, allowing us to examine how multiple pathways contribute to an outcome.

------------------------------------------------------------------------

Several R packages handle multiple mediation models:

-   **manymome**: A flexible package for multiple mediation modeling.
    -   [Vignette](https://cran.r-project.org/web/packages/manymome/vignettes/med_lm.html)

```{r}
library(manymome)
```

-   **mma**: Used for multiple mediator models.

    -   [Package PDF](https://cran.r-project.org/web/packages/mma/mma.pdf)

    -   [Vignette](https://cran.r-project.org/web/packages/mma/vignettes/MMAvignette.html)

```{r}
library(mma)
```

#### Multiple Mediators: Structural Equation Modeling Approach

A popular method for estimating multiple mediation models is **Structural Equation Modeling** using **lavaan**.

To test multiple mediation, we first **simulate data** where two mediators ($M_1$ and $M_2$) contribute to the outcome ($Y$).

```{r}
# Load required packages
library(MASS)  # For mvrnorm (generating correlated errors)
library(lavaan)

# Function to generate synthetic data
generate_data <- function(n = 10000, a1 = 0.5, a2 = -0.35, 
                          b1 = 0.7, b2 = 0.48, 
                          corr = TRUE, correlation_value = 0.7) {
    set.seed(12345)
    X <- rnorm(n)  # Independent variable
    
    # Generate correlated errors for mediators
    if (corr) {
        Sigma <- matrix(c(1, correlation_value, correlation_value, 1), nrow = 2)
        errors <- mvrnorm(n, mu = c(0, 0), Sigma = Sigma) 
    } else {
        errors <- mvrnorm(n, mu = c(0, 0), Sigma = diag(2)) 
    }
    
    M1 <- a1 * X + errors[, 1]
    M2 <- a2 * X + errors[, 2]
    Y  <- b1 * M1 + b2 * M2 + rnorm(n)  # Outcome variable

    return(data.frame(X = X, M1 = M1, M2 = M2, Y = Y))
}
```

We analyze the indirect effects through both mediators ($M_1$ and $M_2$).

1.  Correctly Modeling Correlated Mediators

```{r}
# Generate data with correlated mediators
Data_corr <- generate_data(n = 10000, corr = TRUE, correlation_value = 0.7)

# Define SEM model for multiple mediation
model_corr <- '
  Y ~ b1 * M1 + b2 * M2 + c * X
  M1 ~ a1 * X
  M2 ~ a2 * X
  M1 ~~ M2  # Correlated mediators (modeling correlation correctly)
'

# Fit SEM model
fit_corr <- sem(model_corr, data = Data_corr)

# Extract parameter estimates
parameterEstimates(fit_corr)[, c("lhs", "rhs", "est", "se", "pvalue")]

```

2\. Incorrectly Ignoring Correlation Between Mediators

```{r}
# Define SEM model without modeling mediator correlation
model_uncorr <- '
  Y ~ b1 * M1 + b2 * M2 + c * X
  M1 ~ a1 * X
  M2 ~ a2 * X
'

# Fit incorrect model
fit_uncorr <- sem(model_uncorr, data = Data_corr)

# Compare parameter estimates
parameterEstimates(fit_uncorr)[, c("lhs", "rhs", "est", "se", "pvalue")]

```

**Comparison of Model Fits**

To check whether modeling correlation matters, we compare AIC and RMSEA.

```{r}
# Extract model fit statistics
fit_measures <- function(fit) {
  fitMeasures(fit, c("aic", "bic", "rmsea", "chisq"))
}

# Compare model fits
fit_measures(fit_corr)  # Correct model (correlated mediators)
fit_measures(fit_uncorr)  # Incorrect model (ignores correlation)

```

-   If AIC and RMSEA are lower in the correlated model, it suggests that accounting for correlated errors improves fit.

After fitting the model, we assess:

1.  **Direct Effect**: The effect of $X$ on $Y$ **after accounting for both mediators** ($c'$).

2.  **Indirect Effects**:

    -   $a_1 \times b_1$: Effect of $X \to M_1 \to Y$.

    -   $a_2 \times b_2$: Effect of $X \to M_2 \to Y$.

3.  **Total Effect**: Sum of direct and indirect effects.

```{r}
# Extract indirect and direct effects
parameterEstimates(fit_corr, standardized = TRUE)

```

If $c'$ is reduced (but still significant), we have partial mediation. If $c' \approx 0$, it suggests full mediation.

```{r full code}
# Load required packages
library(MASS)  # for mvrnorm
library(lavaan)

# Function to generate synthetic data with correctly correlated errors for mediators
generate_data <-
  function(n = 10000,
           a1 = 0.5,
           a2 = -0.35,
           b1 = 0.7,
           b2 = 0.48,
           corr = TRUE,
           correlation_value = 0.7) {
    set.seed(12345)
    X <- rnorm(n)
    
    # Generate correlated errors using a multivariate normal distribution
    if (corr) {
        Sigma <- matrix(c(1, correlation_value, correlation_value, 1), nrow = 2)  # Higher covariance matrix for errors
        errors <- mvrnorm(n, mu = c(0, 0), Sigma = Sigma)  # Generate correlated errors
    } else {
        errors <- mvrnorm(n, mu = c(0, 0), Sigma = diag(2))  # Independent errors
    }
    
    M1 <- a1 * X + errors[, 1]
    M2 <- a2 * X + errors[, 2]
    Y <- b1 * M1 + b2 * M2 + rnorm(n)  # Y depends on M1 and M2
    
    data.frame(X = X, M1 = M1, M2 = M2, Y = Y)
}

# Ground truth for comparison
ground_truth <- data.frame(Parameter = c("b1", "b2"), GroundTruth = c(0.7, 0.48))

# Function to extract relevant estimates, standard errors, and model fit
extract_estimates_b1_b2 <- function(fit) {
    estimates <- parameterEstimates(fit)
    estimates <- estimates[estimates$lhs == "Y" & estimates$rhs %in% c("M1", "M2"), c("rhs", "est", "se")]
    estimates$Parameter <- ifelse(estimates$rhs == "M1", "b1", "b2")
    estimates <- estimates[, c("Parameter", "est", "se")]
    fit_stats <- fitMeasures(fit, c("aic", "bic", "rmsea", "chisq"))
    return(list(estimates = estimates, fit_stats = fit_stats))
}

# Case 1: Correlated errors for mediators (modeled correctly)
Data_corr <- generate_data(n = 10000, corr = TRUE, correlation_value = 0.7)
model_corr <- '
  Y ~ b1 * M1 + b2 * M2 + c * X
  M1 ~ a1 * X
  M2 ~ a2 * X
  M1 ~~ M2  # Correlated mediators (errors)
'
fit_corr <- sem(model = model_corr, data = Data_corr)
results_corr <- extract_estimates_b1_b2(fit_corr)

# Case 2: Uncorrelated errors for mediators (modeled correctly)
Data_uncorr <- generate_data(n = 10000, corr = FALSE)
model_uncorr <- '
  Y ~ b1 * M1 + b2 * M2 + c * X
  M1 ~ a1 * X
  M2 ~ a2 * X
'
fit_uncorr <- sem(model = model_uncorr, data = Data_uncorr)
results_uncorr <- extract_estimates_b1_b2(fit_uncorr)

# Case 3: Correlated errors, but not modeled as correlated
fit_corr_incorrect <- sem(model = model_uncorr, data = Data_corr)
results_corr_incorrect <- extract_estimates_b1_b2(fit_corr_incorrect)

# Case 4: Uncorrelated errors, but modeled as correlated
fit_uncorr_incorrect <- sem(model = model_corr, data = Data_uncorr)
results_uncorr_incorrect <- extract_estimates_b1_b2(fit_uncorr_incorrect)

# Combine all estimates for comparison
estimates_combined <- list(
    "Correlated (Correct)" = results_corr$estimates,
    "Uncorrelated (Correct)" = results_uncorr$estimates,
    "Correlated (Incorrect)" = results_corr_incorrect$estimates,
    "Uncorrelated (Incorrect)" = results_uncorr_incorrect$estimates
)

# Combine all into a single table
comparison_table <- do.call(rbind, lapply(names(estimates_combined), function(case) {
    df <- estimates_combined[[case]]
    df$Case <- case
    df
}))

# Merge with ground truth for final comparison
comparison_table <- merge(comparison_table, ground_truth, by = "Parameter")

# Display the comparison table
comparison_table

# Display model fit statistics for each case
fit_stats_combined <- list(
    "Correlated (Correct)" = results_corr$fit_stats,
    "Uncorrelated (Correct)" = results_uncorr$fit_stats,
    "Correlated (Incorrect)" = results_corr_incorrect$fit_stats,
    "Uncorrelated (Incorrect)" = results_uncorr_incorrect$fit_stats
)

fit_stats_combined
```

### Multiple Treatments in Mediation

In some cases, **multiple independent variables** ($X_1$, $X_2$) influence the same mediators. This is called **multiple treatments mediation** [@hayes2014statistical].

For an example in **PROCESS (SPSS/R)**, see:\
[Process Mediation with Multiple Treatments](https://core.ecu.edu/wuenschk/MV/multReg/Mediation_Multicategorical.pdf).

## Causal Inference Approach to Mediation

Traditional mediation models assume that regression-based estimates provide valid causal inference. However, **causal mediation analysis** (CMA) extends beyond traditional models by explicitly defining mediation in terms of **potential outcomes and counterfactuals**.

------------------------------------------------------------------------

### Example: Traditional Mediation Analysis {#sec-example-traditional-mediation-analysis}

We begin with a classic **three-step mediation approach**.

```{r, message=FALSE}
# Load data
myData <- read.csv("data/mediationData.csv")

# Step 1 (Total Effect: X → Y) [No longer required]
model.0 <- lm(Y ~ X, data = myData)
summary(model.0)

# Step 2 (Effect of X on M)
model.M <- lm(M ~ X, data = myData)
summary(model.M)

# Step 3 (Effect of M on Y, controlling for X)
model.Y <- lm(Y ~ X + M, data = myData)
summary(model.Y)

# Step 4: Bootstrapping for ACME
library(mediation)
results <-
    mediate(
        model.M,
        model.Y,
        treat = 'X',
        mediator = 'M',
        boot = TRUE,
        sims = 500
    )
summary(results)
```

-   **Total Effect**: $\hat{c} = 0.3961$ → effect of $X$ on $Y$ **without** controlling for $M$.
-   **Direct Effect (ADE)**: $\hat{c'} = 0.0396$ → effect of $X$ on $Y$ **after** accounting for $M$.
-   **ACME (Average Causal Mediation Effect)**:
    -   ACME = $\hat{c} - \hat{c'} = 0.3961 - 0.0396 = 0.3565$

    -   Equivalent to **product of paths**: $\hat{a} \times \hat{b} = 0.56102 \times 0.6355 = 0.3565$.

These calculations do not rely on strong causal assumptions. For a causal interpretation, we need a more rigorous framework.

### Two Approaches in Causal Mediation Analysis

The `mediation` package [@imai2010general, @imai2010identification] enables **causal mediation analysis**. It supports two inference types:

1.  **Model-Based Inference**

    -   Assumptions:

        -   Treatment is randomized (or approximated via matching).

        -   Sequential Ignorability: No unobserved confounding in:

            1.  Treatment → Mediator

            2.  Treatment → Outcome

            3.  Mediator → Outcome

        -   This assumption is hard to justify in observational studies.

2.  **Design-Based Inference**

    -   Relies on experimental design to isolate the causal mechanism.

**Notation**

We follow the standard potential outcomes framework:

-   $M_i(t)$ = mediator under treatment condition $t$

-   $T_i \in {0,1}$ = treatment assignment

-   $Y_i(t, m)$ = outcome under treatment $t$ and mediator value $m$

-   $X_i$ = observed pre-treatment covariates

The **treatment effect** for an individual $i$: $$
\tau_i = Y_i(1,M_i(1)) - Y_i (0,M_i(0))
$$ which decomposes into:

1.  **Causal Mediation Effect (ACME)**:

$$
\delta_i (t) = Y_i (t,M_i(1)) - Y_i(t,M_i(0))
$$

2.  **Direct Effect (ADE)**:

$$
\zeta_i (t) = Y_i (1, M_i(1)) - Y_i(0, M_i(0))
$$

Summing up:

$$
\tau_i = \delta_i (t) + \zeta_i (1-t)
$$

**Sequential Ignorability Assumption**

For CMA to be valid, we assume:

$$
\begin{aligned}
\{ Y_i (t', m), M_i (t) \} &\perp T_i |X_i = x\\
Y_i(t',m) &\perp M_i(t) | T_i = t, X_i = x
\end{aligned}
$$

-   First condition is the standard strong ignorability condition where treatment assignment is random conditional on pre-treatment confounders.

-   Second condition is stronger where the mediators is also random given the observed treatment and pre-treatment confounders. This condition is satisfied only when there is no unobserved pre-treatment confounders, and post-treatment confounders, and multiple mediators that are correlated.

**Key Challenge**

⚠️ **Sequential Ignorability is not testable.** Researchers should conduct **sensitivity analysis**.

We now fit a **causal mediation model** using `mediation`.

```{r}
library(mediation)
set.seed(2014)
data("framing", package = "mediation")

# Step 1: Fit mediator model (M ~ T, X)
med.fit <-
    lm(emo ~ treat + age + educ + gender + income, data = framing)

# Step 2: Fit outcome model (Y ~ M, T, X)
out.fit <-
    glm(
        cong_mesg ~ emo + treat + age + educ + gender + income,
        data = framing,
        family = binomial("probit")
    )

# Step 3: Causal Mediation Analysis (Quasi-Bayesian)
med.out <-
    mediate(
        med.fit,
        out.fit,
        treat = "treat",
        mediator = "emo",
        robustSE = TRUE,
        sims = 100
    )  # Use sims = 10000 in practice
summary(med.out)

```

**Alternative: Nonparametric Bootstrap**

```{r}
med.out <-
    mediate(
        med.fit,
        out.fit,
        boot = TRUE,
        treat = "treat",
        mediator = "emo",
        sims = 100,
        boot.ci.type = "bca"
    )
summary(med.out)

```

If we suspect **moderation**, we include an **interaction term**.

```{r}
med.fit <-
    lm(emo ~ treat + age + educ + gender + income, data = framing)

out.fit <-
    glm(
        cong_mesg ~ emo * treat + age + educ + gender + income,
        data = framing,
        family = binomial("probit")
    )

med.out <-
    mediate(
        med.fit,
        out.fit,
        treat = "treat",
        mediator = "emo",
        robustSE = TRUE,
        sims = 100
    )

summary(med.out)
test.TMint(med.out, conf.level = .95)  # Tests for interaction effect

```

Since **sequential ignorability** is untestable, we examine how unmeasured confounding affects ACME estimates.

```{r}
# Load required package
library(mediation)

# Simulate some example data
set.seed(123)
n <- 100
data <- data.frame(
  treat = rbinom(n, 1, 0.5),       # Binary treatment
  med = rnorm(n),                   # Continuous mediator
  outcome = rnorm(n)                 # Continuous outcome
)

# Fit the mediator model (med ~ treat)
med_model <- lm(med ~ treat, data = data)

# Fit the outcome model (outcome ~ treat + med)
outcome_model <- lm(outcome ~ treat + med, data = data)

# Perform mediation analysis
med_out <- mediate(med_model, 
                   outcome_model, 
                   treat = "treat", 
                   mediator = "med", 
                   sims = 100)

# Conduct sensitivity analysis
sens_out <- medsens(med_out, sims = 100)

# Print and plot results
summary(sens_out)
plot(sens_out)
```

-   If ACME confidence intervals contain 0, the effect is not robust to confounding.

Alternatively, using $R^2$ interpretation, we need to specify the direction of confounder that affects the mediator and outcome variables in `plot` using `sign.prod = "positive"` (i.e., same direction) or `sign.prod = "negative"` (i.e., opposite direction).

```{r, message=FALSE, eval = FALSE}
plot(sens.out, sens.par = "R2", r.type = "total", sign.prod = "positive")
```

**Summary: Causal Mediation vs. Traditional Mediation**

| **Aspect**                | **Traditional Mediation** | **Causal Mediation**         |
|-----------------------|-----------------------|--------------------------|
| **Model Assumption**      | Linear regressions        | Potential outcomes framework |
| **Assumptions Needed**    | No omitted confounders    | Sequential ignorability      |
| **Inference Method**      | Product of coefficients   | Counterfactual reasoning     |
| **Bootstrapping?**        | Common                    | Essential                    |
| **Sensitivity Analysis?** | Rarely used               | Strongly recommended         |