08.Rmd

# Extending the Fundamental Principles of Moderation Analysis

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

"In this chapter, [we'll see] how [the] principles of moderation analysis are applied when the moderator is dichotomous (rather than a continuum, as in the previous chapter) as well as when both focal antecedent and moderator are continuous" [@hayesIntroductionMediationModeration2018, p. 267].

## Moderation with a dichotomous moderator

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

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

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

glimpse(disaster)
```

Regardless of whether the antecedent variables are continuous or binary, the equation for the simple moderation is still

$$Y = i_Y + b_1 X + b_2 W + b_3 XW + e_Y.$$

We can use that equation to fit our first moderation model with a binary $W$ (i.e., `frame`) like so.

```{r model8.1}
model8.1 <- brm(
  data = disaster, 
  family = gaussian,
  justify ~ 1 + skeptic + frame + frame:skeptic,
  cores = 4,
  file = "fits/model08.01")
```

Check the summary.

```{r}
print(model8.1, digits = 3)
```

We'll compute our Bayesian $R^2$ in the typical way.

```{r}
bayes_R2(model8.1) %>% round(digits = 3)
```

> This model should look familiar to you because it is exactly the same model estimated in the analysis presented in Chapter 7 (see Table 7.4, model 3). The only differences between these two analyses are how the corresponding question is framed, meaning which variable is deemed the focal antecedent and which is the moderator, and how these variables are symbolically labeled as $X$ and $W$. In the analysis in Chapter 7, the focal antecedent variable was a dichotomous variable coding the framing of the cause of the disaster (labeled $X$ then, but $W$ now), whereas in this analysis, the focal antecedent is a continuous variable placing each person on a continuum of climate change skepticism (labeled $W$ then, but $X$ now), with the moderator being a dichotomous variable coding experimental condition.
>
> So this example illustrates the symmetry property of interactions introduced in [section 7.1][Conditional and unconditional effects]. (p. 272)

### Visualizing and probing the interaction.

For the plots in this chapter, we'll take our color palette from the [**ochRe** package](https://github.com/ropenscilabs/ochRe), which provides Australia-inspired colors. We'll also use a few theme settings from good-old [**ggthemes**](https://CRAN.R-project.org/package=ggthemes). As in the last chapter, we'll save our adjusted theme settings as an object, `theme_08`.

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

theme_08 <- theme_minimal() +
  theme(panel.grid.minor = element_blank(),
        plot.background = element_rect(fill = ochre_palettes[["olsen_seq"]][8],
                                       color = "transparent"))
```

Happily, the **ochRe** package has a handy convenience function, `viz_palette()`, that makes it easy to preview the colors available in a given palette. We'll be using "olsen_qual" and "olsen_seq".

```{r, fig.width = 8, fig.height = 3}
viz_palette(ochre_palettes[["olsen_qual"]])
viz_palette(ochre_palettes[["olsen_seq"]])
```

Behold our Figure 8.3.

```{r, fig.width = 5.5, fig.height = 4.75}
# these will come in handy with `geom_text()`, below
green_slope <- (fixef(model8.1)["skeptic", 1] + fixef(model8.1)[4, 1]) %>% round(digits = 3)
blue_slope  <-  fixef(model8.1)["skeptic", 1] %>% round(digits = 3)

# define the new data
nd <- crossing(frame   = 0:1,
               skeptic = c(0, 7))

# extract the `fitted()` draws
fitted(model8.1, 
       newdata = nd, 
       summary = F) %>% 
  # wrangle
  data.frame() %>% 
  set_names(mutate(nd, name = str_c(frame, "_", skeptic)) %>% pull(name)) %>% 
  mutate(iter = 1:n()) %>% 
  pivot_longer(-iter) %>% 
  separate(name, into = c("frame", "skeptic"), sep = "_") %>% 
  mutate(skeptic = as.double(skeptic)) %>% 
  
  # plot!
  ggplot(aes(x = skeptic, y = value, 
             group = interaction(frame, iter))) +
  geom_line(aes(color = frame),
            linewidth = 1/6, alpha = 1/25) +
  annotate(geom = "text",
           x = c(4, 4.6), y = c(3.5, 2.6),
           label = c(paste("the slope when frame = 1 is about", green_slope),
                     paste("the slope when frame = 0 is about", blue_slope)),
           angle = c(28, 5), color = "black") +
  scale_color_manual(NULL, values = ochre_palettes[["olsen_qual"]][(5:6)]) +
  scale_x_continuous(breaks = 1:6) +
  coord_cartesian(xlim = c(1, 6),
                  ylim = c(2, 4)) +
  labs(title = "The moderation effect in a spaghetti plot",
       subtitle = "Using brms defaults, we have 4,000 posterior draws. Much like we have\n4,000 draws for each model parameter, when we combine those draws\nacross parameters, we end up with 4,000 model equations. The\nconsequence is we can also express the line plot as 4,000\nsemitransparent lines, one for each level of our moderator, frame.",
       x = expression("Climate Change Skepticism "*(italic(X))),
       y = "Strength of Justification for Withholding Aid") +
  theme_08 +
  theme(legend.position = "none")
```

