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()))
}
#### Code to analyse additional R&R ATME, LPM, and descriptive-slope requests.
#### Builds data from raw files and saves exact-p/FDR and slope outputs.
### Packages and file paths----

local_lib <- file.path(getwd(), "r-lib")
.libPaths(c(local_lib, .libPaths()))

library(dplyr)
library(tidyr)
library(tibble)
library(ggplot2)
library(sandwich)
library(lmtest)
library(texreg)
library(psych)
library(lavaan)
library(here)

set.seed(2025)

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

### Load raw data and code measures----

us <- read.csv(file.path(project_root, "data", "democ_cultures_us.csv"))
sk <- read.csv(file.path(project_root, "data", "democ_cultures_sk.csv"))
pl <- read.csv(file.path(project_root, "data", "democ_cultures_pl.csv"))

raw_n <- data.frame(
  country = c("United States", "South Korea", "Poland"),
  raw_n = c(nrow(us), nrow(sk), nrow(pl))
)

rt_med <- median(pl$Response_duration_seconds, na.rm = TRUE)
rt_cut <- 0.40 * rt_med
drop_ids <- pl$RespondentID[pl$Response_duration_seconds < rt_cut]
pl <- dplyr::filter(pl, !RespondentID %in% drop_ids)

rt_med <- median(sk$Response_duration_seconds, na.rm = TRUE)
rt_cut <- 0.40 * rt_med
drop_ids <- sk$RespondentID[sk$Response_duration_seconds < rt_cut]
sk <- dplyr::filter(sk, !RespondentID %in% drop_ids)

rt_med <- median(us$Response_duration_seconds, na.rm = TRUE)
rt_cut <- 0.40 * rt_med
drop_ids <- us$RespondentID[us$Response_duration_seconds < rt_cut]
us <- dplyr::filter(us, !RespondentID %in% drop_ids)

us$country <- "US"
sk$country <- "KR"
pl$country <- "PL"

us$gender_male <- ifelse(us$gender == 1, 1, 0)
sk$gender_male <- ifelse(sk$gender == 1, 1, 0)
pl$gender_male <- ifelse(pl$gender == 1, 1, 0)

us <- dplyr::mutate(us, educ_5f = dplyr::case_when(
  education %in% 1:4 ~ "less high school",
  education == 5 ~ "high school grad",
  education == 6 ~ "some college/2yr",
  education == 7 ~ "bachelors",
  education == 8 ~ "postgrad",
  TRUE ~ NA_character_
))
sk <- dplyr::mutate(sk, educ_5f = dplyr::case_when(
  education %in% 1:4 ~ "less high school",
  education == 5 ~ "high school grad",
  education == 6 ~ "some college/2yr",
  education == 7 ~ "bachelors",
  education == 8 ~ "postgrad",
  TRUE ~ NA_character_
))
pl <- dplyr::mutate(pl, educ_5f = dplyr::case_when(
  education %in% 1:4 ~ "less high school",
  education == 5 ~ "high school grad",
  education == 6 ~ "some college/2yr",
  education == 7 ~ "bachelors",
  education == 8 ~ "postgrad",
  TRUE ~ NA_character_
))

us <- dplyr::mutate(us, income_cat = dplyr::case_when(
  hh_income == 1 ~ "Lowest",
  hh_income == 2 ~ "Low",
  hh_income == 3 ~ "Middle",
  hh_income == 4 ~ "High",
  hh_income == 5 ~ "Highest",
  hh_income == 98 ~ "Refused",
  TRUE ~ NA_character_
))
pl <- dplyr::mutate(pl, income_cat = dplyr::case_when(
  hh_income == 6 ~ "Lowest",
  hh_income == 7 ~ "Low",
  hh_income == 8 ~ "Middle",
  hh_income == 9 ~ "High",
  hh_income == 10 ~ "Highest",
  hh_income == 98 ~ "Refused",
  TRUE ~ NA_character_
))
sk <- dplyr::mutate(sk, income_cat = dplyr::case_when(
  hh_income == 11 ~ "Lowest",
  hh_income == 12 ~ "Low",
  hh_income == 13 ~ "Middle",
  hh_income == 14 ~ "High",
  hh_income == 15 ~ "Highest",
  hh_income == 98 ~ "Refused",
  TRUE ~ NA_character_
))

