# 5 More Than One Mediator

In this chapter we’ll explore

models with more than one mediator. [We will] focus on two forms of the multiple mediator model defined by whether the mediators are linked together in a causal chain (the serial multiple mediator model) or are merely allowed to correlate bot not causally influence another mediator in the model (the parallel multiple mediator model). [We’ll] also discuss models that blend parallel and serial processes. (Andrew F. Hayes, 2018, p. 149, emphasis in the original)

## 5.1 The parallel multiple mediator model

Going from one to multiple mediators can be a big step up, conceptually. But from a model fitting perspective, it often isn’t that big of a deal. We just have more parameters.

### 5.1.1 Direct and indirect effects in a parallel multiple mediator model.

With multiple mediators, we use the language of specific indirect effects. We also add the notion of a total indirect effect, following the form

$\text{Total indirect effect of } X \text{ on } Y = \sum_{i = 1}^k a_i b_i,$

where $$k$$ is the number of mediator variables. Thus, the total effect of $$X$$ on $$Y$$ is

$c = c' + \sum_{i = 1}^k a_i b_i.$

## 5.2 Example using the presumed media influence study

Here we load a couple necessary packages, load the data, and take a glimpse().

library(tidyverse)

glimpse(pmi)
## Rows: 123
## Columns: 6
## $cond <dbl> 1, 0, 1, 0, 0, 0, 0, 1, 0, 0, 1, 0, 1, 0, 1, 1, 0, 0, 1, 0, 1, 1, 1, 1, 1, 0, 1, 0, 1, 0, 0… ##$ pmi      <dbl> 7.0, 6.0, 5.5, 6.5, 6.0, 5.5, 3.5, 6.0, 4.5, 7.0, 1.0, 6.0, 5.0, 7.0, 7.0, 7.0, 4.5, 3.5, 7…
## $import <dbl> 6, 1, 6, 6, 5, 1, 1, 6, 6, 6, 3, 3, 4, 7, 1, 6, 3, 3, 2, 4, 4, 6, 7, 4, 5, 4, 6, 5, 5, 7, 4… ##$ reaction <dbl> 5.25, 1.25, 5.00, 2.75, 2.50, 1.25, 1.50, 4.75, 4.25, 6.25, 1.25, 2.75, 3.75, 5.00, 4.00, 5…
## $gender <dbl> 1, 1, 1, 0, 1, 1, 0, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 0, 0, 1, 1, 0, 1, 0, 0, 0, 1, 0… ##$ age      <dbl> 51.0, 40.0, 26.0, 21.0, 27.0, 25.0, 23.0, 25.0, 22.0, 24.0, 22.0, 21.0, 23.0, 21.0, 22.0, 2…

library(brms)

Bayesian correlations, recall, just take an intercepts-only multivariate model.

model5.1 <-
brm(data = pmi,
family = gaussian,
mvbind(pmi, import) ~ 1,
cores = 4,
file = "fits/model05.01")

A little indexing with the posterior_summary() function will get us the Bayesian correlation with its posterior $$SD$$ and intervals.

posterior_summary(model5.1)["rescor__pmi__import", ] %>% round(digits = 3)
##  Estimate Est.Error      Q2.5     Q97.5
##     0.278     0.082     0.115     0.431

As with single mediation models, the multiple mediation model requires we carefully construct the formula for each criterion. Here we’ll use the multiple bf() approach from Chapter 3.

m1_model <- bf(import   ~ 1 + cond)
m2_model <- bf(pmi      ~ 1 + cond)
y_model  <- bf(reaction ~ 1 + import + pmi + cond)

And now we fit the model.

