# Bootstrap Mediation Analyses

Below is the R code for conducting mediation analyses using the bootstrap method detailed in Preacher and Hayes (2004). Please note that my function uses the bias corrected and accelarated (BCa) confidence intervals, whereas their SPSS and SAS macros use the bias corrected (BC) confidence intervals. However, they should yield similar results.

The `med()` function takes the following arguments:

• `iv` = independent variable (X). Use: `iv = "iv`
• `dv` = dependent variable (Y). Use: `dv = "dv"`
• `mv` = mediating variable (M). Use: `mv = "mv"`
• `cv` = covariate(s) if any. Use: `cv = "cv.1"` for single covariate or `cv = c("cv.1", "cv.2", ..., "cv.n")` for multiple covariates.
• `dec` = number of decimal places (default = 4)
• `rep` = number of bootstrap replications (default = 1000)

## Function

```require(boot)
require(car)
require(MASS)
require(QuantPsyc) # lm.beta() function

#---- Bootstrap Mediation Function ----

boot.med <- function(data, formula.m, formula.y, iv, mv, i) {
data.resamp <- data[i,]
reg.m <- lm(formula.m, data=data.resamp)
reg.y <- lm(formula.y, data=data.resamp)
# Make sure that you pull the correct coefficients from each regression equation
# With interactions, quadratics, etc, R may reorder the equation so you need to check it.
a <- coefficients(reg.m)[iv] # the a path coefficient
b <- coefficients(reg.y)[mv] # the b path coefficient
ind <- a*b # the indirect effect (product of a*b)
dir <- coefficients(reg.y)[iv] # the direct effect (c')
total <- ind + dir # teturns total effect
boot.med <- return(c(ind, dir, total))
}

#---- Mediation Analyses ----
#
# data = data frame
# iv = independent variable (X)
# dv = dependent variable (Y)
# mv = mediating variable (M)
# cv = covariate(s) if any. Use: cv = c("var.1", "var.2", "var.3", ..., "var.n")
# reps = number of bootstrap replications (default = 1000)
#

med <- function(data, iv, dv, mv, cv, dec = 4, reps = 1000) {
# Formulas for Regressions
if(missing(cv)){
iv.formula <- as.formula(paste(mv, paste("~"), paste(iv)))
mv.formula <- as.formula(paste(dv, paste("~"), paste(mv)))
dv.formula <- as.formula(paste(dv, paste("~"), paste(iv)))
dvmv.formula <- as.formula(paste(dv, paste("~"), paste(iv), paste("+"), paste(mv)))
} else {
iv.formula <- as.formula(paste(mv, paste("~"), paste(iv), paste("+"), paste(cv, collapse="+")))
mv.formula <- as.formula(paste(dv, paste("~"), paste(mv), paste("+"), paste(cv, collapse="+")))
dv.formula <- as.formula(paste(dv, paste("~"), paste(iv), paste("+"), paste(cv, collapse="+")))
dvmv.formula <- as.formula(paste(dv, paste("~"), paste(iv), paste("+"), paste(mv), paste("+"), paste(cv, collapse="+")))
}
## Regressions

# IV on MV
iv.out <- lm(iv.formula, data = data)
iv.coef <- coef(iv.out) # Coefficients
iv.se <- summary(iv.out)\$coefficients[,2] # SEs
iv.pr <- 2 * pt(abs(iv.coef/iv.se), df=nrow(iv.out\$model)-2, lower.tail = FALSE) # p-values
iv.ll <- iv.coef - 1.96 * iv.se # lower 95% CIs
iv.ul <- iv.coef + 1.96 * iv.se # upper 95% CIs
# Standardised Betas
iv.beta <- lm.beta(iv.out)
iv.beta <- cbind(NA, iv.beta)
iv.beta <- rbind(NA, iv.beta)[,2]
iv.rsquared <- summary(iv.out)\$r.squared #R2
iv.rsquared.pr <- pf(summary(iv.out)\$fstatistic[1], summary(iv.out)\$fstatistic[2],
summary(iv.out)\$fstatistic[3], lower.tail = F) # p-value for R2
# Output Table
iv.table <- cbind("Estimate" = round(iv.coef, dec),
"SE" = round(iv.se, dec),
"Pr(>|z|)" = round(2 * pt(abs(iv.coef/iv.se), df=nrow(iv.out\$model)-2,
lower.tail = FALSE), dec),
"LL" = round(iv.coef - 1.96 * iv.se, dec),
"UL" = round(iv.coef + 1.96 * iv.se, dec),
"Beta" = round(iv.beta, dec))

# MV on DV
mv.out <- lm(mv.formula, data = data)
mv.coef <- coef(mv.out) # Coefficients
mv.se <- summary(mv.out)\$coefficients[,2] # SEs
mv.pr <- 2 * pt(abs(mv.coef/mv.se), df=nrow(mv.out\$model)-2, lower.tail = FALSE) # p-values
mv.ll <- mv.coef - 1.96 * mv.se # lower 95% CIs
mv.ul <- mv.coef + 1.96 * mv.se # upper 95% CIs
# Standardised Betas
mv.beta <- lm.beta(mv.out)
mv.beta <- cbind(NA, mv.beta)
mv.beta <- rbind(NA, mv.beta)[,2]
mv.rsquared <- summary(mv.out)\$r.squared #R2
mv.rsquared.pr <- pf(summary(mv.out)\$fstatistic[1], summary(mv.out)\$fstatistic[2],
summary(mv.out)\$fstatistic[3], lower.tail = F) # p-value for R2
# Output Table
mv.table <- cbind("Estimate" = round(mv.coef, dec),
"SE" = round(mv.se, dec),
"Pr(>|z|)" = round(2 * pt(abs(mv.coef/mv.se), df=nrow(mv.out\$model)-2,
lower.tail = FALSE), dec),
"LL" = round(mv.coef - 1.96 * mv.se, dec),
"UL" = round(mv.coef + 1.96 * mv.se, dec),
"Beta" = round(mv.beta, dec))

# IV on DV
dv.out <- lm(dv.formula, data = data)
dv.coef <- coef(dv.out) # Coefficients
dv.se <- summary(dv.out)\$coefficients[,2] # SEs
dv.pr <- 2 * pt(abs(dv.coef/dv.se), df=nrow(dv.out\$model)-2, lower.tail = FALSE) # p-values
dv.ll <- dv.coef - 1.96 * dv.se # lower 95% CIs
dv.ul <- dv.coef + 1.96 * dv.se # upper 95% CIs
# Standardised Betas
dv.beta <- lm.beta(dv.out)
dv.beta <- cbind(NA, dv.beta)
dv.beta <- rbind(NA, dv.beta)[,2]
dv.rsquared <- summary(dv.out)\$r.squared #R2
dv.rsquared.pr <- pf(summary(dv.out)\$fstatistic[1], summary(dv.out)\$fstatistic[2],
summary(dv.out)\$fstatistic[3], lower.tail = F) # p-value for R2
# Output Table
dv.table <- cbind("Estimate" = round(dv.coef, dec),
"SE" = round(dv.se, dec),
"Pr(>|z|)" = round(2 * pt(abs(dv.coef/dv.se), df=nrow(dv.out\$model)-2,
lower.tail = FALSE), dec),
"LL" = round(dv.coef - 1.96 * dv.se, dec),
"UL" = round(dv.coef + 1.96 * dv.se, dec),
"Beta" = round(dv.beta, dec))

# IV and MV on DV
dvmv.out <- lm(dvmv.formula, data = data)
dvmv.coef <- coef(dvmv.out)
dvmv.se <- summary(dvmv.out)\$coefficients[,2]
dvmv.pr <- 2 * pt(abs(dvmv.coef/dvmv.se), df=nrow(dvmv.out\$model)-2, lower.tail = FALSE)
dvmv.ll <- dvmv.coef - 1.96 * dvmv.se
dvmv.ul <- dvmv.coef + 1.96 * dvmv.se
# Standardised Betas
dvmv.beta <- lm.beta(dvmv.out)
dvmv.beta <- cbind(NA, dvmv.beta)
dvmv.beta <- rbind(NA, dvmv.beta)[,2]
dvmv.rsquared <- summary(dvmv.out)\$r.squared
dvmv.rsquared.pr <- pf(summary(dvmv.out)\$fstatistic[1], summary(dvmv.out)\$fstatistic[2],
summary(dvmv.out)\$fstatistic[3], lower.tail = F)
# Output Table
dvmv.table <- cbind("Estimate" = round(dvmv.coef, dec),
"SE" = round(dvmv.se, dec),
"Pr(>|z|)" = round(2 * pt(abs(dvmv.coef/dvmv.se), df=nrow(data)-2,
lower.tail = FALSE), dec),
"LL" = round(dvmv.coef - 1.96 * dvmv.se, dec),
"UL" = round(dvmv.coef + 1.96 * dvmv.se, dec),
"Beta" = round(dvmv.beta, dec))
## Mediation Bootstrap
medboot.out <- boot(statistic = boot.med,
data = data,
formula.m = iv.formula,
formula.y = dvmv.formula,
iv = iv,
mv = mv,
R = reps,
parallel="multicore")
indirect.b <- round(summary(medboot.out)\$original[1], dec) # Indirect B
indirect.se <- round(summary(medboot.out)\$bootSE[1], dec) # Indirect SE
indirect.bca <- boot.ci(medboot.out, type="bca", index = 1) # Indirect BCA
indirect.bca.l <- round(indirect.bca\$bca[4], dec) # Extract lower BCa CI
indirect.bca.u <- round(indirect.bca\$bca[5], dec) # Extract upper BCa CI
direct.b <- round(summary(medboot.out)\$original[2], dec) # Direct B
direct.se <- round(summary(medboot.out)\$bootSE[2], dec) # Direct SE
direct.bca <- boot.ci(medboot.out, type="bca", index = 2) # indirect
direct.bca.l <- round(direct.bca\$bca[4], dec) # Extract lower BCa CI
direct.bca.u <- round(direct.bca\$bca[5], dec) # Extract upper BCa CI
total.b <- round(summary(medboot.out)\$original[3], dec) # Direct B
total.se <- round(summary(medboot.out)\$bootSE[3], dec) # Direct SE
total.bca <- boot.ci(medboot.out, type="bca", index = 3) # Total Effect
total.bca.l <- round(total.bca\$bca[4], dec) # Extract lower BCa CI
total.bca.u <- round(total.bca\$bca[5], dec) # Extract upper BCa CI
## Output

cat(paste("\n",
"Independent Variable: ", iv,
"\n", "Dependent Variable: ", dv,
"\n", "Mediator Variable: ", mv,
"\n",
"\n",
"IV (and CVs) on MV: Path a",
"\n",
sep = ""))

print(iv.table)

cat(paste("\n",
"R-Squared: ", round(iv.rsquared, dec), ", p-value: ", round(iv.rsquared.pr, dec),
"\n",
sep=""))

cat(paste("\n",
"\n",
"MV (and CVs) on DV: Path b",
"\n",
sep = ""))
print(mv.table)

cat(paste("\n",
"R-Squared: ", round(mv.rsquared, dec), ", p-value: ", round(mv.rsquared.pr, dec),
"\n",
sep=""))

cat(paste("\n",
"\n",
"IV (and CVs) on DV: Path c",
"\n",
sep = ""))
print(dv.table)
cat(paste("\n",
"R-Squared: ", round(dv.rsquared, dec), ", p-value: ", round(dv.rsquared.pr, dec),
sep=""))

cat(paste("\n",
"\n",
"IV and MV (and CVs) on DV",
"\n",
sep = ""))
print(dvmv.table)
cat(paste("\n",
"R-Squared: ", round(dvmv.rsquared, dec), ", p-value: ", round(dvmv.rsquared.pr, dec),
"\n",
sep=""))
cat("\n",
"============================================", "\n",
"Bootstrap Mediation Effects with 95% BCa CIs", "\n",
"============================================", "\n",
"\n",
"Indirect effect", "\n",
"---------------", "\n",
paste("B = ", indirect.b, ", SE = ", indirect.se, ", 95% CI [", indirect.bca.l, ", ", indirect.bca.u, "]", sep=""),
"\n", "\n",
"Direct effect", "\n",
"-------------", "\n",
paste("B = ", direct.b, ", SE = ", direct.se, ", 95% CI [", direct.bca.l, ", ", direct.bca.u, "]", sep=""),
"\n", "\n",
"Total effect", "\n",
"-------------", "\n",
paste("B = ", total.b, ", SE = ", total.se, ", 95% CI [", total.bca.l, ", ", total.bca.u, "]", sep=""),
"\n",
sep ="")
}
```