us$region_pool <- us$cens_division
sk$region_pool <- sk$region_kr
pl$region_pool <- pl$region_pl

lvl_income <- c("Lowest", "Low", "Middle", "High", "Highest", "Refused")
lvl_educ <- c("less high school", "high school grad", "some college/2yr",
              "bachelors", "postgrad")
for (nm in c("us", "sk", "pl")) {
  df <- get(nm)
  df$income_cat <- factor(df$income_cat, levels = lvl_income)
  df$educ_5f <- factor(df$educ_5f, levels = lvl_educ)
  assign(nm, df)
}

for (nm in c("pl", "sk", "us")) {
  df <- get(nm)
  df$demaut_pref <- scale(df$sup_dem1 + df$sup_dem3 - (df$sup_dem2 + df$sup_dem4))[, 1]
  df$sup_ldem <- scale(rowSums(df[, c("sup_ld_1", "sup_ld_2", "sup_ld_3",
                                      "sup_ld_4", "sup_ld_5", "sup_ld_6",
                                      "sup_ld_7")], na.rm = FALSE))[, 1] * -1
  df$sup_ldem_loo_media <- scale(rowSums(df[, c("sup_ld_2", "sup_ld_3",
                                                "sup_ld_4", "sup_ld_5",
                                                "sup_ld_6", "sup_ld_7")],
                                         na.rm = FALSE))[, 1] * -1
  df$sup_ldem_loo_vote <- scale(rowSums(df[, c("sup_ld_1", "sup_ld_2",
                                               "sup_ld_4", "sup_ld_5",
                                               "sup_ld_6", "sup_ld_7")],
                                        na.rm = FALSE))[, 1] * -1
  df$sup_ldem_loo_jud <- scale(rowSums(df[, c("sup_ld_1", "sup_ld_2",
                                              "sup_ld_3", "sup_ld_4",
                                              "sup_ld_6", "sup_ld_7")],
                                       na.rm = FALSE))[, 1] * -1
  df$sup_ldem_no_ov <- scale(rowSums(df[, c("sup_ld_2", "sup_ld_4",
                                            "sup_ld_6", "sup_ld_7")],
                                        na.rm = FALSE))[, 1] * -1
  df$inst_trust <- scale(rowSums(df[, c("trust_leg", "trust_parl",
                                        "trust_civ", "trust_part")],
                                    na.rm = FALSE))[, 1]
  df$pol_tolerance <- scale(rowSums(df[, c("pol_tol_1", "pol_tol_2",
                                           "pol_tol_3")], na.rm = FALSE))[, 1] * -1

  val_names <- c("val_conf_1", "val_trad_1", "val_selfd_1", "val_stim_1",
                 "val_conf_2", "val_trad_2", "val_selfd_2", "val_stim_2")
  res_mean <- rowMeans(df[, val_names], na.rm = TRUE)
  ips <- sapply(df[, val_names], function(x) x - res_mean)
  colnames(ips) <- paste0(val_names, "_ips")
  for (ips_nm in colnames(ips)) df[[ips_nm]] <- ips[, ips_nm]
  df$conformity <- rowMeans(ips[, c("val_conf_1_ips", "val_trad_1_ips",
                                    "val_conf_2_ips", "val_trad_2_ips")],
                            na.rm = TRUE)
  df$autonomy <- rowMeans(ips[, c("val_selfd_1_ips", "val_stim_1_ips",
                                  "val_selfd_2_ips", "val_stim_2_ips")],
                          na.rm = TRUE)
  df$auth_values <- scale(df$conformity - df$autonomy)[, 1]
  df$autonomy_values <- scale(df$autonomy - df$conformity)[, 1]

  df$child_no_obed <- 1 - df$child_obed
  df$autonomy_z <- as.numeric(scale(df$child_ind + df$child_imag + df$child_no_obed))
  df$equality_rev <- rowMeans(df[, c("gend_eq_1", "gend_eq_2", "gend_eq_3")],
                              na.rm = TRUE)
  df$equality_z <- as.numeric(scale(6 - df$equality_rev))
  df$choice_z <- as.numeric(scale(df$gay_just + df$abort_just + df$divorc_just))
  v1 <- dplyr::case_when(
    (df$post_mat_11 == 2 & df$post_mat_12 == 4) | (df$post_mat_11 == 4 & df$post_mat_12 == 2) ~ 3,
    (df$post_mat_11 == 2 & df$post_mat_12 != 4) | (df$post_mat_11 == 4 & df$post_mat_12 != 2) ~ 2,
    (df$post_mat_11 != 2 & df$post_mat_12 == 4) | (df$post_mat_11 != 4 & df$post_mat_12 == 2) ~ 1,
    (df$post_mat_11 != 2 & df$post_mat_11 != 4) | (df$post_mat_12 != 4 & df$post_mat_12 != 2) ~ 0,
    TRUE ~ NA_real_
  )
  v2 <- dplyr::case_when(
    (df$post_mat_21 == 3 & df$post_mat_22 == 4) | (df$post_mat_21 == 4 & df$post_mat_22 == 3) ~ 3,
    (df$post_mat_21 == 3 & df$post_mat_22 != 4) | (df$post_mat_21 == 4 & df$post_mat_22 != 3) ~ 2,
    (df$post_mat_21 != 3 & df$post_mat_22 == 4) | (df$post_mat_21 != 4 & df$post_mat_22 == 3) ~ 1,
    (df$post_mat_21 != 3 & df$post_mat_21 != 4) | (df$post_mat_22 != 4 & df$post_mat_22 != 3) ~ 0,
    TRUE ~ NA_real_
  )
  df$postmat_voice_1 <- v1
  df$postmat_voice_2 <- v2
  df$voice_z <- as.numeric(scale(v1 + v2))
  df$emancip_val <- as.numeric(scale(rowMeans(df[, c("autonomy_z", "equality_z",
                                                     "choice_z", "voice_z")],
                                             na.rm = TRUE)))
  df$gend_pol_left <- scale(df$gend_pol_1 + df$gend_pol_2)[, 1]
  df$redist_pol_left <- scale(df$redist_1 + df$redist_2)[, 1]
  assign(nm, df)
}