In addition to our fancy Australia-inspired colors, we'll also play around a bit with spaghetti plots in this chapter. To my knowledge, this use of spaghetti plots is uniquely Bayesian. If you're trying to wrap your head around what on earth we just did, take a look at the structure of our `as_draws_df()` object, `draws`.

```{r}
draws <- as_draws_df(model8.1)

head(draws)
```

The `head()` function returned six rows, each one corresponding to one of the posterior draws. There are 7 primary columns and 3 additional metadata columns indexing the HMC draws. The final 2 primary columns, `lprior` and `lp__`, are beyond the scope of this ebook. You might think of the fifth primary column, `sigma`, as the Bayesian analogue to what the OLS folks often refer to as *error* or the *residual variance*. Hayes doesn't tend to emphasize it in this text, but it's something you'll want to pay increasing attention to as you move along in your Bayesian career. All the columns starting with `b_` are the regression parameters, the *model coefficients* or the *fixed effects*. But anyways, notice that those `b_` columns correspond to the four parameter values in Formula 8.2 on page 270. Here they are, but reformatted to more closely mimic the text:

1. $\hat{Y}$ = `r round(draws[1, 1], 3)` + `r round(draws[1, 2], 3)`$X$ + `r round(draws[1, 3], 3)`$W$ + `r round(draws[1, 4], 3)`*XW*
2. $\hat{Y}$ = `r round(draws[2, 1], 3)` + `r round(draws[2, 2], 3)`$X$ + `r round(draws[2, 3], 3)`$W$ + `r round(draws[2, 4], 3)`*XW*
3. $\hat{Y}$ = `r round(draws[3, 1], 3)` + `r round(draws[3, 2], 3)`$X$ + `r round(draws[3, 3], 3)`$W$ + `r round(draws[3, 4], 3)`*XW*
4. $\hat{Y}$ = `r round(draws[4, 1], 3)` + `r round(draws[4, 2], 3)`$X$ + `r round(draws[4, 3], 3)`$W$ + `r round(draws[4, 4], 3)`*XW*
5. $\hat{Y}$ = `r round(draws[5, 1], 3)` + `r round(draws[5, 2], 3)`$X$ + `r round(draws[5, 3], 3)`$W$ + `r round(draws[5, 4], 3)`*XW*
6. $\hat{Y}$ = `r round(draws[6, 1], 3)` + `r round(draws[6, 2], 3)`$X$ + `r round(draws[6, 3], 3)`$W$ + `r round(draws[6, 4], 3)`*XW*

Each row of `draws`, each iteration or posterior draw, yields a full model equation that is a credible description of the data--or at least as credible as we can get within the limits of the model we have specified, our priors (which we typically cop out on and just use defaults in this ebook), and how well those fit when applied to the data at hand. So when we use **brms** convenience functions like `fitted()`, we pass specific predictor values through those 4,000 unique model equations, which returns 4,000 similar but distinct expected $Y$-values. So although a nice way to summarize those 4,000 values is with summaries such as the posterior mean/median and 95% intervals, another way is to just plot an individual regression line for each of the iterations. That is what's going on when we depict out models with a spaghetti plot.

The thing I like about spaghetti plots is that they give a three-dimensional sense of the posterior. Note that each individual line is very skinny and semitransparent. When you pile a whole bunch of them atop each other, the peaked or most credible regions of the posterior are the most saturated in color. Less credible posterior regions almost seamlessly merge into the background. Also, note how the combination of many similar but distinct straight lines results in a bowtie shape. Hopefully this clarifies where that shape's been coming from when we use `geom_ribbon()` to plot the 95% intervals.

Back to the text, on the bottom of page 274, Hayes pointed out the conditional effect of `skeptic` when `frame == 1` is $b_1 + b_3 = 0.306$. We can show that with a little arithmetic followed up with `tidybayes::mean_qi()`.

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

draws <- draws %>% 
  mutate(`b1 + b3` = b_skeptic + `b_skeptic:frame`)

draws %>% 
  mean_qi(`b1 + b3`) %>% 
  mutate_if(is.double, round, digits = 3)
```

But anyways, you could recode `frame` in a number of ways, including `if_else()` or, in this case, by simple arithmetic.

```{r}
disaster <- disaster %>% 
  mutate(frame_ep = 1 - frame)
