H₂: Derogations and human rights

library(tidyverse)
library(targets)
library(broom)
library(modelsummary)
library(kableExtra)
library(brms)
library(marginaleffects)
library(tidybayes)
library(ggdist)
library(bayesplot)
library(patchwork)
library(colorspace)
library(scales)

# Bayes stuff
BAYES_SEED <- 1234
options(mc.cores = 4,
        brms.backend = "cmdstanr")

# Generated via random.org
set.seed(2385)

tar_config_set(store = here::here('_targets'),
               script = here::here('_targets.R'))

# Load targets 
# Final data
tar_load(daily_panel)
tar_load(weekly_panel)
tar_load(year_week_lookup)

# Models
tar_load(m_human_rights)
tar_load(human_rights_plot_data)
tar_load(models_tbl_human_rights)

# Plotting functions
invisible(list2env(tar_read(graphic_functions), .GlobalEnv))
invisible(list2env(tar_read(diagnostic_functions), .GlobalEnv))
invisible(list2env(tar_read(modelsummary_functions), .GlobalEnv))

# Tell bayesplot to use the Lakota palette for things like pp_check()
bayesplot::color_scheme_set(clrs[c(1, 2, 4, 5, 7, 8)])

weekly_panel <- weekly_panel %>% 
  mutate(across(c(new_cases, new_deaths, cumulative_cases, cumulative_deaths),
                list(z = ~as.numeric(scale(.)))))

Model and priors

Official formal model

\[ \begin{aligned} \text{COVID policy}_{i, t} \sim&\ \operatorname{Bernoulli}(\pi_{i, t}) \\[0.75em] \operatorname{logit}(\pi_{i, t}) =&\ \alpha_{\text{region} [i]} + \beta_1\ \text{Derogation in effect}_{i}\ + \\ &\ \beta_2\ \text{New cases}_{i}\ + \beta_3\ \text{Cumulative cases}_{i}\ + \\ &\ \beta_4\ \text{New deaths}_{i}\ + \beta_5\ \text{Cumulative deaths}_{i}\ + \\ &\ \beta_6\ \text{Past ICCPR derogation}_{i}\ + \beta_7\ \text{Past ICCPR action}_{i}\ + \\ &\ \beta_8\ \text{Rule of law index}_{i}\ + \beta_9\ \text{Civil liberties index}_{i}\ + \\ &\ \beta_{10}\ \text{Core civil society index}_{i}\ + \beta_{11}\ \text{Week number}_{i} \\ \alpha_j \sim&\ \mathcal{N}(\bar{\alpha}, \sigma_\alpha) \text{, for } j \text{ in } 1 \dots J \\[0.75em] \bar{\alpha} \sim&\ \operatorname{Student t}(\nu = 1, \mu = 0, \sigma = 3) \\ \beta_{1 \dots 11} \sim&\ \operatorname{Student t}(\nu = 1, \mu = 0, \sigma = 3) \\ \sigma_a \sim&\ \operatorname{Cauchy}(x = 0, \gamma = 1) \text{, lower bound} = 0 \end{aligned} \]

Priors

p1 <- ggplot() +
  stat_function(geom = "area", 
                fun = ~extraDistr::dlst(., df = 1, mu = 0, sigma = 3), 
                fill = clrs[2]) +
  xlim(c(-20, 20)) +
  annotate(geom = "label", x = 0, y = 0.02, label = "Student t(1, 0, 3)") +
  labs(x = "α and βs") +
  theme_pandem(prior = TRUE)

p2 <- ggplot() +
  stat_function(geom = "area", 
                fun = ~dcauchy(., 0, 1), 
                fill = clrs[4]) +
  xlim(c(0, 10)) +
  annotate(geom = "label", x = 5, y = 0.063, label = "Cauchy(0, 1)") +
  labs(x = "σ") +
  theme_pandem(prior = TRUE)

p1 | p2

Prior simulation

ologit_priors <- c(prior(student_t(1, 0, 3), class = Intercept),
                   prior(student_t(1, 0, 3), class = b),
                   prior(cauchy(0, 1), class = sd, lb = 0))