party_labels <- c("US_Dem", "US_Rep", "KR_DPK", "KR_PPP",
                  "PL_PiS", "PL_Kon", "PL_KO", "PL_Lew")

eval_mat_pl <- cbind(PL_PiS = pl$eval_pis_pl, PL_Kon = pl$eval_kon_pl,
                     PL_KO = pl$eval_ko_pl, PL_Lew = pl$eval_lew_pl)
labels_pl <- colnames(eval_mat_pl)
all_na_row <- apply(eval_mat_pl, 1, function(x) all(is.na(x)))
top_party_pl <- rep(NA_character_, nrow(pl))
bottom_party_pl <- rep(NA_character_, nrow(pl))
top_party_pl[!all_na_row] <- labels_pl[max.col(ifelse(is.na(eval_mat_pl), -Inf, eval_mat_pl))[!all_na_row]]
bottom_party_pl[!all_na_row] <- labels_pl[max.col(ifelse(is.na(eval_mat_pl), -Inf, -eval_mat_pl))[!all_na_row]]
pl$top_party <- factor(top_party_pl, levels = party_labels)
pl$bottom_party <- factor(bottom_party_pl, levels = party_labels)

eval_mat_sk <- cbind(KR_DPK = sk$eval_dpk_kr, KR_PPP = sk$eval_ppp_kr)
labels_sk <- colnames(eval_mat_sk)
all_na_row <- apply(eval_mat_sk, 1, function(x) all(is.na(x)))
top_party_sk <- rep(NA_character_, nrow(sk))
bottom_party_sk <- rep(NA_character_, nrow(sk))
top_party_sk[!all_na_row] <- labels_sk[max.col(ifelse(is.na(eval_mat_sk), -Inf, eval_mat_sk))[!all_na_row]]
bottom_party_sk[!all_na_row] <- labels_sk[max.col(ifelse(is.na(eval_mat_sk), -Inf, -eval_mat_sk))[!all_na_row]]
sk$top_party <- factor(top_party_sk, levels = party_labels)
sk$bottom_party <- factor(bottom_party_sk, levels = party_labels)

