# 7 Fundamentals of Moderation Analysis

The effect of \(X\) on some variable \(Y\) is moderated by \(W\) if its size, sign, or strength depends on or can be predicted by \(W\). In that case, \(W\) is said to be a

moderatorof \(X\)’s effect on \(Y\), or that \(W\) and \(X\) interact in their influence on \(Y\). Identifying a moderator of an effect helps to establish the boundary conditions of that effect or the circumstances, stimuli, or type of people for which the effect is large versus small, present versus absent, positive versus negative, and so forth. (Andrew F. Hayes, 2018, p. 220,emphasisin the original)

## 7.1 Conditional and unconditional effects

If \(X\)’s effect on \(Y\) is moderated by another variable in the model, that means \(X\)’s effect depends on that other variable. But this model constrains \(X\)’s effect to be unconditional on \(W\), meaning that it is invariant across all values of \(W\). (p. 224)

Hayes’s Table 7.1 and the related Figure 7.2 showcase this well. You might reproduce the data for both like this.

```
library(tidyverse)
<-
d crossing(x = -1:2,
w = 0:2) %>%
mutate(A = 4 + 1 * x + 2 * w,
B = 4 + 1 * x + 2 * w + 1.5 * x * w)
d
```

```
## # A tibble: 12 x 4
## x w A B
## <int> <int> <dbl> <dbl>
## 1 -1 0 3 3
## 2 -1 1 5 3.5
## 3 -1 2 7 4
## 4 0 0 4 4
## 5 0 1 6 6
## 6 0 2 8 8
## 7 1 0 5 5
## 8 1 1 7 8.5
## 9 1 2 9 12
## 10 2 0 6 6
## 11 2 1 8 11
## 12 2 2 10 16
```

In previous chapters, we altered our plot themes using either built-in settings from **ggplot2** or extra themes from the **ggthemes** package. When we wanted to change them further, we did so with extra `theme()`

arguments. One, it’s nice to be know how to make one’s own custom theme and, two, it’d be handy to condense our code a bit. Here we’ll take a few cues from the Building a New Theme section from Peng, Kross, and Anderson’s (2017) text, *Mastering software development in R*. We’ll make our own theme, `theme_07`

, by saving a handful of augmentations from the default `theme_gray()`

theme.

```
library(dutchmasters)
<-
theme_07 theme_gray() +
theme(axis.text = element_text(color = dutchmasters$little_street[6]),
legend.background = element_rect(fill = "transparent"),
legend.key = element_rect(fill = "transparent", color = "transparent"),
panel.background = element_rect(fill = dutchmasters$little_street[2]),
panel.grid = element_blank(),
plot.background = element_rect(fill = dutchmasters$little_street[7]),
strip.background = element_rect(fill = alpha(dutchmasters$little_street[5], 2/3),
color = "transparent"),
text = element_text(family = "Courier",
color = dutchmasters$little_street[6],
size = 14))
```

This chapter’s color palette comes from the **dutchmasters** package (Thoen, 2019), which was itself based of Vermeer’s (1657) *The little street*. To get a quick glance at the full palette, we’ll also use `viz_palette()`

, a convenience function from the **ochRe** package (Allan et al., 2017).

```
# devtools::install_github("ropenscilabs/ochRe")
library(ochRe)
viz_palette(dutchmasters$little_street)
```

With our new `theme_07`

in hand, we’re ready to make our version of Figure 7.2.

```
library(directlabels)
%>%
d pivot_longer(A:B, values_to = "y") %>%
mutate(label = str_c("W = ", w),
w = factor(w),
strip = if_else(name == "A",
"A: where W is a covariate",
"B: where W is a moderator")) %>%
ggplot(aes(x = x, y = y, group = w, size = w, label = label)) +
geom_line(color = dutchmasters$little_street[9]) +
geom_dl(method = list(dl.trans(x = x + 0.2),
"last.points", cex = 1,
color = dutchmasters$little_street[9],
fontfamily = "Courier")) +
scale_size_manual(values = c(.5, .75, 1)) +
coord_cartesian(xlim = c(-1, 2.6),
ylim = c(0, 16)) +
labs(x = expression(italic(X)),
y = expression(hat(italic(Y)))) +
facet_wrap(~ strip) +
+
theme_07 # because we didn't alter the `legend.position` argument in our `theme_07`,
# we need to use `theme()` to remove it
theme(legend.position = "none")
```

We borrowed `geom_dl()`

form the **directlabels** package (Hocking, 2021), which made it easy to insert the “W = \(i\)” labels to the right of the lines.

I played around with the annotation in Figure 7.4 for a while and it was just way more trouble than it was worth. If you’re ever inspired to work it out, please share your code. I’m moving on.

### 7.1.1 Eliminating the constraint of unconditionality.

We can write the generic moderation model as

\[Y = i_Y + f(W) X + b_2 W + e_Y,\]

where \(f(W)\) is a function of \(W\). Consider a simple function of the form \(f(W) = b_1 + b_3 W\). This function of \(W\) looks like a simple linear regression model where \(b_1\) is the constant and \(b_3\) is the regression coefficient for \(W\), except that rather than estimating some consequent variable from \(W\), it is a model of the effect of \(X\) on \(Y\). (p. 226)

If we use that definition of \(f(W)\), we can update our equation to

\[\begin{align*} Y & = i_Y + (b_1 + b_3 W) X + b_2 W + e_Y, \text{or} \\ Y & = i_Y + b_1 X + b_2 W + b_3 XW + e_Y. \end{align*}\]

Hayes called this the *simple linear moderation model*. Out of this equation we can extract \(\theta_{X \rightarrow Y}\), the *conditional effect of* \(X\) *on* \(Y\), which takes the form

\[\theta_{X \rightarrow Y} = b_1 + b_3 W.\]

This is what Figure 7.2 was all about.

### 7.1.2 Symmetry in moderation.

It turns out that, mathematically, there’s no difference in speaking about \(X\) moderating \(W\) than speaking about \(W\) moderating \(X\). These are all just capital letters standing in for variables and perhaps you’ll recall from grade school that \(a + b + ab\) is just the same as \(b + a + ba\). Thus we can rearrange the simple moderation equations above to

\[Y = i_Y + b_1 X + (b_2 + b_3 X) W + e_Y\]

and speak instead about the *conditional effect of* \(W\) *on* \(Y\),

\[\theta_{w \rightarrow Y} = b_2 + b_3 X.\]

### 7.1.3 Interpretation of the regression coefficients.

With all this wacky stuff going on, we should clarify some of our terms:

- \(i_Y\) is still the expected value for \(Y\) with both \(X = 0\) and \(W = 0\),
- \(b_1\) is the conditional effect of \(X\) on \(Y\) when \(W = 0\), and
- \(b_2\) is the conditional effect of \(W\) on \(Y\) when \(X = 0\).

It turns out \(b_3\) requires more words, which we’ll provide in the next subsection.

### 7.1.4 The importance of \(b_3\) wen asking about mediation.

The simple moderation model allows \(X\)’s effect on \(Y\) to be a linear function of \(W\). Of course, allowing that effect to depend on \(W\) doesn’t mean that it actually does in reality. In most any sample of data, [the point estimate (i.e., posterior mean or median) for] \(b_3\) will be different from zero even when \(X\)’s effect on \(Y\) is independent of \(W\). Of interest when [modeling] a moderation hypothesis is not just allowing \(X\)’s effect to be contingent on \(W\), but also determining whether \(b_3\) deviates too far from zero than would be expected given that \(b_3\), like any statistic, is subject to sampling variance. (p. 231)

## 7.2 An example: Climate change disasters and humanitarianism

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

.

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

