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 R&R background, reliability, CFA, correlations, and sample diagnostics.
#### Builds data from raw files and saves background 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)"
)

country_samples <- c("us", "sk", "pl")
sample_labels <- c(us = "United States", sk = "South Korea", pl = "Poland",
                   pool = "Pooled")
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)
}

### Scale reliability----

reverse_to_observed_range <- function(x) {
  rng <- range(x, na.rm = TRUE)
  sum(rng) - x
}

reliability_alpha <- function(mat) {
  suppressWarnings(tryCatch(
    psych::alpha(mat, check.keys = TRUE)$total$std.alpha,
    error = function(e) NA_real_
  ))
}

reliability_alpha_poly <- function(mat) {
  mat_complete <- mat[stats::complete.cases(mat), , drop = FALSE]
  if (nrow(mat_complete) < 10) return(NA_real_)
  suppressWarnings(tryCatch({
    poly_cor <- psych::polychoric(mat_complete, max.cat = 20,
                                  progress = FALSE)$rho
    psych::alpha(poly_cor, n.obs = nrow(mat_complete),
                 check.keys = TRUE)$total$std.alpha
  }, error = function(e) NA_real_))
}

reliability <- NULL
scale_specs <- list(
  "Democracy-autocracy preference" = c("sup_dem1", "sup_dem2", "sup_dem3", "sup_dem4"),
  "Support for liberal democracy" = c("sup_ld_1", "sup_ld_2", "sup_ld_3",
                                      "sup_ld_4", "sup_ld_5", "sup_ld_6",
                                      "sup_ld_7"),
  "Institutional trust" = c("trust_leg", "trust_parl", "trust_civ", "trust_part"),
  "Political tolerance" = c("pol_tol_1", "pol_tol_2", "pol_tol_3"),
  "Emancipative values" =
    c("autonomy_z", "equality_z", "choice_z", "voice_z"),
  "Authoritarian values" = c("val_conf_1_ips",
                             "val_trad_1_ips",
                             "val_selfd_1_ips",
                             "val_stim_1_ips",
                             "val_conf_2_ips",
                             "val_trad_2_ips",
                             "val_selfd_2_ips",
                             "val_stim_2_ips")
)
ordinal_alpha_scales <- c("Democracy-autocracy preference",
                          "Support for liberal democracy",
                          "Institutional trust",
                          "Political tolerance")
for (s in country_samples) {
  df <- get(s)
  for (measure in names(scale_specs)) {
    items <- scale_specs[[measure]]
    mat <- df[, items, drop = FALSE]
    if (measure == "Democracy-autocracy preference") {
      mat$sup_dem2 <- reverse_to_observed_range(mat$sup_dem2)
      mat$sup_dem4 <- reverse_to_observed_range(mat$sup_dem4)
    }
    if (measure == "Support for liberal democracy") {
      for (item in names(mat)) mat[[item]] <- reverse_to_observed_range(mat[[item]])
    }
    if (measure == "Political tolerance") {
      for (item in names(mat)) mat[[item]] <- reverse_to_observed_range(mat[[item]])
    }
    if (measure == "Authoritarian values") {
      mat$val_selfd_1_ips <- -mat$val_selfd_1_ips
      mat$val_stim_1_ips <- -mat$val_stim_1_ips
      mat$val_selfd_2_ips <- -mat$val_selfd_2_ips
      mat$val_stim_2_ips <- -mat$val_stim_2_ips
    }
    valid_n <- sum(stats::complete.cases(mat))
    pct_missing <- mean(!stats::complete.cases(mat)) * 100
    alpha_method <- ifelse(measure %in% ordinal_alpha_scales,
                           "Polychoric correlations",
                           "Pearson correlations")
    alpha_val <- if (measure %in% ordinal_alpha_scales) {
      reliability_alpha_poly(mat)
    } else {
      reliability_alpha(mat)
    }
    reliability <- dplyr::bind_rows(reliability, tibble::tibble(
      country = sample_labels[[s]],
      measure = measure,
      number_of_items = length(items),
      cronbach_alpha = alpha_val,
      alpha_method = alpha_method,
      pct_missing = pct_missing,
      valid_n = valid_n
    ))
  }
}
write.csv(reliability, file.path(project_root, "outputs", "scale_reliability_v3.csv"),
          row.names = FALSE)

