Usage example: modeling revolutionary onset

This article needs the vdemdata, brglm2, and ggplot2 packages, which are not all installed. Install them and rebuild the vignette to see the full walkthrough.

This walks through a complete workflow with contentiousR: construct a country-year panel from several sources, check its coverage, and estimate a rare-events logit of revolutionary onset. It is rather a demonstration of package mechanics than a substantive research claim — the model below is illustrative and has not been vetted as a serious specification.

Construct the panel

library(contentiousR)
library(dplyr)

panel <- build_states_panel(1900, 2010, coding_system = "cow") |>
  add_gdp(dataset = "fariss") |>
  add_population(dataset = "wpp") |>
  add_vdem(vars = c("v2x_polyarchy", "v2x_libdem", "v2x_execorr")) |>
  add_leader_data(dataset = "archigos") |>
  add_military_expenditure(dataset = "nmc") |>
  add_conflict(dataset = "beissinger", aggregate = TRUE) |>
  select(
    cow, year, country, fariss_gdppc, wpp_pop, nmc_milex,
    v2x_polyarchy, v2x_execorr, leader_tenure, beissinger_onset
  ) |>
  mutate(across(starts_with("beissinger_"), ~ tidyr::replace_na(.x, 0))) |> # fill non-onsets with zeros
  add_spell_duration() |>
  add_lag(
    vars = c(
      "fariss_gdppc", "wpp_pop", "v2x_polyarchy",
      "v2x_execorr", "leader_tenure", "nmc_milex"
    ),
    n = 1
  )

dim(panel)

head(panel) |>
  mutate(across(where(is.numeric), ~ round(.x, 2))) |>
  knitr::kable()

add_conflict(dataset = "beissinger", aggregate = TRUE) is used rather than aggregate = FALSE: Revolutionary Episodes Dataset is an episode dataset, and a handful of country-years have more than one recorded episode (e.g. Germany in 1918). Left unaggregated, those country-years would produce duplicate cow-year rows. replace_na(0) treats “not covered” as “no onset” so that add_spell_duration() (which cannot span NA) can compute a peace-years counter (beissinger_spell) across the whole panel. This is a modeling choice, not a neutral default.

Check coverage

plot_coverage(panel, "beissinger_onset", drop_empty = TRUE)

Model onset

Onsets are rare, which biases a regular logit model. The model below uses brglm2’s Jeffreys-prior penalized likelihood (method = "brglmFit", type = "MPL_Jeffreys") instead (see Kosmidis and Firth 2009), with a natural spline on beissinger_spell to flexibly control for time since the last episode and year fixed effects. One-year lagged controls are included in the estimation.

library(brglm2)

model <- glm(
  beissinger_onset ~ log(fariss_gdppc_l + 0.1) + wpp_pop_l +
    v2x_polyarchy_l + I(v2x_polyarchy_l^2) + v2x_execorr_l +
    leader_tenure_l + I(leader_tenure_l^2) + nmc_milex_l +
    splines::bs(beissinger_spell, 3) + as.factor(year),
  family = binomial("logit"),
  method = "brglmFit", type = "MPL_Jeffreys",
  data = panel
)

co <- summary(model)$coefficients
co[!grepl("factor\\(year\\)", rownames(co)), ]

Year fixed effects are omitted from the printed table above for readability.

Reference

Kosmidis, I., and D. Firth. 2009. “Bias Reduction in Exponential Family Nonlinear Models.” Biometrika 96 (4): 793–804. https://doi.org/10.1093/biomet/asp055.