```
## Rows: 211
## Columns: 5
## $ id <dbl> 1, 2, 3, 4, 5, 6, 7, 8, 9, 10, 11, 12, 13, 14, 15, 16, 17, 18, 19, 20, 21, 22, 23, 24, 25, 2…
## $ frame <dbl> 1, 1, 1, 1, 1, 0, 0, 1, 0, 0, 1, 1, 0, 0, 1, 1, 1, 1, 0, 0, 1, 0, 1, 0, 1, 1, 0, 0, 0, 1, 0,…
## $ donate <dbl> 5.6, 4.2, 4.2, 4.6, 3.0, 5.0, 4.8, 6.0, 4.2, 4.4, 5.8, 6.2, 6.0, 4.2, 4.4, 5.8, 5.4, 3.4, 7.…
## $ justify <dbl> 2.95, 2.85, 3.00, 3.30, 5.00, 3.20, 2.90, 1.40, 3.25, 3.55, 1.55, 1.60, 1.65, 2.65, 3.15, 2.…
## $ skeptic <dbl> 1.8, 5.2, 3.2, 1.0, 7.6, 4.2, 4.2, 1.2, 1.8, 8.8, 1.0, 5.4, 2.2, 3.6, 7.8, 1.6, 1.0, 6.4, 3.…
```

Here is how to get the ungrouped mean and \(SD\) values for `justify`

and `skeptic`

, as presented in Table 7.3.

```
%>%
disaster pivot_longer(justify:skeptic) %>%
group_by(name) %>%
summarise(mean = mean(value),
sd = sd(value)) %>%
mutate_if(is.double, round, digits = 3)
```

```
## # A tibble: 2 x 3
## name mean sd
## <chr> <dbl> <dbl>
## 1 justify 2.87 0.93
## 2 skeptic 3.38 2.03
```

And here we get the same summary values, this time grouped by `frame`

.

```
%>%
disaster pivot_longer(justify:skeptic) %>%
group_by(frame, name) %>%
summarise(mean = mean(value),
sd = sd(value)) %>%
mutate_if(is.double, round, digits = 3)
```

```
## # A tibble: 4 x 4
## # Groups: frame [2]
## frame name mean sd
## <dbl> <chr> <dbl> <dbl>
## 1 0 justify 2.80 0.849
## 2 0 skeptic 3.34 2.04
## 3 1 justify 2.94 1.01
## 4 1 skeptic 3.42 2.03
```

Let’s open **brms**.

`library(brms)`

Now fit the simple univariable model.

```
.1 <-
model7brm(data = disaster,
family = gaussian,
~ 1 + frame,
justify cores = 4,
file = "fits/model07.01")
```

`print(model7.1)`

```
## Family: gaussian
## Links: mu = identity; sigma = identity
## Formula: justify ~ 1 + frame
## Data: disaster (Number of observations: 211)
## Samples: 4 chains, each with iter = 2000; warmup = 1000; thin = 1;
## total post-warmup samples = 4000
##
## Population-Level Effects:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept 2.80 0.09 2.63 2.97 1.00 3774 2776
## frame 0.13 0.13 -0.11 0.38 1.00 3963 2724
##
## Family Specific Parameters:
## Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sigma 0.93 0.05 0.85 1.03 1.00 3620 2765
##
## Samples were drawn using sampling(NUTS). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).
```

The ‘Estimate’ (i.e., posterior mean) of the model intercept is the expected `justify`

value for when `frame`

is 0. The ‘Estimate’ for `frame`

is the expected difference when `frame`

is a 1. If all you care about is the posterior mean, you could execute

`fixef(model7.1)["Intercept", 1] + fixef(model7.1)["frame", 1]`

`## [1] 2.934172`

which matches up nicely with the equation on page 233. But this wouldn’t be very Bayesian of us. It’d be more satisfying if we had an expression of the uncertainty in the value. For that, we’ll follow our usual practice of extracting the posterior samples, making nicely-named vectors, and summarizing a bit.

```
<-
post posterior_samples(model7.1) %>%
mutate(when_x_is_0 = b_Intercept,
when_x_is_1 = b_Intercept + b_frame)
%>%
post pivot_longer(when_x_is_0:when_x_is_1) %>%
group_by(name) %>%
summarize(mean = mean(value),
sd = sd(value)) %>%
mutate_if(is.double, round, digits = 3)
```

```
## # A tibble: 2 x 3
## name mean sd
## <chr> <dbl> <dbl>
## 1 when_x_is_0 2.80 0.088
## 2 when_x_is_1 2.93 0.095
```

Hayes referenced a \(t\)-test and accompanying \(p\)-value in the lower part of page 233. We, of course, aren’t going to do that. But we do have the 95% intervals in our `print()`

output, above, which we can also look at with the `brms::posterior_interval()`

function.

`posterior_interval(model7.1)["b_frame", ]`

```
## 2.5% 97.5%
## -0.1104257 0.3847622
```

And we can always plot.

```
%>%
post ggplot(aes(x = b_frame)) +
geom_density(size = 0, fill = dutchmasters$little_street[1]) +
geom_vline(xintercept = 0, color = dutchmasters$little_street[11]) +
scale_x_continuous(breaks = c(-.3, 0, .6)) +
scale_y_continuous(NULL, breaks = NULL) +
+
theme_07 theme(legend.position = "none")
```

We’ll use the `update()`

function to hastily fit `model7.2`

and `model7.3`

.

```
.2 <-
model7update(model7.1,
newdata = disaster,
formula = justify ~ 1 + frame + skeptic,
cores = 4,
file = "fits/model07.02")
.3 <-
model7update(model7.1,
newdata = disaster,
formula = justify ~ 1 + frame + skeptic + frame:skeptic,
cores = 4,
file = "fits/model07.03")
```

Note our use of the `frame:skeptic`

syntax in `model7.3`

. With that syntax we didn’t need to make an interaction variable in the data by hand. The **brms** package just handled it for us. An alternative syntax would have been `frame*skeptic`

. But if you really wanted to make the interaction variable by hand, you’d do this.

```
<-
disaster %>%
disaster mutate(interaction_variable = frame * skeptic)
```

Once you have `interaction_variable`

in the data, you’d specify a model formula within the `brm()`

function like `formula = justify ~ 1 + frame + skeptic + interaction_variable`

. I’m not going to do that, here, but you can play around yourself if so inclined.

Here are the quick and dirty coefficient summaries for our two new models.

`posterior_summary(model7.2)`

```
## Estimate Est.Error Q2.5 Q97.5
## b_Intercept 2.1321750 0.12594154 1.8914550 2.3852145
## b_frame 0.1174051 0.11625236 -0.1110120 0.3459309
## b_skeptic 0.2003137 0.02839508 0.1440296 0.2554813
## sigma 0.8414353 0.04159860 0.7649706 0.9282171
## lp__ -265.6272830 1.42659274 -269.1917073 -263.8615255
```

`posterior_summary(model7.3)`

```
## Estimate Est.Error Q2.5 Q97.5
## b_Intercept 2.4516167 0.14668648 2.17005729 2.7526710
## b_frame -0.5645480 0.21392057 -0.97831354 -0.1458024
## b_skeptic 0.1046097 0.03775143 0.02660806 0.1775666
## b_frame:skeptic 0.2020772 0.05450541 0.09875548 0.3109934
## sigma 0.8181036 0.04037604 0.74317261 0.9020807
## lp__ -259.6144256 1.48807863 -263.07099629 -257.5339074
```

Just focusing on our primary model, `model7.3`

, here’s another way to look at the coefficients.

```
mcmc_plot(model7.3) +
theme_07
```

By default, the `brms::mcmc_plot()`

function returns coefficient plots which depict the parameters of a model by their posterior means (i.e., dots), 50% intervals (i.e., thick horizontal lines), and 95% intervals (i.e., thin horizontal lines). As `mcmc_plot()`