### Measurement equivalence CFA diagnostics----

cfa_fit_indices <- function(fit) {
  fm <- lavaan::fitMeasures(fit)
  pick_fit <- function(preferred, fallback) {
    if (preferred %in% names(fm)) fm[[preferred]] else if (fallback %in% names(fm)) fm[[fallback]] else NA_real_
  }
  c(
    cfi = as.numeric(pick_fit("cfi.robust", "cfi")),
    rmsea = as.numeric(pick_fit("rmsea.robust", "rmsea")),
    srmr = as.numeric(if ("srmr" %in% names(fm)) fm[["srmr"]] else NA_real_),
    chisq = as.numeric(if ("chisq.scaled" %in% names(fm)) fm[["chisq.scaled"]] else fm[["chisq"]]),
    df = as.numeric(if ("df.scaled" %in% names(fm)) fm[["df.scaled"]] else fm[["df"]])
  )
}

run_cfa <- function(dat, sample_name, measure, structure, model, items,
                    estimator = "WLSMV", ordered = TRUE,
                    fallback_model = NULL,
                    fallback_structure = "1-factor CFA with equal loadings") {
  dat_use <- dat[, unique(items), drop = FALSE]
  estimators <- if (ordered) {
    unique(c(estimator, "DWLS", "ULS"))
  } else {
    unique(c(estimator, "MLR", "ML", "ULS"))
  }
  fit <- NULL
  final_model <- model
  final_structure <- structure
  final_estimator <- estimator
  final_note <- NA_character_

  for (est in estimators) {
    fit_try <- tryCatch(
      lavaan::cfa(
        model = model,
        data = dat_use,
        estimator = est,
        ordered = if (ordered) items else NULL,
        std.lv = TRUE,
        missing = if (ordered || est %in% c("ULS", "DWLS", "WLSMV")) "listwise" else "fiml"
      ),
      error = function(e) e
    )
    if (!inherits(fit_try, "error") &&
        isTRUE(lavaan::lavInspect(fit_try, "converged"))) {
      fit <- fit_try
      final_estimator <- est
      break
    }
    final_note <- if (inherits(fit_try, "error")) {
      conditionMessage(fit_try)
    } else {
      "lavaan did not converge"
    }
  }

  if (is.null(fit) && !is.null(fallback_model)) {
    for (est in estimators) {
      fit_try <- tryCatch(
        lavaan::cfa(
          model = fallback_model,
          data = dat_use,
          estimator = est,
          ordered = if (ordered) items else NULL,
          std.lv = TRUE,
          missing = if (ordered || est %in% c("ULS", "DWLS", "WLSMV")) "listwise" else "fiml"
        ),
        error = function(e) e
      )
      if (!inherits(fit_try, "error") &&
          isTRUE(lavaan::lavInspect(fit_try, "converged"))) {
        fit <- fit_try
        final_model <- fallback_model
        final_structure <- fallback_structure
        final_estimator <- est
        final_note <- "Primary CFA failed; equal-loading fallback converged"
        break
      }
      final_note <- if (inherits(fit_try, "error")) {
        conditionMessage(fit_try)
      } else {
        "fallback CFA did not converge"
      }
    }
  }

  if (is.null(fit)) {
    return(tibble::tibble(
      country = sample_name,
      measure = measure,
      model_structure = final_structure,
      estimator = final_estimator,
      ordered_items = if (ordered) paste(items, collapse = "; ") else NA_character_,
      converged = FALSE,
      saturated_or_no_fit_metrics = NA,
      cfi = NA_real_,
      rmsea = NA_real_,
      srmr = NA_real_,
      chisq = NA_real_,
      df = NA_real_,
      valid_n = sum(stats::complete.cases(dat_use)),
      note = final_note
    ))
  }
  converged <- TRUE
  fit_vals <- cfa_fit_indices(fit)
  names(fit_vals) <- c("cfi", "rmsea", "srmr", "chisq", "df")
  tibble::tibble(
    country = sample_name,
    measure = measure,
    model_structure = final_structure,
    estimator = final_estimator,
    ordered_items = if (ordered) paste(items, collapse = "; ") else NA_character_,
    converged = converged,
    saturated_or_no_fit_metrics = converged && !is.na(fit_vals[["df"]]) && fit_vals[["df"]] <= 0,
    cfi = if (is.na(fit_vals[["df"]]) || fit_vals[["df"]] <= 0) NA_real_ else fit_vals[["cfi"]],
    rmsea = if (is.na(fit_vals[["df"]]) || fit_vals[["df"]] <= 0) NA_real_ else fit_vals[["rmsea"]],
    srmr = if (is.na(fit_vals[["df"]]) || fit_vals[["df"]] <= 0) NA_real_ else fit_vals[["srmr"]],
    chisq = fit_vals[["chisq"]],
    df = fit_vals[["df"]],
    valid_n = sum(stats::complete.cases(dat_use)),
    note = final_note
  )
}

