# ==============================================================================
  # FULL REPLICATION CODE (WEIGHTED): Who Should Fix the Climate?
  # ==============================================================================
  # Every model uses Combinwt_new_hh as a weight.
  # Unweighted counterparts are estimated alongside for comparison.
  # ==============================================================================
  
  # --- 0. Setup & Libraries ---
  rm(list = ls())
  library(haven)
  library(dplyr)
  library(tidyr)
  library(tibble)          # <-- provides column_to_rownames()
  library(sandwich)
  library(lmtest)
  library(marginaleffects)
  library(ggplot2)
  library(broom)
  
  # --- LOAD THE DATA ---
  setwd("C:/Users/agbog/Downloads")
  dat <- read_sav("AB_R10.Merge Dataset_25Feb26.38ctry.release.sav")
  
  cat("Raw sample loaded:", nrow(dat), "observations.\n")
  
  # ==============================================================================
  # 1. DATA PREPARATION: CLIMATE-AWARE SAMPLE
  # ==============================================================================
  
  df_lit <- dat %>%
    mutate(
      country       = as_factor(COUNTRY),
      resp_id       = RESPNO,
      climate_aware = case_when(Q60A == 1 ~ 1, TRUE ~ 0),
      
      # --- Attribution Typology ---
      cause_foreign = case_when(
        Q61B == 2 ~ 1,
        Q61B == 1 ~ 0,
        TRUE ~ NA_real_
      ),
      remedy_foreign = case_when(
        Q62D == 4 ~ 1,
        Q62D %in% c(1, 2, 3, 5) ~ 0,
        TRUE ~ NA_real_
      ),
      
      # --- Evaluations of Influence (Rescaled 0-1) ---
      influence_china = case_when(
        Q67B %in% 1:5 ~ (Q67B - 1) / 4,
        Q67B %in% c(8, 9) ~ 0.5,
        TRUE ~ NA_real_
      ),
      influence_us = case_when(
        Q67C %in% 1:5 ~ (Q67C - 1) / 4,
        Q67C %in% c(8, 9) ~ 0.5,
        TRUE ~ NA_real_
      ),
      influence_eu = case_when(
        Q67H %in% 1:5 ~ (Q67H - 1) / 4,
        Q67H %in% c(8, 9) ~ 0.5,
        TRUE ~ NA_real_
      ),
      
      # --- Controls ---
      tertiary = case_when(
        Q96 %in% c(98, 99, -1) ~ NA_real_,
        Q96 >= 6 ~ 1,
        TRUE ~ 0
      ),
      female = case_when(Q101 == 2 ~ 1, Q101 == 1 ~ 0, TRUE ~ NA_real_),
      age    = ifelse(Q1 %in% c(998, 999, -1), NA, Q1),
      urban  = case_when(URBRUR == 1 ~ 1, URBRUR == 2 ~ 0, TRUE ~ NA_real_),
      lived_pov  = LivedPoverty,
      news_count = rowSums(cbind(
        Q65A %in% 1:4, Q65B %in% 1:4, Q65C %in% 1:4,
        Q65D %in% 1:4, Q65E %in% 1:4), na.rm = FALSE),
      
      # --- Typology ---
      attribution_type = case_when(
        cause_foreign == 0 & remedy_foreign == 0 ~ "Consistent internaliser",
        cause_foreign == 1 & remedy_foreign == 1 ~ "Consistent externaliser",
        cause_foreign == 1 & remedy_foreign == 0 ~ "Sovereign adapter",
        cause_foreign == 0 & remedy_foreign == 1 ~ "Inconsistent externaliser",
        TRUE ~ NA_character_
      ),
      attribution_type = factor(attribution_type,
                                levels = c("Consistent internaliser",
                                           "Sovereign adapter",
                                           "Consistent externaliser",
                                           "Inconsistent externaliser")),
      w = Combinwt_new_hh
    ) %>%
    filter(climate_aware == 1, !is.na(attribution_type))
  
  cat("\nAnalytic Sample Size (Climate-Aware):", nrow(df_lit), "\n")
  cat("Missing weights:", sum(is.na(df_lit$w)), "\n")
  cat("Weight summary:\n")
  print(summary(df_lit$w))
  
  # --- Weight concentration diagnostic ---
  w_sorted <- sort(df_lit$w, decreasing = TRUE)
  cat("\nTop 1% of cases carry",
      round(sum(w_sorted[1:ceiling(0.01 * length(w_sorted))]) / sum(w_sorted) * 100, 1),
      "% of total weight\n")
  cat("Top 5% of cases carry",
      round(sum(w_sorted[1:ceiling(0.05 * length(w_sorted))]) / sum(w_sorted) * 100, 1),
      "% of total weight\n")
  
  cat("\nTop 10 weights (country, resp_id, w):\n")
  df_lit %>%
    arrange(desc(w)) %>%
    select(country, resp_id, w) %>%
    head(10) %>%
    print()
  
  # ==============================================================================
  # 2. DESCRIPTIVE STATISTICS (UNWEIGHTED VS WEIGHTED)
  # ==============================================================================
  
  cat("\n--- Table 1: Distribution of Attribution Frameworks ---\n")
  tab_unw <- prop.table(table(df_lit$attribution_type)) * 100
  tab_w   <- prop.table(xtabs(w ~ attribution_type, data = df_lit)) * 100
  print(round(cbind(Unweighted = tab_unw, Weighted = tab_w), 3))
  
  cat("\n--- Cause vs Remedy Asymmetry ---\n")
  cat("Blame Foreign Cause (unweighted):",
      round(mean(df_lit$cause_foreign, na.rm = TRUE) * 100, 1), "%\n")
  cat("Blame Foreign Cause (weighted):  ",
      round(weighted.mean(df_lit$cause_foreign, df_lit$w, na.rm = TRUE) * 100, 1), "%\n")
  cat("Expect Foreign Remedy (unweighted):",
      round(mean(df_lit$remedy_foreign, na.rm = TRUE) * 100, 1), "%\n")
  cat("Expect Foreign Remedy (weighted):  ",
      round(weighted.mean(df_lit$remedy_foreign, df_lit$w, na.rm = TRUE) * 100, 1), "%\n")
  
  # Long format
  stack_geo <- df_lit %>%
    select(resp_id, country, attribution_type, w, tertiary, female, age,
           urban, lived_pov, news_count, influence_china, influence_us, influence_eu) %>%
    pivot_longer(cols = c(influence_china, influence_us, influence_eu),
                 names_to = "power", values_to = "rating") %>%
    filter(!is.na(rating)) %>%
    mutate(power = factor(power,
                          levels = c("influence_china", "influence_us", "influence_eu"),
                          labels = c("China", "US", "EU")))
  
  cat("\n--- Table 2: Mean Evaluations by Power and Type ---\n")
  
  # Unweighted
  tab2_unw <- stack_geo %>%
    group_by(attribution_type, power) %>%
    summarise(mean_u = mean(rating, na.rm = TRUE), .groups = "drop") %>%
    pivot_wider(names_from = power, values_from = mean_u) %>%
    column_to_rownames("attribution_type") %>%
    as.matrix()
  
  # Weighted (summarise keeps rating and w length-matched)
  tab2_w <- stack_geo %>%
    group_by(attribution_type, power) %>%
    summarise(mean_w = weighted.mean(rating, w = w, na.rm = TRUE), .groups = "drop") %>%
    pivot_wider(names_from = power, values_from = mean_w) %>%
    column_to_rownames("attribution_type") %>%
    as.matrix()
  
  cat("\nUnweighted means:\n");  print(round(tab2_unw, 3))
  cat("\nWeighted means:\n");    print(round(tab2_w, 3))
  cat("\nWeighted US-China gap by type:\n")
  print(round(tab2_w[, "US"] - tab2_w[, "China"], 3))
  cat("\nUnweighted US-China gap by type:\n")
  print(round(tab2_unw[, "US"] - tab2_unw[, "China"], 3))
  
  # ==============================================================================
  # 3. MAIN MODELS: CONDITIONAL PENALTY
  # ==============================================================================
  
  # --- Weighted ---
  m_typ_w <- lm(
    rating ~ power * attribution_type + tertiary + female + age + urban +
      lived_pov + news_count + country,
    data = stack_geo, weights = w
  )
  ct_typ_w <- coeftest(m_typ_w, vcov = vcovCL(m_typ_w, cluster = ~ resp_id,
                                              data = stack_geo))
  
  # --- Unweighted counterpart ---
  m_typ_u <- lm(
    rating ~ power * attribution_type + tertiary + female + age + urban +
      lived_pov + news_count + country,
    data = stack_geo
  )
  ct_typ_u <- coeftest(m_typ_u, vcov = vcovCL(m_typ_u, cluster = ~ resp_id,
                                              data = stack_geo))
  
  cat("\n--- Table 3 (WEIGHTED): US-China Gap by Attribution Type ---\n")
  print(ct_typ_w[grep("power", rownames(ct_typ_w)), ])
  
  cat("\n--- Table 3 (UNWEIGHTED, for comparison) ---\n")
  print(ct_typ_u[grep("power", rownames(ct_typ_u)), ])
  
  # --- Weighted marginal means ---
  mm_w <- avg_predictions(
    m_typ_w,
    variables = c("power", "attribution_type"),
    vcov = vcovCL(m_typ_w, cluster = ~ resp_id, data = stack_geo),
    wts = "w"
  )
  cat("\n--- Marginal Means (WEIGHTED) ---\n")
  print(mm_w)
  
  # --- Unweighted marginal means ---
  mm_u <- avg_predictions(
    m_typ_u,
    variables = c("power", "attribution_type"),
    vcov = vcovCL(m_typ_u, cluster = ~ resp_id, data = stack_geo)
  )
  cat("\n--- Marginal Means (UNWEIGHTED) ---\n")
  print(mm_u)
  
  # --- H2: explicit US-China gap by type, weighted ---
  cat("\n--- H2: US-China Gap by Type, Weighted (avg_comparisons) ---\n")
  h2_w <- avg_comparisons(
    m_typ_w,
    variables = list(power = c("China", "US")),
    by = "attribution_type",
    vcov = vcovCL(m_typ_w, cluster = ~ resp_id, data = stack_geo),
    wts = "w"
  )
  print(h2_w)
  
  # ==============================================================================
  # 4. CLIMATE AWARENESS COMPARISON
  # ==============================================================================
  
  df_non_aware <- dat %>%
    filter(Q60A == 0) %>%
    mutate(
      resp_id = RESPNO,
      country = as_factor(COUNTRY),
      influence_china = case_when(Q67B %in% 1:5 ~ (Q67B - 1) / 4,
                                  Q67B %in% c(8, 9) ~ 0.5, TRUE ~ NA_real_),
      influence_us    = case_when(Q67C %in% 1:5 ~ (Q67C - 1) / 4,
                                  Q67C %in% c(8, 9) ~ 0.5, TRUE ~ NA_real_),
      influence_eu    = case_when(Q67H %in% 1:5 ~ (Q67H - 1) / 4,
                                  Q67H %in% c(8, 9) ~ 0.5, TRUE ~ NA_real_),
      tertiary = case_when(Q96 %in% c(98, 99, -1) ~ NA_real_,
                           Q96 >= 6 ~ 1, TRUE ~ 0),
      female = case_when(Q101 == 2 ~ 1, Q101 == 1 ~ 0, TRUE ~ NA_real_),
      age    = ifelse(Q1 %in% c(998, 999, -1), NA, Q1),
      urban  = case_when(URBRUR == 1 ~ 1, URBRUR == 2 ~ 0, TRUE ~ NA_real_),
      lived_pov  = LivedPoverty,
      news_count = rowSums(cbind(Q65A %in% 1:4, Q65B %in% 1:4, Q65C %in% 1:4,
                                 Q65D %in% 1:4, Q65E %in% 1:4), na.rm = FALSE),
      w = Combinwt_new_hh
    )
  
  cat("\nNon-aware N:", nrow(df_non_aware),
      "| Missing weights:", sum(is.na(df_non_aware$w)), "\n")
  
  stack_non_aware <- df_non_aware %>%
    select(resp_id, country, w, tertiary, female, age, urban, lived_pov,
           news_count, influence_china, influence_us, influence_eu) %>%
    pivot_longer(cols = c(influence_china, influence_us, influence_eu),
                 names_to = "power", values_to = "rating") %>%
    filter(!is.na(rating)) %>%
    mutate(power = factor(power,
                          levels = c("influence_china", "influence_us", "influence_eu"),
                          labels = c("China", "US", "EU")))
  
  stack_pooled <- bind_rows(
    stack_non_aware %>% mutate(sample = "Non-aware"),
    stack_geo      %>% mutate(sample = "Climate-aware")
  ) %>%
    mutate(sample = factor(sample, levels = c("Non-aware", "Climate-aware")))
  
  # --- Weighted pooled model ---
  m_pooled_w <- lm(
    rating ~ power * sample + tertiary + female + age + urban + lived_pov +
      news_count + country,
    data = stack_pooled, weights = w
  )
  ct_pooled_w <- coeftest(m_pooled_w, vcov = vcovCL(m_pooled_w,
                                                    cluster = ~ resp_id, data = stack_pooled))
  
  # --- Unweighted pooled model ---
  m_pooled_u <- lm(
    rating ~ power * sample + tertiary + female + age + urban + lived_pov +
      news_count + country,
    data = stack_pooled
  )
  ct_pooled_u <- coeftest(m_pooled_u, vcov = vcovCL(m_pooled_u,
                                                    cluster = ~ resp_id, data = stack_pooled))
  
  cat("\n--- Table 4 (WEIGHTED): Climate Awareness Comparison ---\n")
  print(ct_pooled_w[grep("power|sample", rownames(ct_pooled_w)), ])
  
  cat("\n--- Table 4 (UNWEIGHTED, for comparison) ---\n")
  print(ct_pooled_u[grep("power|sample", rownames(ct_pooled_u)), ])
  
  # --- Baseline and climate-aware totals (weighted) ---
  base_w <- avg_comparisons(
    m_pooled_w,
    variables = list(power = c("China", "US")),
    newdata = datagrid(power = c("China", "US"),
                       sample = "Non-aware",
                       grid_type = "counterfactual"),
    vcov = vcovCL(m_pooled_w, cluster = ~ resp_id, data = stack_pooled),
    wts = "w"
  )
  aware_w <- avg_comparisons(
    m_pooled_w,
    variables = list(power = c("China", "US")),
    newdata = datagrid(power = c("China", "US"),
                       sample = "Climate-aware",
                       grid_type = "counterfactual"),
    vcov = vcovCL(m_pooled_w, cluster = ~ resp_id, data = stack_pooled),
    wts = "w"
  )
  cat("\n--- Baseline US-China Gap (Non-aware, WEIGHTED) ---\n")
  print(base_w)
  cat("\n--- Total US-China Gap (Climate-aware, WEIGHTED) ---\n")
  print(aware_w)
  
  # --- EU placebo, weighted ---
  cat("\n--- EU Placebo by Sample (WEIGHTED) ---\n")
  eu_w <- avg_comparisons(
    m_pooled_w,
    variables = list(power = c("China", "EU")),
    by = "sample",
    vcov = vcovCL(m_pooled_w, cluster = ~ resp_id, data = stack_pooled),
    wts = "w"
  )
  print(eu_w)
  
  # ==============================================================================
  # 5. THREE-WAY COMPARISON (ANOVA), WEIGHTED
  # ==============================================================================
  
  m_simple_w <- lm(
    rating ~ power + tertiary + female + age + urban + lived_pov +
      news_count + country,
    data = stack_geo, weights = w
  )
  m_full_w <- m_typ_w
  
  m_simple_u <- lm(
    rating ~ power + tertiary + female + age + urban + lived_pov +
      news_count + country,
    data = stack_geo
  )
  
  cat("\n--- ANOVA (WEIGHTED): Does Attribution Type Add Explanatory Power? ---\n")
  print(anova(m_simple_w, m_full_w))
  
  cat("\n--- ANOVA (UNWEIGHTED, for comparison) ---\n")
  print(anova(m_simple_u, m_typ_u))
  
  # ==============================================================================
  # 6. ROBUSTNESS CHECK: COMPLETE CASE, WEIGHTED
  # ==============================================================================
  
  df_complete <- dat %>%
    filter(Q60A == 1) %>%
    mutate(
      resp_id = RESPNO,
      country = as_factor(COUNTRY),
      cause_foreign  = case_when(Q61B == 2 ~ 1, Q61B == 1 ~ 0, TRUE ~ NA_real_),
      remedy_foreign = case_when(Q62D == 4 ~ 1,
                                 Q62D %in% c(1, 2, 3, 5) ~ 0, TRUE ~ NA_real_),
      influence_china = case_when(Q67B %in% 1:5 ~ (Q67B - 1) / 4, TRUE ~ NA_real_),
      influence_us    = case_when(Q67C %in% 1:5 ~ (Q67C - 1) / 4, TRUE ~ NA_real_),
      influence_eu    = case_when(Q67H %in% 1:5 ~ (Q67H - 1) / 4, TRUE ~ NA_real_),
      tertiary = case_when(Q96 %in% c(98, 99, -1) ~ NA_real_,
                           Q96 >= 6 ~ 1, TRUE ~ 0),
      female = case_when(Q101 == 2 ~ 1, Q101 == 1 ~ 0, TRUE ~ NA_real_),
      age    = ifelse(Q1 %in% c(998, 999, -1), NA, Q1),
      urban  = case_when(URBRUR == 1 ~ 1, URBRUR == 2 ~ 0, TRUE ~ NA_real_),
      lived_pov  = LivedPoverty,
      news_count = rowSums(cbind(Q65A %in% 1:4, Q65B %in% 1:4, Q65C %in% 1:4,
                                 Q65D %in% 1:4, Q65E %in% 1:4), na.rm = FALSE),
      w = Combinwt_new_hh,
      attribution_type = case_when(
        cause_foreign == 0 & remedy_foreign == 0 ~ "Consistent internaliser",
        cause_foreign == 1 & remedy_foreign == 1 ~ "Consistent externaliser",
        cause_foreign == 1 & remedy_foreign == 0 ~ "Sovereign adapter",
        cause_foreign == 0 & remedy_foreign == 1 ~ "Inconsistent externaliser",
        TRUE ~ NA_character_
      ),
      attribution_type = factor(attribution_type,
                                levels = c("Consistent internaliser",
                                           "Sovereign adapter",
                                           "Consistent externaliser",
                                           "Inconsistent externaliser"))
    ) %>%
    filter(!is.na(attribution_type),
           !is.na(influence_china),
           !is.na(influence_us),
           !is.na(influence_eu))
  
  cat("\n--- Robustness A: Complete Case Sample Size ---\n")
  cat("N =", nrow(df_complete), "| Missing weights:",
      sum(is.na(df_complete$w)), "\n")
  
  stack_complete <- df_complete %>%
    select(resp_id, country, attribution_type, w, tertiary, female, age, urban,
           lived_pov, news_count, influence_china, influence_us, influence_eu) %>%
    pivot_longer(cols = c(influence_china, influence_us, influence_eu),
                 names_to = "power", values_to = "rating") %>%
    filter(!is.na(rating)) %>%
    mutate(power = factor(power,
                          levels = c("influence_china", "influence_us", "influence_eu"),
                          labels = c("China", "US", "EU")))
  
  m_complete_w <- lm(
    rating ~ power * attribution_type + tertiary + female + age + urban +
      lived_pov + news_count + country,
    data = stack_complete, weights = w
  )
  ct_complete_w <- coeftest(m_complete_w,
                            vcov = vcovCL(m_complete_w, cluster = ~ resp_id,
                                          data = stack_complete))
  
  cat("\n--- Robustness A (WEIGHTED, complete case): Power Interactions ---\n")
  print(ct_complete_w[grep("power", rownames(ct_complete_w)), ])
  
  # ==============================================================================
  # 7. SIDE-BY-SIDE SUMMARY: UNWEIGHTED VS WEIGHTED
  # ==============================================================================
  
  cat("\n=== SIDE-BY-SIDE: Table 3 model (power terms) ===\n")
  cmp3 <- data.frame(
    term           = rownames(ct_typ_w),
    est_unweighted = ct_typ_u[, "Estimate"],
    est_weighted   = ct_typ_w[, "Estimate"],
    se_unweighted  = ct_typ_u[, "Std. Error"],
    se_weighted    = ct_typ_w[, "Std. Error"],
    p_unweighted   = ct_typ_u[, "Pr(>|t|)"],
    p_weighted     = ct_typ_w[, "Pr(>|t|)"]
  ) %>%
    filter(grepl("power", term)) %>%
    mutate(across(where(is.numeric), ~ round(.x, 4)))
  print(cmp3, row.names = FALSE)
  
  cat("\n=== SIDE-BY-SIDE: Table 4 model (power x sample terms) ===\n")
  cmp4 <- data.frame(
    term           = rownames(ct_pooled_w),
    est_unweighted = ct_pooled_u[, "Estimate"],
    est_weighted   = ct_pooled_w[, "Estimate"],
    se_unweighted  = ct_pooled_u[, "Std. Error"],
    se_weighted    = ct_pooled_w[, "Std. Error"],
    p_unweighted   = ct_pooled_u[, "Pr(>|t|)"],
    p_weighted     = ct_pooled_w[, "Pr(>|t|)"]
  ) %>%
    filter(grepl("power.*sample", term)) %>%
    mutate(across(where(is.numeric), ~ round(.x, 4)))
  print(cmp4, row.names = FALSE)
  
  cat("\n=== Done. ===\n")