model_human_rights_prior_only <- brm(
  bf(pandem_nolimit ~ derogation_ineffect + 
       new_cases_z + cumulative_cases_z + 
       new_deaths_z + cumulative_deaths_z + 
       prior_iccpr_derogations + prior_iccpr_other_action +
       v2x_rule + v2x_civlib + v2xcs_ccsi +
       year_week_num + (1 | country_name)), 
  data = weekly_panel,
  family = "cumulative",
  prior = ologit_priors,
  sample_prior = "only",
  chains = 4, seed = BAYES_SEED, iter = 5000, refresh = 0,
  threads = threading(2)  # Two CPUs per chain to speed things up
)
## Start sampling
## Running MCMC with 4 parallel chains, with 2 thread(s) per chain...
## 
## Chain 1 finished in 1.6 seconds.
## Chain 2 finished in 1.9 seconds.
## Chain 4 finished in 3.5 seconds.
## Chain 3 finished in 7.8 seconds.
## 
## All 4 chains finished successfully.
## Mean chain execution time: 3.7 seconds.
## Total execution time: 7.9 seconds.
model_human_rights_prior_only
## Warning: Parts of the model have not converged (some Rhats are > 1.05). Be careful when analysing the results! We
## recommend running more iterations and/or setting stronger priors.
##  Family: cumulative 
##   Links: mu = logit; disc = identity 
## Formula: pandem_nolimit ~ derogation_ineffect + new_cases_z + cumulative_cases_z + new_deaths_z + cumulative_deaths_z + prior_iccpr_derogations + prior_iccpr_other_action + v2x_rule + v2x_civlib + v2xcs_ccsi + year_week_num + (1 | country_name) 
##    Data: weekly_panel (Number of observations: 9591) 
##   Draws: 4 chains, each with iter = 5000; warmup = 2500; thin = 1;
##          total post-warmup draws = 10000
## 
## Group-Level Effects: 
## ~country_name (Number of levels: 139) 
##               Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## sd(Intercept)     9.21    401.69     0.04    26.70 1.00    11322     5514
## 
## Population-Level Effects: 
##                              Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## Intercept[1]                    -4.47    899.82 -1435.28  1427.91 1.01     4198     1335
## Intercept[2]                    35.66   1953.60 -1431.40  1545.34 1.01     4231     1371
## derogation_ineffect              3.25     52.64   -37.46    46.31 1.00     2399      876
## new_cases_z                      4.61     46.47   -33.23    61.62 1.01     1008      425
## cumulative_cases_z              82.04    307.76   -44.62  1321.16 1.11       29       13
## new_deaths_z                    -2.58     31.66   -41.65    29.09 1.01     1190      383
## cumulative_deaths_z              0.16     40.38   -41.39    35.78 1.00     4227     1191
## prior_iccpr_derogationsTRUE     -0.61     31.07   -33.36    33.18 1.00     3076     1070
## prior_iccpr_other_actionTRUE    -4.90     46.13   -69.96    35.31 1.01      553      161
## v2x_rule                         0.21     16.79   -25.91    30.00 1.00     4085     1711
## v2x_civlib                       0.73     14.18   -24.48    33.61 1.00     3640     1700
## v2xcs_ccsi                      -1.10     25.21   -34.50    28.82 1.00     2174      855
## year_week_num                    0.10     25.71   -40.80    41.16 1.01     4683     1286
## 
## Family Specific Parameters: 
##      Estimate Est.Error l-95% CI u-95% CI Rhat Bulk_ESS Tail_ESS
## disc     1.00      0.00     1.00     1.00   NA       NA       NA
## 
## Draws were sampled using sample(hmc). For each parameter, Bulk_ESS
## and Tail_ESS are effective sample size measures, and Rhat is the potential
## scale reduction factor on split chains (at convergence, Rhat = 1).

datagrid(model = model_human_rights_prior_only,
         year_week_num = 1:69,
         derogation_ineffect = 0:1) %>% 
  add_epred_draws(model_human_rights_prior_only, ndraws = 100) %>% 
  left_join(year_week_lookup, by = "year_week_num") %>% 
  mutate(derogation_ineffect = factor(derogation_ineffect, 
                                      levels = 0:1,
                                      labels = c("No", "Yes"),
                                      ordered = TRUE)) %>% 
  ggplot(aes(x = year_week_day, y = .epred, color = .category)) +
  geom_line(aes(group = paste(.draw, .category, derogation_ineffect)), alpha = 0.5, linewidth = 0.5) +
  scale_color_manual(values = c(clrs[c(2, 4, 6)])) +
  scale_x_date(date_breaks = "2 months", date_labels = "%b\n%Y") +
  scale_y_continuous(labels = label_percent()) +
  labs(x = NULL, y = "Predicted probability", 
       color = "Severity",
       title = "No Time Limit Measures (priors only)") +
  theme_pandem()

