Model details

library(tidyverse)
library(targets)
library(brms)
library(glue)
library(kableExtra)
library(lubridate)
library(here)

withr::with_dir(here::here(), {
  # Treatment models
  tar_load(c(m_oda_treatment_total, m_oda_treatment_advocacy, 
             m_oda_treatment_entry, m_oda_treatment_funding, 
             m_oda_treatment_ccsi, m_oda_treatment_repress))
  tar_load(c(m_purpose_treatment_total, m_purpose_treatment_advocacy, 
             m_purpose_treatment_entry, m_purpose_treatment_funding, 
             m_purpose_treatment_ccsi, m_purpose_treatment_repress))
  tar_load(c(m_recip_treatment_total_dom, m_recip_treatment_advocacy_dom, 
             m_recip_treatment_entry_dom, m_recip_treatment_funding_dom, 
             m_recip_treatment_ccsi_dom, m_recip_treatment_repress_dom))
  tar_load(c(m_recip_treatment_total_foreign, m_recip_treatment_advocacy_foreign, 
             m_recip_treatment_entry_foreign, m_recip_treatment_funding_foreign, 
             m_recip_treatment_ccsi_foreign, m_recip_treatment_repress_foreign))
  
  # Outcome models
  tar_load(c(m_oda_outcome_total, m_oda_outcome_advocacy, 
             m_oda_outcome_entry, m_oda_outcome_funding,
             m_oda_outcome_ccsi, m_oda_outcome_repress))
  tar_load(c(m_purpose_outcome_total, m_purpose_outcome_advocacy, 
             m_purpose_outcome_entry, m_purpose_outcome_funding,
             m_purpose_outcome_ccsi, m_purpose_outcome_repress))
  tar_load(c(m_recip_outcome_total_dom, m_recip_outcome_advocacy_dom, 
             m_recip_outcome_entry_dom, m_recip_outcome_funding_dom,
             m_recip_outcome_ccsi_dom, m_recip_outcome_repress_dom))
  tar_load(c(m_recip_outcome_total_foreign, m_recip_outcome_advocacy_foreign, 
             m_recip_outcome_entry_foreign, m_recip_outcome_funding_foreign,
             m_recip_outcome_ccsi_foreign, m_recip_outcome_repress_foreign))
})

model_df <- tibble(
  model_name = c("m_oda_treatment_total", "m_oda_treatment_advocacy",
                 "m_oda_treatment_entry", "m_oda_treatment_funding",
                 "m_oda_treatment_ccsi", "m_oda_treatment_repress",
                 "m_oda_outcome_total", "m_oda_outcome_advocacy",
                 "m_oda_outcome_entry", "m_oda_outcome_funding",
                 "m_oda_outcome_ccsi", "m_oda_outcome_repress",
                 "m_purpose_treatment_total", "m_purpose_treatment_advocacy",
                 "m_purpose_treatment_entry", "m_purpose_treatment_funding",
                 # "m_purpose_treatment_ccsi", "m_purpose_treatment_repress",
                 "m_purpose_outcome_total", "m_purpose_outcome_advocacy",
                 "m_purpose_outcome_entry", "m_purpose_outcome_funding",
                 # "m_purpose_outcome_ccsi", "m_purpose_outcome_repress"
                 "m_recip_treatment_total_dom", "m_recip_treatment_advocacy_dom", 
                 "m_recip_treatment_entry_dom", "m_recip_treatment_funding_dom", 
                 # "m_recip_treatment_ccsi_dom", "m_recip_treatment_repress_dom",
                 "m_recip_outcome_total_dom", "m_recip_outcome_advocacy_dom", 
                 "m_recip_outcome_entry_dom", "m_recip_outcome_funding_dom",
                 # "m_recip_outcome_ccsi_dom", "m_recip_outcome_repress_dom",
                 "m_recip_treatment_total_foreign", "m_recip_treatment_advocacy_foreign", 
                 "m_recip_treatment_entry_foreign", "m_recip_treatment_funding_foreign", 
                 # "m_recip_treatment_ccsi_foreign", "m_recip_treatment_repress_foreign",
                 "m_recip_outcome_total_foreign", "m_recip_outcome_advocacy_foreign", 
                 "m_recip_outcome_entry_foreign", "m_recip_outcome_funding_foreign"#,
                 # "m_recip_outcome_ccsi_foreign", "m_recip_outcome_repress_foreign"
  )) %>% 
  mutate(actual_model = map(model_name, ~eval(rlang::sym(.)))) %>% 
  # Put standalone unlisted brms models in a list so that they can be unnested
  mutate(actual_model = ifelse(map_chr(actual_model, ~class(.)) == "brmsfit",
                               map(actual_model, ~list(.)), actual_model)) %>% 
  mutate(list_names = map(actual_model, ~names(.)))

model_details_df <- model_df %>% 
  unnest(c(actual_model, list_names)) %>% 
  mutate(stage = ifelse(str_detect(model_name, "treatment"),
                        "Treatment", "Outcome")) %>% 
  mutate(category = case_when(
    str_detect(list_names, "_num") ~ "Numerator",
    str_detect(list_names, "_denom") ~ "Denominator",
    str_detect(list_names, "\\d+") ~ paste0("IPTW truncated at ", str_extract(list_names, "\\d+")),
    TRUE ~ NA_character_
  )) %>% 
  mutate(outcome = case_when(
    str_detect(model_name, "_oda_") ~ "Total aid",
    str_detect(model_name, "_purpose_") ~ "Aid contentiousness",
    str_detect(model_name, "recip_.*_dom") ~ "Aid recipients (domestic)",
    str_detect(model_name, "recip_.*_foreign") ~ "Aid recipients (foreign)"
  )) %>% 
  mutate(treatment = case_when(
    str_detect(model_name, "_total") ~ "Total barriers",
    str_detect(model_name, "_advocacy") ~ "Barriers to advocacy",
    str_detect(model_name, "_entry") ~ "Barriers to entry",
    str_detect(model_name, "_funding") ~ "Barriers to funding",
    str_detect(model_name, "_ccsi") ~ "Civil society index",
    str_detect(model_name, "_repress") ~ "Civil society repression"
  )) %>% 
  mutate(across(c(stage, category, outcome, treatment), ~fct_inorder(., ordered = TRUE)))

model_details_summarized <- model_details_df %>% 
  mutate(duration = map(actual_model, ~rstan::get_elapsed_time(.$fit)),
         duration = map(duration, ~rownames_to_column(as_tibble(.)))) %>% 
  select(-actual_model) %>% 
  unnest(duration) %>% 
  mutate(model_name = glue("<code>{model_name}</code>")) %>% 
  group_by(Model = model_name, Treatment = treatment, Outcome = outcome,
           Stage = stage, Details = category) %>% 
  summarize(`Total time (i.e. longest chain)` = as.duration(max(warmup + sample))) %>%
  ungroup() %>% 
  arrange(Outcome, Treatment, Stage, Details)

We ran these models on a 2016 MacBook Pro with 4 cores and 16 GB of RAM, using Stan through brms through cmdstanr. In total, it took 12 hours, 34 minutes, and 8 seconds to run everything.

 

Model run times

Models for H₁ (total ODA)

model_details_h1 <- model_details_summarized %>% 
  filter(Outcome == "Total aid")

total_row_h1 <- tibble(
  Treatment = "Total", 
  Model = "All models",
  `Total time (i.e. longest chain)` = 
    as.duration(sum(model_details_h1$`Total time (i.e. longest chain)`))
)

model_time_h1 <- model_details_h1 %>% 
  bind_rows(total_row_h1) 

model_time_h1 %>% 
  select(-Treatment, -Outcome) %>% 
  rename(` ` = Model) %>% 
  kbl(escape = FALSE) %>% 
  pack_rows(index = table(fct_inorder(model_time_h1$Treatment))) %>% 
  kable_styling()

	Stage	Details	Total time (i.e. longest chain)