returns a **ggplot2** object, one can customize the theme and so on.

We’ll extract the \(R^2\) iterations in the usual way once for each model, and then combine them for a plot.

```
# for each of the three models, we create a separare R2 tibble
<-
r2_1 bayes_R2(model7.1, summary = F) %>%
data.frame()
<-
r2_2 bayes_R2(model7.2, summary = F) %>%
data.frame()
<-
r2_3 bayes_R2(model7.3, summary = F) %>%
data.frame()
# here we combine them into one tibble, indexed by `model`
<-
r2 %>%
r2_1 bind_rows(r2_2) %>%
bind_rows(r2_3) %>%
mutate(model = rep(str_c("7.", 1:3), each = 4000))
# now we plot
%>%
r2 ggplot(aes(x = R2, fill = model)) +
geom_density(size = 0, alpha = 2/3) +
scale_fill_manual(values = dutchmasters$little_street[c(3, 4, 8)] %>% as.character()) +
scale_x_continuous(expression(italic(R)^2~distribution), limits = 0:1) +
scale_y_continuous(NULL, breaks = NULL) +
theme_07
```

Here’s the \(\Delta R^2\) distribution for `model7.3`

minus `model7.2`

.

```
%>%
r2_2 bind_cols(r2_3) %>%
set_names(str_c("model7.", 2:3)) %>%
mutate(dif = model7.3 - model7.2) %>%
ggplot(aes(x = dif)) +
geom_density(color = "transparent",
fill = dutchmasters$little_street[9]) +
scale_y_continuous(NULL, breaks = NULL) +
labs(title = expression(The~Delta*italic(R)^2),
subtitle = "Like in the text, the posterior\nmean is about 0.05.",
x = NULL) +
theme_07
```

In addition to the \(R^2\), one can use information criteria to compare the models. Here we’ll use the LOO to compare all three.

```
.1 <- add_criterion(model7.1, "loo")
model7.2 <- add_criterion(model7.2, "loo")
model7.3 <- add_criterion(model7.3, "loo")
model7
loo_compare(model7.1, model7.2, model7.3) %>%
print(simplify = F)
```

```
## elpd_diff se_diff elpd_loo se_elpd_loo p_loo se_p_loo looic se_looic
## model7.3 0.0 0.0 -259.0 10.8 5.4 0.9 518.1 21.7
## model7.2 -5.6 4.1 -264.7 11.2 4.6 0.9 529.3 22.4
## model7.1 -27.1 9.5 -286.2 12.0 3.3 0.8 572.3 24.1
```

The LOO point estimate for both multivariable models were clearly lower than that for `model7.1`

. The point estimate for the moderation model, `model7.3`

, was within the double-digit range lower than that for `model7.2`

, which typically suggests better fit. But notice how wide the standard error was. There’s a lot of uncertainty, there. Hopefully this isn’t surprising. Our \(R^2\) difference was small and uncertain, too. We can also compare them with AIC-type model weighting, which you can learn more about starting at this point in this lecture or the related (2019) vignette, *Bayesian stacking and pseudo-BMA weights using the loo package* (Vehtari & Gabry, 2019). Here we’ll keep things simple and weight with the LOO.

```
model_weights(model7.1, model7.2, model7.3,
weights = "loo") %>%
round(digits = 3)
```

```
## model7.1 model7.2 model7.3
## 0.000 0.004 0.996
```

The `model_weights()`

results put almost all the relative weight on `model7.3`

. This doesn’t mean `model7.3`

is the “true model” or anything like that. It just suggests that it’s the better of the three with respect to the data.

Here are the results of the equations in the second half of page 237.

```
<- posterior_samples(model7.3)
post
%>%
post mutate(if_2 = b_frame + `b_frame:skeptic` * 2,
if_3.5 = b_frame + `b_frame:skeptic` * 3.5,
if_5 = b_frame + `b_frame:skeptic` * 5) %>%
pivot_longer(starts_with("if_")) %>%
group_by(name) %>%
summarise(mean = mean(value),
sd = sd(value)) %>%
mutate_if(is.double, round, digits = 3) %>%
mutate(name = str_remove(name, "if_")) %>%
rename(`conditional effect of X on Y when W =` = name)
```

```
## # A tibble: 3 x 3
## `conditional effect of X on Y when W =` mean sd
## <chr> <dbl> <dbl>
## 1 2 -0.16 0.134
## 2 3.5 0.143 0.113
## 3 5 0.446 0.145
```

### 7.2.1 Estimation using ~~PROCESS~~ brms.

Similar to what Hayes advertised with PROCESS, with our `formula = justify ~ 1 + frame + skeptic + frame:skeptic`

code in `model7.3`

, we didn’t need to hard code an interaction variable into the data. **brms** handled that for us.

### 7.2.2 Interpreting the regression coefficients.

When you add an interaction term into a model, such as \(x_1 \cdot x_2\), this is sometimes called a *higher order term*. The terms that made up the higher order term–\(x_1\) and \(x_2\), in this case–are correspondingly called the *lower order terms*. As Hayes pointed out, these lower order terms are conditional effects. Interpret them with care.

### 7.2.3 Variable scaling and the interpretation of \(b_1\) and \(b_3\).

Making the mean-centered version of our \(W\) variable, `skeptic`

, is a simple `mutate()`

operation. We’ll just call it `skeptic_c`

.

```
<-
disaster %>%
disaster mutate(skeptic_c = skeptic - mean(skeptic))
```

And here’s how we might fit the model.

```
.4 <-
model7update(model7.3,
newdata = disaster,
formula = justify ~ 1 + frame + skeptic_c + frame:skeptic_c,
cores = 4,
file = "fits/model07.04")
```

Here are the summaries of our fixed effects.

`fixef(model7.4)`

```
## Estimate Est.Error Q2.5 Q97.5
## Intercept 2.8073962 0.07701610 2.65744337 2.9562671
## frame 0.1133822 0.11460258 -0.10819717 0.3411017
## skeptic_c 0.1027887 0.03881743 0.02914400 0.1789173
## frame:skeptic_c 0.2040915 0.05496150 0.09173885 0.3086697
```

To practice, `frame`

and `skeptic_c`

are lower order terms and, as such, they are conditional effects. The higher order term is `frame:skeptic_c`

. Anyway, here are the \(R^2\) distributions for `model7.3`

and `model7.4`

. They’re the same within simulation variance.

`bayes_R2(model7.3) %>% round(digits = 3)`

```
## Estimate Est.Error Q2.5 Q97.5
## R2 0.249 0.045 0.158 0.334
```

`bayes_R2(model7.4) %>% round(digits = 3)`

```
## Estimate Est.Error Q2.5 Q97.5
## R2 0.249 0.044 0.161 0.335
```

If you’re bothered by the differences resulting from sampling variation, you might increase the number of HMC iterations from the 2,000-per-chain default. Doing so might look something like this.

```
.3 <-
model7update(model7.3,
chains = 4, cores = 4, warmup = 1000, iter = 10000)
.4 <-
model7update(model7.4,
chains = 4, cores = 4, warmup = 1000, iter = 10000)
```

Before we fit `model7.5`

, we’ll recode `frame`

to a -.5/.5 metric and name it `frame_.5`

.

```
<-
disaster %>%
disaster mutate(frame_.5 = ifelse(frame == 0, -.5, .5))
```

Time to fit `model7.5`

.

```
.5 <-
model7update(model7.4,
newdata = disaster,
formula = justify ~ 1 + frame_.5 + skeptic_c + frame_.5:skeptic_c,
cores = 4,
file = "fits/model07.05")
```

Our posterior summaries match up nicely with the output in Hayes’s Table 7.4.

