## 13.5 Relative conditional direct effects

In order to get the $$R^2$$ difference distribution analogous to the change in $$R^2$$ $$F$$-test Hayes discussed on pages 495–496, we’ll have to first refit the model without the interaction for the $$Y$$ criterion, liking.

m_model <- bf(respappr ~ 1 + D1 + D2 + sexism + D1:sexism + D2:sexism)
y_model <- bf(liking   ~ 1 + D1 + D2 + respappr + sexism)

model3 <-
brm(data = protest, family = gaussian,
m_model + y_model + set_rescor(FALSE),
chains = 4, cores = 4)

Here’s the $$\Delta R^2$$ density.

R2s <-
bayes_R2(model1, resp = "liking", summary = F) %>%
as_tibble() %>%
rename(model1 = R2_liking) %>%
bind_cols(
bayes_R2(model3, resp = "liking", summary = F) %>%
as_tibble() %>%
rename(model3 = R2_liking)
) %>%
mutate(difference = model1 - model3)

R2s %>%
ggplot(aes(x = difference, y = 0)) +

geom_halfeyeh(point_interval = median_qi, .prob = c(0.95, 0.5),
fill = "grey50", color = "white") +
scale_x_continuous(breaks = c(-.5, median(R2s$difference) %>% round(2), .5)) + scale_y_continuous(NULL, breaks = NULL) + coord_cartesian(xlim = c(-.5, .5)) + xlab(expression(paste(Delta, italic(R)^2))) + theme_black() + theme(panel.grid = element_blank()) We’ll also compare the models by their information criteria. loo(model1, model3) ## LOOIC SE ## model1 760.10 29.56 ## model3 761.61 31.10 ## model1 - model3 -1.51 5.47 waic(model1, model3) ## WAIC SE ## model1 759.74 29.47 ## model3 761.39 31.05 ## model1 - model3 -1.65 5.46 As when we went through these steps for resp = "respappr", above, the Bayesian $$R^2$$, the LOO-CV, and the WAIC all suggest there’s little difference between the two models with respect to predictive utility. In such a case, I’d lean on theory to choose between them. If inclined, one could also do Bayesian model averaging. Our approach to plotting the relative conditional direct effects will mirror what we did for the relative conditional indirect effects, above. Here are the brm() parameters that correspond to the parameter names of Hayes’s notation. • $$c_{1}$$ = b_liking_D1 • $$c_{2}$$ = b_liking_D2 • $$c_{4}$$ = b_liking_D1:sexism • $$c_{5}$$ = b_liking_D2:sexism With all clear, we’re off to the races. # c1 + c4W D1_function <- function(w){ post$b_liking_D1 + post$b_liking_D1:sexism*w } # c2 + c5W D2_function <- function(w){ post$b_liking_D2 + post$b_liking_D2:sexism*w } rcde_tibble <- tibble(sexism = seq(from = 3.5, to = 6.5, length.out = 30)) %>% group_by(sexism) %>% mutate(Protest vs. No Protest = map(sexism, D1_function), Collective vs. Individual Protest = map(sexism, D2_function)) %>% unnest() %>% ungroup() %>% mutate(iter = rep(1:4000, times = 30)) %>% gather(direct effect, value, -sexism, -iter) %>% mutate(direct effect = factor(direct effect, levels = c("Protest vs. No Protest", "Collective vs. Individual Protest"))) head(rcde_tibble) ## # A tibble: 6 x 4 ## sexism iter direct effect value ## <dbl> <int> <fct> <dbl> ## 1 3.5 1 Protest vs. No Protest -0.856 ## 2 3.5 2 Protest vs. No Protest -0.482 ## 3 3.5 3 Protest vs. No Protest -1.24 ## 4 3.5 4 Protest vs. No Protest -1.23 ## 5 3.5 5 Protest vs. No Protest -1.06 ## 6 3.5 6 Protest vs. No Protest -0.663 Here is our variant of Figure 13.4, with respect to the relative conditional direct effects. rcde_tibble %>% group_by(direct effect, sexism) %>% summarize(median = median(value), ll = quantile(value, probs = .025), ul = quantile(value, probs = .975)) %>% ggplot(aes(x = sexism, group = direct effect)) + geom_ribbon(aes(ymin = ll, ymax = ul), color = "white", fill = "transparent", linetype = 3) + geom_line(aes(y = median), color = "white") + coord_cartesian(xlim = 4:6, ylim = c(-.6, .8)) + labs(x = expression(paste("Perceived Pervasiveness of Sex Discrimination in Society (", italic(W), ")")), y = "Relative Conditional Effect on Liking") + theme_black() + theme(panel.grid = element_blank(), legend.position = "none") + facet_grid(~ direct effect) Holy smokes, them are some wide 95% CIs! No wonder the information criteria and $$R^2$$ comparisons were so uninspiring. Notice that the y-axis is on the parameter space. In Hayes’s Figure 13.5, the y-axis is on the liking space, instead. When we want things in the parameter space, we work with the output of posterior_samples(); when we want them in the criterion space, we use fitted(). # we need new nd data nd <- tibble(D1 = rep(c(1/3, -2/3, 1/3), each = 30), D2 = rep(c(1/2, 0, -1/2), each = 30), respappr = mean(protest$respappr),
sexism = seq(from = 3.5, to = 6.5, length.out = 30) %>% rep(., times = 3))