Total barriers
m_oda_treatment_total	Treatment	Numerator	189.01s (~3.15 minutes)
m_oda_treatment_total	Treatment	Denominator	151.496s (~2.52 minutes)
m_oda_outcome_total	Outcome		396.215s (~6.6 minutes)
Barriers to advocacy
m_oda_treatment_advocacy	Treatment	Numerator	119.113s (~1.99 minutes)
m_oda_treatment_advocacy	Treatment	Denominator	282.959s (~4.72 minutes)
m_oda_outcome_advocacy	Outcome		279.767s (~4.66 minutes)
Barriers to entry
m_oda_treatment_entry	Treatment	Numerator	109.664s (~1.83 minutes)
m_oda_treatment_entry	Treatment	Denominator	185.992s (~3.1 minutes)
m_oda_outcome_entry	Outcome		302.05s (~5.03 minutes)
Barriers to funding
m_oda_treatment_funding	Treatment	Numerator	141.48s (~2.36 minutes)
m_oda_treatment_funding	Treatment	Denominator	168.247s (~2.8 minutes)
m_oda_outcome_funding	Outcome		392.747s (~6.55 minutes)
Civil society index
m_oda_treatment_ccsi	Treatment	Numerator	192.624s (~3.21 minutes)
m_oda_treatment_ccsi	Treatment	Denominator	568.706s (~9.48 minutes)
m_oda_outcome_ccsi	Outcome	IPTW truncated at 100	249.96s (~4.17 minutes)
m_oda_outcome_ccsi	Outcome	IPTW truncated at 500	439.585s (~7.33 minutes)
m_oda_outcome_ccsi	Outcome	IPTW truncated at 1000	497.64s (~8.29 minutes)
m_oda_outcome_ccsi	Outcome	IPTW truncated at 5000	1004.137s (~16.74 minutes)
Civil society repression
m_oda_treatment_repress	Treatment	Numerator	293.43s (~4.89 minutes)
m_oda_treatment_repress	Treatment	Denominator	268.64s (~4.48 minutes)
m_oda_outcome_repress	Outcome	IPTW truncated at 100	240.046s (~4 minutes)
m_oda_outcome_repress	Outcome	IPTW truncated at 500	571.196s (~9.52 minutes)
m_oda_outcome_repress	Outcome	IPTW truncated at 1000	481.962s (~8.03 minutes)
m_oda_outcome_repress	Outcome	IPTW truncated at 5000	895.046s (~14.92 minutes)
Total
All models			8421.712s (~2.34 hours)

Models for H₂ (aid contentiousness)

model_details_h2 <- model_details_summarized %>% 
  filter(Outcome == "Aid contentiousness")

total_row_h2 <- tibble(
  Treatment = "Total", 
  Model = "All models",
  `Total time (i.e. longest chain)` = 
    as.duration(sum(model_details_h2$`Total time (i.e. longest chain)`))
)

model_time_h2 <- model_details_h2 %>% 
  bind_rows(total_row_h2) 

model_time_h2 %>% 
  select(-Treatment, -Outcome) %>% 
  rename(` ` = Model) %>% 
  kbl(escape = FALSE) %>% 
  pack_rows(index = table(fct_inorder(model_time_h2$Treatment))) %>% 
  kable_styling()

	Stage	Details	Total time (i.e. longest chain)
Total barriers
m_purpose_treatment_total	Treatment	Numerator	105.479s (~1.76 minutes)
m_purpose_treatment_total	Treatment	Denominator	235.396s (~3.92 minutes)
m_purpose_outcome_total	Outcome		509.323s (~8.49 minutes)
Barriers to advocacy
m_purpose_treatment_advocacy	Treatment	Numerator	94.403s (~1.57 minutes)
m_purpose_treatment_advocacy	Treatment	Denominator	387.446s (~6.46 minutes)
m_purpose_outcome_advocacy	Outcome		673.934s (~11.23 minutes)
Barriers to entry
m_purpose_treatment_entry	Treatment	Numerator	119.452s (~1.99 minutes)
m_purpose_treatment_entry	Treatment	Denominator	241.706s (~4.03 minutes)
m_purpose_outcome_entry	Outcome		683.226s (~11.39 minutes)
Barriers to funding
m_purpose_treatment_funding	Treatment	Numerator	76.102s (~1.27 minutes)
m_purpose_treatment_funding	Treatment	Denominator	247.326s (~4.12 minutes)
m_purpose_outcome_funding	Outcome		619.938s (~10.33 minutes)
Total
All models			3993.731s (~1.11 hours)

Models for H₃ (aid recipients)

Domestic NGOs

model_details_h3_dom <- model_details_summarized %>% 
  filter(Outcome == "Aid recipients (domestic)")

total_row_h3_dom <- tibble(
  Treatment = "Total", 
  Model = "All models",
  `Total time (i.e. longest chain)` = 
    as.duration(sum(model_details_h3_dom$`Total time (i.e. longest chain)`))
)

model_time_h3_dom <- model_details_h3_dom %>% 
  bind_rows(total_row_h3_dom) 

model_time_h3_dom %>% 
  select(-Treatment, -Outcome) %>% 
  rename(` ` = Model) %>% 
  kbl(escape = FALSE) %>% 
  pack_rows(index = table(fct_inorder(model_time_h3_dom$Treatment))) %>% 
  kable_styling()

	Stage	Details	Total time (i.e. longest chain)
Total barriers
m_recip_treatment_total_dom	Treatment	Numerator	117.943s (~1.97 minutes)
m_recip_treatment_total_dom	Treatment	Denominator	181.979s (~3.03 minutes)
m_recip_outcome_total_dom	Outcome		542.721s (~9.05 minutes)
Barriers to advocacy
m_recip_treatment_advocacy_dom	Treatment	Numerator	83.917s (~1.4 minutes)
m_recip_treatment_advocacy_dom	Treatment	Denominator	307.042s (~5.12 minutes)
m_recip_outcome_advocacy_dom	Outcome		451.525s (~7.53 minutes)
Barriers to entry
m_recip_treatment_entry_dom	Treatment	Numerator	99.287s (~1.65 minutes)
m_recip_treatment_entry_dom	Treatment	Denominator	253.426s (~4.22 minutes)
m_recip_outcome_entry_dom	Outcome		448.13s (~7.47 minutes)
Barriers to funding
m_recip_treatment_funding_dom	Treatment	Numerator	72.501s (~1.21 minutes)
m_recip_treatment_funding_dom	Treatment	Denominator	192.327s (~3.21 minutes)
m_recip_outcome_funding_dom	Outcome		1206.392s (~20.11 minutes)
Total
All models			3957.19s (~1.1 hours)

Foreign NGOs

model_details_h3_foreign <- model_details_summarized %>% 
  filter(Outcome == "Aid recipients (foreign)")

total_row_h3_foreign <- tibble(
  Treatment = "Total", 
  Model = "All models",
  `Total time (i.e. longest chain)` = 
    as.duration(sum(model_details_h3_foreign$`Total time (i.e. longest chain)`))
)

model_time_h3_foreign <- model_details_h3_foreign %>% 
  bind_rows(total_row_h3_foreign) 

model_time_h3_foreign %>% 
  select(-Treatment, -Outcome) %>% 
  rename(` ` = Model) %>% 
  kbl(escape = FALSE) %>% 
  pack_rows(index = table(fct_inorder(model_time_h3_foreign$Treatment))) %>% 
  kable_styling()

	Stage	Details	Total time (i.e. longest chain)
Total barriers
m_recip_treatment_total_foreign	Treatment	Numerator	87.403s (~1.46 minutes)
m_recip_treatment_total_foreign	Treatment	Denominator	129.949s (~2.17 minutes)
m_recip_outcome_total_foreign	Outcome		24963.38s (~6.93 hours)
Barriers to advocacy
m_recip_treatment_advocacy_foreign	Treatment	Numerator	89.783s (~1.5 minutes)
m_recip_treatment_advocacy_foreign	Treatment	Denominator	442.092s (~7.37 minutes)
m_recip_outcome_advocacy_foreign	Outcome		776.995s (~12.95 minutes)
Barriers to entry
m_recip_treatment_entry_foreign	Treatment	Numerator	100.816s (~1.68 minutes)
m_recip_treatment_entry_foreign	Treatment	Denominator	240.17s (~4 minutes)
m_recip_outcome_entry_foreign	Outcome		1007.792s (~16.8 minutes)
Barriers to funding
m_recip_treatment_funding_foreign	Treatment	Numerator	82.593s (~1.38 minutes)
m_recip_treatment_funding_foreign	Treatment	Denominator	211.803s (~3.53 minutes)
m_recip_outcome_funding_foreign	Outcome		742.711s (~12.38 minutes)
Total
All models			28875.487s (~8.02 hours)

Actual code

All the models are run with a targets pipeline with dataset-specific functions that live in these files:

• R/models_oda.R (for total aid (H₁))
• R/models_purpose.R (for aid contentiousness (H₂))
• R/models_recipients.R (for aid recipients (H₃))
• R/funs_models-iptw.R (for calculating inverse probability weights)

To keep the analysis relatively self contained here in the analysis notebook, and to make it so there’s no need to hunt through the GitHub repository to find the actual code, here’s that code:

R/models_oda.R

# Settings ----------------------------------------------------------------

