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 Modes (LPMs)

# libraries
library(dplyr)
library(tidyr)
library(ggplot2)
library(lmtest)
library(sandwich)
library(texreg)
library(here)

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)

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")

label_map <- c(
  demaut_pref = "Democ-autoc.\npreference",
  sup_ldem = "Support liberal\ndemocracy",
  inst_trust = "Institutional\ntrust",
  pol_tolerance = "Political\ntolerance",
  emancip_val = "Emancipative\nvalues",
  autonomy_values = "Authoritarian\nvalues (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)"
)

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

# 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]]
  
  # Ensure key variables exist / are correctly typed
  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 (H4)
  
  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 (H5)
  
  # Bring evals into long; define party flags by country
  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")
  
  # candidate party flags (use factor labels)
  long2 <- long2 %>%
    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)
      )
    )
  
  # Drop respondents with ambiguous positioning within-country
  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)
      ),
      # Poland: only drop if ALL four parties in, or ALL four out
      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
  
  ## Cross redist ideol (H6)
  
  # Start from long; bring ideology items 
  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(
      # candidate positions (redistribution)
      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_
      ),
      # candidate positions (gender)
      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()
  
  # Redistributive pairs: anti-dem proximate vs pro-dem distal
  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
  
  ## Cross gender ideol (H6)
  
  # Gender pairs: anti-dem proximate vs pro-dem distal
  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")
  
  ## Sample-specific covariates
  covars_base <- setdiff(covars, "top_party")
  covars_s <- switch(
    s,
    "pool" = covars,  
    "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")
  )
  
  ## Overall
  
  dat_choice <- choice_overall_list[[s]]
  set2 <- NULL
  for (v in pc_vars) {
    rhs <- paste(c(v, covars_s), collapse = " + ")
    f   <- as.formula(paste("choose_dem_cand ~", 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("choose_dem_cand ~", 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"))
  }
  res_overall <- 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]] <- res_overall
  
  ## Cross-party
  
  dat_choice <- choice_party_list[[s]]
  set2 <- NULL
  for (v in pc_vars) {
    rhs <- paste(c(v, covars_s), collapse = " + ")
    f   <- as.formula(paste("choose_dem_out_vs_undem_in ~", 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("choose_dem_out_vs_undem_in ~", 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"))
  }
  res_party <- 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_party_list[[s]] <- res_party
  
  ## Redist ideology
  
  dat_choice <- choice_ideol_redist_list[[s]]
  set2 <- NULL
  for (v in pc_vars) {
    rhs <- paste(c(v, covars_s), collapse = " + ")
    f   <- as.formula(paste("choose_dem_dist_vs_undem_prox_redist ~", 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("choose_dem_dist_vs_undem_prox_redist ~", 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"))
  }
  res_ideol_r <- 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_ideol_redist_list[[s]] <- res_ideol_r
  
  ## Gender ideology
  
  dat_choice <- choice_ideol_gender_list[[s]]
  set2 <- NULL
  for (v in pc_vars) {
    rhs <- paste(c(v, covars_s), collapse = " + ")
    f   <- as.formula(paste("choose_dem_dist_vs_undem_prox_gender ~", 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("choose_dem_dist_vs_undem_prox_gender ~", 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"))
  }
  res_ideol_g <- 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_ideol_gender_list[[s]] <- res_ideol_g
}


#### Regression tables----

# Base covariates (no sample-specific items here)
covars_base <- c(
  "gender_male","age_1","educ_5f","income_cat","region_pool",
  "pol_interest","gend_pol_left","redist_pol_left"
)

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

# Map each comparison to its result list, choice data list, dep var & caption stub
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"
  )
)

# Function - sample-specific covariate set (same logic as ATMEs)
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")
  )
}

# Build 7-model LaTeX tables for all samples × comparisons
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_list    <- spec$res_list
    choice_list <- spec$choice_list
    dep_var     <- spec$dep_var
    file_stub   <- spec$file_stub
    caption_stub<- spec$caption_stub
    
    res_s   <- res_list[[s]]
    dchoice <- choice_list[[s]]

    if (is.null(res_s) || is.null(dchoice)) next
    
    ## 1. Extract coefficient estimates & SEs for:
    ##    - Models 1–6: spec == "main cntrls", each PC separately
    ##    - Model 7:    spec == "all cntrls", all PCs together
    
    # rows for individual PC models with main controls
    tab_cov <- res_s %>%
      dplyr::filter(spec == "main cntrls", term %in% pc_vars) %>%
      dplyr::select(term, estimate, se) %>%
      dplyr::slice(match(pc_vars, term))
    
    # rows for the combined PC model (all cntrls spec)
    tab_all <- res_s %>%
      dplyr::filter(spec == "all cntrls", term %in% pc_vars) %>%
      dplyr::select(term, estimate, se) %>%
      dplyr::slice(match(pc_vars, term))
    
    ## 2. Compute N(choices) & N(respondents) per model,
    ##    respecting sample-specific covariates.
    
    n_choices <- list()
    n_resp    <- list()
    
    # Models 1–6: each PC variable separately + covars_s
    for (v in pc_vars) {
      needed <- c(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])
    }
    
    # Model 7: all PC variables together + covars_s
    needed_all <- c(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])
    
    ## 3. Build texreg model objects (7 models)
    
    mods <- list()
    
    # Models 1–6: each PC variable in its own column
    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)
      )
    }
    
    # Model 7: all PC variables together (from tab_all)
    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)
    )
    
  } # end comparison loop
}   # end sample loop



