## ----setup, include = FALSE---------------------------------------------------
knitr::opts_chunk$set(collapse = TRUE, comment = "#>")
has_deps <- requireNamespace("vdemdata", quietly = TRUE) &&
  requireNamespace("brglm2", quietly = TRUE) &&
  requireNamespace("ggplot2", quietly = TRUE)

## ----eval = !has_deps, echo = FALSE, results = "asis"-------------------------
# cat(
#   "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.",
#   sep = ""
# )

## ----eval = has_deps----------------------------------------------------------
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)

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

## ----eval = has_deps, fig.width = 7, fig.height = 6---------------------------
plot_coverage(panel, "beissinger_onset", drop_empty = TRUE)

## ----eval = has_deps----------------------------------------------------------
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)), ]