06.Rmd

# Mediation Analysis with a Multicategorical Antecedent

```{r, echo = FALSE, cache = FALSE}
options(width = 110)
```

"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" [@hayesIntroductionMediationModeration2018, 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)."

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

### 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$.

### 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*.

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

## An example: Sex discrimination in the workplace

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

```{r, warning = F, message = F}
library(tidyverse)

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

glimpse(protest)
```

Here are the ungrouped means and $\textit{SD}$'s for `respappr` and `liking` shown at the bottom of Table 6.1.

```{r, message = F}
protest %>%
  pivot_longer(liking:respappr) %>% 
  group_by(name) %>% 
  summarize(mean = mean(value),
            sd   = sd(value)) %>% 
  mutate_if(is.double, round, digits = 3)
```

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

```{r, warning = F, message = F}
protest %>%
  pivot_longer(liking:respappr) %>% 
  group_by(protest, name) %>% 
  summarize(mean = mean(value),
            sd   = sd(value)) %>% 
  mutate_if(is.double, round, digits = 3)
```
 
It looks like Hayes has a typo in the $\textit{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.

```{r}
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**.

```{r, message = F, warning = F}
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()`.

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

Check the coefficient summaries.

```{r}
fixef(model6.1)
```

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

```{r}
bayes_R2(model6.1)
```

Here's its shape. For the plots in this chapter, we'll take a few formatting cues from Edward Tufte [-@tufteVisualDisplayQuantitative2001], courtesy of the [**ggthemes** package](https://CRAN.R-project.org/package=ggthemes). The `theme_tufte()` function will change the default font and remove some chart junk. We will take our color palette from [Pokemon](https://pokemonbyreview.blogspot.com/2017/02/311-312-plusle-minun.html) via the [**palettetown** package](https://CRAN.R-project.org/package=palettetown) [@R-palettetown].

```{r, fig.width = 4, fig.height = 2.5}
library(ggthemes)
library(palettetown)

bayes_R2(model6.1, summary = F) %>% 
  data.frame() %>% 
  
  ggplot(aes(x = R2)) +
  geom_density(linewidth = 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 draws.

```{r, warning = F}
draws <- as_draws_df(model6.1)

draws %>% 
  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))
```

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.

```{r}
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.

```{r model6.2}
model6.2 <- brm(
  data = protest, 
  family = gaussian,
  mvbf(m_model, y_model, rescor = FALSE),
  cores = 4,
  file = "fits/model06.02")
```

```{r}
print(model6.2)
```

Behold the Bayesian $R^2$ posteriors.

```{r, fig.width = 6, fig.height = 2.5}
bayes_R2(model6.2, summary = F) %>% 
  data.frame() %>% 
  pivot_longer(everything()) %>% 
  
  ggplot(aes(x = value, fill = name)) +
  geom_density(linewidth = 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 `as_draws_df()`, rename a bit, and summarize. Like in the last chapter, here we'll do so with a little help from **tidybayes**.

```{r, warning = F, message = F}
library(tidybayes)

draws <- as_draws_df(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)

draws %>% 
  pivot_longer(a1:i_y) %>% 
  group_by(name) %>% 
  mean_qi(value) %>% 
  mutate_if(is.double, round, digits = 3)
```

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

```{r, warning = F}
draws %>% 
  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))
```

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

```{r, warning = F}
draws <- draws %>% 
  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))

draws %>% 
  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))
```

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.

```{r, fig.width = 6, fig.height = 1, warning = F}
draws %>% 
  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 = 0.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.

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

Here are the competing $R^2$ distributions.

```{r, fig.width = 6, fig.height = 2.5}
# 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(linewidth = 0, alpha = 2/3) +
  annotate("text", x = .15, y = 6.4, label = "model3", color = pokepal(pokemon = "plusle")[7], family = "Times") +
  annotate("text", x = .33, y = 6.9, 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.

```{r, fig.width = 4, fig.height = 2.5}
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(linewidth = 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.

```{r, warning = F}
draws %>% 
  mutate(a1b = a1 * b,
         a2b = a2 * b) %>%
  pivot_longer(a1b:a2b) %>% 
  group_by(name) %>% 
  mean_qi(value) %>% 
  mutate_if(is.double, round, digits = 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.

```{r}
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.

```{r}
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.

```{r model6.4}
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.

```{r}
fixef(model6.4)
```

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

```{r}
posterior_summary(model6.4) %>% 
  data.frame() %>% 
  rownames_to_column("parameter") %>% 
  filter(str_detect(parameter, "b_respappr"))
```

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

```{r, warning = F}
draws <- as_draws_df(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)

draws %>% 
  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))
```

With these in hand, we can compute $a_1$ and $a_2$.

```{r, warning = F}
draws <- draws %>% 
  mutate(a1 = (M_ip + M_cp) / 2 - M_np,
         a2 = M_cp - M_ip)

draws %>% 
  pivot_longer(a1:a2) %>% 
  group_by(name) %>% 
  summarize(mean = mean(value),
            sd   = sd(value))
```

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

```{r, warning = F}
draws <- draws %>% 
  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))

draws %>% 
  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))