# ==============================================================================
# VISUALIZATION CODE FOR FIGURES 1 & 2 (WEIGHTED)
# ==============================================================================

library(ggplot2)
library(dplyr)
library(sf)
library(rnaturalearth)
library(rnaturalearthdata)
library(marginaleffects)
  
  # ==============================================================================
  # FIGURE 1: Interaction Plot (Weighted)
  # ==============================================================================
  
  library(ggplot2)
  library(dplyr)
  library(marginaleffects)
  
  pred_typ <- avg_predictions(
    m_typ_w,
    variables = c("power", "attribution_type"),
    vcov = vcovCL(m_typ_w, cluster = ~ resp_id, data = stack_geo),
    wts = "w"
  )
  
  plot_df <- pred_typ %>%
    mutate(attribution_type = factor(attribution_type,
                                     levels = c("Consistent internaliser", "Sovereign adapter",
                                                "Consistent externaliser", "Inconsistent externaliser")))
  
  p_fig1 <- ggplot(plot_df, aes(x = power, y = estimate,
                                group = attribution_type,
                                color = attribution_type)) +
    geom_ribbon(aes(ymin = conf.low, ymax = conf.high, fill = attribution_type),
                alpha = 0.15, color = NA) +
    geom_line(linewidth = 0.8) +
    geom_point(size = 2) +
    scale_color_manual(values = c(
      "Consistent internaliser"   = "#F8766D",
      "Sovereign adapter"         = "#00BFC4",
      "Consistent externaliser"   = "#7CAE00",
      "Inconsistent externaliser" = "#C77CFF")) +
    scale_fill_manual(values = c(
      "Consistent internaliser"   = "#F8766D",
      "Sovereign adapter"         = "#00BFC4",
      "Consistent externaliser"   = "#7CAE00",
      "Inconsistent externaliser" = "#C77CFF")) +
    labs(x = "Global power",
         y = "Predicted rating (0–1 scale)",
         color = "Attribution type",
         fill = "Attribution type") +
    theme_classic(base_size = 10) +
    theme(legend.position = "right",
          legend.title = element_text(size = 7),
          legend.text = element_text(size = 7),
          axis.text = element_text(size = 7),
          axis.title = element_text(size = 7))
  
  print(p_fig1)
  ggsave("Figure1.pdf", p_fig1, width = 7, height = 4.5)
  
  
  # ==============================================================================
  # FIGURE 2: Map of the US-China Gap (Weighted)
  # ==============================================================================
  
  library(sf)
  library(rnaturalearth)
  library(rnaturalearthdata)
  library(grid)
  
  map_data <- df_lit %>%
    group_by(country) %>%
    summarize(
      mean_us    = weighted.mean(influence_us,    w = w, na.rm = TRUE),
      mean_china = weighted.mean(influence_china, w = w, na.rm = TRUE),
      gap        = mean_us - mean_china,
      .groups = "drop")
  
  africa <- ne_countries(scale = "medium", continent = "Africa", returnclass = "sf") %>%
    mutate(country = case_when(
      name == "Ivory Coast" ~ "Côte d'Ivoire",
      name == "Cape Verde" ~ "Cabo Verde",
      name == "United Republic of Tanzania" ~ "Tanzania",
      name == "Republic of the Congo" ~ "Congo-Brazzaville",
      name == "Democratic Republic of the Congo" ~ "Congo-Kinshasa",
      name == "São Tomé and Principe" ~ "São Tomé and Príncipe",
      name == "Gambia" ~ "Gambia, The",
      name == "Swaziland" ~ "Eswatini",
      TRUE ~ name))
  
  map_df <- africa %>% left_join(map_data, by = "country")
  
  p_fig2 <- ggplot(map_df) +
    geom_sf(aes(fill = gap), color = "#555555", linewidth = 0.2) +
    geom_sf(data = st_union(africa), fill = NA, color = "black", linewidth = 0.8) +
    scale_fill_gradient2(
      low = "#a50026", mid = "white", high = "#313695",
      midpoint = 0, na.value = "#e0e0e0",
      name = "US–China gap\n(0–1 scale)",
      labels = scales::number_format(accuracy = 0.01)) +
    theme_void(base_size = 10) +
    theme(legend.position = "bottom",
          legend.direction = "horizontal",
          legend.title = element_text(size = 7),
          legend.text = element_text(size = 7),
          legend.key.width = unit(2.5, "cm"),
          legend.key.height = unit(0.3, "cm"))
  
  print(p_fig2)
  ggsave("Figure2.pdf", p_fig2, width = 6.5, height = 7)
  
  
  #######################
  
  # Figure 1 — PDF for submission, PNG for working copy
  ggsave("Figure1.pdf", p_fig1, width = 7, height = 4.5)
  ggsave("Figure1.png", p_fig1, width = 7, height = 4.5, dpi = 300)
  
  # Figure 2 — PDF for submission, PNG for working copy
  ggsave("Figure2.pdf", p_fig2, width = 6.5, height = 7)
  ggsave("Figure2.png", p_fig2, width = 6.5, height = 7, dpi = 300)
  
  cat("\n===== FIGURES EXPORTED =====\n")