```

With `frame_ep` in hand, we're ready to refit the model.

```{r model8.2, warning = F, message = F}
model8.2 <- update(
  model8.1, 
  newdata = disaster,
  formula = justify ~ 1 + skeptic + frame_ep + frame_ep:skeptic,
  cores = 4,
  file = "fits/model08.02")
```

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

Our results match nicely with the formula on page 275.

If you want to follow along with Hayes on page 276 and isolate the 95% credible intervals for the `skeptic` parameter, you can use the `posterior_interval()` function.

```{r}
posterior_interval(model8.2)["b_skeptic", ] %>% round(digits = 3)
```

## Interaction between two quantitative variables

Here's the `glbwarm` data.

```{r, message = F}
glbwarm <- read_csv("data/glbwarm/glbwarm.csv")

glimpse(glbwarm)
```

In this section we add three covariates (i.e., $C$ variables) to the basic moderation model. Although Hayes made a distinction between the $X$, $M$, and $C$ variables in the text, that distinction is conceptual and doesn't impact the way we enter them into `brm()`. Rather, the `brm()` formula clarifies they're all just predictors.

```{r model8.3}
model8.3 <- brm(
  data = glbwarm, 
  family = gaussian,
  govact ~ 1 + negemot + age + negemot:age + posemot + ideology + sex,
  cores = 4,
  file = "fits/model08.03")
```

Our results cohere nicely with the Hayes's formula in the middle of page 278 or with the results he displayed in Table 8.2.

```{r}
print(model8.3, digits = 3)
```

Here's the $R^2$ summary.

```{r}
bayes_R2(model8.3) %>% round(digits = 3)
```

As the $R^2$ is a good bit away from the boundaries, it's nicely Gaussian.

```{r, fig.width = 8, fig.height = 1.75}
bayes_R2(model8.3, summary = F) %>% 
  data.frame() %>% 

  ggplot(aes(x = R2)) +
  geom_density(linewidth = 0, fill = ochre_palettes[["olsen_qual"]][(1)]) +
  scale_x_continuous(NULL, limits = 0:1) +
  scale_y_continuous(NULL, breaks = NULL) +
  labs(title = expression("Loot at how Gaussian our "*italic(R)^2*" is!"),
       subtitle = expression("No need to put a "*italic(p)*"-value on that sucker.")) +
  theme_08
```

### Visualizing and probing the interaction.

For our version of Figure 8.5, we'll need to adjust our `nd` data for `fitted()`.

```{r}
(
  nd <- crossing(negemot = c(1, 6),
                 age     = c(30, 50, 70)) %>% 
  mutate(posemot  = mean(glbwarm$posemot),
         ideology = mean(glbwarm$ideology),
         sex      = mean(glbwarm$sex))
)
```

Our `fitted()` and **ggplot2** code will be quite similar to the last spaghetti plot. Only this time we'll use `filter()` to reduce the number of posterior draws we show in the plot.

```{r, fig.width = 8, fig.height = 4.5, warning = F}
# these will come in handy with `geom_text()`, below
slope_30 <- (fixef(model8.3)["negemot", 1] + fixef(model8.3)["negemot:age", 1] * 30) %>% round(digits = 3)
slope_50 <- (fixef(model8.3)["negemot", 1] + fixef(model8.3)["negemot:age", 1] * 50) %>% round(digits = 3)
slope_70 <- (fixef(model8.3)["negemot", 1] + fixef(model8.3)["negemot:age", 1] * 70) %>% round(digits = 3)

# annotation
text <- tibble(
  negemot = 6.1,
  value   = c(5.4, 5.7, 6),
  age     = as.character(c(30, 50, 70)),
  iter    = 0,
  label   = str_c("the slope for ", c(30, 50, 70)," year olds is about ", c(slope_30, slope_50, slope_70)))

# get the fitted draws and wrangle
fitted(model8.3, 
       newdata = nd, 
       summary = F) %>% 
  data.frame() %>% 
  set_names(mutate(nd, name = str_c(negemot, "_", age)) %>% pull(name)) %>% 
  mutate(iter = 1:n()) %>% 
  pivot_longer(-iter) %>% 
  separate(name, into = c("negemot", "age"), sep = "_") %>% 
  mutate(negemot = as.double(negemot)) %>% 
  filter(iter <= 100) %>% 
  
  # plot!
  ggplot(aes(x = negemot, y = value, 
             group = interaction(age, iter),
             color = age)) +
  geom_line(linewidth = 3/4, alpha = 1/8) +
  geom_text(data = text,
            aes(label = label),
            hjust = 0) +
  scale_color_manual(NULL, values = ochre_palettes[["olsen_qual"]][c(4, 1, 3)]) +
  scale_x_continuous(breaks = 1:6, limits = c(1, 9)) +
  scale_y_continuous("Support for Government Action", limits = c(2.9, 6.1)) +
  labs(title = "The moderation effect of age on negemot",
       x = expression("Negative Emotions about Climate Change "*(italic(X)))) +
  theme_08 +
  theme(legend.position = "none")