#### Combined LPM plots----

comp_labels  <- c(overall = "Overall", cross_party = "Cross-party",
                  cross_ideol_redist = "Cross-ideology (redist)",
                  cross_ideol_gender = "Cross-ideology (gend)")
comp_pal <- c(overall = "#000000",
              cross_party = "goldenrod",
              cross_ideol_redist = "dodgerblue3",
              cross_ideol_gender = "dodgerblue3")
comp_shapes  <- c(overall = 21,
                  cross_party = 22,
                  cross_ideol_redist = 24,
                  cross_ideol_gender = 25)
y_offset <- c(overall = +0.24,
              cross_party = +0.08,
              cross_ideol_redist = -0.08,
              cross_ideol_gender = -0.24)

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

for (s in samples) {
  cat("Plotting combined direct LPM effects for", s, "...\n")
  
  res_overall <- res_overall_list[[s]]
  res_party      <- res_party_list[[s]]
  res_ideol_r    <- res_ideol_redist_list[[s]]
  res_ideol_g    <- res_ideol_gender_list[[s]]
  
  ## All controls
  
  block_overall_all <- res_overall %>%
    filter(spec == "all cntrls", term %in% pc_vars) %>%
    mutate(
      var   = factor(recode(term, !!!label_map), levels = rev(unname(label_map[pc_vars]))),
      comp  = "overall",
      y_base = as.numeric(var),
      y     = y_base + y_offset[comp]
    ) %>%
    select(term, var, estimate, se, lo, hi, comp, y_base, y)
  
  block_party_all <- res_party %>%
    filter(spec == "all cntrls", term %in% pc_vars) %>%
    mutate(
      var   = factor(recode(term, !!!label_map), levels = rev(unname(label_map[pc_vars]))),
      comp  = "cross_party",
      y_base = as.numeric(var),
      y     = y_base + y_offset[comp]
    ) %>%
    select(term, var, estimate, se, lo, hi, comp, y_base, y)
  
  block_ideolR_all <- res_ideol_r %>%
    filter(spec == "all cntrls", term %in% pc_vars) %>%
    mutate(
      var   = factor(recode(term, !!!label_map), levels = rev(unname(label_map[pc_vars]))),
      comp  = "cross_ideol_redist",
      y_base = as.numeric(var),
      y     = y_base + y_offset[comp]
    ) %>%
    select(term, var, estimate, se, lo, hi, comp, y_base, y)
  
  block_ideolG_all <- res_ideol_g %>%
    filter(spec == "all cntrls", term %in% pc_vars) %>%
    mutate(
      var   = factor(recode(term, !!!label_map), levels = rev(unname(label_map[pc_vars]))),
      comp  = "cross_ideol_gender",
      y_base = as.numeric(var),
      y     = y_base + y_offset[comp]
    ) %>%
    select(term, var, estimate, se, lo, hi, comp, y_base, y)
  
  combo_all <- bind_rows(block_overall_all, block_party_all,
                         block_ideolR_all, block_ideolG_all)
  
  combo_all$var <- factor(combo_all$var, levels = rev(label_map[pc_order]))
  
  p_all <- ggplot(combo_all, aes(x = estimate, y = y, color = comp, fill = comp, shape = comp)) +
    geom_vline(xintercept = 0, linewidth = 0.4, linetype = 3) +
    geom_errorbar(aes(xmin = lo, xmax = hi), orientation = "y", width = 0, linewidth = 1) +
    geom_point(size = 2.5, stroke = 0.6) +
    scale_y_continuous(
      breaks = sort(unique(combo_all$y_base)),
      labels = levels(combo_all$var)
    ) +
    scale_shape_manual(values = comp_shapes, breaks = names(comp_labels),
                       labels = comp_labels, name = NULL,
                       guide = guide_legend(nrow = 1, byrow = TRUE,
                                            override.aes = list(stroke = 0))) +
    scale_fill_manual(values = comp_pal, breaks = names(comp_labels),
                      labels = comp_labels, name = NULL) +
    scale_color_manual(values = comp_pal, breaks = names(comp_labels),
                       labels = comp_labels, name = NULL) +
    labs(
      title = sample_titles[s],
      x = "LPM effect of a 1 SD increase in culture variable",
      y = NULL
    ) +
    theme_bw(base_size = 14) +
    theme(
      panel.grid.minor  = element_blank(),
      legend.position   = "bottom",
      legend.box.spacing= grid::unit(0.1, "lines"),
      legend.key.height = grid::unit(0.4, "lines"),
      legend.text       = element_text(size = 11)
    )
  
  print(p_all)
  ggsave(file.path(project_root, "outputs", paste0("direct_lpm_combined_allcovar_", s, "_v2.png")),
         p_all, width = 7, height = 5, dpi = 250)
  
  ## Main controls only
  
  block_overall_plus <- res_overall %>%
    filter(spec == "main cntrls", term %in% pc_vars) %>%
    mutate(
      var   = factor(recode(term, !!!label_map), levels = rev(unname(label_map[pc_vars]))),
      comp  = "overall",
      y_base = as.numeric(var),
      y     = y_base + y_offset[comp]
    ) %>%
    select(term, var, estimate, se, lo, hi, comp, y_base, y)
  
  block_party_plus <- res_party %>%
    filter(spec == "main cntrls", term %in% pc_vars) %>%
    mutate(
      var   = factor(recode(term, !!!label_map), levels = rev(unname(label_map[pc_vars]))),
      comp  = "cross_party",
      y_base = as.numeric(var),
      y     = y_base + y_offset[comp]
    ) %>%
    select(term, var, estimate, se, lo, hi, comp, y_base, y)
  
  block_ideolR_plus <- res_ideol_r %>%
    filter(spec == "main cntrls", term %in% pc_vars) %>%
    mutate(
      var   = factor(recode(term, !!!label_map), levels = rev(unname(label_map[pc_vars]))),
      comp  = "cross_ideol_redist",
      y_base = as.numeric(var),
      y     = y_base + y_offset[comp]
    ) %>%
    select(term, var, estimate, se, lo, hi, comp, y_base, y)
  
  block_ideolG_plus <- res_ideol_g %>%
    filter(spec == "main cntrls", term %in% pc_vars) %>%
    mutate(
      var   = factor(recode(term, !!!label_map), levels = rev(unname(label_map[pc_vars]))),
      comp  = "cross_ideol_gender",
      y_base = as.numeric(var),
      y     = y_base + y_offset[comp]
    ) %>%
    select(term, var, estimate, se, lo, hi, comp, y_base, y)
  
  combo_plus <- bind_rows(block_overall_plus, block_party_plus,
                          block_ideolR_plus, block_ideolG_plus)
  
  combo_plus$var <- factor(combo_plus$var, levels = rev(label_map[pc_order]))
  
  p_plus <- ggplot(combo_plus, aes(x = estimate, y = y, color = comp, fill = comp, shape = comp)) +
    geom_vline(xintercept = 0, linewidth = 0.4, linetype = 3) +
    geom_errorbar(aes(xmin = lo, xmax = hi), orientation = "y", width = 0, linewidth = 1) +
    geom_point(size = 2.5, stroke = 0.6) +
    scale_y_continuous(
      breaks = sort(unique(combo_plus$y_base)),
      labels = levels(combo_plus$var)
    ) +
    scale_shape_manual(values = comp_shapes, breaks = names(comp_labels),
                       labels = comp_labels, name = NULL,
                       guide = guide_legend(nrow = 1, byrow = TRUE,
                                            override.aes = list(stroke = 0))) +
    scale_fill_manual(values = comp_pal, breaks = names(comp_labels),
                      labels = comp_labels, name = NULL) +
    scale_color_manual(values = comp_pal, breaks = names(comp_labels),
                       labels = comp_labels, name = NULL) +
    labs(
      title = sample_titles[s],
      x = "LPM effect of a 1 SD increase in culture variable",
      y = NULL
    ) +
    theme_bw(base_size = 14) +
    theme(
      panel.grid.minor  = element_blank(),
      legend.position   = "bottom",
      legend.box.spacing= grid::unit(0.1, "lines"),
      legend.key.height = grid::unit(0.4, "lines"),
      legend.text       = element_text(size = 11)
    )
  
  print(p_plus)
  ggsave(file.path(project_root, "outputs", paste0("direct_lpm_combined_maincovar_", s, "_v2.png")),
         p_plus, width = 7, height = 5, dpi = 250)
}
