Update results of Kim (2021) to include binary variable for hereditary head of state

monarchic_survival.R

### Replicate 
### 
### Kim, Nam Kyu. "Previous Military Rule and Democratic Survival." Journal of Conflict Resolution 65, no. 2-3 (2021): 534-562
###
### including an additional covariate for a hereditary head of state
### 
library(rio)
library(tidyverse)
library(vdemdata)
library(countrycode)
library(brms)
library(modelsummary)
set.seed(1247)

### Cut-down the vdem data
data("vdem")
vdem <- vdem %>%
    dplyr::select(country_name,
                  iso3c = country_text_id,
                  year,
                  v2expathhs)

### Read in the Kim  data
dat <- rio::import("main_data.dta")
dat <- dat %>%
    mutate(iso3c = countrycode(cowcode, "cown", "iso3c"))

### Merge
dat <- merge(dat, vdem,
             by.x = c("iso3c", "year"),
             by.y = c("iso3c", "year"),
             all.x = TRUE,
             all.y = FALSE)


table(dat$v2expathhs, useNA = "always")

dat <- dat %>%
    mutate(is_monarchy = as.numeric(v2expathhs == 4))

### Replicate the original analysis
### specifically model one

f1 <- gwf_fail ~ leadermil_prior +
    reg_Lv2x_polyarchy + log_demyears + Le_migdppcln + Le_migdpgro +
    postcold + poly(gwf_duration, 3)

f2 <- update(f1, . ~ . + is_monarchy)

mod <- glm(f1, data = dat %>% filter(gwf_dem == 1),
           family = binomial(link = "probit"))

confint(mod)

### Restrict it just countries with non-missing values on monarchy
mod2 <- update(mod,
               data = dat %>%
                   filter(gwf_dem == 1) %>%
                   filter(!is.na(is_monarchy)))

### update it now to include is_monarchy
mod3 <- update(mod2,
               formula = f2)

### Does it look any different in a Bayesian model which might deal
### with perfect separation?

mod4 <- brm(f2,
            data = dat %>%
                filter(gwf_dem == 1) %>%
                filter(!is.na(is_monarchy)),
            family = bernoulli(link = "probit"),
            prior = set_prior("normal(0, 1)", coef = "leadermil_prior") +
                set_prior("normal(0, 1)", coef = "Le_migdpgro") +
                set_prior("normal(0, 1)", coef = "Le_migdppcln") +
                set_prior("normal(0, 1)", coef = "log_demyears") +
                set_prior("normal(0, 1)", coef = "postcold") +
                set_prior("normal(0, 1)", coef = "reg_Lv2x_polyarchy"),
            chains = 4,
            cores = 4)

cm <- c('Intercept' = 'Intercept',
        'leadermil_prior' = 'Prior military rule',
        'Le_migdpgro' = 'GDP growth',
        'Le_migdppcln' = 'log(GDP per capita)',
        'log_demyears' = 'log(Previous democracy years)',
        'postcold' = 'Post-Cold War',
        'reg_Lv2x_polyarchy' = 'Regional democracy',
        'b_leadermil_prior' = 'Prior military rule',
        'b_Le_migdpgro' = 'GDP growth',
        'b_Le_migdppcln' = 'log(GDP per capita)',
        'b_log_demyears' = 'log(Previous democracy years)',
        'b_postcold' = 'Post-Cold War',
        'b_reg_Lv2x_polyarchy' = 'Regional democracy',
        'is_monarchy' = 'Hereditary HoS',
        'b_is_monarchy' = 'Hereditary HoS'
        )

modelsummary(list("Original Kim model" = mod,
                  "Model w/ NA  removed" = mod2,
                  "Model w/ monarchy" = mod3,
                  "Bayesian model" = mod4),
             coef_map = cm,
             statistic = "conf.int")