oda_setup <- function() {
  options(worker_options)
  
  # Settings
  CHAINS <- 4
  ITER <- 2000
  WARMUP <- 1000
  BAYES_SEED <- 4045  # From random.org
  
  # Priors
  prior_num <- c(set_prior("normal(0, 10)", class = "Intercept"),
                 set_prior("normal(0, 2.5)", class = "b"),
                 set_prior("cauchy(0, 1)", class = "sd"))
  
  prior_denom <- c(set_prior("normal(0, 10)", class = "Intercept"),
                   set_prior("normal(0, 2.5)", class = "b"),
                   set_prior("normal(0, 2.5)", class = "sd"))
  
  prior_out <- c(set_prior("normal(0, 20)", class = "Intercept"),
                 set_prior("normal(0, 3)", class = "b"),
                 set_prior("cauchy(0, 1)", class = "sd"))
  
  return(list(chains = CHAINS, iter = ITER, warmup = WARMUP, seed = BAYES_SEED,
              prior_num = prior_num, prior_denom = prior_denom, prior_out = prior_out))
}


# Treatment models --------------------------------------------------------

f_oda_treatment_total <- function(dat) {
  oda_settings <- oda_setup()
  
  dat <- dat %>% filter(laws)
  
  model_num <- brm(
    bf(barriers_total ~ barriers_total_lag1 + barriers_total_lag2_cumsum + (1 | gwcode)),
    data = dat,
    family = gaussian(),
    prior = oda_settings$prior_num,
    control = list(adapt_delta = 0.99),
    chains = oda_settings$chains, iter = oda_settings$iter,
    warmup = oda_settings$warmup, seed = oda_settings$seed
  )
  
  model_denom <- brm(
    bf(barriers_total ~ barriers_total_lag1 + barriers_total_lag2_cumsum + total_oda_log_lag1 +
         # Human rights and politics
         v2x_polyarchy + v2x_corr + v2x_rule + v2x_civlib + v2x_clphy + v2x_clpriv +
         # Economics and development
         gdpcap_log + un_trade_pct_gdp + v2peedueq + v2pehealth + e_peinfmor +
         # Conflict and disasters
         internal_conflict_past_5 + natural_dis_count +
         (1 | gwcode)),
    data = dat,
    family = gaussian(),
    prior = oda_settings$prior_num,
    control = list(adapt_delta = 0.9),
    chains = oda_settings$chains, iter = oda_settings$iter,
    warmup = oda_settings$warmup, seed = oda_settings$seed
  )
  
  return(lst(model_num, model_denom))
}

f_oda_treatment_advocacy <- function(dat) {
  oda_settings <- oda_setup()
  
  dat <- dat %>% filter(laws)
  
  model_num <- brm(
    bf(advocacy ~ advocacy_lag1 + advocacy_lag2_cumsum + (1 | gwcode)),
    data = dat,
    family = gaussian(),
    prior = oda_settings$prior_num,
    control = list(adapt_delta = 0.99),
    chains = oda_settings$chains, iter = oda_settings$iter,
    warmup = oda_settings$warmup, seed = oda_settings$seed
  )
  
  model_denom <- brm(
    bf(advocacy ~ advocacy_lag1 + advocacy_lag2_cumsum + total_oda_log_lag1 +
         v2x_polyarchy + v2x_corr + v2x_rule + v2x_civlib + v2x_clphy + v2x_clpriv +
         gdpcap_log + un_trade_pct_gdp + v2peedueq + v2pehealth + e_peinfmor +
         internal_conflict_past_5 + natural_dis_count +
         (1 | gwcode)),
    data = dat,
    family = gaussian(),
    prior = oda_settings$prior_num,
    control = list(adapt_delta = 0.9),
    chains = oda_settings$chains, iter = oda_settings$iter,
    warmup = oda_settings$warmup, seed = oda_settings$seed
  )
  
  return(lst(model_num, model_denom))
}

f_oda_treatment_entry <- function(dat) {
  oda_settings <- oda_setup()
  
  dat <- dat %>% filter(laws)
  
  model_num <- brm(
    bf(entry ~ entry_lag1 + entry_lag2_cumsum + (1 | gwcode)),
    data = dat,
    family = gaussian(),
    prior = oda_settings$prior_num,
    control = list(adapt_delta = 0.99),
    chains = oda_settings$chains, iter = oda_settings$iter,
    warmup = oda_settings$warmup, seed = oda_settings$seed
  )
  
  model_denom <- brm(
    bf(entry ~ entry_lag1 + entry_lag2_cumsum + total_oda_log_lag1 +
         v2x_polyarchy + v2x_corr + v2x_rule + v2x_civlib + v2x_clphy + v2x_clpriv +
         gdpcap_log + un_trade_pct_gdp + v2peedueq + v2pehealth + e_peinfmor +
         internal_conflict_past_5 + natural_dis_count +
         (1 | gwcode)),
    data = dat,
    family = gaussian(),
    prior = oda_settings$prior_num,
    control = list(adapt_delta = 0.9),
    chains = oda_settings$chains, iter = oda_settings$iter,
    warmup = oda_settings$warmup, seed = oda_settings$seed
  )
  
  return(lst(model_num, model_denom))
}

f_oda_treatment_funding <- function(dat) {
  oda_settings <- oda_setup()
  
  dat <- dat %>% filter(laws)
  
  model_num <- brm(
    bf(funding ~ funding_lag1 + funding_lag2_cumsum + (1 | gwcode)),
    data = dat,
    family = gaussian(),
    prior = oda_settings$prior_num,
    control = list(adapt_delta = 0.99),
    chains = oda_settings$chains, iter = oda_settings$iter,
    warmup = oda_settings$warmup, seed = oda_settings$seed
  )
  
  model_denom <- brm(
    bf(funding ~ funding_lag1 + funding_lag2_cumsum + total_oda_log_lag1 +
         v2x_polyarchy + v2x_corr + v2x_rule + v2x_civlib + v2x_clphy + v2x_clpriv +
         gdpcap_log + un_trade_pct_gdp + v2peedueq + v2pehealth + e_peinfmor +
         internal_conflict_past_5 + natural_dis_count +
         (1 | gwcode)),
    data = dat,
    family = gaussian(),
    prior = oda_settings$prior_num,
    control = list(adapt_delta = 0.9),
    chains = oda_settings$chains, iter = oda_settings$iter,
    warmup = oda_settings$warmup, seed = oda_settings$seed
  )
  
  return(lst(model_num, model_denom))
}

f_oda_treatment_ccsi <- function(dat) {
  oda_settings <- oda_setup()
  
  model_num <- brm(
    bf(v2xcs_ccsi ~ v2xcs_ccsi_lag1 + v2xcs_ccsi_lag2_cumsum + (1 | gwcode)),
    data = dat,
    family = gaussian(),
    prior = oda_settings$prior_num,
    control = list(adapt_delta = 0.99),
    chains = oda_settings$chains, iter = oda_settings$iter,
    warmup = oda_settings$warmup, seed = oda_settings$seed
  )
  
  model_denom <- brm(
    bf(v2xcs_ccsi ~ v2xcs_ccsi_lag1 + v2xcs_ccsi_lag2_cumsum + total_oda_log_lag1 +
         v2x_polyarchy + v2x_corr + v2x_rule + v2x_civlib + v2x_clphy + v2x_clpriv +
         gdpcap_log + un_trade_pct_gdp + v2peedueq + v2pehealth + e_peinfmor +
         internal_conflict_past_5 + natural_dis_count +
         (1 | gwcode)),
    data = dat,
    family = gaussian(),
    prior = oda_settings$prior_num,
    control = list(adapt_delta = 0.9),
    chains = oda_settings$chains, iter = oda_settings$iter,
    warmup = oda_settings$warmup, seed = oda_settings$seed
  )
  
  return(lst(model_num, model_denom))
}

f_oda_treatment_repress <- function(dat) {
  oda_settings <- oda_setup()
  
  model_num <- brm(
    bf(v2csreprss ~ v2csreprss_lag1 + v2csreprss_lag2_cumsum + (1 | gwcode)),
    data = dat,
    family = gaussian(),
    prior = oda_settings$prior_num,
    control = list(adapt_delta = 0.99),
    chains = oda_settings$chains, iter = oda_settings$iter,
    warmup = oda_settings$warmup, seed = oda_settings$seed
  )
  
  model_denom <- brm(
    bf(v2csreprss ~ v2csreprss_lag1 + v2csreprss_lag2_cumsum + total_oda_log_lag1 +
         v2x_polyarchy + v2x_corr + v2x_rule + v2x_civlib + v2x_clphy + v2x_clpriv +
         gdpcap_log + un_trade_pct_gdp + v2peedueq + v2pehealth + e_peinfmor +
         internal_conflict_past_5 + natural_dis_count +
         (1 | gwcode)),
    data = dat,
    family = gaussian(),
    prior = oda_settings$prior_num,
    control = list(adapt_delta = 0.9),
    chains = oda_settings$chains, iter = oda_settings$iter,
    warmup = oda_settings$warmup, seed = oda_settings$seed
  )
  
  return(lst(model_num, model_denom))
}


# Outcome models ----------------------------------------------------------