# ==============================================================================
# FULL REPLICATION CODE (WEIGHTED, WITH LITERACY ROBUSTNESS)
# Who Should Fix the Climate?
# ==============================================================================

# --- 0. Setup & Libraries ---
rm(list = ls())
library(haven)
library(dplyr)
library(tidyr)
library(tibble)
library(sandwich)
library(lmtest)
library(marginaleffects)
library(ggplot2)
library(broom)
library(car)

# --- LOAD THE DATA ---
setwd("C:/Users/agbog/Downloads")
dat <- read_sav("AB_R10.Merge Dataset_25Feb26.38ctry.release.sav")
cat("Raw sample loaded:", nrow(dat), "observations.\n")

# ==============================================================================
# 1. DATA PREPARATION: CLIMATE-AWARE SAMPLE + LITERACY FLAG
# ==============================================================================

df_lit <- dat %>%
  mutate(
    country       = as_factor(COUNTRY),
    resp_id       = RESPNO,
    climate_aware = case_when(Q60A == 1 ~ 1, TRUE ~ 0),
    
    # --- Climate Literacy (Simpson et al. 2021) ---
    # Literate = aware AND says human activity is a cause (Q61A = 1 or 3)
    climate_literate = case_when(
      Q60A == 1 & Q61A %in% c(1, 3) ~ 1,
      Q60A == 1 & Q61A %in% c(2, 4, 5, 8, 9) ~ 0,
      TRUE ~ NA_real_
    ),
    
    # --- Attribution Typology ---
    cause_foreign = case_when(
      Q61B == 2 ~ 1,
      Q61B == 1 ~ 0,
      TRUE ~ NA_real_
    ),
    remedy_foreign = case_when(
      Q62D == 4 ~ 1,
      Q62D %in% c(1, 2, 3, 5) ~ 0,
      TRUE ~ NA_real_
    ),
    
    # --- Evaluations of Influence (Rescaled 0-1) ---
    influence_china = case_when(
      Q67B %in% 1:5 ~ (Q67B - 1) / 4,
      Q67B %in% c(8, 9) ~ 0.5,
      TRUE ~ NA_real_
    ),
    influence_us = case_when(
      Q67C %in% 1:5 ~ (Q67C - 1) / 4,
      Q67C %in% c(8, 9) ~ 0.5,
      TRUE ~ NA_real_
    ),
    influence_eu = case_when(
      Q67H %in% 1:5 ~ (Q67H - 1) / 4,
      Q67H %in% c(8, 9) ~ 0.5,
      TRUE ~ NA_real_
    ),
    
    # --- Controls ---
    tertiary = case_when(
      Q96 %in% c(98, 99, -1) ~ NA_real_,
      Q96 >= 6 ~ 1,
      TRUE ~ 0
    ),
    female = case_when(Q101 == 2 ~ 1, Q101 == 1 ~ 0, TRUE ~ NA_real_),
    age    = ifelse(Q1 %in% c(998, 999, -1), NA, Q1),
    urban  = case_when(URBRUR == 1 ~ 1, URBRUR == 2 ~ 0, TRUE ~ NA_real_),
    lived_pov  = LivedPoverty,
    news_count = rowSums(cbind(
      Q65A %in% 1:4, Q65B %in% 1:4, Q65C %in% 1:4,
      Q65D %in% 1:4, Q65E %in% 1:4), na.rm = FALSE),
    
    # --- Typology ---
    attribution_type = case_when(
      cause_foreign == 0 & remedy_foreign == 0 ~ "Consistent internaliser",
      cause_foreign == 1 & remedy_foreign == 1 ~ "Consistent externaliser",
      cause_foreign == 1 & remedy_foreign == 0 ~ "Sovereign adapter",
      cause_foreign == 0 & remedy_foreign == 1 ~ "Inconsistent externaliser",
      TRUE ~ NA_character_
    ),
    attribution_type = factor(attribution_type,
                              levels = c("Consistent internaliser",
                                         "Sovereign adapter",
                                         "Consistent externaliser",
                                         "Inconsistent externaliser")),
    w = Combinwt_new_hh
  ) %>%
  filter(climate_aware == 1, !is.na(attribution_type))

