mediation.Rmd

# Mediation

## [Tutorial](https://ademos.people.uic.edu/Chapter14.html#2_mediation_analyses)

In this example we’ll say we are interested in whether the number of hours since dawn (X) affect the subjective ratings of wakefulness (Y) 100 graduate students through the consumption of coffee (M).

Note that we are intentionally creating a mediation effect here (because statistics is always more fun if we have something to find) and we do so below by creating M so that it is related to X and Y so that it is related to M. This creates the causal chain for our analysis to parse.
```{r}
set.seed(123) #Standardizes the numbers generated by rnorm; see Chapter 5
N <- 100 #Number of participants; graduate students
X <- rnorm(N, 175, 7) #IV; hours since dawn
M <- 0.7*X + rnorm(N, 0, 5) #Suspected mediator; coffee consumption 
Y <- 0.4*M + rnorm(N, 0, 5) #DV; wakefulness
(Meddata <- data.frame(X, M, Y))
```

This package (mediation) uses the more recent bootstrapping method of Preacher & Hayes (2004) to address the power limitations of the Sobel Test. This method computes the point estimate of the indirect effect (ab) over a large number of random sample (typically 1000) so it does not assume that the data are normally distributed and is especially more suitable for small sample sizes than the Barron & Kenny method.

To run the mediate function, we will again need a model of our IV (hours since dawn), predicting our mediator (coffee consumption) like our Path A model above. We will also need a model of the direct effect of our IV (hours since dawn) on our DV (wakefulness), when controlling for our mediator (coffee consumption). When can then use mediate to repeatedly simulate a comparison between these models and to test the significance of the indirect effect of coffee consumption.

```{r}
library(mediation)
?mediate
fitM <- lm(M ~ X,     data=Meddata) #IV on M; Hours since dawn predicting coffee consumption
fitY <- lm(Y ~ X + M, data=Meddata) #IV and M on DV; Hours since dawn and coffee predicting wakefulness
# gvlma(fitM) #data is positively skewed; could log transform (see Chap. 10 on assumptions)
```

```{r}
fitMed <- mediate(fitM, fitY, treat="X", mediator="M")
summary(fitMed)
plot(fitMed)
```

```{r}
#Bootstrap
fitMedBoot <- mediate(fitM, fitY, boot=TRUE, sims=999, treat="X", mediator="M")
summary(fitMedBoot)
```

## Actual analysis for our paper
```{r}
library(tidyverse)
library(arm)
empa <- read_csv("data/cortical_thickness_base.csv")

empa <- empa %>% 
  mutate(
    age_s = rescale(age),
    fs_s = rescale(IRI.FS),
    pt_s = rescale(IRI.PT),
    supre_s = rescale(ERQ.Supre),
    sex = factor(sex),
    group = factor(group)
  )

glimpse(empa)

temp <- lm(r015 ~ IRI.FS + IRI.PT + IRI.EC + IRI.PD + ERQ.Reev + ERQ.Supre + age + sex, data = empa)
temp <- lm(r015 ~ group + IRI.FS + ERQ.Supre + IRI.PT, data = empa)
temp <- lm(r015 ~ IRI.FS + ERQ.Supre + IRI.PT + IRI.EC + age + sex, data = empa)
temp <- lm(r015 ~ IRI.FS + group + age + sex, data = empa)
temp <- lm(r015 ~ IRI.PT + group + age + sex, data = empa)
temp <- lm(r015 ~ supre_s + fs_s + pt_s + group + age + sex, data = empa) # si
# reportar con apatables, editable a word + path graph
temp <- lm(r015 ~ supre_s + fs_s + pt_s + group + age + sex + supre_s:group, data = empa) # si

temp <- lm(r015 ~ ERQ.Supre + age + sex, data = empa) # no
temp <- lm(r015 ~ ERQ.Supre + IRI.PT + IRI.FS + I(ERQ.Supre*IRI.FS) + I(IRI.FS*IRI.PT) + I(IRI.PT*ERQ.Supre) + group + age + sex, data = empa) 
summary(temp)

fitA <- lm(r015 ~ group, data = empa)
fitB <- lm(r015 ~ group + IRI.FS, data = empa)
fitC <- mediate(fitA, fitB, treat = "group", mediator = "IRI.FS")
fitC <- mediate(fitA, fitB, boot = TRUE, sims = 999, treat = "group", mediator = "IRI.FS")
summary(fitC)
plot(fitC)

fitA <- lm(r015 ~ group, data = empa)
fitB <- lm(r015 ~ group + IRI.EC, data = empa)
fitC <- mediate(fitA, fitB, treat = "group", mediator = "IRI.EC", boot = TRUE, sims = 999)
summary(fitC)
plot(fitC)
```


```{r}
library(lavaan)
library(semPlot)
```


```{r}
model <- 'r015 ~ b1 * IRI.FS + b2 * IRI.EC + c * group
IRI.FS ~ a1 * group
IRI.EC ~ a2 * group
indirect1 := a1 * b1
indirect2 := a2 * b2
total := c + (a1 * b1) + (a2 * b2)
IRI.FS ~~ IRI.EC' # probs OK!
```

```{r}
fitL <- sem(model, data = empa)
summary(fitL)

semPaths(fitL, 'std', layout = 'circle')
```

###  Only FS
```{r}
library(lavaan)
library(semPlot)
```

```{r}
model2 <- 'r015 ~ b * IRI.FS + c * group
IRI.FS ~ a * group
indirect := a * b
total := c + (a * b)'
```

```{r}
fitL2 <- sem(model2, data = empa)
summary(fitL2)

semPaths(fitL2, 'std', layout = 'circle')
```

### FS and Supresion
```{r}
model3 <- 'r015 ~ b1 * IRI.FS + b2 * ERQ.Supre + c * group
IRI.FS ~ a1 * group
ERQ.Supre ~ a2 * group
indirect1 := a1 * b1
indirect2 := a2 * b2
total := c + (a1 * b1) + (a2 * b2)
IRI.FS ~~ ERQ.Supre' 
```

### FS and supresion and PT
```{r}
model3 <- 'r015 ~ b1 * IRI.FS + b2 * ERQ.Supre + b3 * IRI.PT + c * group
IRI.FS ~ a1 * group
ERQ.Supre ~ a2 * group
IRI.PT ~ a3 * group
indirect1 := a1 * b1
indirect2 := a2 * b2
indirect3 := a3*b3
total := c + (a1 * b1) + (a2 * b2) + (a3*b3)
IRI.FS ~~ ERQ.Supre'

```

```{r}
model3 <- 'r015 ~ b * supre_s + c * group
supre_s ~ a * group
indirect := a * b
total := c + (a * b)' 
```


### FS and supresion and PT
```{r}
model3 <- 'r015 ~ IRI.FS + ERQ.Supre + IRI.PT'

```

```{r}
fitL3 <- sem(model3, data = empa)
summary(fitL3)

semPaths(fitL3, 'std', layout = 'spring')
```
```{r}
library(mediation)
fitA <- lm(r015 ~ group, data = empa)
fitB <- lm(r015 ~ group + supre_s, data = empa)
fitC <- mediate(fitA, fitB, treat = "group", mediator = "supre_s")
fitC <- mediate(fitA, fitB, boot = TRUE, sims = 999, treat = "group", mediator = "supre_s")
summary(fitC)
plot(fitC)
```