Armin Seimel
465 lines
24 KiB
R
Executable File
465 lines
24 KiB
R
Executable File
#!/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: SPD", "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("SPD 1972","SPD 2021","CDU 1983","CDU 2021","Labour 1983","Labour 1997",
|
||
"Conservatives 1979","Conservatives 2019","Swedish SAP 1994","Swedish SAP 2022",
|
||
"Sweden Democrats 2010","Sweden Democrats 2022","French FN 1988","French FN 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", relig_vparty="V-Party religion",
|
||
welf_vparty="V-Party welfare", culsup_vparty="V-Party cultural superiority",
|
||
immig_vparty="V-Party immigration", lgbt_vparty="V-Party LGBT",
|
||
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.")
|