3  The Simple Mediation Model

Hayes closed the opening with:

Whereas answering questions about when or for whom is the domain of moderation analysis, questions that ask about how pertain to mediation, the focus of this and the next three chapters. In this chapter, [we’ll explore] the simple mediation model and illustrate using OLS regression-based Bayesian path analysis how the effect of an antecedent variable \(X\) on some final consequent \(Y\) can be partitioned into two paths of influence, direct and indirect. (Andrew F. Hayes, 2018, p. 78, emphasis in the original)

3.1 The simple mediation model

Mediation analysis is a statistical method used to evaluate evidence from studies designed to test hypotheses about how some causal antecedent variable \(X\) transmits its effect on a consequent variable \(Y\).

When thinking about whether a phenomenon or theory you are studying could be conceptualized as a mediation process, it is important to keep in mind that mediation is ultimately a causal explanation. It is assumed that the relationships in the system are causal, and, importantly, that \(M\) is causally located between \(X\) and \(Y\). It must be assumed, if not also empirically substantiated, that \(X\) causes \(M\), which in turn causes \(Y\). \(M\) cannot possibly carry \(X\)’s effect on \(Y\) if \(M\) is not located causally between \(X\) and \(Y\). (pp. 78–81, emphasis in the original)

3.2 Estimation of the direct, indirect, and total effects of \(X\)

Given the simple three-term mediation model, the statistical model is expressed in the two equations

\[\begin{align*} M & = i_M + a X + e_M \\ Y & = i_Y + c' X + b M + e_Y. \end{align*}\]

When using OLS software, as Hayes promotes throughout the text, these equations are estimated sequentially. However, the brms package has multivariate capabilities. As such, our results will be from a Bayesian multivariate model that computes both equations at once. They’re both part of a joint model. When we consider more advanced models later in the text, our multivariate models will fit even more than two equations at once. None of this is a problem for brms.

3.3 Example with dichotomous \(X\): The influence of presumed media influence

Here we load a couple necessary packages, load the data, and take a peek.

library(tidyverse)

pmi <- read_csv("data/pmi/pmi.csv")

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, 1, 0, 0, 1, 0, 1, …
$ 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.0, 7.0, 4.5, 6.0, 6…
$ 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, 1, 4, 5, 5, 4, 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.00, 3.50, 2.25, 4.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, 0, 0, 0, 1, 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, 23.0, 23.0, 23.0, 22.…

You can get the male/female split like so.

pmi %>% 
  count(gender)

# A tibble: 2 × 2
  gender     n
   <dbl> <int>
1      0    80
2      1    43

Here is the split by condition.

pmi %>% 
  count(cond)

# A tibble: 2 × 2
   cond     n
  <dbl> <int>
1     0    65
2     1    58

Here is how to get the ungrouped mean and \(\textit{SD}\) values for reaction and pmi, as presented in Table 3.1.

pmi %>% 
  pivot_longer(c(reaction, pmi)) %>% 
  group_by(name) %>% 
  summarise(mean = mean(value),
            sd   = sd(value)) %>% 
  mutate_if(is.double, round, digits = 3)

# A tibble: 2 × 3
  name      mean    sd
  <chr>    <dbl> <dbl>
1 pmi       5.60  1.32
2 reaction  3.48  1.55

You might get the mean and \(\textit{SD}\) values for reaction and pmi grouped by cond like this.

pmi %>% 
  pivot_longer(c(reaction, pmi)) %>% 
  group_by(cond, name) %>% 
  summarise(mean = mean(value),
            sd   = sd(value)) %>% 
  mutate_if(is.double, round, digits = 3)

# A tibble: 4 × 4
# Groups:   cond [2]
   cond name      mean    sd
  <dbl> <chr>    <dbl> <dbl>
1     0 pmi       5.38  1.34
2     0 reaction  3.25  1.61
3     1 pmi       5.85  1.27
4     1 reaction  3.75  1.45

Let’s load our primary statistical package, brms.

library(brms)

Before we begin, I should acknowledge that I greatly benefited by this great blog post on path analysis in brms by Jarrett Byrnes. With brms, we handle mediation models using the multivariate syntax (Bürkner, 2022). There are a few ways to do this. Let’s start simple.

If you look at the path model in Figure 3.3, you’ll note that reaction is predicted by pmi and cond. pmi, in turn, is predicted solely by cond. So we have two regression models, which is just the kind of thing the brms multivariate syntax is for. So first let’s specify both models, which we’ll nest in bf() functions and save as objects.

y_model <- bf(reaction ~ 1 + pmi + cond)
m_model <- bf(pmi ~ 1 + cond)

Now we have our bf() objects in hand, we’ll combine them with the + operator within the brm() function. We’ll also specify set_rescor(FALSE) because we are not interested in adding a residual correlation between reaction and pmi.

model3.1 <- brm(
  data = pmi, 
  family = gaussian,
  y_model + m_model + set_rescor(FALSE),
  cores = 4,
  file = "fits/model03.01")

Here are our results.

print(model3.1)

 Family: MV(gaussian, gaussian) 
  Links: mu = identity
         mu = identity 
Formula: reaction ~ 1 + pmi + cond 
         pmi ~ 1 + cond 
   Data: pmi (Number of observations: 123) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Regression Coefficients:
                   Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
reaction_Intercept     0.53      0.56    -0.55     1.61 1.00     6657     3384
pmi_Intercept          5.38      0.16     5.07     5.70 1.00     5550     2884
reaction_pmi           0.51      0.10     0.31     0.70 1.00     6492     3241
reaction_cond          0.26      0.25    -0.24     0.76 1.00     5018     2931
pmi_cond               0.48      0.23     0.02     0.94 1.00     5820     3049

Further Distributional Parameters:
               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
sigma_reaction     1.40      0.09     1.24     1.60 1.00     5432     2758
sigma_pmi          1.32      0.08     1.16     1.49 1.00     5453     3090

Draws were sampled 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).