model5.2 <-
brm(data = pmi,
family = gaussian,
y_model + m1_model + m2_model + set_rescor(FALSE),
cores = 4,
file = "fits/model05.02")
print(model5.2, digits = 3)
##  Family: MV(gaussian, gaussian, gaussian)
##   Links: mu = identity; sigma = identity
##          mu = identity; sigma = identity
##          mu = identity; sigma = identity
## Formula: reaction ~ 1 + import + pmi + cond
##          import ~ 1 + cond
##          pmi ~ 1 + cond
##    Data: pmi (Number of observations: 123)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup samples = 4000
##
## Population-Level Effects:
##                    Estimate Est.Error l-95% CI u-95% CI  Rhat Bulk_ESS Tail_ESS
## reaction_Intercept   -0.159     0.533   -1.225    0.878 1.001     8229     3166
## import_Intercept      3.906     0.212    3.488    4.340 1.001     9383     2960
## pmi_Intercept         5.379     0.164    5.058    5.696 1.001     8461     3029
## reaction_import       0.325     0.071    0.181    0.460 1.003     7906     3403
## reaction_pmi          0.397     0.092    0.217    0.586 1.001     7868     3379
## reaction_cond         0.106     0.240   -0.360    0.587 1.000     7723     2872
## import_cond           0.628     0.306    0.044    1.224 1.001     9949     3317
## pmi_cond              0.475     0.240    0.003    0.960 1.004     9187     3054
##
## Family Specific Parameters:
##                Estimate Est.Error l-95% CI u-95% CI  Rhat Bulk_ESS Tail_ESS
## sigma_reaction    1.302     0.086    1.148    1.487 1.001     9615     2684
## sigma_import      1.730     0.112    1.524    1.968 1.001     9796     3095
## sigma_pmi         1.317     0.089    1.160    1.506 1.001    10111     2946
##
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

Because we have three criterion variables, we’ll have three Bayesian $$R^2$$ posteriors.

library(ggthemes)

bayes_R2(model5.2, summary = F) %>%
data.frame() %>%
pivot_longer(everything()) %>%
mutate(name = str_remove(name, "R2")) %>%

ggplot(aes(x = value, color = name, fill = name)) +
geom_density(alpha = .5) +
scale_color_ptol() +
scale_fill_ptol() +
scale_x_continuous(NULL, limits = 0:1) +
scale_y_continuous(NULL, breaks = NULL) +
labs(title = expression(Our~italic(R)^{2}~distributions),
subtitle = "The densities for import and pmi are asymmetric, small, and largely overlapping.\nThe density for reaction is approximately Gaussian and more impressive in magnitude.") +
theme_minimal() +
theme(legend.title = element_blank()) It’ll take a bit of data wrangling to rename our model parameters to the $$a$$, $$b$$… configuration. We’ll compute the indirect effects and $$c$$, too.

post <- posterior_samples(model5.2)

post <-
post %>%
mutate(a1 = b_import_cond,
a2 = b_pmi_cond,
b1 = b_reaction_import,
b2 = b_reaction_pmi,
c_prime = b_reaction_cond) %>%
mutate(a1b1 = a1 * b1,
a2b2 = a2 * b2) %>%
mutate(c = c_prime + a1b1 + a2b2)

Next we compute their summaries. Since Bayesians use means, medians, and sometimes the mode to describe the central tendencies of a parameter, this time we’ll mix it up and just use the median.

We’ve been summarizing our posteriors within the summarize() function. This approach gives us a lot of control. It’s also on the verbose side. Another approach is to use a family of functions from the tidybayes package. Here we’ll use median_qi() to give us the posterior medians and quantile-based 95% intervals for our parameters of interest.

library(tidybayes)

post %>%
pivot_longer(a1:c) %>%
group_by(name) %>%
median_qi(value)
## # A tibble: 8 x 7
##   name    value   .lower .upper .width .point .interval
##   <chr>   <dbl>    <dbl>  <dbl>  <dbl> <chr>  <chr>
## 1 a1      0.628  0.0436   1.22    0.95 median qi
## 2 a1b1    0.196  0.0133   0.449   0.95 median qi
## 3 a2      0.476  0.00292  0.960   0.95 median qi
## 4 a2b2    0.180  0.00137  0.435   0.95 median qi
## 5 b1      0.326  0.181    0.460   0.95 median qi
## 6 b2      0.396  0.217    0.586   0.95 median qi
## 7 c       0.496 -0.0437   1.02    0.95 median qi
## 8 c_prime 0.103 -0.360    0.587   0.95 median qi

In the value column, we have our measure of central tendency (i.e., median). The 95% intervals are in the next two columns. With tidybayes, we can ask for different kinds of intervals and different kinds of measures of central tendency, as indicated by the .width and .point columns, respectively. For example, here’s the output for the same variables when we ask for posterior means and 80% intervals.

