## 12.1 Revisiting the disaster framing study

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

library(tidyverse)

glimpse(disaster)
## Observations: 211
## Variables: 5
## $id <int> 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... ##$ frame   <int> 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,...
## $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,... ##$ 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,...
## $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,... Load brms. library(brms) model1 is the simple moderation model. model1 <- brm(data = disaster, family = gaussian, donate ~ 1 + frame + skeptic + frame:skeptic, chains = 4, cores = 4) Our model1 summary matches nicely with the text. print(model1, digits = 3) ## Family: gaussian ## Links: mu = identity; sigma = identity ## Formula: donate ~ 1 + frame + skeptic + frame:skeptic ## 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 Eff.Sample Rhat ## Intercept 5.029 0.228 4.584 5.479 2070 1.002 ## frame 0.681 0.335 0.014 1.333 1969 1.001 ## skeptic -0.139 0.058 -0.255 -0.025 1949 1.003 ## frame:skeptic -0.171 0.085 -0.341 0.003 1816 1.002 ## ## Family Specific Parameters: ## Estimate Est.Error l-95% CI u-95% CI Eff.Sample Rhat ## sigma 1.240 0.062 1.124 1.366 2855 1.001 ## ## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample ## is a crude measure of effective sample size, and Rhat is the potential ## scale reduction factor on split chains (at convergence, Rhat = 1). For the figures in this chapter, we’ll take theme cues from Matthew Kay’s tidybayes package. Otherwise, our Figure 12.2 is business as usual at this point. theme_set(theme_light()) nd <- tibble(frame = rep(0:1, each = 30), skeptic = rep(seq(from = 0, to = 7, length.out = 30), times = 2)) fitted(model1, newdata = nd) %>% as_tibble() %>% bind_cols(nd) %>% mutate(frame = ifelse(frame == 0, str_c("Natural causes (X = ", frame, ")"), str_c("Climate change (X = ", frame, ")"))) %>% ggplot(aes(x = skeptic, y = Estimate)) + geom_ribbon(aes(ymin = Q2.5, ymax = Q97.5, fill = frame), alpha = 1/3) + geom_line(aes(color = frame)) + scale_fill_brewer(type = "qual") + scale_color_brewer(type = "qual") + coord_cartesian(xlim = 1:6, ylim = c(3.5, 5.5)) + labs(x = expression(paste("Climate Change Skepticism (", italic(W), ")")), y = "Willingness to Donate to Victims") + theme(legend.position = "top", legend.direction = "horizontal", legend.title = element_blank()) In Hayes’s Figure 12.2, he emphasized the differences at the three levels of skeptic. If you want the full difference score distributions in a pick-a-point-approach sort of way, you might plot the densities with tidybayes::geom_halfeyeh(), which places coefficient plots at the base of the densities. In this case, we show the posterior medians with the dots, the 50% intervals with the thick horizontal lines, and the 95% intervals with the thinner horizontal lines. library(tidybayes) nd <- tibble(frame = rep(0:1, times = 3), skeptic = rep(quantile(disaster$skeptic, probs = c(.16, .5, .86)),
times = 2))

fitted(model1, summary = F,
newdata = nd) %>%
as_tibble() %>%
gather() %>%
mutate(frame = rep(rep(0:1, times = 3),
each = 4000),
skeptic = rep(rep(quantile(disaster$skeptic, probs = c(.16, .5, .86)), times = 2), each = 4000), iter = rep(1:4000, times = 6)) %>% select(-key) %>% spread(key = frame, value = value) %>% mutate(difference = 1 - 0) %>% ggplot(aes(x = difference, y = skeptic, fill = skeptic %>% as.character())) + geom_halfeyeh(point_interval = median_qi, .prob = c(0.95, 0.5)) + scale_fill_brewer() + scale_y_continuous(breaks = quantile(disaster$skeptic, probs = c(.16, .5, .86)),
labels = quantile(disaster\$skeptic, probs = c(.16, .5, .86)) %>% round(2)) +
theme(legend.position = "none",
panel.grid.minor.y = element_blank())

Here’s our simple mediation model, model2, using the multivariate syntax right in the brm() function.

model2 <-
brm(data = disaster, family = gaussian,
bf(justify ~ 1 + frame) +
bf(donate ~ 1 + frame + justify) +
set_rescor(FALSE),
chains = 4, cores = 4)
print(model2, digits = 3)
##  Family: MV(gaussian, gaussian)
##   Links: mu = identity; sigma = identity
##          mu = identity; sigma = identity
## Formula: justify ~ 1 + frame
##          donate ~ 1 + frame + justify
##    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 Eff.Sample  Rhat
## justify_Intercept    2.801     0.091    2.619    2.981       4000 1.000
## donate_Intercept     7.235     0.233    6.776    7.688       4000 0.999
## justify_frame        0.134     0.127   -0.116    0.379       4000 1.000
## donate_frame         0.211     0.135   -0.060    0.480       4000 1.000
## donate_justify      -0.954     0.076   -1.102   -0.805       4000 0.999
##
## Family Specific Parameters:
##               Estimate Est.Error l-95% CI u-95% CI Eff.Sample  Rhat
## sigma_justify    0.935     0.045    0.850    1.026       4000 1.000
## sigma_donate     0.987     0.049    0.898    1.090       4000 0.999
##
## Samples were drawn using sampling(NUTS). For each parameter, Eff.Sample
## is a crude measure of effective sample size, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

Consider the Bayesian $$R^2$$ summaries.

bayes_R2(model2) %>% round(digits = 3)
##            Estimate Est.Error  Q2.5 Q97.5
## R2_justify    0.010     0.011 0.000 0.040
## R2_donate     0.449     0.039 0.367 0.519

If you want the indirect effect with its intervals, you use posterior_samples() and data wrangle, as usual.

posterior_samples(model2) %>%
mutate(ab = b_justify_frame*b_donate_justify) %>%
summarize(mean = mean(ab),
ll = quantile(ab, probs = .025),
ul = quantile(ab, probs = .975)) %>%
mutate_if(is.double, round, digits = 3)
##     mean     ll   ul
## 1 -0.128 -0.364 0.11

We might also streamline our code a touch using tidybayes::mean_qi() in place of tidyverse::summarize().

posterior_samples(model2) %>%
mutate(ab = b_justify_frame*b_donate_justify) %>%
mean_qi(ab, .prob = .95) %>%
mutate_if(is.double, round, digits = 3)
## # A tibble: 1 x 4
##       ab conf.low conf.high .prob
##    <dbl>    <dbl>     <dbl> <dbl>
## 1 -0.128   -0.364      0.11  0.95

Note that the last column explicates what interval level we used.