party2d/diagnostics/generate_diagnostics.R

#!/usr/bin/env Rscript

suppressPackageStartupMessages(library(tidyverse))

find_repo_root <- function() {
  args <- commandArgs(trailingOnly = FALSE)
  file_arg <- "--file="
  script_arg <- args[startsWith(args, file_arg)][1]
  if (!is.na(script_arg)) {
    return(normalizePath(file.path(dirname(sub(file_arg, "", script_arg)), "..")))
  }
  if (file.exists("data/text_data.csv")) return(normalizePath(getwd()))
  stop("Cannot find repository root. Run from the party2d repository root.")
}

repo_root <- find_repo_root()
setwd(repo_root)

release_version <- Sys.getenv("PARTY2D_RELEASE_VERSION", "v0")
outputs_dir <- Sys.getenv("PARTY2D_OUTPUTS_DIR", "outputs")
if (!grepl("^/", outputs_dir)) outputs_dir <- file.path(repo_root, outputs_dir)
if (!dir.exists(outputs_dir)) {
  stop("Model output directory not found: ", outputs_dir, ". Run estimation/validation first or set PARTY2D_OUTPUTS_DIR.")
}
review_package_dir <- Sys.getenv("PARTY2D_REVIEW_PACKAGE_DIR", file.path(dirname(repo_root), "archive", "review_package"))
if (!grepl("^/", review_package_dir)) review_package_dir <- file.path(repo_root, review_package_dir)

generated_dir <- file.path(repo_root, "diagnostics", "generated")
release_dir <- file.path(repo_root, "data", "releases")
dir.create(generated_dir, recursive = TRUE, showWarnings = FALSE)
dir.create(release_dir, recursive = TRUE, showWarnings = FALSE)

required_inputs <- c(
  "data/text_data.csv",
  "data/expert.csv",
  "data/lr_data.csv",
  "data/union_mapping.csv",
  "data/party_families.csv"
)
missing_inputs <- required_inputs[!file.exists(required_inputs)]
if (length(missing_inputs) > 0) {
  stop("Missing required model input files: ", paste(missing_inputs, collapse = ", "))
}

latest_file <- function(path, pattern) {
  if (!dir.exists(path)) return(NA_character_)
  files <- list.files(path, pattern = pattern, full.names = TRUE)
  if (length(files) == 0) return(NA_character_)
  sort(files)[length(files)]
}

read_if_exists <- function(path) {
  if (is.na(path) || !file.exists(path)) return(tibble())
  readr::read_csv(path, show_col_types = FALSE)
}

review_file <- function(...) {
  path <- file.path(review_package_dir, ...)
  if (file.exists(path)) path else NA_character_
}

public_dimension <- function(x) {
  dplyr::recode(as.character(x),
    economic_lr = "economic left-right",
    galtan = "cultural cosmopolitan--traditionalist",
    Economic = "economic left-right",
    Cultural = "cultural cosmopolitan--traditionalist",
    `Economic Left-Right` = "economic left-right",
    .default = as.character(x)
  )
}

infer_dimension <- function(type_low, type_high) {
  dplyr::case_when(
    type_low %in% c("pro_market", "pro_welfare", "left", "right") |
      type_high %in% c("pro_market", "pro_welfare", "left", "right") ~ "economic left-right",
    type_low %in% c("cosmopolitan", "traditional") |
      type_high %in% c("cosmopolitan", "traditional") ~ "cultural cosmopolitan--traditionalist",
    TRUE ~ "general left-right"
  )
}

is_reversed_for_reporting <- function(type_high) {
  type_high %in% c("pro_welfare", "left", "cosmopolitan")
}

fmt_num <- function(x, digits = 3) {
  ifelse(is.na(x), "NA", formatC(x, digits = digits, format = "f"))
}

display_path <- function(path) {
  if (is.na(path) || !nzchar(path)) return("not available")
  normalized <- normalizePath(path, mustWork = FALSE)
  root_prefix <- paste0(normalizePath(repo_root, mustWork = FALSE), .Platform$file.sep)
  if (startsWith(normalized, root_prefix)) return(sub(root_prefix, "", normalized, fixed = TRUE))
  basename(path)
}