post %>%
pivot_longer(a1:c) %>%
group_by(name) %>%
mean_qi(value, .width = .8) %>%
# for good measure
mutate_if(is_double, round, digits = 3)
## # A tibble: 8 x 7
##   name    value .lower .upper .width .point .interval
##   <chr>   <dbl>  <dbl>  <dbl>  <dbl> <chr>  <chr>
## 1 a1      0.628  0.23   1.02     0.8 mean   qi
## 2 a1b1    0.205  0.068  0.356    0.8 mean   qi
## 3 a2      0.475  0.161  0.784    0.8 mean   qi
## 4 a2b2    0.189  0.058  0.326    0.8 mean   qi
## 5 b1      0.325  0.234  0.414    0.8 mean   qi
## 6 b2      0.397  0.283  0.515    0.8 mean   qi
## 7 c       0.5    0.154  0.847    0.8 mean   qi
## 8 c_prime 0.106 -0.204  0.407    0.8 mean   qi

For more in this family of tidybayes functions, check out the Point summaries and intervals subsection of Kay’s (2020a) vignette, Extracting and visualizing tidy draws from brms models.

In the middle paragraph of page 158, Hayes showed how the mean difference in imprt between the two cond groups multiplied by b1, the coefficient of import predicting reaction, is equal to the a1b1 indirect effect. He did that with simple algebra using the group means and the point estimates, following the formula

$a_1 b_1 = \left \{ \left [\overline M_1 | (X = 1) \right] - \left[\overline M_1 | (X = 0) \right ] \right \} b_1.$

Let’s follow along. First, we’ll get those two group means and save them as numbers to arbitrary precision.

(
import_means <-
pmi %>%
group_by(cond) %>%
summarize(mean = mean(import))
)
## # A tibble: 2 x 2
##    cond  mean
##   <dbl> <dbl>
## 1     0  3.91
## 2     1  4.53
(cond_0_import_mean <- import_means[1, 2] %>% pull())
##  3.907692
(cond_1_import_mean <- import_means[2, 2] %>% pull())
##  4.534483

Here we follow the formula in the last sentence of the paragraph and then compare the results to the posterior for a1b1.

post %>%
# use Hayes's formula to make a new vector, handmade a1b1
mutate(handmade a1b1 = (cond_1_import_mean - cond_0_import_mean) * b1) %>%
# wrangle as usual
pivot_longer(c(a1b1, handmade a1b1)) %>%
group_by(name) %>%
mean_qi(value) %>%
mutate_if(is_double, round, digits = 3)
## # A tibble: 2 x 7
##   name          value .lower .upper .width .point .interval
##   <chr>         <dbl>  <dbl>  <dbl>  <dbl> <chr>  <chr>
## 1 a1b1          0.205  0.013  0.449   0.95 mean   qi
## 2 handmade a1b1 0.204  0.113  0.288   0.95 mean   qi

Yep, Hayes’s formula is good at the mean. But the distributions are distinct with vastly different posterior intervals. I’m no mathematician, so take this with a grain of salt, but I suspect this has to do with how we used fixed values (i.e., the difference of the subsample means) to compute handmade a1b1, but all the components in a1b1 were estimated.

Here we’ll follow the same protocol for a2b2.

(
pmi_means <-
pmi %>%
group_by(cond) %>%
summarize(mean = mean(pmi))
)
## # A tibble: 2 x 2
##    cond  mean
##   <dbl> <dbl>
## 1     0  5.38
## 2     1  5.85
cond_0_pmi_mean <- pmi_means[1, 2] %>% pull()
cond_1_pmi_mean <- pmi_means[2, 2] %>% pull()
post %>%
mutate(handmade a2b2 = (cond_1_pmi_mean - cond_0_pmi_mean) * b2) %>%
pivot_longer(c(a2b2, handmade a2b2)) %>%
group_by(name) %>%
mean_qi(value) %>%
mutate_if(is_double, round, digits = 3)
## # A tibble: 2 x 7
##   name          value .lower .upper .width .point .interval
##   <chr>         <dbl>  <dbl>  <dbl>  <dbl> <chr>  <chr>
## 1 a2b2          0.189  0.001  0.435   0.95 mean   qi
## 2 handmade a2b2 0.189  0.104  0.279   0.95 mean   qi