cfa_results <- NULL
for (s in country_samples) {
  df <- get(s)
  cfa_df <- df
  cfa_df$sup_dem2_r <- reverse_to_observed_range(cfa_df$sup_dem2)
  cfa_df$sup_dem4_r <- reverse_to_observed_range(cfa_df$sup_dem4)
  cfa_df$sup_ld_1_r <- reverse_to_observed_range(cfa_df$sup_ld_1)
  cfa_df$sup_ld_2_r <- reverse_to_observed_range(cfa_df$sup_ld_2)
  cfa_df$sup_ld_3_r <- reverse_to_observed_range(cfa_df$sup_ld_3)
  cfa_df$sup_ld_4_r <- reverse_to_observed_range(cfa_df$sup_ld_4)
  cfa_df$sup_ld_5_r <- reverse_to_observed_range(cfa_df$sup_ld_5)
  cfa_df$sup_ld_6_r <- reverse_to_observed_range(cfa_df$sup_ld_6)
  cfa_df$sup_ld_7_r <- reverse_to_observed_range(cfa_df$sup_ld_7)
  cfa_df$pol_tol_1_r <- reverse_to_observed_range(cfa_df$pol_tol_1)
  cfa_df$pol_tol_2_r <- reverse_to_observed_range(cfa_df$pol_tol_2)
  cfa_df$pol_tol_3_r <- reverse_to_observed_range(cfa_df$pol_tol_3)
  cfa_df$gend_eq_1_r <- reverse_to_observed_range(cfa_df$gend_eq_1)
  cfa_df$gend_eq_2_r <- reverse_to_observed_range(cfa_df$gend_eq_2)
  cfa_df$gend_eq_3_r <- reverse_to_observed_range(cfa_df$gend_eq_3)
  cfa_df$val_selfd_1_ips_r <- -cfa_df$val_selfd_1_ips
  cfa_df$val_stim_1_ips_r <- -cfa_df$val_stim_1_ips
  cfa_df$val_selfd_2_ips_r <- -cfa_df$val_selfd_2_ips
  cfa_df$val_stim_2_ips_r <- -cfa_df$val_stim_2_ips

  cfa_specs <- list(
    list(
      measure = "Democracy-autocracy preference",
      structure = "1-factor CFA",
      model = "demaut =~ sup_dem1 + sup_dem2_r + sup_dem3 + sup_dem4_r",
      items = c("sup_dem1", "sup_dem2_r", "sup_dem3", "sup_dem4_r"),
      estimator = "WLSMV",
      ordered = TRUE
    ),
    list(
      measure = "Support for liberal democracy",
      structure = "1-factor CFA",
      model = "sld =~ sup_ld_1_r + sup_ld_2_r + sup_ld_3_r + sup_ld_4_r + sup_ld_5_r + sup_ld_6_r + sup_ld_7_r",
      items = paste0("sup_ld_", 1:7, "_r"),
      estimator = "WLSMV",
      ordered = TRUE
    ),
    list(
      measure = "Institutional trust",
      structure = "1-factor CFA",
      model = "trust =~ trust_leg + trust_parl + trust_civ + trust_part",
      items = c("trust_leg", "trust_parl", "trust_civ", "trust_part"),
      estimator = "MLR",
      ordered = FALSE
    ),
    list(
      measure = "Political tolerance",
      structure = "loading constrained 1-factor CFA",
      model = "tolerance =~ a*pol_tol_1_r + a*pol_tol_2_r + a*pol_tol_3_r",
      items = c("pol_tol_1_r", "pol_tol_2_r", "pol_tol_3_r"),
      estimator = "WLSMV",
      ordered = TRUE
    ),
    list(
      measure = "Emancipative values",
      structure = "hierarchical 4-factor CFA",
      model = paste(
        "autonomy =~ child_ind + child_imag + child_no_obed",
        "equality =~ gend_eq_1_r + gend_eq_2_r + gend_eq_3_r",
        "choice =~ gay_just + abort_just + divorc_just",
        "voice =~ postmat_voice_1 + postmat_voice_2",
        "emancipative =~ autonomy + equality + choice + voice",
        sep = "\n"
      ),
      items = c("child_ind", "child_imag", "child_no_obed",
                "gend_eq_1_r", "gend_eq_2_r", "gend_eq_3_r",
                "gay_just", "abort_just", "divorc_just",
                "postmat_voice_1", "postmat_voice_2"),
      estimator = "WLSMV",
      ordered = TRUE
    ),
    list(
      measure = "Authoritarian-libertarian values",
      structure = "1-factor CFA on ipsatised items; one redundant ipsatised item omitted for covariance identification",
      model = paste(
        "authlib =~ val_conf_1_ips + val_trad_1_ips + val_conf_2_ips + val_trad_2_ips +",
        "           val_selfd_1_ips_r + val_stim_1_ips_r + val_selfd_2_ips_r",
        sep = "\n"
      ),
      items = c("val_conf_1_ips", "val_trad_1_ips", "val_conf_2_ips", "val_trad_2_ips",
                "val_selfd_1_ips_r", "val_stim_1_ips_r", "val_selfd_2_ips_r"),
      estimator = "ULS",
      ordered = FALSE
    )
  )

  for (spec in cfa_specs) {
    cfa_results <- dplyr::bind_rows(
      cfa_results,
      run_cfa(
        cfa_df,
        sample_labels[[s]],
        spec$measure,
        spec$structure,
        spec$model,
        spec$items,
        spec$estimator,
        spec$ordered,
        fallback_model = paste(
          "equal_factor =~ a*", spec$items[1],
          paste0(" + a*", spec$items[-1], collapse = ""),
          sep = ""
        )
      )
    )
  }
}
cfa_outfile <- file.path(project_root, "outputs", "measurement_equivalence_cfa_v3.csv")
cfa_write_ok <- tryCatch({
  write.csv(cfa_results, cfa_outfile, row.names = FALSE)
  TRUE
}, error = function(e) FALSE)
if (!cfa_write_ok) {
  write.csv(cfa_results,
            file.path(project_root, "outputs", "measurement_equivalence_cfa_v3_unlocked_copy.csv"),
            row.names = FALSE)
}

