#!/usr/bin/env Rscript

options(stringsAsFactors = FALSE)

required_packages <- c("ggplot2", "sandwich", "lme4")
missing_packages <- required_packages[!vapply(required_packages, requireNamespace, logical(1), quietly = TRUE)]
if (length(missing_packages)) stop("Missing R packages: ", paste(missing_packages, collapse = ", "))

library(ggplot2)

script_args <- commandArgs(trailingOnly = FALSE)
file_arg <- script_args[grepl("^--file=", script_args)][1]
script_path <- normalizePath(sub("^--file=", "", file_arg))
repo_root <- normalizePath(file.path(dirname(script_path), ".."))
setwd(repo_root)

output_dir <- file.path("revision", "outputs")
figure_dir <- file.path("revision", "figures")
dir.create(output_dir, recursive = TRUE, showWarnings = FALSE)
dir.create(figure_dir, recursive = TRUE, showWarnings = FALSE)

panel_file <- file.path("data", "releases", "party_2d_election_year_panel_v0.csv.gz")
annual_file <- file.path("data", "releases", "party_2d_annual_model_output_v0.csv.gz")
no_vparty_file <- "/projects/party4d_validation_runs/no_vparty/outputs/estimations/latest/party_positions_2026-06-05_14-57-17.csv"

for (f in c(panel_file, annual_file, no_vparty_file)) {
  if (!file.exists(f)) stop("Required input not found: ", f)
}

panel <- read.csv(panel_file, check.names = FALSE)
annual <- read.csv(annual_file, check.names = FALSE)
no_vparty <- read.csv(no_vparty_file, check.names = FALSE)
families <- read.csv(file.path("data", "party_families.csv"), check.names = FALSE)
names(families)[names(families) == "partyfacts_id"] <- "party_id"

iso_names <- c(
  AL="Albania", AM="Armenia", AR="Argentina", AT="Austria", AU="Australia",
  AZ="Azerbaijan", BA="Bosnia and Herzegovina", BE="Belgium", BG="Bulgaria",
  BO="Bolivia", BR="Brazil", BY="Belarus", CA="Canada", CH="Switzerland",
  CL="Chile", CO="Colombia", CR="Costa Rica", CY="Cyprus", CZ="Czechia",
  DE="Germany", DK="Denmark", DO="Dominican Republic", EC="Ecuador",
  EE="Estonia", ES="Spain", FI="Finland", FR="France", GB="United Kingdom",
  GE="Georgia", GR="Greece", HR="Croatia", HU="Hungary", IE="Ireland",
  IL="Israel", IS="Iceland", IT="Italy", JP="Japan", KR="South Korea",
  LK="Sri Lanka", LT="Lithuania", LU="Luxembourg", LV="Latvia", MD="Moldova",
  ME="Montenegro", MK="North Macedonia", MT="Malta", MX="Mexico",
  NL="Netherlands", NO="Norway", NZ="New Zealand", PA="Panama", PE="Peru",
  PL="Poland", PT="Portugal", RO="Romania", RS="Serbia", RU="Russia",
  SE="Sweden", SI="Slovenia", SK="Slovakia", TR="Türkiye", UA="Ukraine",
  US="United States", UY="Uruguay", ZA="South Africa"
)

region_for <- function(country) {
  europe <- c("AL","AT","BA","BE","BG","BY","CH","CY","CZ","DE","DK","EE","ES","FI","FR","GB","GE","GR","HR","HU","IE","IS","IT","LT","LU","LV","MD","ME","MK","MT","NL","NO","PL","PT","RO","RS","RU","SE","SI","SK","TR","UA")
  latin <- c("AR","BO","BR","CL","CO","CR","DO","EC","MX","PA","PE","UY")
  north_america <- c("CA", "US")
  asia_pacific <- c("AU", "JP", "KR", "LK", "NZ")
  out <- rep("Other", length(country))
  out[country %in% europe] <- "Europe"
  out[country %in% latin] <- "Latin America"
  out[country %in% north_america] <- "North America"
  out[country %in% asia_pacific] <- "Asia-Pacific"
  out
}

