current_dir <- normalizePath(getwd(), winslash = "/", mustWork = TRUE)
project_root <- if (basename(current_dir) == "scripts") {
  dirname(current_dir)
} else if (dir.exists(file.path(current_dir, "scripts")) &&
           dir.exists(file.path(current_dir, "data"))) {
  current_dir
} else {
  dirname(current_dir)
}
setwd(project_root)
dir.create(file.path(project_root, "data"), showWarnings = FALSE)
dir.create(file.path(project_root, "outputs"), showWarnings = FALSE)
local_lib <- file.path(project_root, "r-lib")
if (dir.exists(local_lib)) {
  .libPaths(c(local_lib, .libPaths()))
}
local_lib <- file.path(getwd(), "r-lib")
.libPaths(c(local_lib, .libPaths()))

###### Voting for democracy - Direct Linear Probability Models (LPMs), faceted plots
#### Variable by variable plots, with samples, excluding pooled, in rows

library(dplyr)
library(tidyr)
library(ggplot2)
library(lmtest)
library(sandwich)
library(texreg)
library(here)
library(grid)

setwd(project_root)
dir.create(file.path(project_root, "data"), showWarnings = FALSE)
dir.create(file.path(project_root, "outputs"), showWarnings = FALSE)

# Read analytical datasets
us   <- read.csv(file.path(project_root, "data", "us_analytic_v2.csv"))
sk   <- read.csv(file.path(project_root, "data", "sk_analytic_v2.csv"))
pl   <- read.csv(file.path(project_root, "data", "pl_analytic_v2.csv"))
pool <- read.csv(file.path(project_root, "data", "pool_analytic_v2.csv"))

# Read long conjoint datasets
us_long   <- read.csv(file.path(project_root, "data", "us_long_v2.csv"))
sk_long   <- read.csv(file.path(project_root, "data", "sk_long_v2.csv"))
pl_long   <- read.csv(file.path(project_root, "data", "pl_long_v2.csv"))
pool_long <- read.csv(file.path(project_root, "data", "pool_long_v2.csv"))

# Set up lists for looping
wide_list <- list(pool = pool, us = us, sk = sk, pl = pl)
long_list <- list(pool = pool_long, us = us_long, sk = sk_long, pl = pl_long)
samples   <- names(wide_list)
plot_samples <- c("pl", "sk", "us")

pro_dem_labels <- c("keep polling stations","respect court rulings",
                    "don't prosecute journalists")

pc_vars <- c("demaut_pref","sup_ldem","inst_trust","pol_tolerance",
             "emancip_val","autonomy_values")

covars  <- c("gender_male","age_1","educ_5f","income_cat","region_pool",
             "pol_interest","gend_pol_left","redist_pol_left","top_party")
covars_base <- setdiff(covars, "top_party")

pc_labels <- c(
  demaut_pref = "Democ-autoc. preference",
  sup_ldem = "Support liberal democracy",
  inst_trust = "Institutional trust",
  pol_tolerance = "Political tolerance",
  emancip_val = "Emancipative values",
  autonomy_values = "Authoritarian values (rev)"
)

label_map <- c(
  demaut_pref = "Democ-autoc. preference",
  sup_ldem = "Support liberal democracy",
  inst_trust = "Institutional trust",
  pol_tolerance = "Political tolerance",
  emancip_val = "Emancipative values",
  autonomy_values = "Authoritarian values (rev)"
)

label_reg <- c(
  demaut_pref = "Dem-aut. pref.",
  sup_ldem = "Sup. liberal democ.",
  inst_trust = "Instit. trust",
  pol_tolerance = "Pol. tolerance",
  emancip_val = "Emancip. values",
  autonomy_values = "Auth. values (rev)"
)

sample_labels <- c(
  pl = "Poland",
  sk = "South Korea",
  us = "United States"
)

scenario_labels <- c(
  overall = "Overall",
  cross_party = "Cross-party",
  cross_ideol_redist = "Cross-ideology (redist)",
  cross_ideol_gender = "Cross-ideology (gender)"
)

scenario_shapes <- c(
  "Overall" = 16,
  "Cross-party" = 15,
  "Cross-ideology (redist)" = 17,
  "Cross-ideology (gender)" = 25
)

scenario_colors <- c(
  "Overall" = "#000000",
  "Cross-party" = "goldenrod",
  "Cross-ideology (redist)" = "dodgerblue3",
  "Cross-ideology (gender)" = "darkorchid3"
)