### Correlation matrices----

cor_label_map <- c(
  demaut_pref = "Dem-aut preference",
  sup_ldem = "Support liberal democracy",
  inst_trust = "Institutional trust",
  pol_tolerance = "Political tolerance",
  autonomy_values = "Auth. values (rev)",
  emancip_val = "Emancipative values",
  gend_pol_left = "Gender policy (left)",
  redist_pol_left = "Redistributive policy (left)",
  eval_pis_pl = "Evals: PiS",
  eval_kon_pl = "Evals: Kon",
  eval_ko_pl = "Evals: KO",
  eval_lew_pl = "Evals: Lewica",
  eval_dem_us = "Evals: Democrat",
  eval_rep_us = "Evals: Republican",
  eval_dpk_kr = "Evals: DPK",
  eval_ppp_kr = "Evals: PPP",
  dem_choice_rate = "Pro-democratic choices"
)

cor_specs <- list(
  pool = list(
    data = pool,
    vars = c("demaut_pref", "sup_ldem", "inst_trust", "pol_tolerance",
             "autonomy_values", "emancip_val", "gend_pol_left",
             "redist_pol_left", "dem_choice_rate"),
    caption = "Pooled sample",
    outfile = "cor_tab_pool_v3"
  ),
  pl = list(
    data = pl,
    vars = c("demaut_pref", "sup_ldem", "inst_trust", "pol_tolerance",
             "autonomy_values", "emancip_val", "gend_pol_left",
             "redist_pol_left", "eval_pis_pl", "eval_kon_pl", "eval_ko_pl",
             "eval_lew_pl", "dem_choice_rate"),
    caption = "Polish sample",
    outfile = "cor_tab_pl_v3"
  ),
  us = list(
    data = us,
    vars = c("demaut_pref", "sup_ldem", "inst_trust", "pol_tolerance",
             "autonomy_values", "emancip_val", "gend_pol_left",
             "redist_pol_left", "eval_dem_us", "eval_rep_us",
             "dem_choice_rate"),
    caption = "US sample",
    outfile = "cor_tab_us_v3"
  ),
  sk = list(
    data = sk,
    vars = c("demaut_pref", "sup_ldem", "inst_trust", "pol_tolerance",
             "autonomy_values", "emancip_val", "gend_pol_left",
             "redist_pol_left", "eval_dpk_kr", "eval_ppp_kr",
             "dem_choice_rate"),
    caption = "South Korean sample",
    outfile = "cor_tab_sk_v3"
  )
)