cat("\nAnalytic Sample Size (Climate-Aware):", nrow(df_lit), "\n")
cat("Missing weights:", sum(is.na(df_lit$w)), "\n")
cat("Climate-literate subsample size:",
    sum(df_lit$climate_literate == 1, na.rm = TRUE), "\n")
cat("Climate-literate share of aware sample:",
    round(mean(df_lit$climate_literate == 1, na.rm = TRUE) * 100, 1), "%\n")

# ==============================================================================
# 2. DESCRIPTIVE STATISTICS
# ==============================================================================

cat("\n--- Table 1: Distribution of Attribution Frameworks (aware sample) ---\n")
tab_unw <- prop.table(table(df_lit$attribution_type)) * 100
tab_w   <- prop.table(xtabs(w ~ attribution_type, data = df_lit)) * 100
print(round(cbind(Unweighted = tab_unw, Weighted = tab_w), 3))

cat("\n--- Table 1b: Distribution among climate-literate subsample ---\n")
df_literate <- df_lit %>% filter(climate_literate == 1)
tab_unw_lit <- prop.table(table(df_literate$attribution_type)) * 100
tab_w_lit   <- prop.table(xtabs(w ~ attribution_type, data = df_literate)) * 100
print(round(cbind(Unweighted = tab_unw_lit, Weighted = tab_w_lit), 3))