# Containers for analysis datasets
choice_overall_list <- vector("list", length(samples))
names(choice_overall_list) <- samples
choice_party_list <- vector("list", length(samples))
names(choice_party_list) <- samples
choice_ideol_redist_list <- vector("list", length(samples))
names(choice_ideol_redist_list) <- samples
choice_ideol_gender_list <- vector("list", length(samples))
names(choice_ideol_gender_list) <- samples

#### Loop 1: build analysis datasets for each sample and contrast ----

for (s in samples) {
  cat("Building analysis datasets for", s, "...\n")

  base_df <- wide_list[[s]]
  dat_long <- long_list[[s]]

  if (!("is_cand_dem" %in% names(dat_long))) {
    dat_long <- dat_long %>%
      mutate(is_cand_dem = as.integer(cand_system_change %in% pro_dem_labels))
  }
  dat_long <- dat_long %>%
    mutate(chosen = as.integer(chosen))

  ## (1) Overall
  task_core <- dat_long %>%
    mutate(pro_dem = cand_system_change %in% pro_dem_labels) %>%
    group_by(RespondentID, task) %>%
    summarise(choose_dem_cand = as.integer(any(chosen == 1L & pro_dem)),
              .groups = "drop")

  non_cj_cols <- names(base_df)[!grepl("^CJ_Task_", names(base_df))]

  choice_overall <- task_core %>%
    left_join(base_df %>% select(all_of(non_cj_cols)), by = "RespondentID") %>%
    relocate(country, RespondentID, task, choose_dem_cand)

  choice_overall_list[[s]] <- choice_overall

  ## (2) Cross-party
  eval_cols <- c("eval_dem_us","eval_rep_us",
                 "eval_dpk_kr","eval_ppp_kr",
                 "eval_pis_pl","eval_kon_pl","eval_ko_pl","eval_lew_pl")
  eval_cols <- intersect(eval_cols, names(base_df))

  long2 <- dat_long %>%
    left_join(base_df %>% select(RespondentID, all_of(eval_cols)), by = "RespondentID") %>%
    mutate(
      party_dem_us = country == "US" & cand_party == "Democrat",
      party_rep_us = country == "US" & cand_party == "Republican",
      party_dpk_kr = country == "KR" & cand_party == "DPK",
      party_ppp_kr = country == "KR" & cand_party == "PPP",
      party_pis_pl = country == "PL" & cand_party == "PiS",
      party_kon_pl = country == "PL" & cand_party == "Konfederacja",
      party_ko_pl  = country == "PL" & cand_party == "KO",
      party_lew_pl = country == "PL" & cand_party == "Lewica"
    ) %>%
    mutate(
      in_party = (
        (party_dem_us & eval_dem_us >= 6) | (party_rep_us & eval_rep_us >= 6) |
          (party_dpk_kr & eval_dpk_kr >= 6) | (party_ppp_kr & eval_ppp_kr >= 6) |
          (party_pis_pl & eval_pis_pl >= 6) | (party_kon_pl & eval_kon_pl >= 6) |
          (party_ko_pl  & eval_ko_pl  >= 6) | (party_lew_pl & eval_lew_pl >= 6)
      ),
      out_party = (
        (party_dem_us & eval_dem_us <= 4) | (party_rep_us & eval_rep_us <= 4) |
          (party_dpk_kr & eval_dpk_kr <= 4) | (party_ppp_kr & eval_ppp_kr <= 4) |
          (party_pis_pl & eval_pis_pl <= 4) | (party_kon_pl & eval_kon_pl <= 4) |
          (party_ko_pl  & eval_ko_pl  <= 4) | (party_lew_pl & eval_lew_pl <= 4)
      )
    )

  keep_resp <- base_df %>%
    mutate(
      dem_in_us  = if ("eval_dem_us" %in% names(.))  eval_dem_us >= 6 else NA,
      rep_in_us  = if ("eval_rep_us" %in% names(.))  eval_rep_us >= 6 else NA,
      dem_out_us = if ("eval_dem_us" %in% names(.))  eval_dem_us <= 4 else NA,
      rep_out_us = if ("eval_rep_us" %in% names(.))  eval_rep_us <= 4 else NA,
      dpk_in_kr  = if ("eval_dpk_kr" %in% names(.))  eval_dpk_kr >= 6 else NA,
      ppp_in_kr  = if ("eval_ppp_kr" %in% names(.))  eval_ppp_kr >= 6 else NA,
      dpk_out_kr = if ("eval_dpk_kr" %in% names(.))  eval_dpk_kr <= 4 else NA,
      ppp_out_kr = if ("eval_ppp_kr" %in% names(.))  eval_ppp_kr <= 4 else NA,
      pis_in_pl  = if ("eval_pis_pl" %in% names(.))  eval_pis_pl >= 6 else NA,
      kon_in_pl  = if ("eval_kon_pl" %in% names(.))  eval_kon_pl >= 6 else NA,
      ko_in_pl   = if ("eval_ko_pl"  %in% names(.))  eval_ko_pl  >= 6 else NA,
      lew_in_pl  = if ("eval_lew_pl" %in% names(.))  eval_lew_pl >= 6 else NA,
      pis_out_pl = if ("eval_pis_pl" %in% names(.))  eval_pis_pl <= 4 else NA,
      kon_out_pl = if ("eval_kon_pl" %in% names(.))  eval_kon_pl <= 4 else NA,
      ko_out_pl  = if ("eval_ko_pl"  %in% names(.))  eval_ko_pl  <= 4 else NA,
      lew_out_pl = if ("eval_lew_pl" %in% names(.))  eval_lew_pl <= 4 else NA
    ) %>%
    mutate(
      drop_us = country == "US" & ((dem_in_us  & rep_in_us) | (dem_out_us & rep_out_us)),
      drop_kr = country == "KR" & ((dpk_in_kr  & ppp_in_kr) | (dpk_out_kr & ppp_out_kr)),
      all_in_pl  = rowSums(cbind(pis_in_pl, kon_in_pl, ko_in_pl, lew_in_pl), na.rm = FALSE) == 4,
      all_out_pl = rowSums(cbind(pis_out_pl, kon_out_pl, ko_out_pl, lew_out_pl), na.rm = FALSE) == 4,
      all_in_pl  = tidyr::replace_na(all_in_pl,  FALSE),
      all_out_pl = tidyr::replace_na(all_out_pl, FALSE),
      drop_pl    = country == "PL" & (all_in_pl | all_out_pl),
      drop = (drop_us | drop_kr | drop_pl)
    ) %>%
    filter(!drop) %>%
    select(RespondentID)

  long2_keep <- long2 %>%
    semi_join(keep_resp, by = "RespondentID")

  pair_flags_party <- long2_keep %>%
    group_by(RespondentID, task) %>%
    summarise(
      has_undem_in = any(is_cand_dem == 0 & in_party,  na.rm = TRUE),
      has_dem_out  = any(is_cand_dem == 1 & out_party, na.rm = TRUE),
      .groups = "drop"
    ) %>%
    mutate(keep_pair = has_undem_in & has_dem_out) %>%
    filter(keep_pair) %>%
    select(RespondentID, task)

  party_pairs <- long2_keep %>%
    inner_join(pair_flags_party, by = c("RespondentID","task"))

  party_task_core <- party_pairs %>%
    group_by(RespondentID, task) %>%
    summarise(
      choose_dem_out_vs_undem_in = as.integer(any(chosen == 1L & is_cand_dem == 1L & out_party)),
      .groups = "drop"
    )

  choice_party <- party_task_core %>%
    left_join(base_df %>% select(all_of(non_cj_cols)), by = "RespondentID") %>%
    relocate(country, RespondentID, task, choose_dem_out_vs_undem_in)

  choice_party_list[[s]] <- choice_party

  ## (3) Cross-ideology
  ideonames <- c("redist_1","redist_2","gend_pol_1","gend_pol_2")
  missing_ideo <- setdiff(ideonames, names(dat_long))
  long3 <- dat_long
  if (length(missing_ideo) > 0) {
    long3 <- long3 %>%
      left_join(base_df %>% select(RespondentID, all_of(missing_ideo)), by = "RespondentID")
  }

  long3 <- long3 %>%
    mutate(
      cand_redist_1 = case_when(
        cand_tax_spend == "spend more" ~ 4.5,
        cand_tax_spend == "spend less" ~ 1.5,
        TRUE ~ NA_real_
      ),
      cand_redist_2 = case_when(
        cand_tax_spend == "tax rich more" ~ 4.5,
        cand_tax_spend == "tax rich less" ~ 1.5,
        TRUE ~ NA_real_
      ),
      cand_gend_1 = case_when(
        cand_gend_policy == "require hire pref" ~ 4.5,
        cand_gend_policy == "no hire pref"      ~ 1.5,
        TRUE ~ NA_real_
      ),
      cand_gend_2 = case_when(
        cand_gend_policy == "require quotas" ~ 4.5,
        cand_gend_policy == "no quotas"      ~ 1.5,
        TRUE ~ NA_real_
      ),
      dist_r1 = abs(redist_1   - cand_redist_1),
      dist_r2 = abs(redist_2   - cand_redist_2),
      dist_g1 = abs(gend_pol_1 - cand_gend_1),
      dist_g2 = abs(gend_pol_2 - cand_gend_2)
    ) %>%
    rowwise() %>%
    mutate(
      prox_r = any(c(dist_r1, dist_r2) == 0.5, na.rm = TRUE),
      any_r_far = any(c(dist_r1, dist_r2) %in% c(1.5, 2.5, 3.5), na.rm = TRUE),
      dist_r = (!prox_r) & isTRUE(any_r_far),
      prox_g = any(c(dist_g1, dist_g2) == 0.5, na.rm = TRUE),
      any_g_far = any(c(dist_g1, dist_g2) %in% c(1.5, 2.5, 3.5), na.rm = TRUE),
      dist_g = (!prox_g) & isTRUE(any_g_far)
    ) %>%
    ungroup()

  pair_flags_r <- long3 %>%
    group_by(RespondentID, task) %>%
    summarise(
      has_undem_prox_r = any(is_cand_dem == 0 & prox_r, na.rm = TRUE),
      has_dem_dist_r   = any(is_cand_dem == 1 & dist_r, na.rm = TRUE),
      .groups = "drop"
    ) %>%
    mutate(keep_pair_r = has_undem_prox_r & has_dem_dist_r) %>%
    filter(keep_pair_r) %>%
    select(RespondentID, task)

  ideol_pairs_redist <- long3 %>%
    inner_join(pair_flags_r, by = c("RespondentID","task"))

  ideol_task_core_redist <- ideol_pairs_redist %>%
    group_by(RespondentID, task) %>%
    summarise(
      choose_dem_dist_vs_undem_prox_redist = as.integer(any(chosen == 1L & is_cand_dem == 1L & dist_r)),
      .groups = "drop"
    )

  choice_ideol_redist <- ideol_task_core_redist %>%
    left_join(base_df %>% select(all_of(non_cj_cols)), by = "RespondentID") %>%
    relocate(country, RespondentID, task, choose_dem_dist_vs_undem_prox_redist)

  choice_ideol_redist_list[[s]] <- choice_ideol_redist

  pair_flags_g <- long3 %>%
    group_by(RespondentID, task) %>%
    summarise(
      has_undem_prox_g = any(is_cand_dem == 0 & prox_g, na.rm = TRUE),
      has_dem_dist_g   = any(is_cand_dem == 1 & dist_g, na.rm = TRUE),
      .groups = "drop"
    ) %>%
    mutate(keep_pair_g = has_undem_prox_g & has_dem_dist_g) %>%
    filter(keep_pair_g) %>%
    select(RespondentID, task)

  ideol_pairs_gender <- long3 %>%
    inner_join(pair_flags_g, by = c("RespondentID","task"))

  ideol_task_core_gender <- ideol_pairs_gender %>%
    group_by(RespondentID, task) %>%
    summarise(
      choose_dem_dist_vs_undem_prox_gender = as.integer(any(chosen == 1L & is_cand_dem == 1L & dist_g)),
      .groups = "drop"
    )

  choice_ideol_gender <- ideol_task_core_gender %>%
    left_join(base_df %>% select(all_of(non_cj_cols)), by = "RespondentID") %>%
    relocate(country, RespondentID, task, choose_dem_dist_vs_undem_prox_gender)

  choice_ideol_gender_list[[s]] <- choice_ideol_gender
}