eval_mat_us <- cbind(US_Dem = us$eval_dem_us, US_Rep = us$eval_rep_us)
labels_us <- colnames(eval_mat_us)
all_na_row <- apply(eval_mat_us, 1, function(x) all(is.na(x)))
top_party_us <- rep(NA_character_, nrow(us))
bottom_party_us <- rep(NA_character_, nrow(us))
top_party_us[!all_na_row] <- labels_us[max.col(ifelse(is.na(eval_mat_us), -Inf, eval_mat_us))[!all_na_row]]
bottom_party_us[!all_na_row] <- labels_us[max.col(ifelse(is.na(eval_mat_us), -Inf, -eval_mat_us))[!all_na_row]]
us$top_party <- factor(top_party_us, levels = party_labels)
us$bottom_party <- factor(bottom_party_us, levels = party_labels)

pool <- dplyr::bind_rows(us, sk, pl)
pool$country <- factor(pool$country, levels = c("US", "KR", "PL"))
pool$region_pool <- factor(pool$region_pool)
pool$top_party <- factor(pool$top_party, levels = party_labels)

write.csv(us, file.path(project_root, "data", "us_analytic_v3.csv"), row.names = FALSE)
write.csv(sk, file.path(project_root, "data", "sk_analytic_v3.csv"), row.names = FALSE)
write.csv(pl, file.path(project_root, "data", "pl_analytic_v3.csv"), row.names = FALSE)
write.csv(pool, file.path(project_root, "data", "pool_analytic_v3.csv"), row.names = FALSE)

### Reshape conjoint data and build comparison samples----

dataset_names <- c("pool", "us", "sk", "pl")
for (nm in dataset_names) {
  df <- get(nm)
  long_list <- vector("list", length = 24L)
  idx <- 1L
  for (tk in 1:12) {
    choice_var <- sprintf("CJ_Task_%d_choice", tk)
    c1_vars <- c("RespondentID", choice_var,
                 sprintf("CJ_Task_%d_Cand_1_Att%d", tk, 1:7))
    c1 <- df[, c1_vars, drop = FALSE]
    names(c1)[-(1:2)] <- paste0("Att", 1:7)
    c1$task <- tk
    c1$cand <- 1L
    c1$chosen <- as.integer(c1[[choice_var]] == 1L)
    c1 <- c1[, c("RespondentID", "task", "cand", "chosen", paste0("Att", 1:7))]
    long_list[[idx]] <- c1
    idx <- idx + 1L

    c2_vars <- c("RespondentID", choice_var,
                 sprintf("CJ_Task_%d_Cand_2_Att%d", tk, 1:7))
    c2 <- df[, c2_vars, drop = FALSE]
    names(c2)[-(1:2)] <- paste0("Att", 1:7)
    c2$task <- tk
    c2$cand <- 2L
    c2$chosen <- as.integer(c2[[choice_var]] == 2L)
    c2 <- c2[, c("RespondentID", "task", "cand", "chosen", paste0("Att", 1:7))]
    long_list[[idx]] <- c2
    idx <- idx + 1L
  }
  df_long <- dplyr::bind_rows(long_list)
  df_long <- dplyr::left_join(df_long, dplyr::select(df, RespondentID, country),
                              by = "RespondentID")
  names(df_long)[names(df_long) == "Att1"] <- "cand_age"
  names(df_long)[names(df_long) == "Att2"] <- "cand_gender"
  names(df_long)[names(df_long) == "Att3"] <- "cand_background"
  names(df_long)[names(df_long) == "Att4"] <- "cand_party"
  names(df_long)[names(df_long) == "Att5"] <- "cand_tax_spend"
  names(df_long)[names(df_long) == "Att6"] <- "cand_gend_policy"
  names(df_long)[names(df_long) == "Att7"] <- "cand_system_change"
  df_long$cand_age <- factor(df_long$cand_age, levels = 1:6,
                             labels = c("35", "44", "51", "59", "64", "68"))
  df_long$cand_gender <- factor(df_long$cand_gender, levels = 1:2,
                                labels = c("male", "female"))
  df_long$cand_background <- factor(df_long$cand_background, levels = 1:7,
                                    labels = c("lawyer", "politician", "journalist",
                                               "founder", "military", "athlete",
                                               "scientist"))
  df_long$cand_tax_spend <- factor(df_long$cand_tax_spend, levels = 1:4,
                                   labels = c("spend more", "spend less",
                                              "tax rich more", "tax rich less"))
  df_long$cand_gend_policy <- factor(df_long$cand_gend_policy, levels = 1:4,
                                     labels = c("require hire pref", "no hire pref",
                                                "require quotas", "no quotas"))
  df_long$cand_system_change <- factor(df_long$cand_system_change, levels = 1:6,
                                       labels = c("cut polling stations",
                                                  "keep polling stations",
                                                  "ignore court rulings",
                                                  "respect court rulings",
                                                  "prosecute journalists",
                                                  "don't prosecute journalists"))
  party_map <- tibble::tribble(
    ~country, ~cand_party, ~party_label,
    "US", 1L, "Democrat",
    "US", 2L, "Republican",
    "KR", 1L, "DPK",
    "KR", 2L, "PPP",
    "PL", 1L, "PiS",
    "PL", 2L, "KO",
    "PL", 3L, "Konfederacja",
    "PL", 4L, "Lewica"
  )
  df_long$cand_party <- as.integer(df_long$cand_party)
  df_long <- dplyr::left_join(df_long, party_map, by = c("country", "cand_party"))
  df_long$cand_party <- factor(df_long$party_label,
                               levels = c("Democrat", "Republican", "DPK", "PPP",
                                          "PiS", "KO", "Konfederacja", "Lewica"))
  df_long$party_label <- NULL
  assign(paste0(nm, "_long"), df_long)
  write.csv(df_long, file.path(project_root, "data", paste0(nm, "_long_v3.csv")),
            row.names = FALSE)
}