for (nm in names(cor_specs)) {
  cor_df <- cor_specs[[nm]]$data
  cor_vars <- intersect(cor_specs[[nm]]$vars, names(cor_df))
  cor_matrix <- stats::cor(cor_df[, cor_vars, drop = FALSE],
                           use = "pairwise.complete.obs")
  cor_rounded <- round(cor_matrix, 2)
  cor_labels <- cor_label_map[cor_vars]
  cor_labels[is.na(cor_labels)] <- cor_vars[is.na(cor_labels)]
  colnames(cor_rounded) <- cor_labels
  rownames(cor_rounded) <- cor_labels
  write.csv(cor_rounded,
            file.path(project_root, "outputs", paste0(cor_specs[[nm]]$outfile, ".csv")))

}

### Sample composition diagnostics----

sample_demo <- NULL
for (s in country_samples) {
  df <- get(s)
  if (s == "us") {
    df$age_group <- dplyr::case_when(
      df$age_1 < 35 ~ "18-34",
      df$age_1 < 50 ~ "35-49",
      df$age_1 < 65 ~ "50-64",
      df$age_1 >= 65 ~ "65+",
      TRUE ~ NA_character_
    )
  } else {
    df$age_group <- dplyr::case_when(
      df$age_1 >= 20 & df$age_1 < 30 ~ "20-29",
      df$age_1 >= 30 & df$age_1 < 45 ~ "30-44",
      df$age_1 >= 45 & df$age_1 < 60 ~ "45-59",
      df$age_1 >= 60 ~ "60+",
      TRUE ~ NA_character_
    )
  }
  df$gender_cat <- dplyr::case_when(
    df$gender == 1 ~ "Male",
    df$gender == 2 ~ "Female",
    df$gender == 9 ~ "Other",
    df$gender == 98 ~ "Refused",
    TRUE ~ NA_character_
  )
  df$gender_cat <- factor(df$gender_cat,
                          levels = c("Male", "Female", "Other", "Refused"))
  gender_tab <- as.data.frame(prop.table(table(df$gender_cat, useNA = "no")))
  names(gender_tab) <- c("category", "sample_share")
  gender_tab <- dplyr::filter(gender_tab, sample_share > 0)
  gender_tab$variable <- "Gender"
  gender_tab$country <- sample_labels[[s]]
  age_tab <- as.data.frame(prop.table(table(df$age_group, useNA = "no")))
  names(age_tab) <- c("category", "sample_share")
  age_tab$variable <- "Age group"
  age_tab$country <- sample_labels[[s]]
  region_tab <- as.data.frame(prop.table(table(df$region_pool, useNA = "no")))
  names(region_tab) <- c("category", "sample_share")
  region_tab$variable <- "Region"
  region_tab$country <- sample_labels[[s]]
  sample_demo <- dplyr::bind_rows(
    sample_demo,
    dplyr::select(gender_tab, country, variable, category, sample_share),
    dplyr::select(age_tab, country, variable, category, sample_share),
    dplyr::select(region_tab, country, variable, category, sample_share)
  )
  party_vars <- switch(s,
    "us" = c("eval_dem_us", "eval_rep_us"),
    "sk" = c("eval_dpk_kr", "eval_ppp_kr"),
    "pl" = c("eval_pis_pl", "eval_kon_pl", "eval_ko_pl", "eval_lew_pl")
  )
  party_names <- switch(s,
    "us" = c("Democrat", "Republican"),
    "sk" = c("DPK", "PPP"),
    "pl" = c("PiS", "Konfederacja", "KO", "Lewica")
  )
  for (i in seq_along(party_vars)) {
    vals <- df[[party_vars[i]]]
    party_row <- tibble::tibble(
      country = sample_labels[[s]],
      variable = "Partisan support",
      category = party_names[i],
      sample_share = mean(vals >= 6, na.rm = TRUE)
    )
    sample_demo <- dplyr::bind_rows(sample_demo, party_row)
  }
}
write.csv(sample_demo, file.path(project_root, "outputs", "sample_demographic_margins_v3.csv"),
          row.names = FALSE)