md_table <- function(df, n = Inf) {
  if (nrow(df) == 0) return("_No rows available._\n")
  df <- head(df, n)
  df <- mutate(df, across(everything(), as.character))
  header <- paste0("| ", paste(names(df), collapse = " | "), " |")
  sep <- paste0("| ", paste(rep("---", ncol(df)), collapse = " | "), " |")
  rows <- apply(df, 1, function(x) paste0("| ", paste(x, collapse = " | "), " |"))
  paste(c(header, sep, rows), collapse = "\n")
}

text_data <- read_csv("data/text_data.csv", show_col_types = FALSE)
expert <- read_csv("data/expert.csv", show_col_types = FALSE)
lr_data <- read_csv("data/lr_data.csv", show_col_types = FALSE)
union_mapping <- read_csv("data/union_mapping.csv", show_col_types = FALSE)
party_families <- read_csv("data/party_families.csv", show_col_types = FALSE)

excluded_poldem <- text_data %>%
  filter(project == "PolDem", str_detect(str_to_lower(var), "reform"))
if (nrow(excluded_poldem) > 0) {
  stop("Excluded PolDem reform item is present in data/text_data.csv. Final-model diagnostics must be regenerated after removing it.")
}

annual_release <- file.path(release_dir, paste0("party_2d_annual_model_output_", release_version, ".csv.gz"))
panel_release <- file.path(release_dir, paste0("party_2d_election_year_panel_", release_version, ".csv.gz"))
model_positions_file <- latest_file(file.path(outputs_dir, "estimations", "latest"), "^party_positions_.*\\.csv$")
if (is.na(model_positions_file)) {
  stop("No post-estimation party-position output found under ", outputs_dir, ". Run model estimation/post-estimation first, or set PARTY2D_OUTPUTS_DIR to an outputs directory.")
}
positions <- read_csv(model_positions_file, show_col_types = FALSE)

item_rows <- bind_rows(
  text_data %>%
    mutate(source_file = "text_data.csv") %>%
    group_by(source_file, item = var, source = project, type_low, type_high) %>%
    summarise(observations = n(), party_years = n_distinct(party, year), parties = n_distinct(party), countries = n_distinct(country), min_year = min(year), max_year = max(year), .groups = "drop"),
  expert %>%
    mutate(source_file = "expert.csv") %>%
    group_by(source_file, item = var, source = project, type_low, type_high) %>%
    summarise(observations = n(), party_years = n_distinct(party, year), parties = n_distinct(party), countries = n_distinct(country), min_year = min(year), max_year = max(year), .groups = "drop"),
  lr_data %>%
    mutate(source_file = "lr_data.csv", type_low = NA_character_, type_high = NA_character_) %>%
    group_by(source_file, item = var, source = project, type_low, type_high) %>%
    summarise(observations = n(), party_years = n_distinct(party, year), parties = n_distinct(party), countries = n_distinct(country), min_year = min(year), max_year = max(year), .groups = "drop")
) %>%
  mutate(
    dimension = infer_dimension(type_low, type_high),
    higher_values_indicate = if_else(is.na(type_high), "source-coded left-right", type_high),
    reversed_for_reporting = if_else(is_reversed_for_reporting(type_high), "yes", "no")
  ) %>%
  select(source_file, item, source, dimension, type_low, type_high, higher_values_indicate, reversed_for_reporting, observations, party_years, parties, countries, min_year, max_year) %>%
  arrange(source_file, source, dimension, item)

source_coverage <- bind_rows(
  text_data %>% transmute(source_file = "text_data.csv", source = project, item = var, party, country, year),
  expert %>% transmute(source_file = "expert.csv", source = project, item = var, party, country, year),
  lr_data %>% transmute(source_file = "lr_data.csv", source = project, item = var, party, country, year)
) %>%
  group_by(source_file, source) %>%
  summarise(items = n_distinct(item), observations = n(), party_years = n_distinct(party, year), parties = n_distinct(party), countries = n_distinct(country), min_year = min(year), max_year = max(year), .groups = "drop") %>%
  arrange(source_file, source)

party_year_source_coverage <- bind_rows(
  text_data %>% distinct(party, country, year) %>% mutate(has_text = TRUE, has_expert = FALSE, has_general_lr = FALSE),
  expert %>% distinct(party, country, year) %>% mutate(has_text = FALSE, has_expert = TRUE, has_general_lr = FALSE),
  lr_data %>% distinct(party, country, year) %>% mutate(has_text = FALSE, has_expert = FALSE, has_general_lr = TRUE)
) %>%
  group_by(party, country, year) %>%
  summarise(has_text = any(has_text), has_expert = any(has_expert), has_general_lr = any(has_general_lr), n_source_types = has_text + has_expert + has_general_lr, .groups = "drop") %>%
  arrange(year, country, party)

