Add diagnostics report workflow

This commit is contained in:
aseimel
2026-06-15 16:33:19 +02:00
parent 91403e2bdf
commit bc2d616605
24 changed files with 7896 additions and 0 deletions
+387
View File
@@ -0,0 +1,387 @@
#!/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.")
}
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)
}
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(relationship, name = "value") %>% transmute(metric = "mappings_by_relationship", category = relationship, 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) %>%
mutate(source_file = convergence_summary_file)
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), source_file = convergence_detail_file) %>%
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), source_file = external_validation_file) %>%
arrange(dimension, var)
construct_family_positions <- read_if_exists(construct_families_file) %>%
rename(mean_cultural = mean_galtan, sd_cultural = sd_galtan) %>%
mutate(source_file = construct_families_file) %>%
arrange(mean_economic)
construct_temporal_stability <- read_if_exists(construct_unstable_file) %>%
mutate(dimension = public_dimension(dimension), source_file = construct_unstable_file) %>%
arrange(desc(annual_change))
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(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), "`"),
"",
"## Overview",
"",
"This report collects the large technical diagnostics that support the Scientific Data Data Descriptor. It follows the structure of the earlier technical appendix: data and item coverage, coding and scale orientation, party-union harmonization, model convergence, and validation. The report 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)))),
"",
"## Generated tables",
"",
paste0("- `", list.files(generated_dir, pattern = "\\.csv$"), "`"),
""
)
report_path <- file.path(generated_dir, "diagnostics_report.md")
writeLines(report_lines, report_path)
release_report <- file.path(release_dir, paste0("party_2d_diagnostics_report_", release_version, ".md"))
invisible(file.copy(report_path, release_report, overwrite = 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 main Markdown report is `diagnostics_report.md`; the same report 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_report)
)
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 report written to ", release_report)