#### Loop 2: fit LPMs and collect results ----

res_overall_list <- vector("list", length(samples)); names(res_overall_list) <- samples
res_party_list <- vector("list", length(samples)); names(res_party_list) <- samples
res_ideol_redist_list <- vector("list", length(samples)); names(res_ideol_redist_list) <- samples
res_ideol_gender_list <- vector("list", length(samples)); names(res_ideol_gender_list) <- samples

spec_levels <- c("main cntrls","all cntrls")

for (s in samples) {
  cat("Estimating LPMs for", s, "...\n")

  covars_s <- switch(
    s,
    "pool" = c(covars_base, "top_party"),
    "us"   = c(covars_base, "eval_dem_us", "eval_rep_us"),
    "pl"   = c(covars_base, "eval_pis_pl", "eval_kon_pl", "eval_ko_pl", "eval_lew_pl"),
    "sk"   = c(covars_base, "eval_dpk_kr", "eval_ppp_kr")
  )

  fit_one_set <- function(dat_choice, dep_var) {
    set2 <- NULL
    for (v in pc_vars) {
      rhs <- paste(c(v, covars_s), collapse = " + ")
      f   <- as.formula(paste(dep_var, "~", rhs))
      req <- unique(c(all.vars(f), "RespondentID"))
      rows <- complete.cases(dat_choice[, req, drop = FALSE])
      m  <- lm(f, data = dat_choice[rows, , drop = FALSE])
      V  <- vcovCL(m, cluster = dat_choice$RespondentID[rows])
      ct <- coeftest(m, vcov = V)
      est <- if (v %in% rownames(ct)) unname(ct[v,1]) else NA_real_
      se  <- if (v %in% rownames(ct)) unname(ct[v,2]) else NA_real_
      set2 <- bind_rows(set2, tibble(term = v, estimate = est, se = se, spec = "main cntrls"))
    }

    rhs_all <- paste(c(pc_vars, covars_s), collapse = " + ")
    f_all   <- as.formula(paste(dep_var, "~", rhs_all))
    req_all <- unique(c(all.vars(f_all), "RespondentID"))
    rows_all <- complete.cases(dat_choice[, req_all, drop = FALSE])
    m_all   <- lm(f_all, data = dat_choice[rows_all, , drop = FALSE])
    V_all   <- vcovCL(m_all, cluster = dat_choice$RespondentID[rows_all])
    ct_all  <- coeftest(m_all, vcov = V_all)
    set3 <- NULL
    for (v in pc_vars) {
      est <- if (v %in% rownames(ct_all)) unname(ct_all[v,1]) else NA_real_
      se  <- if (v %in% rownames(ct_all)) unname(ct_all[v,2]) else NA_real_
      set3 <- bind_rows(set3, tibble(term = v, estimate = est, se = se, spec = "all cntrls"))
    }

    bind_rows(set2, set3) %>%
      mutate(
        var  = factor(recode(term, !!!label_map), levels = rev(unname(label_map[pc_vars]))),
        lo   = estimate - 1.96*se,
        hi   = estimate + 1.96*se,
        spec = factor(spec, levels = spec_levels)
      )
  }

  res_overall_list[[s]] <- fit_one_set(choice_overall_list[[s]], "choose_dem_cand")
  res_party_list[[s]] <- fit_one_set(choice_party_list[[s]], "choose_dem_out_vs_undem_in")
  res_ideol_redist_list[[s]] <- fit_one_set(choice_ideol_redist_list[[s]], "choose_dem_dist_vs_undem_prox_redist")
  res_ideol_gender_list[[s]] <- fit_one_set(choice_ideol_gender_list[[s]], "choose_dem_dist_vs_undem_prox_gender")
}