alliance_union_harmonization <- bind_rows(
  tibble(metric = "constituent_mappings", category = "all", value = nrow(union_mapping)),
  tibble(metric = "unique_union_or_alliance_ids", category = "all", value = n_distinct(union_mapping$manifesto_pf_id)),
  tibble(metric = "unique_constituent_party_ids", category = "all", value = n_distinct(union_mapping$expert_pf_id)),
  union_mapping %>% count(country, name = "value") %>% transmute(metric = "mappings_by_country", category = country, value),
  union_mapping %>% count(status, name = "value") %>% transmute(metric = "mappings_by_status", category = status, value)
)

party_col <- if ("party_id" %in% names(positions)) "party_id" else "party"
party_family_coverage <- positions %>%
  transmute(partyfacts_id = .data[[party_col]], country, year) %>%
  inner_join(party_families, by = "partyfacts_id") %>%
  group_by(family) %>%
  summarise(parties = n_distinct(partyfacts_id), party_years = n(), countries = n_distinct(country), min_year = min(year), max_year = max(year), .groups = "drop") %>%
  arrange(desc(party_years))

convergence_summary_file <- latest_file(file.path(outputs_dir, "diagnostics"), "^convergence_summary_.*\\.csv$")
convergence_detail_file <- latest_file(file.path(outputs_dir, "diagnostics"), "^convergence_diagnostics_.*\\.csv$")
if (is.na(convergence_summary_file) || is.na(convergence_detail_file)) {
  stop("Convergence diagnostics not found under ", outputs_dir, ". Run the model diagnostics before generating the report.")
}
model_convergence_summary <- read_if_exists(convergence_summary_file) %>%
  identity()
model_convergence_by_dimension <- read_if_exists(convergence_detail_file) %>%
  group_by(dimension) %>%
  summarise(parameters = n(), mean_rhat = mean(rhat, na.rm = TRUE), max_rhat = max(rhat, na.rm = TRUE), min_ess_bulk = min(ess_bulk, na.rm = TRUE), mean_ess_bulk = mean(ess_bulk, na.rm = TRUE), .groups = "drop") %>%
  mutate(dimension = public_dimension(dimension)) %>%
  arrange(dimension)

convergent_summary_file <- latest_file(file.path(outputs_dir, "validation", "latest"), "^convergent_summary_.*\\.csv$")
discriminant_summary_file <- latest_file(file.path(outputs_dir, "validation", "latest"), "^discriminant_summary_.*\\.csv$")
uncertainty_summary_file <- latest_file(file.path(outputs_dir, "validation", "latest"), "^uncertainty_cic_summary_.*\\.csv$")
external_validation_file <- latest_file(file.path(outputs_dir, "validation", "latest"), "^external_validation_.*\\.csv$")
construct_families_file <- latest_file(file.path(outputs_dir, "validation", "latest"), "^construct_families_.*\\.csv$")
construct_unstable_file <- latest_file(file.path(outputs_dir, "validation", "latest"), "^construct_unstable_.*\\.csv$")
if (any(is.na(c(convergent_summary_file, discriminant_summary_file, uncertainty_summary_file, external_validation_file, construct_families_file, construct_unstable_file)))) {
  stop("Validation diagnostics not found under ", outputs_dir, ". Run validation before generating the report.")
}

convergent_summary <- read_if_exists(convergent_summary_file) %>%
  mutate(diagnostic = "convergent validity", dimension = public_dimension(dimension))
discriminant_summary <- read_if_exists(discriminant_summary_file) %>%
  mutate(diagnostic = "discriminant validity", model_dim = public_dimension(model_dim), expert_dim = public_dimension(expert_dim))
uncertainty_summary <- read_if_exists(uncertainty_summary_file) %>%
  mutate(diagnostic = "posterior predictive coverage", dimension = public_dimension(dimension))

external_validation_correlations <- read_if_exists(external_validation_file) %>%
  group_by(var, dimension) %>%
  summarise(n = n(), pearson_r = cor(expert_val, model_val, use = "complete.obs"), mean_absolute_error = mean(abs_error, na.rm = TRUE), coverage_95 = mean(covered_95, na.rm = TRUE), .groups = "drop") %>%
  mutate(dimension = public_dimension(dimension)) %>%
  arrange(dimension, var)