Posterior

notes <- c(
  "Estimates are median posterior log odds from ordered logistic and binary logistic regression models; 95% credible intervals (highest density posterior interval, or HDPI) in brackets.", 
  "Total \\(R^2\\) considers the variance of both population and group effects;", 
  "marginal \\(R^2\\) only takes population effects into account."
)

modelsummary(models_tbl_human_rights,
             estimate = "{estimate}",
             statistic = "[{conf.low}, {conf.high}]",
             coef_map = coef_map,
             gof_map = gof_map,
             output = "kableExtra",
             fmt =  list(estimate = 2, conf.low = 2, conf.high = 2)) %>% 
  kable_styling(htmltable_class = "table-sm light-border") %>% 
  footnote(general = notes, footnote_as_chunk = TRUE)

	Discriminatory Policy	Non-Derogable Rights	No Time Limit Measures	Abusive Enforcement
Derogation in effect	−0.20	−0.32	−1.84	0.22
	[−0.63, 0.19]	[−0.80, 0.18]	[−2.17, −1.48]	[0.04, 0.41]
New cases (standardized)	0.04	0.003	−0.13	0.03
	[−0.05, 0.13]	[−0.30, 0.25]	[−0.32, 0.04]	[−0.04, 0.11]
Cumulative cases (standardized)	−0.05	0.21	−0.45	0.23
	[−0.23, 0.14]	[−0.22, 0.59]	[−0.68, −0.24]	[0.11, 0.36]
New deaths (standardized)	−0.05	−0.13	0.06	0.08
	[−0.17, 0.07]	[−0.40, 0.09]	[−0.06, 0.18]	[−0.008, 0.16]
Cumulative deaths (standardized)	0.16	−0.40	0.28	−0.22
	[−0.06, 0.36]	[−0.78, −0.02]	[0.12, 0.44]	[−0.35, −0.09]
Past ICCPR derogation	0.94	0.24	0.08	0.38
	[0.78, 1.10]	[0.003, 0.47]	[−0.05, 0.21]	[0.26, 0.50]
Past ICCPR action	0.28	1.48	−0.52	0.07
	[0.13, 0.44]	[1.30, 1.68]	[−0.65, −0.40]	[−0.04, 0.19]
Rule of law	0.99	1.30	0.87	−0.48
	[0.53, 1.42]	[0.71, 1.89]	[0.59, 1.15]	[−0.76, −0.20]
Civil liberties	1.60	−5.47	−1.99	0.07
	[0.75, 2.42]	[−6.60, −4.33]	[−2.62, −1.40]	[−0.49, 0.69]
Core civil society index	−2.61	1.57	0.06	−0.47
	[−3.11, −2.04]	[0.83, 2.28]	[−0.37, 0.48]	[−0.89, −0.08]
Constant		−1.17
		[−2.07, −0.17]
Cut 1	1.43		0.43	−0.41
	[0.07, 2.64]		[−0.59, 1.39]	[−1.06, 0.17]
Cut 2	2.18		0.49	0.60
	[0.73, 3.30]		[−0.52, 1.45]	[−0.04, 1.19]
Cut 3	2.38			2.02
	[0.94, 3.50]			[1.41, 2.63]
Year-week	−0.01	−0.01	0.002	−0.03
	[−0.02, −0.01]	[−0.01, −0.005]	[−0.0002, 0.005]	[−0.03, −0.03]
Region random effects σ	1.38	0.96	1.13	0.63
	[0.71, 2.65]	[0.48, 1.80]	[0.57, 2.14]	[0.31, 1.18]
N	9591	9591	9591	9591
\(R^2\) (total)	0.15	0.12	0.07	0.13
\(R^2\) (marginal)	0.07	0.09	0.04	0.07
Note: Estimates are median posterior log odds from ordered logistic and binary logistic regression models; 95% credible intervals (highest density posterior interval, or HDPI) in brackets. Total \(R^2\) considers the variance of both population and group effects; marginal \(R^2\) only takes population effects into account.

Model evaluation

params_to_show <- c("b_derogation_ineffect", "b_year_week_num", 
                    "b_v2xcs_ccsi", "sd_who_region__Intercept")

Discriminatory Policy

