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. (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().

## Observations: 123
## Variables: 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, …
## $ 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, …
## $ 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, …
## $ 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, …
## $ 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, …
## $ 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, …

Now load brms.

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

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

##  Estimate Est.Error      Q2.5     Q97.5 
##     0.274     0.084     0.102     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.

And now we fit the model.

##  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.155     0.529   -1.220    0.897 1.001     7000     3193
## import_Intercept      3.907     0.217    3.478    4.334 1.001     8850     3188
## pmi_Intercept         5.378     0.161    5.062    5.697 1.001     9396     2705
## reaction_import       0.325     0.071    0.183    0.461 1.000     6348     2918
## reaction_pmi          0.398     0.094    0.216    0.582 1.002     6464     3227
## reaction_cond         0.099     0.238   -0.369    0.566 1.002     7431     3107
## import_cond           0.630     0.312    0.014    1.256 1.001     8770     2883
## pmi_cond              0.474     0.231    0.021    0.928 1.000     7921     2470
## 
## Family Specific Parameters: 
##                Estimate Est.Error l-95% CI u-95% CI  Rhat Bulk_ESS Tail_ESS
## sigma_reaction    1.303     0.087    1.147    1.482 1.003     6724     2801
## sigma_import      1.734     0.112    1.535    1.958 1.003     8273     3088
## sigma_pmi         1.318     0.087    1.159    1.502 1.002     8969     3099
## 
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
## is a crude measure of effective sample size, 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.

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.

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 summarising 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.

## # A tibble: 8 x 7
##   name    value   .lower .upper .width .point .interval
##   <chr>   <dbl>    <dbl>  <dbl>  <dbl> <chr>  <chr>    
## 1 a1      0.628  0.0143   1.26    0.95 median qi       
## 2 a1b1    0.197  0.00490  0.456   0.95 median qi       
## 3 a2      0.477  0.0206   0.928   0.95 median qi       
## 4 a2b2    0.180  0.00726  0.415   0.95 median qi       
## 5 b1      0.325  0.183    0.461   0.95 median qi       
## 6 b2      0.396  0.216    0.582   0.95 median qi       
## 7 c       0.495 -0.0329   1.01    0.95 median qi       
## 8 c_prime 0.101 -0.369    0.566   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.

## # A tibble: 8 x 7
##   name    value .lower .upper .width .point .interval
##   <chr>   <dbl>  <dbl>  <dbl>  <dbl> <chr>  <chr>    
## 1 a1      0.63   0.229  1.03     0.8 mean   qi       
## 2 a1b1    0.205  0.07   0.35     0.8 mean   qi       
## 3 a2      0.474  0.174  0.765    0.8 mean   qi       
## 4 a2b2    0.188  0.062  0.323    0.8 mean   qi       
## 5 b1      0.325  0.235  0.413    0.8 mean   qi       
## 6 b2      0.398  0.279  0.52     0.8 mean   qi       
## 7 c       0.492  0.149  0.844    0.8 mean   qi       
## 8 c_prime 0.099 -0.207  0.396    0.8 mean   qi

For more in this family of tidybayes functions, check out the Point summaries and intervals subsection of Kay’s helpful document, 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 = \{[\overline M_1 | (X = 1)] - [\overline M_1 | (X = 0)]\} b_1.\]

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

## # A tibble: 2 x 2
##    cond  mean
##   <dbl> <dbl>
## 1     0  3.91
## 2     1  4.53
## [1] 3.907692
## [1] 4.534483

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

## # A tibble: 2 x 7
##   name          value .lower .upper .width .point .interval
##   <chr>         <dbl>  <dbl>  <dbl>  <dbl> <chr>  <chr>    
## 1 a1b1          0.205  0.005  0.456   0.95 mean   qi       
## 2 handmade a1b1 0.203  0.115  0.289   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.

## # A tibble: 2 x 2
##    cond  mean
##   <dbl> <dbl>
## 1     0  5.38
## 2     1  5.85
## # A tibble: 2 x 7
##   name          value .lower .upper .width .point .interval
##   <chr>         <dbl>  <dbl>  <dbl>  <dbl> <chr>  <chr>    
## 1 a2b2          0.188  0.007  0.415   0.95 mean   qi       
## 2 handmade a2b2 0.189  0.103  0.277   0.95 mean   qi

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