To get the total indirect effect as discussed on page 160, we simply add the a1b1 and a2b2 columns.

post <-
post %>%
mutate(total_indirect_effect = a1b1 + a2b2)

post %>%
mean_qi(total_indirect_effect) %>%
mutate_if(is_double, round, digits = 3)
## # A tibble: 1 x 6
##   total_indirect_effect .lower .upper .width .point .interval
##                   <dbl>  <dbl>  <dbl>  <dbl> <chr>  <chr>
## 1                 0.394   0.11  0.707   0.95 mean   qi

To use the equations on the top of page 161, we’ll just work directly with the original vectors in post.

post %>%
mutate(Y_bar_given_X_1 = b_import_Intercept + b_reaction_cond * 1 + b_reaction_import * b_import_Intercept + b_reaction_pmi * b_pmi_Intercept,
Y_bar_given_X_0 = b_import_Intercept + b_reaction_cond * 0 + b_reaction_import * b_import_Intercept + b_reaction_pmi * b_pmi_Intercept) %>%
mutate(c_prime by hand = Y_bar_given_X_1 - Y_bar_given_X_0) %>%
pivot_longer(c(c_prime, c_prime by hand)) %>%
group_by(name) %>%
mean_qi(value)
## # A tibble: 2 x 7
##   name            value .lower .upper .width .point .interval
##   <chr>           <dbl>  <dbl>  <dbl>  <dbl> <chr>  <chr>
## 1 c_prime         0.106 -0.360  0.587   0.95 mean   qi
## 2 c_prime by hand 0.106 -0.360  0.587   0.95 mean   qi

We computed c a while ago.

post %>%
mean_qi(c)
## # A tibble: 1 x 6
##       c  .lower .upper .width .point .interval
##   <dbl>   <dbl>  <dbl>  <dbl> <chr>  <chr>
## 1 0.500 -0.0437   1.02   0.95 mean   qi

And c minus c_prime is straight subtraction.

post %>%
mutate(c minus c_prime = c - c_prime) %>%
mean_qi(c minus c_prime)
## # A tibble: 1 x 6
##   c minus c_prime .lower .upper .width .point .interval
##               <dbl>  <dbl>  <dbl>  <dbl> <chr>  <chr>
## 1             0.394  0.110  0.707   0.95 mean   qi

## 5.3 Statistical inference

We’ve been focusing on this all along with our posterior intervals.

### 5.3.1 Inference about the direct and total effects.

We’re not going to bother with $$p$$-values and we’ve already computed the 95% Bayesian credible intervals, above. But we can examine our parameters with a density plot.

post %>%
pivot_longer(c(c, c_prime)) %>%

ggplot(aes(x = value, fill = name, color = name)) +
geom_vline(xintercept = 0, color = "black") +
geom_density(alpha = .5) +
scale_color_ptol(NULL) +
scale_fill_ptol(NULL) +
scale_y_continuous(NULL, breaks = NULL) +
labs(title = expression("It appears zero is more credible for the direct effect"~italic(c)*"', than it is the total effect, "*italic(c)*"."),
x = NULL) +
coord_cartesian(xlim = -c(-1.5, 1.5)) +
theme_minimal() ### 5.3.2 Inference about specific indirect effects.

Again, no need to worry about bootstrapping within the Bayesian paradigm. We can compute high-quality percentile-based intervals with our HMC-based posterior samples.

post %>%
pivot_longer(c(a1b1, a2b2)) %>%
group_by(name) %>%
median_qi(value) %>%
mutate_if(is.double, round, digits = 3)
## # A tibble: 2 x 7
##   name  value .lower .upper .width .point .interval
##   <chr> <dbl>  <dbl>  <dbl>  <dbl> <chr>  <chr>
## 1 a1b1  0.196  0.013  0.449   0.95 median qi
## 2 a2b2  0.18   0.001  0.435   0.95 median qi

### 5.3.3 Pairwise comparisons between specific indirect effects.

Within the Bayesian paradigm, it’s straightforward to compare indirect effects. All one has to do is compute a difference score and summarize it somehow. Here it is, a1b1 minus a2b2.