cat("\n--- Cause vs Remedy Asymmetry ---\n")
cat("Blame Foreign Cause (aware, weighted):    ",
    round(weighted.mean(df_lit$cause_foreign, df_lit$w, na.rm = TRUE) * 100, 1), "%\n")
cat("Expect Foreign Remedy (aware, weighted):  ",
    round(weighted.mean(df_lit$remedy_foreign, df_lit$w, na.rm = TRUE) * 100, 1), "%\n")
cat("Blame Foreign Cause (literate, weighted): ",
    round(weighted.mean(df_literate$cause_foreign, df_literate$w, na.rm = TRUE) * 100, 1), "%\n")
cat("Expect Foreign Remedy (literate, weighted):",
    round(weighted.mean(df_literate$remedy_foreign, df_literate$w, na.rm = TRUE) * 100, 1), "%\n")

# Long format
stack_geo <- df_lit %>%
  select(resp_id, country, attribution_type, w, tertiary, female, age,
         urban, lived_pov, news_count, influence_china, influence_us, influence_eu,
         climate_literate) %>%
  pivot_longer(cols = c(influence_china, influence_us, influence_eu),
               names_to = "power", values_to = "rating") %>%
  filter(!is.na(rating)) %>%
  mutate(power = factor(power,
                        levels = c("influence_china", "influence_us", "influence_eu"),
                        labels = c("China", "US", "EU")))