```

When we reduce the number of lines depicted in the plot, we lose some of the three-dimensional illusion. It's nice, however, to get a closer look to each individual line. To each their own.

We'll continue with our spaghetti plot approach for Figure 8.7. Again, when we made the JN technique plot for Chapter 7, we computed values for the posterior mean and the 95% intervals. Because the intervals follow a bowtie shape, we had to compute the $Y$-values for many values across the $x$-axis in order to make the curve look smooth. But as long as we stick with the spaghetti plot approach, all we need are the values at the endpoints of each iteration. Although each line is straight, the combination of many lines is what produces the bowtie effect.

```{r, fig.width = 6, fig.height = 4.5, warning = F}
# here is our primary data object
draws <- as_draws_df(model8.3) %>% 
  transmute(at_15 = b_negemot + `b_negemot:age` * 15,
            at_90 = b_negemot + `b_negemot:age` * 90,
            .draw = .draw) %>% 
  pivot_longer(-.draw, names_to = "age") %>% 
  mutate(age = str_remove(age, "at_") %>% as.double())

# here we compute the points for the posterior mean
post_means <- draws %>% 
  group_by(age) %>% 
  summarize(value = mean(value)) %>% 
  mutate(.draw = 0)

# plot!
draws %>% 
  filter(.draw < 501) %>% 

  ggplot(aes(x = age, y = value, group = .draw)) +
  geom_line(color = ochre_palettes[["olsen_qual"]][1],
            alpha = 1/12) +
  geom_line(data = post_means,
            color = ochre_palettes[["olsen_qual"]][3],
            linewidth = 1.1) +
  scale_y_continuous(breaks = seq(from = -.25, to = 1, by = .25)) +
  coord_cartesian(xlim = c(20, 85),
                  ylim = c(-.25, 1)) +
  labs(subtitle = "Each orange line is the consequence of one of 1,000 posterior draws. All are\ncredible regression lines for the data, but because they are semitransparent,\nthe most credible regions are the ones with the darkest color. The line\ncorresponding to the posterior mean is in red-orange.",
       x = expression(Age~(italic(W))),
       y = "Conditional Effects of Negative Emotions on\nSupport for Government Action") +
  theme_08
```

In other words, each of those orange lines is a credible expression of $\theta_{X \rightarrow Y}$ (i.e., $b_1 + b_3 W$) across a continuous range of $W$ values.

## Hierarchical versus simultaneous entry

> Many investigators test a moderation hypothesis in regression analysis using a method that on the surface seems different than the procedure described thus far. This alternative approach is to build a regression model by adding the product of $X$ and $W$ to a model already containing $X$ and $W$. This procedure is sometimes called *hierarchical regression* or *hierarchical variable entry* (and easily confused by name with *hierarchical linear modeling*, which is an entirely different thing). The goal using this method is to determine whether allowing $X$'s effect to be contingent on $W$ produces a better fitting model than one in which the effect of $X$ is constrained to be unconditional on $W$. According to the logic of hierarchical entry, if the contingent model accounts for more of the variation in Y than the model that forces $X$'s effect to be independent of $W$, then the better model is one in which $W$ is allowed to moderate $X$'s effect. Although this approach works, it is a widely believed myth that it is *necessary* to use this approach in order to test a moderation hypothesis. (p. 289, *emphasis* in the original)

Although this method is not necessary, it can be handy to slowly build your model. This method can also serve nice rhetorical purposes in a paper. Anyway, here's our multivariable but non-moderation model, `model8.4`.

```{r model8.4}
model8.4 <- update(
  model8.1,
  formula = justify ~ 1 + skeptic + frame,
  cores = 4,
  file = "fits/model08.04")
```

Here we'll compute the corresponding $R^2$ and compare it with the one for the original interaction model with a difference score.

```{r, warning = F}
r2 <- tibble(moderation_model    = bayes_R2(model8.1, summary = F)[, 1],
             multivariable_model = bayes_R2(model8.4, summary = F)[, 1]) %>% 
  # we'll need a difference score
  mutate(difference = moderation_model - multivariable_model) %>% 
  # putting the data in the long format and grouping will make summarizing easier
  pivot_longer(everything(), names_to = "R2")

r2 %>% 
  group_by(R2) %>% 
  median_qi(value) %>% 
  mutate_if(is.double, round, digits = 3)
