6 Mediation Analysis with a Multicategorical Antecedent

“Historically, investigators interested in doing a mediation analysis with a multicategorical antecedents \(X\) have resorted to some less optimal strategies than the one [Hayes] discuss[ed] in this chapter” (Andrew F. Hayes, 2018, p. 188). Happily, the approach outlined in this chapter avoids such gaffs. Hayes’s procedure “does not require discarding any data; the entire sample is analyzed simultaneously. Furthermore, the multicategorical nature of \(X\) is respected and retained (p. 189).”

6.1 Relative total, direct, and indirect effects

In review of regression analysis in Chapter 2, we saw that a multicategorical antecedent variable with \(g\) categories can be used as an antecedent variable in a regression model if it is represented by \(g - 1\) variables using some kind of group coding system (see section 2.7). [Hayes] described indicator or dummy coding as one such system, where groups are represented with \(g - 1\) variables set to either zero or one (see Table 2.1). With indicator coding, one of the \(g\) groups is chosen as the reference group. Cases in the reference group receive a zero on all \(g - 1\) variables coding \(X\). Each of the remaining \(g - 1\) groups gets its own indicator variable that is set to 1 for cases in that group, with all other cases set to zero. Using such a system, which of the \(g\) groups a case is in is represented by its pattern of zeros and ones on the \(g - 1\) indicator variables. These \(g - 1\) indicator variables are then used as antecedent variables in a regression model as a stand-in for \(X\). (pp. 189–190, emphasis in the original)

6.1.1 Relative indirect effects.

When our \(X\) is multicategorical, we end up with \(g - 1\) \(a\) coefficients. Presuming the \(M\) variable is continuous or binary, this will yield \(g - 1\) relative indirect effects, \(a_j b\).

6.1.2 Relative direct effects.