# we feed nd into fitted()
model1_fitted <-
fitted(model1,
newdata = nd,
resp = "liking",
summary = T) %>%
as_tibble() %>%
bind_cols(nd) %>%
mutate(condition = ifelse(D2 == 0, "No Protest",
ifelse(D2 == -1/2, "Individual Protest", "Collective Protest"))) %>%
mutate(condition = factor(condition, levels = c("No Protest", "Individual Protest", "Collective Protest")))

model1_fitted %>%
ggplot(aes(x = sexism, group = condition)) +
geom_ribbon(aes(ymin = Q2.5, ymax = Q97.5),
linetype = 3, color = "white", fill = "transparent") +
geom_line(aes(y = Estimate), color = "white") +
geom_point(data = protest, aes(x = sexism, y = liking),
color = "red", size = 2/3) +
coord_cartesian(xlim = 4:6,
ylim = 4:7) +
labs(x = expression(paste("Perceived Pervasiveness of Sex Discrimination in Society (", italic(W), ")")),
y = expression(paste("Evaluation of the Attorney (", italic(Y), ")"))) +
theme_black() +
theme(panel.grid = element_blank()) +
facet_wrap(~condition) We expanded the range of the y-axis, a bit, to show more of that data (and there’s even more data outside of our expanded range). Also note how after doing so and after including the 95% CI bands, the crossing regression line effect in Hayes’s Figure 13.5 isn’t as impressive looking any more.

On pages 497–498, Hayes discussed more omnibus $$F$$-tests. Much like with the $$M$$ criterion, we won’t come up with Bayesian $$F$$-tests, but we might go ahead and make pairwise comparisons at the three percentiles Hayes prefers.

# we need new nd data
nd <-
tibble(D1 = rep(c(1/3, -2/3, 1/3), each = 3),
D2 = rep(c(1/2, 0, -1/2), each = 3),
respappr = mean(protest\$respappr),
sexism = rep(c(4.250, 5.120, 5.896), times = 3))

# this tie we'll use summary = F
model1_fitted <-
fitted(model1,
newdata = nd,
resp = "liking",
summary = F) %>%
as_tibble() %>%
gather() %>%
mutate(condition = rep(c("Collective Protest", "No Protest", "Individual Protest"),
each = 3*4000),
sexism = rep(c(4.250, 5.120, 5.896), times = 3) %>% rep(., each = 4000),
iter = rep(1:4000, times = 9)) %>%
select(-key) %>%
spread(key = condition, value = value) %>%
mutate(Individual Protest - No Protest = Individual Protest - No Protest,
Collective Protest - No Protest = Collective Protest - No Protest,
Collective Protest - Individual Protest = Collective Protest - Individual Protest)