If you compare our model summary with the coefficients in the path model in Figure 3.3, you’ll see our coefficients are the same. The brms summary also includes intercepts and residual variances, which are typically omitted in path diagrams even though they’re still part of the model.

If you’re getting lost in all the model output, try taking out the constant and error terms.

fixef(model3.1)[3:5, ] %>% round(digits = 3)

              Estimate Est.Error   Q2.5 Q97.5
reaction_pmi     0.506     0.099  0.313 0.702
reaction_cond    0.258     0.250 -0.244 0.756
pmi_cond         0.485     0.233  0.018 0.940

In his Table 3.2, Hayes included the \(R^2\) values. Here are ours.

bayes_R2(model3.1) %>% round(digits = 3)

           Estimate Est.Error Q2.5 Q97.5
R2reaction    0.209     0.057  0.1 0.317
R2pmi         0.040     0.030  0.0 0.113

It’s worth it to actually plot the \(R^2\) distributions. We’ll take our color palette from the ggthemes package (Arnold, 2021).

library(ggthemes)

bayes_R2(model3.1, summary = F) %>% 
  data.frame() %>% 
  pivot_longer(everything()) %>% 
  mutate(name = str_remove(name, "R2")) %>% 
  
  ggplot(aes(x = value, fill = name)) +
  geom_density(color = "transparent", alpha = 2/3) +
  scale_fill_colorblind(NULL) +  # we got this color palette from the ggthemes package
  coord_cartesian(xlim = 0:1) +
  labs(title = expression(paste("The ", italic("R")^{2}, " distributions for model3.1")),
       x = NULL) +
  theme_classic()

We went through the trouble of plotting the \(R^2\) distributions because it’s useful to understand that they won’t often be symmetric when they’re near their logical boundaries (i.e., 0 and 1). This is where asymmetric Bayesian credible intervals can really shine.

Let’s get down to business and examine the indirect effect, the \(ab\) pathway. In our model,

• \(a\) = pmi_cond and
• \(b\) = reaction_pmi.

You can isolate them with fixef()[i, ].

fixef(model3.1)[5, ]

  Estimate  Est.Error       Q2.5      Q97.5 
0.48464685 0.23299553 0.01787585 0.93982559 

fixef(model3.1)[3, ]

  Estimate  Est.Error       Q2.5      Q97.5 
0.50620027 0.09875084 0.31318150 0.70186569 

So the naive approach would be to just multiply them.

(fixef(model3.1)[5, ] * fixef(model3.1)[3, ]) %>% round(digits = 3)

 Estimate Est.Error      Q2.5     Q97.5 
    0.245     0.023     0.006     0.660 