post <-
post %>%
mutate(difference = a1b1 - a2b2)

post %>%
mean_qi(difference) %>%
mutate_if(is.double, round, digits = 3)
## # A tibble: 1 x 6
##   difference .lower .upper .width .point .interval
##        <dbl>  <dbl>  <dbl>  <dbl> <chr>  <chr>
## 1      0.017 -0.303  0.339   0.95 mean   qi

Why not plot?

post %>%
ggplot(aes(x = difference)) +
geom_vline(xintercept = 0, color = "black", linetype = 2) +
geom_density(color = "black", fill = "black", alpha = .5) +
scale_x_continuous(NULL, limits = c(-1, 1)) +
scale_y_continuous(NULL, breaks = NULL) +
labs(title = "The difference score between the indirect effects",
subtitle = expression(No~italic(p)*"-value or 95% intervals needed for this one.")) +
theme_minimal() Although note well this does not mean their difference is exactly zero. The shape of the posterior distribution testifies our uncertainty in their difference. Our best bet is that the difference is approximately zero, but it could easily be plus or minus a quarter of a point or more.

### 5.3.4 Inference about the total indirect effect.

Here’s the plot.

post %>%
ggplot(aes(x = total_indirect_effect, fill = factor(0), color = factor(0))) +
geom_density(alpha = .5) +
scale_color_ptol() +
scale_fill_ptol() +
scale_y_continuous(NULL, breaks = NULL) +
labs(title = "The total indirect effect of condition on reaction",
subtitle = expression("This is the sum of"~italic(a)*italic(b)~and~italic(a)*italic(b)*". It's wide and uncertain."),
x = NULL) +
theme_minimal() +
theme(legend.position = "none") ## 5.4 The serial multiple mediator model

Examples of the parallel multiple mediator model like that described in the prior section are in abundance in the literature. A distinguishing feature of this model is the assumption that no mediator causally influences another. In practice, mediators will be correlated, but this model specified that they are not causally so. In the serial multiple mediator model, the assumption of no causal association between two or more mediators is not only relaxed, it is rejected outright a priori. The goal when an investigator estimates a serial multiple mediator model is to investigate the direct and indirect effects of $$X$$ on $$Y$$ while modeling a process in which $$X$$ causes $$M_1$$, which in turn causes $$M_2$$, and so forth, concluding with $$Y$$ as the final consequent. (p. 167, emphasis in the original)

### 5.4.1 Direct and indirect effects in a serial multiple mediator model.

In a serial multiple mediator model, the total effect of $$X$$ on $$Y$$ partitions into direct and indirect components, just as it does in the simple and parallel multiple mediator models. Regardless of the number of mediators in the model, the direct effect is $$c'$$ and interpreted as always–the estimated difference in $$Y$$ between two cases that differ by one unit on $$X$$ but that are equal on all mediators in the model. The indirect effects, of which there may be many depending on the number of mediators in the model, are all constructed by multiplying the regression weights corresponding to each step in an indirect pathway. And they are all interpreted as the estimated difference in $$Y$$ between two cases that differ by one unit on $$X$$ through the causal sequence from $$X$$ to mediator(s) to $$Y$$. Regardless of the number of mediators, the sum of all the specific indirect effects is the total indirect effect of $$X$$, and the direct and indirect effects sum to the total effect of $$X$$. (p. 170)

### 5.4.2 Statistical inference.

“In principle, Monte Carlo confidence intervals can be constructed for all indirect effects in a serial multiple mediator model” (p. 172). I’m pretty sure Hayes didn’t intend this to refer to Bayesian estimation, but I couldn’t resist the quote.

### 5.4.3 Example from the presumed media influence study.

The model syntax is similar to the earlier multiple mediator model. All we change is adding import to the list of predictors in the submodel for m2_model. But this time, let’s take the approach from last chapter where we define our bf() formulas all within brm().