```

Note that the Bayesian $R^2$ performed differently than the $F$-test in the text.

```{r, fig.width = 8, fig.height = 2.5}
r2 %>% 
  filter(R2 == "difference") %>% 

  ggplot(aes(x = value, fill = model)) +
  geom_density(linewidth = 0, fill = ochre_palettes[["olsen_seq"]][14]) +
  scale_x_continuous(NULL, limits = c(-.4, .4)) +
  scale_y_continuous(NULL, breaks = NULL) +
  labs(title = expression(paste("The Bayesian ", Delta, italic(R)^2, " distribution")),
       subtitle = "Although most of the posterior mass is positive--suggesting the moderation model accounted for more variance than\nthe simple multivariable model--, a substantial portion of the postrior is within the negative parameter space. Sure,\nif we had to bet, the safer bet is on the moderation model. But that bet wouled be quite uncertain and we might well\nloose our shirts. Also, note the width of the distribution; credible values range from -0.1 to nearly 0.2.") +
  theme_08
```

We can also compare these with the LOO, which, as is typical of information criteria, corrects for model complexity. First, we compute them and attach the results to the model fit objects.

```{r}
model8.1 <- add_criterion(model8.1, "loo")
model8.4 <- add_criterion(model8.4, "loo")
```

Now use the `loo_compare()` function to compare them directly.

```{r}
loo_compare(model8.1, model8.4) %>% 
  print(simplify = F)
```

As a reminder, we generally prefer models with lower information criteria, which in this case is clearly the moderation model (i.e., `model8.1`). However, the standard error value (i.e., `se_diff`) for the difference (i.e., `elpd_diff`) is quite large, which suggests that the model with the lowest value isn't the clear winner. Happily, these results match nicely with the Bayesian $R^2$ difference score. The moderation model appears somewhat better than the multivariable model, but its superiority is hardly decisive.

## The equivalence between moderated regression analysis and a 2 X 2 factorial analysis of variance

I'm just not going to encourage ANOVA $F$-testing methodology. However, I will show the Bayesian regression model. First, here are the data.

```{r, message = F}
caskets <- read_csv("data/caskets/caskets.csv")

glimpse(caskets)
```

Fit the moderation model.

```{r model8.5}
model8.5 <- brm(
  data = caskets, 
  family = gaussian,
  interest ~ 1 + policy + kerry + policy:kerry,
  cores = 4,
  file = "fits/model08.05")
```

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

Those results don't look anything like what Hayes reported in Tables 8.3 or 8.4. However, a little deft manipulation of the posterior samples can yield equivalent results to Hayes's Table 8.3.

```{r}
draws <- as_draws_df(model8.5) %>% 
  mutate(Y_bar_1 = b_Intercept + b_policy * 0 + b_kerry * 0 + `b_policy:kerry` * 0 * 0,
         Y_bar_3 = b_Intercept + b_policy * 0 + b_kerry * 1 + `b_policy:kerry` * 0 * 1,
         Y_bar_2 = b_Intercept + b_policy * 1 + b_kerry * 0 + `b_policy:kerry` * 1 * 0,
         Y_bar_4 = b_Intercept + b_policy * 1 + b_kerry * 1 + `b_policy:kerry` * 1 * 1,
         
         Y_bar_12 = b_Intercept + b_policy * .5 + b_kerry * 0  + `b_policy:kerry` * .5 * 0,
         Y_bar_34 = b_Intercept + b_policy * .5 + b_kerry * 1  + `b_policy:kerry` * .5 * 1,
         Y_bar_13 = b_Intercept + b_policy * 0  + b_kerry * .5 + `b_policy:kerry` * 0  * .5,
         Y_bar_24 = b_Intercept + b_policy * 1  + b_kerry * .5 + `b_policy:kerry` * 1  * .5)
```

Here are the cell-specific means in Table 8.3.

```{r, warning = F}
draws %>% 
  pivot_longer(Y_bar_1:Y_bar_4) %>% 
  group_by(name) %>% 
  median_qi(value) %>% 
  mutate_if(is.double, round, digits = 3)
```

And here are the marginal means from Table 8.3.

```{r, warning = F, message = F}
draws %>% 
  pivot_longer(Y_bar_12:Y_bar_24) %>% 
  group_by(name) %>% 
  median_qi(value) %>% 
  mutate_if(is.double, round, digits = 3)
```

For kicks and giggles, here are what the cell-specific means look like in box plots.

```{r, fig.width = 3.25, fig.height = 3, warning = F}
draws %>% 
  pivot_longer(Y_bar_1:Y_bar_4) %>% 

  ggplot(aes(x = name, y = value, fill = name)) +
  geom_boxplot(linewidth = 1/3) +
  scale_fill_manual(values = ochre_palettes[["olsen_qual"]][c(5, 6, 4, 3)]) +
  labs(title = "Cell-specific effects",
       x = NULL, 
       y = "interest") +
  theme_08 +
  theme(legend.position = "none")