## Example

Assessing whether body image inflexibility (`biaaq.total`) mediates the relationship between negative body appearance evaluation (`mbsrq.ae`) and avoidance behaviours (`bicsi.avoid`) while controlling for age (`age`), gender (`gender`), and body image investment (`asir.comp`).

```>med(data = cbm.bdd, dv = "bicsi.avoid", iv = "mbsrq.ae", mv = "biaaq.total", cv = c("gender", "asir.comp"), dec = 3, reps = 1000)

Independent Variable: mbsrq.ae
Dependent Variable: bicsi.avoid
Mediator Variable: biaaq.total

IV (and CVs) on MV
Estimate    SE Pr(>|z|)      LL      UL   Beta
(Intercept)   79.286 6.112        0  67.306  91.265     NA
mbsrq.ae       6.926 0.890        0   5.181   8.671  0.394
genderMale    -6.301 1.774        0  -9.779  -2.824 -0.171
asir.comp    -12.881 1.277        0 -15.384 -10.378 -0.518

R-Squared: 0.547, p-value: 0

MV (and CVs) on DV
Estimate    SE Pr(>|z|)     LL     UL   Beta
(Intercept)    2.761 0.268    0.000  2.236  3.286     NA
biaaq.total   -0.021 0.002    0.000 -0.025 -0.017 -0.737
genderMale    -0.272 0.059    0.000 -0.388 -0.157 -0.256
asir.comp     -0.173 0.051    0.001 -0.273 -0.073 -0.241

R-Squared: 0.417, p-value: 0

IV (and CVs) on DV
Estimate    SE Pr(>|z|)     LL     UL   Beta
(Intercept)    1.230 0.233    0.000  0.773  1.688     NA
mbsrq.ae      -0.181 0.034    0.000 -0.248 -0.114 -0.357
genderMale    -0.133 0.068    0.051 -0.266  0.000 -0.125
asir.comp      0.084 0.049    0.086 -0.011  0.180  0.117

R-Squared: 0.208, p-value: 0

IV and MV (and CVs) on DV
Estimate    SE Pr(>|z|)     LL     UL   Beta
(Intercept)    2.806 0.269    0.000  2.278  3.334     NA
mbsrq.ae      -0.043 0.033    0.191 -0.108  0.022 -0.086
biaaq.total   -0.020 0.002    0.000 -0.024 -0.015 -0.688
genderMale    -0.258 0.060    0.000 -0.375 -0.141 -0.243
asir.comp     -0.172 0.051    0.001 -0.272 -0.072 -0.239

R-Squared: 0.422, p-value: 0

============================================
Bootstrap Mediation Effects with 95% BCa CIs
============================================

Indirect effect
---------------
B = -0.138, SE = 0.024, 95% CI [-0.194, -0.098]

Direct effect
-------------
B = -0.043, SE = 0.038, 95% CI [-0.126, 0.025]

Total effect
-------------
B = -0.181, SE = 0.038, 95% CI [-0.259, -0.108]
```

## Write up of results in APA format

Controlling for the gender, age, and body image investment, the mediation analysis showed that appearance evaluation had a significant total effect on appearance-fixing (B = -0.18, 95% CI [-0.26, -0.11], p < .05), a non-significant direct effect (B = -0.04, 95% CI [-0.13, 0.03], p = n.s.), and a significant indirect effect (B = -0.14, 95% CI [-0.19, -0.10], p < .05). This pattern of results suggests that body image inflexibility fully mediates the relationship between appearance evaluation and avoidance behaviours.

## Reference

Preacher, K.J., & Hayes, A.F. (2004). SPSS and SAS procedures for estimating indirect effects in simple Mediation models. Behavior Research Methods, Instruments, & Computers, 36, 717-731.

## One thought on “Bootstrap Mediation Analyses”

This site uses Akismet to reduce spam. Learn how your comment data is processed.