model5.3 <-
brm(data = pmi,
family = gaussian,
bf(import ~ 1 + cond) +
bf(pmi ~ 1 + import + cond) +
bf(reaction ~ 1 + import + pmi + cond) +
set_rescor(FALSE),
cores = 4,
file = "fits/model05.03")
print(model5.3)
##  Family: MV(gaussian, gaussian, gaussian)
##   Links: mu = identity; sigma = identity
##          mu = identity; sigma = identity
##          mu = identity; sigma = identity
## Formula: import ~ 1 + cond
##          pmi ~ 1 + import + cond
##          reaction ~ 1 + import + pmi + cond
##    Data: pmi (Number of observations: 123)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
##          total post-warmup samples = 4000
##
## Population-Level Effects:
##                    Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## import_Intercept       3.91      0.22     3.50     4.33 1.00     9837     3064
## pmi_Intercept          4.61      0.31     3.99     5.21 1.00     8070     2924
## reaction_Intercept    -0.14      0.55    -1.22     0.90 1.00     8998     3012
## import_cond            0.63      0.31     0.03     1.24 1.00     9589     3063
## pmi_import             0.20      0.07     0.06     0.33 1.00     7796     2979
## pmi_cond               0.35      0.23    -0.11     0.81 1.00     7838     3162
## reaction_import        0.32      0.07     0.18     0.47 1.00     6835     3181
## reaction_pmi           0.40      0.10     0.20     0.59 1.00     7755     3038
## reaction_cond          0.11      0.24    -0.37     0.58 1.00     7839     3026
##
## Family Specific Parameters:
##                Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma_import       1.73      0.11     1.53     1.97 1.00     9188     2670
## sigma_pmi          1.28      0.08     1.12     1.45 1.00     8057     2837
## sigma_reaction     1.30      0.08     1.15     1.48 1.00     7701     3114
##
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

Behold the $$R^2$$ posterior densities.

bayes_R2(model5.3, summary = F) %>%
data.frame() %>%
pivot_longer(everything()) %>%
mutate(name = str_remove(name, "R2")) %>%

ggplot(aes(x = value, color = name, fill = name)) +
geom_density(alpha = .5) +
scale_color_ptol() +
scale_fill_ptol() +
scale_x_continuous(NULL, limits = 0:1) +
scale_y_continuous(NULL, breaks = NULL) +
labs(title = expression(The~italic(R)^2*" distributions for model3, the serial multiple mediator model"),
subtitle = "The density for reaction hasn't changed from model5.2. However, look how the pmi density separated from import.") +
theme_minimal() +
theme(legend.title = element_blank()) As before, here we’ll save the posterior samples into a data frame and rename the parameters a bit to match Hayes’s nomenclature.

post <-
posterior_samples(model5.3) %>%
mutate(a1 = b_import_cond,
a2 = b_pmi_cond,
b1 = b_reaction_import,
b2 = b_reaction_pmi,
c_prime = b_reaction_cond,
d21 = b_pmi_import)

Here are the parameter summaries for the pathways depicted in Figure 5.6.

post %>%
pivot_longer(a1:d21) %>%
group_by(name) %>%
mean_qi(value) %>%
mutate_if(is_double, round, digits = 3)
## # A tibble: 6 x 7
##   name    value .lower .upper .width .point .interval
##   <chr>   <dbl>  <dbl>  <dbl>  <dbl> <chr>  <chr>
## 1 a1      0.626  0.028  1.24    0.95 mean   qi
## 2 a2      0.351 -0.106  0.808   0.95 mean   qi
## 3 b1      0.324  0.182  0.47    0.95 mean   qi
## 4 b2      0.395  0.204  0.589   0.95 mean   qi
## 5 c_prime 0.107 -0.374  0.581   0.95 mean   qi
## 6 d21     0.196  0.062  0.33    0.95 mean   qi

To get our version of the parameter summaries in Table 5.2, all you have to do is add the summaries for the intercepts to what we did above.

post %>%
rename(im1 = b_import_Intercept,
im2 = b_pmi_Intercept,
iy  = b_reaction_Intercept) %>%
pivot_longer(c(a1:d21, starts_with("i"))) %>%
group_by(name) %>%
mean_qi(value) %>%
# simplify the output
select(name:.upper) %>%
mutate_if(is_double, round, digits = 3)
## # A tibble: 9 x 4
##   name     value .lower .upper
##   <chr>    <dbl>  <dbl>  <dbl>
## 1 a1       0.626  0.028  1.24
## 2 a2       0.351 -0.106  0.808
## 3 b1       0.324  0.182  0.47
## 4 b2       0.395  0.204  0.589
## 5 c_prime  0.107 -0.374  0.581
## 6 d21      0.196  0.062  0.33
## 7 im1      3.91   3.50   4.33
## 8 im2      4.61   3.99   5.21
## 9 iy      -0.143 -1.22   0.895