Similar to above, when our \(X\) is multicategorical, we end up with \(g - 1\) \(c'\) coefficients, each of which is a relative direct effects.

6.1.3 Relative total effects.

With the two prior subsections in mind, when our \(X\) is multicategorical, we end up with \(g - 1\) \(c\) coefficients, each of which is a relative total effect. These follow the form

\[c_j = c_j' + a_j b,\]

where \(j\) indexes a given group.

6.2 An example: Sex discrimination in the workplace

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

library(tidyverse)

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

glimpse(protest)
## Rows: 129
## Columns: 6
## $ subnum   <dbl> 209, 44, 124, 232, 30, 140, 27, 64, 67, 182, 85, 109, 122, 69, 45, 28, 170, 66, 168, 97, 7,…
## $ protest  <dbl> 2, 0, 2, 2, 2, 1, 2, 0, 0, 0, 2, 2, 0, 1, 1, 0, 1, 2, 2, 1, 2, 1, 1, 2, 2, 0, 1, 1, 0, 1, 1…
## $ sexism   <dbl> 4.87, 4.25, 5.00, 5.50, 5.62, 5.75, 5.12, 6.62, 5.75, 4.62, 4.75, 6.12, 4.87, 5.87, 4.87, 4…
## $ angry    <dbl> 2, 1, 3, 1, 1, 1, 2, 1, 6, 1, 2, 5, 2, 1, 1, 1, 2, 1, 3, 4, 1, 1, 1, 5, 1, 5, 1, 1, 2, 1, 1…
## $ liking   <dbl> 4.83, 4.50, 5.50, 5.66, 6.16, 6.00, 4.66, 6.50, 1.00, 6.83, 5.00, 5.66, 5.83, 6.50, 4.50, 4…
## $ respappr <dbl> 4.25, 5.75, 4.75, 7.00, 6.75, 5.50, 5.00, 6.25, 3.00, 5.75, 5.25, 7.00, 4.50, 6.25, 5.00, 5…

Here are the ungrouped means and \(SD\)s for respappr and liking shown at the bottom of Table 6.1.

protest %>%
  pivot_longer(liking:respappr) %>% 
  group_by(name) %>% 
  summarize(mean = mean(value),
            sd   = sd(value)) %>% 
  mutate_if(is.double, round, digits = 3)
## # A tibble: 2 x 3
##   name      mean    sd
##   <chr>    <dbl> <dbl>
## 1 liking    5.64  1.05
## 2 respappr  4.87  1.35

We compute the summaries for respappr and liking, grouped by protest, like so.

protest %>%
  pivot_longer(liking:respappr) %>% 
  group_by(protest, name) %>% 
  summarize(mean = mean(value),
            sd   = sd(value)) %>% 
  mutate_if(is.double, round, digits = 3)
## # A tibble: 6 x 4
## # Groups:   protest [3]
##   protest name      mean    sd
##     <dbl> <chr>    <dbl> <dbl>
## 1       0 liking    5.31 1.30 
## 2       0 respappr  3.88 1.46 
## 3       1 liking    5.83 0.819
## 4       1 respappr  5.14 1.08 
## 5       2 liking    5.75 0.936
## 6       2 respappr  5.49 0.936

It looks like Hayes has a typo in the \(SD\) for liking when protest == 0. It seems he accidentally entered the value for when protest == 1 in that slot.

You’ll have to wait a minute to see where the adjusted \(Y\) values came from.

With a little if_else(), computing the dummies d1 and d2 is easy enough.

protest <-
  protest %>% 
  mutate(d1 = if_else(protest == 1, 1, 0),
         d2 = if_else(protest == 2, 1, 0))

We’re almost ready to fit the model. Let’s load brms.

library(brms)

This is the first time we’ve had a simple univariate regression model in a while–no special mvbind() syntax or multiple bf() formulas, just straight up brms::brm().

model6.1 <-
  brm(data = protest, 
      family = gaussian,
      liking ~ 1 + d1 + d2,
      cores = 4,
      file = "fits/model06.01")

Check the coefficient summaries.

fixef(model6.1)
##            Estimate Est.Error         Q2.5     Q97.5
## Intercept 5.3105221 0.1646877  4.987615216 5.6267716
## d1        0.5183924 0.2306089  0.071900531 0.9716942
## d2        0.4445530 0.2275941 -0.004987925 0.8885953

Our \(R^2\) differences a bit from the OLS version in the text. This shouldn’t be surprising when it’s near the boundary.

bayes_R2(model6.1)
##      Estimate  Est.Error        Q2.5     Q97.5
## R2 0.05873881 0.03607054 0.006018393 0.1416788

Here’s its shape. For the plots in this chapter, we’ll take a few formatting cues from Edward Tufte (2001), courtesy of the ggthemes package. The theme_tufte() function will change the default font and remove some chart junk. We will take our color palette from Pokemon via the palettetown package (Lucas, 2016).

library(ggthemes)
library(palettetown)

bayes_R2(model6.1, summary = F) %>% 
  data.frame() %>% 
  
  ggplot(aes(x = R2)) +
  geom_density(size = 0, fill = pokepal(pokemon = "plusle")[2]) +
  scale_y_continuous(NULL, breaks = NULL) +
  coord_cartesian(xlim = 0:1) +
  xlab(expression(italic(R)^2)) +
  theme_tufte() +
  theme(legend.title = element_blank(),
        plot.background = element_rect(fill = pokepal(pokemon = "plusle")[8]))

To use the model-implied equations to compute the means for each group on the criterion, we’ll extract the posterior samples.

post <- posterior_samples(model6.1)

post %>% 
  mutate(Y_np = b_Intercept + b_d1 * 0 + b_d2 * 0,
         Y_ip = b_Intercept + b_d1 * 1 + b_d2 * 0,
         Y_cp = b_Intercept + b_d1 * 0 + b_d2 * 1) %>% 
  pivot_longer(contains("Y_")) %>%
  # this line will order our output the same way Hayes did in the text (p. 197)
  mutate(name = factor(name, levels = c("Y_np", "Y_ip", "Y_cp"))) %>% 
  group_by(name) %>% 
  summarize(mean = mean(value),
            sd   = sd(value))
## # A tibble: 3 x 3
##   name   mean    sd
##   <fct> <dbl> <dbl>
## 1 Y_np   5.31 0.165
## 2 Y_ip   5.83 0.159
## 3 Y_cp   5.76 0.154

What Hayes called the “relative total effects” \(c_1\) and \(c_2\) are the d1 and d2 lines in our fixef() output, above.

Here are the sub-models for the mediation model.

m_model <- bf(respappr ~ 1 + d1 + d2)
y_model <- bf(liking   ~ 1 + d1 + d2 + respappr)

There’s a third way to fit multivariate models in brms. It uses either the mvbrmsformula() function, or its abbreviated version, mvbf(). With these, we first define our submodels in br() statements like before. We then combine them within mvbf(), separated with a comma. If we’d like to avoid estimating a residual correlation, which we do in this project–, we then set rescore = FALSE. Here’s how it looks like for our second model.

model6.2 <-
  brm(data = protest, 
      family = gaussian,
      mvbf(m_model, y_model, rescor = FALSE),
      cores = 4,
      file = "fits/model06.02")
print(model6.2)
##  Family: MV(gaussian, gaussian) 
##   Links: mu = identity; sigma = identity
##          mu = identity; sigma = identity 
## Formula: respappr ~ 1 + d1 + d2 
##          liking ~ 1 + d1 + d2 + respappr 
##    Data: protest (Number of observations: 129) 
## 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
## respappr_Intercept     3.89      0.18     3.54     4.24 1.00     5949     3113
## liking_Intercept       3.71      0.30     3.15     4.33 1.00     6332     3358
## respappr_d1            1.26      0.26     0.76     1.76 1.00     6080     3328
## respappr_d2            1.60      0.25     1.11     2.08 1.00     5515     3108
## liking_d1             -0.00      0.22    -0.43     0.42 1.00     4014     3295
## liking_d2             -0.22      0.23    -0.67     0.23 1.00     3916     3034
## liking_respappr        0.41      0.07     0.27     0.55 1.00     4690     3390
## 
## Family Specific Parameters: 
##                Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma_respappr     1.18      0.08     1.04     1.34 1.00     8008     3137
## sigma_liking       0.93      0.06     0.82     1.05 1.00     6155     3168
## 
## 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 Bayesian \(R^2\) posteriors.

bayes_R2(model6.2, summary = F) %>% 
  data.frame() %>% 
  pivot_longer(everything()) %>% 
  
  ggplot(aes(x = value, fill = name)) +
  geom_density(size = 0, alpha = 2/3) +
  annotate("text", x = .2, y = 7, label = "liking", color = pokepal(pokemon = "plusle")[2], family = "Times") +
  annotate("text", x = .355, y = 6, label = "respappr", color = pokepal(pokemon = "plusle")[6], family = "Times") +
  scale_fill_manual(values = pokepal(pokemon = "plusle")[c(2, 6)]) +
  scale_x_continuous(NULL, limits = c(0:1)) +
  scale_y_continuous(NULL, breaks = NULL) +
  labs(title = expression(The~italic(R)^2*" densities overlap near perfectly, both hovering around .25.")) +
  theme_tufte() +
  theme(legend.position = "none",
        plot.background = element_rect(fill = pokepal(pokemon = "plusle")[8]))

To get the model summaries as presented in the second two columns in Table 6.2, we use posterior_samples(), rename a bit, and summarize. Like in the last chapter, here we’ll do so with a little help from tidybayes.

library(tidybayes)

post <-
  posterior_samples(model6.2) %>% 
  mutate(a1       = b_respappr_d1,
         a2       = b_respappr_d2,
         b        = b_liking_respappr,
         c1_prime = b_liking_d1,
         c2_prime = b_liking_d2,
         i_m      = b_respappr_Intercept,
         i_y      = b_liking_Intercept)

post %>% 
  pivot_longer(a1:i_y) %>% 
  group_by(name) %>% 
  mean_qi(value) %>% 
  mutate_if(is.double, round, digits = 3)
## # A tibble: 7 x 7
##   name      value .lower .upper .width .point .interval
##   <chr>     <dbl>  <dbl>  <dbl>  <dbl> <chr>  <chr>    
## 1 a1        1.26   0.759  1.76    0.95 mean   qi       
## 2 a2        1.60   1.11   2.08    0.95 mean   qi       
## 3 b         0.412  0.272  0.548   0.95 mean   qi       
## 4 c1_prime -0.003 -0.43   0.418   0.95 mean   qi       
## 5 c2_prime -0.22  -0.67   0.226   0.95 mean   qi       
## 6 i_m       3.89   3.54   4.24    0.95 mean   qi       
## 7 i_y       3.71   3.15   4.33    0.95 mean   qi

Working with the \(\overline M_{ij}\) formulas in page 199 is quite similar to what we did above.

post %>% 
  mutate(M_np = b_respappr_Intercept + b_respappr_d1 * 0 + b_respappr_d2 * 0,
         M_ip = b_respappr_Intercept + b_respappr_d1 * 1 + b_respappr_d2 * 0,
         M_cp = b_respappr_Intercept + b_respappr_d1 * 0 + b_respappr_d2 * 1) %>% 
  pivot_longer(starts_with("M_")) %>%
  # this line will order our output the same way Hayes did in the text (p. 199)
  mutate(name = factor(name, levels = c("M_np", "M_ip", "M_cp"))) %>% 
  group_by(name) %>% 
  summarize(mean = mean(value),
            sd   = sd(value))
## # A tibble: 3 x 3
##   name   mean    sd
##   <fct> <dbl> <dbl>
## 1 M_np   3.89 0.180
## 2 M_ip   5.15 0.181
## 3 M_cp   5.49 0.177

The \(\overline Y^*_{ij}\) formulas are more of the same.

post <-
  post %>% 
  mutate(Y_np = b_liking_Intercept + b_liking_d1 * 0 + b_liking_d2 * 0 + b_liking_respappr * mean(protest$respappr),
         Y_ip = b_liking_Intercept + b_liking_d1 * 1 + b_liking_d2 * 0 + b_liking_respappr * mean(protest$respappr),
         Y_cp = b_liking_Intercept + b_liking_d1 * 0 + b_liking_d2 * 1 + b_liking_respappr * mean(protest$respappr))

post %>% 
  pivot_longer(starts_with("Y_")) %>% 
  mutate(name = factor(name, levels = c("Y_np", "Y_ip", "Y_cp"))) %>% 
  group_by(name) %>% 
  summarize(mean = mean(value),
            sd   = sd(value))
## # A tibble: 3 x 3
##   name   mean    sd
##   <fct> <dbl> <dbl>
## 1 Y_np   5.72 0.161
## 2 Y_ip   5.71 0.141
## 3 Y_cp   5.50 0.144

Note, these are where the adjusted \(Y\) values came from in Table 6.1.

This is as fine a spot as any to introduce coefficient plots. The brms, tidybayes, and bayesplot packages all offer convenience functions for coefficient plots. Before we get all lazy using convenience functions, it’s good to know how to make coefficient plots by hand. Here’s ours for those last three \(\overline Y^*_{ij}\)-values.

post %>% 
  pivot_longer(starts_with("Y_")) %>% 
  
  ggplot(aes(x = value, y = name, color = name)) +
  stat_summary(geom = "pointrange",
               fun = median,
               fun.min = function(x) {quantile(x, probs = .025)},
               fun.max = function(x) {quantile(x, probs = .975)},
               size = .75) +
  stat_summary(geom = "linerange",
               fun.min = function(x) {quantile(x, probs = .25)},
               fun.max = function(x) {quantile(x, probs = .75)},
               size = 1.5) +
  scale_color_manual(values = pokepal(pokemon = "plusle")[c(3, 7, 9)]) +
  labs(x = NULL, y = NULL) +
  theme_tufte() +
  theme(axis.ticks.y = element_blank(),
        legend.position = "none",
        plot.background = element_rect(fill = pokepal(pokemon = "plusle")[8]))

The points are the posterior medians, the thick inner lines the 50% intervals, and the thinner outer lines the 95% intervals. For kicks, we distinguished the three values by color.

If we want to examine \(R^2\) change for dropping the dummy variables, we’ll first fit a model that omits them.

model6.3 <-
  brm(data = protest, 
      family = gaussian,
      liking ~ 1 + respappr,
      cores = 4,
      file = "fits/model06.03")

Here are the competing \(R^2\) distributions.

# get the R2 draws and wrangle
r2 <-
  rbind(bayes_R2(model6.2, resp = "liking", summary = F),
        bayes_R2(model6.3, summary = F)) %>% 
  data.frame() %>% 
  set_names("R2") %>% 
  mutate(fit = rep(c("model6.2", "model6.3"), each = n() / 2))

# plot!
r2 %>% 
  ggplot(aes(x = R2, fill = fit)) +
  geom_density(size = 0, alpha = 2/3) +
  annotate("text", x = .18, y = 6.75, label = "model3", color = pokepal(pokemon = "plusle")[7], family = "Times") +
  annotate("text", x = .35, y = 6.75, label = "model2", color = pokepal(pokemon = "plusle")[6], family = "Times") +
  scale_fill_manual(values = pokepal(pokemon = "plusle")[c(6, 7)]) +
  scale_x_continuous(NULL, limits = 0:1) +
  scale_y_continuous(NULL, breaks = NULL) +
  ggtitle(expression(The~italic(R)^2*" densities for LIKING overlap a lot.")) +
  theme_tufte() +
  theme(legend.position = "none",
        plot.background = element_rect(fill = pokepal(pokemon = "plusle")[8]))

If you want to compare then with a change score, do something like this.

r2 %>%
  mutate(iter = rep(1:4000, times = 2)) %>% 
  pivot_wider(names_from = fit, values_from = R2) %>% 

  ggplot(aes(x = model6.2 - model6.3)) +
  geom_density(size = 0, fill = pokepal(pokemon = "plusle")[4]) +
  geom_vline(xintercept = 0, color = pokepal(pokemon = "plusle")[8]) +
  scale_y_continuous(NULL, breaks = NULL) +
  labs(title = expression(The~Delta*italic(R)^2~distribution),
       subtitle = "Doesn't appear we have a lot of change.",
       x = NULL) +
  theme_tufte() +
  theme(legend.title = element_blank(),
        plot.background = element_rect(fill = pokepal(pokemon = "plusle")[8]))

Now compute the posterior means and 95% intervals for \(a_1 b\) and \(a_2 b\), the conditional indirect effects.

post %>% 
  mutate(a1b = a1 * b,
         a2b = a2 * b) %>%
  pivot_longer(a1b:a2b) %>% 
  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 a1b   0.517  0.27   0.807   0.95 mean   qi       
## 2 a2b   0.659  0.381  0.982   0.95 mean   qi

6.3 Using a different group coding system

Here we’ll make our alternative dummies, what we’ll call d_1 and d_2, with orthogonal contrast coding.

protest <-
  protest %>% 
  mutate(d_1 = if_else(protest == 0, -2/3, 1/3),
         d_2 = if_else(protest == 0, 0, 
                       if_else(protest == 1, -1/2, 1/2)))

Here are the sub-models.

m_model <- bf(respappr ~ 1 + d_1 + d_2)
y_model <- bf(liking   ~ 1 + d_1 + d_2 + respappr)

Now we fit using the mvbf() approach.

model6.4 <-
  brm(data = protest, 
      family = gaussian,
      mvbf(m_model, y_model, rescor = FALSE),
      cores = 4,
      file = "fits/model06.04")

Here are our intercepts and regression coefficient summaries.

fixef(model6.4)
##                      Estimate  Est.Error       Q2.5     Q97.5
## respappr_Intercept  4.8408491 0.10400545  4.6391636 5.0524578
## liking_Intercept    3.6363534 0.35318446  2.9437136 4.3215543
## respappr_d_1        1.4376291 0.21964377  1.0104644 1.8696691
## respappr_d_2        0.3478888 0.24998359 -0.1359316 0.8558661
## liking_d_1         -0.1120374 0.20567008 -0.5092694 0.2919660
## liking_d_2         -0.2193080 0.19541186 -0.6081167 0.1656851
## liking_respappr     0.4117305 0.07136229  0.2730765 0.5551897

It’s important to note that these will not correspond to the “TOTAL EFFECT MODEL” section of the PROCESS output of Figure 6.3. Hayes’s PROCESS has the mcx=3 command which tells the program to reparametrize the orthogonal contrasts. brms doesn’t have such a command.

For now, we’ll have to jump to Equation 6.8 towards the bottom of page 207. Those parameters are evident in our output. For good measure, here we’ll practice with posterior_summary().

posterior_summary(model6.4) %>% 
  data.frame() %>% 
  rownames_to_column("parameter") %>% 
  filter(str_detect(parameter, "b_respappr"))
##              parameter  Estimate Est.Error       Q2.5     Q97.5
## 1 b_respappr_Intercept 4.8408491 0.1040055  4.6391636 5.0524578
## 2       b_respappr_d_1 1.4376291 0.2196438  1.0104644 1.8696691
## 3       b_respappr_d_2 0.3478888 0.2499836 -0.1359316 0.8558661

Thus it’s easy to get the \(\overline M_{ij}\) means with a little posterior manipulation.

post <-
  posterior_samples(model6.4) %>% 
  mutate(M_np = b_respappr_Intercept + b_respappr_d_1 * -2/3 + b_respappr_d_2 *    0,
         M_ip = b_respappr_Intercept + b_respappr_d_1 *  1/3 + b_respappr_d_2 * -1/2,
         M_cp = b_respappr_Intercept + b_respappr_d_1 *  1/3 + b_respappr_d_2 *  1/2)

post %>% 
  pivot_longer(starts_with("M_")) %>% 
  mutate(name = factor(name, levels = c("M_np", "M_ip", "M_cp"))) %>% 
  group_by(name) %>% 
  summarize(mean = mean(value),
            sd   = sd(value))
## # A tibble: 3 x 3
##   name   mean    sd
##   <fct> <dbl> <dbl>
## 1 M_np   3.88 0.182
## 2 M_ip   5.15 0.181
## 3 M_cp   5.49 0.173

With these in hand, we can compute \(a_1\) and \(a_2\).

post <-
  post %>% 
  mutate(a1 = (M_ip + M_cp) / 2 - M_np,
         a2 = M_cp - M_ip)

post %>% 
  pivot_longer(a1:a2) %>% 
  group_by(name) %>% 
  summarize(mean = mean(value),
            sd   = sd(value))
## # A tibble: 2 x 3
##   name   mean    sd
##   <chr> <dbl> <dbl>
## 1 a1    1.44  0.220
## 2 a2    0.348 0.250

Happily, our model output will allow us to work with Hayes’s \(\overline Y^*_{ij}\) equations in the middle of page 210.

post <-
  post %>% 
  mutate(Y_np = b_liking_Intercept + b_liking_d_1 * -2/3 + b_liking_d_2 *    0 + b_liking_respappr * mean(protest$respappr),
         Y_ip = b_liking_Intercept + b_liking_d_1 *  1/3 + b_liking_d_2 * -1/2 + b_liking_respappr * mean(protest$respappr),
         Y_cp = b_liking_Intercept + b_liking_d_1 *  1/3 + b_liking_d_2 *  1/2 + b_liking_respappr * mean(protest$respappr))

post %>% 
  pivot_longer(starts_with("Y_")) %>% 
  mutate(name = factor(name, levels = c("Y_np", "Y_ip", "Y_cp"))) %>% 
  group_by(name) %>%
  summarize(mean = mean(value),
            sd   = sd(value))
## # A tibble: 3 x 3
##   name   mean    sd
##   <fct> <dbl> <dbl>
## 1 Y_np   5.71 0.164
## 2 Y_ip   5.71 0.141
## 3 Y_cp   5.49 0.144

And with these in hand, we can compute \(c'_1\) and \(c'_2\).

post <-
  post %>% 
  mutate(c1_prime = (Y_ip + Y_cp) / 2 - Y_np,
         c2_prime = Y_cp - Y_ip)

post %>% 
  pivot_longer(c1_prime:c2_prime) %>% 
  group_by(name) %>% 
  summarize(mean = mean(value),
            sd   = sd(value))
## # A tibble: 2 x 3
##   name       mean    sd
##   <chr>     <dbl> <dbl>
## 1 c1_prime -0.112 0.206
## 2 c2_prime -0.219 0.195

It appears Hayes has a typo in the formula for \(c'_2\) on page 211. The value he has down for \(\overline Y^*_{IP}\), 5.145, is incorrect. It’s not the one he displayed at the bottom of the previous page and it also contradicts the analyses herein. So it goes… These things happen.

We haven’t spelled it out, but the \(b\) parameter is currently labeled b_liking_respappr in our post object. Here we’ll make a b column to make things easier. While we’re at it, we’ll compute the indirect effects, too.

post <-
  post %>%
  mutate(b = b_liking_respappr) %>% 
  mutate(a1b = a1 * b,
         a2b = a2 * b)

post %>% 
  pivot_longer(a1b:a2b) %>% 
  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 a1b   0.591  0.351  0.875   0.95 mean   qi       
## 2 a2b   0.143 -0.051  0.368   0.95 mean   qi

Now we can compute and summarize() our \(c_1\) and \(c_2\).

post <-
  post %>% 
  mutate(c1 = c1_prime + a1b,
         c2 = c2_prime + a2b)

post %>% 
  pivot_longer(c1:c2) %>% 
  group_by(name) %>% 
  summarize(mean = mean(value),
            sd   = sd(value))
## # A tibble: 2 x 3
##   name     mean    sd
##   <chr>   <dbl> <dbl>
## 1 c1     0.479  0.201
## 2 c2    -0.0761 0.217

6.4 Some miscellaneous issues [unrelated to those Hayes covered in the text]

Do you recall how way back in Chapter 2 we covered an alternative way to fit models with multicategorical grouping variables? Well, we did. The basic strategy is to save our grouping variable as a factor and then enter it into the model with the special 0 + syntax, which removes the typical intercept. Since this chapter is all about multicategorical variables, it might make sense to explore what happens when we use this approach. For our first step, well prepare the data.

protest <-
  protest %>% 
  mutate(group = factor(protest,
                        levels = 0:2,
                        labels = c("none", "individual", "collective")))

protest %>% 
  select(protest, group)
## # A tibble: 129 x 2
##    protest group     
##      <dbl> <fct>     
##  1       2 collective
##  2       0 none      
##  3       2 collective
##  4       2 collective
##  5       2 collective
##  6       1 individual
##  7       2 collective
##  8       0 none      
##  9       0 none      
## 10       0 none      
## # … with 119 more rows

Before we fit a full mediation model, we should warm up. Here we fit a univariable model for liking. This is an alternative to what we did way back with model6.1.

model6.5 <-
  brm(data = protest, 
      family = gaussian,
      liking ~ 0 + group,
      cores = 4,
      file = "fits/model06.05")

Check the summary.

print(model6.5)
##  Family: gaussian 
##   Links: mu = identity; sigma = identity 
## Formula: liking ~ 0 + group 
##    Data: protest (Number of observations: 129) 
## 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
## groupnone           5.31      0.16     4.99     5.63 1.00     4289     2749
## groupindividual     5.82      0.16     5.51     6.14 1.00     4466     2830
## groupcollective     5.75      0.15     5.46     6.04 1.00     4168     3341
## 
## Family Specific Parameters: 
##       Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma     1.04      0.07     0.92     1.18 1.00     5139     3391
## 
## 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).

There’s no conventional intercept parameter. Rather, each of the each of the levels of group get its own conditional intercept. To get a sense of what this model is, let’s practice our coefficient plotting skills. This time we’ll compute the necessary values before plugging them into ggplot2.

# compute the means for `liking` by `group`
group_means <-
  protest %>% 
  group_by(group) %>% 
  summarize(mu_liking = mean(liking))

# pull the posterior summaries and wrangle
fixef(model6.5) %>% 
  data.frame() %>% 
  rownames_to_column("parameter") %>% 
  mutate(group = str_remove(parameter, "group")) %>% 
  
  # plot!
  ggplot(aes(y = group)) +
  # this is the main function for our coefficient plots
  geom_pointrange(aes(x = Estimate, xmin = Q2.5, xmax = Q97.5, color = group),
                  size = 3/4) +
  geom_point(data = group_means,
             aes(x = mu_liking)) +
  scale_color_manual(values = pokepal(pokemon = "plusle")[c(3, 7, 9)]) +
  labs(x = NULL, y = NULL) +
  theme_tufte() +
  theme(axis.text.y = element_text(hjust = 0),
        axis.ticks.y = element_blank(),
        legend.position = "none",
        plot.background = element_rect(fill = pokepal(pokemon = "plusle")[8]))

The results from the model are in colored point ranges. The black dots in the foreground are the empirical means. It looks like our model did a good job estimating the group means for liking.

Let’s see how this coding approach works when you fit a full mediation model. First, define the sub-models with two bf() lines.

m_model <- bf(respappr ~ 0 + group)
y_model <- bf(liking   ~ 0 + group + respappr)

Now fit model6 using the mvbf() approach.

model6.6 <-
  brm(data = protest, 
      family = gaussian,
      mvbf(m_model, y_model, rescor = FALSE),
      cores = 4,
      file = "fits/model06.06")

What will the summary hold?

print(model6.6)
##  Family: MV(gaussian, gaussian) 
##   Links: mu = identity; sigma = identity
##          mu = identity; sigma = identity 
## Formula: respappr ~ 0 + group 
##          liking ~ 0 + group + respappr 
##    Data: protest (Number of observations: 129) 
## 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
## respappr_groupnone           3.89      0.18     3.52     4.23 1.00     4039     2601
## respappr_groupindividual     5.14      0.18     4.77     5.51 1.00     3599     2467
## respappr_groupcollective     5.50      0.17     5.15     5.84 1.00     3960     2741
## liking_groupnone             3.71      0.31     3.11     4.34 1.00     1673     2291
## liking_groupindividual       3.70      0.39     2.93     4.48 1.00     1658     1967
## liking_groupcollective       3.49      0.41     2.70     4.28 1.00     1687     2011
## liking_respappr              0.41      0.07     0.27     0.55 1.00     1590     1950
## 
## Family Specific Parameters: 
##                Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma_respappr     1.18      0.07     1.04     1.34 1.00     3862     2463
## sigma_liking       0.93      0.06     0.82     1.05 1.00     4106     2867
## 
## 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).

If you flip back to page 199, you’ll notice the posterior mean in the first three rows (i.e., respappr_groupnone through respappr_groupcollective) correspond to the estimates for \(\overline M_{NP}\), \(\overline M_{IP}\), and \(\overline M_{CP}\), respectively.

Let’s get to the \(a_j b\) estimates.

post <- 
  posterior_samples(model6.6) %>% 
  mutate(a1 = b_respappr_groupnone,
         a2 = b_respappr_groupindividual,
         a3 = b_respappr_groupcollective,
         b  = b_liking_respappr) %>% 
  mutate(a1b = a1 * b,
         a2b = a2 * b,
         a3b = a3 * b)

post %>% 
  pivot_longer(a1b:a3b) %>% 
  group_by(name) %>% 
  mean_qi(value) %>% 
  mutate_if(is.double, round, digits = 3)
## # A tibble: 3 x 7
##   name  value .lower .upper .width .point .interval
##   <chr> <dbl>  <dbl>  <dbl>  <dbl> <chr>  <chr>    
## 1 a1b    1.6    1.05   2.15   0.95 mean   qi       
## 2 a2b    2.12   1.39   2.85   0.95 mean   qi       
## 3 a3b    2.26   1.49   3.04   0.95 mean   qi

With this parameterization, it’s a little difficult to interpret the \(a_j b\) estimates. None of them are in comparison to anything. However, this approach is quite useful once you compute their various difference scores.

# compute the difference scores
post <- 
  post %>% 
  mutate(diff_individual_minus_none       = a2b - a1b,
         diff_collective_minus_none       = a3b - a1b,
         diff_collective_minus_individual = a3b - a2b)

# summarize
post %>% 
  pivot_longer(contains("diff")) %>% 
  group_by(name) %>% 
  mean_qi(value) %>% 
  mutate_if(is.double, round, digits = 3)
## # A tibble: 3 x 7
##   name                             value .lower .upper .width .point .interval
##   <chr>                            <dbl>  <dbl>  <dbl>  <dbl> <chr>  <chr>    
## 1 diff_collective_minus_individual 0.147 -0.053  0.364   0.95 mean   qi       
## 2 diff_collective_minus_none       0.663  0.378  0.994   0.95 mean   qi       
## 3 diff_individual_minus_none       0.517  0.272  0.82    0.95 mean   qi

If you look back to our results from model6.2, you’ll see that diff_individual_minus_none and diff_collective_minus_none correspond to a1b and a2b, respectively. But with our model6.6 approach, we get the additional information of what kind of indirect effect we might have yielded had we used a different coding scheme for our original set of dummy variables. That is, we get diff_collective_minus_individual.

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:
## [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_3.0.0   palettetown_0.1.1 ggthemes_4.2.4    brms_2.15.0       Rcpp_1.0.6       
##  [6] forcats_0.5.1     stringr_1.4.0     dplyr_1.0.6       purrr_0.3.4       readr_1.4.0      
## [11] tidyr_1.1.3       tibble_3.1.2      ggplot2_3.3.5     tidyverse_1.3.1  
## 
## loaded via a namespace (and not attached):
##   [1] readxl_1.3.1         backports_1.2.1      plyr_1.8.6           igraph_1.2.6         svUnit_1.0.3        
##   [6] splines_4.0.4        crosstalk_1.1.0.1    TH.data_1.0-10       rstantools_2.1.1     inline_0.3.17       
##  [11] digest_0.6.27        htmltools_0.5.1.1    rsconnect_0.8.16     fansi_0.4.2          checkmate_2.0.0     
##  [16] magrittr_2.0.1       modelr_0.1.8         RcppParallel_5.0.2   matrixStats_0.57.0   xts_0.12.1          
##  [21] sandwich_3.0-0       prettyunits_1.1.1    colorspace_2.0-0     rvest_1.0.1          ggdist_3.0.0        
##  [26] haven_2.3.1          xfun_0.23            callr_3.7.0          crayon_1.4.1         jsonlite_1.7.2      
##  [31] lme4_1.1-25          survival_3.2-10      zoo_1.8-8            glue_1.4.2           gtable_0.3.0        
##  [36] emmeans_1.5.2-1      V8_3.4.0             distributional_0.2.2 pkgbuild_1.2.0       rstan_2.21.2        
##  [41] abind_1.4-5          scales_1.1.1         mvtnorm_1.1-1        DBI_1.1.0            miniUI_0.1.1.1      
##  [46] xtable_1.8-4         stats4_4.0.4         StanHeaders_2.21.0-7 DT_0.16              htmlwidgets_1.5.3   
##  [51] httr_1.4.2           threejs_0.3.3        arrayhelpers_1.1-0   posterior_1.0.1      ellipsis_0.3.2      
##  [56] farver_2.1.0         pkgconfig_2.0.3      loo_2.4.1            sass_0.3.1           dbplyr_2.1.1        
##  [61] utf8_1.2.1           tidyselect_1.1.1     labeling_0.4.2       rlang_0.4.11         reshape2_1.4.4      
##  [66] later_1.2.0          munsell_0.5.0        cellranger_1.1.0     tools_4.0.4          cli_3.0.1           
##  [71] generics_0.1.0       broom_0.7.6          ggridges_0.5.3       evaluate_0.14        fastmap_1.1.0       
##  [76] processx_3.5.2       knitr_1.33           fs_1.5.0             nlme_3.1-152         mime_0.10           
##  [81] projpred_2.0.2       xml2_1.3.2           compiler_4.0.4       bayesplot_1.8.0      shinythemes_1.1.2   
##  [86] rstudioapi_0.13      gamm4_0.2-6          curl_4.3             reprex_2.0.0         statmod_1.4.35      
##  [91] bslib_0.2.4          stringi_1.6.2        highr_0.9            ps_1.6.0             Brobdingnag_1.2-6   
##  [96] lattice_0.20-41      Matrix_1.3-2         nloptr_1.2.2.2       markdown_1.1         tensorA_0.36.2      
## [101] shinyjs_2.0.0        vctrs_0.3.8          pillar_1.6.1         lifecycle_1.0.0      jquerylib_0.1.4     
## [106] bridgesampling_1.0-0 estimability_1.3     httpuv_1.6.0         R6_2.5.0             bookdown_0.22       
## [111] promises_1.2.0.1     gridExtra_2.3        codetools_0.2-18     boot_1.3-26          colourpicker_1.1.0  
## [116] MASS_7.3-53          gtools_3.8.2         assertthat_0.2.1     withr_2.4.2          shinystan_2.5.0     
## [121] multcomp_1.4-16      mgcv_1.8-33          parallel_4.0.4       hms_1.1.0            grid_4.0.4          
## [126] coda_0.19-4          minqa_1.2.4          rmarkdown_2.8        shiny_1.6.0          lubridate_1.7.10    
## [131] 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
Lucas, T. (2016). palettetown: Use Pokemon inspired colour palettes [Manual]. https://CRAN.R-project.org/package=palettetown
Tufte, E. R. (2001). The visual display of quantitative information (Second Edition). Graphics Press. https://www.edwardtufte.com/tufte/books_vdqi