Now, this does get us the correct ‘Estimate’ (i.e., posterior mean). However, the posterior \(\textit{SD}\) and 95% intervals are off. If you want to do this properly, you need to work with the poster samples themselves. We do that with the as_draws_df() function.

draws <- as_draws_df(model3.1)

glimpse(draws)

Rows: 4,000
Columns: 14
$ b_reaction_Intercept <dbl> 0.884817256, 0.709373041, 1.246341949, -0.184538213, 0.297656358, -0.313754147, 1.130191959, -0.388…
$ b_pmi_Intercept      <dbl> 5.344391, 5.152021, 5.611714, 5.102691, 5.236361, 5.616774, 5.256399, 5.200104, 5.512576, 5.350389,…
$ b_reaction_pmi       <dbl> 0.4467195, 0.4015191, 0.3969717, 0.5977755, 0.5274852, 0.6683733, 0.3662360, 0.6694417, 0.6382508, …
$ b_reaction_cond      <dbl> -0.0411748546, 0.5160009658, 0.2991970528, 0.3947024041, 0.3395527263, -0.0131784151, 0.5578033005,…
$ b_pmi_cond           <dbl> 0.2946610, 0.5237199, 0.4640553, 0.5611707, 0.4044745, 0.2125948, 0.6780048, 0.4123121, 0.2576005, …
$ sigma_reaction       <dbl> 1.603663, 1.470381, 1.358391, 1.447387, 1.454218, 1.463181, 1.450798, 1.395807, 1.415815, 1.377325,…
$ sigma_pmi            <dbl> 1.375654, 1.436679, 1.311643, 1.290451, 1.315679, 1.319538, 1.331141, 1.446826, 1.435405, 1.202506,…
$ Intercept_reaction   <dbl> 3.367757, 3.201850, 3.611114, 3.350097, 3.412545, 3.424009, 3.444738, 3.392543, 3.506240, 3.470382,…
$ Intercept_pmi        <dbl> 5.483337, 5.398978, 5.830536, 5.367308, 5.427089, 5.717022, 5.576109, 5.394528, 5.634047, 5.634196,…
$ lprior               <dbl> -6.763072, -6.747718, -6.666736, -6.709087, -6.712288, -6.690134, -6.701261, -6.734492, -6.718296, …
$ lp__                 <dbl> -430.8661, -432.1431, -429.2444, -429.6944, -427.8512, -429.2601, -428.2888, -430.5593, -428.7798, …
$ .chain               <int> 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, …
$ .iteration           <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, …
$ .draw                <int> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 26, 27, …

Here we compute the indirect effect, ab.

draws <- draws %>% 
  mutate(ab = b_pmi_cond * b_reaction_pmi)

Now we have ab as a properly computed vector, we can summarize it with the quantile() function.

quantile(draws$ab, probs = c(0.5, 0.025, 0.975)) %>% 
  round(digits = 3)

  50%  2.5% 97.5% 
0.239 0.008 0.520 

And we can even visualize it as a density.

draws %>% 
  ggplot(aes(x = ab)) +
  geom_density(color = "transparent", 
               fill = colorblind_pal()(3)[3]) +
  scale_y_continuous(NULL, breaks = NULL) +
  labs(title = expression(paste("Our indirect effect, the ", italic("ab"), " pathway")),
       x = NULL) +
  theme_classic()

It’s also worth pointing out that as the indirect effect isn’t perfectly symmetric, its mean and median aren’t quite the same.

draws %>% 
  summarize(mean   = mean(ab),
            median = median(ab)) %>% 
  round(digits = 3)

# A tibble: 1 × 2
   mean median
  <dbl>  <dbl>
1 0.246  0.239

Their magnitudes are similar, but this asymmetry will be a source of contrast to our estimates and the OLS estimates Hayes reported throughout his text. This is also something to consider when reporting on central tendency. When the indirect effect–or any other parameter, for that matter–is quite asymmetric, you might prefer reporting the median rather than the mean.

On page 90, Hayes computed the adjusted means for \(Y\). For both cond == 1 and cond == 0, he computed the expected values for reaction when pmi was at its mean. A natural way to do that in brms is with fitted(). First we’ll put our input values for cond and pmi in a tibble, which we’ll call nd. Then we’ll feed nd into the newdata argument within the fitted() function.

