Testing indirect effects/mediation in R

What is mediation?

There are many ways to define mediation and mediators. Here’s one way: Mediation is the process by which one variable transmits an effect onto another through one or more mediating variables. For example, as room temperature increases, people get thirstier, and then they drink more water. In this case, thirst transmits the effect of room temperature on water drinking.

What is an indirect effect?

The indirect effect quantifies a mediation effect, if such an effect exists. Referring to the thirst example above, in statistical terms, the indirect effect quantifies the extent to which room temperature is associated with water drinking indirectly through thirstiness. If you’re familiar with interpreting regression coefficients and the idea of controlling for other variables, then you might find it intuitive to think of the indirect effect as the decrease in the relationship between room temperature and water drinking after you’ve partialed out the association between room temperature and thirtiness. In other words, how much does the coefficient for room temperature decrease when you control for thirstiness?

Model and Conceptual Assumptions

Correct functional form. Your model variables share linear relationships and don’t interact with eachother.

No omitted influences. This one is hard: Your model accounts for all relevant influences on the variables included. All models are wrong, but how wrong is yours?

Accurate measurement. Your measurements are valid and reliable. Note that unreliable measures can’t be valid, and reliable measures don’t necessairly measure just one construct or even your construct.

Well-behaved residuals. Residuals (i.e., prediction errors) aren’t correlated with predictor variables or eachother, and residuals have constant variance across values of your predictor variables. Also, residual error terms aren’t correlated across regression equations. This could happen if, for example, some omitted variable causes both thirst and water drinking.

Libraries

library(tidyverse)
library(knitr)
library(lavaan)
library(psych)
library(MBESS)

Data

I took the data from Table 3.1 in Mackinnon (2008, p. 56) [.csv]

thirst_dat <- "data/mackinnon_2008_t3.1.csv" %>% read_csv()

Print first and last five observations

thirst_dat %>% 
  headTail() %>% 
  kable()

id	room_temp	thirst	consume
1	70	4	3
2	71	4	3
3	69	1	3
4	70	1	3
…	…	…	…
47	71	4	4
48	71	3	2
49	72	4	5
50	70	2	2

Visualize relationships

It’s always a good idea to look at your data. Check some assumptions.

thirst_dat %>% 
  select(room_temp, thirst, consume) %>% 
  pairs.panels()

Conceptual Diagram

Figure 1. Room temperature is associated with water drinking indirectly through thirstiness

Write model to test indirect effect using `sem()` from lavaan

~ = Regress onto …

Within the regression models, I label coefficients with the astrix.

:= = Define a new parameter. Note when you define new parameter with :=, you can use the astrix to multiply values

For more details about lavaan syntax, see the tutorials tab at the lavaan website (linked in Resources below)

mod1 <- "# a path
         thirst ~ a * room_temp

         # b path
         consume ~ b * thirst

         # c prime path 
         consume ~ cp * room_temp

         # indirect and total effects
         ab := a * b
         total := cp + ab"

Set random seed so results can be reproduced

set.seed(1234)

Fit model

You must specify bootstrapping in the sem() function

fsem1 <- sem(mod1, data = thirst_dat, se = "bootstrap", bootstrap = 10000)

Summarize model

standardized = TRUE adds standardized estimate to the model output. Also, see help("standardizedsolution")

summary(fsem1, standardized = TRUE)

## lavaan (0.5-23.1097) converged normally after  11 iterations
## 
##   Number of observations                            50
## 
##   Estimator                                         ML
##   Minimum Function Test Statistic                0.000
##   Degrees of freedom                                 0
##   Minimum Function Value               0.0000000000000
## 
## Parameter Estimates:
## 
##   Information                                 Observed
##   Standard Errors                            Bootstrap
##   Number of requested bootstrap draws            10000
##   Number of successful bootstrap draws           10000
## 
## Regressions:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##   thirst ~                                                              
##     room_temp  (a)    0.339    0.101    3.362    0.001    0.339    0.371
##   consume ~                                                             
##     thirst     (b)    0.451    0.149    3.025    0.002    0.451    0.413
##     room_temp (cp)    0.208    0.130    1.599    0.110    0.208    0.208
## 
## Variances:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##    .thirst            0.911    0.153    5.966    0.000    0.911    0.862
##    .consume           0.912    0.155    5.892    0.000    0.912    0.723
## 
## Defined Parameters:
##                    Estimate  Std.Err  z-value  P(>|z|)   Std.lv  Std.all
##     ab                0.153    0.064    2.396    0.017    0.153    0.153
##     total             0.360    0.116    3.098    0.002    0.360    0.361