`fixef(model7.5)`

```
## Estimate Est.Error Q2.5 Q97.5
## Intercept 2.8661008 0.05693406 2.75424820 2.9800025
## frame_.5 0.1162708 0.11283424 -0.10375391 0.3302863
## skeptic_c 0.2059356 0.02710684 0.15156950 0.2584645
## frame_.5:skeptic_c 0.2006627 0.05580013 0.09389709 0.3151269
```

Here’s a summary of the Bayesian \(R^2\).

`bayes_R2(model7.5)`

```
## Estimate Est.Error Q2.5 Q97.5
## R2 0.2493167 0.04306566 0.1624235 0.3303375
```

## 7.3 Visualizing moderation

A regression model with the product to two antecedent variables in an abstract mathematical representation of one’s data that can be harder to interpret than a model without such a produce. As described earlier, the coefficients for \(X\) and \(W\) are conditional effects that may not have any substantive interpretation, and the coefficient for \(XW\) is interpreted as a difference between differences that can be hard to make sense of without more information. (pp. 223–224)

This is why we plot. To get quick plots for the interaction effect with **brms**, you might use the `conditional_effects()`

function.

`conditional_effects(model7.3)`

By default, `conditional_effects()`

will show three levels of the variable on the right side of the interaction term. The formula in `model7.3`

was `justify ~ frame + skeptic + frame:skeptic`

, with `frame:skeptic`

as the interaction term and `skeptic`

making up the right hand side of the term. The three levels of `skeptic`

in the plot, above, are the mean \(\pm 1\) standard deviation. See the **brms** reference manual for details on the `conditional_effects()`

function (Bürkner, 2021e).

On page 244, Hayes discussed using the 16th, 50th, and 84th percentiles for the moderator variable. We can compute those with `quantile()`

.

`quantile(disaster$skeptic, probs = c(.16, .5, .84))`

```
## 16% 50% 84%
## 1.6 2.8 5.2
```

The first two columns in Hayes’s Table 7.5 contain the values he combined with the point estimates of his model to get the \(\widehat Y\) column. The way we’ll push those values through `model7.3`

’s posterior is with `brms::fitted()`

. As a preparatory step, we’ll put the predictor values in a data object, `nd`

.

```
(<-
nd crossing(frame = 0:1,
skeptic = quantile(disaster$skeptic,
probs = c(.16, .5, .84)))
)
```

```
## # A tibble: 6 x 2
## frame skeptic
## <int> <dbl>
## 1 0 1.6
## 2 0 2.8
## 3 0 5.2
## 4 1 1.6
## 5 1 2.8
## 6 1 5.2
```

Now we’ve got our `nd`

, we’ll get our posterior estimates for \(Y\) with `fitted()`

.

`fitted(model7.3, newdata = nd)`

```
## Estimate Est.Error Q2.5 Q97.5
## [1,] 2.618992 0.10059981 2.424506 2.822780
## [2,] 2.744524 0.07975865 2.589767 2.901872
## [3,] 2.995587 0.10562745 2.785480 3.196689
## [4,] 2.377768 0.10841435 2.164639 2.588615
## [5,] 2.745792 0.08410312 2.583013 2.908011
## [6,] 3.481841 0.10967583 3.262311 3.702851
```

When using the default `summary = TRUE`

settings in `fitted()`

, the function returns posterior means, \(SD\)s and 95% intervals for \(Y\) based on each row in the `nd`

data we specified in the `newdata = nd`

argument. You don’t have to name your newdata `nd`

or anything like that; it’s just my convention.

Here’ a quick plot of what those values imply.

```
fitted(model7.3, newdata = nd) %>%
data.frame() %>%
bind_cols(nd) %>%
ggplot(aes(x = skeptic)) +
geom_ribbon(aes(ymin = Q2.5, ymax = Q97.5, fill = frame %>% as.character()),
alpha = 1/3) +
geom_line(aes(y = Estimate, color = frame %>% as.character())) +
scale_fill_manual("frame",
values = dutchmasters$little_street[c(10, 5)] %>% as.character()) +
scale_color_manual("frame",
values = dutchmasters$little_street[c(10, 5)] %>% as.character()) +
theme_07
```

That plot is okay, but we can do better.

In order to plot the model-implied effects across the full range of `skeptic`

values presented in Figure 7.7, you need to change the range of those values in the `nd`

data. Also, although the effect is subtle in the above example, 95% intervals often follow a bowtie shape. In order to insure the contours of that shape are smooth, it’s often helpful to specify 30 or so evenly-spaced values in the variable on the \(x\)-axis, `skeptic`

in this case. We’ll employ the `seq()`

function for that and specify `length.out = 30`

. In addition, we add a few other flourishes to make our plot more closely resemble the one in the text.

Here’s our Figure 7.7.

```
<-
nd crossing(frame = 0:1,
skeptic = seq(from = 0, to = 7, length.out = 30))
fitted(model7.3, newdata = nd) %>%
as_tibble() %>%
bind_cols(nd) %>%
ggplot(aes(x = skeptic)) +
geom_ribbon(aes(ymin = Q2.5, ymax = Q97.5, fill = frame %>% as.character()),
alpha = 1/3) +
geom_line(aes(y = Estimate, color = frame %>% as.character())) +
scale_fill_manual("frame",
values = dutchmasters$little_street[c(10, 5)] %>% as.character()) +
scale_color_manual("frame",
values = dutchmasters$little_street[c(10, 5)] %>% as.character()) +
scale_x_continuous(breaks = 1:6) +
coord_cartesian(xlim = c(1, 6),
ylim = c(2, 4)) +
labs(x = expression(paste("Climate Change Skepticism (", italic("W"), ")")),
y = "Strength of Justification\nfor Withholding Aid") +
+
theme_07 theme(legend.position = "top")
```

Do you see that subtle bowtie shape?

## 7.4 Probing an interaction

As with other modeling contexts, there the results from our moderation models, like depicted in the last plot, carry uncertainty with them.

To deal with the uncertainty, it is common to follow up a test of interaction with a set of additional inferential tests to establish where in the distribution of the moderator \(X\) has an effect on \(Y\) that is different from zero and where it does not. This exercise is commonly known as “probing” an interaction, like you might squeeze an avocado or a mango in the produce section of the grocery store to assess its ripeness. The goal is to ascertain where in the distribution of the moderator \(X\) is related to \(Y\) and where it is not in an attempt to better discern the substantive interpretation of the interaction. In this section [we will cover] two approaches to probing an interaction. (p. 249)

As in earlier sections and chapters, our approach will largely follow Hayes’s, but we will not be emphasizing statistical tests, \(p\)-values, and so on. You’ll see.

### 7.4.1 The pick-a-point approach.

The pick-a-point approach (Bauer & Curran, 2005; Rogosa, 1980), sometimes called an

analysis of simple slopesor aspotlight analysis, is perhaps the most popular approach to probing an interaction and is described in most discussions of multiple regression with interactions (e.g., Aiken & West, 1991; Cohen et al., 2003; Darlington & Hayes, 2017; Andrew F. Hayes, 2005; Jaccard & Turrisi, 2003; Spiller et al., 2013). This procedure involves selecting a value or values of the moderator \(W\), calculating the conditional effect of \(X\) on \(Y (\theta_{X \rightarrow Y})\) at that value or values, and then conducting an inferential test or generating a confidence interval. (p. 249,emphasisin the original)

Though we will not be using inferential tests, we will use various ways of expressing the uncertainty in our Bayesian models, such as with 95% credible intervals and other ways of summarizing the posterior.

#### 7.4.1.1 The pick-a-point approach ~~implimented by regression centering~~ working directly with the posterior.

