This tutorial shows you how to do mediation analysis in R from scratch. As we know, quantitatively, mediation analysis is about if the product of a*b is statistically significant (see figure below).
If you assume that a*b is normally distributed, you can just directly test the level of significance. If you do not assume normal distribution, you could use quantile or other methods.
Mediation Analysis in R from Scratch
The following is the key, complete R code for the mediation analysis.
library(boot)
set.seed(123)
Mediation_function<-function(data_used,i)
{
# Sample a data
data_temp=data_used[i,]
# a path
result_a<-lm(M~X, data = data_temp)
a<-result_a$coefficients[2]
# b path
result_b<-lm(Y~M+X, data = data_temp)
b<-result_b$coefficients[2]
#calculating the indirect effect
indirect_effect<-a*b
return(indirect_effect)
}
# use boot() to do bootstrapping mediation analysis
boot_mediation <- boot(Mediation_data, Mediation_function, R=5000)
# read data
Mediation_data <-read.csv('https://raw.githubusercontent.com/TidyPython/Mediation_analysis/main/mediation_hypothetical_data.csv')
# show the head of the data
head(Mediation_data)
# print out confidence intervals
boot.ci(boot.out = boot_mediation, type = c("norm", "perc"))
The following is the output. As you can see, it produces the confidence intervals based on normal distribution assumptions and percentile as well. Neither of the confidence intervals includes zero. Thus, we can conclude that the mediation model works. (It should be noted that, the normal distribution assumption is for the 5000 bootstrapping indirect effects.)
BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
Based on 5000 bootstrap replicates
CALL :
boot.ci(boot.out = boot_mediation, type = c("norm", "perc"))
Intervals :
Level Normal Percentile
95% ( 0.0377, 0.1570 ) ( 0.0379, 0.1559 )
Calculations and Intervals on Original Scale
We can also plot the 5000 indirect effects. As we can see, it is pretty normally distributed.
# plot the 5000 indirect effects
plot(boot_mediation)
SE and SD in Bootstrapping Mediation Analysis
Note that, standard error (SE) of a statistic is the standard deviation (SD) of its sampling distribution or an estimate of that standard deviation (based on Wikipedia).
Thus, we can calculate the SE of the indirect effects by calculating the SD of the bootstrapped sample. The following is the R code to show how.
# calculate the standard deviation
sd(boot_mediation$t)
The following is the output:
## [1] 0.03044034
We can double check by printing out the the CI again.
# print out the bootstrapping output again
boot_mediation
The following is the output.
ORDINARY NONPARAMETRIC BOOTSTRAP
Call:
boot(data = Mediation_data, statistic = Mediation_function, R = 5000)
Bootstrap Statistics :
original bias std. error
t1* 0.09748237 0.0001266305 0.03044034
Thus, we can see that SE of the indirect effect is the same as the sd of the bootstrapped sample of indirect effect.
Understand “bias” in the output
The difference between the mean of the bootstrap estimates and the original sample estimate is the bias.
# the mean of 5000 indirect effect
print(mean(boot_mediation$t))
The following is the output:
[1] 0.097609
From the original sample, the following is the output.
#the original sample estimate
print(boot_mediation$t0)
The following is the output:
X 0.09748237
The difference between these two numbers is the bias.
# the difference = the bias
bias=mean(boot_mediation$t)-boot_mediation$t0
print(bias)
The output:
> print(bias)
X
0.0001266305
Bias used in Confidence Interval
We can apply this bias in the calculation of CI based on normal distribution assumption.
# calculate the confidence interval from scratch
print(boot_mediation$t0-bias + c(-1, 1) * 1.96 * sd(boot_mediation$t))
The following is the output.
[1] 0.03769268 0.15701880
We can print again the CI from the boot() function.
# print out confidence interval based on normal distribution
boot.ci(boot.out = boot_mediation, type = c("norm"))
The following is the output:
BOOTSTRAP CONFIDENCE INTERVAL CALCULATIONS
Based on 5000 bootstrap replicates
CALL :
boot.ci(boot.out = boot_mediation, type = c("norm"))
Intervals :
Level Normal
95% ( 0.0377, 0.1570 )
Calculations and Intervals on Original Scale
We can see that (0.03769268 0.15701880) is the same as ( 0.0377, 0.1570 ), after some rounding up. That means, we need to apply bias when calculate CI based on normal distribution assumption.