construct_family_positions <- read_if_exists(construct_families_file) %>%
  rename(mean_cultural = mean_galtan, sd_cultural = sd_galtan) %>%
  arrange(mean_economic)
construct_temporal_stability <- read_if_exists(construct_unstable_file) %>%
  mutate(dimension = public_dimension(dimension)) %>%
  arrange(desc(annual_change))
source_composition_balance <- read_if_exists(review_file("validation", "source_composition_balance.csv")) %>%
  mutate(dimension = public_dimension(dimension))
robustness_sensitivity <- read_if_exists(review_file("validation", "table10_sensitivity.csv")) %>%
  mutate(
    dimension = public_dimension(dimension),
    across(everything(), ~ na_if(as.character(.x), "[INSERT VALUE]"))
  ) %>%
  select(specification, ablated_source, dimension, matched_n, correlation_with_production,
         mean_abs_difference, median_abs_difference, p95_abs_difference,
         mean_interval_width_production, mean_interval_width_ablation)

posterior_uncertainty <- positions %>%
  summarise(
    rows = n(),
    parties = n_distinct(.data[[party_col]]),
    countries = n_distinct(country),
    min_year = min(year),
    max_year = max(year),
    mean_economic_se = mean(economic_lr_se, na.rm = TRUE),
    median_economic_se = median(economic_lr_se, na.rm = TRUE),
    mean_cultural_se = mean(galtan_se, na.rm = TRUE),
    median_cultural_se = median(galtan_se, na.rm = TRUE)
  )

write_csv(item_rows, file.path(generated_dir, "item_coverage.csv"))
write_csv(source_coverage, file.path(generated_dir, "source_coverage.csv"))
write_csv(party_year_source_coverage, file.path(generated_dir, "party_year_source_coverage.csv"))
write_csv(item_rows, file.path(generated_dir, "item_coding_orientation.csv"))
write_csv(filter(item_rows, reversed_for_reporting == "yes"), file.path(generated_dir, "reversed_items.csv"))
write_csv(alliance_union_harmonization, file.path(generated_dir, "alliance_union_harmonization.csv"))
write_csv(party_family_coverage, file.path(generated_dir, "party_family_coverage.csv"))
write_csv(model_convergence_summary, file.path(generated_dir, "model_convergence_summary.csv"))
write_csv(model_convergence_by_dimension, file.path(generated_dir, "model_convergence_by_dimension.csv"))
write_csv(convergent_summary, file.path(generated_dir, "posterior_validation_convergent_summary.csv"))
write_csv(discriminant_summary, file.path(generated_dir, "posterior_validation_discriminant_summary.csv"))
write_csv(uncertainty_summary, file.path(generated_dir, "posterior_validation_uncertainty_summary.csv"))
write_csv(external_validation_correlations, file.path(generated_dir, "external_validation_correlations.csv"))
write_csv(construct_family_positions, file.path(generated_dir, "construct_family_positions.csv"))
write_csv(construct_temporal_stability, file.path(generated_dir, "construct_temporal_stability_flags.csv"))
write_csv(source_composition_balance, file.path(generated_dir, "source_composition_balance.csv"))
write_csv(robustness_sensitivity, file.path(generated_dir, "robustness_sensitivity.csv"))
write_csv(posterior_uncertainty, file.path(generated_dir, "posterior_uncertainty_summary.csv"))

item_counts <- item_rows %>% count(source_file, dimension, name = "items")
source_counts <- source_coverage %>% select(source_file, source, items, observations, party_years, parties, countries, min_year, max_year)
reversed_items <- filter(item_rows, reversed_for_reporting == "yes") %>% select(item, source, dimension, higher_values_indicate, observations, min_year, max_year)
conv_display <- model_convergence_summary %>% select(-any_of("source_file"))
conv_dim_display <- model_convergence_by_dimension %>% select(-any_of("source_file")) %>% mutate(across(where(is.numeric), ~ round(.x, 3)))
val_display <- bind_rows(
  convergent_summary %>% transmute(diagnostic, dimension, n, pearson_r = round(r_pearson, 3), spearman_r = round(r_spearman, 3), mae = round(mae, 3), coverage = NA_real_),
  uncertainty_summary %>% transmute(diagnostic, dimension, n, pearson_r = NA_real_, spearman_r = NA_real_, mae = NA_real_, coverage = round(cic, 3))
)