Here we compute the four indirect effects.

post <-
post %>%
mutate(a1b1    = a1 * b1,
a2b2    = a2 * b2,
a1d21b2 = a1 * d21 * b2) %>%
mutate(total_indirect_effect = a1b1 + a2b2 + a1d21b2)

Anticipating the skew typical of indirect effects, we’ll summarize these posteriors with medians rather than means.

post %>%
pivot_longer(a1b1:total_indirect_effect) %>%
group_by(name) %>%
median_qi(value) %>%
select(name:.upper) %>%
mutate_if(is_double, round, digits = 3)
## # A tibble: 4 x 4
##   name                  value .lower .upper
##   <chr>                 <dbl>  <dbl>  <dbl>
## 1 a1b1                  0.194  0.008  0.447
## 2 a1d21b2               0.043  0.001  0.13
## 3 a2b2                  0.13  -0.041  0.355
## 4 total_indirect_effect 0.383  0.091  0.734

To get the contrasts Hayes presented in page 179, we just do a little subtraction.

post %>%
mutate(c1 = a1b1 - a2b2,
c2 = a1b1 - a1d21b2,
c3 = a2b2 - a1d21b2) %>%
pivot_longer(c1:c3) %>%
group_by(name) %>%
median_qi(value) %>%
select(name:.upper) %>%
mutate_if(is_double, round, digits = 3)
## # A tibble: 3 x 4
##   name  value .lower .upper
##   <chr> <dbl>  <dbl>  <dbl>
## 1 c1    0.065 -0.244  0.368
## 2 c2    0.141  0.004  0.375
## 3 c3    0.085 -0.113  0.313

And just because it’s fun, we may as well plot our indirect effects.

# this will help us save a little space with the plot code
my_labels <- c(expression(italic(a)*italic(b)),
expression(italic(a)*italic(d)*italic(b)),
expression(italic(a)*italic(b)),
"total indirect effect")
# wrangle
post %>%
pivot_longer(a1b1:total_indirect_effect) %>%

# plot!
ggplot(aes(x = value, fill = name, color = name)) +
geom_density(alpha = .5) +
scale_color_ptol(NULL, labels = my_labels,
guide = guide_legend(label.hjust = 0)) +
scale_fill_ptol(NULL, labels = my_labels,
guide = guide_legend(label.hjust = 0)) +
scale_y_continuous(NULL, breaks = NULL) +
labs(title = "The four indirect effects of the serial multiple mediator model",
x = NULL) +
theme_minimal() ## 5.5 Models with parallel and serial mediation properties

In a model with two mediators, the only difference between a serial and a parallel multiple mediator model is the inclusion of a causal path from $$M_1$$ to $$M_2$$. The serial model estimates this effect, whereas the parallel model assumes it is zero, which is equivalent to leaving it out of the model entirely. With more than three mediators, a model can be a blend of parallel and serial mediation processes, depending on which paths between mediators are estimated and which are fixed to zero through their exclusion from the model. (p. 180)

## 5.6 Complementarity and competition among mediators

This chapter has been dedicated to mediation models containing more than one mediator. At this point, the benefits of estimating multiple mechanisms of influence in a single model are no doubt apparent. But the inclusion of more than one mediator in a model does entail certain risks as well, and at times the results of multiple mediator model may appear to contradict the results obtained when estimating a simpler model with a single mediator. Some of the risks, paradoxes, and contradictions that sometimes can occur are worth some acknowledgement and discussion. (p. 183)

## Session info