Print all model parameters

In the boot.ci.type argument, I ask for bia-corrected and accelerated confidence intervals.

parameterestimates(fsem1, boot.ci.type = "bca.simple", standardized = TRUE) %>% 
  kable()

lhs	op	rhs	label	est	se	z	pvalue	ci.lower	ci.upper	std.lv	std.all	std.nox
thirst	~	room_temp	a	0.3385926	0.1006968	3.362496	0.0007724	0.1303043	0.5290935	0.3385926	0.3708979	0.3294301
consume	~	thirst	b	0.4510391	0.1491257	3.024556	0.0024900	0.1633510	0.7466368	0.4510391	0.4126006	0.4126006
consume	~	room_temp	cp	0.2076475	0.1298255	1.599436	0.1097238	-0.0568067	0.4465309	0.2076475	0.2080748	0.1848112
thirst	~~	thirst		0.9110761	0.1527217	5.965598	0.0000000	0.6618918	1.2978137	0.9110761	0.8624347	0.8624347
consume	~~	consume		0.9124389	0.1548648	5.891840	0.0000000	0.6772054	1.3154000	0.9124389	0.7227811	0.7227811
room_temp	~~	room_temp		1.2676000	0.0000000	NA	NA	1.2676000	1.2676000	1.2676000	1.0000000	1.2676000
ab	:=	a*b	ab	0.1527185	0.0637271	2.396444	0.0165550	0.0557988	0.3153882	0.1527185	0.1530327	0.1359231
total	:=	cp+ab	total	0.3603660	0.1163348	3.097664	0.0019505	0.1186142	0.5777528	0.3603660	0.3611075	0.3207343

Interpretation

Every 1&degF increase in room temperature was associated with an a = 0.339 (S.E. = 0.101) increase in thirstiness units. Adjusting for room temperature, every 1-unit increase in thirstiness was associated with drinking b = 0.451 (S.E. = 0.149) more deciliters of water. Increases in room temperature were associated with increases in water drinking indirectly through increases in thirstiness. Specifically, for every a = 0.339 unit increase in the association between room temperature and thirstiness, there was an ab = 0.153 (S.E. = 0.064) increase in deciliters of water people drank. Importatnly, a bias-corrected bootstrapped confidence interval with 10,000 samples was above zero, 95% CI [0.06, 0.32]. Last, there was no sufficient evidence that room temperature was associated with how many deciliters of water people drank independent of its association with thirstiness, c’ = 0.208 (S.E. = 0.130).

Test same model using `mediation()` from MBESS

The syntax for mediation() doesn’t have as steep a learning curve as lavaan, but lavaan (and SEM in general) has a gazillion-fold more flexability in specifying more involved models.

with(thirst_dat, mediation(x = room_temp, mediator = thirst, dv = consume, bootstrap = TRUE, which.boot = "BCa", B = 10000))

## [1] "Bootstrap resampling has begun. This process may take a considerable amount of time if the number of replications is large, which is optimal for the bootstrap procedure."

##                                          Estimate CI.Lower_BCa
## Indirect.Effect                        0.15271850  0.054774249
## Indirect.Effect.Partially.Standardized 0.13455697  0.047500499
## Index.of.Mediation                     0.15303271  0.056483230
## R2_4.5                                 0.09305941  0.017104506
## R2_4.6                                 0.02322601  0.003224223
## R2_4.7                                 0.08378220  0.018697373
## Ratio.of.Indirect.to.Total.Effect      0.42378715  0.114709549
## Ratio.of.Indirect.to.Direct.Effect     0.73546979 -0.139295208
## Success.of.Surrogate.Endpoint          1.06430568  0.356034288
## Residual.Based_Gamma                   0.05094366 -0.012101027
## Residual.Based.Standardized_gamma      0.04865370 -0.015292266
## SOS                                    0.71365334  0.250959300
##                                        CI.Upper_BCa
## Indirect.Effect                          0.31703756
## Indirect.Effect.Partially.Standardized   0.26609057
## Index.of.Mediation                       0.32588089
## R2_4.5                                   0.24214258
## R2_4.6                                   0.08866451
## R2_4.7                                   0.23564896
## Ratio.of.Indirect.to.Total.Effect        1.46825891
## Ratio.of.Indirect.to.Direct.Effect      17.49103913
## Success.of.Surrogate.Endpoint            2.63708443
## Residual.Based_Gamma                     0.16891177
## Residual.Based.Standardized_gamma        0.16425113
## SOS                                      0.99749314