party_therm <- tibble::tibble(
  country = c(rep("United States", 2), rep("South Korea", 2), rep("Poland", 4)),
  party = c("Democrat", "Republican", "DPK", "PPP", "PiS", "Konfederacja", "KO", "Lewica"),
  mean_thermometer = c(mean(us$eval_dem_us, na.rm = TRUE),
                       mean(us$eval_rep_us, na.rm = TRUE),
                       mean(sk$eval_dpk_kr, na.rm = TRUE),
                       mean(sk$eval_ppp_kr, na.rm = TRUE),
                       mean(pl$eval_pis_pl, na.rm = TRUE),
                       mean(pl$eval_kon_pl, na.rm = TRUE),
                       mean(pl$eval_ko_pl, na.rm = TRUE),
                       mean(pl$eval_lew_pl, na.rm = TRUE)),
  valid_n = c(sum(!is.na(us$eval_dem_us)), sum(!is.na(us$eval_rep_us)),
              sum(!is.na(sk$eval_dpk_kr)), sum(!is.na(sk$eval_ppp_kr)),
              sum(!is.na(pl$eval_pis_pl)), sum(!is.na(pl$eval_kon_pl)),
              sum(!is.na(pl$eval_ko_pl)), sum(!is.na(pl$eval_lew_pl)))
)
write.csv(party_therm, file.path(project_root, "outputs", "party_thermometer_means_v3.csv"),
          row.names = FALSE)

# Combine sample margins with already-downloaded population benchmarks.
benchmark_margins <- read.csv(file.path(project_root, "data", "population_benchmark_margins_v3.csv"))
sample_margins <- read.csv(file.path(project_root, "outputs", "sample_demographic_margins_v3.csv"))
sample_benchmark <- dplyr::left_join(
  sample_margins,
  benchmark_margins,
  by = c("country", "variable", "category")
)
sample_benchmark$source[
  sample_benchmark$variable == "Partisan support" &
    is.na(sample_benchmark$source)
] <- "Survey thermometer rating only; no external partisan benchmark used."
sample_benchmark$benchmark_note[
  sample_benchmark$variable == "Partisan support" &
    is.na(sample_benchmark$benchmark_note)
] <- "Sample share is percent with thermometer rating 6-10."
sample_benchmark$sample_percent <- round(sample_benchmark$sample_share * 100, 1)
sample_benchmark$population_percent <- round(sample_benchmark$benchmark_share * 100, 1)
sample_benchmark <- dplyr::select(
  sample_benchmark,
  country, variable, category, sample_percent, population_percent,
  source, benchmark_note
)
write.csv(sample_benchmark,
          file.path(project_root, "outputs", "sample_benchmark_comparison_all_v3.csv"),
          row.names = FALSE, na = "")

for (cn in unique(sample_benchmark$country)) {
  one_sample <- dplyr::filter(sample_benchmark, country == cn)
  file_stub <- tolower(gsub("[^A-Za-z0-9]+", "_", cn))
  write.csv(one_sample,
            file.path(project_root, "outputs", paste0("sample_benchmark_", file_stub, "_v3.csv")),
            row.names = FALSE, na = "")
}