#### Regression tables ----

sample_titles <- c(
  pool = "Pooled samples",
  us   = "United States",
  sk   = "South Korea",
  pl   = "Poland"
)

comparison_specs <- list(
  overall = list(
    res_list      = res_overall_list,
    choice_list   = choice_overall_list,
    dep_var       = "choose_dem_cand",
    file_stub     = "lpm_overall",
    caption_stub  = "LPM: Choosing the pro-democratic candidate"
  ),
  cross_party = list(
    res_list      = res_party_list,
    choice_list   = choice_party_list,
    dep_var       = "choose_dem_out_vs_undem_in",
    file_stub     = "lpm_crossparty",
    caption_stub  = "LPM: Choosing pro-democratic out-party vs anti-democratic in-party candidate"
  ),
  cross_ideol_redist = list(
    res_list      = res_ideol_redist_list,
    choice_list   = choice_ideol_redist_list,
    dep_var       = "choose_dem_dist_vs_undem_prox_redist",
    file_stub     = "lpm_cross_ideol_redist",
    caption_stub  = "LPM: Choosing pro-democratic ideologically distal (redistribution) vs anti-democratic proximate candidate"
  ),
  cross_ideol_gender = list(
    res_list      = res_ideol_gender_list,
    choice_list   = choice_ideol_gender_list,
    dep_var       = "choose_dem_dist_vs_undem_prox_gender",
    file_stub     = "lpm_cross_ideol_gender",
    caption_stub  = "LPM: Choosing pro-democratic ideologically distal (gender) vs anti-democratic proximate candidate"
  )
)