Plot the mediation effect

with(thirst_dat, mediation.effect.plot(x = room_temp, mediator = thirst, dv = consume, ylab = "Water Drank (dl)", xlab = "Thirstiness (1/5 = Not at all thirty/Very thirsty)"))

The plot above depicts the relationship between the proposed mediator (thirstiness) and outcome variable (water drank, in dl) at different levels of the proposed antecedent (room temperature, in &degF). The plot doesn’t label this, but if check out the right triangle formed in between the vertical lines marking the a coefficient, you’ll see the indirect effect, which is the height of this triangle.

Test the same model using `mediate()` from psych

mediate(consume ~ room_temp + thirst, data = thirst_dat, n.iter = 10000) %>% print(short = FALSE)

## Call: mediate(y = consume ~ room_temp + thirst, data = thirst_dat, 
##     n.iter = 10000)
## 
## The DV (Y) was  consume . The IV (X) was  room_temp . The mediating variable(s) =  thirst .
## 
## Total Direct effect(c) of  room_temp  on  consume  =  0.36   S.E. =  0.13  t direct =  2.68   with probability =  0.01
## Direct effect (c') of  room_temp  on  consume  removing  thirst  =  0.21   S.E. =  0.13  t direct =  1.56   with probability =  0.13
## Indirect effect (ab) of  room_temp  on  consume  through  thirst   =  0.15 
## Mean bootstrapped indirect effect =  0.15  with standard error =  0.06  Lower CI =  0.04    Upper CI =  0.29
## R2 of model =  0.28
## 
##  Full output  
## 
##  Total effect estimates (c) 
##           consume   se    t    Prob
## room_temp    0.36 0.13 2.68 0.00998
## 
## Direct effect estimates     (c') 
##           consume   se    t    Prob
## room_temp    0.21 0.13 1.56 0.12600
## thirst       0.45 0.15 3.09 0.00336
## 
##  'a'  effect estimates 
##        room_temp   se    t    Prob
## thirst      0.34 0.12 2.77 0.00801
## 
##  'b'  effect estimates 
##        consume   se    t    Prob
## thirst    0.45 0.15 3.09 0.00336
## 
##  'ab'  effect estimates 
##           consume boot   sd lower upper
## room_temp    0.15 0.15 0.06  0.04  0.29

Resources

MacKinnon, D. P. (2008). Introduction to statistical mediation analysis. New York, NY: Lawrence Erlbaum Associates.

Revelle, W. (2017) How to use the psych package for mediation/moderation/regression analysis. [.pdf]

Rosseel, Y. (2012). Lavaan: An R package for structural equation modeling and more. Version 0.5–12 (BETA). Journal of statistical software, 48(2), 1-36. [website]

Rucker, D. D., Preacher, K. J., Tormala, Z. L., & Petty, R. E. (2011). Mediation analysis in social psychology: Current practices and new recommendations. Social and Personality Psychology Compass, 5(6), 359-371. [.pdf]

General word of caution

Above, I listed resources prepared by experts on these and related topics. Although I generally do my best to write accurate posts, don’t assume my posts are 100% accurate or that they apply to your data or research questions. Trust statistics and methodology experts, not blog posts.

Testing indirect effects/mediation in R

What is mediation?

What is an indirect effect?

Model and Conceptual Assumptions

Libraries

Data

Print first and last five observations

Visualize relationships

Conceptual Diagram

Write model to test indirect effect using sem() from lavaan

Set random seed so results can be reproduced

Fit model

Summarize model

Print all model parameters

Interpretation

Test same model using mediation() from MBESS

Plot the mediation effect

Test the same model using mediate() from psych

Resources

General word of caution

Write model to test indirect effect using `sem()` from lavaan

Test same model using `mediation()` from MBESS

Test the same model using `mediate()` from psych