```

And with these in hand, we can compute $c'_1$ and $c'_2$.

```{r, warning = F}
draws <- draws %>% 
  mutate(c1_prime = (Y_ip + Y_cp) / 2 - Y_np,
         c2_prime = Y_cp - Y_ip)

draws %>% 
  pivot_longer(c1_prime:c2_prime) %>% 
  group_by(name) %>% 
  summarize(mean = mean(value),
            sd   = sd(value))
```

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 `draws` object. Here we'll make a `b` column to make things easier. While we're at it, we'll compute the indirect effects, too.

```{r, warning = F}
draws <- draws %>%
  mutate(b = b_liking_respappr) %>% 
  mutate(a1b = a1 * b,
         a2b = a2 * b)

draws %>% 
  pivot_longer(a1b:a2b) %>% 
  group_by(name) %>% 
  mean_qi(value) %>% 
  mutate_if(is.double, round, digits = 3)
```

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

```{r, warning = F}
draws <- draws %>% 
  mutate(c1 = c1_prime + a1b,
         c2 = c2_prime + a2b)

draws %>% 
  pivot_longer(c1:c2) %>% 
  group_by(name) %>% 
  summarize(mean = mean(value),
            sd   = sd(value))
```

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

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

protest %>% 
  select(protest, group)
```

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

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

Check the summary.

```{r}
print(model6.5)
```

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

```{r, fig.width = 6, fig.height = 1, message = F}
# 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),
                  linewidth = 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.

```{r}
m_model <- bf(respappr ~ 0 + group)
y_model <- bf(liking   ~ 0 + group + respappr)
```

Now fit `model6` using the `mvbf()` approach.

```{r model6.6}
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?

```{r}
print(model6.6)
```

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.

```{r, warning = F}
draws <- as_draws_df(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)

draws %>% 
  pivot_longer(a1b:a3b) %>% 
  group_by(name) %>% 
  mean_qi(value) %>% 
  mutate_if(is.double, round, digits = 3)
```

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.

```{r, warning = F}
# compute the difference scores
draws <- draws %>% 
  mutate(diff_individual_minus_none       = a2b - a1b,
         diff_collective_minus_none       = a3b - a1b,
         diff_collective_minus_individual = a3b - a2b)

# summarize
draws %>% 
  pivot_longer(contains("diff")) %>% 
  group_by(name) %>% 
  mean_qi(value) %>% 
  mutate_if(is.double, round, digits = 3)
```

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 {-}

```{r}
sessionInfo()
```

```{r, echo = F, message = F, warning = F, results = "hide"}
pacman::p_unload(pacman::p_loaded(), character.only = TRUE)
```