cat("\n--- Table 2: Weighted Mean Ratings by Power and Type (aware sample) ---\n")
tab2_w <- stack_geo %>%
  group_by(attribution_type, power) %>%
  summarise(mean_w = weighted.mean(rating, w = w, na.rm = TRUE), .groups = "drop") %>%
  pivot_wider(names_from = power, values_from = mean_w) %>%
  column_to_rownames("attribution_type") %>%
  as.matrix()
print(round(tab2_w, 3))

cat("\n--- Table 2b: Weighted Mean Ratings (climate-literate subsample) ---\n")
tab2_w_lit <- stack_geo %>%
  filter(climate_literate == 1) %>%
  group_by(attribution_type, power) %>%
  summarise(mean_w = weighted.mean(rating, w = w, na.rm = TRUE), .groups = "drop") %>%
  pivot_wider(names_from = power, values_from = mean_w) %>%
  column_to_rownames("attribution_type") %>%
  as.matrix()
print(round(tab2_w_lit, 3))

# ==============================================================================
# 3. MAIN MODELS: CONDITIONAL PENALTY (AWARE SAMPLE)
# ==============================================================================

m_typ_w <- lm(
  rating ~ power * attribution_type + tertiary + female + age + urban +
    lived_pov + news_count + country,
  data = stack_geo, weights = w
)
ct_typ_w <- coeftest(m_typ_w, vcov = vcovCL(m_typ_w, cluster = ~ resp_id,
                                            data = stack_geo))

m_typ_u <- lm(
  rating ~ power * attribution_type + tertiary + female + age + urban +
    lived_pov + news_count + country,
  data = stack_geo
)
ct_typ_u <- coeftest(m_typ_u, vcov = vcovCL(m_typ_u, cluster = ~ resp_id,
                                            data = stack_geo))

cat("\n--- Table 3 (WEIGHTED): US-China Gap by Attribution Type ---\n")
print(ct_typ_w[grep("power", rownames(ct_typ_w)), ])

cat("\n--- Table 3 (UNWEIGHTED, comparison) ---\n")
print(ct_typ_u[grep("power", rownames(ct_typ_u)), ])