report_lines <- c(
  "# Diagnostics report",
  "",
  paste0("Generated: ", format(Sys.time(), "%Y-%m-%d %H:%M:%S %Z")),
  paste0("Release: ", release_version),
  paste0("Model positions source: `", display_path(model_positions_file), "`"),
  "",
  "## Purpose",
  "",
  "The purpose of this report is to provide the appendix-style diagnostics that document how the released party-position estimates are constructed, checked, and validated. The main article reports the central validation evidence; this report keeps the larger technical tables with the release so readers can inspect item coverage, source coverage, coding orientation, harmonization, convergence, posterior uncertainty, and validation details in one reproducible place.",
  "",
  "## Overview",
  "",
  "This report follows the structure of the technical appendix material: data and item coverage, coding and scale orientation, party-union harmonization, construct checks, model convergence, and validation. It is generated from the model-ready inputs and completed model outputs; it is not part of the raw-data setup workflow.",
  "",
  "## Data and item coverage",
  "",
  "The model combines text-coded item counts, dimension-specific expert placements, and general left-right expert placements. Text items enter as positive/sample counts, expert items enter as aggregated ratings with scale and expert-count information, and general left-right ratings anchor the relationship between the two dimensions.",
  "",
  md_table(item_counts),
  "",
  "### Source coverage",
  "",
  md_table(source_counts),
  "",
  "## Data coding and item orientation",
  "",
  "All indicators are oriented toward the two reported dimensions: economic left-right and cultural cosmopolitan--traditionalist. For interpretability, generated diagnostics report whether higher observed values point toward the public high pole or are reversed for reporting. Original source variable names are preserved in the tables.",
  "",
  "### Reversed items",
  "",
  md_table(reversed_items),
  "",
  "## Party unions and electoral coalitions",
  "",
  "Alliance and union labels are handled through constituent mappings so the released party identifiers represent individual parties. Shared text evidence can inform constituent parties through the union mapping while expert data continue to constrain individual parties directly.",
  "",
  md_table(alliance_union_harmonization %>% head(30)),
  "",
  "## Party-family coverage",
  "",
  "Party-family classifications are used for construct-validity diagnostics and coverage summaries. The table below reports coverage in the completed model output by family code.",
  "",
  md_table(party_family_coverage),
  "",
  "### Construct-validity family means",
  "",
  "Substantive party-family means provide a construct-validity check: families should follow the expected ordering on the economic left-right and cultural cosmopolitan--traditionalist dimensions.",
  "",
  md_table(construct_family_positions %>% select(family_name, n_parties, n_obs, mean_economic, sd_economic, mean_cultural, sd_cultural) %>% mutate(across(where(is.numeric), ~ round(.x, 3)))),
  "",
  "### Temporal-stability flags",
  "",
  "The model permits movement through random walks, but unusually large one-year changes are flagged for inspection rather than treated as automatic errors.",
  "",
  md_table(construct_temporal_stability %>% select(party_id, country, dimension, year_from, year_to, val_from, val_to, annual_change) %>% mutate(across(where(is.numeric), ~ round(.x, 3))), n = 20),
  "",
  "## Model convergence diagnostics",
  "",
  if (nrow(model_convergence_summary) > 0) "Convergence is assessed using split R-hat and effective sample size over monitored parameters." else "Convergence summary files were not found in the configured outputs directory.",
  "",
  md_table(conv_display),
  "",
  "### Convergence by parameter group",
  "",
  md_table(conv_dim_display),
  "",
  "## Posterior uncertainty",
  "",
  "The completed party-position output reports posterior standard errors and interval endpoints for both dimensions. These summaries describe the typical uncertainty in the release file used by the report.",
  "",
  md_table(posterior_uncertainty %>% mutate(across(where(is.numeric), ~ round(.x, 3)))),
  "",
  "## Validation diagnostics",
  "",
  "The validation diagnostics combine convergent and discriminant comparisons with expert surveys, posterior predictive coverage, construct checks, and out-of-sample validation when the corresponding outputs are available.",
  "",
  md_table(val_display),
  "",
  "### Discriminant validity",
  "",
  md_table(discriminant_summary %>% select(-any_of("source_file")) %>% mutate(across(where(is.numeric), ~ round(.x, 3)))),
  "",
  "### External validation correlations",
  "",
  md_table(external_validation_correlations %>% select(-any_of("source_file")) %>% mutate(across(where(is.numeric), ~ round(.x, 3)))),
  "",
  "## Evidence-composition balance",
  "",
  "Evidence-composition balance checks whether estimates informed by different nearby source combinations are systematically shifted relative to rows with both text and expert evidence. The reported differences are adjusted differences on the unit scale relative to the overlapping text-and-expert reference category.",
  "",
  md_table(source_composition_balance),
  "",
  "## Robustness and sensitivity checks",
  "",
  "Sensitivity checks compare the released election-year estimates with source-ablation, segmentation-threshold, and item-screening variants where available. Correlations near one and small absolute differences indicate that the released estimates are stable to the corresponding design choice.",
  "",
  md_table(robustness_sensitivity),
  "",
  "## Generated tables",
  "",
  paste0("- `", list.files(generated_dir, pattern = "\\.csv$"), "`"),
  ""
)

