Armin
385 lines
20 KiB
R
Executable File
385 lines
20 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)
|
|
|
|
# ---------------------------------------------------------------------------
|
|
# Historically anchored scale examples
|
|
# ---------------------------------------------------------------------------
|
|
|
|
source(file.path("validation", "plot_scale_examples.R"), local = TRUE)
|
|
|
|
# ---------------------------------------------------------------------------
|
|
# 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.")
|