family_names <- c(
  com="Communist/Far Left", eco="Green/Ecological", soc="Social Democratic",
  chr="Christian Democratic", agr="Agrarian", lib="Liberal", con="Conservative",
  right="Radical Right", spec="Special issue", other="Other", nofam="Unclassified"
)

theme_revision <- function() {
  theme_minimal(base_size = 9) +
    theme(
      panel.grid.minor = element_blank(),
      strip.text = element_text(face = "bold"),
      legend.position = "bottom",
      plot.title.position = "plot"
    )
}

write_clean_csv <- function(x, path) {
  write.csv(x, path, row.names = FALSE, na = "")
  message("Wrote ", path, " (", nrow(x), " rows)")
}

# ---------------------------------------------------------------------------
# Versioned country coverage
# ---------------------------------------------------------------------------

country_rows <- lapply(sort(unique(panel$country)), function(cc) {
  d <- panel[panel$country == cc, ]
  tab <- table(d$source_support_class)
  get_n <- function(label) if (label %in% names(tab)) unname(tab[label]) else 0
  data.frame(
    release = "v0", iso2 = cc, country = unname(iso_names[cc]),
    first_year = min(d$year), last_year = max(d$year),
    parties = length(unique(d$party_id)), election_year_rows = nrow(d),
    both_text_expert = get_n("both_direct_or_nearby"),
    text_only = get_n("text_only_direct_or_nearby"),
    expert_only = get_n("expert_only_direct_or_nearby"),
    temporal_propagation = get_n("temporal_propagation")
  )
})
country_coverage <- do.call(rbind, country_rows)
write_clean_csv(country_coverage, file.path("metadata", "country_coverage_v0.csv"))

# ---------------------------------------------------------------------------
# Source-support balance with uncertainty and hierarchical heterogeneity
# ---------------------------------------------------------------------------

balance <- panel[is.na(panel$pervote) | panel$pervote > 0, ]
balance$pervote[is.na(balance$pervote)] <- 0
balance$decade <- factor(floor(balance$year / 10) * 10)
balance$region <- factor(region_for(balance$country))
balance$country <- factor(balance$country)
balance$text_only <- as.integer(balance$source_support_class == "text_only_direct_or_nearby")
balance$expert_only <- as.integer(balance$source_support_class == "expert_only_direct_or_nearby")
balance$temporal <- as.integer(balance$source_support_class == "temporal_propagation")

class_labels <- c(
  text_only="Text only", expert_only="Expert only", temporal="Temporal propagation"
)

pooled_rows <- list()
interaction_rows <- list()
row_i <- 1
interaction_i <- 1

for (dimension in c("economic_lr", "galtan")) {
  form <- as.formula(paste(dimension, "~ text_only + expert_only + temporal + factor(country) + decade + log1p(pervote)"))
  model <- lm(form, data = balance)
  cluster_vcov <- sandwich::vcovCL(model, cluster = balance$country, type = "HC1")
  for (term in names(class_labels)) {
    estimate <- unname(coef(model)[term])
    se <- sqrt(cluster_vcov[term, term])
    category <- switch(term,
      text_only="text_only_direct_or_nearby",
      expert_only="expert_only_direct_or_nearby",
      temporal="temporal_propagation")
    pooled_rows[[row_i]] <- data.frame(
      dimension = dimension, comparison = class_labels[term],
      category_n = sum(balance$source_support_class == category),
      model_n = nobs(model), estimate = estimate, clustered_se = se,
      ci_lower = estimate - 1.96 * se, ci_upper = estimate + 1.96 * se,
      p_value = 2 * pnorm(abs(estimate / se), lower.tail = FALSE)
    )
    row_i <- row_i + 1
  }

  # Partially pool the well-populated text-only contrast across region and
  # decade. Expert-only and temporal-propagation contrasts remain pooled
  # because their samples are too sparse for stable varying slopes.
  varying_form <- as.formula(paste(
    dimension,
    "~ expert_only + temporal + text_only + log1p(pervote) +",
    "(1 | country) + (0 + text_only | region) + (0 + text_only | decade)"
  ))
  varying_model <- lme4::lmer(varying_form, data = balance, REML = TRUE,
                              control = lme4::lmerControl(check.nobs.vs.nRE = "ignore"))
  fixed_text <- lme4::fixef(varying_model)["text_only"]
  fixed_var <- as.matrix(vcov(varying_model))["text_only", "text_only"]
  re <- lme4::ranef(varying_model, condVar = TRUE)

  for (group_type in c("region", "decade")) {
    group_values <- rownames(re[[group_type]])
    post_var <- attr(re[[group_type]], "postVar")
    for (group_index in seq_along(group_values)) {
      g <- group_values[group_index]
      template <- balance[as.character(balance[[group_type]]) == g, , drop = FALSE]
      estimate <- fixed_text + re[[group_type]][g, "text_only"]
      # Approximate conditional interval. Adding fixed-effect and conditional
      # random-effect variances is conservative because their covariance is
      # not exposed by ranef().
      se <- sqrt(fixed_var + post_var[1, 1, group_index])
      interaction_rows[[interaction_i]] <- data.frame(
        dimension = dimension, group_type = group_type, group = as.character(g),
        comparison = "Text only", group_n = nrow(template),
        text_only_n = sum(template$source_support_class == "text_only_direct_or_nearby"),
        estimate = unname(estimate), approximate_se = unname(se),
        ci_lower = unname(estimate - 1.96 * se), ci_upper = unname(estimate + 1.96 * se),
        model_n = nobs(varying_model), singular_fit = lme4::isSingular(varying_model)
      )
      interaction_i <- interaction_i + 1
    }
  }
}