pdf_source <- file.path(generated_dir, "diagnostics_report.Rmd")
pdf_file <- file.path(generated_dir, "diagnostics_report.pdf")
release_pdf <- file.path(release_dir, paste0("party_2d_diagnostics_report_", release_version, ".pdf"))
pdf_lines <- c(
  "---",
  "title: \"Diagnostics report\"",
  paste0("date: \"", format(Sys.time(), "%Y-%m-%d"), "\""),
  "output:",
  "  pdf_document:",
  "    toc: true",
  "    number_sections: true",
  "    latex_engine: pdflatex",
  "geometry: margin=0.75in",
  "fontsize: 10pt",
  "header-includes:",
  "  - \\usepackage{booktabs}",
  "  - \\usepackage{longtable}",
  "  - \\usepackage{array}",
  "  - \\usepackage{pdflscape}",
  "  - \\setlength{\\tabcolsep}{4pt}",
  "  - \\renewcommand{\\arraystretch}{1.12}",
  "---",
  "",
  "```{r setup, include=FALSE}",
  "knitr::opts_chunk$set(echo = FALSE, message = FALSE, warning = FALSE)",
  "print_table <- function(x, n = Inf, size = 'footnotesize') {",
  "  if (nrow(x) == 0) { cat('No rows available.\\n'); return(invisible(NULL)) }",
  "  x <- head(x, n)",
  "  x <- mutate(x, across(everything(), as.character))",
  "  x[is.na(x)] <- ''",
  "  names(x) <- gsub('_', ' ', names(x), fixed = TRUE)",
  "  cat(paste0(\"\\n\\\\begingroup\\\\\", size, \"\\n\"))",
  "  print(knitr::kable(x, format = 'latex', booktabs = TRUE, longtable = FALSE, digits = 3))",
  "  cat(\"\\n\\\\endgroup\\n\")",
  "}",
  "short_dim <- function(x) dplyr::recode(as.character(x), 'cultural cosmopolitan--traditionalist' = 'cultural', 'economic left-right' = 'economic', .default = as.character(x))",
  "```",
  "",
  paste0("Generated: ", format(Sys.time(), "%Y-%m-%d %H:%M:%S %Z")),
  "",
  paste0("Release: ", release_version),
  "",
  paste0("Model positions source: `", display_path(model_positions_file), "`"),
  "",
  "# Purpose",
  "",
  "The purpose of this report is to provide the appendix-style diagnostics that document how the released party-position estimates are constructed, checked, and validated. The main article reports the central validation evidence; this report keeps the larger technical tables with the release so readers can inspect item coverage, source coverage, coding orientation, harmonization, convergence, posterior uncertainty, and validation details in one reproducible place.",
  "",
  "# Overview",
  "",
  "This report follows the structure of the technical appendix material: data and item coverage, coding and scale orientation, party-union harmonization, construct checks, model convergence, and validation. It is generated from the model-ready inputs and completed model outputs; it is not part of the raw-data setup workflow.",
  "",
  "# Data and item coverage",
  "",
  "The model combines text-coded item counts, dimension-specific expert placements, and general left-right expert placements. Text items enter as positive/sample counts, expert items enter as aggregated ratings with scale and expert-count information, and general left-right ratings anchor the relationship between the two dimensions.",
  "",
  "```{r item-counts, results='asis'}",
  "print_table(item_counts %>% mutate(dimension = short_dim(dimension)))",
  "```",
  "",
  "## Source coverage",
  "",
  "```{r source-coverage, results='asis'}",
  "print_table(source_counts %>% transmute(file = recode(source_file, text_data.csv = 'text', expert.csv = 'expert', lr_data.csv = 'general LR'), source, items, obs = observations, party_years, parties, countries, years = paste0(min_year, '--', max_year)), size = 'scriptsize')",
  "```",
  "",
  "Full item-level coverage and coding-orientation details are provided as generated CSV files listed at the end of this report.",
  "",
  "# Data coding and item orientation",
  "",
  "All indicators are oriented toward the two reported dimensions: economic left-right and cultural cosmopolitan--traditionalist. For interpretability, generated diagnostics report whether higher observed values point toward the public high pole or are reversed for reporting. Original source variable names are preserved in the tables.",
  "",
  "## Reversed items",
  "",
  "```{r reversed-items, results='asis'}",
  "print_table(reversed_items %>% count(source, dimension = short_dim(dimension), higher_values_indicate, name = 'items'))",
  "```",
  "",
  "# Party unions and electoral coalitions",
  "",
  "Alliance and union labels are handled through constituent mappings so the released party identifiers represent individual parties. Shared text evidence can inform constituent parties through the union mapping while expert data continue to constrain individual parties directly.",
  "",
  "```{r union-summary, results='asis'}",
  "print_table(alliance_union_harmonization %>% transmute(metric, category, value), n = 40)",
  "```",
  "",
  "# Party-family and construct coverage",
  "",
  "Party-family classifications are used for construct-validity diagnostics and coverage summaries. The table below reports coverage in the completed model output by family code.",
  "",
  "```{r family-coverage, results='asis'}",
  "print_table(party_family_coverage %>% transmute(family, parties, party_years, countries, years = paste0(min_year, '--', max_year)))",
  "```",
  "",
  "## Construct-validity family means",
  "",
  "Substantive party-family means provide a construct-validity check: families should follow the expected ordering on the economic left-right and cultural cosmopolitan--traditionalist dimensions.",
  "",
  "```{r construct-family, results='asis'}",
  "print_table(construct_family_positions %>% transmute(family = family_name, parties = n_parties, obs = n_obs, econ_mean = round(mean_economic, 3), econ_sd = round(sd_economic, 3), cult_mean = round(mean_cultural, 3), cult_sd = round(sd_cultural, 3)), size = 'scriptsize')",
  "```",
  "",
  "## Temporal-stability flags",
  "",
  "The model permits movement through random walks, but unusually large one-year changes are flagged for inspection rather than treated as automatic errors.",
  "",
  "```{r temporal-stability, results='asis'}",
  "print_table(construct_temporal_stability %>% transmute(party = party_id, country, dim = short_dim(dimension), from = year_from, to = year_to, start = round(val_from, 3), end = round(val_to, 3), annual_change = round(annual_change, 3)), n = 12, size = 'scriptsize')",
  "```",
  "",
  "# Model convergence diagnostics",
  "",
  "Convergence is assessed using split R-hat and effective sample size over monitored parameters.",
  "",
  "```{r convergence-summary, results='asis'}",
  "print_table(conv_display)",
  "```",
  "",
  "## Convergence by parameter group",
  "",
  "```{r convergence-dim, results='asis'}",
  "print_table(conv_dim_display %>% mutate(dimension = short_dim(dimension)))",
  "```",
  "",
  "# Posterior uncertainty",
  "",
  "The completed party-position output reports posterior standard errors and interval endpoints for both dimensions. These summaries describe the typical uncertainty in the release file used by the report.",
  "",
  "```{r posterior-uncertainty, results='asis'}",
  "print_table(posterior_uncertainty %>% transmute(rows, parties, countries, years = paste0(min_year, '--', max_year), mean_econ_se = round(mean_economic_se, 3), median_econ_se = round(median_economic_se, 3), mean_cult_se = round(mean_cultural_se, 3), median_cult_se = round(median_cultural_se, 3)), size = 'scriptsize')",
  "```",
  "",
  "# Validation diagnostics",
  "",
  "The validation diagnostics combine convergent and discriminant comparisons with expert surveys, posterior predictive coverage, construct checks, and out-of-sample validation when the corresponding outputs are available.",
  "",
  "```{r validation-summary, results='asis'}",
  "print_table(val_display %>% mutate(dimension = short_dim(dimension)), size = 'scriptsize')",
  "```",
  "",
  "## Discriminant validity",
  "",
  "```{r discriminant, results='asis'}",
  "print_table(discriminant_summary %>% transmute(type, model = short_dim(model_dim), expert = expert_dim, n, pearson = round(r_pearson, 3), spearman = round(r_spearman, 3)))",
  "```",
  "",
  "## External validation correlations",
  "",
  "```{r external-validation, results='asis'}",
  "print_table(external_validation_correlations %>% transmute(item = var, dim = short_dim(dimension), n, r = round(pearson_r, 3), mae = round(mean_absolute_error, 3), coverage = round(coverage_95, 3)))",
  "```",
  "",
  "# Evidence-composition balance",
  "",
  "Evidence-composition balance checks whether estimates informed by different nearby source combinations are systematically shifted relative to rows with both text and expert evidence. The reported differences are adjusted differences on the unit scale relative to the overlapping text-and-expert reference category.",
  "",
  "```{r source-balance, results='asis'}",
  "print_table(source_composition_balance %>% transmute(dim = short_dim(dimension), evidence = recode(source_composition_class, both_direct_or_nearby = 'both', text_only_direct_or_nearby = 'text only', expert_only_direct_or_nearby = 'expert only', temporal_propagation = 'temporal'), ref = recode(reference_class, both_direct_or_nearby = 'both'), n, adj_diff = round(adjusted_difference, 3)))",
  "```",
  "",
  "# Robustness and sensitivity checks",
  "",
  "Sensitivity checks compare the released election-year estimates with source-ablation, segmentation-threshold, and item-screening variants where available. Correlations near one and small absolute differences indicate that the released estimates are stable to the corresponding design choice.",
  "",
  "```{r robustness-sensitivity, results='asis'}",
  "print_table(robustness_sensitivity %>% transmute(spec = specification, source = ablated_source, dim = short_dim(dimension), n = matched_n, r = round(as.numeric(correlation_with_production), 3), mean_abs = round(as.numeric(mean_abs_difference), 3), median_abs = round(as.numeric(median_abs_difference), 3), p95_abs = round(as.numeric(p95_abs_difference), 3)), size = 'scriptsize')",
  "```",
  "",
  "# Generated tables",
  "",
  paste0("- `", list.files(generated_dir, pattern = "\\.csv$"), "`")
)
writeLines(pdf_lines, pdf_source)
if (!requireNamespace("rmarkdown", quietly = TRUE)) {
  stop("The rmarkdown package is required to render the diagnostics PDF.")
}
rmarkdown::render(
  input = pdf_source,
  output_format = rmarkdown::pdf_document(toc = TRUE, number_sections = TRUE),
  output_file = basename(pdf_file),
  output_dir = generated_dir,
  quiet = TRUE,
  envir = environment()
)
invisible(file.copy(pdf_file, release_pdf, overwrite = TRUE))
unlink(c(
  pdf_source,
  file.path(generated_dir, "diagnostics_report.log"),
  file.path(generated_dir, "diagnostics_report.aux"),
  file.path(generated_dir, "diagnostics_report.out"),
  file.path(generated_dir, "diagnostics_report.toc"),
  file.path(generated_dir, "diagnostics_report.tex")
), force = TRUE)