f_oda_outcome_total <- function(dat) {
  oda_settings <- oda_setup()
  
  dat <- dat %>% filter(laws)
  
  model <- brm(
    bf(total_oda_log_lead1 | weights(iptw) ~ 
         barriers_total + barriers_total_lag1_cumsum + 
         (1 | gwcode) + (1 | year)),
    data = dat,
    family = gaussian(),
    prior = oda_settings$prior_out,
    chains = oda_settings$chains, iter = oda_settings$iter * 2,
    warmup = oda_settings$warmup, seed = oda_settings$seed
  )
  
  return(model)
}

f_oda_outcome_advocacy <- function(dat) {
  oda_settings <- oda_setup()
  
  dat <- dat %>% filter(laws)
  
  model <- brm(
    bf(total_oda_log_lead1 | weights(iptw) ~ 
         advocacy + advocacy_lag1_cumsum + 
         (1 | gwcode) + (1 | year)),
    data = dat,
    family = gaussian(),
    prior = oda_settings$prior_out,
    chains = oda_settings$chains, iter = oda_settings$iter * 2,
    warmup = oda_settings$warmup, seed = oda_settings$seed
  )
  
  return(model)
}

f_oda_outcome_entry <- function(dat) {
  oda_settings <- oda_setup()
  
  dat <- dat %>% filter(laws)
  
  model <- brm(
    bf(total_oda_log_lead1 | weights(iptw) ~ 
         entry + entry_lag1_cumsum + 
         (1 | gwcode) + (1 | year)),
    data = dat,
    family = gaussian(),
    prior = oda_settings$prior_out,
    chains = oda_settings$chains, iter = oda_settings$iter * 2,
    warmup = oda_settings$warmup, seed = oda_settings$seed
  )
  
  return(model)
}

f_oda_outcome_funding <- function(dat) {
  oda_settings <- oda_setup()
  
  dat <- dat %>% filter(laws)
  
  model <- brm(
    bf(total_oda_log_lead1 | weights(iptw) ~ 
         funding + funding_lag1_cumsum + 
         (1 | gwcode) + (1 | year)),
    data = dat,
    family = gaussian(),
    prior = oda_settings$prior_out,
    chains = oda_settings$chains, iter = oda_settings$iter * 2,
    warmup = oda_settings$warmup, seed = oda_settings$seed
  )
  
  return(model)
}

f_oda_outcome_ccsi <- function(dat) {
  oda_settings <- oda_setup()
  
  dat <- dat %>%
    mutate(across(v2xcs_ccsi_lag1_cumsum, my_scale))
  
  dat_100 <- dat %>% mutate(iptw = ifelse(iptw > 100, 100, iptw))
  dat_500 <- dat %>% mutate(iptw = ifelse(iptw > 500, 500, iptw))
  dat_1000 <- dat %>% mutate(iptw = ifelse(iptw > 1000, 1000, iptw))
  dat_5000 <- dat %>% mutate(iptw = ifelse(iptw > 5000, 5000, iptw))
  
  model_100 <- brm(
    bf(total_oda_log_lead1 | weights(iptw) ~ 
         v2xcs_ccsi + v2xcs_ccsi_lag1_cumsum + 
         (1 | gwcode)),
    data = dat_100,
    family = gaussian(),
    prior = oda_settings$prior_out,
    chains = oda_settings$chains, iter = oda_settings$iter * 2,
    warmup = oda_settings$warmup, seed = oda_settings$seed
  )
  
  model_500 <- brm(
    bf(total_oda_log_lead1 | weights(iptw) ~ 
         v2xcs_ccsi + v2xcs_ccsi_lag1_cumsum + 
         (1 | gwcode)),
    data = dat_500,
    family = gaussian(),
    prior = oda_settings$prior_out,
    chains = oda_settings$chains, iter = oda_settings$iter * 2,
    warmup = oda_settings$warmup, seed = oda_settings$seed
  )
  
  model_1000 <- brm(
    bf(total_oda_log_lead1 | weights(iptw) ~ 
         v2xcs_ccsi + v2xcs_ccsi_lag1_cumsum + 
         (1 | gwcode)),
    data = dat_1000,
    family = gaussian(),
    prior = oda_settings$prior_out,
    chains = oda_settings$chains, iter = oda_settings$iter * 2,
    warmup = oda_settings$warmup, seed = oda_settings$seed
  )
  
  model_5000 <- brm(
    bf(total_oda_log_lead1 | weights(iptw) ~ 
         v2xcs_ccsi + v2xcs_ccsi_lag1_cumsum + 
         (1 | gwcode)),
    data = dat_5000,
    family = gaussian(),
    prior = oda_settings$prior_out,
    chains = oda_settings$chains, iter = oda_settings$iter * 2,
    warmup = oda_settings$warmup, seed = oda_settings$seed
  )
  
  return(lst(model_100, model_500, model_1000, model_5000))
}

f_oda_outcome_repress <- function(dat) {
  oda_settings <- oda_setup()
  
  dat <- dat %>%
    mutate(across(v2csreprss_lag1_cumsum, my_scale))
  
  dat_100 <- dat %>% mutate(iptw = ifelse(iptw > 100, 100, iptw))
  dat_500 <- dat %>% mutate(iptw = ifelse(iptw > 500, 500, iptw))
  dat_1000 <- dat %>% mutate(iptw = ifelse(iptw > 1000, 1000, iptw))
  dat_5000 <- dat %>% mutate(iptw = ifelse(iptw > 5000, 5000, iptw))
  
  model_100 <- brm(
    bf(total_oda_log_lead1 | weights(iptw) ~ 
         v2csreprss + v2csreprss_lag1_cumsum + 
         (1 | gwcode)),
    data = dat_100,
    family = gaussian(),
    prior = oda_settings$prior_out,
    chains = oda_settings$chains, iter = oda_settings$iter * 2,
    warmup = oda_settings$warmup, seed = oda_settings$seed
  )
  
  model_500 <- brm(
    bf(total_oda_log_lead1 | weights(iptw) ~ 
         v2csreprss + v2csreprss_lag1_cumsum + 
         (1 | gwcode)),
    data = dat_500,
    family = gaussian(),
    prior = oda_settings$prior_out,
    chains = oda_settings$chains, iter = oda_settings$iter * 2,
    warmup = oda_settings$warmup, seed = oda_settings$seed
  )
  
  model_1000 <- brm(
    bf(total_oda_log_lead1 | weights(iptw) ~ 
         v2csreprss + v2csreprss_lag1_cumsum + 
         (1 | gwcode)),
    data = dat_1000,
    family = gaussian(),
    prior = oda_settings$prior_out,
    chains = oda_settings$chains, iter = oda_settings$iter * 2,
    warmup = oda_settings$warmup, seed = oda_settings$seed
  )
  
  model_5000 <- brm(
    bf(total_oda_log_lead1 | weights(iptw) ~ 
         v2csreprss + v2csreprss_lag1_cumsum + 
         (1 | gwcode)),
    data = dat_5000,
    family = gaussian(),
    prior = oda_settings$prior_out,
    chains = oda_settings$chains, iter = oda_settings$iter * 2,
    warmup = oda_settings$warmup, seed = oda_settings$seed
  )
  
  return(lst(model_100, model_500, model_1000, model_5000))
}


# Example models for demonstrating lognormal interpretation ---------------

f_example_normal <- function(dat) {
  oda_settings <- oda_setup()
  
  dat <- dat %>% filter(laws)
  
  model <- brm(
    bf(total_oda_log_lead1 | weights(iptw) ~ 
         barriers_total + barriers_total_lag1_cumsum + 
         (1 | gwcode) + (1 | year)),
    data = dat,
    family = gaussian(), 
    chains = oda_settings$chains, iter = oda_settings$iter,
    warmup = oda_settings$warmup, seed = oda_settings$seed
  )
  
  return(model)
}

f_example_hurdle <- function(dat) {
  oda_settings <- oda_setup()
  
  dat <- dat %>% filter(laws)
  
  model <- brm(
    bf(total_oda_lead1 | weights(iptw) ~ 
         barriers_total + barriers_total_lag1_cumsum + 
         (1 | gwcode) + (1 | year), 
       hu ~ 1),
    data = dat,
    family = hurdle_lognormal(), backend = "rstan",
    chains = oda_settings$chains, iter = oda_settings$iter * 2,
    warmup = oda_settings$warmup, seed = oda_settings$seed
  )
  
  return(model)
}

R/models_purpose.R

# Settings ----------------------------------------------------------------

