## 7.6 Artificial categorization and subgroups

There are multiple ways to dichotomize the data by skeptic. A quick simple way is to use ifelse() to make a skeptic_hi dummy.

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

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

model6_low <-
brm(data = disaster %>% filter(skeptic_hi == 0),
family = gaussian,
justify ~ 1 + frame,
chains = 4, cores = 4)

model6_high <-
update(model6_low, newdata = disaster %>% filter(skeptic_hi == 1),
chains = 4, cores = 4)

The coefficient summaries:

fixef(model6_low) %>% round(digits = 3)
##           Estimate Est.Error   Q2.5 Q97.5
## Intercept    2.626     0.099  2.431 2.824
## frame       -0.108     0.145 -0.390 0.173
fixef(model6_high) %>% round(digits = 3)
##           Estimate Est.Error  Q2.5 Q97.5
## Intercept    3.070     0.144 2.786 3.348
## frame        0.477     0.209 0.071 0.891

You can use fitted() to get the posterior means and other summaries for the two frame groups, by model.

fitted(model6_low,
newdata = tibble(frame = 0:1)) %>%
round(digits = 3)
##      Estimate Est.Error  Q2.5 Q97.5
## [1,]    2.626     0.099 2.431 2.824
## [2,]    2.518     0.105 2.309 2.724
fitted(model6_high,
newdata = tibble(frame = 0:1)) %>%
round(digits = 3)
##      Estimate Est.Error  Q2.5 Q97.5
## [1,]    3.070     0.144 2.786 3.348
## [2,]    3.548     0.149 3.256 3.836

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 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 model6_low and model6_high. Here we’ll use bind_rows() to combine their posterior samples and then plot.

posterior_samples(model6_low) %>%
select(b_frame) %>%
bind_rows(
posterior_samples(model6_high) %>%
select(b_frame)
) %>%
mutate(model = rep(c("model6_low", "model6_high"), 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 levels,\ncolor 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.