pc_vars <- c("demaut_pref", "sup_ldem", "inst_trust", "pol_tolerance",
             "emancip_val", "autonomy_values")
pc_labels <- c(
  demaut_pref = "Democracy-autocracy preference",
  sup_ldem = "Support for liberal democracy",
  inst_trust = "Institutional trust",
  pol_tolerance = "Political tolerance",
  emancip_val = "Emancipative values",
  autonomy_values = "Authoritarian values (rev)"
)
pc_plot_labels <- 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)"
)
covars <- c("gender_male", "age_1", "educ_5f", "income_cat", "region_pool",
            "pol_interest", "gend_pol_left", "redist_pol_left", "top_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")
pro_dem_labels <- c("keep polling stations", "respect court rulings",
                    "don't prosecute journalists")

for (nm in dataset_names) {
  base_df <- get(nm)
  df_long <- get(paste0(nm, "_long"))
  for (cl in eval_cols) {
    if (!cl %in% names(base_df)) base_df[[cl]] <- NA_real_
  }
  covars_base <- setdiff(covars, "top_party")
  covars_s <- switch(nm,
    "pool" = c(covars_base, "country", "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")
  )
  join_vars <- setdiff(
    unique(c("RespondentID", pc_vars, "sup_ldem_loo_media",
             "sup_ldem_loo_vote", "sup_ldem_loo_jud",
             "sup_ldem_no_ov", covars_s,
             eval_cols, "redist_1", "redist_2", "gend_pol_1", "gend_pol_2")),
    "country"
  )
  to_join <- dplyr::select(base_df, dplyr::all_of(join_vars))
  dat_long <- dplyr::mutate(
    df_long,
    is_cand_dem = as.integer(cand_system_change %in% pro_dem_labels),
    chosen = as.integer(chosen)
  )
  dat_long <- dplyr::left_join(dat_long, to_join, by = "RespondentID")

  long2 <- dplyr::mutate(
    dat_long,
    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"
  )
  long2 <- dplyr::mutate(
    long2,
    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)
    )
  )
  pair_flags <- long2 %>%
    dplyr::group_by(RespondentID, task) %>%
    dplyr::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"
    ) %>%
    dplyr::filter(has_undem_in & has_dem_out) %>%
    dplyr::select(RespondentID, task)
  dat_h2 <- dplyr::inner_join(long2, pair_flags, by = c("RespondentID", "task"))

  long3 <- dplyr::mutate(
    dat_long,
    cand_redist_1 = dplyr::case_when(cand_tax_spend == "spend more" ~ 4.5,
                                     cand_tax_spend == "spend less" ~ 1.5,
                                     TRUE ~ NA_real_),
    cand_redist_2 = dplyr::case_when(cand_tax_spend == "tax rich more" ~ 4.5,
                                     cand_tax_spend == "tax rich less" ~ 1.5,
                                     TRUE ~ NA_real_),
    cand_gend_1 = dplyr::case_when(cand_gend_policy == "require hire pref" ~ 4.5,
                                   cand_gend_policy == "no hire pref" ~ 1.5,
                                   TRUE ~ NA_real_),
    cand_gend_2 = dplyr::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)
  )
  long3 <- dplyr::rowwise(long3)
  long3 <- dplyr::mutate(
    long3,
    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)
  )
  long3 <- dplyr::ungroup(long3)
  pair_flags_r <- long3 %>%
    dplyr::group_by(RespondentID, task) %>%
    dplyr::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"
    ) %>%
    dplyr::filter(has_undem_prox_r & has_dem_dist_r) %>%
    dplyr::select(RespondentID, task)
  dat_h3_r <- dplyr::inner_join(long3, pair_flags_r, by = c("RespondentID", "task"))
  pair_flags_g <- long3 %>%
    dplyr::group_by(RespondentID, task) %>%
    dplyr::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"
    ) %>%
    dplyr::filter(has_undem_prox_g & has_dem_dist_g) %>%
    dplyr::select(RespondentID, task)
  dat_h3_g <- dplyr::inner_join(long3, pair_flags_g, by = c("RespondentID", "task"))

  assign(paste0("dat_long_", nm), dat_long)
  assign(paste0("dat_h2_", nm), dat_h2)
  assign(paste0("dat_h3_r_", nm), dat_h3_r)
  assign(paste0("dat_h3_g_", nm), dat_h3_g)
}