pooled <- do.call(rbind, pooled_rows)
interactions <- do.call(rbind, interaction_rows)
write_clean_csv(pooled, file.path(output_dir, "source_support_pooled.csv"))
write_clean_csv(interactions, file.path(output_dir, "source_support_interactions.csv"))

interactions$dimension_label <- ifelse(interactions$dimension == "economic_lr", "Economic", "Cultural")
interactions$group_type_label <- ifelse(interactions$group_type == "region", "Region", "Decade")
p_balance <- ggplot(interactions, aes(x = estimate, y = reorder(group, estimate))) +
  geom_vline(xintercept = 0, colour = "grey40", linewidth = 0.4, linetype = 2) +
  geom_errorbar(aes(xmin = ci_lower, xmax = ci_upper), orientation = "y", width = 0, linewidth = 0.4, colour = "grey20") +
  geom_point(size = 1.6, colour = "black") +
  facet_wrap(~dimension_label + group_type_label, scales = "free_y", ncol = 2) +
  labs(x = "Adjusted difference from overlapping text-and-expert support", y = NULL) +
  theme_revision()
ggsave(file.path(figure_dir, "source_support_heterogeneity.pdf"), p_balance,
       width = 8.2, height = 6.8, device = cairo_pdf)

# ---------------------------------------------------------------------------
# V-Party ablation: case and subgroup sensitivity
# ---------------------------------------------------------------------------

keys <- c("party_id", "country", "year")
vparty <- merge(panel, no_vparty, by = keys, suffixes = c("_production", "_no_vparty"))
vparty <- merge(vparty, families, by = "party_id", all.x = TRUE)
vparty$family[is.na(vparty$family)] <- "nofam"
vparty$family_label <- unname(family_names[vparty$family])
vparty$region <- region_for(vparty$country)
vparty$decade <- floor(vparty$year / 10) * 10

for (dimension in c("economic_lr", "galtan")) {
  prod <- vparty[[paste0(dimension, "_production")]]
  ablated <- vparty[[paste0(dimension, "_no_vparty")]]
  vparty[[paste0(dimension, "_difference")]] <- ablated - prod
  vparty[[paste0(dimension, "_abs_difference")]] <- abs(ablated - prod)
  width_prod <- vparty[[paste0(dimension, "_q975_production")]] - vparty[[paste0(dimension, "_q025_production")]]
  width_ablation <- vparty[[paste0(dimension, "_q975_no_vparty")]] - vparty[[paste0(dimension, "_q025_no_vparty")]]
  vparty[[paste0(dimension, "_interval_width_production")]] <- width_prod
  vparty[[paste0(dimension, "_interval_width_no_vparty")]] <- width_ablation
  vparty[[paste0(dimension, "_interval_width_change")]] <- width_ablation - width_prod
}