## # A tibble: 1 x 6
##   total_indirect_effect .lower .upper .width .point .interval
##                   <dbl>  <dbl>  <dbl>  <dbl> <chr>  <chr>    
## 1                 0.393  0.119  0.708   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.

## # A tibble: 2 x 7
##   name             value .lower .upper .width .point .interval
##   <chr>            <dbl>  <dbl>  <dbl>  <dbl> <chr>  <chr>    
## 1 c_prime         0.0991 -0.369  0.566   0.95 mean   qi       
## 2 c_prime by hand 0.0991 -0.369  0.566   0.95 mean   qi

We computed c a while ago.

## # A tibble: 1 x 6
##       c  .lower .upper .width .point .interval
##   <dbl>   <dbl>  <dbl>  <dbl> <chr>  <chr>    
## 1 0.492 -0.0329   1.01   0.95 mean   qi

And c minus c_prime is straight subtraction.

## # A tibble: 1 x 6
##   `c minus c_prime` .lower .upper .width .point .interval
##               <dbl>  <dbl>  <dbl>  <dbl> <chr>  <chr>    
## 1             0.393  0.119  0.708   0.95 mean   qi

5.3 Statistical inference

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

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.

## # A tibble: 2 x 7
##   name  value .lower .upper .width .point .interval
##   <chr> <dbl>  <dbl>  <dbl>  <dbl> <chr>  <chr>    
## 1 a1b1  0.197  0.005  0.456   0.95 median qi       
## 2 a2b2  0.18   0.007  0.415   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.

## # A tibble: 1 x 6
##   difference .lower .upper .width .point .interval
##        <dbl>  <dbl>  <dbl>  <dbl> <chr>  <chr>    
## 1      0.016 -0.289  0.337   0.95 mean   qi

Why not plot?

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.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().

##  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.46     4.33 1.00     8200     3008
## pmi_Intercept          4.62      0.30     4.01     5.20 1.00     7898     3279
## reaction_Intercept    -0.15      0.54    -1.19     0.88 1.00     7148     2965
## import_cond            0.63      0.32     0.01     1.25 1.00     8085     3363
## pmi_import             0.20      0.07     0.06     0.33 1.00     7493     3110
## pmi_cond               0.35      0.23    -0.09     0.80 1.00     7073     3145
## reaction_import        0.32      0.07     0.18     0.47 1.00     8088     3161
## reaction_pmi           0.40      0.09     0.21     0.58 1.00     6706     3482
## reaction_cond          0.10      0.25    -0.38     0.58 1.01     9372     2636
## 
## Family Specific Parameters: 
##                Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma_import       1.73      0.11     1.52     1.98 1.00     8014     2541
## sigma_pmi          1.28      0.08     1.12     1.45 1.00     8945     3206
## sigma_reaction     1.30      0.08     1.15     1.48 1.00     8368     2379
## 
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample 
## is a crude measure of effective sample size, and Rhat is the potential 
## scale reduction factor on split chains (at convergence, Rhat = 1).

Behold the \(R^2\) posterior densities.

As before, here we’ll save the posterior samples into a data frame and rename the parameters a bit to match Hayes’s nomenclature.

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

## # A tibble: 6 x 7
##   name    value .lower .upper .width .point .interval
##   <chr>   <dbl>  <dbl>  <dbl>  <dbl> <chr>  <chr>    
## 1 a1      0.628  0.009  1.25    0.95 mean   qi       
## 2 a2      0.355 -0.089  0.8     0.95 mean   qi       
## 3 b1      0.325  0.184  0.466   0.95 mean   qi       
## 4 b2      0.396  0.211  0.576   0.95 mean   qi       
## 5 c_prime 0.103 -0.38   0.580   0.95 mean   qi       
## 6 d21     0.195  0.061  0.327   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.

## # A tibble: 9 x 4
##   name     value .lower .upper
##   <chr>    <dbl>  <dbl>  <dbl>
## 1 a1       0.628  0.009  1.25 
## 2 a2       0.355 -0.089  0.8  
## 3 b1       0.325  0.184  0.466
## 4 b2       0.396  0.211  0.576
## 5 c_prime  0.103 -0.38   0.580
## 6 d21      0.195  0.061  0.327
## 7 im1      3.91   3.46   4.33 
## 8 im2      4.62   4.01   5.20 
## 9 iy      -0.147 -1.19   0.884

Here we compute the four indirect effects.

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