plot_trace(m_human_rights$model[[1]], params_to_show)

plot_trank(m_human_rights$model[[1]], params_to_show)

plot_pp(m_human_rights$model[[1]])

Non-Derogable Rights

plot_trace(m_human_rights$model[[2]], params_to_show)

plot_trank(m_human_rights$model[[2]], params_to_show)

plot_pp(m_human_rights$model[[2]])

No Time Limit Measures

plot_trace(m_human_rights$model[[3]], params_to_show)

plot_trank(m_human_rights$model[[3]], params_to_show)

plot_pp(m_human_rights$model[[3]])

Abusive Enforcement

plot_trace(m_human_rights$model[[4]], params_to_show)

plot_trank(m_human_rights$model[[4]], params_to_show)

plot_pp(m_human_rights$model[[4]])

Interpretation and marginal effects

Predictions

human_rights_plots <- human_rights_plot_data %>% 
  # Predicted probabilities
  mutate(pred_plot = pmap(list(pred_data, nice, family), \(.data, .nice_name, .family) {
    if (.family == "cumulative") {
      .data %>% 
        ggplot(aes(
          x = year_week_day, y = .epred, ymin = .lower, ymax = .upper,
          color = .category, fill = after_scale(color),
          group = forcats::fct_rev(ordered(.width))
        )) +
        geom_lineribbon(alpha = 0.3) +
        scale_color_manual(values = c(clrs[c(2, 4, 6, 7)])) +
        scale_x_date(date_breaks = "2 months", date_labels = "%b\n%Y") +
        scale_y_continuous(labels = label_percent()) +
        labs(x = NULL, y = "Predicted probability", 
             color = "Severity", fill = "Severity",
             title = .nice_name) +
        facet_wrap(vars(derogation_ineffect)) +
        theme_pandem()
    } else {
      .data %>% 
        ggplot(aes(
          x = year_week_day, y = .epred, ymin = .lower, ymax = .upper,
          color = derogation_ineffect, fill = after_scale(color),
          group = forcats::fct_rev(ordered(.width))
        )) +
        geom_lineribbon(alpha = 0.3) +
        scale_color_manual(values = c(clrs[1], clrs[8])) +
        scale_x_date(date_breaks = "2 months", date_labels = "%b\n%Y") +
        scale_y_continuous(labels = label_percent()) +
        labs(x = NULL, y = "Predicted probability", 
             color = "Derogation in effect", fill = "Derogation in effect",
             title = .nice_name) +
        theme_pandem()
    }
  })) %>% 
  # Marginal effects / contrasts
  mutate(mfx_plot = pmap(list(mfx_data, nice, family), \(.data, .nice_name, .family) {
    if (.family == "cumulative") {
      .data %>% 
        mutate(across(c(draw, .lower, .upper), ~ . * 100)) %>%
        ggplot(aes(
          x = year_week_day, y = draw, ymin = .lower, ymax = .upper,
          color = group, fill = after_scale(color),
          group = forcats::fct_rev(ordered(.width))
        )) +
        geom_lineribbon(alpha = 0.3) +
        geom_hline(yintercept = 0, linewidth = 0.5, color = "grey50") +
        scale_color_manual(values = darken(clrs[c(2, 4, 6, 7)], 0.5, space = "combined"), guide = "none") +
        scale_x_date(date_breaks = "2 months", date_labels = "%b\n%Y") +
        labs(x = NULL,
             y = "Percentage point difference\nin probability of outcome level",
             title = .nice_name) +
        facet_wrap(vars(group), ncol = 2) +
        theme_pandem()
    } else {
      .data %>% 
        mutate(across(c(draw, .lower, .upper), ~ . * 100)) %>%
        ggplot(aes(
          x = year_week_day, y = draw, ymin = .lower, ymax = .upper,
          group = forcats::fct_rev(ordered(.width))
        )) +
        geom_lineribbon(alpha = 0.3, fill = clrs[5], color = clrs[5]) +
        geom_hline(yintercept = 0, linewidth = 0.5, color = "grey50") +
        scale_x_date(date_breaks = "2 months", date_labels = "%b\n%Y") +
        labs(x = NULL,
             y = "Percentage point difference\nin probability of outcome",
             title = .nice_name) +
        theme_pandem()
    }
  }))

wrap_plots(human_rights_plots$pred_plot)

Marginal effects

wrap_plots(human_rights_plots$mfx_plot)