Yes, if you wanted to use the regression centering approach, you could do that in **brms**. Just center the necessary variables in the way Hayes described in the text, refit the model, and `summarize()`

. I suspect this would be particularly approachable for someone new to **R** and to the ins and outs of data wrangling. But I’m going leave that as an exercise for the interested reader.

Now that we’ve already got a posterior for our model, we can just either algebraically manipulate the vectors yielded by `posterior_samples()`

or push predictions through `fitted()`

. To give a sense, we’ll start off with the 16^{th} percentile for `skeptic`

. Recall we can get that with the `quantile()`

function.

`quantile(disaster$skeptic, probs = .16) `

```
## 16%
## 1.6
```

Now we just need to feed that value and different values of `frame`

into the posterior samples of the model coefficients. We then create a `difference`

score for the model-implied estimates given `frame`

is either 0 or 1 and then plot that `difference`

.

```
%>%
post # use the model formula
mutate(y_given_frame_0_skeptic_1.6 = b_Intercept + b_frame * 0 + b_skeptic * 1.6 + `b_frame:skeptic` * 0 * 1.6,
y_given_frame_1_skeptic_1.6 = b_Intercept + b_frame * 1 + b_skeptic * 1.6 + `b_frame:skeptic` * 1 * 1.6) %>%
# compute the difference
mutate(difference = y_given_frame_1_skeptic_1.6 - y_given_frame_0_skeptic_1.6) %>%
# plot!
ggplot(aes(x = difference)) +
geom_density(color = "transparent",
fill = dutchmasters$little_street[9]) +
geom_vline(xintercept = 0, color = dutchmasters$little_street[7], linetype = 2) +
scale_y_continuous(NULL, breaks = NULL) +
labs(subtitle = "The effect of frame on justify\ngiven skeptic = 1.6",
x = NULL) +
theme_07
```

Note how nicely that distribution corresponds to the output in the lower left corner of Hayes’s Figure 7.8. If we wanted the values for other values of `skeptic`

(e.g., 2.8 and 5.2 as in the text), we’d just rinse, wash, and repeat. A nice quality of this method is it requires you to work explicitly with the model formula. But it’s also clunky if you want to do this over many values. The `fitted()`

function offers an alternative approach.

Recall how the default `fitted()`

settings are to return summaries of a model’s \(Y\)-variable given values of the predictor variables. In the previous section we put our preferred `frame`

and `skeptic`

values into a data object named `nd`

and used the `newdata`

argument to push those values through `fitted()`

. Buy default, this yielded the typical posterior means, \(SD\)s, and 95% intervals for the predictions. However, if one sets `summary = F`

, the output will differ. First. Let’s revisit what `nd`

looks like.

```
(<-
nd crossing(frame = 0:1,
skeptic = quantile(disaster$skeptic,
probs = c(.16, .5, .84)))
)
```

```
## # A tibble: 6 x 2
## frame skeptic
## <int> <dbl>
## 1 0 1.6
## 2 0 2.8
## 3 0 5.2
## 4 1 1.6
## 5 1 2.8
## 6 1 5.2
```

Here’s what happens when we use `summary = F`

.

```
<-
f fitted(model7.3, newdata = nd, summary = F)
%>% str() f
```

`## num [1:4000, 1:6] 2.54 2.64 2.56 2.58 2.76 ...`

`%>% head() f `

```
## [,1] [,2] [,3] [,4] [,5] [,6]
## [1,] 2.544529 2.605685 2.727995 2.592793 2.922759 3.582689
## [2,] 2.642688 2.753064 2.973817 2.416421 2.770837 3.479670
## [3,] 2.560774 2.672912 2.897188 2.472279 2.835495 3.561927
## [4,] 2.580077 2.698479 2.935283 2.449071 2.807149 3.523305
## [5,] 2.755264 2.800871 2.892086 2.153414 2.598460 3.488551
## [6,] 2.601511 2.692850 2.875529 2.391907 2.748605 3.462002
```

With `summary = F`

, `fitted()`

returned a matrix of 4,000 rows (i.e., one for each posterior draw) and 6 vectors (i.e., one for each row in our `nd`

data). So now instead of summary information, we have a full expression of the uncertainty in terms of 4,000 draws. If you prefer working within the **tidyverse** and plotting with **ggplot2**, matrices aren’t the most useful data type. Let’s wrangle a bit.

```
<-
f %>%
f 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 = "_") %>%
pivot_wider(names_from = frame, values_from = value) %>%
mutate(difference = `1` - `0`,
# this isn't necessary, but will help with the facet labels
skeptic = str_c("skeptic = ", skeptic))
%>%
f head()
```

```
## # A tibble: 6 x 5
## iter skeptic `0` `1` difference
## <int> <chr> <dbl> <dbl> <dbl>
## 1 1 skeptic = 1.6 2.54 2.59 0.0483
## 2 1 skeptic = 2.8 2.61 2.92 0.317
## 3 1 skeptic = 5.2 2.73 3.58 0.855
## 4 2 skeptic = 1.6 2.64 2.42 -0.226
## 5 2 skeptic = 2.8 2.75 2.77 0.0178
## 6 2 skeptic = 5.2 2.97 3.48 0.506
```

Now we have our draws in a nice structure, we’re ready to plot.

```
%>%
f ggplot(aes(x = difference)) +
geom_density(color = "transparent", fill = dutchmasters$little_street[9]) +
geom_vline(xintercept = 0, color = dutchmasters$little_street[7], linetype = 2) +
scale_y_continuous(NULL, breaks = NULL) +
labs(subtitle = "The effect of frame on justify given three different values of skeptic",
x = NULL) +
+
theme_07 facet_wrap(~ skeptic)
```

If you prefer summary information instead of plots, you might just use `tidybayes::median_qi()`

.

```
library(tidybayes)
%>%
f group_by(skeptic) %>%
median_qi(difference) %>%
mutate_if(is.double, round, digits = 3)
```

```
## # A tibble: 3 x 7
## skeptic difference .lower .upper .width .point .interval
## <chr> <dbl> <dbl> <dbl> <dbl> <chr> <chr>
## 1 skeptic = 1.6 -0.243 -0.53 0.051 0.95 median qi
## 2 skeptic = 2.8 -0.003 -0.228 0.229 0.95 median qi
## 3 skeptic = 5.2 0.486 0.193 0.786 0.95 median qi
```

### 7.4.2 The Johnson-Neyman technique.

The JN technique generalizes this approach over many values of \(W\) (i.e., `skeptic`

in this example) in order to get a sense of the trend and summarize regions of the trend in terms of \(p\)-value thresholds. Since we’re emphasizing modeling and deemphasizing null-hypothesis testing in this project, I’ll show a Bayesian version of the approach without the \(p\)-values.

#### 7.4.2.1 Implementation in ~~PROCESS~~ brms.

Since Figure 7.9 had `skeptic`

values ranging from 1 to 6 with ticks on the 0.5’s, we’ll use a similar approach for our version. We will display posterior samples with `fitted()`

for `skeptic`

values ranging from .5 to 6.5, one for each 0.5–13 in total. But since we have two levels of `frame`

(i.e., 0 and 1), that really gives us 26. And we don’t just want 26 summaries; we want full posterior distributions for each of those 26.

We’ve got a lot of moving parts in the code, below. To help make sure everything adds up, we’ll save several important values as **R** objects.