purpose_setup <- function() {
  options(worker_options)
  
  # Settings
  CHAINS <- 4
  ITER <- 2000
  WARMUP <- 1000
  BAYES_SEED <- 3246  # From random.org
  
  # Priors
  prior_num <- c(set_prior("normal(0, 10)", class = "Intercept"),
                 set_prior("normal(0, 2.5)", class = "b"),
                 set_prior("cauchy(0, 1)", class = "sd"))
  
  prior_denom <- c(set_prior("normal(0, 10)", class = "Intercept"),
                   set_prior("normal(0, 2.5)", class = "b"),
                   set_prior("normal(0, 2.5)", class = "sd"))
  
  prior_out <- c(set_prior("normal(0, 10)", class = "Intercept"),
                 set_prior("normal(0, 3)", class = "b"),
                 set_prior("cauchy(0, 1)", class = "sd"),
                 set_prior("logistic(-1.5, 0.5)", class = "Intercept", dpar = "zi"),
                 set_prior("gamma(0.01, 0.01)", class = "phi"))
  
  return(list(chains = CHAINS, iter = ITER, warmup = WARMUP, seed = BAYES_SEED,
              prior_num = prior_num, prior_denom = prior_denom, prior_out = prior_out))
}


# Treatment models --------------------------------------------------------

f_purpose_treatment_total <- function(dat) {
  purpose_settings <- purpose_setup()
  
  dat <- dat %>% 
    filter(laws) %>% 
    mutate(across(barriers_total_lag2_cumsum, my_scale))
  
  model_num <- brm(
    bf(barriers_total ~ barriers_total_lag1 + barriers_total_lag2_cumsum + (1 | gwcode)),
    data = dat,
    family = gaussian(),
    prior = purpose_settings$prior_num,
    control = list(adapt_delta = 0.99),
    chains = purpose_settings$chains, iter = purpose_settings$iter,
    warmup = purpose_settings$warmup, seed = purpose_settings$seed
  )
  
  model_denom <- brm(
    bf(barriers_total ~ barriers_total_lag1 + barriers_total_lag2_cumsum + prop_contentious_lag1 +
         # Human rights and politics
         v2x_polyarchy + v2x_corr + v2x_rule + v2x_civlib + v2x_clphy + v2x_clpriv +
         # Economics and development
         gdpcap_log + un_trade_pct_gdp + v2peedueq + v2pehealth + e_peinfmor +
         # Conflict and disasters
         internal_conflict_past_5 + natural_dis_count +
         (1 | gwcode)),
    data = dat,
    family = gaussian(),
    prior = purpose_settings$prior_num,
    control = list(adapt_delta = 0.9),
    chains = purpose_settings$chains, iter = purpose_settings$iter,
    warmup = purpose_settings$warmup, seed = purpose_settings$seed
  )
  
  return(lst(model_num, model_denom))
}

f_purpose_treatment_logit_total <- function(dat) {
  purpose_settings <- purpose_setup()
  
  dat <- dat %>% 
    filter(laws) %>% 
    mutate(across(barriers_total_lag2_cumsum, my_scale))
  
  model_num <- brm(
    bf(barriers_total ~ barriers_total_lag1 + barriers_total_lag2_cumsum + (1 | gwcode)),
    data = dat,
    family = gaussian(),
    prior = purpose_settings$prior_num,
    control = list(adapt_delta = 0.99),
    chains = purpose_settings$chains, iter = purpose_settings$iter,
    warmup = purpose_settings$warmup, seed = purpose_settings$seed
  )
  
  model_denom <- brm(
    bf(barriers_total ~ barriers_total_lag1 + barriers_total_lag2_cumsum + prop_contentious_logit_lag1 +
         # Human rights and politics
         v2x_polyarchy + v2x_corr + v2x_rule + v2x_civlib + v2x_clphy + v2x_clpriv +
         # Economics and development
         gdpcap_log + un_trade_pct_gdp + v2peedueq + v2pehealth + e_peinfmor +
         # Conflict and disasters
         internal_conflict_past_5 + natural_dis_count +
         (1 | gwcode)),
    data = dat,
    family = gaussian(),
    prior = purpose_settings$prior_num,
    control = list(adapt_delta = 0.9),
    chains = purpose_settings$chains, iter = purpose_settings$iter,
    warmup = purpose_settings$warmup, seed = purpose_settings$seed
  )
  
  return(lst(model_num, model_denom))
}

f_purpose_treatment_advocacy <- function(dat) {
  purpose_settings <- purpose_setup()
  
  dat <- dat %>% 
    filter(laws) %>% 
    mutate(across(advocacy_lag2_cumsum, my_scale))
  
  model_num <- brm(
    bf(advocacy ~ advocacy_lag1 + advocacy_lag2_cumsum + (1 | gwcode)),
    data = dat,
    family = gaussian(),
    prior = purpose_settings$prior_num,
    control = list(adapt_delta = 0.99),
    chains = purpose_settings$chains, iter = purpose_settings$iter,
    warmup = purpose_settings$warmup, seed = purpose_settings$seed
  )
  
  model_denom <- brm(
    bf(advocacy ~ advocacy_lag1 + advocacy_lag2_cumsum + prop_contentious_lag1 +
         v2x_polyarchy + v2x_corr + v2x_rule + v2x_civlib + v2x_clphy + v2x_clpriv +
         gdpcap_log + un_trade_pct_gdp + v2peedueq + v2pehealth + e_peinfmor +
         internal_conflict_past_5 + natural_dis_count +
         (1 | gwcode)),
    data = dat,
    family = gaussian(),
    prior = purpose_settings$prior_num,
    control = list(adapt_delta = 0.9),
    chains = purpose_settings$chains, iter = purpose_settings$iter,
    warmup = purpose_settings$warmup, seed = purpose_settings$seed
  )
  
  return(lst(model_num, model_denom))
}

f_purpose_treatment_logit_advocacy <- function(dat) {
  purpose_settings <- purpose_setup()
  
  dat <- dat %>% 
    filter(laws) %>% 
    mutate(across(advocacy_lag2_cumsum, my_scale))
  
  model_num <- brm(
    bf(advocacy ~ advocacy_lag1 + advocacy_lag2_cumsum + (1 | gwcode)),
    data = dat,
    family = gaussian(),
    prior = purpose_settings$prior_num,
    control = list(adapt_delta = 0.99),
    chains = purpose_settings$chains, iter = purpose_settings$iter,
    warmup = purpose_settings$warmup, seed = purpose_settings$seed
  )
  
  model_denom <- brm(
    bf(advocacy ~ advocacy_lag1 + advocacy_lag2_cumsum + prop_contentious_logit_lag1 +
         v2x_polyarchy + v2x_corr + v2x_rule + v2x_civlib + v2x_clphy + v2x_clpriv +
         gdpcap_log + un_trade_pct_gdp + v2peedueq + v2pehealth + e_peinfmor +
         internal_conflict_past_5 + natural_dis_count +
         (1 | gwcode)),
    data = dat,
    family = gaussian(),
    prior = purpose_settings$prior_num,
    control = list(adapt_delta = 0.9),
    chains = purpose_settings$chains, iter = purpose_settings$iter,
    warmup = purpose_settings$warmup, seed = purpose_settings$seed
  )
  
  return(lst(model_num, model_denom))
}

f_purpose_treatment_entry <- function(dat) {
  purpose_settings <- purpose_setup()
  
  dat <- dat %>% 
    filter(laws) %>% 
    mutate(across(entry_lag2_cumsum, my_scale))
  
  model_num <- brm(
    bf(entry ~ entry_lag1 + entry_lag2_cumsum + (1 | gwcode)),
    data = dat,
    family = gaussian(),
    prior = purpose_settings$prior_num,
    control = list(adapt_delta = 0.99),
    chains = purpose_settings$chains, iter = purpose_settings$iter,
    warmup = purpose_settings$warmup, seed = purpose_settings$seed
  )
  
  model_denom <- brm(
    bf(entry ~ entry_lag1 + entry_lag2_cumsum + prop_contentious_lag1 +
         v2x_polyarchy + v2x_corr + v2x_rule + v2x_civlib + v2x_clphy + v2x_clpriv +
         gdpcap_log + un_trade_pct_gdp + v2peedueq + v2pehealth + e_peinfmor +
         internal_conflict_past_5 + natural_dis_count +
         (1 | gwcode)),
    data = dat,
    family = gaussian(),
    prior = purpose_settings$prior_num,
    control = list(adapt_delta = 0.9),
    chains = purpose_settings$chains, iter = purpose_settings$iter,
    warmup = purpose_settings$warmup, seed = purpose_settings$seed
  )
  
  return(lst(model_num, model_denom))
}