summary_by <- function(data, group_var, dimension) {
  split_data <- split(data, data[[group_var]], drop = TRUE)
  rows <- lapply(names(split_data), function(g) {
    d <- split_data[[g]]
    absdiff <- d[[paste0(dimension, "_abs_difference")]]
    diff <- d[[paste0(dimension, "_difference")]]
    wprod <- d[[paste0(dimension, "_interval_width_production")]]
    wabl <- d[[paste0(dimension, "_interval_width_no_vparty")]]
    data.frame(
      dimension = dimension, group_type = group_var, group = g,
      n = nrow(d), parties = length(unique(d$party_id)),
      mean_abs_difference = mean(absdiff), median_abs_difference = median(absdiff),
      p90_abs_difference = unname(quantile(absdiff, .90)),
      p95_abs_difference = unname(quantile(absdiff, .95)),
      mean_signed_difference = mean(diff),
      mean_interval_width_production = mean(wprod),
      mean_interval_width_no_vparty = mean(wabl),
      mean_interval_width_change = mean(wabl - wprod)
    )
  })
  do.call(rbind, rows)
}

vparty_summaries <- do.call(rbind, lapply(c("country", "decade", "region", "family_label", "source_support_class"), function(g) {
  do.call(rbind, lapply(c("economic_lr", "galtan"), function(d) summary_by(vparty, g, d)))
}))

write_clean_csv(vparty, file.path(output_dir, "vparty_sensitivity_matched_rows.csv"))
write_clean_csv(vparty_summaries, file.path(output_dir, "vparty_sensitivity_groups.csv"))

top_cases <- do.call(rbind, lapply(c("economic_lr", "galtan"), function(dimension) {
  ord <- order(vparty[[paste0(dimension, "_abs_difference")]], decreasing = TRUE)
  d <- vparty[head(ord, 30), ]
  data.frame(
    dimension = dimension, party_id = d$party_id,
    party_name = d$party_name_english, party_short = d$party_name_short,
    country = d$country, year = d$year,
    production = d[[paste0(dimension, "_production")]],
    no_vparty = d[[paste0(dimension, "_no_vparty")]],
    signed_difference = d[[paste0(dimension, "_difference")]],
    absolute_difference = d[[paste0(dimension, "_abs_difference")]],
    production_interval_width = d[[paste0(dimension, "_interval_width_production")]],
    no_vparty_interval_width = d[[paste0(dimension, "_interval_width_no_vparty")]],
    source_support_class = d$source_support_class,
    family = d$family_label
  )
}))
write_clean_csv(top_cases, file.path(output_dir, "vparty_sensitivity_top_cases.csv"))

country_plot <- vparty_summaries[vparty_summaries$group_type == "country" & vparty_summaries$n >= 10, ]
country_plot <- do.call(rbind, lapply(split(country_plot, country_plot$dimension), function(d) {
  head(d[order(d$mean_abs_difference, decreasing = TRUE), ], 15)
}))
country_plot$dimension_label <- ifelse(country_plot$dimension == "economic_lr", "Economic", "Cultural")
p_vparty <- ggplot(country_plot, aes(x = mean_abs_difference, y = reorder(group, mean_abs_difference))) +
  geom_point(aes(size = n), colour = "black") +
  facet_wrap(~dimension_label, scales = "free_y") +
  labs(x = "Mean absolute change after removing V-Party", y = "Country", size = "Party-years") +
  theme_revision()
ggsave(file.path(figure_dir, "vparty_sensitivity_countries.pdf"), p_vparty,
       width = 7.2, height = 5.8, device = cairo_pdf)

# ---------------------------------------------------------------------------
# Illustrative trajectories and landmark cases
# ---------------------------------------------------------------------------

trajectory_ids <- c(383, 1567, 379, 409)
traj <- panel[panel$party_id %in% trajectory_ids, ]
traj$party_label <- factor(traj$party_id, levels = trajectory_ids,
  labels = c("Germany: Social Democrats", "United Kingdom: Conservatives", "Denmark: Social Democrats", "Sweden: Sweden Democrats"))