### ATME and LPM result tables----

country_samples <- c("us", "sk", "pl")
sample_labels <- c(us = "United States", sk = "South Korea", pl = "Poland",
                   pool = "Pooled")
scenario_names <- c("Overall ATME", "Cross-party ATME",
                    "Cross-ideology ATME, redistribution",
                    "Cross-ideology ATME, gender")
scenario_objects <- c("dat_long_", "dat_h2_", "dat_h3_r_", "dat_h3_g_")
hyp_fam <- c("H1", "H2", "H3", "H3")

atme_table <- NULL
for (control_spec in c("Full controls", "Basic controls")) {
  for (s in c(country_samples, "pool")) {
    covars_base <- setdiff(covars, "top_party")
    covars_s <- switch(s,
      "pool" = c(covars_base, "country", "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 (j in seq_along(scenario_names)) {
      dat_use <- get(paste0(scenario_objects[j], s))
      for (v in pc_vars) {
        preds_model <- if (control_spec == "Full controls") c(pc_vars, covars_s) else c(v, covars_s)
        dsub_all <- dat_use %>%
          dplyr::select(chosen, RespondentID, is_cand_dem, dplyr::all_of(preds_model)) %>%
          tidyr::drop_na()
        d1 <- dplyr::filter(dsub_all, is_cand_dem == 1)
        d0 <- dplyr::filter(dsub_all, is_cand_dem == 0)
        m1 <- lm(reformulate(preds_model, "chosen"), data = d1)
        m0 <- lm(reformulate(preds_model, "chosen"), data = d0)
        ct1 <- lmtest::coeftest(m1, vcov = sandwich::vcovCL(m1, cluster = d1$RespondentID))
        ct0 <- lmtest::coeftest(m0, vcov = sandwich::vcovCL(m0, cluster = d0$RespondentID))
        est <- unname(ct1[v, 1] - ct0[v, 1])
        se <- sqrt(unname(ct1[v, 2])^2 + unname(ct0[v, 2])^2)
        pval <- 2 * pnorm(abs(est / se), lower.tail = FALSE)
        atme_table <- dplyr::bind_rows(atme_table, tibble::tibble(
          control_spec = control_spec,
          country = sample_labels[[s]],
          hypothesis_family = hyp_fam[j],
          scenario = scenario_names[j],
          political_culture_measure = pc_labels[[v]],
          estimate = est,
          standard_error = se,
          conf_low = est - 1.96 * se,
          conf_high = est + 1.96 * se,
          exact_p_value = pval,
          n_respondents = length(unique(dsub_all$RespondentID)),
          sample = s,
          pc_var = v
        ))
      }
    }
  }
}

atme_table$fdr_q_value <- NA_real_
for (control_spec in unique(atme_table$control_spec)) {
  for (v in pc_vars) {
    fdr_rows <- atme_table$control_spec == control_spec &
      atme_table$pc_var == v &
      atme_table$sample %in% country_samples
    atme_table$fdr_q_value[fdr_rows] <- p.adjust(
      atme_table$exact_p_value[fdr_rows],
      method = "BH"
    )
  }
}

atme_out <- dplyr::select(
  atme_table,
  control_spec, political_culture_measure, country, hypothesis_family,
  scenario, estimate, standard_error, conf_low, conf_high,
  exact_p_value, fdr_q_value, n_respondents, sample, pc_var
)
atme_full_out <- dplyr::filter(atme_out, control_spec == "Full controls")
atme_basic_out <- dplyr::filter(atme_out, control_spec == "Basic controls")
atme_pooled <- dplyr::filter(atme_full_out, sample == "pool")

write.csv(dplyr::select(atme_full_out, -sample, -pc_var),
          file.path(project_root, "outputs", "atme_exact_p_fdr_full_v3.csv"), row.names = FALSE,
          na = "")
write.csv(dplyr::select(atme_basic_out, -sample, -pc_var),
          file.path(project_root, "outputs", "atme_exact_p_fdr_basic_v3.csv"), row.names = FALSE,
          na = "")
write.csv(dplyr::select(atme_full_out, -control_spec, -sample, -pc_var),
          file.path(project_root, "outputs", "atme_exact_p_fdr_v3.csv"), row.names = FALSE,
          na = "")
write.csv(dplyr::select(atme_pooled, -control_spec, -sample, -pc_var),
          file.path(project_root, "outputs", "pooled_atme_exact_p_v3.csv"), row.names = FALSE,
          na = "")

lpm_table <- NULL
for (control_spec in c("Full controls", "Basic controls")) {
  for (s in c(country_samples, "pool")) {
    covars_base <- setdiff(covars, "top_party")
    covars_s <- switch(s,
      "pool" = c(covars_base, "country", "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 (j in seq_along(scenario_names)) {
      dat_rows <- get(paste0(scenario_objects[j], s))
      choice_dat <- dat_rows %>%
        dplyr::group_by(RespondentID, task) %>%
        dplyr::summarise(y = as.integer(any(chosen == 1 & is_cand_dem == 1)),
                         .groups = "drop")
      cov_join <- get(s) %>%
        dplyr::select(RespondentID, dplyr::all_of(c(pc_vars, covars_s)))
      choice_dat <- dplyr::left_join(choice_dat, cov_join, by = "RespondentID")
      for (v in pc_vars) {
        preds_model <- if (control_spec == "Full controls") c(pc_vars, covars_s) else c(v, covars_s)
        dsub <- choice_dat %>%
          dplyr::select(y, RespondentID, dplyr::all_of(preds_model)) %>%
          tidyr::drop_na()
        m <- lm(reformulate(preds_model, "y"), data = dsub)
        ct <- lmtest::coeftest(m, vcov = sandwich::vcovCL(m, cluster = dsub$RespondentID))
        lpm_table <- dplyr::bind_rows(lpm_table, tibble::tibble(
          control_spec = control_spec,
          country = sample_labels[[s]],
          hypothesis_family = hyp_fam[j],
          scenario = scenario_names[j],
          political_culture_measure = pc_labels[[v]],
          estimate = unname(ct[v, 1]),
          standard_error = unname(ct[v, 2]),
          conf_low = unname(ct[v, 1]) - 1.96 * unname(ct[v, 2]),
          conf_high = unname(ct[v, 1]) + 1.96 * unname(ct[v, 2]),
          exact_p_value = unname(ct[v, 4]),
          n_respondents = length(unique(dsub$RespondentID)),
          sample = s,
          pc_var = v
        ))
      }
    }
  }
}

lpm_table$fdr_q_value <- NA_real_
for (control_spec in unique(lpm_table$control_spec)) {
  for (v in pc_vars) {
    fdr_rows <- lpm_table$control_spec == control_spec &
      lpm_table$pc_var == v &
      lpm_table$sample %in% country_samples
    lpm_table$fdr_q_value[fdr_rows] <- p.adjust(
      lpm_table$exact_p_value[fdr_rows],
      method = "BH"
    )
  }
}

lpm_out <- dplyr::select(
  lpm_table,
  control_spec, political_culture_measure, country, hypothesis_family,
  scenario, estimate, standard_error, conf_low, conf_high,
  exact_p_value, fdr_q_value, n_respondents, sample, pc_var
)
lpm_full_out <- dplyr::filter(lpm_out, control_spec == "Full controls")
lpm_basic_out <- dplyr::filter(lpm_out, control_spec == "Basic controls")
lpm_pooled <- dplyr::filter(lpm_full_out, sample == "pool")

write.csv(dplyr::select(lpm_full_out, -sample, -pc_var),
          file.path(project_root, "outputs", "lpm_exact_p_fdr_full_v3.csv"), row.names = FALSE,
          na = "")
write.csv(dplyr::select(lpm_basic_out, -sample, -pc_var),
          file.path(project_root, "outputs", "lpm_exact_p_fdr_basic_v3.csv"), row.names = FALSE,
          na = "")
write.csv(dplyr::select(lpm_full_out, -control_spec, -sample, -pc_var),
          file.path(project_root, "outputs", "lpm_exact_p_fdr_v3.csv"), row.names = FALSE,
          na = "")
write.csv(dplyr::select(lpm_pooled, -control_spec, -sample, -pc_var),
          file.path(project_root, "outputs", "pooled_lpm_exact_p_v3.csv"), row.names = FALSE,
          na = "")
model_control_audit <- tibble::tibble(
  analysis = c("Country-specific ATME", "Country-specific LPM",
               "Pooled ATME", "Pooled LPM", "SLD overlap ATME"),
  political_culture_controls = c(
    "All six political-culture variables in each model.",
    "All six political-culture variables in each model.",
    "All six political-culture variables in each model.",
    "All six political-culture variables in each model.",
    "Reduced support for liberal democracy plus the five other political-culture variables."
  ),
  non_culture_controls = c(
    "gender_male, age_1, educ_5f, income_cat, region_pool, pol_interest, gend_pol_left, redist_pol_left, and country-specific party thermometer evaluations.",
    "gender_male, age_1, educ_5f, income_cat, region_pool, pol_interest, gend_pol_left, redist_pol_left, and country-specific party thermometer evaluations.",
    "gender_male, age_1, educ_5f, income_cat, region_pool, pol_interest, gend_pol_left, redist_pol_left, country fixed effects, and top_party.",
    "gender_male, age_1, educ_5f, income_cat, region_pool, pol_interest, gend_pol_left, redist_pol_left, country fixed effects, and top_party.",
    "gender_male, age_1, educ_5f, income_cat, region_pool, pol_interest, gend_pol_left, redist_pol_left, and country-specific party thermometer evaluations."
  )
)
write.csv(model_control_audit, file.path(project_root, "outputs", "model_control_audit_v3.csv"),
          row.names = FALSE)

### Figure 1 descriptive regression slopes----

slope_table <- NULL
fig1_panel_names <- c("Overall", "Cross-party",
                      "Cross-ideology redistribution",
                      "Cross-ideology gender")
for (s in c(country_samples, "pool")) {
  for (j in seq_along(fig1_panel_names)) {
    dat_rows <- get(paste0(scenario_objects[j], s))
    choice_dat <- dat_rows %>%
      dplyr::group_by(RespondentID, task) %>%
      dplyr::summarise(y = as.integer(any(chosen == 1 & is_cand_dem == 1)),
                       .groups = "drop")
    cov_join <- get(s) %>% dplyr::select(RespondentID, dplyr::all_of(pc_vars))
    choice_dat <- dplyr::left_join(choice_dat, cov_join, by = "RespondentID")
    for (v in pc_vars) {
      dsub <- choice_dat %>%
        dplyr::select(y, RespondentID, dplyr::all_of(v)) %>%
        tidyr::drop_na()
      m <- lm(reformulate(v, "y"), data = dsub)
      ct <- lmtest::coeftest(m, vcov = sandwich::vcovCL(m, cluster = dsub$RespondentID))
      slope <- unname(ct[v, 1])
      se <- unname(ct[v, 2])
      slope_table <- dplyr::bind_rows(slope_table, tibble::tibble(
        panel = fig1_panel_names[j],
        country = sample_labels[[s]],
        political_culture_measure = pc_labels[[v]],
        slope = slope,
        standard_error = se,
        conf_low = slope - 1.96 * se,
        conf_high = slope + 1.96 * se,
        exact_p_value = unname(ct[v, 4]),
        n_respondents = length(unique(dsub$RespondentID))
      ))
    }
  }
}
write.csv(slope_table, file.path(project_root, "outputs", "figure1_descriptive_slopes_v3.csv"),
          row.names = FALSE)