```
         
And here are the same for the marginal means. This time we'll show the shapes of the posteriors with violin plots with horizontal lines depicting the median and interquartile ranges.

```{r, fig.width = 3.25, fig.height = 3, warning = F}
draws %>% 
  pivot_longer(Y_bar_12:Y_bar_24) %>% 

  ggplot(aes(x = name, y = value, fill = name)) +
  geom_violin(draw_quantiles = c(.25, .5, .75),
              color = ochre_palettes[["olsen_seq"]][8]) +
  scale_fill_manual(values = ochre_palettes[["olsen_qual"]][c(5, 6, 4, 3)]) +
  labs(title = "Marginal means",
       x = NULL, 
       y = "interest") +
  theme_08 +
  theme(legend.position = "none")
```

On page 294, Hayes used point estimates to compute the *simple effect* of policy information among Kerry supporters and then the same thing among Bush supporters. Here's how we'd do that when working with the full vector of posterior draws.

```{r, warning = F}
draws %>% 
  mutate(simple_effect_Kerry = Y_bar_4 - Y_bar_3,
         simple_effect_Bush  = Y_bar_2 - Y_bar_1) %>% 
  pivot_longer(contains("simple")) %>% 
  group_by(name) %>% 
  mean_qi(value) %>% 
  mutate_if(is.double, round, digits = 3)
```

So then computing the main effect for policy information using the simple effects is little more than an extension of those steps.

```{r, warning = F}
draws %>% 
  transmute(main_effect = ((Y_bar_4 - Y_bar_3) + (Y_bar_2 - Y_bar_1)) / 2) %>% 
  mean_qi(main_effect) %>% 
  mutate_if(is.double, round, digits = 3)
```

And we get the same results by strategically subtracting the marginal means.

```{r, warning = F}
draws %>% 
  transmute(main_effect = Y_bar_24 - Y_bar_13) %>% 
  mean_qi(main_effect) %>% 
  mutate_if(is.double, round, digits = 3)
```

The main effect of for candidate is similarly computed using either approach.

```{r, warning = F}
draws %>% 
  transmute(main_effect = ((Y_bar_4 - Y_bar_2) + (Y_bar_3 - Y_bar_1)) / 2) %>% 
  mean_qi(main_effect) %>% 
  mutate_if(is.double, round, digits = 3)

draws %>% 
  transmute(main_effect = Y_bar_34 - Y_bar_12) %>% 
  mean_qi(main_effect) %>% 
  mutate_if(is.double, round, digits = 3)
```

We don't have an $F$-test for our Bayesian moderation model. But we do have an interaction term. Here's its distribution.

```{r, fig.width = 3.25, fig.height = 3}
draws %>% 
  ggplot(aes(x = `b_policy:kerry`)) +
  geom_density(linewidth = 0, fill = ochre_palettes[["olsen_qual"]][2]) +
  geom_vline(xintercept = fixef(model8.5)["policy:kerry", c(1, 3, 4)],
             color = ochre_palettes[["olsen_seq"]][8], linetype = c(1, 2, 2)) +
  scale_x_continuous(NULL,
                     breaks = fixef(model8.5)["policy:kerry", c(1, 3, 4)],
                     labels = fixef(model8.5)["policy:kerry", c(1, 3, 4)] %>% round(digits = 2)) +
  scale_y_continuous(NULL, breaks = NULL) +
  labs(title = "The interaction term, `policy:kerry`",
       subtitle = "The solid vertical line is the posterior mean\nand the dashed lines to either end denote the\npercentile-based 95% intervals.") +
  theme_08 +
  theme(legend.position = "none")
```

Following Hayes's work on the bottom of page 295, here's how you'd reproduce that by manipulating our $\overline Y$ vectors.

```{r, warning = F}
draws %>% 
  transmute(reproduced_interaction_term = (Y_bar_4 - Y_bar_3) - (Y_bar_2 - Y_bar_1)) %>% 
  mean_qi(reproduced_interaction_term) %>% 
  mutate_if(is.double, round, digits = 2)
```

Extending that logic, we also get the answer this way.

```{r, warning = F}
draws %>% 
  transmute(reproduced_interaction_term = (Y_bar_4 - Y_bar_2) - (Y_bar_3 - Y_bar_1)) %>% 
  mean_qi(reproduced_interaction_term) %>% 
  mutate_if(is.double, round, digits = 2)