get_covars_s <- function(s) {
  switch(
    s,
    "pool" = c(covars_base, "top_party"),
    "us"   = c(covars_base, "eval_dem_us", "eval_rep_us"),
    "pl"   = c(covars_base, "eval_pis_pl", "eval_kon_pl", "eval_ko_pl", "eval_lew_pl"),
    "sk"   = c(covars_base, "eval_dpk_kr", "eval_ppp_kr")
  )
}

for (s in samples) {
  cat("Building regression tables for sample", s, "...\n")

  covars_s <- get_covars_s(s)

  for (comp_name in names(comparison_specs)) {
    spec <- comparison_specs[[comp_name]]
    res_s   <- spec$res_list[[s]]
    dchoice <- spec$choice_list[[s]]
    if (is.null(res_s) || is.null(dchoice)) next

    tab_cov <- res_s %>%
      dplyr::filter(spec == "main cntrls", term %in% pc_vars) %>%
      dplyr::select(term, estimate, se) %>%
      dplyr::slice(match(pc_vars, term))

    tab_all <- res_s %>%
      dplyr::filter(spec == "all cntrls", term %in% pc_vars) %>%
      dplyr::select(term, estimate, se) %>%
      dplyr::slice(match(pc_vars, term))

    n_choices <- list(); n_resp <- list()
    for (v in pc_vars) {
      needed <- c(spec$dep_var, "RespondentID", v, covars_s)
      needed <- intersect(needed, names(dchoice))
      rows   <- complete.cases(dchoice[, needed, drop = FALSE])
      n_choices[[v]] <- sum(rows)
      n_resp[[v]]    <- dplyr::n_distinct(dchoice$RespondentID[rows])
    }

    needed_all <- c(spec$dep_var, "RespondentID", pc_vars, covars_s)
    needed_all <- intersect(needed_all, names(dchoice))
    rows_all   <- complete.cases(dchoice[, needed_all, drop = FALSE])
    n_choices_all <- sum(rows_all)
    n_resp_all    <- dplyr::n_distinct(dchoice$RespondentID[rows_all])

    mods <- list()
    for (i in seq_along(pc_vars)) {
      v <- pc_vars[i]
      co <- rep(0, length(pc_vars)); names(co) <- pc_vars
      se <- rep(0, length(pc_vars)); names(se) <- pc_vars
      pv <- rep(NA_real_, length(pc_vars)); names(pv) <- pc_vars
      est_v <- tab_cov$estimate[tab_cov$term == v]
      se_v  <- tab_cov$se[tab_cov$term == v]
      co[v] <- est_v
      se[v] <- se_v
      if (!is.na(est_v) && !is.na(se_v) && se_v > 0) {
        pv[v] <- 2 * pnorm(abs(est_v / se_v), lower.tail = FALSE)
      }
      mods[[length(mods) + 1L]] <- createTexreg(
        coef.names = pc_vars,
        coef       = co,
        se         = se,
        pvalues    = pv,
        gof.names  = c("N (choices)", "N (respondents)"),
        gof        = c(n_choices[[v]], n_resp[[v]]),
        gof.decimal = c(FALSE, FALSE)
      )
    }

    co_all <- tab_all$estimate; names(co_all) <- pc_vars
    se_all <- tab_all$se;       names(se_all) <- pc_vars
    pv_all <- ifelse(!is.na(co_all) & !is.na(se_all) & se_all > 0,
                     2 * pnorm(abs(co_all / se_all), lower.tail = FALSE),
                     NA_real_)

    mods[[length(mods) + 1L]] <- createTexreg(
      coef.names = pc_vars,
      coef       = co_all,
      se         = se_all,
      pvalues    = pv_all,
      gof.names  = c("N (choices)", "N (respondents)"),
      gof        = c(n_choices_all, n_resp_all),
      gof.decimal = c(FALSE, FALSE)
    )

  }
}