traj_long <- rbind(
  data.frame(traj[c("party_id","party_label","country","year","source_support_class")],
             dimension="Economic", estimate=traj$economic_lr,
             lower=traj$economic_lr_q025, upper=traj$economic_lr_q975),
  data.frame(traj[c("party_id","party_label","country","year","source_support_class")],
             dimension="Cultural", estimate=traj$galtan,
             lower=traj$galtan_q025, upper=traj$galtan_q975)
)
write_clean_csv(traj_long, file.path(output_dir, "party_trajectory_plot_data.csv"))
p_traj <- ggplot(traj_long, aes(x = year, y = estimate)) +
  geom_ribbon(aes(ymin = lower, ymax = upper), fill = "grey85", alpha = 0.5) +
  geom_line(colour = "black", linewidth = .65) +
  geom_point(colour = "black", size = .9) +
  facet_grid(dimension ~ party_label) +
  scale_y_continuous(limits = c(0, 1), breaks = c(0, .5, 1)) +
  labs(x = "Election year", y = "Posterior position (0–1)") +
  theme_revision() +
  theme(axis.text.x = element_text(angle = 45, hjust = 1))
ggsave(file.path(figure_dir, "party_trajectories.pdf"), p_traj,
       width = 10.2, height = 4.9, device = cairo_pdf)

landmark_spec <- data.frame(
  party_id = c(383,383,1375,1375,1516,1516,1567,1567,487,487,409,409,433,433,1545,432,809),
  target_year = c(1972,2021,1983,2021,1983,1997,1979,2019,1994,2022,2010,2022,1988,2022,2021,2020,2020),
  label = c("German Social Democrats 1972","German Social Democrats 2021","German Christian Democrats 1983","German Christian Democrats 2021","Labour 1983","Labour 1997",
            "Conservatives 1979","Conservatives 2019","Swedish Social Democrats 1994","Swedish Social Democrats 2022",
            "Sweden Democrats 2010","Sweden Democrats 2022","French National Front 1988","French National Front 2022",
            "The Left 2021","US Democrats 2020","US Republicans 2020")
)

landmark_rows <- lapply(seq_len(nrow(landmark_spec)), function(i) {
  spec <- landmark_spec[i, ]
  d <- panel[panel$party_id == spec$party_id & abs(panel$year - spec$target_year) <= 3, ]
  if (!nrow(d)) return(data.frame(spec, status="missing", actual_year=NA, economic_lr=NA, galtan=NA))
  d <- d[which.min(abs(d$year - spec$target_year)), ]
  data.frame(spec, status=ifelse(d$year == spec$target_year, "exact", "nearest within 3 years"),
             actual_year=d$year, party_name=d$party_name_english,
             country=d$country, economic_lr=d$economic_lr, economic_lower=d$economic_lr_q025,
             economic_upper=d$economic_lr_q975, galtan=d$galtan,
             galtan_lower=d$galtan_q025, galtan_upper=d$galtan_q975,
             source_support_class=d$source_support_class)
})
landmarks <- do.call(rbind, landmark_rows)
landmarks$era <- ifelse(landmarks$target_year < 2000, "Historical & Cold War Era (1970–1999)", "Contemporary Era (2000–2022)")
landmarks$era <- factor(landmarks$era, levels = c("Historical & Cold War Era (1970–1999)", "Contemporary Era (2000–2022)"))

write_clean_csv(landmarks, file.path(output_dir, "party_landmark_plot_data.csv"))
plot_landmarks <- landmarks[!is.na(landmarks$economic_lr), ]

label_offsets <- data.frame(
  label = landmark_spec$label,
  dx = c(-.060,-.060,-.055,.045,-.040,-.060,.035,-.055,-.005,.080,-.060,.060,.030,-.030,.035,.080,.025),
  dy = c(.070,.040,.055,-.035,-.070,-.040,-.025,.030,.040,.060,-.035,.035,.035,.050,.035,-.060,.035)
)
plot_landmarks <- merge(plot_landmarks, label_offsets, by = "label", all.x = TRUE, sort = FALSE)
plot_landmarks$label_x <- pmin(.98, pmax(.02, plot_landmarks$economic_lr + plot_landmarks$dx))
plot_landmarks$label_y <- pmin(.98, pmax(.02, plot_landmarks$galtan + plot_landmarks$dy))