# Marginal means (weighted)
mm_w <- avg_predictions(
  m_typ_w,
  variables = c("power", "attribution_type"),
  vcov = vcovCL(m_typ_w, cluster = ~ resp_id, data = stack_geo),
  wts = "w"
)
cat("\n--- Marginal Means (WEIGHTED) ---\n")
print(mm_w)

# H2: explicit US-China gap by type
h2_w <- avg_comparisons(
  m_typ_w,
  variables = list(power = c("China", "US")),
  by = "attribution_type",
  vcov = vcovCL(m_typ_w, cluster = ~ resp_id, data = stack_geo),
  wts = "w"
)
cat("\n--- H2: US-China Gap by Type (Weighted) ---\n")
print(h2_w)

# ==============================================================================
# 4. CLIMATE AWARENESS COMPARISON (AWARE vs NON-AWARE)
# ==============================================================================

df_non_aware <- dat %>%
  filter(Q60A == 0) %>%
  mutate(
    resp_id = RESPNO,
    country = as_factor(COUNTRY),
    influence_china = case_when(Q67B %in% 1:5 ~ (Q67B - 1) / 4,
                                Q67B %in% c(8, 9) ~ 0.5, TRUE ~ NA_real_),
    influence_us    = case_when(Q67C %in% 1:5 ~ (Q67C - 1) / 4,
                                Q67C %in% c(8, 9) ~ 0.5, TRUE ~ NA_real_),
    influence_eu    = case_when(Q67H %in% 1:5 ~ (Q67H - 1) / 4,
                                Q67H %in% c(8, 9) ~ 0.5, TRUE ~ NA_real_),
    tertiary = case_when(Q96 %in% c(98, 99, -1) ~ NA_real_,
                         Q96 >= 6 ~ 1, TRUE ~ 0),
    female = case_when(Q101 == 2 ~ 1, Q101 == 1 ~ 0, TRUE ~ NA_real_),
    age    = ifelse(Q1 %in% c(998, 999, -1), NA, Q1),
    urban  = case_when(URBRUR == 1 ~ 1, URBRUR == 2 ~ 0, TRUE ~ NA_real_),
    lived_pov  = LivedPoverty,
    news_count = rowSums(cbind(Q65A %in% 1:4, Q65B %in% 1:4, Q65C %in% 1:4,
                               Q65D %in% 1:4, Q65E %in% 1:4), na.rm = FALSE),
    w = Combinwt_new_hh
  )

cat("\nNon-aware N:", nrow(df_non_aware),
    "| Missing weights:", sum(is.na(df_non_aware$w)), "\n")

stack_non_aware <- df_non_aware %>%
  select(resp_id, country, w, tertiary, female, age, urban, lived_pov,
         news_count, influence_china, influence_us, influence_eu) %>%
  pivot_longer(cols = c(influence_china, influence_us, influence_eu),
               names_to = "power", values_to = "rating") %>%
  filter(!is.na(rating)) %>%
  mutate(power = factor(power,
                        levels = c("influence_china", "influence_us", "influence_eu"),
                        labels = c("China", "US", "EU")))

stack_pooled <- bind_rows(
  stack_non_aware %>% mutate(sample = "Non-aware"),
  stack_geo      %>% mutate(sample = "Climate-aware")
) %>%
  mutate(sample = factor(sample, levels = c("Non-aware", "Climate-aware")))

m_pooled_w <- lm(
  rating ~ power * sample + tertiary + female + age + urban + lived_pov +
    news_count + country,
  data = stack_pooled, weights = w
)
ct_pooled_w <- coeftest(m_pooled_w, vcov = vcovCL(m_pooled_w,
                                                  cluster = ~ resp_id, data = stack_pooled))

cat("\n--- Table 5 (WEIGHTED): Climate Awareness Comparison ---\n")
print(ct_pooled_w[grep("power|sample", rownames(ct_pooled_w)), ])

base_w <- avg_comparisons(
  m_pooled_w,
  variables = list(power = c("China", "US")),
  newdata = datagrid(power = c("China", "US"),
                     sample = "Non-aware",
                     grid_type = "counterfactual"),
  vcov = vcovCL(m_pooled_w, cluster = ~ resp_id, data = stack_pooled),
  wts = "w"
)
aware_w <- avg_comparisons(
  m_pooled_w,
  variables = list(power = c("China", "US")),
  newdata = datagrid(power = c("China", "US"),
                     sample = "Climate-aware",
                     grid_type = "counterfactual"),
  vcov = vcovCL(m_pooled_w, cluster = ~ resp_id, data = stack_pooled),
  wts = "w"
)
cat("\n--- Baseline US-China Gap (Non-aware, WEIGHTED) ---\n")
print(base_w)
cat("\n--- Total US-China Gap (Climate-aware, WEIGHTED) ---\n")
print(aware_w)

cat("\n--- EU Placebo by Sample (WEIGHTED) ---\n")
eu_w <- avg_comparisons(
  m_pooled_w,
  variables = list(power = c("China", "EU")),
  by = "sample",
  vcov = vcovCL(m_pooled_w, cluster = ~ resp_id, data = stack_pooled),
  wts = "w"
)
print(eu_w)

# ==============================================================================
# 5. LITERACY ROBUSTNESS CHECK
# ==============================================================================

cat("\n=== LITERACY ROBUSTNESS CHECK ===\n")

# Re-estimate Table 3 model on the literate subsample
m_lit_w <- lm(
  rating ~ power * attribution_type + tertiary + female + age + urban +
    lived_pov + news_count + country,
  data = stack_geo %>% filter(climate_literate == 1),
  weights = w
)
ct_lit_w <- coeftest(m_lit_w, vcov = vcovCL(m_lit_w, cluster = ~ resp_id,
                                            data = stack_geo %>% filter(climate_literate == 1)))

cat("\n--- Table 3b (WEIGHTED, literate subsample): Power Interactions ---\n")
print(ct_lit_w[grep("power", rownames(ct_lit_w)), ])