```

### Simple effects parameterization.

We might reacquaint ourselves with the `formula` from `model8.5`.

```{r}
model8.5$formula
```

The results cohere nicely with the "Model 1" results at the top of Table 8.5.

```{r}
fixef(model8.5) %>% round(digits = 3)
```

The Bayesian $R^2$ portion looks on point, too.

```{r}
bayes_R2(model8.5) %>% round(digits = 3)
```

Our various `Y_bar` transformations from before continue to cohere with the coefficients, above, just like in the text. E.g., the `policy` coefficient may be returned like so.

```{r, warning = F}
draws %>% 
  transmute(b1                  = b_policy,
            `Y_bar_2 - Y_bar_1` = Y_bar_2 - Y_bar_1) %>% 
  pivot_longer(everything()) %>% 
  group_by(name) %>% 
  summarize(mean = mean(value),
            sd   = sd(value)) %>% 
  mutate_if(is.double, round, digits = 3)
```

We can continue to use Hayes's `Y_bar` transformations to return the `kerry` coefficient, too.

```{r, warning = F}
draws %>% 
  transmute(b2                  = b_kerry,
            `Y_bar_3 - Y_bar_1` = Y_bar_3 - Y_bar_1) %>% 
  pivot_longer(everything()) %>% 
  group_by(name) %>% 
  summarize(mean = mean(value),
            sd   = sd(value)) %>% 
  mutate_if(is.double, round, digits = 3)
```

Here we compute $b_3$ with the difference between the simple effects of $X$ at levels of $W$.

```{r, warning = F}
draws %>% 
  transmute(b3                                          = `b_policy:kerry`,
            `(Y_bar_4 - Y_bar_3) - (Y_bar_2 - Y_bar_1)` = (Y_bar_4 - Y_bar_3) - (Y_bar_2 - Y_bar_1)) %>% 
  pivot_longer(everything()) %>% 
  group_by(name) %>% 
  summarize(mean = mean(value),
            sd   = sd(value)) %>% 
  mutate_if(is.double, round, digits = 3)
```

And now $b_{3}$ with the difference between the simple effects of $W$ at levels of $X$.

```{r, warning = F}
draws %>% 
  transmute(b3                                          = `b_policy:kerry`,
            `(Y_bar_4 - Y_bar_2) - (Y_bar_3 - Y_bar_1)` = (Y_bar_4 - Y_bar_2) - (Y_bar_3 - Y_bar_1)) %>% 
  pivot_longer(everything()) %>% 
  group_by(name) %>% 
  summarize(mean = mean(value),
            sd   = sd(value)) %>% 
  mutate_if(is.double, round, digits = 3)
```

### Main effects parameterization.

A nice feature of **brms** is you can transform your data right within the `brm()` or `update()` functions. Here we'll make our two new main-effects-coded variables, `policy_me` and `kerry_me`, with the `mutate()` function right within `update()`.

```{r model8.6}
model8.6 <- update(
  model8.5,
  newdata = caskets %>%
    mutate(policy_me = policy - .5,
           kerry_me  = kerry  - .5), 
  family = gaussian,
  interest ~ 1 + policy_me + kerry_me + policy_me:kerry_me,
  cores = 4,
  file = "fits/model08.06")
```

Transforming your data within the **brms** functions won't change the original data structure. However, **brms** will save the data used to fit the model within the `brm()` object. You can access that data like so.

```{r}
model8.6$data %>%
  head()
```

But we digress. Here's our analogue to the "Model 2" portion of Table 8.5.

```{r}
rbind(fixef(model8.6), bayes_R2(model8.6)) %>% 
  round(digits = 3)
```

Like with `model8.6`, above, we'll need a bit of algebra to compute our $\overline Y_i$ vectors.

```{r}
draws <- as_draws_df(model8.6) %>% 
  mutate(Y_bar_1 = b_Intercept + b_policy_me * -.5 + b_kerry_me * -.5 + `b_policy_me:kerry_me` * -.5 * -.5,
         Y_bar_3 = b_Intercept + b_policy_me * -.5 + b_kerry_me *  .5 + `b_policy_me:kerry_me` * -.5 *  .5,
         Y_bar_2 = b_Intercept + b_policy_me *  .5 + b_kerry_me * -.5 + `b_policy_me:kerry_me` *  .5 * -.5,
         Y_bar_4 = b_Intercept + b_policy_me *  .5 + b_kerry_me *  .5 + `b_policy_me:kerry_me` *  .5 *  .5)
```

With our `draws` for `fit5` in hand, we'll follow the formulas at the top of page 298 to compute our $b_1$ and $b_2$ distributions.

```{r, warning = F}
draws %>% 
  transmute(b1 = ((Y_bar_4 - Y_bar_3) + (Y_bar_2 - Y_bar_1)) / 2,
            b2 = ((Y_bar_4 - Y_bar_2) + (Y_bar_3 - Y_bar_1)) / 2) %>% 
  pivot_longer(b1:b2) %>% 
  group_by(name) %>% 
  summarize(mean = mean(value),
            sd   = sd(value)) %>% 
  mutate_if(is.double, round, digits = 3)