# Country shapes for clear grayscale distinction
country_shapes <- c(DE = 16, GB = 17, SE = 15, FR = 18, US = 8)

p_land <- ggplot(plot_landmarks, aes(x = economic_lr, y = galtan)) +
  geom_errorbar(aes(ymin = galtan_lower, ymax = galtan_upper), width = 0, alpha = .35, colour = "grey30") +
  geom_errorbar(aes(xmin = economic_lower, xmax = economic_upper), orientation = "y", width = 0, alpha = .35, colour = "grey30") +
  geom_segment(aes(xend = label_x, yend = label_y), linewidth = .2, colour = "grey50") +
  geom_point(aes(shape = country), size = 2.2, fill = "black", colour = "black") +
  geom_text(aes(x = label_x, y = label_y, label = label), size = 2.4,
            check_overlap = TRUE, show.legend = FALSE) +
  facet_wrap(~era, ncol = 2) +
  scale_shape_manual(values = country_shapes) +
  scale_x_continuous(limits = c(0,1), breaks = c(0,.25,.5,.75,1)) +
  scale_y_continuous(limits = c(0,1), breaks = c(0,.25,.5,.75,1)) +
  labs(x = "Economic: left (0) to right (1)",
       y = "Cultural: cosmopolitan (0) to traditionalist (1)", shape = "Country") +
  theme_revision()
ggsave(file.path(figure_dir, "party_landmarks.pdf"), p_land,
       width = 9.2, height = 5.2, device = cairo_pdf)

# ---------------------------------------------------------------------------
# Predictive-interval calibration and concentration of misses
# ---------------------------------------------------------------------------