## # A tibble: 4 x 4
##   name                  value .lower .upper
##   <chr>                 <dbl>  <dbl>  <dbl>
## 1 a1b1                  0.195  0.003  0.46 
## 2 a1d21b2               0.042  0      0.128
## 3 a2b2                  0.136 -0.035  0.353
## 4 total_indirect_effect 0.385  0.089  0.737

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

## # A tibble: 3 x 4
##   name  value .lower .upper
##   <chr> <dbl>  <dbl>  <dbl>
## 1 c1    0.06  -0.243  0.38 
## 2 c2    0.143 -0.001  0.379
## 3 c3    0.088 -0.103  0.31

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

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)

Tread carefully, friends.

Session info

## R version 3.6.0 (2019-04-26)
## Platform: x86_64-apple-darwin15.6.0 (64-bit)
## Running under: macOS High Sierra 10.13.6
## 
## Matrix products: default
## BLAS:   /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRblas.0.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/3.6/Resources/lib/libRlapack.dylib
## 
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
## 
## attached base packages:
## [1] stats     graphics  grDevices utils     datasets  methods   base     
## 
## other attached packages:
##  [1] tidybayes_1.1.0 ggthemes_4.2.0  brms_2.10.3     Rcpp_1.0.2      forcats_0.4.0   stringr_1.4.0  
##  [7] dplyr_0.8.3     purrr_0.3.3     readr_1.3.1     tidyr_1.0.0     tibble_2.1.3    ggplot2_3.2.1  
## [13] tidyverse_1.2.1
## 
## loaded via a namespace (and not attached):
##  [1] colorspace_1.4-1          ggridges_0.5.1            rsconnect_0.8.15          ggstance_0.3.2           
##  [5] markdown_1.1              base64enc_0.1-3           rstudioapi_0.10           rstan_2.19.2             
##  [9] svUnit_0.7-12             DT_0.9                    fansi_0.4.0               lubridate_1.7.4          
## [13] xml2_1.2.0                bridgesampling_0.7-2      knitr_1.23                shinythemes_1.1.2        
## [17] zeallot_0.1.0             bayesplot_1.7.0           jsonlite_1.6              broom_0.5.2              
## [21] shiny_1.3.2               compiler_3.6.0            httr_1.4.0                backports_1.1.5          
## [25] assertthat_0.2.1          Matrix_1.2-17             lazyeval_0.2.2            cli_1.1.0                
## [29] later_1.0.0               htmltools_0.4.0           prettyunits_1.0.2         tools_3.6.0              
## [33] igraph_1.2.4.1            coda_0.19-3               gtable_0.3.0              glue_1.3.1.9000          
## [37] reshape2_1.4.3            cellranger_1.1.0          vctrs_0.2.0               nlme_3.1-139             
## [41] crosstalk_1.0.0           xfun_0.10                 ps_1.3.0                  rvest_0.3.4              
## [45] mime_0.7                  miniUI_0.1.1.1            lifecycle_0.1.0           gtools_3.8.1             
## [49] zoo_1.8-6                 scales_1.0.0              colourpicker_1.0          hms_0.4.2                
## [53] promises_1.1.0            Brobdingnag_1.2-6         parallel_3.6.0            inline_0.3.15            
## [57] shinystan_2.5.0           gridExtra_2.3             loo_2.1.0                 StanHeaders_2.19.0       
## [61] stringi_1.4.3             dygraphs_1.1.1.6          pkgbuild_1.0.5            rlang_0.4.1              
## [65] pkgconfig_2.0.3           matrixStats_0.55.0        evaluate_0.14             lattice_0.20-38          
## [69] rstantools_2.0.0          htmlwidgets_1.5           labeling_0.3              tidyselect_0.2.5         
## [73] processx_3.4.1            plyr_1.8.4                magrittr_1.5              R6_2.4.0                 
## [77] generics_0.0.2            pillar_1.4.2              haven_2.1.0               withr_2.1.2              
## [81] xts_0.11-2                abind_1.4-5               modelr_0.1.4              crayon_1.3.4             
## [85] arrayhelpers_1.0-20160527 utf8_1.1.4                rmarkdown_1.13            grid_3.6.0               
## [89] readxl_1.3.1              callr_3.3.2               threejs_0.3.1             digest_0.6.21            
## [93] xtable_1.8-4              httpuv_1.5.2              stats4_3.6.0              munsell_0.5.0            
## [97] shinyjs_1.0