# Introduction to mediation analysis

In [Peters et al.](https://www.nature.com/articles/ncomms9570) where the authors developed a transcriptomics ageing clock, the authors hypothesized that chrononological age can be predicted by the expression level.

![](images/clipboard-4135855152.png)

The authors write,

> In addition, this is the first large-scale study testing the hypothesis that changes in gene expression with chronological age are epigenetically mediated by changes of methylation levels at specific loci.

In other words, the change in expression drive changes in methylation level which further drive the chronological age. This is a typical case of "mediation" analysis where methylation is a mediator that explains the chronological age.

![](images/clipboard-3338854842.png)

We will do a toy example trying to understand how to test if a variable has a "mediating" effect.

# Analyzing mediation effects

The goal of mediation analysis is **not to draw causality** but to show that possibly the independent variable influences the **mediator** variable that influences the final variable.

There are multiple steps involved in doing a mediation analysis. here are three steps of regression involved: $X \longrightarrow Y$ , $X \longrightarrow M$ and $X+M \longrightarrow Y$.

## Step 1

We try to estimate the total effect $Y = \gamma_1 + \tau X + \epsilon_1$

We will read the data first:

```{r}
data <- read.csv("mediation_data.csv")
head(data)
```

```{r}
library(ggplot2)
library(patchwork)
p1 <- ggplot(data, aes(x = X, y = Y)) +
  geom_point() +
  theme_minimal()
p1
```

```{r}
p2 <- ggplot(data, aes(x = X, y = M)) +
  geom_point() +
  theme_minimal()
p2
```

```{r}
p3 <- ggplot(data, aes(x = M, y = Y)) +
  geom_point() +
  theme_minimal()
p3
```

We can visualize them side by side
```{r, fig.width=12, fig.height=4}
p1 | p2 | p3
```

```{r}
model.1 <- lm(Y ~ X, data)
summary(model.1)
```

## Step 2

In step 2, we estimate the following relationship: $M = \gamma_2 + \alpha X + \epsilon_2$


```{r}
model.M <- lm(M ~ X, data)
summary(model.M)
```


## Step 3
In step 3, we estimate, if the mediator $M$ has an effect on $Y$ after controlling for $X$: $Y=\gamma_3 + \tau' X + \beta M + \epsilon_3$

```{r}
model.Y <- lm(Y ~ X + M, data)
summary(model.Y)
```


After the three steps, if the effect of $X$ on $Y$ is reduced after including $M$ in the model, then we can say that $M$ is a mediator. We can infer this by comparing the coefficients of $X$ in model 1 and model 3. To test the significance of the mediation effect, we can use the Sobel test or preferably the bootstrapping method.

We will use the `mediate()` method to do the bootstrapping method (over the Sobel method which is more conservative). `mediate()` takes two arguments, the two models $X \longrightarrow M$ and $X +M \longrightarrow Y$ along with the treatment variable and mediator variable.

```{r}
library(mediation)
results <- mediate(
  model.m = model.M,
  model.y = model.Y,
  treat = "X", mediator = "M",
  boot = TRUE, sims = 500
)
summary(results)
```

Total effect of X on Y is 0.3961 (ACME)

Direct effect of X on Y after taking account into M is 0.0396 (ADE), this is the coefficient $\tau'$ in the Step 3 model.

The mediation effect (ACME) is the total effect minus the direct effect (0.3961 - 0.0396 = 0.3565), which  is also equals to the product of coeffiecient $\alpha$ of $X$ in the second step and coefficient of $\beta$ of  $M$ in the step  or 0.56102 * 0.6355 = 0.3565.

We can also run Sobel's test using the multilevel package:

```{r}
library(multilevel)
fit <- sobel(pred = data$X, med = data$M, out = data$Y)
fit
```


or the bda package:
```{r}
library(bda)
fit <- mediation.test(mv = data$M, iv = data$Y, dv = data$X)
fit
```