```

Hayes pointed out that the interaction effect, $b_3$, is the same across models his OLS Models 1 and 2. This is largely true for our Bayesian HMC `model8.5` and `model8.6` models.

```{r}
fixef(model8.5)[4, ] %>% round(digits = 3)
fixef(model8.6)[4, ] %>% round(digits = 3)
```

However, the results aren't exactly the same because of simulation error. If you were working on a project requiring high precision, increase the number of posterior iterations. To demonstrate, here we'll increase each chain's post-warmup iteration count by an order of magnitude, resulting in 80,000 post-warmup iterations rather than the default 4,000.

```{r model8.7, warning = F, message = F}
model8.7 <- update(
  model8.5,
  chains = 4, cores = 4, warmup = 1000, iter = 21000,
  file = "fits/model08.07")

model8.8 <- update(
  model8.6,
  chains = 4, cores = 4, warmup = 1000, iter = 21000,
  file = "fits/model08.08")
```

Now they're quite a bit closer. 

```{r}
fixef(model8.7)[4, ] %>% round(digits = 3)
fixef(model8.8)[4, ] %>% round(digits = 3)
```

And before you get fixate on how there are still differences after 80,000 iterations, each, consider comparing the two density plots.

```{r, fig.width = 4.5, fig.height = 3.25, warning = F}
bind_rows(
  as_draws_df(model8.7) %>% transmute(iteraction = `b_policy:kerry`),
  as_draws_df(model8.8) %>% transmute(iteraction = `b_policy_me:kerry_me`)
) %>% 
  mutate(model = rep(str_c("8.", 7:8), times = n() / 2)) %>% 
  
  ggplot(aes(x = iteraction, fill = model)) +
  geom_density(linewidth = 0, alpha = 1/2) +
  scale_fill_manual(values = ochre_palettes[["olsen_qual"]][c(3, 6)]) +
  scale_y_continuous(NULL, breaks = NULL) +
  labs(title = "The interaction densities, by model",
       subtitle = "Yes, they are indeed different. And yet that difference is so\ntrivial that we'd expect greater variability from measurement\nerror than we still have from MCMC simulation error.",
       x = NULL) +
  theme_08 
```

### Conducting ~~a $2 \times 2$ between-participants factorial ANOVA using PROCESS~~ another regression model with brms.

Since we're square in single-level regression land with our **brms** approach, there's no direct analogue for us, here. However, notice the post-ANOVA $t$-tests Hayes presented on page 300. If we just want to consider the $2 \times 2$ structure of our two dummy variables as indicative of four groups, we have one more coding system available for the job. With the handy `str_c()` function, we'll concatenate the `policy` and `kerry` values into a nominal variable, `policy_kerry`. Here's what that looks like:

```{r}
caskets <- caskets %>% 
  mutate(policy_kerry = str_c(policy, kerry))

head(caskets)
```

Now check out what happens if we reformat our formula to `interest ~ 0 + policy_kerry`.

```{r model8.9}
model8.9 <- brm(
  data = caskets, 
  family = gaussian,
  interest ~ 0 + policy_kerry,
  cores = 4,
  file = "fits/model08.09")
```

The `brm()` function recognized `policy_kerry` as a character vector and treated it as a nominal variable. The `0 + ` part of the function removed the model intercept. Here's how that effects the output.

```{r}
print(model8.9)
```

Without the typical intercept, `brm()` estimated the means for each of the four `policy_kerry` groups. It's kinda like an intercept-only model, but with four intercepts. Here's what their densities look like.

```{r, fig.width = 8, fig.height = 2.5, warning = F}
draws <- as_draws_df(model8.9)

draws %>% 
  pivot_longer(b_policy_kerry00:b_policy_kerry11) %>% 
  mutate(name = str_remove(name, "b_")) %>% 
  
  ggplot(aes(x = value, fill = name)) +
  geom_density(color = "transparent", alpha = 2/3) +
  scale_fill_manual(NULL,
                    values = ochre_palettes[["olsen_qual"]][c(5, 6, 4, 3)]) +
  scale_y_continuous(NULL, breaks = NULL) +
  labs(title = "Group means",
       x = NULL) +
  theme_08
```

Since each of the four primary vectors in our `draws` object is of a group mean, it's trivial to compute difference scores. To compute the difference score analogous to Hayes's two $t$-tests, we'd do the following.

```{r, warning = F}
draws %>% 
  mutate(difference_1 = b_policy_kerry10 - b_policy_kerry00,
         difference_2 = b_policy_kerry11 - b_policy_kerry01) %>% 
  pivot_longer(starts_with("difference_")) %>% 
  group_by(name) %>% 
  median_qi(value) %>% 
  mutate_if(is.double, round, digits = 3)
```

## Session info {-}

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

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