#### Faceted LPM plots ----

build_lpm_facet_plot <- function(spec_to_plot = c("all cntrls", "main cntrls")) {
  spec_to_plot <- match.arg(spec_to_plot)

  plot_df <- bind_rows(
    bind_rows(lapply(plot_samples, function(s) {
      res_overall_list[[s]] %>% filter(spec == spec_to_plot) %>% mutate(sample = s, scenario = "Overall")
    })),
    bind_rows(lapply(plot_samples, function(s) {
      res_party_list[[s]] %>% filter(spec == spec_to_plot) %>% mutate(sample = s, scenario = "Cross-party")
    })),
    bind_rows(lapply(plot_samples, function(s) {
      res_ideol_redist_list[[s]] %>% filter(spec == spec_to_plot) %>% mutate(sample = s, scenario = "Cross-ideology (redist)")
    })),
    bind_rows(lapply(plot_samples, function(s) {
      res_ideol_gender_list[[s]] %>% filter(spec == spec_to_plot) %>% mutate(sample = s, scenario = "Cross-ideology (gender)")
    }))
  ) %>%
    mutate(
      sample = factor(sample, levels = c("pl", "sk", "us")),
      sample_label = factor(sample_labels[as.character(sample)], levels = c("Poland", "South Korea", "United States")),
      scenario = factor(scenario,
                        levels = c("Overall", "Cross-party", "Cross-ideology (redist)", "Cross-ideology (gender)")),
      pc_label = factor(recode(term, !!!label_map), levels = unname(pc_labels[pc_vars]))
    )

  dodge_pos <- position_dodge(width = 0.55, reverse = TRUE)

  ggplot(
    plot_df,
    aes(x = estimate, y = sample_label, color = scenario, fill = scenario, shape = scenario)
  ) +
    geom_vline(xintercept = 0, linetype = "dotted", color = "gray40") +
    geom_errorbar(
      aes(xmin = lo, xmax = hi),
      orientation = "y",
      width = 0.12,
      linewidth = 0.6,
      position = dodge_pos
    ) +
    geom_point(size = 2.7, stroke = 0.7, position = dodge_pos) +
    facet_wrap(~ pc_label, ncol = 2) +
    scale_color_manual(values = scenario_colors) +
    scale_fill_manual(values = scenario_colors) +
    scale_shape_manual(values = scenario_shapes) +
    scale_x_continuous(limits = c(-0.10, 0.20), breaks = seq(-0.10, 0.20, by = 0.05)) +
    labs(
      x = "LPM coefficient",
      y = NULL,
      color = NULL,
      fill = NULL
    ) +
    guides(
      color = guide_legend(
        override.aes = list(shape = unname(scenario_shapes), fill = unname(scenario_colors), size = 3.2),
        nrow = 1,
        byrow = TRUE
      ),
      fill = "none",
      shape = "none"
    ) +
    theme_bw(base_size = 11) +
    theme(
      panel.grid.minor = element_blank(),
      panel.grid.major.x = element_line(color = "grey85", linewidth = 0.3),
      panel.grid.major.y = element_line(color = "grey85", linewidth = 0.3),
      strip.background = element_blank(),
      strip.text = element_text(size = 10, face = "plain"),
      legend.position = "bottom",
      legend.title = element_blank(),
      legend.margin = margin(t = 2, r = 0, b = 2, l = 0),
      legend.box.spacing = grid::unit(0.35, "lines"),
      legend.key.height = grid::unit(0.32, "lines"),
      legend.text = element_text(size = 10),
      panel.spacing.x = grid::unit(0.8, "lines"),
      panel.spacing.y = grid::unit(0.8, "lines"),
      plot.margin = margin(4, 4, 10, 4)
    )
}

p_lpm_all <- build_lpm_facet_plot("all cntrls")
print(p_lpm_all)
ggsave(file.path(project_root, "outputs", "LPMs_faceted_allcovar_v2.png"), p_lpm_all, width = 8.6, height = 7.6, dpi = 250)

p_lpm_main <- build_lpm_facet_plot("main cntrls")
print(p_lpm_main)
ggsave(file.path(project_root, "outputs", "LPMs_faceted_maincovar_v2.png"), p_lpm_main, width = 8.6, height = 7.6, dpi = 250)