```
# this number comes from the total number of post-warmup posterior iterations from `brm()`
<- 4000
iter
# we're choosing 13 in this example to follow some of the sensibilities in Figure 7.9. You'll see.
<- 13
n_skeptic_values
# as before, we'll make `nd` to feed in to `fitted()`
<-
nd crossing(frame = 0:1,
skeptic = seq(from = .5, to = 6.5, length.out = n_skeptic_values))
# after the initial `fitted()` action, we need a few steps to wrangle the data into a useful format
<-
f fitted(model7.3, newdata = nd, summary = F) %>%
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 = "_") %>%
pivot_wider(names_from = frame, values_from = value) %>%
mutate(difference = `1` - `0`,
skeptic = as.double(skeptic))
# finally, here's the plot
%>%
f ggplot(aes(x = skeptic %>% as.character, y = difference)) +
geom_hline(yintercept = 0, color = dutchmasters$little_street[7]) +
geom_violin(size = 0, fill = dutchmasters$little_street[6]) +
stat_summary(fun = median,
fun.min = function(x){quantile(x, probs = .025)},
fun.max = function(x){quantile(x, probs = .975)},
color = dutchmasters$little_street[5]) +
labs(x = expression(paste("Climate Change Skepticism (", italic(W), ")")),
y = expression(atop(theta[paste(italic(X), " on ", italic(Y))], paste("Conditional Effect of Disaster Frame")))) +
theme_07
```

[Note. I got the `atop()`

trick for the label for the \(y\)-axis from Drew Steen’s answer to this stackoverflow question.]

This isn’t quite our version of Figure 7.9, but I’m hoping it’ll add some pedagogical value for what we’re doing. Since we specified `summary = F`

within `fitted()`

, we got full posterior distributions for each of our 26 conditions. *Because Figure 7.9 is all about differences between each frame pair across the various values of skeptic*, we needed to make a

`difference`

score for each pair; this is what we did with the last `mutate()`

line before the plot code. This initial version of the plot shows the full posterior distribution for each `difference`

score. The posteriors are depicted with violin plots, which are density plots set on their side and symmetrically reflected as if by a mirror to give a pleasing leaf- or violin-like shape (though beware). The light dots and vertical lines are the posterior medians and 95% intervals for each.Going from left to right, it appears we have a clearly emerging trend. We can more simply express the trend by summarizing each posterior with medians and 95% intervals.

```
%>%
f group_by(skeptic) %>%
median_qi(difference) %>%
ggplot(aes(x = skeptic)) +
geom_hline(yintercept = 0, color = dutchmasters$little_street[7]) +
geom_vline(xintercept = c(1.171, 3.934), color = dutchmasters$little_street[7]) +
geom_ribbon(aes(ymin = .lower, ymax = .upper),
fill = dutchmasters$little_street[5],
alpha = 1/2) +
geom_line(aes(y = difference),
color = dutchmasters$little_street[5], size = 1) +
scale_x_continuous(breaks = 1:6) +
coord_cartesian(xlim = c(1, 6),
ylim = c(-1, 1.5)) +
labs(x = expression(paste("Climate Change Skepticism (", italic(W), ")")),
y = expression(atop(theta[paste(italic(X), " on ", italic(Y))], paste("Conditional Effect of Disaster Frame")))) +
theme_07
```

Notice how the contour boundaries of the 95% intervals are a little clunky. That’s because our bowtie-shape is based on only 13 \(x\)-axis values. If you wanted a smoother shape, you’d specify more `skeptic`

values in the data object you feed into `fitted()`

’s `newdata`

argument. For linear effects, I find 30 or so usually does it.

Anyway, I got the values for the two vertical lines directly out of the text. It’s not clear to me how one might elegantly determine those values within the paradigm we’ve been using. But that leads to an important digression. The two vertical lines are quite \(p\)-value centric. They are an attempt to separate the \(x\)-axis into areas where the `difference`

trend either is or is not statistically-significantly different from zero. That is, we’re dichotomizing–or “trichotomizing,” depending on how you look at it–a continuous phenomenon. This is somewhat at odds with the sensibilities of the Bayesians associated with Stan and **brms** (e.g., here).

On page 259, Hayes wrote: “Although the JN technique eliminates the need to select arbitrary values of \(W\) when probing an interaction, it does not eliminate your need to keep your brain turned into the task and thinking critically about the answer the method gives you.” I think this is valuable advice, particularly when working within the Bayesian paradigm. Our version of Figure 7.9 gives some interesting insights into the moderation model, `model7.3`

. I’m just not so sure I’d want to encourage people to interpret a continuous phenomenon by heuristically dividing it into discrete regions.

## 7.5 The difference between testing for moderation and probing it

This is another section where the NHST-type paradigm contrasts with many within the contemporary Bayesian paradigm. E.g., Hayes opened the section with: “We test for evidence of moderation when we want to know whether the relationship between \(X\) and \(Y\) varies systematically as a function of a proposed moderator \(W\).” His use of “whether” suggests we are talking about a binary answer–either there is an effect or there isn’t. But, as Gelman (2015) has argued, the default presumption in social science [and warning, I’m a psychologist and thus biased towards thinking in terms of social science] is that treatment effects–and more generally, causal effects–vary across contexts^{3} . As such, asking “whether” there’s a difference or an interaction effect isn’t really the right question. Rather, we should presume variation at the outset and ask instead what the magnitude of that variation is and how much accounting for it matters for our given purposes. If the variation–read *interaction effect*–is tiny and of little theoretical interest, perhaps we might just ignore it and not include it in the model. Alternatively, if the variation is large or of theoretical interest, we might should include it in the model regardless of statistical significance.

Another way into this topic is posterior predictive checking. We’ve already done a bit of this in previous chapters. The basic idea, recall, is that better models should give us a better sense of the patterns in the data. In the plot below, we continue to show the interaction effect with two regression lines, but this time we separate them into their own panels by `frame`

. In addition, we add the original data which we also separate and color code by `frame`

.

```
<-
nd crossing(frame = 0:1,
skeptic = seq(from = 0, to = 10, length.out = 30))
fitted(model7.3, newdata = nd) %>%
data.frame() %>%
bind_cols(nd) %>%
mutate(frame = frame %>% as.character()) %>%
ggplot(aes(x = skeptic)) +
geom_ribbon(aes(ymin = Q2.5, ymax = Q97.5, fill = frame),
alpha = 1/3) +
geom_line(aes(y = Estimate, color = frame)) +
geom_point(data = disaster %>% mutate(frame = frame %>% as.character()),
aes(y = justify, color = frame),
alpha = 3/4) +
scale_fill_manual("frame",
values = dutchmasters$little_street[c(10, 5)] %>% as.character()) +
scale_color_manual("frame",
values = dutchmasters$little_street[c(10, 5)] %>% as.character()) +
scale_x_continuous(breaks = c(1, 9)) +
coord_cartesian(xlim = c(1, 9)) +
labs(title = "model 7.3, the interaction model",
x = expression(paste("Climate Change Skepticism (", italic("W"), ")")),
y = "Strength of Justification\nfor Withholding Aid") +
+
theme_07 theme(legend.position = "top") +
facet_wrap(~ frame)
```

When we separate out the data this way, it really does appear that when `frame == 1`

, the `justify`

values do increase as the `skeptic`

values increase, but not so much when `frame == 0`

. We can use the same plotting approach, but this time with the results from the non-interaction multivariable model, `model7.2`

.

```
fitted(model7.2, newdata = nd) %>%
data.frame() %>%
bind_cols(nd) %>%
mutate(frame = frame %>% as.character()) %>%
ggplot(aes(x = skeptic)) +
geom_ribbon(aes(ymin = Q2.5, ymax = Q97.5, fill = frame),
alpha = 1/3) +
geom_line(aes(y = Estimate, color = frame %>% as.character())) +
geom_point(data = disaster %>% mutate(frame = frame %>% as.character()),
aes(y = justify, color = frame),
alpha = 3/4) +
scale_fill_manual("frame",
values = dutchmasters$little_street[c(10, 5)] %>% as.character()) +
scale_color_manual("frame",
values = dutchmasters$little_street[c(10, 5)] %>% as.character()) +
scale_x_continuous(breaks = c(1, 9)) +
coord_cartesian(xlim = c(1, 9)) +
labs(title = "model 7.2, the multivariable model",
x = expression(paste("Climate Change Skepticism (", italic("W"), ")")),
y = "Strength of Justification\nfor Withholding Aid") +
+
theme_07 theme(legend.position = "top") +
facet_wrap(~ frame)
```