sessionInfo()
## R version 4.0.4 (2021-02-15)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Catalina 10.15.7
##
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
##
## locale:
##  en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
##  stats     graphics  grDevices utils     datasets  methods   base
##
## other attached packages:
##   tidybayes_3.0.0 ggthemes_4.2.4  brms_2.15.0     Rcpp_1.0.6      forcats_0.5.1   stringr_1.4.0
##   dplyr_1.0.6     purrr_0.3.4     readr_1.4.0     tidyr_1.1.3     tibble_3.1.2    ggplot2_3.3.5
##  tidyverse_1.3.1
##
## loaded via a namespace (and not attached):
##    readxl_1.3.1         backports_1.2.1      plyr_1.8.6           igraph_1.2.6         svUnit_1.0.3
##    splines_4.0.4        crosstalk_1.1.0.1    TH.data_1.0-10       rstantools_2.1.1     inline_0.3.17
##   digest_0.6.27        htmltools_0.5.1.1    rsconnect_0.8.16     fansi_0.4.2          checkmate_2.0.0
##   magrittr_2.0.1       modelr_0.1.8         RcppParallel_5.0.2   matrixStats_0.57.0   xts_0.12.1
##   sandwich_3.0-0       prettyunits_1.1.1    colorspace_2.0-0     rvest_1.0.1          ggdist_3.0.0
##   haven_2.3.1          xfun_0.23            callr_3.7.0          crayon_1.4.1         jsonlite_1.7.2
##   lme4_1.1-25          survival_3.2-10      zoo_1.8-8            glue_1.4.2           gtable_0.3.0
##   emmeans_1.5.2-1      V8_3.4.0             distributional_0.2.2 pkgbuild_1.2.0       rstan_2.21.2
##   abind_1.4-5          scales_1.1.1         mvtnorm_1.1-1        DBI_1.1.0            miniUI_0.1.1.1
##   xtable_1.8-4         stats4_4.0.4         StanHeaders_2.21.0-7 DT_0.16              htmlwidgets_1.5.3
##   httr_1.4.2           threejs_0.3.3        arrayhelpers_1.1-0   posterior_1.0.1      ellipsis_0.3.2
##   pkgconfig_2.0.3      loo_2.4.1            farver_2.1.0         sass_0.3.1           dbplyr_2.1.1
##   utf8_1.2.1           tidyselect_1.1.1     labeling_0.4.2       rlang_0.4.11         reshape2_1.4.4
##   later_1.2.0          munsell_0.5.0        cellranger_1.1.0     tools_4.0.4          cli_3.0.1
##   generics_0.1.0       broom_0.7.6          ggridges_0.5.3       evaluate_0.14        fastmap_1.1.0
##   processx_3.5.2       knitr_1.33           fs_1.5.0             nlme_3.1-152         mime_0.10
##   projpred_2.0.2       xml2_1.3.2           compiler_4.0.4       bayesplot_1.8.0      shinythemes_1.1.2
##   rstudioapi_0.13      gamm4_0.2-6          curl_4.3             reprex_2.0.0         statmod_1.4.35
##   bslib_0.2.4          stringi_1.6.2        highr_0.9            ps_1.6.0             Brobdingnag_1.2-6
##   lattice_0.20-41      Matrix_1.3-2         nloptr_1.2.2.2       markdown_1.1         tensorA_0.36.2
##  shinyjs_2.0.0        vctrs_0.3.8          pillar_1.6.1         lifecycle_1.0.0      jquerylib_0.1.4
##  bridgesampling_1.0-0 estimability_1.3     httpuv_1.6.0         R6_2.5.0             bookdown_0.22
##  promises_1.2.0.1     gridExtra_2.3        codetools_0.2-18     boot_1.3-26          colourpicker_1.1.0
##  MASS_7.3-53          gtools_3.8.2         assertthat_0.2.1     withr_2.4.2          shinystan_2.5.0
##  multcomp_1.4-16      mgcv_1.8-33          parallel_4.0.4       hms_1.1.0            grid_4.0.4
##  coda_0.19-4          minqa_1.2.4          rmarkdown_2.8        shiny_1.6.0          lubridate_1.7.10
##  base64enc_0.1-3      dygraphs_1.1.1.6

### References

Hayes, Andrew F. (2018). Introduction to mediation, moderation, and conditional process analysis: A regression-based approach (Second edition). The Guilford Press. https://www.guilford.com/books/Introduction-to-Mediation-Moderation-and-Conditional-Process-Analysis/Andrew-Hayes/9781462534654
Kay, M. (2020a). Extracting and visualizing tidy draws from brms models. https://mjskay.github.io/tidybayes/articles/tidy-brms.html