f_purpose_treatment_funding <- function(dat) {
  purpose_settings <- purpose_setup()
  
  dat <- dat %>% 
    filter(laws) %>% 
    mutate(across(funding_lag2_cumsum, my_scale))
  
  model_num <- brm(
    bf(funding ~ funding_lag1 + funding_lag2_cumsum + (1 | gwcode)),
    data = dat,
    family = gaussian(),
    prior = purpose_settings$prior_num,
    control = list(adapt_delta = 0.99),
    chains = purpose_settings$chains, iter = purpose_settings$iter,
    warmup = purpose_settings$warmup, seed = purpose_settings$seed
  )
  
  model_denom <- brm(
    bf(funding ~ funding_lag1 + funding_lag2_cumsum + prop_contentious_lag1 +
         v2x_polyarchy + v2x_corr + v2x_rule + v2x_civlib + v2x_clphy + v2x_clpriv +
         gdpcap_log + un_trade_pct_gdp + v2peedueq + v2pehealth + e_peinfmor +
         internal_conflict_past_5 + natural_dis_count +
         (1 | gwcode)),
    data = dat,
    family = gaussian(),
    prior = purpose_settings$prior_num,
    control = list(adapt_delta = 0.9),
    chains = purpose_settings$chains, iter = purpose_settings$iter,
    warmup = purpose_settings$warmup, seed = purpose_settings$seed
  )
  
  return(lst(model_num, model_denom))
}


# Outcome models ----------------------------------------------------------

f_purpose_outcome_total <- function(dat) {
  purpose_settings <- purpose_setup()
  
  dat <- dat %>% 
    filter(laws) %>% 
    mutate(across(barriers_total_lag1_cumsum, my_scale))
  
  model <- brm(
    bf(prop_contentious_lead1 | weights(iptw) ~ 
         barriers_total + barriers_total_lag1_cumsum + 
         (1 | gwcode),
       zi ~ 1),
    data = dat,
    family = zero_inflated_beta(),
    inits = 0,
    prior = purpose_settings$prior_out,
    chains = purpose_settings$chains, iter = purpose_settings$iter * 2,
    warmup = purpose_settings$warmup, seed = purpose_settings$seed
  )
  
  return(model)
}

f_purpose_outcome_logit_total <- function(dat) {
  purpose_settings <- purpose_setup()
  
  dat <- dat %>% 
    filter(laws) %>% 
    mutate(across(barriers_total_lag1_cumsum, my_scale))
  
  model <- brm(
    bf(prop_contentious_logit_lead1 | weights(iptw) ~ 
         barriers_total + barriers_total_lag1_cumsum + 
         (1 | gwcode)),
    data = dat,
    family = gaussian(),
    inits = 0,
    prior = c(set_prior("normal(0, 10)", class = "Intercept"),
              set_prior("normal(0, 2.5)", class = "b"),
              set_prior("normal(0, 2.5)", class = "sd")),
    chains = purpose_settings$chains, iter = purpose_settings$iter * 2,
    warmup = purpose_settings$warmup, seed = purpose_settings$seed
  )
  
  return(model)
}

f_purpose_outcome_advocacy <- function(dat) {
  purpose_settings <- purpose_setup()
  
  dat <- dat %>% 
    filter(laws) %>% 
    mutate(across(advocacy_lag1_cumsum, my_scale))
  
  model <- brm(
    bf(prop_contentious_lead1 | weights(iptw) ~ 
         advocacy + advocacy_lag1_cumsum + 
         (1 | gwcode) + (1 | year),
       zi ~ 1),
    data = dat,
    family = zero_inflated_beta(),
    inits = 0,
    prior = purpose_settings$prior_out,
    chains = purpose_settings$chains, iter = purpose_settings$iter * 2,
    warmup = purpose_settings$warmup, seed = purpose_settings$seed
  )
  
  return(model)
}

f_purpose_outcome_logit_advocacy <- function(dat) {
  purpose_settings <- purpose_setup()
  
  dat <- dat %>% 
    filter(laws) %>% 
    mutate(across(advocacy_lag1_cumsum, my_scale))
  
  model <- brm(
    bf(prop_contentious_logit_lead1 | weights(iptw) ~ 
         advocacy + advocacy_lag1_cumsum + 
         (1 | gwcode) + (1 | year)),
    data = dat,
    family = gaussian(),
    inits = 0,
    prior = c(set_prior("normal(0, 10)", class = "Intercept"),
              set_prior("normal(0, 2.5)", class = "b"),
              set_prior("normal(0, 2.5)", class = "sd")),
    chains = purpose_settings$chains, iter = purpose_settings$iter * 2,
    warmup = purpose_settings$warmup, seed = purpose_settings$seed
  )
  
  return(model)
}

f_purpose_outcome_entry <- function(dat) {
  purpose_settings <- purpose_setup()
  
  dat <- dat %>% 
    filter(laws) %>% 
    mutate(across(entry_lag1_cumsum, my_scale))
  
  model <- brm(
    bf(prop_contentious_lead1 | weights(iptw) ~ 
         entry + entry_lag1_cumsum + 
         (1 | gwcode) + (1 | year),
       zi ~ 1),
    data = dat,
    family = zero_inflated_beta(),
    inits = 0,
    prior = purpose_settings$prior_out,
    chains = purpose_settings$chains, iter = purpose_settings$iter * 2,
    warmup = purpose_settings$warmup, seed = purpose_settings$seed
  )
  
  return(model)
}

f_purpose_outcome_funding <- function(dat) {
  purpose_settings <- purpose_setup()
  
  dat <- dat %>% 
    filter(laws) %>% 
    mutate(across(funding_lag1_cumsum, my_scale))
  
  model <- brm(
    bf(prop_contentious_lead1 | weights(iptw) ~ 
         funding + funding_lag1_cumsum + 
         (1 | gwcode) + (1 | year),
       zi ~ 1),
    data = dat,
    family = zero_inflated_beta(),
    inits = 0,
    prior = purpose_settings$prior_out,
    chains = purpose_settings$chains, iter = purpose_settings$iter * 2,
    warmup = purpose_settings$warmup, seed = purpose_settings$seed
  )
  
  return(model)
}

R/models_recipients.R

# Settings ----------------------------------------------------------------

recip_setup <- function() {
  options(worker_options)
  
  # Settings
  CHAINS <- 4
  ITER <- 2000
  WARMUP <- 1000
  BAYES_SEED <- 4045  # From random.org
  
  # Priors
  prior_num <- c(set_prior("normal(0, 10)", class = "Intercept"),
                 set_prior("normal(0, 2.5)", class = "b"),
                 set_prior("cauchy(0, 1)", class = "sd"))
  
  prior_denom <- c(set_prior("normal(0, 10)", class = "Intercept"),
                   set_prior("normal(0, 2.5)", class = "b"),
                   set_prior("normal(0, 2.5)", class = "sd"))
  
  prior_out <- c(set_prior("normal(0, 10)", class = "Intercept"),
                 set_prior("normal(0, 3)", class = "b"),
                 set_prior("cauchy(0, 1)", class = "sd"),
                 set_prior("logistic(-0.5, 0.35)", class = "Intercept", dpar = "zi"),
                 set_prior("gamma(0.01, 0.01)", class = "phi"))
  
  return(list(chains = CHAINS, iter = ITER, warmup = WARMUP, seed = BAYES_SEED,
              prior_num = prior_num, prior_denom = prior_denom, prior_out = prior_out))
}


# Treatment models --------------------------------------------------------

f_recip_treatment_total_dom <- function(dat) {
  recip_settings <- recip_setup()
  
  dat <- dat %>% 
    filter(laws) %>% 
    mutate(across(barriers_total_lag2_cumsum, my_scale))
  
  model_num <- brm(
    bf(barriers_total ~ barriers_total_lag1 + barriers_total_lag2_cumsum + (1 | gwcode)),
    data = dat,
    family = gaussian(),
    prior = recip_settings$prior_num,
    control = list(adapt_delta = 0.99),
    chains = recip_settings$chains, iter = recip_settings$iter,
    warmup = recip_settings$warmup, seed = recip_settings$seed
  )
  
  model_denom <- brm(
    bf(barriers_total ~ barriers_total_lag1 + barriers_total_lag2_cumsum + prop_ngo_dom_lag1 +
         # Human rights and politics
         v2x_polyarchy + v2x_corr + v2x_rule + v2x_civlib + v2x_clphy + v2x_clpriv +
         # Economics and development
         gdpcap_log + un_trade_pct_gdp + v2peedueq + v2pehealth + e_peinfmor +
         # Conflict and disasters
         internal_conflict_past_5 + natural_dis_count +
         (1 | gwcode)),
    data = dat,
    family = gaussian(),
    prior = recip_settings$prior_num,
    control = list(adapt_delta = 0.9),
    chains = recip_settings$chains, iter = recip_settings$iter,
    warmup = recip_settings$warmup, seed = recip_settings$seed
  )
  
  return(lst(model_num, model_denom))
}