ppc_dir <- file.path(output_dir, "ppc")
ppc_item_file <- file.path(ppc_dir, "posterior_predictive_by_dimension_item.csv")
ppc_calibration_file <- file.path(ppc_dir, "posterior_predictive_calibration_curve.csv")
if (file.exists(ppc_item_file) && file.exists(ppc_calibration_file)) {
  ppc_items <- read.csv(ppc_item_file, check.names = FALSE)
  ppc_items$dimension <- ifelse(ppc_items$dim_idx == 1, "Economic", "Cultural")
  ppc_items$item_label <- c(
    gender_vparty="V-Party gender-policy item", relig_vparty="V-Party religion-policy item",
    welf_vparty="V-Party welfare-position item", culsup_vparty="V-Party cultural-superiority item",
    immig_vparty="V-Party immigration-policy item", lgbt_vparty="V-Party LGBT-policy item",
    lrecon_vparty="V-Party economic LR", galtan_ches="CHES GAL--TAN",
    lrecon_ches="CHES economic LR", libcon_gps="GPS cultural",
    lrecon_gps="GPS economic LR", lrecon_poppa="POPPA economic LR"
  )[ppc_items$var]
  ppc_items$item_label[is.na(ppc_items$item_label)] <- ppc_items$var[is.na(ppc_items$item_label)]
  ppc_items$item_label <- factor(ppc_items$item_label,
    levels = rev(ppc_items$item_label[order(ppc_items$observed_coverage)]))
  p_ppc <- ggplot(ppc_items, aes(x = observed_coverage, y = item_label)) +
    geom_vline(xintercept = .95, linetype = 2, colour = "grey40", linewidth = .35) +
    geom_errorbar(aes(xmin = coverage_ci_lower, xmax = coverage_ci_upper),
                  orientation = "y", width = 0, linewidth = .35, colour = "grey20") +
    geom_point(aes(size = n, shape = project), colour = "black", fill = "gray30") +
    facet_wrap(~dimension, scales = "free_y") +
    scale_x_continuous(limits = c(.70, 1), breaks = c(.70,.80,.90,.95,1),
                       labels = function(x) paste0(round(100*x), "%")) +
    labs(x = "Observed coverage of nominal 95% predictive interval", y = NULL,
         shape = "Source project", size = "Ratings") +
    theme_revision()
  ggsave(file.path(figure_dir, "predictive_coverage_by_item.pdf"), p_ppc,
         width = 7.8, height = 5.4, device = cairo_pdf)

  ppc_cal <- read.csv(ppc_calibration_file, check.names = FALSE)
  ppc_cal$dimension <- ifelse(ppc_cal$dim_idx == 1, "Economic", "Cultural")
  p_cal <- ggplot(ppc_cal, aes(x = nominal_level, y = observed_coverage, linetype = dimension, shape = dimension)) +
    geom_abline(slope = 1, intercept = 0, linetype = 3, colour = "grey50") +
    geom_line(linewidth = .6, colour = "black") + geom_point(size = 1.8, colour = "black") +
    scale_x_continuous(limits = c(.45, 1), breaks = c(.5,.8,.9,.95,1),
                       labels = function(x) paste0(round(100*x), "%")) +
    scale_y_continuous(limits = c(.4, 1), breaks = c(.4,.5,.6,.7,.8,.9,1),
                       labels = function(x) paste0(round(100*x), "%")) +
    labs(x = "Nominal interval level", y = "Observed coverage", linetype = "Dimension", shape = "Dimension") +
    coord_equal() + theme_revision()
  ggsave(file.path(figure_dir, "predictive_calibration_curve.pdf"), p_cal,
         width = 5.7, height = 5.0, device = cairo_pdf)


  ppc_obs_file <- file.path(ppc_dir, "posterior_predictive_observations.csv")
  if (file.exists(ppc_obs_file)) {
    ppc_obs <- read.csv(ppc_obs_file, check.names = FALSE)
    ppc_obs$covered_95 <- tolower(as.character(ppc_obs$covered_95)) == "true"
    expert_input <- read.csv(file.path("data", "expert.csv"), check.names = FALSE)
    key <- function(d) paste(d$party, d$country, d$year, d$var, sprintf("%.10f", d$val), sep = "|")
    ppc_obs$n_experts <- expert_input$n_experts[match(key(ppc_obs), key(expert_input))]
    ppc_obs$interval_width <- ppc_obs$pred_upper - ppc_obs$pred_lower
    ppc_obs$fitted_bin <- cut(ppc_obs$pred_median,
      breaks = c(-Inf,.10,.25,.50,.75,.90,Inf), right = FALSE,
      labels = c("<0.10","0.10--0.25","0.25--0.50","0.50--0.75","0.75--0.90",">=0.90"))
    ppc_obs$expert_count_bin <- cut(ppc_obs$n_experts,
      breaks = c(-Inf,1,3,5,10,Inf), right = TRUE,
      labels = c("1","2--3","4--5","6--10",">10"))
    width_breaks <- unique(quantile(ppc_obs$interval_width, seq(0, 1, .25), na.rm = TRUE))
    ppc_obs$interval_width_quartile <- cut(ppc_obs$interval_width, breaks = width_breaks,
      include.lowest = TRUE, labels = paste0("Q", seq_len(length(width_breaks)-1)))

    summarize_ppc_bins <- function(variable) {
      pieces <- split(ppc_obs, interaction(ppc_obs$dim_idx, ppc_obs[[variable]], drop = TRUE))
      do.call(rbind, lapply(pieces, function(d) {
        data.frame(
          dim_idx = d$dim_idx[1], group_type = variable,
          group = as.character(d[[variable]][1]), n = nrow(d),
          covered = sum(d$covered_95), observed_coverage = mean(d$covered_95),
          mean_interval_width = mean(d$interval_width),
          mean_fitted_value = mean(d$pred_median), mean_expert_count = mean(d$n_experts, na.rm = TRUE)
        )
      }))
    }
    ppc_residual <- do.call(rbind, lapply(
      c("fitted_bin", "expert_count_bin", "interval_width_quartile"), summarize_ppc_bins))
    write_clean_csv(ppc_residual, file.path(ppc_dir, "posterior_predictive_residual_patterns.csv"))
  }
}

message("Post-processing complete.")