generated_readme <- c(
  "# Generated diagnostics",
  "",
  "These files are generated by:",
  "",
  "```bash",
  "Rscript diagnostics/generate_diagnostics.R",
  "```",
  "",
  "The command requires completed model/post-estimation outputs. If those outputs are outside the repository root, set `PARTY2D_OUTPUTS_DIR` before running the script.",
  "",
  "The report file is `diagnostics_report.pdf`; the same PDF is copied into `data/releases/` for the release files."
)
writeLines(generated_readme, file.path(generated_dir, "README.md"))

sha_file <- file.path(release_dir, "SHA256SUMS")
release_files_for_sha <- c(
  paste0("party_2d_election_year_panel_", release_version, ".csv.gz"),
  paste0("party_2d_annual_model_output_", release_version, ".csv.gz"),
  basename(release_pdf)
)
existing_release_files <- release_files_for_sha[file.exists(file.path(release_dir, release_files_for_sha))]
sha_lines <- vapply(existing_release_files, function(f) {
  old <- getwd()
  on.exit(setwd(old), add = TRUE)
  setwd(release_dir)
  system2("sha256sum", f, stdout = TRUE)
}, character(1))
writeLines(sha_lines, sha_file)

message("Diagnostics written to diagnostics/generated")
message("Release diagnostics PDF written to ", release_pdf)