f_recip_treatment_total_foreign <- function(dat) {
  recip_settings <- recip_setup()
  
  dat <- dat %>% 
    filter(laws) %>% 
    mutate(across(barriers_total_lag2_cumsum, my_scale))
  
  model_num <- brm(
    bf(barriers_total ~ barriers_total_lag1 + barriers_total_lag2_cumsum + (1 | gwcode)),
    data = dat,
    family = gaussian(),
    prior = recip_settings$prior_num,
    control = list(adapt_delta = 0.99),
    chains = recip_settings$chains, iter = recip_settings$iter,
    warmup = recip_settings$warmup, seed = recip_settings$seed
  )
  
  model_denom <- brm(
    bf(barriers_total ~ barriers_total_lag1 + barriers_total_lag2_cumsum + prop_ngo_foreign_lag1 +
         v2x_polyarchy + v2x_corr + v2x_rule + v2x_civlib + v2x_clphy + v2x_clpriv +
         gdpcap_log + un_trade_pct_gdp + v2peedueq + v2pehealth + e_peinfmor +
         internal_conflict_past_5 + natural_dis_count +
         (1 | gwcode)),
    data = dat,
    family = gaussian(),
    prior = recip_settings$prior_num,
    control = list(adapt_delta = 0.9),
    chains = recip_settings$chains, iter = recip_settings$iter,
    warmup = recip_settings$warmup, seed = recip_settings$seed
  )
  
  return(lst(model_num, model_denom))
}

f_recip_treatment_advocacy_dom <- function(dat) {
  recip_settings <- recip_setup()
  
  dat <- dat %>% 
    filter(laws) %>% 
    mutate(across(advocacy_lag2_cumsum, my_scale))
  
  model_num <- brm(
    bf(advocacy ~ advocacy_lag1 + advocacy_lag2_cumsum + (1 | gwcode)),
    data = dat,
    family = gaussian(),
    prior = recip_settings$prior_num,
    control = list(adapt_delta = 0.99),
    chains = recip_settings$chains, iter = recip_settings$iter,
    warmup = recip_settings$warmup, seed = recip_settings$seed
  )
  
  model_denom <- brm(
    bf(advocacy ~ advocacy_lag1 + advocacy_lag2_cumsum + prop_ngo_dom_lag1 +
         v2x_polyarchy + v2x_corr + v2x_rule + v2x_civlib + v2x_clphy + v2x_clpriv +
         gdpcap_log + un_trade_pct_gdp + v2peedueq + v2pehealth + e_peinfmor +
         internal_conflict_past_5 + natural_dis_count +
         (1 | gwcode)),
    data = dat,
    family = gaussian(),
    prior = recip_settings$prior_num,
    control = list(adapt_delta = 0.9),
    chains = recip_settings$chains, iter = recip_settings$iter,
    warmup = recip_settings$warmup, seed = recip_settings$seed
  )
  
  return(lst(model_num, model_denom))
}

f_recip_treatment_advocacy_foreign <- function(dat) {
  recip_settings <- recip_setup()
  
  dat <- dat %>% 
    filter(laws) %>% 
    mutate(across(advocacy_lag2_cumsum, my_scale))
  
  model_num <- brm(
    bf(advocacy ~ advocacy_lag1 + advocacy_lag2_cumsum + (1 | gwcode)),
    data = dat,
    family = gaussian(),
    prior = recip_settings$prior_num,
    control = list(adapt_delta = 0.99),
    chains = recip_settings$chains, iter = recip_settings$iter,
    warmup = recip_settings$warmup, seed = recip_settings$seed
  )
  
  model_denom <- brm(
    bf(advocacy ~ advocacy_lag1 + advocacy_lag2_cumsum + prop_ngo_foreign_lag1 +
         v2x_polyarchy + v2x_corr + v2x_rule + v2x_civlib + v2x_clphy + v2x_clpriv +
         gdpcap_log + un_trade_pct_gdp + v2peedueq + v2pehealth + e_peinfmor +
         internal_conflict_past_5 + natural_dis_count +
         (1 | gwcode)),
    data = dat,
    family = gaussian(),
    prior = recip_settings$prior_num,
    control = list(adapt_delta = 0.9),
    chains = recip_settings$chains, iter = recip_settings$iter,
    warmup = recip_settings$warmup, seed = recip_settings$seed
  )
  
  return(lst(model_num, model_denom))
}

f_recip_treatment_entry_dom <- function(dat) {
  recip_settings <- recip_setup()
  
  dat <- dat %>% 
    filter(laws) %>% 
    mutate(across(entry_lag2_cumsum, my_scale))
  
  model_num <- brm(
    bf(entry ~ entry_lag1 + entry_lag2_cumsum + (1 | gwcode)),
    data = dat,
    family = gaussian(),
    prior = recip_settings$prior_num,
    control = list(adapt_delta = 0.99),
    chains = recip_settings$chains, iter = recip_settings$iter,
    warmup = recip_settings$warmup, seed = recip_settings$seed
  )
  
  model_denom <- brm(
    bf(entry ~ entry_lag1 + entry_lag2_cumsum + prop_ngo_dom_lag1 +
         v2x_polyarchy + v2x_corr + v2x_rule + v2x_civlib + v2x_clphy + v2x_clpriv +
         gdpcap_log + un_trade_pct_gdp + v2peedueq + v2pehealth + e_peinfmor +
         internal_conflict_past_5 + natural_dis_count +
         (1 | gwcode)),
    data = dat,
    family = gaussian(),
    prior = recip_settings$prior_num,
    control = list(adapt_delta = 0.9),
    chains = recip_settings$chains, iter = recip_settings$iter,
    warmup = recip_settings$warmup, seed = recip_settings$seed
  )
  
  return(lst(model_num, model_denom))
}

f_recip_treatment_entry_foreign <- function(dat) {
  recip_settings <- recip_setup()
  
  dat <- dat %>% 
    filter(laws) %>% 
    mutate(across(entry_lag2_cumsum, my_scale))
  
  model_num <- brm(
    bf(entry ~ entry_lag1 + entry_lag2_cumsum + (1 | gwcode)),
    data = dat,
    family = gaussian(),
    prior = recip_settings$prior_num,
    control = list(adapt_delta = 0.99),
    chains = recip_settings$chains, iter = recip_settings$iter,
    warmup = recip_settings$warmup, seed = recip_settings$seed
  )
  
  model_denom <- brm(
    bf(entry ~ entry_lag1 + entry_lag2_cumsum + prop_ngo_foreign_lag1 +
         v2x_polyarchy + v2x_corr + v2x_rule + v2x_civlib + v2x_clphy + v2x_clpriv +
         gdpcap_log + un_trade_pct_gdp + v2peedueq + v2pehealth + e_peinfmor +
         internal_conflict_past_5 + natural_dis_count +
         (1 | gwcode)),
    data = dat,
    family = gaussian(),
    prior = recip_settings$prior_num,
    control = list(adapt_delta = 0.9),
    chains = recip_settings$chains, iter = recip_settings$iter,
    warmup = recip_settings$warmup, seed = recip_settings$seed
  )
  
  return(lst(model_num, model_denom))
}

f_recip_treatment_funding_dom <- function(dat) {
  recip_settings <- recip_setup()
  
  dat <- dat %>% 
    filter(laws) %>% 
    mutate(across(funding_lag2_cumsum, my_scale))
  
  model_num <- brm(
    bf(funding ~ funding_lag1 + funding_lag2_cumsum + (1 | gwcode)),
    data = dat,
    family = gaussian(),
    prior = recip_settings$prior_num,
    control = list(adapt_delta = 0.99),
    chains = recip_settings$chains, iter = recip_settings$iter,
    warmup = recip_settings$warmup, seed = recip_settings$seed
  )
  
  model_denom <- brm(
    bf(funding ~ funding_lag1 + funding_lag2_cumsum + prop_ngo_dom_lag1 +
         v2x_polyarchy + v2x_corr + v2x_rule + v2x_civlib + v2x_clphy + v2x_clpriv +
         gdpcap_log + un_trade_pct_gdp + v2peedueq + v2pehealth + e_peinfmor +
         internal_conflict_past_5 + natural_dis_count +
         (1 | gwcode)),
    data = dat,
    family = gaussian(),
    prior = recip_settings$prior_num,
    control = list(adapt_delta = 0.9),
    chains = recip_settings$chains, iter = recip_settings$iter,
    warmup = recip_settings$warmup, seed = recip_settings$seed
  )
  
  return(lst(model_num, model_denom))
}

