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 ="")
} 

Download Code

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.

Category:
Mediation Analyses
Tags:
, ,

Join the conversation! 1 Comment

Leave a Reply

Please log in using one of these methods to post your comment:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out /  Change )

Google+ photo

You are commenting using your Google+ account. Log Out /  Change )

Twitter picture

You are commenting using your Twitter account. Log Out /  Change )

Facebook photo

You are commenting using your Facebook account. Log Out /  Change )

w

Connecting to %s

%d bloggers like this: