10  Multicategorical Focal Antecedents and Moderators

In this chapter, [Hayes] extend[ed] the principles of moderation analysis described in Chapters 7 and 8 to testing interaction involving a multicategorical focal antecedent variable or moderator. As you will see, the principles discussed in those chapters generalize quite readily, although the model necessarily requires more than one product to capture an interaction between two variables. This makes the formulas a bit more complex, and the visualizing and probing process a bit more involved. But with comfort with the fundamentals described so far, you should not find it difficult to master this extension of multiple regression analysis. (Hayes, 2018, p. 350)

10.1 Moderation of the effect of a multicategorical antecedent variable

Take the case of a continuous or dichotomous moderator \(W\) and a multicategorical \(X\) “with \(g\) groups, include \(g − 1\) variables coding membership in the groups, the moderator variable \(W\), and \(g − 1\) products between the \(g − 1\) group codes and moderator \(W\) in a regression model” (p. 351) following the form

\[ Y = i_Y + \sum_{i = 1}^{g - 1} b_i D_i + b_g W + \sum_{j = g + 1}^{2g - 1} b_j D_{j - g} W + e_Y, \]

where \(D_i\) denotes the \(i\)^th dummy variable. Given the case where \(g = 4\), that formula can be re-expressed as

\[\begin{align*} Y & = i_Y + b_1 D_1 + b_2 D_2 + b_3 D_3 + b_4 W + b_5 D_1 W + b_6 D_2 W + b_7 D_3 W + e_Y, \;\;\;\text{or} \\ & = i_Y + (b_1 + b_5 W) D_1 + (b_2 + b_6 W) D_2 + (b_3 + b_7 W) D_3 + b_4 W + e_Y. \end{align*}\]

10.2 An example from the sex disrimination in the workplace study

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, 29, 116, 35, 13, 12…
$ 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, 1, 0, 2, 1, 2, 2, …
$ 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.62, 4.50, 5.50, 4.2…
$ 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, 1, 4, 1, 2, 6, 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.83, 6.16, 4.33, 4.8…
$ 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.25, 5.25, 5.00, 4.5…

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

Load brms.

library(brms)

With model10.1 and model10.2 we fit the multicategorical multivariable model and the multicategorical moderation models, respectively.

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

model10.2 <- update(
  model10.1,
  newdata = protest,
  liking ~ 1 + d1 + d2 + sexism + d1:sexism + d2:sexism,
  cores = 4,
  file = "fits/model10.02")

Behold the \(R^2\) summaries.

r2 <- tibble(`Model 10.1` = bayes_R2(model10.1, summary = F)[, 1],
             `Model 10.2` = bayes_R2(model10.2, summary = F)[, 1]) %>% 
  mutate(`The R2 difference` = `Model 10.2` - `Model 10.1`)
  
r2 %>% 
  pivot_longer(everything())  %>% 
  # this line isn't necessary, but it sets the order the summaries appear in
  mutate(name = factor(name, levels = c("Model 10.1", "Model 10.2", "The R2 difference"))) %>% 
  group_by(name) %>% 
  summarize(mean   = mean(value),
            median = median(value),
            ll     = quantile(value, probs = 0.025),
            ul     = quantile(value, probs = 0.975)) %>% 
  mutate_if(is.double, round, digits = 3)

# A tibble: 3 × 5
  name               mean median     ll    ul
  <fct>             <dbl>  <dbl>  <dbl> <dbl>
1 Model 10.1        0.069  0.065  0.012 0.152
2 Model 10.2        0.155  0.154  0.067 0.252
3 The R2 difference 0.086  0.085 -0.034 0.205

Interestingly, even though our posterior means and medians for the model-specific \(R^2\) values differed some from the OLS estimates in the text, their difference corresponded quite nicely to the one in the text. Let’s take a look at their distributions.

r2 %>% 
  pivot_longer(everything()) %>% 
  ggplot(aes(x = value)) +
  geom_density(linewidth = 0, fill = "grey33") +
  scale_y_continuous(NULL, breaks = NULL) +
  facet_wrap(~ name, scales = "free_y") +
  theme_minimal()

The model coefficient summaries cohere well with those in Table 10.1.

print(model10.1, digits = 3)

 Family: gaussian 
  Links: mu = identity 
Formula: liking ~ 1 + d1 + d2 + sexism 
   Data: protest (Number of observations: 129) 
  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
Intercept    4.763     0.606    3.576    5.927 1.000     4439     3314
d1           0.491     0.226    0.047    0.923 1.000     3765     3088
d2           0.443     0.225   -0.002    0.893 1.001     3959     3022
sexism       0.109     0.116   -0.117    0.333 1.000     4250     3313

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI  Rhat Bulk_ESS Tail_ESS
sigma    1.043     0.067    0.926    1.184 1.000     4641     2727

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

print(model10.2, digits = 3)

 Family: gaussian 
  Links: mu = identity 
Formula: liking ~ d1 + d2 + sexism + d1:sexism + d2:sexism 
   Data: protest (Number of observations: 129) 
  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
Intercept    7.706     1.043    5.685    9.700 1.001     1661     2023
d1          -4.114     1.510   -7.004   -1.176 1.003     1647     2027
d2          -3.467     1.405   -6.242   -0.723 1.001     1572     1861
sexism      -0.472     0.203   -0.862   -0.083 1.001     1644     2106
d1:sexism    0.898     0.289    0.338    1.457 1.003     1636     1988
d2:sexism    0.773     0.274    0.234    1.314 1.001     1541     1937

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI  Rhat Bulk_ESS Tail_ESS
sigma    1.005     0.063    0.893    1.138 1.001     2185     2441

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

10.3 Visualizing the model

To get our version of the values in Table 10.2, we’ll first recreate columns for \(d_1\) through \(W\) (SEXISM) and save then as a tibble, nd.

(
  nd <- tibble(d1 = c(0, 1, 0),
               d2 = c(0, 0, 1)) %>% 
    expand_grid(sexism = quantile(protest$sexism, probs = c(0.16, 0.5, 0.84)))
)

# A tibble: 9 × 3
     d1    d2 sexism
  <dbl> <dbl>  <dbl>
1     0     0   4.31
2     0     0   5.12
3     0     0   5.87
4     1     0   4.31
5     1     0   5.12
6     1     0   5.87
7     0     1   4.31
8     0     1   5.12
9     0     1   5.87

With nd in hand, we’ll feed the predictor values into fitted() for the typical posterior summaries.

fitted(model10.2, newdata = nd) %>% round(digits = 3)

      Estimate Est.Error  Q2.5 Q97.5
 [1,]    5.671     0.222 5.238 6.101
 [2,]    5.287     0.157 4.973 5.581
 [3,]    4.933     0.223 4.476 5.366
 [4,]    5.424     0.247 4.938 5.918
 [5,]    5.769     0.155 5.460 6.080
 [6,]    6.088     0.204 5.675 6.490
 [7,]    5.534     0.203 5.144 5.927
 [8,]    5.778     0.152 5.485 6.077
 [9,]    6.003     0.217 5.579 6.438

The values in our Estimate column correspond to those in the \(\hat Y\) column in the table. We, of course, add summaries of uncertainty to the point estimates.

If we want to make a decent line plot for our version of Figure 10.3, we’ll need many more values for sexism, which will appear on the \(x\)-axis.

nd <- tibble(d1 = c(0, 1, 0),
             d2 = c(0, 0, 1)) %>% 
  expand_grid(sexism = seq(from = 3.5, to = 6.5, length.out = 30))

This time we’ll save the results from fitted() as a tlbble and wrangle a bit to get ready for the plot.

f <- fitted(model10.2, 
            newdata = nd, 
            probs = c(0.025, 0.25, 0.75, 0.975)) %>% 
  data.frame() %>% 
  bind_cols(nd) %>% 
  mutate(condition = if_else(d1 == 1, "Individual Protest",
                             if_else(d2 == 1, "Collective Protest", "No Protest"))) %>% 
  # this line is not necessary, but it will help order the facets of the plot
  mutate(condition = factor(condition, levels = c("No Protest", "Individual Protest", "Collective Protest"))) 

glimpse(f)

Rows: 90
Columns: 10
$ Estimate  <dbl> 6.052223, 6.003355, 5.954486, 5.905618, 5.856749, 5.807881, 5.759013, 5.710144, 5.661276, 5.612407, 5.563539, …
$ Est.Error <dbl> 0.3572760, 0.3385261, 0.3200545, 0.3019124, 0.2841628, 0.2668840, 0.2501737, 0.2341535, 0.2189751, 0.2048256, …
$ Q2.5      <dbl> 5.361808, 5.347695, 5.335067, 5.322182, 5.302876, 5.283948, 5.264961, 5.252240, 5.233252, 5.209949, 5.187643, …
$ Q25       <dbl> 5.808220, 5.773004, 5.737041, 5.702387, 5.664708, 5.627959, 5.589772, 5.553440, 5.514996, 5.475050, 5.437265, …
$ Q75       <dbl> 6.294548, 6.233033, 6.171710, 6.114005, 6.052796, 5.994043, 5.935402, 5.875317, 5.812195, 5.750999, 5.690777, …
$ Q97.5     <dbl> 6.725991, 6.643750, 6.564150, 6.482485, 6.394971, 6.316943, 6.238490, 6.158867, 6.085771, 6.010054, 5.936342, …
$ d1        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 1, 1, 1, 1, 1, 1, 1,…
$ d2        <dbl> 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,…
$ sexism    <dbl> 3.500000, 3.603448, 3.706897, 3.810345, 3.913793, 4.017241, 4.120690, 4.224138, 4.327586, 4.431034, 4.534483, …
$ condition <fct> No Protest, No Protest, No Protest, No Protest, No Protest, No Protest, No Protest, No Protest, No Protest, No…

For Figure 10.3 and many to follow for this chapter, we’ll superimpose 50% intervals on top of 95% intervals.

# this will help us add the original data points to the plot
protest <- protest %>% 
  mutate(condition = ifelse(protest == 0, "No Protest",
                            ifelse(protest == 1, "Individual Protest", 
                                   "Collective Protest"))) %>% 
  mutate(condition = factor(condition, levels = c("No Protest", "Individual Protest", "Collective Protest"))) 
  

# this will help us with the x-axis
breaks <- tibble(values = quantile(protest$sexism, probs = c(0.16, 0.5, 0.84))) %>% 
  mutate(labels = values %>% round(digits = 2) %>% as.character())

# Here we plot
f %>% 
  ggplot(aes(x = sexism)) +
  geom_ribbon(aes(ymin = Q2.5, ymax = Q97.5),
              alpha = 1/3) +
  geom_ribbon(aes(ymin = Q25, ymax = Q75),
              alpha = 1/3) +
  geom_line(aes(y = Estimate)) +
  geom_point(data = protest,
             aes(y = liking),
             size = 2/3) +
  scale_x_continuous(breaks = breaks$values,
                     labels = breaks$labels) +
  coord_cartesian(xlim = c(4, 6),
                  ylim = c(2.5, 7.2)) +
  labs(x = expression(paste("Perceived Pervasiveness of Sex Discrimination in Society (", italic(W), ")")),
       y = "Evaluation of the Attorney") +
  theme_minimal() +
  facet_wrap(~ condition)

By adding the data to the plots, they are both more informative and now serve as a posterior predictive check.

10.4 Probing the interaction

These will involve both omnibus tests and pairwise comparisons.

10.4.1 The pick-a-point approach

“The pick-a-point approach requires you to choose values of the moderator \(W\) and then estimate the conditional effect of \(X\) on \(Y\) at those values and conduct an inferential test” [evaluate the posterior distribution] (p. 368).

10.4.1.1 Omnibus inference

Hayes used the omnibus testing framework to assess how important coefficients \(b_1\) and \(b_2\) were to our interaction model, model1. Before fitting the models, he discussed why he preferred to fit models after centering sexism (i.e., \(W\)) to 4.25. Here we’ll call our centered variable sexism_p, where _p stands in for “prime”.

protest <- protest %>% 
  mutate(sexism_p = sexism - 4.25)

From here on, model10.3 is the moderation model without the lower-order d1 and d2 terms; model10.4 is the full moderation model. But we’re going to be fitting both these models three different ways, based on how we centersexism. So for this first set where we centered sexism on 4.25, we’ll give them the suffix a.

# the model without d1 + d2
model10.3a <- update(
  model10.2,
  newdata = protest,
  liking ~ 1 + sexism_p + d1:sexism_p + d2:sexism_p,
  cores = 4,
  file = "fits/model10.03a")

# the full model with d1 + d2
model10.4a <- update(
  model10.2,
  newdata = protest,
  liking ~ 1 + d1 + d2 + sexism_p + d1:sexism_p + d2:sexism_p,
  cores = 4,
  file = "fits/model10.04a")

The coefficient summaries for model10.4a correspond to the top section of Table 10.3 (p. 373).

fixef(model10.4a) %>% round(digits = 3)

            Estimate Est.Error   Q2.5  Q97.5
Intercept      5.697     0.234  5.249  6.155
d1            -0.297     0.351 -0.981  0.392
d2            -0.184     0.314 -0.786  0.459
sexism_p      -0.472     0.213 -0.916 -0.058
d1:sexism_p    0.897     0.298  0.313  1.486
d2:sexism_p    0.775     0.283  0.226  1.335

We can compare their Bayesian \(R^2\) distributions like we usually do.

library(tidybayes)

r2 <- tibble(`Model without d1 + d2` = bayes_R2(model10.3a, summary = F)[, 1],
             `Model with d1 + d2`    = bayes_R2(model10.4a, summary = F)[, 1]) %>% 
  mutate(`The R2 difference` = `Model with d1 + d2` - `Model without d1 + d2`)
  
r2 %>% 
  pivot_longer(everything()) %>% 
  mutate(name = factor(name, levels = c("Model without d1 + d2", "Model with d1 + d2", "The R2 difference"))) %>% 
  group_by(name) %>% 
  median_qi(value) %>% 
  mutate_if(is.double, round, digits = 3)

# A tibble: 3 × 7
  name                  value .lower .upper .width .point .interval
  <fct>                 <dbl>  <dbl>  <dbl>  <dbl> <chr>  <chr>    
1 Model without d1 + d2 0.14   0.051  0.241   0.95 median qi       
2 Model with d1 + d2    0.154  0.063  0.255   0.95 median qi       
3 The R2 difference     0.013 -0.122  0.153   0.95 median qi       

Our results differ a bit from those in the text (p. 370), but the substantive interpretation is the same. The d1 and d2 parameters added little predictive power to the model in terms of \(R^2\). We can also use information criteria to compare the models. Here are the results from using the LOO-CV.

model10.3a <- add_criterion(model10.3a, criterion = "loo")
model10.4a <- add_criterion(model10.4a, criterion = "loo")

loo_compare(model10.3a, model10.4a) %>% 
  print(simplify = F)

           elpd_diff se_diff elpd_loo se_elpd_loo p_loo  se_p_loo looic  se_looic
model10.3a    0.0       0.0  -186.1     11.1         6.7    2.1    372.2   22.2  
model10.4a   -1.4       0.9  -187.5     10.9         8.0    2.0    375.0   21.8  

The LOO-CV difference between the two models was pretty small. Thus, the LOO-CV gives the same general message as the \(R^2\). The d1 and d2 parameters were sufficiently small and uncertain enough that constraining them to zero did little in terms of reducing the explanatory power of the statistical model.

Here’s the same thing all over again, but this time after centering sexism on 5.120.

protest <- protest %>% 
  mutate(sexism_p = sexism - 5.120)

Now fit the models.

# the model without d1 + d2
model10.3b <- update(
  model10.2,
  newdata = protest,
  liking ~ 1 + sexism_p + d1:sexism_p + d2:sexism_p,
  cores = 4,
  file = "fits/model10.03b")

# the full model with d1 + d2
model10.4b <- update(
  model10.2,
  newdata = protest,
  liking ~ 1 + d1 + d2 + sexism_p + d1:sexism_p + d2:sexism_p,
  cores = 4,
  file = "fits/model10.04b")

These coefficient summaries correspond to the middle section of Table 10.3 (p. 373).

fixef(model10.4b) %>% round(digits = 3)

            Estimate Est.Error   Q2.5  Q97.5
Intercept      5.289     0.156  4.989  5.597
d1             0.480     0.222  0.046  0.904
d2             0.486     0.216  0.060  0.895
sexism_p      -0.471     0.211 -0.890 -0.053
d1:sexism_p    0.903     0.298  0.329  1.486
d2:sexism_p    0.777     0.282  0.213  1.328

Here are the Bayesian \(R^2\) summaries and the summary for their difference.

r2 <- tibble(`Model without d1 + d2` = bayes_R2(model10.3b, summary = F)[, 1],
             `Model with d1 + d2`    = bayes_R2(model10.4b, summary = F)[, 1]) %>%
  mutate(`The R2 difference` = `Model with d1 + d2` - `Model without d1 + d2`)

r2 %>% 
  pivot_longer(everything()) %>% 
  mutate(name = factor(name, levels = c("Model without d1 + d2", "Model with d1 + d2", "The R2 difference"))) %>% 
  group_by(name) %>% 
  median_qi(value) %>% 
  mutate_if(is.double, round, digits = 3)

# A tibble: 3 × 7
  name                  value .lower .upper .width .point .interval
  <fct>                 <dbl>  <dbl>  <dbl>  <dbl> <chr>  <chr>    
1 Model without d1 + d2 0.101  0.028  0.199   0.95 median qi       
2 Model with d1 + d2    0.154  0.064  0.255   0.95 median qi       
3 The R2 difference     0.052 -0.074  0.18    0.95 median qi       

This time, our \(\Delta R^2\) distribution was more similar to the results Hayes reported in the text (p. 370, toward the bottom).

Here’s the updated LOO-CV.

model10.3b <- add_criterion(model10.3b, criterion = "loo")
model10.4b <- add_criterion(model10.4b, criterion = "loo")

loo_compare(model10.3b, model10.4b) %>% 
  print(simplify = F)

           elpd_diff se_diff elpd_loo se_elpd_loo p_loo  se_p_loo looic  se_looic
model10.4b    0.0       0.0  -187.4     10.9         8.0    2.0    374.8   21.8  
model10.3b   -1.2       3.0  -188.6     11.9         6.0    1.9    377.3   23.7  

Here again our Bayesian \(R^2\) and loo() results cohere, both suggesting the d1 and d2 parameters were of little predictive utility. Note how this differs a little from the second \(F\)-test on page 370.

Here’s what happens when we center sexism on 5.896. First center.

protest <- protest %>% 
  mutate(sexism_p = sexism - 5.896)

Fit the models.

# the model without d1 + d2
model10.3c <- update(
  model10.2,
  newdata = protest,
  liking ~ 1 + sexism_p + d1:sexism_p + d2:sexism_p,
  cores = 4,
  file = "fits/model10.03c")

# the full model with d1 + d2
model10.4c <- update(
  model10.2,
  newdata = protest,
  liking ~ 1 + d1 + d2 + sexism_p + d1:sexism_p + d2:sexism_p,
  cores = 4,
  file = "fits/model10.04c")

These coefficient summaries correspond to the lower section of Table 10.3 (p. 373).

fixef(model10.4c) %>% round(digits = 3)

            Estimate Est.Error   Q2.5  Q97.5
Intercept      4.921     0.233  4.464  5.379
d1             1.178     0.315  0.560  1.801
d2             1.096     0.323  0.452  1.729
sexism_p      -0.468     0.205 -0.878 -0.066
d1:sexism_p    0.895     0.292  0.325  1.460
d2:sexism_p    0.774     0.276  0.229  1.299

Again, compute the \(R^2\) distributions and their difference-score distribution.

r2 <- tibble(`Model without d1 + d2` = bayes_R2(model10.3c, summary = F)[, 1],
             `Model with d1 + d2`    = bayes_R2(model10.4c, summary = F)[, 1]) %>%
  mutate(`The R2 difference` = `Model with d1 + d2` - `Model without d1 + d2`)

r2 %>% 
  pivot_longer(everything()) %>% 
  mutate(name = factor(name, levels = c("Model without d1 + d2", "Model with d1 + d2", "The R2 difference"))) %>% 
  group_by(name) %>% 
  median_qi(value) %>% 
  mutate_if(is.double, round, digits = 3)

# A tibble: 3 × 7
  name                  value .lower .upper .width .point .interval
  <fct>                 <dbl>  <dbl>  <dbl>  <dbl> <chr>  <chr>    
1 Model without d1 + d2 0.028  0.003  0.091   0.95 median qi       
2 Model with d1 + d2    0.154  0.064  0.256   0.95 median qi       
3 The R2 difference     0.122  0.015  0.232   0.95 median qi       

That \(\Delta R^2\) distribution matches up nicely with the one Hayes reported at the bottom of page 370. Now compare the models with the LOO.

model10.3c <- add_criterion(model10.3c, criterion = "loo")
model10.4c <- add_criterion(model10.4c, criterion = "loo")

loo_compare(model10.3c, model10.4c) %>% 
  print(simplify = F)

           elpd_diff se_diff elpd_loo se_elpd_loo p_loo  se_p_loo looic  se_looic
model10.4c    0.0       0.0  -187.3     11.0         7.9    2.0    374.7   22.0  
model10.3c   -6.5       5.6  -193.9     13.5         5.9    2.0    387.7   27.0  

Although our Bayesian \(R^2\) difference is now predominantly positive, the LOO-CV difference for the two models remains uncertain. Here’s a look at the two parameters in question using a handmade coefficient plot.

as_draws_df(model10.4c) %>% 
  pivot_longer(b_d1:b_d2) %>% 
  mutate(name = str_remove(name, "b_")) %>% 
  
  ggplot(aes(x = value, y = name)) +
  stat_summary(fun = median,
               fun.min = function(i) quantile(i, probs = 0.025),
               fun.max = function(i) quantile(i, probs = 0.975),
               color = "grey33") +
  stat_summary(geom = "linerange",
               fun.min = function(i) quantile(i, probs = 0.25),
               fun.max = function(i) quantile(i, probs = 0.75),
               color = "grey33",
               linewidth = 1.25) +
  ylab(NULL) +
  coord_cartesian(xlim = c(0, 2)) +
  theme_minimal()

For Figure 10.4, we’ll drop our faceting approach and just make one big plot. Heads up: I’m going to drop the 50% intervals from this plot. They’d just make it too busy.

f %>% 
  ggplot(aes(x = sexism, y = Estimate, ymin = Q2.5, ymax = Q97.5, alpha = condition)) +
  geom_ribbon() +
  geom_line() +
  scale_alpha_manual(values = c(0.2, 0.5, 0.8)) +
  scale_x_continuous(breaks = breaks$values,
                     labels = breaks$labels) +
  coord_cartesian(xlim = c(4, 6),
                  ylim = c(4.5, 6.7)) +
  labs(x = expression("Perceived Pervasiveness of Sex Discrimination in Society "*(italic(W))),
       y = "Evaluation of the Attorney") +
  theme_minimal() +
  theme(legend.direction = "vertical",
        legend.position = "top",
        legend.title = element_blank())

10.4.1.2 Pairwise inference

Hayes continues to reference Table 10.3. In the last subsection, we reproduced those results one model at a time. Why not practice doing it altogether? There are a lot of ways you could do this. A good first try is to extend the fixef() approach from before with a little help from bind_rows().

bind_rows(
  # start with `model4a`
  fixef(model10.4a) %>% 
    data.frame() %>% 
    rownames_to_column("parameter"),
  # add `model4b`
  fixef(model10.4b) %>% 
    data.frame() %>% 
    rownames_to_column("parameter"),
  # add `model4c`
  fixef(model10.4c) %>% 
    data.frame() %>% 
    rownames_to_column("parameter")
) %>% 
  # wrangle a bit
  mutate(`w'` = str_c("w - ", c(4.25, 5.12, 5.896)) %>% rep(., each = 6)) %>% 
  select(`w'`, everything()) %>% 
  mutate_if(is.double, round, digits = 3)

          w'   parameter Estimate Est.Error   Q2.5  Q97.5
1   w - 4.25   Intercept    5.697     0.234  5.249  6.155
2   w - 4.25          d1   -0.297     0.351 -0.981  0.392
3   w - 4.25          d2   -0.184     0.314 -0.786  0.459
4   w - 4.25    sexism_p   -0.472     0.213 -0.916 -0.058
5   w - 4.25 d1:sexism_p    0.897     0.298  0.313  1.486
6   w - 4.25 d2:sexism_p    0.775     0.283  0.226  1.335
7   w - 5.12   Intercept    5.289     0.156  4.989  5.597
8   w - 5.12          d1    0.480     0.222  0.046  0.904
9   w - 5.12          d2    0.486     0.216  0.060  0.895
10  w - 5.12    sexism_p   -0.471     0.211 -0.890 -0.053
11  w - 5.12 d1:sexism_p    0.903     0.298  0.329  1.486
12  w - 5.12 d2:sexism_p    0.777     0.282  0.213  1.328
13 w - 5.896   Intercept    4.921     0.233  4.464  5.379
14 w - 5.896          d1    1.178     0.315  0.560  1.801
15 w - 5.896          d2    1.096     0.323  0.452  1.729
16 w - 5.896    sexism_p   -0.468     0.205 -0.878 -0.066
17 w - 5.896 d1:sexism_p    0.895     0.292  0.325  1.460
18 w - 5.896 d2:sexism_p    0.774     0.276  0.229  1.299

This code works okay, but it’s redundant. Here’s a streamlined approach where we use a combination of nested tibbles and the purrr::map() function to work with our three model fits–model10.4a, model10.4b, and model10.4c–in bulk.

t <- tibble(`w'` = str_c("w - ", c(4.25, 5.12, 5.896)),
            name = str_c("model10.4", letters[1:3])) %>% 
  mutate(fit = map(name, get)) %>% 
  mutate(s = map(fit, ~fixef(.) %>% 
                   data.frame() %>% 
                   rownames_to_column("parameter"))) %>% 
  unnest(s) %>% 
  select(`w'`, parameter:Q97.5)

t %>% 
  mutate_if(is.double, round, digits = 3)

# A tibble: 18 × 6
   `w'`      parameter   Estimate Est.Error   Q2.5  Q97.5
   <chr>     <chr>          <dbl>     <dbl>  <dbl>  <dbl>
 1 w - 4.25  Intercept      5.70      0.234  5.25   6.16 
 2 w - 4.25  d1            -0.297     0.351 -0.981  0.392
 3 w - 4.25  d2            -0.184     0.314 -0.786  0.459
 4 w - 4.25  sexism_p      -0.472     0.213 -0.916 -0.058
 5 w - 4.25  d1:sexism_p    0.897     0.298  0.313  1.49 
 6 w - 4.25  d2:sexism_p    0.775     0.283  0.226  1.34 
 7 w - 5.12  Intercept      5.29      0.156  4.99   5.60 
 8 w - 5.12  d1             0.48      0.222  0.046  0.904
 9 w - 5.12  d2             0.486     0.216  0.06   0.895
10 w - 5.12  sexism_p      -0.471     0.211 -0.89  -0.053
11 w - 5.12  d1:sexism_p    0.903     0.298  0.329  1.49 
12 w - 5.12  d2:sexism_p    0.777     0.282  0.213  1.33 
13 w - 5.896 Intercept      4.92      0.233  4.46   5.38 
14 w - 5.896 d1             1.18      0.315  0.56   1.80 
15 w - 5.896 d2             1.10      0.323  0.452  1.73 
16 w - 5.896 sexism_p      -0.468     0.205 -0.878 -0.066
17 w - 5.896 d1:sexism_p    0.895     0.292  0.325  1.46 
18 w - 5.896 d2:sexism_p    0.774     0.276  0.229  1.30 

Summary tables like this are precise and very common in the literature. But you can get lost in all those numbers. A coefficient plot can be better. This first version is pretty close to the Table 10.3 format.

t %>% 
  # this will help us order our y-axis
  mutate(parameter = factor(parameter,
                            levels = c("d2:sexism_p", "d1:sexism_p", 
                                       "sexism_p", "d2", "d1", "Intercept")),
         # this is just for aesthetics
         `w'` = str_c("w' = ", `w'`)) %>% 
  
  # plot!
  ggplot(aes(x = Estimate, xmin = Q2.5, xmax = Q97.5, y = parameter)) +
  geom_pointrange() +
  labs(x = NULL,
       y = NULL) +
  theme_minimal() +
  theme(axis.text.y = element_text(hjust = 0)) +
  facet_wrap(~ `w'`, nrow = 1)

Notice how this arrangement makes it easiest to compare coefficients within models. If we wanted to make it easier to compare coefficients across models, we might arrange the plot like so.

t %>% 
  # this will help us order our y-axis
  mutate(parameter = factor(parameter,
                            levels = c("Intercept", "d1", "d2", "sexism_p", "d1:sexism_p", "d2:sexism_p"))) %>% 
  
  # plot!
  ggplot(aes(x = Estimate, xmin = Q2.5, xmax = Q97.5, y = `w'`)) +
  geom_pointrange() +
  labs(x = NULL,
       y = NULL) +
  theme_minimal() +
  theme(axis.text.y = element_text(hjust = 0)) +
  facet_wrap(~ parameter, ncol = 1)

Oh man–with sweet plots like these, who needs tables! This makes it much easier to see what happened as we changed values we centered sexism on. In the middle paragraph on page 374, Hayes pointed out “that \(b_1\) and \(b_2\) differ in these analyses, but \(b_3\), \(b_4\), and \(b_5\) are unaffected by the centering”. Our coefficient plot clarified that in a way I don’t think a table ever could. But before we move on, let’s back up a little in the text.

“To make this more concrete, consider the effect of Catherine’s behavior on how she is perceived among people who are relatively high in their perceptions of the pervasiveness of sex discrimination in society” (p. 372). For this, Hayes defined “relatively high” as \(W = 5.896\). To get those estimates for each condition, we’ll use fitted(). Since the number of unique predictor values is small for this example, we’ll just plug them directly into the newdata argument rather than first saving them as a nd object.

fitted(model10.2,
       newdata = tibble(d1 = c(0, 1, 0),
                        d2 = c(0, 0, 1),
                        sexism = 5.896)) %>% 
  round(digits = 3)

     Estimate Est.Error  Q2.5 Q97.5
[1,]    4.920     0.227 4.457 5.359
[2,]    6.099     0.207 5.678 6.508
[3,]    6.011     0.221 5.579 6.451

Those posterior summaries match up nicely with the point estimates Hayes presented at the bottom of page 372. Hayes further expounded:

So by using the regression centering strategy described earlier in the context of an omnibus test of equality of values of \(\hat Y\), the regression coefficients \(b_1\) and \(b_2\) provide pairwise inferences consistent with the coding system used to represent the three groups, conditioned on the value that \(W\) is centered around.

In the next few sentences, he focused on what happened when \(W = 4.250\) (i.e., in model4a). Recall that the two coefficients in question, \(b_1\) and \(b_2\), are named d1 and d2 when we pull their summaries with fixef().

fixef(model10.4a)[c("d1", "d2"), ] %>% 
  round(digits = 3)

   Estimate Est.Error   Q2.5 Q97.5
d1   -0.297     0.351 -0.981 0.392
d2   -0.184     0.314 -0.786 0.459

Hayes then clarified that in this model

\[\begin{align*} b_1 & = \theta_{D_1 \rightarrow Y} \mid (W = 4.250) = 5.400 - 5.698 = -0.299 \;\;\; \text{ and} \\ b_2 & = \theta_{D_2 \rightarrow Y} \mid (W = 4.250) = 5.513 - 5.698 = -0.185. \end{align*}\]

That is, it is the same as a difference score of each of the experimental conditions minus the “No protest” condition. To further show the difference-score quality of these coefficients, we can continue using fitted() in conjunction with the original model10.2 to get the group comparisons for when \(W = 4.250\). Since these involve computing difference scores, we’ll have to use summary = F and do some wrangling.

fitted(model10.2,
       newdata = tibble(d1     = c(0, 1, 0),
                        d2     = c(0, 0, 1),
                        sexism = 4.25),
       summary = F) %>% 
  data.frame() %>% 
  set_names("No Protest", "Individual Protest", "Collective Protest") %>% 
  mutate(difference_a = `Individual Protest` - `No Protest`,
         difference_b = `Collective Protest` - `No Protest`) %>% 
  pivot_longer(everything()) %>% 
  mutate(name = factor(name, levels = c("No Protest", "Individual Protest", "Collective Protest", 
                                        "difference_a", "difference_b"))) %>% 
  group_by(name) %>% 
  mean_qi(value) %>%
  select(name:.upper) %>% 
  mutate_if(is.double, round, digits = 3)

# A tibble: 5 × 4
  name                value .lower .upper
  <fct>               <dbl>  <dbl>  <dbl>
1 No Protest          5.70   5.25   6.14 
2 Individual Protest  5.40   4.89   5.92 
3 Collective Protest  5.52   5.11   5.92 
4 difference_a       -0.299 -0.96   0.396
5 difference_b       -0.182 -0.781  0.45 

Within simulation variance, difference_a is the same as \(b_{1 \mid \text{model10.4a}}\) and difference_b is the same as \(b_{2 \mid \text{model10.4a}}\). Here’s the same thing for when \(W = 5.120\).

fitted(model10.2,
       newdata = tibble(d1     = c(0, 1, 0),
                        d2     = c(0, 0, 1),
                        sexism = 5.120),
       summary = F) %>% 
  data.frame() %>% 
  set_names("No Protest", "Individual Protest", "Collective Protest") %>% 
  mutate(difference_a = `Individual Protest` - `No Protest`,
         difference_b = `Collective Protest` - `No Protest`) %>% 
  pivot_longer(everything()) %>% 
  mutate(name = factor(name, levels = c("No Protest", "Individual Protest", "Collective Protest",
                                        "difference_a", "difference_b"))) %>% 
  group_by(name) %>% 
  mean_qi(value) %>%
  select(name:.upper) %>% 
  mutate_if(is.double, round, digits = 3)

# A tibble: 5 × 4
  name               value .lower .upper
  <fct>              <dbl>  <dbl>  <dbl>
1 No Protest         5.29   4.97   5.58 
2 Individual Protest 5.77   5.46   6.08 
3 Collective Protest 5.78   5.48   6.08 
4 difference_a       0.482  0.049  0.909
5 difference_b       0.491  0.076  0.928

Finally, here it is for when \(W = 5.986\).

fitted(model10.2,
       newdata = tibble(d1     = c(0, 1, 0),
                        d2     = c(0, 0, 1),
                        sexism = 5.986),
       summary = F) %>% 
  data.frame() %>% 
  set_names("No Protest", "Individual Protest", "Collective Protest") %>% 
  mutate(difference_a = `Individual Protest` - `No Protest`,
         difference_b = `Collective Protest` - `No Protest`) %>% 
  pivot_longer(everything()) %>% 
  mutate(name = factor(name, levels = c("No Protest", "Individual Protest", "Collective Protest", 
                                        "difference_a", "difference_b"))) %>% 
  group_by(name) %>% 
  mean_qi(value) %>%
  select(name:.upper) %>% 
  mutate_if(is.double, round, digits = 3)

# A tibble: 5 × 4
  name               value .lower .upper
  <fct>              <dbl>  <dbl>  <dbl>
1 No Protest          4.88  4.39    5.34
2 Individual Protest  6.14  5.69    6.56
3 Collective Protest  6.04  5.58    6.50
4 difference_a        1.26  0.628   1.90
5 difference_b        1.16  0.532   1.82

10.4.2 The Johnson-Neyman technique

As discussed in section 7.4, a problem with the pick-a-point approach to probing an interaction is having to choose values of the moderator. When the moderator is a continuum, you may not have any basis for choosing some values rather than others, and the choice you make will certainly influence the results of the probing exercise to some extent… Actively choosing a different system or con- vention, such as using the sample mean of \(W\), a standard deviation below the mean, and a standard deviation above the mean also does not eliminate the problem. But the Johnson–Neyman (JN) technique avoids this problem entirely. (p. 376)

10.4.2.1 Omnibus inference

Consider the first sentence of the section:

Applied to probing an interaction between a multicategorical \(X\) and a continuous \(W\), an omnibus version of the JM technique involves finding the value or values of \(W\) where their \(F\)-ratio comparing the \(g\) estimated values of \(Y\) is just statistically significant. (p. 376)

Since we’re not using \(F\)-tests with our approach to Bayesian modeling, the closest we might have is a series of \(R^2\) difference tests, which would require refitting the model multiple times over many ways of centering the \(W\)-variable, sexism. I suppose you could do this if you wanted, but it just seems silly, to me. I’ll leave this one up to the interested reader.

10.4.2.2 Pairwise inference

Hayes didn’t make plots for this section, but if you’re careful constructing your nd and with the subsequent wrangling, you can make the usual plots. Since we have two conditions we’d like to compare with No Protest, we’ll make two plots. Here’s the comparison using Individual Protest, first.

# the transition value Hayes identified in the text
Hayes_value <- 5.065

# we need some new data
nd <- tibble(d1 = 0:1,
             d2 = 0) %>% 
  expand_grid(sexism = seq(from = 3.5, to = 6.5, length.out = 30)) %>% 
  mutate(row = 1:n())

# plug those data into `fitted()`
fitted(model10.2,
       newdata = nd,
       summary = F) %>% 
  # wrangle
  data.frame() %>% 
  set_names(pull(nd, row)) %>%
  mutate(draw = 1:n()) %>% 
  pivot_longer(-draw, values_to = "estimate") %>% 
  mutate(row = as.double(name)) %>%
  left_join(nd, by = "row") %>% 
  mutate(condition = if_else(d1 == 0, "No Protest", "Individual Protest")) %>% 
  select(estimate, sexism, draw, condition) %>% 
  pivot_wider(names_from = condition, values_from = estimate) %>% 
  mutate(difference = `Individual Protest` - `No Protest`)%>% 
  
  # plot!
  ggplot(aes(x = sexism, y = difference)) +
  stat_summary(geom = "ribbon",
               fun.min = function(i) quantile(i, probs = 0.025),
               fun.max = function(i) quantile(i, probs = 0.975),
               alpha = 1/3) +
  stat_summary(geom = "ribbon",
               fun.min = function(i) quantile(i, probs = 0.25),
               fun.max = function(i) quantile(i, probs = 0.75),
               alpha = 1/3) +
  stat_summary(geom = "line",
               fun = median) +
  scale_x_continuous(breaks = c(4, Hayes_value, 6),
                     labels = c("4", Hayes_value, "6")) +
  coord_cartesian(xlim = c(4, 6)) +
  labs(subtitle = expression("Our JN-technique plot for "*italic(Individual~Protest)*" compared with "*italic(No~Protest))) +
  theme_minimal()

Now we’re ready to compare No Protest to Collective Protest. The main data difference is which values we assigned to the d1 and d2 columns in nd. For kicks, we should practice another way to get the median line and interval ribbons. The stat_summary() approach from above works great, but it’s verbose. The tidybayes::stat_lineribbon() function will give us the same results with fewer lines of code.

# the transition value Hayes identified in the text
Hayes_value <- 5.036

# new data
nd <- tibble(d1 = 0,
             d2 = 0:1) %>% 
  expand_grid(sexism = seq(from = 3.5, to = 6.5, length.out = 30)) %>% 
  mutate(row = 1:n())

# this part is the same as before
fitted(model10.2,
       newdata = nd,
       summary = F) %>% 
  data.frame() %>% 
  set_names(pull(nd, row)) %>%
  mutate(draw = 1:n()) %>% 
  pivot_longer(-draw, values_to = "estimate") %>% 
  mutate(row = as.double(name)) %>%
  left_join(nd, by = "row") %>% 
  # there are some mild differences in the two mutate() lines
  mutate(condition = if_else(d2 == 0, "No Protest", "Collective Protest")) %>% 
  select(estimate, sexism, draw, condition) %>% 
  pivot_wider(names_from = condition, values_from = estimate) %>% 
  mutate(difference = `Collective Protest` - `No Protest`) %>% 

  # plot!
  ggplot(aes(x = sexism, y = difference)) +
  # look how compact this is!
  stat_lineribbon(.width = c(0.5, 0.95),
                  alpha = 1/3, fill = "black") +
  scale_x_continuous(breaks = c(4, Hayes_value, 6),
                     labels = c("4", Hayes_value, "6")) +
  coord_cartesian(xlim = c(4, 6)) +
  labs(subtitle = expression("Our JN-technique plot for "*italic(Collective~Protest)*" compared with "*italic(No~Protest))) +
  theme_minimal()

And here we do it one last time between the two active protest conditions. For good measure, we will continue experimenting with different ways of plotting the results. This time well first summarize the posterior median and intervals with tidybayes::median_qi() before plotting. We’ll then feed those results into our plot with the aid of tidybayes::geom_lineribbon() and a follow-up scale_fill_manual() line.

nd <- tibble(d1 = 1:0,
             d2 = 0:1) %>%
  expand_grid(sexism = seq(from = 3.5, to = 6.5, length.out = 30)) %>% 
  mutate(row = 1:n())

fitted(model10.2,
       newdata = nd,
       summary = F) %>% 
  data.frame() %>% 
  set_names(pull(nd, row)) %>%
  mutate(draw = 1:n()) %>% 
  pivot_longer(-draw, values_to = "estimate") %>% 
  mutate(row = as.double(name)) %>%
  left_join(nd, by = "row") %>% 
  # there are some mild differences in the two mutate() lines
  mutate(condition = if_else(d1 == 0, "Individual Protest", "Collective Protest")) %>% 
  select(estimate, sexism, draw, condition) %>% 
  pivot_wider(names_from = condition, values_from = estimate) %>% 
  mutate(difference = `Collective Protest` - `Individual Protest`) %>% 
  # group and summarise, here
  group_by(sexism) %>%
  median_qi(difference, .width = c(0.5, 0.95))  %>% 
  
  # plot!
  ggplot(aes(x = sexism, y = difference, ymin = .lower, ymax = .upper)) +
  # look how simple these two lines are
  geom_lineribbon(show.legend = F) +
  scale_fill_manual(values = c("grey75", "grey50")) +
  coord_cartesian(xlim = c(4, 6)) +
  labs(subtitle = expression("Our JN-technique plot for "*italic(Collective~Protest)*" compared with "*italic(Individual~Protest))) +
  theme_minimal()

Little difference between those conditions.

10.5 When the moderator is multicategorical

From a substantive standpoint the combination of

• a multicategorical variable \(X\) and a dichotomous or continuous moderator \(W\) versus
• a dichotomous or continuous variable \(X\) and a multicategorical moderator \(W\)

might seem different. From a modeling perspective, the difference is trivial. As Hayes pointed out, “when we claim from a statistical test of moderation that \(X\)’s effect is moderated by \(W\), then it is also true that \(W\)’s effect is moderated by \(X\). This is the symmetry property of interactions” (p. 381). This symmetry holds when we’re not using the hypothesis-testing framework, too.

10.5.1 An example

Just as a refresher, here’s the print() output for model2.

print(model10.2, digits = 3)

 Family: gaussian 
  Links: mu = identity 
Formula: liking ~ d1 + d2 + sexism + d1:sexism + d2:sexism 
   Data: protest (Number of observations: 129) 
  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
Intercept    7.706     1.043    5.685    9.700 1.001     1661     2023
d1          -4.114     1.510   -7.004   -1.176 1.003     1647     2027
d2          -3.467     1.405   -6.242   -0.723 1.001     1572     1861
sexism      -0.472     0.203   -0.862   -0.083 1.001     1644     2106
d1:sexism    0.898     0.289    0.338    1.457 1.003     1636     1988
d2:sexism    0.773     0.274    0.234    1.314 1.001     1541     1937

Further Distributional Parameters:
      Estimate Est.Error l-95% CI u-95% CI  Rhat Bulk_ESS Tail_ESS
sigma    1.005     0.063    0.893    1.138 1.001     2185     2441

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 Bayesian \(R^2\):

bayes_R2(model10.2) %>% round(digits = 3)

   Estimate Est.Error  Q2.5 Q97.5
R2    0.155     0.049 0.067 0.252

And the \(R^2\) difference between this and the model excluding the interaction terms:

tibble(`Model 10.1` = bayes_R2(model10.1, summary = F)[, 1],
       `Model 10.2` = bayes_R2(model10.2, summary = F)[, 1])  %>% 
  transmute(difference = `Model 10.2` - `Model 10.1`) %>% 
  mean_qi(difference) %>% 
  mutate_if(is.double, round, digits = 3)

# A tibble: 1 × 6
  difference .lower .upper .width .point .interval
       <dbl>  <dbl>  <dbl>  <dbl> <chr>  <chr>    
1      0.086 -0.034  0.205   0.95 mean   qi       

Much like in the text, our Figure 10.7 is just a little different from what we did with Figure 10.3.

# this will help us with the `geom_text()` annotation
slopes <- tibble(
  slope = c(fixef(model10.2)["sexism", "Estimate"] + fixef(model10.2)["d1:sexism", "Estimate"],
            fixef(model10.2)["sexism", "Estimate"] + fixef(model10.2)["d2:sexism", "Estimate"],
            fixef(model10.2)["sexism", "Estimate"]),
  x = c(4.8, 4.6, 5),
  y = c(6.37, 6.25, 4.5),
  condition = factor(c("Individual Protest", "Collective Protest", "No Protest"),
                     levels = c("No Protest", "Individual Protest", "Collective Protest"))) %>% 
  mutate(label = str_c("This slope is about ", slope %>% round(digits = 3)))

# now we plot
f %>% 
  ggplot(aes(x = sexism)) +
  geom_ribbon(aes(ymin = Q2.5, ymax = Q97.5),
              alpha = 1/3) +
  geom_ribbon(aes(ymin = Q25, ymax = Q75),
              alpha = 1/3) +
  geom_line(aes(y = Estimate)) +
  geom_text(data = slopes,
            aes(x = x, y = y, label = label)) +
  coord_cartesian(xlim = c(4, 6)) +
  labs(x = expression(paste("Perceived Pervasiveness of Sex Discrimination in Society (", italic(X), ")")),
       y = "Evaluation of the Attorney") +
  facet_wrap(~ condition) +
  theme_minimal()

10.5.2 Probing the interaction and interpreting the regression coefficients

We computed the posterior means for the slopes when prepping for the figure, above. Here’s how we might get more complete posterior summaries. Much like in the text, our Figure 10.7 is just a little different from what we did with Figure 10.3.

draws <- as_draws_df(model10.2) %>% 
  transmute(`No Protest`         = b_sexism + `b_d1:sexism` * 0 + `b_d2:sexism` * 0,
            `Individual Protest` = b_sexism + `b_d1:sexism` * 1 + `b_d2:sexism` * 0,
            `Collective Protest` = b_sexism + `b_d1:sexism` * 0 + `b_d2:sexism` * 1)

draws %>% 
  pivot_longer(everything()) %>% 
  mutate(name = factor(name, levels = c("No Protest", "Individual Protest", "Collective Protest"))) %>% 
  group_by(name) %>% 
  mean_qi(value) %>% 
  mutate_if(is.double, round, digits = 3)

# A tibble: 3 × 7
  name                value .lower .upper .width .point .interval
  <fct>               <dbl>  <dbl>  <dbl>  <dbl> <chr>  <chr>    
1 No Protest         -0.472 -0.862 -0.083   0.95 mean   qi       
2 Individual Protest  0.425  0.022  0.817   0.95 mean   qi       
3 Collective Protest  0.301 -0.056  0.657   0.95 mean   qi       

Here are the differences among the three protest groups.

draws %>% 
  transmute(`Individual Protest - No Protest`         = `Individual Protest` - `No Protest`,
            `Collective Protest - No Protest`         = `Collective Protest` - `No Protest`,
            `Individual Protest - Collective Protest` = `Individual Protest` - `Collective Protest`) %>% 
  pivot_longer(everything()) %>% 
  # again, not necessary, but useful for reordering the summaries
  mutate(name = factor(name, levels = c("Individual Protest - No Protest", "Collective Protest - No Protest", "Individual Protest - Collective Protest"))) %>% 
  group_by(name) %>% 
  mean_qi(value) %>% 
  mutate_if(is.double, round, digits = 3)

# A tibble: 3 × 7
  name                                    value .lower .upper .width .point .interval
  <fct>                                   <dbl>  <dbl>  <dbl>  <dbl> <chr>  <chr>    
1 Individual Protest - No Protest         0.898  0.338  1.46    0.95 mean   qi       
2 Collective Protest - No Protest         0.773  0.234  1.31    0.95 mean   qi       
3 Individual Protest - Collective Protest 0.125 -0.43   0.665   0.95 mean   qi       

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

Comments

Hayes, A. 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