This time when we allowed the intercept but not the slope to vary by `frame`

, it appears the regression lines are missing part of the story. They look okay, but it appears that the red line on the left is sloping up to quickly and that the cream line on the right isn’t sloping steeply enough. We have missed an insight.

Now imagine scenarios in which the differences by `frame`

are more or less pronounced. Imagine those scenarios fall along a continuum. It’s not so much that you can say with certainty where on such a continuous an interaction effect would exist or not, but rather, such a continuum suggests it would appear more or less important, of greater or smaller magnitude. It’s not that the effect exists or is non-zero. It’s that it is orderly enough and of a large enough magnitude, and perhaps of theoretical interest, that it appears to matter in terms of explaining the data.

**And none of this is to serve as a harsh criticism of Andrew Hayes**. His text is a fine effort to teach mediation and moderation from a frequentist OLS perspective. I have benefited tremendously from his work. Yet I’d also like to connect his work to some other sensibilities.

Building further, consider this sentence from the text:

Rather, probing moderation involves ascertaining whether the conditional effect of \(X\) on \(Y\) is different from zero at certain specified values of \(W\) (if using the pick-a-point approach) or exploring where in the distribution of \(W\) the conditional effect of \(X\) on \(Y\) transitions between statistically significant and non-significant (if using the Johnson-Neyman technique). (pp. 259–260)

From an NHST/frequentist perspective, this makes clear sense. But we’re dealing with an entire posterior distribution. Consider again a figure from above.

```
<-
nd crossing(frame = 0:1,
skeptic = quantile(disaster$skeptic,
probs = c(.16, .5, .84)))
fitted(model7.3, newdata = nd, summary = F) %>%
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 = "_") %>%
pivot_wider(names_from = frame, values_from = value) %>%
mutate(difference = `1` - `0`,
# this isn't necessary, but will help with the facet labels
skeptic = str_c("skeptic = ", skeptic)) %>%
ggplot(aes(x = difference)) +
geom_density(color = "transparent", fill = dutchmasters$little_street[9]) +
scale_y_continuous(NULL, breaks = NULL) +
labs(subtitle = "The effect of frame on justify given three different values of skeptic",
x = NULL) +
+
theme_07 facet_wrap(~ skeptic)
```

With the pick pick-a-point approach one could fixate on whether zero was a credible value within the posterior, given a particular `skeptic`

value. And yet zero is just one point in the parameter space. One might also focus on the whole shapes of the posteriors of these three `skeptic`

values. You could focus on where the most credible values (i.e., those at and around their peaks) are on the number line (i.e., the effect sizes) and you could also focus on the relative widths of the distributions (i.e., the precision with which the effect sizes are estimated). These sensibilities can apply to the JN technique, as well. Sure, we might be interested in how credible zero is. But there’s a lot more to notice, too.

Now consider a modified version of our JN technique plot, from above.

```
%>%
f ggplot(aes(x = skeptic, y = difference)) +
# the main action is in `tidybayes::stat_lineribbon()`
stat_lineribbon(.width = c(.1, .2, .3, .4, .5, .6, .7, .8, .9, .99),
fill = dutchmasters$little_street[5],
alpha = 1/5, size = 0) +
scale_x_continuous(breaks = 1:6) +
coord_cartesian(xlim = c(1, 6),
ylim = c(-1, 1.5)) +
labs(x = expression(paste("Climate Change Skepticism (", italic(W), ")")),
y = expression(atop(theta[paste(italic(X), " on ", italic(Y))], paste("Conditional Effect of Disaster Frame")))) +
theme_07
```

This time we emphasized the shape of the posterior with stacked semitransparent 10, 20, 30, 40, 50, 60, 70, 80, 90, and 99% intervals. We also deemphasized the central tendency–our analogue to the OLS point estimate–by removing the median line. Yes, one could focus on where the 95% intervals cross zero. And yes one could request we emphasize central tendency. But such focuses miss a lot of information about the shape–the entire smooth, seamless distribution of credible values.

I suppose you could consider this our version of Figure 7.10.

## 7.6 Artificial categorization and subgroups

There are multiple ways to dichotomize the data by `skeptic`

. A quick simple way is to use `if_else()`

to make a `skeptic_hi`

dummy.

```
<-
disaster %>%
disaster mutate(skeptic_hi = if_else(skeptic >= mean(skeptic), 1, 0))
```

With our dummy in hand, we’re ready to fit the two models.

```
.6 <-
model7brm(data = disaster %>% filter(skeptic_hi == 0),
family = gaussian,
~ 1 + frame,
justify cores = 4,
file = "fits/model07.06")
.7 <-
model7update(model7.6, newdata = disaster %>% filter(skeptic_hi == 1),
cores = 4,
file = "fits/model07.07")
```

Behold the coefficient summaries.

`fixef(model7.6) %>% round(digits = 3)`

```
## Estimate Est.Error Q2.5 Q97.5
## Intercept 2.622 0.096 2.434 2.806
## frame -0.104 0.141 -0.382 0.173
```

`fixef(model7.7) %>% round(digits = 3)`

```
## Estimate Est.Error Q2.5 Q97.5
## Intercept 3.068 0.143 2.790 3.348
## frame 0.472 0.205 0.077 0.871
```

You can use `fitted()`

to get the posterior means and other summaries for the two `frame`

groups, by model.

```
fitted(model7.6,
newdata = tibble(frame = 0:1)) %>%
round(digits = 3)
```

```
## Estimate Est.Error Q2.5 Q97.5
## [1,] 2.622 0.096 2.434 2.806
## [2,] 2.518 0.103 2.315 2.714
```

```
fitted(model7.7,
newdata = tibble(frame = 0:1)) %>%
round(digits = 3)
```

```
## Estimate Est.Error Q2.5 Q97.5
## [1,] 3.068 0.143 2.790 3.348
## [2,] 3.541 0.147 3.243 3.827
```

Do note that though **brms** ‘Est.Error’ is the posterior \(SD\) for the coefficient, it is not the same thing as descriptive statistic \(SD\) of a subset of the data. Thus, although our ‘Estimates’ correspond nicely to the mean values Hayes reported in the middle of page 264, his \(SD\)s will not match up with our ‘Est.Error’ values, and nor should they.

Anyway, our results don’t yield \(t\)-tests. But you don’t need those anyway. We’re working within the Bayesian regression paradigm! But if you’re really interested in the sub-model-implied differences between the two levels of `frame`

by `skeptic_hi`

subgroup, all you need is the `frame`

coefficient of `model7.6`

and `model7.7`

. Here we’ll use `bind_rows()`

to combine their posterior samples and then plot.

```
bind_rows(
posterior_samples(model7.6) %>% select(b_frame),
posterior_samples(model7.7) %>% select(b_frame)
%>%
) mutate(model = rep(c("model7.6", "model7.7"), each = 4000)) %>%
ggplot(aes(x = b_frame, fill = model)) +
geom_density(size = 0, alpha = .8) +
scale_fill_manual(NULL,
values = dutchmasters$little_street[c(1, 9)] %>% as.character()) +
scale_y_continuous(NULL, breaks = NULL) +
labs(subtitle = "The difference score distributions between frame\nlevels, color coded by mean-split skeptic",
x = NULL) +
theme_07
```

