Comparing OLS w/ interactions to RE w/ and w/o interactions with oversight data

comparing.R

# load packages
library(tidyverse)
library(marginaleffects)
library(rstanarm); options(mc.cores = parallel::detectCores())

# load data
data <- read_rds("https://www.dropbox.com/s/9grn8kkb5yzwagx/data.rds?raw=1") %>%
  filter(passed_mvc1) %>%
  glimpse()

# simplest model
fit1 <- lm(cong_overall ~ failure*amplify, data = data)

# adding partisan strength and partisanship
fit2 <- lm(cong_overall ~ failure*amplify*pid_strength, data = data)

# adding partisan strength and partisanship
fit3 <- lm(cong_overall ~ failure*amplify*pid_strength*pid2, data = data)

# AIC likes the increasingly complicated models better
AIC(fit1, fit2, fit3)

c <- comparisons(fit3, 
            variables = list(amplify = "reference"), 
            newdata = datagrid(pid_strength = unique, 
                               pid2 = unique, 
                               failure = unique)) %>%
  glimpse()

gg1 <- ggplot(c, aes(x = pid_strength, y = estimate, 
              ymin = conf.low, ymax = conf.high, 
              color = failure)) + 
  geom_errorbar(width = 0, position = position_dodge(width = 0.2)) + 
  geom_point(position = position_dodge(width = 0.2)) + 
  facet_wrap(vars(pid2)) + 
  labs(title = "OLS With Interactions", 
       x = "PID Strength", 
       y = "Effect of Congress Amplifying\nWhistleblower on Their Approval", 
       color = "Outcome of Operation")


# adding partisan strength and partisanship

hier_fit <- stan_lmer(cong_overall ~ failure*amplify + 
                        (failure*amplify | pid_strength:pid2) + 
                        (failure*amplify | pid_strength) + 
                        (failure*amplify | pid2), data = data)


c <- comparisons(hier_fit, 
                 variables = list(amplify = "reference"), 
                 newdata = datagrid(pid_strength = unique, 
                                    pid2 = unique, 
                                    failure = unique)) %>%
  glimpse()

gg2 <- ggplot(c, aes(x = pid_strength, y = estimate, 
              ymin = conf.low, ymax = conf.high, 
              color = failure)) + 
  geom_errorbar(width = 0, position = position_dodge(width = 0.2)) + 
  geom_point(position = position_dodge(width = 0.2)) + 
  facet_wrap(vars(pid2)) + 
  labs(title = "Random Effects Without Group-Level Predictors", 
       x = "PID Strength", 
       y = "Effect of Congress Amplifying\nWhistleblower on Their Approval", 
       color = "Outcome of Operation")

# adding partisan strength and partisanship

hier_fit2 <- stan_lmer(cong_overall ~ failure*amplify*pid_strength*pid2 + 
                        (failure*amplify | pid_strength:pid2) + 
                        (failure*amplify | pid_strength) + 
                        (failure*amplify | pid2), data = data)


c <- comparisons(hier_fit2, 
                 variables = list(amplify = "reference"), 
                 newdata = datagrid(pid_strength = unique, 
                                    pid2 = unique, 
                                    failure = unique)) %>%
  glimpse()

gg3 <- ggplot(c, aes(x = pid_strength, y = estimate, 
              ymin = conf.low, ymax = conf.high, 
              color = failure)) + 
  geom_errorbar(width = 0, position = position_dodge(width = 0.2)) + 
  geom_point(position = position_dodge(width = 0.2)) + 
  facet_wrap(vars(pid2)) + 
  labs(title = "RE With Group-Level Predictors", 
       x = "PID Strength", 
       y = "Effect of Congress Amplifying\nWhistleblower on Their Approval", 
       color = "Outcome of Operation")

library(patchwork)

gg1 / gg2 / gg3