nd <- tibble(cond = 1:0,
             pmi  = mean(pmi$pmi))

fitted(model3.1, newdata = nd)

, , reaction

     Estimate Est.Error     Q2.5    Q97.5
[1,] 3.621997 0.1782274 3.265353 3.975661
[2,] 3.364489 0.1735926 3.015155 3.704067

, , pmi

     Estimate Est.Error    Q2.5    Q97.5
[1,] 5.861110 0.1736487 5.52749 6.204191
[2,] 5.376463 0.1600364 5.07328 5.695525

Because model3.1 is a multivariate model, fitted() returned the model-implied summaries for both reaction and pmi. If you just want the adjusted means for reaction, you can use the resp argument within fitted().

fitted(model3.1, newdata = nd, resp = "reaction") %>% round(digits = 3)

     Estimate Est.Error  Q2.5 Q97.5
[1,]    3.622     0.178 3.265 3.976
[2,]    3.364     0.174 3.015 3.704

Note how this is where the two values in the \(Y\) adjusted column in Table 3.1 came from.

However, if we want to reproduce how Hayes computed the total effect (i.e., \(c' + ab\)), we’ll need to work with the posterior draws themselves, draws. Recall, we’ve already saved the indirect effect as a vector, ab. The direct effect, \(c'\), is labeled b_reaction_cond within draws. To get the total effect, \(c\), all we need to is add those vectors together.

draws <- draws %>% 
  mutate(total_effect = b_reaction_cond + ab)

Here are the posterior mean with its quantile-based 95% intervals.

draws %>% 
  summarize(mean = mean(total_effect),
            ll   = quantile(total_effect, prob = 0.025),
            ul   = quantile(total_effect, prob = 0.975))

# A tibble: 1 × 3
   mean      ll    ul
  <dbl>   <dbl> <dbl>
1 0.503 -0.0371  1.04

3.3.1 Estimation of the model in PROCESS for SPSS and SAS

Nothing new for us, here.

3.4 Statistical inference

Our approach will not match up neatly with Hayes’s on this topic.

3.4.1 Inference about the total effect of \(X\) on \(Y\)

As we mentioned in Chapter 2, we can indeed focus on rejecting \(H_0\) when using Bayesian statistics. I, however, am not a fan of that approach and I will not be focusing on Bayesian \(p\)-values. But throughout this project, we will make great efforts to express the (un)certainty in our models with various plots of posterior distributions and summary statistics, such as measures of central tendency (e.g., means) and spread (e.g., percentile-based 95% intervals).

So instead of \(t\)- and \(p\)-values for \(c'\), we are going to focus on the distribution. We already gave the mean and 95% intervals, above. Here’s a look at the density.

 draws %>% 
  ggplot(aes(x = total_effect)) +
  geom_density(color = "transparent", 
               fill = colorblind_pal()(3)[2]) +
  scale_y_continuous(NULL, breaks = NULL) +
  xlab(expression(paste(italic(c)," (i.e., the total effect)"))) +
  theme_classic()

3.4.2 Inference about the direct effect of \(X\) on \(Y\)

Like in the last section, we will just look at the posterior distribution for the direct effect (i.e., \(c'\), b_reaction_cond).

draws %>% 
  ggplot(aes(x = b_reaction_cond)) +
  geom_density(color = "transparent", 
               fill = colorblind_pal()(4)[4]) +
  geom_vline(xintercept = 0, color = "white", linetype = 2) +
  scale_y_continuous(NULL, breaks = NULL) +
  labs(title = expression(paste("Yep, 0 is a credible value for ", italic("c"), ".")),
       x = NULL) +
  theme_classic()

If you wanted to quantify what proportion of the density was less than 0, you could execute something like this.

draws %>% 
  summarize(proportion_below_zero = mean(b_reaction_cond < 0))

# A tibble: 1 × 1
  proportion_below_zero
                  <dbl>
1                 0.151

This is something like a Bayesian \(p\)-value. But of course, you could always just look at the posterior intervals.

posterior_interval(model3.1)["b_reaction_cond", ]

      2.5%      97.5% 
-0.2439048  0.7558042 

3.4.3 Inference about the indirect of \(X\) on \(Y\) through \(M\)

The indirect effect quantifies how much two cases that differ by a unit on \(X\) are estimated to differ on \(Y\) as a result of \(X\)’s influence on \(M\), which in turn influences \(Y\). The indirect effect is relevant as to [what extent] \(X\)’s effect on \(Y\) can be said to be transmitted through the mechanism represented by the \(X \rightarrow M \rightarrow Y\) causal chain of events. (p. 95)

3.4.3.1 The normal theory approach

This is not our approach.

3.4.3.2 Bootstrap confidence interval

This is not our approach.

However, Markov chain Monte Carlo (MCMC) methods are iterative and share some characteristics with boostrapping. On page 98, Hayes outlined 6 steps for constructing the \(ab\) bootstrap confidence interval. Here are our responses to those steps w/r/t Bayes with MCMC–or in our case HMC (i.e., Hamiltonian Monte Carlo).

If HMC or MCMC, in general, are new to you, you might check out this lecture or the Stan Reference Manual if you’re more technically oriented.

Anyway, Hayes’s 6 steps:

3.4.3.2.1 Step 1

With HMC we do not take random samples of the data themselves. Rather, we take random draws from the posterior distribution. The posterior distribution is the joint probability distribution of our model.

3.4.3.2.2 Step 2

After we fit our model with the brm() function and save our posterior draws in a data frame (i.e., draws <- as_draws_df(my_model_fit)), we then make a new column (a.k.a. vector, variable) that is the product of our coefficients for the \(a\) and \(b\) pathways. In the example above, this looked like draws %>% mutate(ab = b_pmi_cond * b_reaction_pmi). Let’s take a look at those columns.

draws %>% 
  select(b_pmi_cond, b_reaction_pmi, ab) %>% 
  slice(1:10)

# A tibble: 10 × 3
   b_pmi_cond b_reaction_pmi    ab
        <dbl>          <dbl> <dbl>
 1      0.295          0.447 0.132
 2      0.524          0.402 0.210
 3      0.464          0.397 0.184
 4      0.561          0.598 0.335
 5      0.404          0.527 0.213
 6      0.213          0.668 0.142
 7      0.678          0.366 0.248
 8      0.412          0.669 0.276
 9      0.258          0.638 0.164
10      0.602          0.366 0.220

Our draws data frame has 4,000 rows. Why 4,000? By default, brms runs 4 HMC chains. Each chain has 2,000 iterations, 1,000 of which are warmups, which we always discard. As such, there are 1,000 good iterations left in each chain and \(1{,}000 \times 4 = 4{,}000\). We can change these defaults as needed.

Each row in draws contains the parameter values based on one of those draws. And again, these are draws from the posterior distribution. They are not draws from the data.

3.4.3.2.3 Step 3

We don’t refit the model \(k\) times based on the samples from the data. We take a number of draws from the posterior distribution. Hayes likes to take 5,000 samples when he bootstraps. Happily, that number is quite similar to our default 4,000 HMC draws. Whether 5,000, 4,000 or 10,000, these are all large enough numbers that the distributions become fairly stable. With HMC, however, you might want to increase the number of iterations if either measure of effective sample size, ‘Bulk_ESS’ and ‘Tail_ESS’ in the print() output, are substantially smaller than the number of iterations.

3.4.3.2.4 Step 4

When we use the quantile() function to compute our Bayesian credible intervals, we’ve sorted. Conceptually, we’ve done this.

draws %>% 
  select(ab) %>% 
  arrange(ab) %>% 
  slice(1:10)

# A tibble: 10 × 1
        ab
     <dbl>
 1 -0.230 
 2 -0.154 
 3 -0.145 
 4 -0.140 
 5 -0.114 
 6 -0.113 
 7 -0.0954
 8 -0.0932
 9 -0.0901
10 -0.0804

3.4.3.2.5 Step 5

Yes, this is what we do, too.

ci <- 95

0.5 * (100 - ci)

[1] 2.5

3.4.3.2.6 Step 6

This is also what we do.

ci <- 95

(100 - 0.5 * (100 - ci))

[1] 97.5

Also, notice the headers in the rightmost two columns in our posterior_summary() output:

posterior_summary(model3.1)

                         Estimate  Est.Error          Q2.5        Q97.5
b_reaction_Intercept    0.5289444 0.55759464   -0.55309143    1.6079577
b_pmi_Intercept         5.3764629 0.16003640    5.07327992    5.6955248
b_reaction_pmi          0.5062003 0.09875084    0.31318150    0.7018657
b_reaction_cond         0.2575082 0.25043123   -0.24390475    0.7558042
b_pmi_cond              0.4846469 0.23299553    0.01787585    0.9398256
sigma_reaction          1.4038404 0.09056935    1.23587566    1.5987953
sigma_pmi               1.3153359 0.08415618    1.16488606    1.4894074
Intercept_reaction      3.4859156 0.12359260    3.24843622    3.7269914
Intercept_pmi           5.6049956 0.11927172    5.37620847    5.8419304
lprior                 -6.6883996 0.03594460   -6.76542468   -6.6242599
lp__                 -429.5784108 1.89073892 -434.18604665 -426.9614705

Those .025 and .975 quantiles from above are just what brms is giving us in our 95% Bayesian credible intervals.

Here’s our version of Figure 3.5.

# these will come in handy in the subtitle
ll <- quantile(draws$ab, probs = 0.025) %>% round(digits = 3)
ul <- quantile(draws$ab, probs = 0.975) %>% round(digits = 3)

draws %>% 
  ggplot(aes(x = ab)) +
  geom_histogram(color = "white", linewidth = 0.25, 
                 fill = colorblind_pal()(5)[5],
                 binwidth = 0.025, boundary = 0) +
  geom_vline(xintercept = quantile(draws$ab, probs = c(0.025, 0.975)),
             linetype = 3, color = colorblind_pal()(6)[6]) +
  labs(subtitle = str_c("95% of the posterior draws are between ", ll, " and ", ul),
       x = expression(Indirect~effect~(italic(ab))),
       y = "Frequency in 4,000 HMC posterior draws") +
  theme_classic()

Again, as Hayes discussed how to specify different types of intervals in PROCESS on page 102, you can ask for different kinds of intervals in your print() or summary() output with the probs argument, just as you can with quantile() when working directly with the posterior draws.

Hayes discussed setting the seed in PROCESS on page 104. You can do this with the seed argument in the brm() function, too.

3.4.3.3 Alternative “asymmetric” confidence interval approaches

This section does not quite refer to us. I’m a little surprised Hayes didn’t at least dedicate a paragraph or two on Bayesian estimation. Sure, he mentioned Monte Carlo, but not within the context of Bayes. So it goes… But if you’re interested, you can read about Bayesian intervals for mediation models in Yuan and MacKinnon’s (2009) Bayesian mediation analysis. And yes, Hayes is aware of this. He has cited it in his work (e.g., Andrew F. Hayes, 2015).

3.5 An example with continuous \(X\): Economic stress among small-business owners

Here’s the estress data.

estress <- read_csv("data/estress/estress.csv")

glimpse(estress)

Rows: 262
Columns: 7
$ tenure   <dbl> 1.67, 0.58, 0.58, 2.00, 5.00, 9.00, 0.00, 2.50, 0.50, 0.58, 9.00, 1.92, 2.00, 1.42, 0.92, 20.00, 2.00, 10.00, 0…
$ estress  <dbl> 6.0, 5.0, 5.5, 3.0, 4.5, 6.0, 5.5, 3.0, 5.5, 6.0, 5.5, 4.0, 3.0, 2.5, 3.5, 6.0, 4.0, 6.0, 3.5, 4.0, 2.5, 4.0, 2…
$ affect   <dbl> 2.60, 1.00, 2.40, 1.16, 1.00, 1.50, 1.00, 1.16, 1.33, 3.00, 3.00, 2.00, 1.83, 1.16, 1.16, 1.00, 1.33, 1.50, 1.3…
$ withdraw <dbl> 3.00, 1.00, 3.66, 4.66, 4.33, 3.00, 1.00, 1.00, 2.00, 4.00, 4.33, 1.00, 5.00, 1.66, 4.00, 1.33, 1.66, 1.66, 2.0…
$ sex      <dbl> 1, 0, 1, 1, 1, 1, 0, 0, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 1, 1, 1, 0, 1, 0, 0, 0, 1, 0, 1, 1, 1, 0, 1, 1, …
$ age      <dbl> 51, 45, 42, 50, 48, 48, 51, 47, 40, 43, 57, 36, 33, 29, 33, 48, 40, 45, 37, 42, 54, 57, 37, 49, 40, 30, 49, 31,…
$ ese      <dbl> 5.33, 6.05, 5.26, 4.35, 4.86, 5.05, 3.66, 6.13, 5.26, 4.00, 2.53, 6.60, 5.20, 5.66, 5.66, 5.40, 6.00, 6.13, 5.4…

The model set up is just like before. There are no complications switching from a binary \(X\) variable to a continuous one.

y_model <- bf(withdraw ~ 1 + estress + affect)
m_model <- bf(affect ~ 1 + estress)

With our y_model and m_model defined, we’re ready to fit.

model3.2 <- brm(
  data = estress, 
  family = gaussian,
  y_model + m_model + set_rescor(FALSE),
  cores = 4,
  file = "fits/model03.02")

Let’s take a look at the results.

print(model3.2, digits = 3)

 Family: MV(gaussian, gaussian) 
  Links: mu = identity
         mu = identity 
Formula: withdraw ~ 1 + estress + affect 
         affect ~ 1 + estress 
   Data: estress (Number of observations: 262) 
  Draws: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
         total post-warmup draws = 4000

Regression Coefficients:
                   Estimate Est.Error l-95% CI u-95% CI  Rhat Bulk_ESS Tail_ESS
withdraw_Intercept    1.443     0.254    0.953    1.942 1.000     6799     2866
affect_Intercept      0.802     0.143    0.528    1.083 1.002     5511     2931
withdraw_estress     -0.076     0.051   -0.176    0.025 1.001     6085     3307
withdraw_affect       0.769     0.102    0.564    0.968 1.000     4736     3170
affect_estress        0.172     0.030    0.114    0.230 1.003     5317     2851

Further Distributional Parameters:
               Estimate Est.Error l-95% CI u-95% CI  Rhat Bulk_ESS Tail_ESS
sigma_withdraw    1.137     0.049    1.045    1.239 1.002     6125     3206
sigma_affect      0.685     0.030    0.631    0.746 1.001     5435     3025

Draws were sampled 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).

The ‘Rhat’, ‘Bulk_ESS’, and ‘Tail_ESS’ values look great. Happily, the values in our summary cohere well with those Hayes reported in Table 3.5. Here are our \(R^2\) values.

bayes_R2(model3.2)

            Estimate  Est.Error       Q2.5     Q97.5
R2withdraw 0.1823843 0.03789610 0.10751336 0.2577403
R2affect   0.1163203 0.03423755 0.05327693 0.1859644

These are also quite similar to those in the text. Here’s our indirect effect.

# putting the posterior draws into a data frame
draws <- as_draws_df(model3.2)

# computing the ab coefficient with multiplication
draws <- draws %>% 
  mutate(ab = b_affect_estress*b_withdraw_affect)

# getting the posterior median and 95% intervals with `quantile()`
quantile(draws$ab, probs = c(0.5, 0.025, 0.975)) %>% round(digits = 3)

  50%  2.5% 97.5% 
0.131 0.080 0.193 

We can visualize its shape, median, and 95% intervals in a density plot.

draws %>% 
  ggplot(aes(x = ab)) +
  geom_density(color = "transparent", 
               fill = colorblind_pal()(7)[7]) +
  geom_vline(xintercept = quantile(draws$ab, probs = c(0.025, 0.5, 0.975)), 
             color = "white", linetype = c(2, 1, 2), linewidth = c(0.5, 0.8, 0.5)) +
  scale_x_continuous(breaks = quantile(draws$ab, probs = c(0.025, 0.5, 0.975)),
                     labels = quantile(draws$ab, probs = c(0.025, 0.5, 0.975)) %>% round(2) %>% as.character()) +
  scale_y_continuous(NULL, breaks = NULL) +
  labs(title = expression(Behold~our~italic("ab")*'!'),
       x = NULL) +
  theme_classic()

Here’s \(c'\), the direct effect of esterss predicting withdraw.

posterior_summary(model3.2)["b_withdraw_estress", ]

   Estimate   Est.Error        Q2.5       Q97.5 
-0.07571124  0.05094141 -0.17600229  0.02457440 

It has wide flapping intervals which do straddle zero. A little addition will give us the direct effect, \(c\).

draws <- draws %>% 
  mutate(c = b_withdraw_estress + ab)

quantile(draws$c, probs = c(0.5, 0.025, 0.975)) %>% round(digits = 3)

   50%   2.5%  97.5% 
 0.056 -0.050  0.167 

Session info

sessionInfo()

R version 4.5.1 (2025-06-13)
Platform: aarch64-apple-darwin20
Running under: macOS Ventura 13.4

Matrix products: default
BLAS:   /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRblas.0.dylib 
LAPACK: /Library/Frameworks/R.framework/Versions/4.5-arm64/Resources/lib/libRlapack.dylib;  LAPACK version 3.12.1

locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8

time zone: America/Chicago
tzcode source: internal

attached base packages:
[1] stats     graphics  grDevices utils     datasets  methods   base     

other attached packages:
 [1] ggthemes_5.1.0  brms_2.23.0     Rcpp_1.1.0      lubridate_1.9.4 forcats_1.0.1   stringr_1.6.0   dplyr_1.1.4    
 [8] purrr_1.2.1     readr_2.1.5     tidyr_1.3.2     tibble_3.3.1    ggplot2_4.0.1   tidyverse_2.0.0

loaded via a namespace (and not attached):
 [1] tidyselect_1.2.1        farver_2.1.2            loo_2.9.0.9000          S7_0.2.1                fastmap_1.2.0          
 [6] TH.data_1.1-4           tensorA_0.36.2.1        digest_0.6.39           estimability_1.5.1      timechange_0.3.0       
[11] lifecycle_1.0.5         StanHeaders_2.36.0.9000 survival_3.8-3          magrittr_2.0.4          posterior_1.6.1.9000   
[16] compiler_4.5.1          rlang_1.1.7             tools_4.5.1             utf8_1.2.6              knitr_1.51             
[21] labeling_0.4.3          bridgesampling_1.2-1    htmlwidgets_1.6.4       curl_7.0.0              bit_4.6.0              
[26] pkgbuild_1.4.8          plyr_1.8.9              RColorBrewer_1.1-3      abind_1.4-8             multcomp_1.4-29        
[31] withr_3.0.2             grid_4.5.1              stats4_4.5.1            xtable_1.8-4            inline_0.3.21          
[36] emmeans_1.11.2-8        scales_1.4.0            MASS_7.3-65             cli_3.6.5               mvtnorm_1.3-3          
[41] rmarkdown_2.30          crayon_1.5.3            generics_0.1.4          RcppParallel_5.1.11-1   rstudioapi_0.17.1      
[46] reshape2_1.4.5          tzdb_0.5.0              rstan_2.36.0.9000       splines_4.5.1           bayesplot_1.15.0.9000  
[51] parallel_4.5.1          matrixStats_1.5.0       vctrs_0.7.0             V8_8.0.1                Matrix_1.7-3           
[56] sandwich_3.1-1          jsonlite_2.0.0          hms_1.1.4               bit64_4.6.0-1           glue_1.8.0             
[61] codetools_0.2-20        distributional_0.5.0    stringi_1.8.7           gtable_0.3.6            QuickJSR_1.8.1         
[66] pillar_1.11.1           htmltools_0.5.9         Brobdingnag_1.2-9       R6_2.6.1                vroom_1.6.6            
[71] evaluate_1.0.5          lattice_0.22-7          backports_1.5.0         rstantools_2.5.0.9000   coda_0.19-4.1          
[76] gridExtra_2.3           nlme_3.1-168            checkmate_2.3.3         xfun_0.55               zoo_1.8-14             
[81] pkgconfig_2.0.3        

Comments

Arnold, J. B. (2021). ggthemes: Extra themes, scales and geoms for ’ggplot2’. https://CRAN.R-project.org/package=ggthemes

Bürkner, P.-C. (2022). Estimating multivariate models with brms. https://CRAN.R-project.org/package=brms/vignettes/brms_multivariate.html

Hayes, Andrew F. (2015). An index and test of linear moderated mediation. Multivariate Behavioral Research, 50(1), 1–22. https://doi.org/10.1080/00273171.2014.962683

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

Yuan, Y., & MacKinnon, D. P. (2009). Bayesian mediation analysis. Psychological Methods, 14(4), 301–322. https://doi.org/10.1037/a0016972