f_recip_treatment_funding_foreign <- function(dat) {
  recip_settings <- recip_setup()
  
  dat <- dat %>% 
    filter(laws) %>% 
    mutate(across(funding_lag2_cumsum, my_scale))
  
  model_num <- brm(
    bf(funding ~ funding_lag1 + funding_lag2_cumsum + (1 | gwcode)),
    data = dat,
    family = gaussian(),
    prior = recip_settings$prior_num,
    control = list(adapt_delta = 0.99),
    chains = recip_settings$chains, iter = recip_settings$iter,
    warmup = recip_settings$warmup, seed = recip_settings$seed
  )
  
  model_denom <- brm(
    bf(funding ~ funding_lag1 + funding_lag2_cumsum + prop_ngo_foreign_lag1 +
         v2x_polyarchy + v2x_corr + v2x_rule + v2x_civlib + v2x_clphy + v2x_clpriv +
         gdpcap_log + un_trade_pct_gdp + v2peedueq + v2pehealth + e_peinfmor +
         internal_conflict_past_5 + natural_dis_count +
         (1 | gwcode)),
    data = dat,
    family = gaussian(),
    prior = recip_settings$prior_num,
    control = list(adapt_delta = 0.9),
    chains = recip_settings$chains, iter = recip_settings$iter,
    warmup = recip_settings$warmup, seed = recip_settings$seed
  )
  
  return(lst(model_num, model_denom))
}


# Outcome models ----------------------------------------------------------

f_recip_outcome_total_dom <- function(dat) {
  recip_settings <- recip_setup()
  
  dat <- dat %>% 
    filter(laws) %>% 
    mutate(across(barriers_total_lag1_cumsum, my_scale))
  
  model <- brm(
    bf(prop_ngo_dom_lead1 | weights(iptw) ~ 
         barriers_total + barriers_total_lag1_cumsum + 
         (1 | gwcode) + (1 | year),
       zi ~ 1),
    data = dat,
    family = zero_inflated_beta(),
    prior = recip_settings$prior_out,
    chains = recip_settings$chains, iter = recip_settings$iter * 2,
    warmup = recip_settings$warmup, seed = recip_settings$seed
  )
  
  return(model)
}

f_recip_outcome_total_foreign <- function(dat) {
  recip_settings <- recip_setup()
  
  dat <- dat %>% 
    filter(laws) %>% 
    mutate(across(barriers_total_lag1_cumsum, my_scale))
  
  model <- brm(
    bf(prop_ngo_foreign_lead1 | weights(iptw) ~ 
         barriers_total + barriers_total_lag1_cumsum + 
         (1 | gwcode) + (1 | year),
       zi ~ 1),
    data = dat,
    family = zero_inflated_beta(),
    prior = recip_settings$prior_out,
    chains = recip_settings$chains, iter = recip_settings$iter * 2,
    warmup = recip_settings$warmup, seed = recip_settings$seed
  )
  
  return(model)
}

f_recip_outcome_advocacy_dom <- function(dat) {
  recip_settings <- recip_setup()
  
  dat <- dat %>% 
    filter(laws) %>% 
    mutate(across(advocacy_lag1_cumsum, my_scale))
  
  model <- brm(
    bf(prop_ngo_dom_lead1 | weights(iptw) ~ 
         advocacy + advocacy_lag1_cumsum + 
         (1 | gwcode) + (1 | year),
       zi ~ 1),
    data = dat,
    family = zero_inflated_beta(),
    prior = recip_settings$prior_out,
    chains = recip_settings$chains, iter = recip_settings$iter * 2,
    warmup = recip_settings$warmup, seed = recip_settings$seed
  )
  
  return(model)
}

f_recip_outcome_advocacy_foreign <- function(dat) {
  recip_settings <- recip_setup()
  
  dat <- dat %>% 
    filter(laws) %>% 
    mutate(across(advocacy_lag1_cumsum, my_scale))
  
  model <- brm(
    bf(prop_ngo_foreign_lead1 | weights(iptw) ~ 
         advocacy + advocacy_lag1_cumsum + 
         (1 | gwcode) + (1 | year),
       zi ~ 1),
    data = dat,
    family = zero_inflated_beta(),
    prior = recip_settings$prior_out,
    chains = recip_settings$chains, iter = recip_settings$iter * 2,
    warmup = recip_settings$warmup, seed = recip_settings$seed
  )
  
  return(model)
}

f_recip_outcome_entry_dom <- function(dat) {
  recip_settings <- recip_setup()
  
  dat <- dat %>% 
    filter(laws) %>% 
    mutate(across(entry_lag1_cumsum, my_scale))
  
  model <- brm(
    bf(prop_ngo_dom_lead1 | weights(iptw) ~ 
         entry + entry_lag1_cumsum + 
         (1 | gwcode) + (1 | year),
       zi ~ 1),
    data = dat,
    family = zero_inflated_beta(),
    prior = recip_settings$prior_out,
    chains = recip_settings$chains, iter = recip_settings$iter * 2,
    warmup = recip_settings$warmup, seed = recip_settings$seed
  )
  
  return(model)
}

f_recip_outcome_entry_foreign <- function(dat) {
  recip_settings <- recip_setup()
  
  dat <- dat %>% 
    filter(laws) %>% 
    mutate(across(entry_lag1_cumsum, my_scale))
  
  model <- brm(
    bf(prop_ngo_foreign_lead1 | weights(iptw) ~ 
         entry + entry_lag1_cumsum + 
         (1 | gwcode) + (1 | year),
       zi ~ 1),
    data = dat,
    family = zero_inflated_beta(),
    prior = recip_settings$prior_out,
    chains = recip_settings$chains, iter = recip_settings$iter * 2,
    warmup = recip_settings$warmup, seed = recip_settings$seed
  )
  
  return(model)
}

f_recip_outcome_funding_dom <- function(dat) {
  recip_settings <- recip_setup()
  
  dat <- dat %>% 
    filter(laws) %>% 
    mutate(across(funding_lag1_cumsum, my_scale))
  
  model <- brm(
    bf(prop_ngo_dom_lead1 | weights(iptw) ~ 
         funding + funding_lag1_cumsum + 
         (1 | gwcode) + (1 | year),
       zi ~ 1),
    data = dat,
    family = zero_inflated_beta(),
    prior = recip_settings$prior_out,
    chains = recip_settings$chains, iter = recip_settings$iter * 2,
    warmup = recip_settings$warmup, seed = recip_settings$seed
  )
  
  return(model)
}

f_recip_outcome_funding_foreign <- function(dat) {
  recip_settings <- recip_setup()
  
  dat <- dat %>% 
    filter(laws) %>% 
    mutate(across(funding_lag1_cumsum, my_scale))
  
  model <- brm(
    bf(prop_ngo_foreign_lead1 | weights(iptw) ~ 
         funding + funding_lag1_cumsum + 
         (1 | gwcode) + (1 | year),
       zi ~ 1),
    data = dat,
    family = zero_inflated_beta(),
    prior = recip_settings$prior_out,
    chains = recip_settings$chains, iter = recip_settings$iter * 2,
    warmup = recip_settings$warmup, seed = recip_settings$seed
  )
  
  return(model)
}

R/funs_models-iptw.R

cumprod_na <- function(x) {
  x[is.na(x)] <- 1
  return(cumprod(x))
}

cumsum_na <- function(x) {
  x[is.na(x)] <- 0
  return(cumsum(x))
}

# Via https://stackoverflow.com/a/55323097/120898
lhs <- function(x) {
  if (attr(terms(as.formula(x)), which = "response")) {
    all.vars(x)[1]
  } else {
    NULL
  }
}

# Base R's scale() outputs a matrix, which requires a bunch of shenanigans to
# wrangle into a column, like mutate(x = as.numeric(scale(x))). So instead we
# just use our own here
my_scale <- function(x) {
  (x - mean(x, na.rm = TRUE)) / sd(x, na.rm = TRUE)
}

create_iptws <- function(dat, wt_model) {
  pred_num <- predict(wt_model$model_num, newdata = dat, allow_new_levels = TRUE)
  resid_num <- residuals(wt_model$model_num, newdata = dat, allow_new_levels = TRUE)
  lhs_num <- lhs(wt_model$model_num$formula$formula)
  
  num_actual <- dnorm(dat[[lhs_num]],
                      pred_num[,1],
                      sd(resid_num[,1], na.rm = TRUE))
  
  pred_denom <- predict(wt_model$model_denom, newdata = dat, allow_new_levels = TRUE)
  resid_denom <- residuals(wt_model$model_denom, newdata = dat, allow_new_levels = TRUE)
  lhs_denom <- lhs(wt_model$model_denom$formula$formula)
  
  denom_actual <- dnorm(dat[[lhs_denom]],
                        pred_denom[,1],
                        sd(resid_denom[,1], na.rm = TRUE))
  
  dat <- dat %>% 
    mutate(weights_sans_time = num_actual / denom_actual) %>% 
    group_by(gwcode) %>% 
    mutate(iptw = cumprod_na(weights_sans_time)) %>% 
    ungroup()
  
  return(dat)
}

Copyright © 2017–2021 Suparna Chaudhry and Andrew Heiss. Licensed under the MIT License. Code available on GitHub.