# Re-estimate Table 5 model with literate as the "treated" sample
# (pooled: non-aware vs literate)
stack_pooled_lit <- bind_rows(
  stack_non_aware %>% mutate(sample = "Non-aware"),
  stack_geo %>% filter(climate_literate == 1) %>% mutate(sample = "Climate-literate")
) %>%
  mutate(sample = factor(sample, levels = c("Non-aware", "Climate-literate")))

m_pooled_lit_w <- lm(
  rating ~ power * sample + tertiary + female + age + urban + lived_pov +
    news_count + country,
  data = stack_pooled_lit, weights = w
)
ct_pooled_lit_w <- coeftest(m_pooled_lit_w, vcov = vcovCL(m_pooled_lit_w,
                                                          cluster = ~ resp_id, data = stack_pooled_lit))

cat("\n--- Table 5b (WEIGHTED, literate vs non-aware): Power Interactions ---\n")
print(ct_pooled_lit_w[grep("power|sample", rownames(ct_pooled_lit_w)), ])

# ==============================================================================
# 6. THREE-WAY COMPARISON (ANOVA), WEIGHTED
# ==============================================================================

m_simple_w <- lm(
  rating ~ power + tertiary + female + age + urban + lived_pov +
    news_count + country,
  data = stack_geo, weights = w
)
m_full_w <- m_typ_w

cat("\n--- ANOVA (WEIGHTED): Does Attribution Type Add Explanatory Power? ---\n")
print(anova(m_simple_w, m_full_w))

# ==============================================================================
# 7. ROBUSTNESS CHECK: COMPLETE CASE, WEIGHTED
# ==============================================================================

df_complete <- dat %>%
  filter(Q60A == 1) %>%
  mutate(
    resp_id = RESPNO,
    country = as_factor(COUNTRY),
    cause_foreign  = case_when(Q61B == 2 ~ 1, Q61B == 1 ~ 0, TRUE ~ NA_real_),
    remedy_foreign = case_when(Q62D == 4 ~ 1,
                               Q62D %in% c(1, 2, 3, 5) ~ 0, TRUE ~ NA_real_),
    influence_china = case_when(Q67B %in% 1:5 ~ (Q67B - 1) / 4, TRUE ~ NA_real_),
    influence_us    = case_when(Q67C %in% 1:5 ~ (Q67C - 1) / 4, TRUE ~ NA_real_),
    influence_eu    = case_when(Q67H %in% 1:5 ~ (Q67H - 1) / 4, TRUE ~ NA_real_),
    tertiary = case_when(Q96 %in% c(98, 99, -1) ~ NA_real_,
                         Q96 >= 6 ~ 1, TRUE ~ 0),
    female = case_when(Q101 == 2 ~ 1, Q101 == 1 ~ 0, TRUE ~ NA_real_),
    age    = ifelse(Q1 %in% c(998, 999, -1), NA, Q1),
    urban  = case_when(URBRUR == 1 ~ 1, URBRUR == 2 ~ 0, TRUE ~ NA_real_),
    lived_pov  = LivedPoverty,
    news_count = rowSums(cbind(Q65A %in% 1:4, Q65B %in% 1:4, Q65C %in% 1:4,
                               Q65D %in% 1:4, Q65E %in% 1:4), na.rm = FALSE),
    w = Combinwt_new_hh,
    attribution_type = case_when(
      cause_foreign == 0 & remedy_foreign == 0 ~ "Consistent internaliser",
      cause_foreign == 1 & remedy_foreign == 1 ~ "Consistent externaliser",
      cause_foreign == 1 & remedy_foreign == 0 ~ "Sovereign adapter",
      cause_foreign == 0 & remedy_foreign == 1 ~ "Inconsistent externaliser",
      TRUE ~ NA_character_
    ),
    attribution_type = factor(attribution_type,
                              levels = c("Consistent internaliser",
                                         "Sovereign adapter",
                                         "Consistent externaliser",
                                         "Inconsistent externaliser"))
  ) %>%
  filter(!is.na(attribution_type),
         !is.na(influence_china),
         !is.na(influence_us),
         !is.na(influence_eu))

cat("\n--- Robustness A: Complete Case Sample Size ---\n")
cat("N =", nrow(df_complete), "| Missing weights:",
    sum(is.na(df_complete$w)), "\n")

stack_complete <- df_complete %>%
  select(resp_id, country, attribution_type, w, tertiary, female, age, urban,
         lived_pov, news_count, influence_china, influence_us, influence_eu) %>%
  pivot_longer(cols = c(influence_china, influence_us, influence_eu),
               names_to = "power", values_to = "rating") %>%
  filter(!is.na(rating)) %>%
  mutate(power = factor(power,
                        levels = c("influence_china", "influence_us", "influence_eu"),
                        labels = c("China", "US", "EU")))

m_complete_w <- lm(
  rating ~ power * attribution_type + tertiary + female + age + urban +
    lived_pov + news_count + country,
  data = stack_complete, weights = w
)
ct_complete_w <- coeftest(m_complete_w,
                          vcov = vcovCL(m_complete_w, cluster = ~ resp_id,
                                        data = stack_complete))

cat("\n--- Robustness A (WEIGHTED, complete case): Power Interactions ---\n")
print(ct_complete_w[grep("power", rownames(ct_complete_w)), ])

# ==============================================================================
# 8. SUMMARY: AWARE vs LITERATE COMPARISON
# ==============================================================================

cat("\n=== SUMMARY: Table 3 power terms, aware vs literate ===\n")

power_terms <- grep("power", rownames(ct_typ_w))

cmp <- data.frame(
  term         = rownames(ct_typ_w)[power_terms],
  aware_est    = ct_typ_w[power_terms, "Estimate"],
  aware_p      = ct_typ_w[power_terms, "Pr(>|t|)"],
  literate_est = ct_lit_w[power_terms, "Estimate"],
  literate_p   = ct_lit_w[power_terms, "Pr(>|t|)"],
  stringsAsFactors = FALSE
) %>%
  mutate(across(where(is.numeric), ~ round(.x, 4)))

print(cmp, row.names = FALSE)