# a tiny bit more wrangling and we're ready to plot the difference distributions
model1_fitted %>%
select(sexism, contains("-")) %>%
gather(key, value, -sexism) %>%

ggplot(aes(x = value)) +
geom_halfeyeh(aes(y = 0), fill = "grey50", color = "white",
point_interval = median_qi, .prob = 0.95) +
geom_vline(xintercept = 0, color = "grey25", linetype = 2) +
scale_y_continuous(NULL, breaks = NULL) +
facet_grid(sexism~key) +
theme_black() +
theme(panel.grid = element_blank()) Now we have model1_fitted, it’s easy to get the typical numeric summaries for the differences.

model1_fitted %>%
select(sexism, contains("-")) %>%
gather(key, value, -sexism) %>%
group_by(key, sexism) %>%
summarize(mean = mean(value),
ll = quantile(value, probs = .025),
ul = quantile(value, probs = .975)) %>%
mutate_if(is.double, round, digits = 3)
## # A tibble: 9 x 5
## # Groups:   key 
##   key                                     sexism   mean     ll    ul
##   <chr>                                    <dbl>  <dbl>  <dbl> <dbl>
## 1 Collective Protest - Individual Protest   4.25 -0.11  -0.713 0.487
## 2 Collective Protest - Individual Protest   5.12 -0.147 -0.539 0.254
## 3 Collective Protest - Individual Protest   5.90 -0.179 -0.72  0.361
## 4 Collective Protest - No Protest           4.25 -0.542 -1.11  0.036
## 5 Collective Protest - No Protest           5.12 -0.109 -0.569 0.338
## 6 Collective Protest - No Protest           5.90  0.277 -0.381 0.904
## 7 Individual Protest - No Protest           4.25 -0.432 -1.05  0.18
## 8 Individual Protest - No Protest           5.12  0.038 -0.373 0.465
## 9 Individual Protest - No Protest           5.90  0.456 -0.147 1.06

We don’t have $$p$$-values, but all the differences are small in magnitude and have wide 95% intervals straddling zero.

To get the difference scores Hayes presented on pages 498–500, one might:

post %>%
mutate(Difference in liking between being told she protested or not when W is 4.250 = b_liking_D1 + b_liking_D1:sexism*4.250,
Difference in liking between being told she protested or not when W is 5.120 = b_liking_D1 + b_liking_D1:sexism*5.120,
Difference in liking between being told she protested or not when W is 5.896 = b_liking_D1 + b_liking_D1:sexism*5.896,

Difference in liking between collective vs. individual protest when W is 4.250 = b_liking_D2 + b_liking_D2:sexism*4.250,
Difference in liking between collective vs. individual protest when W is 5.120 = b_liking_D2 + b_liking_D2:sexism*5.120,
Difference in liking between collective vs. individual protest when W is 5.896 = b_liking_D2 + b_liking_D2:sexism*5.896) %>%
select(contains("Difference in liking")) %>%
gather() %>%
group_by(key) %>%
summarize(mean = mean(value),
ll = quantile(value, probs = .025),
ul = quantile(value, probs = .975)) %>%
mutate_if(is.double, round, digits = 3)
## # A tibble: 6 x 4
##   key                                                                             mean     ll    ul
##   <chr>                                                                          <dbl>  <dbl> <dbl>
## 1 Difference in liking between being told she protested or not when W is 4.250  -0.487 -1.00  0.017
## 2 Difference in liking between being told she protested or not when W is 5.120  -0.036 -0.43  0.35
## 3 Difference in liking between being told she protested or not when W is 5.896   0.367 -0.207 0.923
## 4 Difference in liking between collective vs. individual protest when W is 4.2… -0.11  -0.713 0.487
## 5 Difference in liking between collective vs. individual protest when W is 5.1… -0.147 -0.539 0.254
## 6 Difference in liking between collective vs. individual protest when W is 5.8… -0.179 -0.72  0.361