As within the frequentist paradigm, please don’t mean split as a Bayesian. When possible, use all available data and use the regression formula to model theoretically-meaningful variables in your analyses.

## Session info

`sessionInfo()`

```
## R version 4.0.4 (2021-02-15)
## Platform: x86_64-apple-darwin17.0 (64-bit)
## Running under: macOS Catalina 10.15.7
##
## Matrix products: default
## BLAS: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRblas.dylib
## LAPACK: /Library/Frameworks/R.framework/Versions/4.0/Resources/lib/libRlapack.dylib
##
## locale:
## [1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
##
## attached base packages:
## [1] stats graphics grDevices utils datasets methods base
##
## other attached packages:
## [1] tidybayes_3.0.0 brms_2.15.0 Rcpp_1.0.6 directlabels_2021.1.13
## [5] ochRe_1.0.0 dutchmasters_0.1.0 forcats_0.5.1 stringr_1.4.0
## [9] dplyr_1.0.6 purrr_0.3.4 readr_1.4.0 tidyr_1.1.3
## [13] tibble_3.1.2 ggplot2_3.3.5 tidyverse_1.3.1
##
## loaded via a namespace (and not attached):
## [1] readxl_1.3.1 backports_1.2.1 plyr_1.8.6 igraph_1.2.6 svUnit_1.0.3
## [6] sp_1.4-4 splines_4.0.4 crosstalk_1.1.0.1 TH.data_1.0-10 rstantools_2.1.1
## [11] inline_0.3.17 digest_0.6.27 htmltools_0.5.1.1 rsconnect_0.8.16 fansi_0.4.2
## [16] checkmate_2.0.0 magrittr_2.0.1 modelr_0.1.8 RcppParallel_5.0.2 matrixStats_0.57.0
## [21] sandwich_3.0-0 xts_0.12.1 prettyunits_1.1.1 colorspace_2.0-0 rvest_1.0.1
## [26] ggdist_3.0.0 haven_2.3.1 xfun_0.23 callr_3.7.0 crayon_1.4.1
## [31] jsonlite_1.7.2 lme4_1.1-25 survival_3.2-10 zoo_1.8-8 glue_1.4.2
## [36] gtable_0.3.0 emmeans_1.5.2-1 V8_3.4.0 distributional_0.2.2 pkgbuild_1.2.0
## [41] rstan_2.21.2 abind_1.4-5 scales_1.1.1 mvtnorm_1.1-1 DBI_1.1.0
## [46] miniUI_0.1.1.1 viridisLite_0.4.0 xtable_1.8-4 stats4_4.0.4 StanHeaders_2.21.0-7
## [51] DT_0.16 htmlwidgets_1.5.3 httr_1.4.2 threejs_0.3.3 arrayhelpers_1.1-0
## [56] posterior_1.0.1 ellipsis_0.3.2 pkgconfig_2.0.3 loo_2.4.1 farver_2.1.0
## [61] sass_0.3.1 dbplyr_2.1.1 utf8_1.2.1 tidyselect_1.1.1 labeling_0.4.2
## [66] rlang_0.4.11 reshape2_1.4.4 later_1.2.0 munsell_0.5.0 cellranger_1.1.0
## [71] tools_4.0.4 cli_3.0.1 generics_0.1.0 broom_0.7.6 ggridges_0.5.3
## [76] evaluate_0.14 fastmap_1.1.0 processx_3.5.2 knitr_1.33 fs_1.5.0
## [81] nlme_3.1-152 mime_0.10 projpred_2.0.2 xml2_1.3.2 compiler_4.0.4
## [86] bayesplot_1.8.0 shinythemes_1.1.2 rstudioapi_0.13 gamm4_0.2-6 curl_4.3
## [91] reprex_2.0.0 statmod_1.4.35 bslib_0.2.4 stringi_1.6.2 highr_0.9
## [96] ps_1.6.0 Brobdingnag_1.2-6 lattice_0.20-41 Matrix_1.3-2 nloptr_1.2.2.2
## [101] markdown_1.1 tensorA_0.36.2 shinyjs_2.0.0 vctrs_0.3.8 pillar_1.6.1
## [106] lifecycle_1.0.0 jquerylib_0.1.4 bridgesampling_1.0-0 estimability_1.3 raster_3.4-5
## [111] httpuv_1.6.0 R6_2.5.0 bookdown_0.22 promises_1.2.0.1 gridExtra_2.3
## [116] codetools_0.2-18 boot_1.3-26 colourpicker_1.1.0 MASS_7.3-53 gtools_3.8.2
## [121] assertthat_0.2.1 withr_2.4.2 shinystan_2.5.0 multcomp_1.4-16 mgcv_1.8-33
## [126] parallel_4.0.4 hms_1.1.0 quadprog_1.5-8 grid_4.0.4 coda_0.19-4
## [131] minqa_1.2.4 rmarkdown_2.8 shiny_1.6.0 lubridate_1.7.10 base64enc_0.1-3
## [136] dygraphs_1.1.1.6
```

## Footnote

### References

*Multiple regression: Testing and interpreting interactions*. SAGE.

*ochRe: Australia-themed colour palettes*[Manual]. https://github.com/ropenscilabs/ochRe

*Multivariate Behavioral Research*,

*40*(3), 373–400. https://doi.org/10.1207/s15327906mbr4003_5

*Journal of Experimental Psychology: General*,

*148*(4), 601–618. https://doi.org/10.1037/xge0000558

*brms reference manual, Version 2.15.0*. https://CRAN.R-project.org/package=brms/brms.pdf

*Social Psychological and Personality Science*,

*7*(1), 13–20. https://doi.org/10.1177/1948550615590448

*Applied multiple regression/correlation analysis for the behavioral sciences*(3rd edition). Lawrence Erlbaum Associates. https://doi.org/10.4324/9780203774441

*Regression analysis and linear models: Concepts, applications, and implementation*. Guilford Press. https://www.guilford.com/books/Regression-Analysis-and-Linear-Models/Darlington-Hayes/9781462521135/reviews

*Journal of Management*,

*41*(2), 632–643. https://doi.org/10.1177/0149206314525208

*Handbook of research methods for studying daily life*(pp. 43–61). The Guilford Press. https://www.guilford.com/books/Handbook-of-Research-Methods-for-Studying-Daily-Life/Mehl-Conner/9781462513055

*Statistical methods for communication science*. Routledge. https://doi.org/10.4324/9781410613707

*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

*Directlabels: Direct labels for multicolor plots*[Manual]. https://CRAN.R-project.org/package=directlabels

*Interaction effects in multiple regression*(2nd edition). Sage Publications. https://dx.doi.org/10.4135/9781412984522

*Mastering software development in {}R{}*. https://github.com/rdpeng/RProgDA

*Psychological Bulletin*,

*88*(2), 307–321. https://doi.org/10.1037/0033-2909.88.2.307

*Journal of Marketing Research*,

*50*(2), 277–288. https://doi.org/10.1509/jmr.12.0420

*dutchmasters*[Manual]. https://github.com/EdwinTh/dutchmasters

*Bayesian stacking and pseudo-BMA weights using the loo package*. https://CRAN.R-project.org/package=loo/vignettes/loo2-weights.html

*The little street*.

If you’re an experimental psychologist, check out the great (2019) paper by Bolger, Zee, Rossignac-Milon, and Hassin,

*Causal processes in psychology are heterogeneous*. The rest of you social scientists aren’t off the hook. Check out Ellen Hamaker’s excellent (2012) chapter,*Why researchers should think “within-person”: A paradigmatic rationale*. Both works suggest researchers might do well to switch out their fixed-effects models for multilevel models. Does this mean the methods we’ve been covering are of no use? No. But don’t stop here, friends. Keep learning!↩︎