Extensions of Bell and Banks (2018)

library(tidyverse)
library(haven)
library(lme4)
library(broom)
library(broom.mixed)
library(emmeans)
library(brms)
library(tidybayes)
library(modelsummary)
library(kableExtra)
library(patchwork)
library(here)

options(
  mc.cores = 4,
  brms.backend = "cmdstanr"
)

ssq <- read_stata(here("data", "Final.9.26.2016.dta"))

Bayesian versions of original models

Table 3

table_3_b <- brm(
  bf(f_Overall3P ~ wingos + whrnoposmgoldtotal + 
       (wingos * whrnoposmgoldtotal) + 
       cedaw_rat_1 + log_gdp + log_pop + polity2 +
       overallglobalizationindex + (1 | ccode),
     decomp = "QR"), 
  data = ssq,
  family = "gaussian",
  chains = 4)

tidy(table_3_b) %>% 
  kbl() %>% 
  kable_styling()

effect	component	group	term	estimate	std.error	conf.low	conf.high
fixed	cond	NA	(Intercept)	0.2863300	2.3472551	-4.4235070	4.6888592
fixed	cond	NA	wingos	0.0175739	0.0061941	0.0055428	0.0296754
fixed	cond	NA	whrnoposmgoldtotal	-0.0173717	0.0106123	-0.0379116	0.0033147
fixed	cond	NA	cedaw_rat_1	1.0270168	0.5979503	-0.1969347	2.2089430
fixed	cond	NA	log_gdp	-0.2149541	0.1703789	-0.5533549	0.1163722
fixed	cond	NA	log_pop	0.2010756	0.1175592	-0.0241090	0.4330191
fixed	cond	NA	polity2	0.0695767	0.0245524	0.0197258	0.1169868
fixed	cond	NA	overallglobalizationindex	0.0945391	0.0168001	0.0622672	0.1281644
fixed	cond	NA	wingos:whrnoposmgoldtotal	0.0002117	0.0001606	-0.0000983	0.0005308
ran_pars	cond	ccode	sd__(Intercept)	1.6011981	0.1256395	1.3736490	1.8682823
ran_pars	cond	Residual	sd__Observation	1.4891239	0.0390635	1.4152690	1.5680955

Marginal effects

ame1 <- table_3_b %>% 
  emtrends(~ whrnoposmgoldtotal,
           var = "wingos",
           at = list(whrnoposmgoldtotal = seq(0, 107, 1)),
           epred = TRUE) %>% 
  gather_emmeans_draws()

plot_ame1 <- ggplot(ame1, aes(x = whrnoposmgoldtotal, y = .value)) +
  stat_lineribbon() +
  geom_hline(yintercept = 0) +
  labs(x = "Women's INGO shaming", y = "Marginal effect of INGO presence",
       title = "Panel A: 3P index", fill = "Credible interval") +
  coord_cartesian(ylim = c(-0.05, 0.1)) +
  theme_bw() +
  theme(panel.grid.major.x = element_blank(),
        panel.grid.minor.x = element_blank(),
        panel.grid.minor.y = element_blank(),
        legend.position = "bottom")

ame2 <- table_3_b %>% 
  emtrends(~ wingos,
           var = "whrnoposmgoldtotal",
           at = list(wingos = seq(0, 107, 1)),
           epred = TRUE) %>% 
  gather_emmeans_draws()

plot_ame2 <- ggplot(ame2, aes(x = wingos, y = .value)) +
  stat_lineribbon() +
  geom_hline(yintercept = 0) +
  guides(fill = "none") +
  labs(x = "Women's INGO presence", y = "Marginal effect of shaming",
       title = "Panel B: 3P index") +
  coord_cartesian(ylim = c(-0.05, 0.05)) +
  theme_bw() +
  theme(panel.grid.major.x = element_blank(),
        panel.grid.minor.x = element_blank(),
        panel.grid.minor.y = element_blank())

plot_ame1 | plot_ame2

Binary treatment

Basic MSMs

Official Bayesian MSMs