From c811159c29f3e4a99add2f9af5155d10110df2a7 Mon Sep 17 00:00:00 2001 From: cm401 Date: Thu, 7 May 2026 11:04:04 +0100 Subject: [PATCH 01/18] Add data-format ablation analysis script Fits three model variants per eligible pathogen (individual-level only, summary-statistics only, federated) and compares posterior predictive summaries, interval ratios, W2/JS/OVL distributional metrics, tail probabilities, and study-level PSIS-LOO ELPD. Co-Authored-By: Claude Sonnet 4.6 --- analysis/data_format_ablation.R | 784 ++++++++++++++++++++++++++++++++ 1 file changed, 784 insertions(+) create mode 100644 analysis/data_format_ablation.R diff --git a/analysis/data_format_ablation.R b/analysis/data_format_ablation.R new file mode 100644 index 0000000..44a3275 --- /dev/null +++ b/analysis/data_format_ablation.R @@ -0,0 +1,784 @@ +# ============================================================================= +# data_format_ablation.R +# ----------------------------------------------------------------------------- +# "Data-format ablation" experiment on the real dataset. +# +# For every pathogen that has at least MIN_DATASETS_PER_ARM datasets of +# individual-level format (summary types 4/5: exact or interval-censored +# frequency tables) AND at least MIN_DATASETS_PER_ARM datasets of summary- +# statistic format (types 1/2/3: median+range, median+IQR, mean+SD), the +# script fits three model versions: +# +# A. Individual-level only — summary types 4 and 5 only +# B. Summary-statistics only — summary types 1, 2, and 3 only +# C. Federated / full model — all available data formats jointly +# +# Comparison metrics (per pathogen × distribution): +# +# Posterior summaries (all three arms) +# pred_median posterior predictive median (days): median (2.5%, 97.5%) +# pred_q95 posterior predictive 95th percentile (days) +# tau between-study SD: median (2.5%, 97.5%); "— (n<5)" when sparse +# +# Uncertainty +# width 95% CrI width of pred_median (posterior uncertainty) +# interval_ratio A_width / C_width and B_width / C_width; > 1 = federated is +# tighter +# norm_shift (C_median − arm_median) / arm_width; normalised shift in units +# of arm's own uncertainty (|shift| > 0.5 = meaningful pull) +# +# Distributional discrepancy (posterior-median predictive CDFs) +# W2 Wasserstein-2 distance (days); quantile-function L2 over [0.1%, 99.9%] +# JS Jensen-Shannon divergence (bits, base-2); symmetric, bounded [0, 1] +# OVL Overlap coefficient ∫ min(f_A, f_C) dx; 1 = identical, 0 = disjoint +# +# Tail probabilities P(T > d | arm) at TAIL_THRESHOLDS (days) +# ptail_d_A, ptail_d_B, ptail_d_C (one triplet per threshold) +# +# LOO/ELPD +# elpd / elpd_se study-level PSIS-LOO ELPD (NA when n_datasets < 4) +# n_k_high number of studies with Pareto k̂ > 0.7 (LOO reliability flag) +# +# Differences (federated minus single-arm, posterior-median point estimates) +# delta_CA_median, delta_CA_q95 +# delta_CB_median, delta_CB_q95 +# +# LOO/ELPD notes +# -------------- +# log_lik[n_datasets * 3] is generated by the Stan model. For types 1/2, the +# three entries per study are the individual order-statistic contributions; for +# types 3/4/5 only the first entry is used (unused slots are 0). Summing each +# triplet gives a per-study log-likelihood matrix for study-level PSIS-LOO. +# ELPD values are comparable within an arm across distributions but NOT across +# arms — each arm conditions on a different number and type of observations. +# +# Distributional-metric notes +# --------------------------- +# W2, JS, and OVL are computed from the posterior-median predictive CDF, i.e. +# the CDF at the draw-wise median of the predictive distribution. This is a +# single representative distribution, not a posterior average; it is consistent +# with how pred_median and pred_q95 are reported. All three metrics compare +# two arms pairwise: (A, C), (B, C), and (A, B). +# +# Output +# ------ +# results/data_format_ablation.rds — nested list of stanfit objects + tables +# results/data_format_ablation.csv — main comparison table (wide format) +# +# Results file format: +# $ablation_fits[[pathogen]][[arm]][[dist]] +# $fit rstan stanfit object (NULL if sampling failed) +# $stan_data list passed to Stan +# $datasets dataset list used for this arm +# $skipped TRUE + $reason when deliberately skipped +# $comparison_tbl data.frame with all comparison columns +# ============================================================================= + +devtools::load_all(here::here(), quiet = TRUE) + +local({ + pkg_datasets <- grep("^datasets_", ls(asNamespace("ddsynth")), value = TRUE) + stale <- intersect(ls(envir = .GlobalEnv), pkg_datasets) + if (length(stale) > 0L) { + message("Removing ", length(stale), " stale dataset object(s) from .GlobalEnv: ", + paste(stale, collapse = ", ")) + rm(list = stale, envir = .GlobalEnv) + } +}) + +library(rstan) +rstan_options(auto_write = TRUE) +options(mc.cores = parallel::detectCores()) + + +# ── 1. Settings ─────────────────────────────────────────────────────────────── + +CHAINS <- 4 +ITER <- 12000 +WARMUP <- 2000 +THIN <- 1 +SEED <- 123 +CONTROL <- list(adapt_delta = 0.999, max_treedepth = 12, stepsize = 0.01) + +DIST_CODES <- c(lognormal = 1L, gamma = 2L, weibull = 3L, burr = 4L, gengamma = 5L) + +# Minimum datasets required in both arms before fitting. +MIN_DATASETS_PER_ARM <- 2L + +# Tail-probability thresholds (days) for P(T > d | arm). +TAIL_THRESHOLDS <- c(7, 14, 21) + +# Set FALSE to skip distributional metrics (W2/JS/OVL/tail) after fitting — +# useful for a quick tabular run when the CDFs are expensive to compute. +COMPUTE_CDF_METRICS <- TRUE + +# Number of posterior draws used inside compute_predictive_cdf(). +# Higher = more accurate; 200 is a reasonable trade-off for a comparison run. +N_CDF_DRAWS <- 200 + +# Pathogens to wipe and refit from scratch. +FORCE_RERUN <- character(0) + +CI_PROBS <- c(0.025, 0.5, 0.975) + + +# ── 2. Output paths ─────────────────────────────────────────────────────────── + +OUTPUT_DIR <- here::here("results") +OUTPUT_FILE <- file.path(OUTPUT_DIR, "data_format_ablation.rds") +CSV_FILE <- file.path(OUTPUT_DIR, "data_format_ablation.csv") +dir.create(OUTPUT_DIR, showWarnings = FALSE) + + +# ── 3. Pathogen registry ────────────────────────────────────────────────────── + +pathogen_registry <- list( + Nipah = datasets_Nipah, + MVD = datasets_MVD, + EVD = datasets_EVD, + Lassa = datasets_Lassa, + SARS = datasets_SARS, + MERS = datasets_MERS, + Zika = datasets_Zika, + Measles = datasets_Measles, + Mpox = datasets_Mpox, + Cholera = datasets_Cholera, + RVF = datasets_RVF, + CCHF = datasets_CCHF_extended, + COVID_19 = datasets_COVID_19, + Dengue = datasets_Dengue, + YFV = datasets_YFV, + Typhoid = datasets_typhoid, + Smallpox = datasets_Smallpox, + Flu = datasets_flu +) + + +# ── 4. Dataset classification helpers ──────────────────────────────────────── + +.get_summary_type <- function(d) { + if (!is.null(d$freq_lower) && !is.null(d$freq_upper) && !is.null(d$freq_count)) 5L + else if (!is.null(d$freq_value) && !is.null(d$freq_count)) 4L + else if (!is.null(d$mean) && !is.null(d$sd)) 3L + else if (!is.null(d$median) && !is.null(d$Q1) && !is.null(d$Q3)) 2L + else if (!is.null(d$median) && (!is.null(d$min) || !is.null(d$max))) 1L + else NA_integer_ +} + +.make_arms <- function(datasets) { + types <- vapply(datasets, .get_summary_type, integer(1L)) + list( + individual_only = datasets[!is.na(types) & types %in% c(4L, 5L)], + summary_only = datasets[!is.na(types) & types %in% c(1L, 2L, 3L)], + federated = datasets + ) +} + +# compute_predictive_cdf() uses "gg" for generalised gamma; all others match. +.cdf_dist_name <- function(dist_name) { + if (dist_name == "gengamma") "gg" else dist_name +} + + +# ── 5. Identify eligible pathogens ─────────────────────────────────────────── + +eligible_pathogens <- Filter(function(ds) { + arms <- .make_arms(ds) + length(arms$individual_only) >= MIN_DATASETS_PER_ARM && + length(arms$summary_only) >= MIN_DATASETS_PER_ARM +}, pathogen_registry) + +if (length(eligible_pathogens) == 0L) + stop("No pathogens meet the eligibility criterion (>= ", MIN_DATASETS_PER_ARM, + " datasets per arm). Reduce MIN_DATASETS_PER_ARM.") + +message("\n", strrep("=", 70)) +message("Eligible pathogens for data-format ablation") +message("(require >= ", MIN_DATASETS_PER_ARM, " datasets per arm):") +message(sprintf(" %-12s %9s %10s %7s", + "Pathogen", "indiv (A)", "sumstat (B)", "total (C)")) +message(strrep("-", 46)) +for (p in names(eligible_pathogens)) { + arms <- .make_arms(eligible_pathogens[[p]]) + message(sprintf(" %-12s %9d %10d %7d", + p, length(arms$individual_only), + length(arms$summary_only), + length(arms$federated))) +} +message(strrep("=", 70)) + + +# ── 6. Compile Stan model ───────────────────────────────────────────────────── + +stan_file <- system.file("stan", + "hierarchical_data_synthesis_summary_stats.stan", + package = "ddsynth") +if (!nzchar(stan_file)) + stop("Stan file 'hierarchical_data_synthesis_summary_stats.stan' not found.") +stan_model <- rstan::stan_model(stan_file) + + +# ── 7. Fitting helper ───────────────────────────────────────────────────────── + +.fit_one <- function(datasets, dist_name, arm, pathogen) { + dist_type <- DIST_CODES[[dist_name]] + + stan_data <- tryCatch( + prepare_stan_data_from_datasets(datasets, dist_type = dist_type), + error = function(e) { + message(" [SKIP] prepare_stan_data failed for ", + pathogen, "/", arm, "/", dist_name, ": ", conditionMessage(e)) + NULL + } + ) + if (is.null(stan_data)) return(NULL) + stan_data <- update_phi_prior(stan_data, datasets) + + fit <- tryCatch( + rstan::sampling( + stan_model, + data = stan_data, + chains = CHAINS, + iter = ITER, + warmup = WARMUP, + thin = THIN, + control = CONTROL, + seed = SEED, + refresh = 200 + ), + error = function(e) { + message(" [FAIL] sampling failed for ", + pathogen, "/", arm, "/", dist_name, ": ", conditionMessage(e)) + NULL + } + ) + + list(fit = fit, stan_data = stan_data, datasets = datasets) +} + + +# ── 8. Load stored results ─────────────────────────────────────────────────── + +ablation_results <- if (file.exists(OUTPUT_FILE)) { + message("Resuming from existing results: ", OUTPUT_FILE) + readRDS(OUTPUT_FILE)$ablation_fits +} else { + list() +} + +if (length(FORCE_RERUN) > 0L) { + unknown <- setdiff(FORCE_RERUN, names(eligible_pathogens)) + if (length(unknown) > 0L) + warning("FORCE_RERUN contains unrecognised pathogen(s): ", + paste(unknown, collapse = ", "), call. = FALSE) + for (p in intersect(FORCE_RERUN, names(ablation_results))) { + message("Force-rerun: clearing stored results for '", p, "'.") + ablation_results[[p]] <- NULL + } +} + + +# ── 9. Main fitting loop ────────────────────────────────────────────────────── + +ARM_LABELS <- c("individual_only", "summary_only", "federated") + +for (pathogen in names(eligible_pathogens)) { + + message("\n", strrep("=", 70)) + message("PATHOGEN: ", pathogen) + message(strrep("=", 70)) + + if (is.null(ablation_results[[pathogen]])) + ablation_results[[pathogen]] <- list() + + arms <- .make_arms(eligible_pathogens[[pathogen]]) + + all_done <- all(vapply(ARM_LABELS, function(arm) { + all(vapply(names(DIST_CODES), function(d) { + !is.null(ablation_results[[pathogen]][[arm]][[d]]) + }, logical(1L))) + }, logical(1L))) + + if (all_done) { + message(" All fits already complete — nothing new to do.") + next + } + + for (arm in ARM_LABELS) { + + arm_datasets <- arms[[arm]] + + if (is.null(ablation_results[[pathogen]][[arm]])) + ablation_results[[pathogen]][[arm]] <- list() + + n_arm <- length(arm_datasets) + if (n_arm < MIN_DATASETS_PER_ARM) { + message("\n [SKIP] ", arm, ": only ", n_arm, + " dataset(s) (min = ", MIN_DATASETS_PER_ARM, ").") + next + } + + message("\n Arm: ", arm, " (n_datasets = ", n_arm, ")") + + for (dist_name in names(DIST_CODES)) { + + if (!is.null(ablation_results[[pathogen]][[arm]][[dist_name]])) { + message(" [SKIP] ", pathogen, " / ", arm, " / ", dist_name, + " — already in results.") + next + } + + if (dist_name == "gengamma" && + !should_attempt_gg(arm_datasets, verbose = FALSE)) { + message(" [SKIP] ", pathogen, " / ", arm, + " / gengamma — identifiability heuristic.") + ablation_results[[pathogen]][[arm]][["gengamma"]] <- + list(skipped = TRUE, reason = "identifiability_heuristic", + datasets = arm_datasets) + next + } + + if (dist_name == "gamma" && + !gamma_type2_reliable(arm_datasets, verbose = FALSE)) { + message(" [SKIP] ", pathogen, " / ", arm, + " / gamma — gamma_type2_reliable check.") + ablation_results[[pathogen]][[arm]][["gamma"]] <- + list(skipped = TRUE, reason = "gamma_type2_reliable", + datasets = arm_datasets) + next + } + + message(" Fitting: ", pathogen, " | ", arm, " | ", dist_name) + ablation_results[[pathogen]][[arm]][[dist_name]] <- + .fit_one(arm_datasets, dist_name, arm, pathogen) + } + } + + saveRDS(list(ablation_fits = ablation_results), OUTPUT_FILE) + message("\n Saved to: ", OUTPUT_FILE) +} + +message("\n", strrep("=", 70)) +message("Fitting complete.") +message(strrep("=", 70)) + + +# ── 10. Summary extraction helpers ─────────────────────────────────────────── + +# Study-level PSIS-LOO ELPD using the log_lik triplets from Stan. +# Returns list(elpd, se, n_k_high) — all NA/0 on failure or when n < 4. +.compute_loo_elpd <- function(fit, n_datasets) { + na_res <- list(elpd = NA_real_, se = NA_real_, n_k_high = NA_integer_) + if (is.null(fit) || n_datasets < 4L) return(na_res) + + ll_raw <- tryCatch( + rstan::extract(fit, pars = "log_lik")$log_lik, + error = function(e) NULL + ) + if (is.null(ll_raw) || ncol(ll_raw) != n_datasets * 3L) return(na_res) + + n_draws <- nrow(ll_raw) + ll_per_study <- matrix(0, nrow = n_draws, ncol = n_datasets) + for (d in seq_len(n_datasets)) { + base <- (d - 1L) * 3L + 1L + ll_per_study[, d] <- ll_raw[, base] + ll_raw[, base + 1L] + ll_raw[, base + 2L] + } + ll_per_study[!is.finite(ll_per_study)] <- -700 + + loo_res <- tryCatch( + loo::loo(ll_per_study, r_eff = rep(1, n_datasets)), + error = function(e) NULL + ) + if (is.null(loo_res)) return(na_res) + + n_k_high <- sum(loo_res$diagnostics$pareto_k > 0.7, na.rm = TRUE) + + list( + elpd = loo_res$estimates["elpd_loo", "Estimate"], + se = loo_res$estimates["elpd_loo", "SE"], + n_k_high = n_k_high + ) +} + +# Format a posterior sample vector as "median (lo, hi)". +.fmt <- function(x, digits = 2L) { + x <- x[is.finite(x)] + if (length(x) == 0L) return(NA_character_) + q <- quantile(x, CI_PROBS) + sprintf("%.*f (%.*f, %.*f)", digits, q[2L], digits, q[1L], digits, q[3L]) +} + +# Extract posterior summaries and LOO from one arm's result slot. +.extract_summaries <- function(res) { + if (is.null(res) || isTRUE(res$skipped) || is.null(res$fit)) return(NULL) + sims <- tryCatch(rstan::extract(res$fit), error = function(e) NULL) + if (is.null(sims)) return(NULL) + + n_datasets <- length(res$datasets) + pm <- sims$pred_median[is.finite(sims$pred_median)] + pq <- sims$pred_q95[is.finite(sims$pred_q95)] + tau_raw <- if (n_datasets >= 5L) sims$tau else NULL + loo_res <- .compute_loo_elpd(res$fit, n_datasets) + + list( + pred_median_fmt = .fmt(pm, 1L), + pred_q95_fmt = .fmt(pq, 1L), + pred_median_med = median(pm), + pred_q95_med = median(pq), + pred_median_width = diff(quantile(pm, c(0.025, 0.975))), + tau_fmt = if (!is.null(tau_raw)) .fmt(tau_raw, 2L) else "— (n<5)", + tau_med = if (!is.null(tau_raw)) median(tau_raw) else NA_real_, + elpd = loo_res$elpd, + elpd_se = loo_res$se, + n_k_high = loo_res$n_k_high, + n_datasets = n_datasets + ) +} + +# Compute distributional metrics between two posterior-median predictive CDFs. +# +# Uses the posterior-median CDF (the draw whose CDF is closest to the +# median across all draws) as a single representative distribution. W2 is +# in the same units as the outcome (days); JS is in bits (base-2). OVL is +# the overlap coefficient ∫ min(f_ref, f_arm) dx ∈ [0, 1]. +# +# Also returns tail probabilities P(T > d) at each threshold in thr_days. +# Returns NULL when either fit is absent or CDF computation fails. +.compute_dist_metrics <- function(fit_ref, fit_arm, dist_name, + thr_days = TAIL_THRESHOLDS) { + if (is.null(fit_ref) || is.null(fit_arm)) return(NULL) + + cdf_name <- .cdf_dist_name(dist_name) + + # Determine x upper bound from the upper tail of pred_q95 across both fits. + .q99_pred_q95 <- function(f) { + s <- tryCatch(rstan::extract(f, pars = "pred_q95")$pred_q95, error = function(e) NULL) + if (is.null(s) || length(s) == 0L) return(30) + quantile(s[is.finite(s)], 0.99, na.rm = TRUE) + } + x_max <- max(.q99_pred_q95(fit_ref), .q99_pred_q95(fit_arm)) * 1.1 + x_max <- max(x_max, max(thr_days) * 1.2, 14) + + x_seq <- seq(0.001, x_max, length.out = 1000L) + + cdf_ref_obj <- tryCatch( + compute_predictive_cdf(fit_ref, cdf_name, x_seq = x_seq, n_draws = N_CDF_DRAWS), + error = function(e) NULL + ) + cdf_arm_obj <- tryCatch( + compute_predictive_cdf(fit_arm, cdf_name, x_seq = x_seq, n_draws = N_CDF_DRAWS), + error = function(e) NULL + ) + if (is.null(cdf_ref_obj) || is.null(cdf_arm_obj)) return(NULL) + + # Posterior-median CDF: clamp to [ε, 1−ε] for numerical stability. + F_ref <- pmax(pmin(cdf_ref_obj$summary$median, 1 - 1e-9), 1e-9) + F_arm <- pmax(pmin(cdf_arm_obj$summary$median, 1 - 1e-9), 1e-9) + + dx <- diff(x_seq)[1L] + + # PDF via forward differences; renormalise to guard against CDF + # non-monotonicity from interpolation artefacts. + .cdf_to_pdf <- function(F) { + f <- c(diff(F), 0) / dx + f <- pmax(f, 1e-12) + f / (sum(f) * dx) + } + f_ref <- .cdf_to_pdf(F_ref) + f_arm <- .cdf_to_pdf(F_arm) + + # Jensen-Shannon divergence (bits, symmetric, bounded [0, 1]) + f_mix <- 0.5 * (f_ref + f_arm) + JS <- 0.5 * sum(f_ref * log2(f_ref / f_mix) * dx, na.rm = TRUE) + + 0.5 * sum(f_arm * log2(f_arm / f_mix) * dx, na.rm = TRUE) + JS <- pmax(JS, 0) # guard against tiny negative values from numerics + + # Overlap coefficient + OVL <- sum(pmin(f_ref, f_arm) * dx, na.rm = TRUE) + + # Wasserstein-2 via quantile-function L2 norm + # Q(p) = inf{x : F(x) >= p} approximated by linear interpolation. + p_grid <- seq(0.001, 0.999, length.out = 500L) + .quantile_fun <- function(F_vals) { + sapply(p_grid, function(p) { + idx <- which(F_vals >= p)[1L] + if (is.na(idx)) return(x_seq[length(x_seq)]) + if (idx == 1L) return(x_seq[1L]) + # Linear interpolation between grid points + x_seq[idx - 1L] + + (p - F_vals[idx - 1L]) / (F_vals[idx] - F_vals[idx - 1L]) * + (x_seq[idx] - x_seq[idx - 1L]) + }) + } + Q_ref <- .quantile_fun(F_ref) + Q_arm <- .quantile_fun(F_arm) + W2 <- sqrt(mean((Q_ref - Q_arm)^2, na.rm = TRUE)) + + # Tail probabilities P(T > d) = 1 - F(d) at each threshold. + .tail_prob <- function(F_vals, d) { + idx <- which(x_seq >= d)[1L] + if (is.na(idx)) return(0) + 1 - F_vals[idx] + } + tail_ref <- setNames(sapply(thr_days, .tail_prob, F_vals = F_ref), paste0("d", thr_days)) + tail_arm <- setNames(sapply(thr_days, .tail_prob, F_vals = F_arm), paste0("d", thr_days)) + + list( + W2 = round(W2, 3L), + JS = round(JS, 4L), + OVL = round(OVL, 4L), + tail_ref = tail_ref, + tail_arm = tail_arm + ) +} + + +# ── 11. Build comparison table ──────────────────────────────────────────────── + +dist_labels <- c( + lognormal = "Log-normal", gamma = "Gamma", weibull = "Weibull", + burr = "Burr XII", gengamma = "Gen. gamma" +) + +rows <- list() + +for (pathogen in names(eligible_pathogens)) { + for (dist_name in names(DIST_CODES)) { + + res_A <- ablation_results[[pathogen]][["individual_only"]][[dist_name]] + res_B <- ablation_results[[pathogen]][["summary_only"]][[dist_name]] + res_C <- ablation_results[[pathogen]][["federated"]][[dist_name]] + + s_A <- .extract_summaries(res_A) + s_B <- .extract_summaries(res_B) + s_C <- .extract_summaries(res_C) + + if (is.null(s_C)) next # arm C is the reference; skip row if absent + + # ── Distributional metrics (if CDF computation is enabled) ──────────────── + dm_CA <- dm_CB <- dm_AB <- NULL + if (COMPUTE_CDF_METRICS && + !is.null(res_C) && !isTRUE(res_C$skipped) && !is.null(res_C$fit)) { + + message(" Computing CDF metrics: ", pathogen, " / ", dist_name) + + if (!is.null(res_A) && !isTRUE(res_A$skipped) && !is.null(res_A$fit)) + dm_CA <- .compute_dist_metrics(res_C$fit, res_A$fit, dist_name) + + if (!is.null(res_B) && !isTRUE(res_B$skipped) && !is.null(res_B$fit)) + dm_CB <- .compute_dist_metrics(res_C$fit, res_B$fit, dist_name) + + if (!is.null(res_A) && !isTRUE(res_A$skipped) && !is.null(res_A$fit) && + !is.null(res_B) && !isTRUE(res_B$skipped) && !is.null(res_B$fit)) + dm_AB <- .compute_dist_metrics(res_A$fit, res_B$fit, dist_name) + } + + # ── Helpers ─────────────────────────────────────────────────────────────── + .n <- function(res) { + if (is.null(res) || isTRUE(res$skipped)) NA_integer_ + else length(res$datasets) + } + .delta <- function(s_ref, s_arm, field) { + if (is.null(s_ref) || is.null(s_arm)) NA_real_ + else round(s_ref[[field]] - s_arm[[field]], 2L) + } + .ratio <- function(s_arm, s_ref) { + # arm_width / ref_width; > 1 means arm is wider (federated is tighter) + if (is.null(s_arm) || is.null(s_ref)) NA_real_ + else round(s_arm$pred_median_width / s_ref$pred_median_width, 3L) + } + .norm_shift <- function(s_arm, s_ref) { + # (ref_median − arm_median) / arm_width + if (is.null(s_arm) || is.null(s_ref)) NA_real_ + else round((s_ref$pred_median_med - s_arm$pred_median_med) / + s_arm$pred_median_width, 3L) + } + .dm_field <- function(dm, field) { + if (is.null(dm)) NA_real_ else dm[[field]] + } + .tail <- function(dm, side, threshold) { + # side = "tail_ref" (arm C) or "tail_arm" (the other arm) + key <- paste0("d", threshold) + if (is.null(dm) || is.null(dm[[side]])) NA_real_ + else round(dm[[side]][[key]], 4L) + } + + # ── Tail probability columns ─────────────────────────────────────────────── + tail_cols <- list() + for (d in TAIL_THRESHOLDS) { + tail_cols[[paste0("ptail_A_", d)]] <- .tail(dm_CA, "tail_arm", d) + tail_cols[[paste0("ptail_B_", d)]] <- .tail(dm_CB, "tail_arm", d) + tail_cols[[paste0("ptail_C_", d)]] <- .tail(dm_CA, "tail_ref", d) + } + + row <- c( + list( + pathogen = pathogen, + dist = dist_labels[[dist_name]], + n_indiv = .n(res_A), + n_sumstat = .n(res_B), + n_total = .n(res_C), + + # ── Arm A: individual-level only ───────────────────────────────────── + A_pred_median = if (!is.null(s_A)) s_A$pred_median_fmt else NA_character_, + A_pred_q95 = if (!is.null(s_A)) s_A$pred_q95_fmt else NA_character_, + A_tau = if (!is.null(s_A)) s_A$tau_fmt else NA_character_, + A_width = if (!is.null(s_A)) round(s_A$pred_median_width, 3L) else NA_real_, + A_elpd = if (!is.null(s_A) && !is.na(s_A$elpd)) round(s_A$elpd, 1L) else NA_real_, + A_elpd_se = if (!is.null(s_A) && !is.na(s_A$elpd_se)) round(s_A$elpd_se, 1L) else NA_real_, + A_n_k_high = if (!is.null(s_A)) s_A$n_k_high else NA_integer_, + + # ── Arm B: summary-statistics only ────────────────────────────────── + B_pred_median = if (!is.null(s_B)) s_B$pred_median_fmt else NA_character_, + B_pred_q95 = if (!is.null(s_B)) s_B$pred_q95_fmt else NA_character_, + B_tau = if (!is.null(s_B)) s_B$tau_fmt else NA_character_, + B_width = if (!is.null(s_B)) round(s_B$pred_median_width, 3L) else NA_real_, + B_elpd = if (!is.null(s_B) && !is.na(s_B$elpd)) round(s_B$elpd, 1L) else NA_real_, + B_elpd_se = if (!is.null(s_B) && !is.na(s_B$elpd_se)) round(s_B$elpd_se, 1L) else NA_real_, + B_n_k_high = if (!is.null(s_B)) s_B$n_k_high else NA_integer_, + + # ── Arm C: federated ───────────────────────────────────────────────── + C_pred_median = s_C$pred_median_fmt, + C_pred_q95 = s_C$pred_q95_fmt, + C_tau = s_C$tau_fmt, + C_width = round(s_C$pred_median_width, 3L), + C_elpd = if (!is.na(s_C$elpd)) round(s_C$elpd, 1L) else NA_real_, + C_elpd_se = if (!is.na(s_C$elpd_se)) round(s_C$elpd_se, 1L) else NA_real_, + C_n_k_high = s_C$n_k_high, + + # ── Differences: federated minus each arm ─────────────────────────── + delta_CA_median = .delta(s_C, s_A, "pred_median_med"), + delta_CA_q95 = .delta(s_C, s_A, "pred_q95_med"), + delta_CB_median = .delta(s_C, s_B, "pred_median_med"), + delta_CB_q95 = .delta(s_C, s_B, "pred_q95_med"), + + # ── Interval ratio: arm_width / C_width (> 1 = federated is tighter) ── + ratio_CA = .ratio(s_A, s_C), + ratio_CB = .ratio(s_B, s_C), + + # ── Normalised shift: (C_median − arm_median) / arm_width ────────── + norm_shift_CA = .norm_shift(s_A, s_C), + norm_shift_CB = .norm_shift(s_B, s_C), + + # ── Distributional metrics (A vs C) ────────────────────────────────── + W2_CA = .dm_field(dm_CA, "W2"), + JS_CA = .dm_field(dm_CA, "JS"), + OVL_CA = .dm_field(dm_CA, "OVL"), + + # ── Distributional metrics (B vs C) ────────────────────────────────── + W2_CB = .dm_field(dm_CB, "W2"), + JS_CB = .dm_field(dm_CB, "JS"), + OVL_CB = .dm_field(dm_CB, "OVL"), + + # ── Distributional metrics (A vs B) ────────────────────────────────── + W2_AB = .dm_field(dm_AB, "W2"), + JS_AB = .dm_field(dm_AB, "JS"), + OVL_AB = .dm_field(dm_AB, "OVL") + ), + tail_cols + ) + + rows[[length(rows) + 1L]] <- as.data.frame(row, stringsAsFactors = FALSE) + } +} + +comparison_tbl <- if (length(rows) > 0L) { + tbl <- do.call(rbind, rows) + rownames(tbl) <- NULL + tbl +} else { + message("WARNING: no rows in comparison table — check fitting completed.") + data.frame() +} + + +# ── 12. Save and report ─────────────────────────────────────────────────────── + +saveRDS( + list(ablation_fits = ablation_results, comparison_tbl = comparison_tbl), + OUTPUT_FILE +) +message("\nRDS saved to: ", OUTPUT_FILE) + +if (nrow(comparison_tbl) > 0L) { + write.csv(comparison_tbl, CSV_FILE, row.names = FALSE) + message("CSV saved to: ", CSV_FILE) +} + +message("\n", strrep("=", 80)) +message("DATA-FORMAT ABLATION — summary tables") +message(strrep("=", 80)) +message( + "\nArms:\n", + " A individual_only — types 4/5: exact + interval-censored freq tables\n", + " B summary_only — types 1/2/3: median+range, median+IQR, mean+SD\n", + " C federated — all data formats jointly\n" +) + +if (nrow(comparison_tbl) > 0L) { + + # ── Table 1: Predictive summaries ─────────────────────────────────────────── + message("\n=== Table 1: Predictive summaries (days, median [95% CrI]) ===\n") + cols1 <- c("pathogen", "dist", "n_indiv", "n_sumstat", "n_total", + "A_pred_median", "B_pred_median", "C_pred_median", + "delta_CA_median", "delta_CB_median") + print(comparison_tbl[, intersect(cols1, names(comparison_tbl))], + row.names = FALSE, right = FALSE) + + message("\n=== Table 2: Predictive 95th percentile (days) ===\n") + cols2 <- c("pathogen", "dist", + "A_pred_q95", "B_pred_q95", "C_pred_q95", + "delta_CA_q95", "delta_CB_q95") + print(comparison_tbl[, intersect(cols2, names(comparison_tbl))], + row.names = FALSE, right = FALSE) + + # ── Table 2: Uncertainty and tau ───────────────────────────────────────────── + message("\n=== Table 3: Posterior uncertainty and heterogeneity tau ===\n") + message(" width = 95% CrI width of pred_median\n", + " ratio_CA = A_width / C_width (> 1 = federated is tighter)\n", + " norm_shift = (C_median − arm_median) / arm_width\n") + cols3 <- c("pathogen", "dist", + "A_width", "B_width", "C_width", + "ratio_CA", "ratio_CB", + "norm_shift_CA", "norm_shift_CB", + "A_tau", "B_tau", "C_tau") + print(comparison_tbl[, intersect(cols3, names(comparison_tbl))], + row.names = FALSE, right = FALSE) + + # ── Table 3: Distributional metrics ────────────────────────────────────────── + if (COMPUTE_CDF_METRICS) { + message("\n=== Table 4: Distributional discrepancy metrics ===\n") + message(" W2 Wasserstein-2 (days); CA = arm C vs A, CB = C vs B, AB = A vs B\n", + " JS Jensen-Shannon divergence (bits); symmetric, bounded [0, 1]\n", + " OVL Overlap coefficient; 1 = identical distributions\n") + cols4 <- c("pathogen", "dist", + "W2_CA", "JS_CA", "OVL_CA", + "W2_CB", "JS_CB", "OVL_CB", + "W2_AB", "JS_AB", "OVL_AB") + print(comparison_tbl[, intersect(cols4, names(comparison_tbl))], + row.names = FALSE, right = FALSE) + + message("\n=== Table 5: Tail probabilities P(T > d days) ===\n") + tail_col_names <- c("pathogen", "dist", + unlist(lapply(TAIL_THRESHOLDS, function(d) + c(paste0("ptail_A_", d), + paste0("ptail_B_", d), + paste0("ptail_C_", d))))) + print(comparison_tbl[, intersect(tail_col_names, names(comparison_tbl))], + row.names = FALSE, right = FALSE) + } + + # ── Table 4: LOO / ELPD ─────────────────────────────────────────────────────── + message("\n=== Table 6: Study-level PSIS-LOO ELPD ===\n") + message(" n_k_high = studies with Pareto k̂ > 0.7 (LOO reliability flag)\n") + cols6 <- c("pathogen", "dist", + "A_elpd", "A_elpd_se", "A_n_k_high", + "B_elpd", "B_elpd_se", "B_n_k_high", + "C_elpd", "C_elpd_se", "C_n_k_high") + print(comparison_tbl[, intersect(cols6, names(comparison_tbl))], + row.names = FALSE, right = FALSE) +} + +message("\n", strrep("=", 80)) +message("Done. Results written to: ", OUTPUT_FILE) +message(strrep("=", 80)) From 7587da73858054f15693dbaec4888f34f470d2ea Mon Sep 17 00:00:00 2001 From: cm401 Date: Sun, 10 May 2026 10:32:19 +0800 Subject: [PATCH 02/18] Add plotting script for data-format ablation results Three-figure suite: (1) three-arm forest plot of pred_median with 95% CrI, (2) federated-gain metrics strip (interval ratio, JS divergence, OVL), (3) posterior predictive density overlays for top pathogens by JS divergence. Supplementary version of figure 2 faceted across all distributions. Co-Authored-By: Claude Sonnet 4.6 --- analysis/plot_data_format_ablation.R | 554 +++++++++++++++++++++++++++ 1 file changed, 554 insertions(+) create mode 100644 analysis/plot_data_format_ablation.R diff --git a/analysis/plot_data_format_ablation.R b/analysis/plot_data_format_ablation.R new file mode 100644 index 0000000..652fcd5 --- /dev/null +++ b/analysis/plot_data_format_ablation.R @@ -0,0 +1,554 @@ +# ============================================================================= +# plot_data_format_ablation.R +# ----------------------------------------------------------------------------- +# Visualises the data-format ablation results produced by +# data_format_ablation.R. +# +# Three figures, each answering a distinct question: +# +# Figure 1 — "Do the arms agree?" +# Three-arm forest plot. Pathogens on the y-axis, pred_median on a +# log-scaled x-axis. Arms A, B, C dodged vertically per pathogen, with +# horizontal 95% CrI bars. A vertical reference line marks the federated +# (arm C) point estimate. Faceted by distribution family. +# +# Figure 2 — "How much does federation help?" +# Side-by-side metrics strip using the same pathogen ordering. +# Left panel : interval ratio (A_width / C_width, B_width / C_width); +# reference line at 1 (equal precision). +# Middle panel: Jensen-Shannon divergence (JS_CA, JS_CB; bits). +# Right panel : overlap coefficient (OVL_CA, OVL_CB). +# All three panels share the y-axis (pathogens) and are assembled with +# patchwork. Restricted to one distribution (PRIMARY_DIST) for clarity; +# a supplementary version facets over all distributions. +# +# Figure 3 — "What does the gain look like in practice?" +# Posterior predictive density overlays for the TOP_N_DENSITY pathogens +# with the largest JS_CA divergence (i.e. where individual-level data +# pulls the federated estimate furthest from the summary-stat arm). +# Three shaded ribbons (A, B, C) per pathogen panel, computed from +# compute_predictive_cdf(). +# +# Output +# ------ +# results/figures/fig_ablation_forest.pdf (+ .png) +# results/figures/fig_ablation_gain.pdf (+ .png) +# results/figures/fig_ablation_densities.pdf (+ .png) +# ============================================================================= + +devtools::load_all(here::here(), quiet = TRUE) +library(ggplot2) +library(dplyr) +library(tidyr) +library(patchwork) + +# ── 1. Settings ─────────────────────────────────────────────────────────────── + +# Distribution shown in Figure 2 (gain strip) and Figure 3 (densities). +PRIMARY_DIST <- "Log-normal" + +# Number of pathogens shown in Figure 3, chosen by largest JS_CA divergence. +TOP_N_DENSITY <- 6L + +# Save dimensions (inches) and resolution. +FIG_WIDTH_WIDE <- 14 +FIG_WIDTH_HALF <- 7 +FIG_HEIGHT_ROW <- 0.55 # height per pathogen row in Figures 1 & 2 +FIG_DPI <- 300 + +# Arm colours (Wong colorblind-safe palette). +ARM_COLOURS <- c( + "A" = "#0072B2", # blue + "B" = "#D55E00", # vermilion + "C" = "#009E73" # green +) +ARM_LABELS <- c( + "A" = "Individual-level only (A)", + "B" = "Summary-statistics only (B)", + "C" = "Federated (C)" +) +ARM_SHAPES <- c("A" = 19, "B" = 17, "C" = 18) # circle, triangle, diamond + +OUTPUT_DIR <- here::here("results", "figures") +dir.create(OUTPUT_DIR, showWarnings = FALSE, recursive = TRUE) + + +# ── 2. Load results ─────────────────────────────────────────────────────────── + +rds_file <- here::here("results", "data_format_ablation.rds") +if (!file.exists(rds_file)) + stop("Results file not found: ", rds_file, + "\nRun analysis/data_format_ablation.R first.") + +stored <- readRDS(rds_file) +comparison_tbl <- stored$comparison_tbl +ablation_fits <- stored$ablation_fits + +if (nrow(comparison_tbl) == 0L) + stop("comparison_tbl is empty — check that data_format_ablation.R completed.") + + +# ── 3. Shared helpers ───────────────────────────────────────────────────────── + +# Parse "9.1 (8.2, 10.3)" -> named vector c(med, lo, hi). +.parse_cri <- function(s) { + if (is.na(s) || grepl("—|NA", s)) + return(c(med = NA_real_, lo = NA_real_, hi = NA_real_)) + nums <- suppressWarnings( + as.numeric(unlist(regmatches(s, gregexpr("[0-9]+\\.?[0-9]*", s)))) + ) + if (length(nums) >= 3L) c(med = nums[1L], lo = nums[2L], hi = nums[3L]) + else c(med = NA_real_, lo = NA_real_, hi = NA_real_) +} + +# Consistent theme for all figures. +theme_ablation <- function(base_size = 10) { + theme_minimal(base_size = base_size) + + theme( + panel.grid.major.y = element_blank(), + panel.grid.minor = element_blank(), + strip.text = element_text(face = "bold"), + legend.position = "bottom", + legend.title = element_blank(), + axis.title.y = element_blank() + ) +} + +# Order pathogens by federated pred_median (log-normal, descending). +.pathogen_order <- function(tbl, dist = PRIMARY_DIST) { + tbl |> + filter(.data$dist == !!dist) |> + rowwise() |> + mutate(C_med = .parse_cri(C_pred_median)[["med"]]) |> + ungroup() |> + arrange(C_med) |> + pull(pathogen) |> + unique() +} + +pathogen_order <- .pathogen_order(comparison_tbl) + + +# ── 4. Figure 1: Three-arm forest plot ─────────────────────────────────────── +# +# For each arm, parse the formatted CrI string into med/lo/hi and reshape to +# long format. Dodge three arms vertically within each pathogen × distribution +# cell. Log-scale x-axis with a vertical reference line at arm C's median. + +.parse_arm_col <- function(tbl, col, arm_label) { + tbl |> + select(pathogen, dist, cri_str = all_of(col)) |> + rowwise() |> + mutate(parsed = list(.parse_cri(cri_str))) |> + ungroup() |> + mutate( + med = vapply(parsed, `[[`, numeric(1L), "med"), + lo = vapply(parsed, `[[`, numeric(1L), "lo"), + hi = vapply(parsed, `[[`, numeric(1L), "hi"), + arm = arm_label + ) |> + select(pathogen, dist, arm, med, lo, hi) +} + +forest_long <- bind_rows( + .parse_arm_col(comparison_tbl, "A_pred_median", "A"), + .parse_arm_col(comparison_tbl, "B_pred_median", "B"), + .parse_arm_col(comparison_tbl, "C_pred_median", "C") +) |> + filter(!is.na(med)) |> + mutate( + pathogen = factor(pathogen, levels = pathogen_order), + arm = factor(arm, levels = c("A", "B", "C")) + ) + +# Federated reference line (arm C) per facet cell. +ref_lines <- forest_long |> + filter(arm == "C") |> + select(pathogen, dist, ref_med = med) + +dodge_height <- 0.55 # vertical spread per pathogen row + +fig1 <- ggplot(forest_long, + aes(x = med, y = pathogen, + colour = arm, shape = arm)) + + # Reference line at arm C's point estimate + geom_vline(data = ref_lines, + aes(xintercept = ref_med), + colour = ARM_COLOURS[["C"]], linewidth = 0.25, + linetype = "dashed", inherit.aes = FALSE) + + # CrI bars + geom_errorbarh(aes(xmin = lo, xmax = hi), + height = 0, + linewidth = 0.5, + position = position_dodge(height = dodge_height)) + + # Point estimates + geom_point(size = 1.8, + position = position_dodge(height = dodge_height)) + + scale_x_log10( + breaks = c(1, 2, 5, 10, 20, 50), + labels = c("1", "2", "5", "10", "20", "50") + ) + + scale_colour_manual(values = ARM_COLOURS, labels = ARM_LABELS) + + scale_shape_manual(values = ARM_SHAPES, labels = ARM_LABELS) + + facet_wrap(~ dist, nrow = 1L, scales = "free_x") + + labs( + x = "Posterior predictive median (days, log scale)", + caption = "Horizontal bars: 95% credible intervals. Dashed line: federated (C) point estimate." + ) + + theme_ablation() + + guides(colour = guide_legend(override.aes = list(size = 2.5))) + +n_pathogens <- length(pathogen_order) +fig1_height <- max(4, n_pathogens * FIG_HEIGHT_ROW + 1.5) + +ggsave(file.path(OUTPUT_DIR, "fig_ablation_forest.pdf"), + fig1, width = FIG_WIDTH_WIDE, height = fig1_height, device = "pdf") +ggsave(file.path(OUTPUT_DIR, "fig_ablation_forest.png"), + fig1, width = FIG_WIDTH_WIDE, height = fig1_height, dpi = FIG_DPI) +message("Figure 1 saved.") + + +# ── 5. Figure 2: Federated-gain metrics strip ───────────────────────────────── +# +# Restricted to PRIMARY_DIST. Three ggplot objects assembled side-by-side with +# patchwork. Each panel shares the y-axis; only the left panel shows pathogen +# labels. + +gain_df <- comparison_tbl |> + filter(dist == PRIMARY_DIST) |> + mutate(pathogen = factor(pathogen, levels = pathogen_order)) |> + select(pathogen, + ratio_CA, ratio_CB, + JS_CA, JS_CB, + OVL_CA, OVL_CB) |> + pivot_longer( + cols = -pathogen, + names_to = c("metric", "comparison"), + names_sep = "_", + values_to = "value" + ) |> + filter(!is.na(value)) |> + mutate( + comparison = factor(comparison, + levels = c("CA", "CB"), + labels = c("vs A (individual-level)", "vs B (summary-stats)")), + metric = factor(metric, + levels = c("ratio", "JS", "OVL"), + labels = c("Interval ratio\n(arm / federated)", + "Jensen-Shannon\ndivergence (bits)", + "Overlap\ncoefficient")) + ) + +# Shared y scale +y_scale <- scale_y_discrete(drop = FALSE) + +# Helper: one panel +.gain_panel <- function(df, metric_label, x_lab, ref_val = NULL, + show_y = TRUE) { + d <- filter(df, metric == metric_label) + p <- ggplot(d, aes(x = value, y = pathogen, + colour = comparison, shape = comparison)) + + geom_point(size = 2.2, alpha = 0.85) + + y_scale + + labs(x = x_lab) + + scale_colour_manual( + values = c("vs A (individual-level)" = ARM_COLOURS[["A"]], + "vs B (summary-stats)" = ARM_COLOURS[["B"]]) + ) + + scale_shape_manual( + values = c("vs A (individual-level)" = 19, + "vs B (summary-stats)" = 17) + ) + + theme_ablation() + + theme(legend.position = if (show_y) "none" else "bottom") + + if (!is.null(ref_val)) + p <- p + geom_vline(xintercept = ref_val, + linetype = "dashed", colour = "grey50", + linewidth = 0.4) + if (!show_y) + p <- p + theme(axis.text.y = element_blank()) + p +} + +p2a <- .gain_panel(gain_df, + "Interval ratio\n(arm / federated)", + "Interval ratio (> 1 = federated is tighter)", + ref_val = 1, + show_y = TRUE) + +p2b <- .gain_panel(gain_df, + "Jensen-Shannon\ndivergence (bits)", + "JS divergence (bits)", + ref_val = 0, + show_y = FALSE) + +p2c <- .gain_panel(gain_df, + "Overlap\ncoefficient", + "Overlap coefficient (1 = identical)", + ref_val = 1, + show_y = FALSE) + + theme(legend.position = "bottom") + +fig2 <- (p2a | p2b | p2c) + + plot_layout(guides = "collect") & + theme(legend.position = "bottom") + +ggsave(file.path(OUTPUT_DIR, "fig_ablation_gain.pdf"), + fig2, width = FIG_WIDTH_WIDE * 0.7, height = fig1_height, + device = "pdf") +ggsave(file.path(OUTPUT_DIR, "fig_ablation_gain.png"), + fig2, width = FIG_WIDTH_WIDE * 0.7, height = fig1_height, + dpi = FIG_DPI) +message("Figure 2 saved.") + + +# ── 6. Figure 3: Posterior predictive density overlays ─────────────────────── +# +# Selects the TOP_N_DENSITY pathogens with the largest JS_CA (arm A vs +# federated) for PRIMARY_DIST. Calls compute_predictive_cdf() for each of the +# three arms and overlays the posterior-median density with a ±95% CrI ribbon. +# Falls back gracefully if a fit is absent. + +dist_code_map <- c( + "Log-normal" = "lognormal", "Gamma" = "gamma", "Weibull" = "weibull", + "Burr XII" = "burr", "Gen. gamma" = "gg" +) + +# Identify the top pathogens to plot. +top_pathogens <- comparison_tbl |> + filter(dist == PRIMARY_DIST, !is.na(JS_CA)) |> + arrange(desc(JS_CA)) |> + slice_head(n = TOP_N_DENSITY) |> + pull(pathogen) + +if (length(top_pathogens) == 0L) { + message("No JS_CA values available — skipping Figure 3.") +} else { + + cdf_dist_name <- dist_code_map[[PRIMARY_DIST]] + + # Determine a common x upper bound for each pathogen using the maximum + # pred_q95 upper CrI bound across all three arms. + .x_upper <- function(pathogen_name) { + row <- filter(comparison_tbl, + pathogen == pathogen_name & dist == PRIMARY_DIST) + vals <- sapply(c("A_pred_q95", "B_pred_q95", "C_pred_q95"), function(col) { + p <- tryCatch(.parse_cri(row[[col]])[["hi"]], error = function(e) NA_real_) + p + }) + v <- max(vals, na.rm = TRUE) + if (!is.finite(v)) 50 else v * 1.25 + } + + # Build density data for one arm of one pathogen. + .get_density_df <- function(fit, arm_label, pathogen_name) { + if (is.null(fit)) return(NULL) + x_upper <- .x_upper(pathogen_name) + x_seq <- seq(0.001, x_upper, length.out = 400L) + cdf_obj <- tryCatch( + compute_predictive_cdf(fit, cdf_dist_name, x_seq = x_seq, n_draws = 300L), + error = function(e) NULL + ) + if (is.null(cdf_obj)) return(NULL) + + dx <- diff(x_seq)[1L] + .to_pdf <- function(F_vals) { + f <- c(diff(pmax(pmin(F_vals, 1 - 1e-9), 1e-9)), 0) / dx + f[f < 0] <- 0 + f + } + + # Convert CDF summary columns to PDF by forward-differencing. + data.frame( + x = x_seq, + pdf_med = .to_pdf(cdf_obj$summary$median), + pdf_lo = .to_pdf(cdf_obj$summary$low), # CrI on density + pdf_hi = .to_pdf(cdf_obj$summary$high), + arm = arm_label, + pathogen = pathogen_name + ) + } + + density_list <- list() + + for (pg in top_pathogens) { + arm_names <- c("individual_only", "summary_only", "federated") + arm_labels <- c("A", "B", "C") + primary_key <- names(dist_code_map[dist_code_map == cdf_dist_name & + names(dist_code_map) == PRIMARY_DIST]) + # Map PRIMARY_DIST label back to DIST_CODES key + dist_key_map <- c( + "Log-normal" = "lognormal", "Gamma" = "gamma", "Weibull" = "weibull", + "Burr XII" = "burr", "Gen. gamma" = "gengamma" + ) + dk <- dist_key_map[[PRIMARY_DIST]] + + for (i in seq_along(arm_names)) { + res <- ablation_fits[[pg]][[arm_names[i]]][[dk]] + if (is.null(res) || isTRUE(res$skipped) || is.null(res$fit)) next + df <- .get_density_df(res$fit, arm_labels[i], pg) + if (!is.null(df)) density_list[[length(density_list) + 1L]] <- df + } + } + + if (length(density_list) == 0L) { + message("No density data could be computed — skipping Figure 3.") + } else { + + density_df <- bind_rows(density_list) |> + mutate( + arm = factor(arm, levels = c("A", "B", "C")), + pathogen = factor(pathogen, levels = top_pathogens) + ) + + # Also annotate each panel with the JS_CA and JS_CB values. + js_labels <- comparison_tbl |> + filter(pathogen %in% top_pathogens, dist == PRIMARY_DIST) |> + mutate( + label = sprintf( + "JS[CA] == %.3f~~~~~JS[CB] == %.3f", + round(JS_CA, 3L), round(JS_CB, 3L) + ), + pathogen = factor(pathogen, levels = top_pathogens) + ) |> + select(pathogen, label) + + # x position for label: 70th percentile of x range per pathogen + label_x <- density_df |> + group_by(pathogen) |> + summarise(x_pos = quantile(x, 0.72), .groups = "drop") + y_pos <- density_df |> + group_by(pathogen) |> + summarise(y_pos = max(pdf_hi, na.rm = TRUE) * 0.97, .groups = "drop") + js_labels <- left_join(js_labels, label_x, by = "pathogen") |> + left_join(y_pos, by = "pathogen") + + # Ribbon alpha by arm: federated slightly more prominent. + ribbon_alpha <- c("A" = 0.20, "B" = 0.20, "C" = 0.25) + line_lwd <- c("A" = 0.6, "B" = 0.6, "C" = 0.9) + + fig3 <- ggplot(density_df, + aes(x = x, group = arm, fill = arm, colour = arm)) + + geom_ribbon(aes(ymin = pdf_lo, ymax = pdf_hi), + alpha = 0.18, colour = NA) + + geom_line(aes(y = pdf_med, linewidth = arm)) + + geom_text(data = js_labels, + aes(x = x_pos, y = y_pos, label = label), + inherit.aes = FALSE, + size = 2.6, hjust = 0, parse = TRUE, colour = "grey30") + + scale_fill_manual(values = ARM_COLOURS, labels = ARM_LABELS) + + scale_colour_manual(values = ARM_COLOURS, labels = ARM_LABELS) + + scale_linewidth_manual(values = line_lwd, labels = ARM_LABELS) + + facet_wrap(~ pathogen, scales = "free", ncol = 2L) + + labs( + x = "Incubation period (days)", + y = "Density", + caption = paste0( + "Top ", TOP_N_DENSITY, " pathogens by JS divergence (arm A vs federated). ", + PRIMARY_DIST, " distribution. ", + "Shaded bands: 95% posterior credible intervals." + ) + ) + + theme_ablation() + + theme( + axis.title.y = element_text(), + legend.position = "bottom" + ) + + guides( + fill = guide_legend(override.aes = list(alpha = 0.4)), + linewidth = guide_legend(override.aes = list(linewidth = 1)) + ) + + n_rows_fig3 <- ceiling(length(top_pathogens) / 2L) + fig3_height <- max(4, n_rows_fig3 * 3 + 1.5) + + ggsave(file.path(OUTPUT_DIR, "fig_ablation_densities.pdf"), + fig3, width = FIG_WIDTH_HALF * 1.5, height = fig3_height, + device = "pdf") + ggsave(file.path(OUTPUT_DIR, "fig_ablation_densities.png"), + fig3, width = FIG_WIDTH_HALF * 1.5, height = fig3_height, + dpi = FIG_DPI) + message("Figure 3 saved.") + } +} + + +# ── 7. Supplementary: gain strip faceted over all distributions ─────────────── +# +# Same as Figure 2 but faceted by distribution so the reader can check whether +# the federated-gain pattern is consistent across distributional assumptions. + +gain_all_dists <- comparison_tbl |> + mutate(pathogen = factor(pathogen, levels = pathogen_order)) |> + select(pathogen, dist, + ratio_CA, ratio_CB, + JS_CA, JS_CB, + OVL_CA, OVL_CB) |> + pivot_longer( + cols = c(ratio_CA, ratio_CB, JS_CA, JS_CB, OVL_CA, OVL_CB), + names_to = c("metric", "comparison"), + names_sep = "_", + values_to = "value" + ) |> + filter(!is.na(value)) |> + mutate( + comparison = factor(comparison, + levels = c("CA", "CB"), + labels = c("vs A (individual-level)", + "vs B (summary-stats)")), + metric = factor(metric, + levels = c("ratio", "JS", "OVL"), + labels = c("Interval ratio", + "JS divergence (bits)", + "Overlap coefficient")) + ) + +ref_vals <- c("Interval ratio" = 1, "JS divergence (bits)" = 0, + "Overlap coefficient" = 1) + +fig_supp <- ggplot(gain_all_dists, + aes(x = value, y = pathogen, + colour = comparison, shape = comparison)) + + geom_point(size = 1.8, alpha = 0.8) + + geom_vline( + data = data.frame( + metric = names(ref_vals), + ref_val = unname(ref_vals) + ), + aes(xintercept = ref_val), + colour = "grey50", + linetype = "dashed", + linewidth = 0.35, + inherit.aes = FALSE + ) + + scale_colour_manual( + values = c("vs A (individual-level)" = ARM_COLOURS[["A"]], + "vs B (summary-stats)" = ARM_COLOURS[["B"]]) + ) + + scale_shape_manual( + values = c("vs A (individual-level)" = 19, + "vs B (summary-stats)" = 17) + ) + + facet_grid(dist ~ metric, scales = "free_x") + + labs(x = NULL) + + theme_ablation() + + theme( + legend.position = "bottom", + strip.text.y = element_text(angle = 0, hjust = 0) + ) + +ggsave(file.path(OUTPUT_DIR, "fig_ablation_gain_all_dists.pdf"), + fig_supp, + width = FIG_WIDTH_WIDE * 0.7, + height = max(5, n_pathogens * FIG_HEIGHT_ROW * 5 + 1.5), + device = "pdf") +ggsave(file.path(OUTPUT_DIR, "fig_ablation_gain_all_dists.png"), + fig_supp, + width = FIG_WIDTH_WIDE * 0.7, + height = max(5, n_pathogens * FIG_HEIGHT_ROW * 5 + 1.5), + dpi = FIG_DPI) +message("Supplementary gain figure (all distributions) saved.") + + +message("\nAll figures written to: ", OUTPUT_DIR) From 84bd0b81a6f84f5e8bf9b7f15b2af86299b79bc1 Mon Sep 17 00:00:00 2001 From: cm401 Date: Sun, 10 May 2026 10:37:58 +0800 Subject: [PATCH 03/18] Fix position_dodge and geom_vline arguments in ablation plot position_dodge() takes width= not height=; geom_vline() does not accept inherit.aes. Co-Authored-By: Claude Sonnet 4.6 --- analysis/plot_data_format_ablation.R | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) diff --git a/analysis/plot_data_format_ablation.R b/analysis/plot_data_format_ablation.R index 652fcd5..3586ef3 100644 --- a/analysis/plot_data_format_ablation.R +++ b/analysis/plot_data_format_ablation.R @@ -175,15 +175,15 @@ fig1 <- ggplot(forest_long, geom_vline(data = ref_lines, aes(xintercept = ref_med), colour = ARM_COLOURS[["C"]], linewidth = 0.25, - linetype = "dashed", inherit.aes = FALSE) + + linetype = "dashed") + # CrI bars geom_errorbarh(aes(xmin = lo, xmax = hi), height = 0, linewidth = 0.5, - position = position_dodge(height = dodge_height)) + + position = position_dodge(width = dodge_height)) + # Point estimates geom_point(size = 1.8, - position = position_dodge(height = dodge_height)) + + position = position_dodge(width = dodge_height)) + scale_x_log10( breaks = c(1, 2, 5, 10, 20, 50), labels = c("1", "2", "5", "10", "20", "50") From c6bab00f756ef9c9aace5a5e696e90a8d107734e Mon Sep 17 00:00:00 2001 From: cm401 Date: Sun, 10 May 2026 19:33:44 +0800 Subject: [PATCH 04/18] Add combined multi-panel figure (fig_ablation_combined) Three-column patchwork layout: (A) dumbbell forest plot with grey segment connecting arm A to arm B point estimates and prominent arm C overlay, (B) interval ratio strip, (C) JS divergence strip. Also fixes stray inherit.aes warning in supplementary figure geom_vline call. Co-Authored-By: Claude Sonnet 4.6 --- analysis/plot_data_format_ablation.R | 179 ++++++++++++++++++++++++++- 1 file changed, 178 insertions(+), 1 deletion(-) diff --git a/analysis/plot_data_format_ablation.R b/analysis/plot_data_format_ablation.R index 3586ef3..720ab11 100644 --- a/analysis/plot_data_format_ablation.R +++ b/analysis/plot_data_format_ablation.R @@ -520,7 +520,6 @@ fig_supp <- ggplot(gain_all_dists, colour = "grey50", linetype = "dashed", linewidth = 0.35, - inherit.aes = FALSE ) + scale_colour_manual( values = c("vs A (individual-level)" = ARM_COLOURS[["A"]], @@ -551,4 +550,182 @@ ggsave(file.path(OUTPUT_DIR, "fig_ablation_gain_all_dists.png"), message("Supplementary gain figure (all distributions) saved.") +# ── 8. Figure 4 (combined): Multi-panel federated gain summary ─────────────── +# +# A single three-column figure suitable for the main paper body. +# +# Panel A (wide) — Dumbbell forest plot (PRIMARY_DIST only) +# A grey segment connects the arm A and arm B point estimates; its length +# encodes how much the two single-format arms disagree. Three +# point + CrI layers (A, B, C) are drawn on top. Arm C is rendered more +# prominently (larger point, heavier error bar) because it is the focus. +# A log-scaled x-axis accommodates the range across pathogens. +# +# Panel B (narrow) — Interval ratio (arm_width / C_width) +# One dot per arm × pathogen pair (CA in blue, CB in orange). +# Reference line at 1: points to the right mean the single-format arm is +# wider (less precise) than the federated model. +# +# Panel C (narrow) — Jensen-Shannon divergence (bits) +# JS_CA and JS_CB on the same scale. Larger values signal a bigger +# distributional shift between that arm and the federated result. +# Reference line at 0 (identical distributions). +# +# Panels A/B/C share the pathogen y-axis via patchwork alignment. +# Panel labels A/B/C are added automatically by plot_annotation(). +# The combined legend is collected at the bottom. +# +# Output: results/figures/fig_ablation_combined.{pdf,png} + +# ── 8a. Wide data for PRIMARY_DIST ─────────────────────────────────────────── + +wide_ln <- comparison_tbl |> + filter(dist == PRIMARY_DIST) |> + mutate(pathogen = factor(pathogen, levels = pathogen_order)) |> + rowwise() |> + mutate( + A_med = .parse_cri(A_pred_median)[["med"]], + A_lo = .parse_cri(A_pred_median)[["lo"]], + A_hi = .parse_cri(A_pred_median)[["hi"]], + B_med = .parse_cri(B_pred_median)[["med"]], + B_lo = .parse_cri(B_pred_median)[["lo"]], + B_hi = .parse_cri(B_pred_median)[["hi"]], + C_med = .parse_cri(C_pred_median)[["med"]], + C_lo = .parse_cri(C_pred_median)[["lo"]], + C_hi = .parse_cri(C_pred_median)[["hi"]] + ) |> + ungroup() + +# Long format for the point + CrI layers. +forest_ln <- wide_ln |> + select(pathogen, A_med, A_lo, A_hi, B_med, B_lo, B_hi, C_med, C_lo, C_hi) |> + pivot_longer( + cols = -pathogen, + names_to = c("arm", ".value"), + names_pattern = "^(.)_(.*)" + ) |> + filter(!is.na(med)) |> + mutate(arm = factor(arm, levels = c("A", "B", "C"))) + +# Segments connecting A to B point estimates (length = arm disagreement). +dumbbell_segs <- wide_ln |> + filter(!is.na(A_med), !is.na(B_med)) + +# Arm C rendered more prominently than A/B. +ARM_SIZES_C <- c("A" = 1.8, "B" = 1.8, "C" = 2.8) +ARM_EBW_C <- c("A" = 0.35, "B" = 0.35, "C" = 0.75) + +# ── 8b. Panel A: dumbbell forest plot ──────────────────────────────────────── + +pA <- ggplot(forest_ln, + aes(x = med, y = pathogen, colour = arm, shape = arm)) + + geom_segment( + data = dumbbell_segs, + aes(x = A_med, xend = B_med, y = pathogen, yend = pathogen), + colour = "grey78", + linewidth = 0.9, + lineend = "round" + ) + + geom_errorbarh( + aes(xmin = lo, xmax = hi, linewidth = arm), + height = 0 + ) + + geom_point(aes(size = arm)) + + scale_x_log10( + breaks = c(1, 2, 5, 10, 20, 50), + labels = c("1", "2", "5", "10", "20", "50") + ) + + scale_colour_manual(values = ARM_COLOURS, labels = ARM_LABELS) + + scale_shape_manual( values = ARM_SHAPES, labels = ARM_LABELS) + + scale_size_manual( values = ARM_SIZES_C, labels = ARM_LABELS) + + scale_linewidth_manual(values = ARM_EBW_C, labels = ARM_LABELS) + + labs(x = "Posterior predictive median (days, log scale)", + subtitle = PRIMARY_DIST) + + theme_ablation() + + guides( + colour = guide_legend(override.aes = list(size = 2.5)), + size = "none", + linewidth = "none", + shape = "none" + ) + +# ── 8c. Panels B and C: gain metrics ───────────────────────────────────────── + +# Shared data for the two metric panels. +COMP_COLOURS <- c( + "vs A (individual-level)" = ARM_COLOURS[["A"]], + "vs B (summary-stats)" = ARM_COLOURS[["B"]] +) +COMP_SHAPES <- c("vs A (individual-level)" = 19, "vs B (summary-stats)" = 17) + +gain_ln <- comparison_tbl |> + filter(dist == PRIMARY_DIST) |> + mutate(pathogen = factor(pathogen, levels = pathogen_order)) |> + select(pathogen, ratio_CA, ratio_CB, JS_CA, JS_CB) |> + pivot_longer( + cols = -pathogen, + names_to = c("metric", "comparison"), + names_pattern = "^(.*)_(C[AB])$" + ) |> + filter(!is.na(value)) |> + mutate( + comparison = factor( + comparison, + levels = c("CA", "CB"), + labels = c("vs A (individual-level)", "vs B (summary-stats)") + ) + ) + +# Shared y scale for panels B and C — same factor, labels suppressed. +y_shared <- list( + scale_y_discrete(drop = FALSE), + theme_ablation(), + theme( + axis.text.y = element_blank(), + axis.ticks.y = element_blank(), + axis.title.y = element_blank() + ) +) + +pB <- ggplot( + filter(gain_ln, metric == "ratio"), + aes(x = value, y = pathogen, colour = comparison, shape = comparison) + ) + + geom_vline(xintercept = 1, linetype = "dashed", + colour = "grey50", linewidth = 0.4) + + geom_point(size = 2.2) + + scale_colour_manual(values = COMP_COLOURS) + + scale_shape_manual( values = COMP_SHAPES) + + labs(x = "Interval ratio\n(arm width / C width)") + + y_shared + +pC <- ggplot( + filter(gain_ln, metric == "JS"), + aes(x = value, y = pathogen, colour = comparison, shape = comparison) + ) + + geom_vline(xintercept = 0, linetype = "dashed", + colour = "grey50", linewidth = 0.4) + + geom_point(size = 2.2) + + scale_colour_manual(values = COMP_COLOURS) + + scale_shape_manual( values = COMP_SHAPES) + + labs(x = "JS divergence (bits)") + + y_shared + +# ── 8d. Assemble and save ───────────────────────────────────────────────────── + +fig4 <- (pA | pB | pC) + + plot_layout(widths = c(3, 1, 1), guides = "collect") + + plot_annotation(tag_levels = "A") & + theme(legend.position = "bottom", + plot.tag = element_text(face = "bold", size = 10)) + +fig4_height <- max(4, n_pathogens * FIG_HEIGHT_ROW + 1.8) + +ggsave(file.path(OUTPUT_DIR, "fig_ablation_combined.pdf"), + fig4, width = FIG_WIDTH_WIDE, height = fig4_height, device = "pdf") +ggsave(file.path(OUTPUT_DIR, "fig_ablation_combined.png"), + fig4, width = FIG_WIDTH_WIDE, height = fig4_height, dpi = FIG_DPI) +message("Figure 4 (combined) saved.") + + message("\nAll figures written to: ", OUTPUT_DIR) From d6c7c92e8fbdf8a0be2099beb7a420731c32f0a5 Mon Sep 17 00:00:00 2001 From: cm401 Date: Sun, 10 May 2026 19:43:33 +0800 Subject: [PATCH 05/18] Use best-fitting distribution per pathogen in combined figure; add position_dodge to Panel A Panel A of the combined figure (fig_ablation_combined) now selects the distribution with highest LOO-ELPD for arm C per pathogen rather than always using the primary (log-normal) distribution. Tie-breaking follows lognormal > gamma > Weibull > Burr XII > gen. gamma; pathogens with no valid ELPD fall back to PRIMARY_DIST. Y-axis labels show the selected distribution in parentheses. Arms A/B/C are now separated with position_dodge(width = 0.5) on both the errorbar and point layers so estimates do not overlap. Panels B and C use the same per-pathogen best-distribution data (gain_best). Co-Authored-By: Claude Sonnet 4.6 --- analysis/plot_data_format_ablation.R | 97 +++++++++++++++++++++------- 1 file changed, 72 insertions(+), 25 deletions(-) diff --git a/analysis/plot_data_format_ablation.R b/analysis/plot_data_format_ablation.R index 720ab11..bb9c44b 100644 --- a/analysis/plot_data_format_ablation.R +++ b/analysis/plot_data_format_ablation.R @@ -577,10 +577,42 @@ message("Supplementary gain figure (all distributions) saved.") # # Output: results/figures/fig_ablation_combined.{pdf,png} -# ── 8a. Wide data for PRIMARY_DIST ─────────────────────────────────────────── +# ── 8a. Best-fitting distribution per pathogen ─────────────────────────────── +# +# "Best" = highest study-level PSIS-LOO ELPD for arm C (federated). +# Tie-breaking preference: lognormal > gamma > weibull > burr > gengamma. +# Falls back to PRIMARY_DIST for pathogens where all ELPD values are NA +# (fewer than 4 datasets in arm C). + +DIST_PREFERENCE <- c("Log-normal" = 1L, "Gamma" = 2L, "Weibull" = 3L, + "Burr XII" = 4L, "Gen. gamma" = 5L) + +best_dist_tbl <- comparison_tbl |> + select(pathogen, dist, C_elpd) |> + filter(!is.na(C_elpd)) |> + mutate(pref = DIST_PREFERENCE[dist]) |> + group_by(pathogen) |> + arrange(desc(C_elpd), pref, .by_group = TRUE) |> + slice_head(n = 1L) |> + ungroup() |> + select(pathogen, best_dist = dist) + +# Pathogens with no valid ELPD fall back to PRIMARY_DIST. +fallback_pathogens <- setdiff(unique(comparison_tbl$pathogen), + best_dist_tbl$pathogen) +if (length(fallback_pathogens) > 0L) { + best_dist_tbl <- bind_rows( + best_dist_tbl, + data.frame(pathogen = fallback_pathogens, + best_dist = PRIMARY_DIST, + stringsAsFactors = FALSE) + ) +} -wide_ln <- comparison_tbl |> - filter(dist == PRIMARY_DIST) |> +# ── 8b. Wide data using best distribution per pathogen ─────────────────────── + +wide_best <- comparison_tbl |> + inner_join(best_dist_tbl, by = c("pathogen", "dist" = "best_dist")) |> mutate(pathogen = factor(pathogen, levels = pathogen_order)) |> rowwise() |> mutate( @@ -596,29 +628,39 @@ wide_ln <- comparison_tbl |> ) |> ungroup() -# Long format for the point + CrI layers. -forest_ln <- wide_ln |> +# y-axis labels: "Pathogen (Distribution)". +dist_label_lookup <- setNames( + paste0(wide_best$pathogen, "\n(", wide_best$dist, ")"), + as.character(wide_best$pathogen) +) + +# Long format for point + CrI layers. +forest_best <- wide_best |> select(pathogen, A_med, A_lo, A_hi, B_med, B_lo, B_hi, C_med, C_lo, C_hi) |> pivot_longer( - cols = -pathogen, - names_to = c("arm", ".value"), + cols = -pathogen, + names_to = c("arm", ".value"), names_pattern = "^(.)_(.*)" ) |> filter(!is.na(med)) |> mutate(arm = factor(arm, levels = c("A", "B", "C"))) -# Segments connecting A to B point estimates (length = arm disagreement). -dumbbell_segs <- wide_ln |> +# Segments connecting A to B (length = arm disagreement). +dumbbell_segs <- wide_best |> filter(!is.na(A_med), !is.na(B_med)) +# Dodging width for the three arms along the discrete y-axis. +COMB_DODGE <- 0.5 + # Arm C rendered more prominently than A/B. ARM_SIZES_C <- c("A" = 1.8, "B" = 1.8, "C" = 2.8) ARM_EBW_C <- c("A" = 0.35, "B" = 0.35, "C" = 0.75) -# ── 8b. Panel A: dumbbell forest plot ──────────────────────────────────────── +# ── 8c. Panel A: dumbbell forest plot ──────────────────────────────────────── -pA <- ggplot(forest_ln, +pA <- ggplot(forest_best, aes(x = med, y = pathogen, colour = arm, shape = arm)) + + # Grey segment from A to B (drawn at the un-dodged y position). geom_segment( data = dumbbell_segs, aes(x = A_med, xend = B_med, y = pathogen, yend = pathogen), @@ -626,21 +668,28 @@ pA <- ggplot(forest_ln, linewidth = 0.9, lineend = "round" ) + + # CrI bars, dodged so the three arms sit at distinct y positions. geom_errorbarh( aes(xmin = lo, xmax = hi, linewidth = arm), - height = 0 + height = 0, + position = position_dodge(width = COMB_DODGE) + ) + + # Point estimates, same dodge. + geom_point( + aes(size = arm), + position = position_dodge(width = COMB_DODGE) ) + - geom_point(aes(size = arm)) + scale_x_log10( breaks = c(1, 2, 5, 10, 20, 50), labels = c("1", "2", "5", "10", "20", "50") ) + - scale_colour_manual(values = ARM_COLOURS, labels = ARM_LABELS) + - scale_shape_manual( values = ARM_SHAPES, labels = ARM_LABELS) + - scale_size_manual( values = ARM_SIZES_C, labels = ARM_LABELS) + - scale_linewidth_manual(values = ARM_EBW_C, labels = ARM_LABELS) + + scale_y_discrete(labels = dist_label_lookup) + + scale_colour_manual(values = ARM_COLOURS, labels = ARM_LABELS) + + scale_shape_manual( values = ARM_SHAPES, labels = ARM_LABELS) + + scale_size_manual( values = ARM_SIZES_C, labels = ARM_LABELS) + + scale_linewidth_manual(values = ARM_EBW_C, labels = ARM_LABELS) + labs(x = "Posterior predictive median (days, log scale)", - subtitle = PRIMARY_DIST) + + subtitle = "Best-fitting distribution per pathogen (highest LOO-ELPD)") + theme_ablation() + guides( colour = guide_legend(override.aes = list(size = 2.5)), @@ -649,22 +698,20 @@ pA <- ggplot(forest_ln, shape = "none" ) -# ── 8c. Panels B and C: gain metrics ───────────────────────────────────────── +# ── 8d. Panels B and C: gain metrics ───────────────────────────────────────── -# Shared data for the two metric panels. +# Shared colours/shapes for the two comparison directions. COMP_COLOURS <- c( "vs A (individual-level)" = ARM_COLOURS[["A"]], "vs B (summary-stats)" = ARM_COLOURS[["B"]] ) COMP_SHAPES <- c("vs A (individual-level)" = 19, "vs B (summary-stats)" = 17) -gain_ln <- comparison_tbl |> - filter(dist == PRIMARY_DIST) |> - mutate(pathogen = factor(pathogen, levels = pathogen_order)) |> +gain_best <- wide_best |> select(pathogen, ratio_CA, ratio_CB, JS_CA, JS_CB) |> pivot_longer( - cols = -pathogen, - names_to = c("metric", "comparison"), + cols = -pathogen, + names_to = c("metric", "comparison"), names_pattern = "^(.*)_(C[AB])$" ) |> filter(!is.na(value)) |> From ece11079b771d8ad0ab3fb7b77ba5caf0e5c9ab9 Mon Sep 17 00:00:00 2001 From: cm401 Date: Sun, 10 May 2026 19:49:20 +0800 Subject: [PATCH 06/18] Use main-analysis model weights to select best distribution in combined figure Section 8a now loads results/model_weights.rds (pre-computed by compute_pathogen_model_bayes_factors() from main_results.rds) and picks the highest pseudo-Bayes-factor distribution per pathogen, matching the criterion used in plot_main_figure() and make_supplementary_figures.R. Previously, best-distribution was determined from the ablation arm C ELPD, which was inconsistent with the rest of the paper. Internal Stan names (lognormal, gamma, weibull, burr, gengamma) are mapped to display labels via MAIN_DIST_TO_DISPLAY before joining to comparison_tbl. Falls back to PRIMARY_DIST if a pathogen is absent from model_weights or its best distribution was not fitted in the ablation. Also fixes pB and pC panels which were referencing gain_ln instead of gain_best. Co-Authored-By: Claude Sonnet 4.6 --- analysis/plot_data_format_ablation.R | 86 +++++++++++++++++----------- 1 file changed, 54 insertions(+), 32 deletions(-) diff --git a/analysis/plot_data_format_ablation.R b/analysis/plot_data_format_ablation.R index bb9c44b..60aae20 100644 --- a/analysis/plot_data_format_ablation.R +++ b/analysis/plot_data_format_ablation.R @@ -577,36 +577,58 @@ message("Supplementary gain figure (all distributions) saved.") # # Output: results/figures/fig_ablation_combined.{pdf,png} -# ── 8a. Best-fitting distribution per pathogen ─────────────────────────────── +# ── 8a. Best-fitting distribution per pathogen (from main analysis) ────────── # -# "Best" = highest study-level PSIS-LOO ELPD for arm C (federated). -# Tie-breaking preference: lognormal > gamma > weibull > burr > gengamma. -# Falls back to PRIMARY_DIST for pathogens where all ELPD values are NA -# (fewer than 4 datasets in arm C). - -DIST_PREFERENCE <- c("Log-normal" = 1L, "Gamma" = 2L, "Weibull" = 3L, - "Burr XII" = 4L, "Gen. gamma" = 5L) - -best_dist_tbl <- comparison_tbl |> - select(pathogen, dist, C_elpd) |> - filter(!is.na(C_elpd)) |> - mutate(pref = DIST_PREFERENCE[dist]) |> - group_by(pathogen) |> - arrange(desc(C_elpd), pref, .by_group = TRUE) |> - slice_head(n = 1L) |> - ungroup() |> - select(pathogen, best_dist = dist) - -# Pathogens with no valid ELPD fall back to PRIMARY_DIST. -fallback_pathogens <- setdiff(unique(comparison_tbl$pathogen), - best_dist_tbl$pathogen) -if (length(fallback_pathogens) > 0L) { - best_dist_tbl <- bind_rows( - best_dist_tbl, - data.frame(pathogen = fallback_pathogens, - best_dist = PRIMARY_DIST, - stringsAsFactors = FALSE) - ) +# Load model weights pre-computed from main_results.rds via +# compute_pathogen_model_bayes_factors(). The best distribution is the +# first entry in each pathogen's weight vector (highest pseudo-Bayes factor +# weight). This is the same criterion used in plot_main_figure() and +# make_supplementary_figures.R, ensuring consistency across the paper. +# +# Falls back to PRIMARY_DIST for pathogens absent from model_weights or +# whose best distribution was not fitted in the ablation. + +# Map internal Stan dist names -> display labels used in comparison_tbl$dist +MAIN_DIST_TO_DISPLAY <- c( + lognormal = "Log-normal", + gamma = "Gamma", + weibull = "Weibull", + burr = "Burr XII", + gengamma = "Gen. gamma" +) + +model_weights_file <- here::here("results", "model_weights.rds") +if (!file.exists(model_weights_file)) + stop("model_weights.rds not found at: ", model_weights_file, + "\nPre-compute via: mw <- compute_pathogen_model_bayes_factors(all_results); ", + "saveRDS(mw, 'results/model_weights.rds')") + +main_model_weights <- readRDS(model_weights_file) + +# For each ablation pathogen, extract the highest-weight distribution from +# the main analysis. +ablation_pathogens <- unique(comparison_tbl$pathogen) +best_internal <- vapply(ablation_pathogens, function(p) { + w <- main_model_weights[[p]] + if (is.null(w) || length(w) == 0L) return(NA_character_) + names(w)[1L] # best = highest pseudo-BF weight +}, character(1L)) + +best_dist_tbl <- data.frame( + pathogen = ablation_pathogens, + best_dist = MAIN_DIST_TO_DISPLAY[best_internal], + stringsAsFactors = FALSE +) + +# Fall back to PRIMARY_DIST where the main-analysis best dist is unavailable +# in the ablation results (e.g. Gen. gamma was not fitted in arm B/C). +available_dists <- unique(comparison_tbl$dist) +needs_fallback <- is.na(best_dist_tbl$best_dist) | + !(best_dist_tbl$best_dist %in% available_dists) +if (any(needs_fallback)) { + message("Note: falling back to PRIMARY_DIST for: ", + paste(best_dist_tbl$pathogen[needs_fallback], collapse = ", ")) + best_dist_tbl$best_dist[needs_fallback] <- PRIMARY_DIST } # ── 8b. Wide data using best distribution per pathogen ─────────────────────── @@ -689,7 +711,7 @@ pA <- ggplot(forest_best, scale_size_manual( values = ARM_SIZES_C, labels = ARM_LABELS) + scale_linewidth_manual(values = ARM_EBW_C, labels = ARM_LABELS) + labs(x = "Posterior predictive median (days, log scale)", - subtitle = "Best-fitting distribution per pathogen (highest LOO-ELPD)") + + subtitle = "Best-fitting distribution per pathogen (main analysis model weights)") + theme_ablation() + guides( colour = guide_legend(override.aes = list(size = 2.5)), @@ -735,7 +757,7 @@ y_shared <- list( ) pB <- ggplot( - filter(gain_ln, metric == "ratio"), + filter(gain_best, metric == "ratio"), aes(x = value, y = pathogen, colour = comparison, shape = comparison) ) + geom_vline(xintercept = 1, linetype = "dashed", @@ -747,7 +769,7 @@ pB <- ggplot( y_shared pC <- ggplot( - filter(gain_ln, metric == "JS"), + filter(gain_best, metric == "JS"), aes(x = value, y = pathogen, colour = comparison, shape = comparison) ) + geom_vline(xintercept = 0, linetype = "dashed", From d2a34d13c135cf22a757737a0f5851ef081d4ea9 Mon Sep 17 00:00:00 2001 From: cm401 Date: Sun, 10 May 2026 20:10:26 +0800 Subject: [PATCH 07/18] Add P95 facet to combined figure; unify arm labels; widen forest panel Combined figure (fig_ablation_combined) changes: - Panel A is now a faceted dumbbell forest with two side-by-side columns: 'Median (P50)' and '95th percentile (P95)', sharing the y-axis with independent x-scales (log10). Both facets use the same position_dodge, colour/shape encoding and dumbbell backbone. - patchwork widths changed from c(3,1,1) to c(2,1,1). - ARM_LABELS now reads 'Individual-level only', 'Summary-statistics only', 'Federated' (letter suffixes removed). COMP_COLOURS/COMP_SHAPES in section 8 use the same strings so the collected legend is unified: blue/orange/green mean the same entity in all three panels with no duplicate entries. - P95 CrI bounds parsed from A_pred_q95/B_pred_q95/C_pred_q95 columns in the same rowwise() block as the median. Co-Authored-By: Claude Sonnet 4.6 --- analysis/plot_data_format_ablation.R | 128 ++++++++++++++++++--------- 1 file changed, 87 insertions(+), 41 deletions(-) diff --git a/analysis/plot_data_format_ablation.R b/analysis/plot_data_format_ablation.R index 60aae20..ada4d46 100644 --- a/analysis/plot_data_format_ablation.R +++ b/analysis/plot_data_format_ablation.R @@ -63,9 +63,9 @@ ARM_COLOURS <- c( "C" = "#009E73" # green ) ARM_LABELS <- c( - "A" = "Individual-level only (A)", - "B" = "Summary-statistics only (B)", - "C" = "Federated (C)" + "A" = "Individual-level only", + "B" = "Summary-statistics only", + "C" = "Federated" ) ARM_SHAPES <- c("A" = 19, "B" = 17, "C" = 18) # circle, triangle, diamond @@ -638,15 +638,26 @@ wide_best <- comparison_tbl |> mutate(pathogen = factor(pathogen, levels = pathogen_order)) |> rowwise() |> mutate( - A_med = .parse_cri(A_pred_median)[["med"]], - A_lo = .parse_cri(A_pred_median)[["lo"]], - A_hi = .parse_cri(A_pred_median)[["hi"]], - B_med = .parse_cri(B_pred_median)[["med"]], - B_lo = .parse_cri(B_pred_median)[["lo"]], - B_hi = .parse_cri(B_pred_median)[["hi"]], - C_med = .parse_cri(C_pred_median)[["med"]], - C_lo = .parse_cri(C_pred_median)[["lo"]], - C_hi = .parse_cri(C_pred_median)[["hi"]] + # Median (P50) posterior predictive estimates + 95 % CrI bounds + A_med = .parse_cri(A_pred_median)[["med"]], + A_lo = .parse_cri(A_pred_median)[["lo"]], + A_hi = .parse_cri(A_pred_median)[["hi"]], + B_med = .parse_cri(B_pred_median)[["med"]], + B_lo = .parse_cri(B_pred_median)[["lo"]], + B_hi = .parse_cri(B_pred_median)[["hi"]], + C_med = .parse_cri(C_pred_median)[["med"]], + C_lo = .parse_cri(C_pred_median)[["lo"]], + C_hi = .parse_cri(C_pred_median)[["hi"]], + # 95th percentile (P95) posterior predictive estimates + 95 % CrI bounds + A_q95 = .parse_cri(A_pred_q95)[["med"]], + A_q95_lo = .parse_cri(A_pred_q95)[["lo"]], + A_q95_hi = .parse_cri(A_pred_q95)[["hi"]], + B_q95 = .parse_cri(B_pred_q95)[["med"]], + B_q95_lo = .parse_cri(B_pred_q95)[["lo"]], + B_q95_hi = .parse_cri(B_pred_q95)[["hi"]], + C_q95 = .parse_cri(C_pred_q95)[["med"]], + C_q95_lo = .parse_cri(C_pred_q95)[["lo"]], + C_q95_hi = .parse_cri(C_pred_q95)[["hi"]] ) |> ungroup() @@ -656,41 +667,74 @@ dist_label_lookup <- setNames( as.character(wide_best$pathogen) ) -# Long format for point + CrI layers. -forest_best <- wide_best |> +# Long format for point + CrI layers — both P50 and P95 as facets. +.to_forest_long <- function(wb, med_col, lo_col, hi_col, metric_label) { + wb |> + select(pathogen, + A_med = {{ med_col }}, A_lo = {{ lo_col }}, A_hi = {{ hi_col }}, + B_med = B_med, B_lo = B_lo, B_hi = B_hi, + C_med = C_med, C_lo = C_lo, C_hi = C_hi) |> + pivot_longer( + cols = -pathogen, + names_to = c("arm", ".value"), + names_pattern = "^(.)_(.*)" + ) |> + filter(!is.na(med)) |> + mutate(arm = factor(arm, levels = c("A", "B", "C")), + metric = metric_label) +} + +# Manually build each half so column references are unambiguous. +forest_p50 <- wide_best |> select(pathogen, A_med, A_lo, A_hi, B_med, B_lo, B_hi, C_med, C_lo, C_hi) |> - pivot_longer( - cols = -pathogen, - names_to = c("arm", ".value"), - names_pattern = "^(.)_(.*)" - ) |> + pivot_longer(-pathogen, names_to = c("arm", ".value"), + names_pattern = "^(.)_(.*)") |> filter(!is.na(med)) |> - mutate(arm = factor(arm, levels = c("A", "B", "C"))) + mutate(arm = factor(arm, levels = c("A", "B", "C")), metric = "Median (P50)") -# Segments connecting A to B (length = arm disagreement). -dumbbell_segs <- wide_best |> - filter(!is.na(A_med), !is.na(B_med)) +forest_p95 <- wide_best |> + select(pathogen, + A_med = A_q95, A_lo = A_q95_lo, A_hi = A_q95_hi, + B_med = B_q95, B_lo = B_q95_lo, B_hi = B_q95_hi, + C_med = C_q95, C_lo = C_q95_lo, C_hi = C_q95_hi) |> + pivot_longer(-pathogen, names_to = c("arm", ".value"), + names_pattern = "^(.)_(.*)") |> + filter(!is.na(med)) |> + mutate(arm = factor(arm, levels = c("A", "B", "C")), metric = "95th percentile (P95)") + +forest_both <- bind_rows(forest_p50, forest_p95) |> + mutate(metric = factor(metric, levels = c("Median (P50)", "95th percentile (P95)"))) + +# Dumbbell backbone segments for both metrics. +dumbbell_segs <- bind_rows( + wide_best |> + filter(!is.na(A_med), !is.na(B_med)) |> + transmute(pathogen, x = A_med, xend = B_med, metric = "Median (P50)"), + wide_best |> + filter(!is.na(A_q95), !is.na(B_q95)) |> + transmute(pathogen, x = A_q95, xend = B_q95, metric = "95th percentile (P95)") +) |> + mutate(metric = factor(metric, levels = c("Median (P50)", "95th percentile (P95)"))) -# Dodging width for the three arms along the discrete y-axis. +# Dodging width and size scales. COMB_DODGE <- 0.5 - -# Arm C rendered more prominently than A/B. ARM_SIZES_C <- c("A" = 1.8, "B" = 1.8, "C" = 2.8) ARM_EBW_C <- c("A" = 0.35, "B" = 0.35, "C" = 0.75) -# ── 8c. Panel A: dumbbell forest plot ──────────────────────────────────────── +# ── 8c. Panel A: faceted dumbbell forest (P50 | P95) ───────────────────────── -pA <- ggplot(forest_best, +pA <- ggplot(forest_both, aes(x = med, y = pathogen, colour = arm, shape = arm)) + - # Grey segment from A to B (drawn at the un-dodged y position). + # Grey dumbbell backbone (A-to-B range), drawn at the un-dodged y position. geom_segment( data = dumbbell_segs, - aes(x = A_med, xend = B_med, y = pathogen, yend = pathogen), + aes(x = x, xend = xend, y = pathogen, yend = pathogen), colour = "grey78", linewidth = 0.9, - lineend = "round" + lineend = "round", + inherit.aes = FALSE ) + - # CrI bars, dodged so the three arms sit at distinct y positions. + # 95 % posterior CrI bars, dodged vertically. geom_errorbarh( aes(xmin = lo, xmax = hi, linewidth = arm), height = 0, @@ -701,16 +745,17 @@ pA <- ggplot(forest_best, aes(size = arm), position = position_dodge(width = COMB_DODGE) ) + + facet_wrap(~ metric, nrow = 1L, scales = "free_x") + scale_x_log10( - breaks = c(1, 2, 5, 10, 20, 50), - labels = c("1", "2", "5", "10", "20", "50") + breaks = c(1, 2, 5, 10, 20, 50, 100), + labels = c("1", "2", "5", "10", "20", "50", "100") ) + scale_y_discrete(labels = dist_label_lookup) + scale_colour_manual(values = ARM_COLOURS, labels = ARM_LABELS) + scale_shape_manual( values = ARM_SHAPES, labels = ARM_LABELS) + scale_size_manual( values = ARM_SIZES_C, labels = ARM_LABELS) + scale_linewidth_manual(values = ARM_EBW_C, labels = ARM_LABELS) + - labs(x = "Posterior predictive median (days, log scale)", + labs(x = "Days (log scale)", subtitle = "Best-fitting distribution per pathogen (main analysis model weights)") + theme_ablation() + guides( @@ -722,12 +767,13 @@ pA <- ggplot(forest_best, # ── 8d. Panels B and C: gain metrics ───────────────────────────────────────── -# Shared colours/shapes for the two comparison directions. +# Colours/shapes for the two comparison directions, using the same names as +# ARM_LABELS so the collected legend is unified across all three panels. COMP_COLOURS <- c( - "vs A (individual-level)" = ARM_COLOURS[["A"]], - "vs B (summary-stats)" = ARM_COLOURS[["B"]] + "Individual-level only" = ARM_COLOURS[["A"]], + "Summary-statistics only" = ARM_COLOURS[["B"]] ) -COMP_SHAPES <- c("vs A (individual-level)" = 19, "vs B (summary-stats)" = 17) +COMP_SHAPES <- c("Individual-level only" = 19, "Summary-statistics only" = 17) gain_best <- wide_best |> select(pathogen, ratio_CA, ratio_CB, JS_CA, JS_CB) |> @@ -741,7 +787,7 @@ gain_best <- wide_best |> comparison = factor( comparison, levels = c("CA", "CB"), - labels = c("vs A (individual-level)", "vs B (summary-stats)") + labels = c("Individual-level only", "Summary-statistics only") ) ) @@ -783,7 +829,7 @@ pC <- ggplot( # ── 8d. Assemble and save ───────────────────────────────────────────────────── fig4 <- (pA | pB | pC) + - plot_layout(widths = c(3, 1, 1), guides = "collect") + + plot_layout(widths = c(2, 1, 1), guides = "collect") + plot_annotation(tag_levels = "A") & theme(legend.position = "bottom", plot.tag = element_text(face = "bold", size = 10)) From 34a61d5dd695d836e612f0a5d67e53f94bd70071 Mon Sep 17 00:00:00 2001 From: cm401 Date: Sun, 10 May 2026 21:00:41 +0800 Subject: [PATCH 08/18] Add mu0 and tau panels to decompose opposing forces in combined figure MIME-Version: 1.0 Content-Type: text/plain; charset=UTF-8 Content-Transfer-Encoding: 8bit The combined figure (fig_ablation_combined) now has four panels (2:1:1:1): Panel A — faceted dumbbell forest (Median P50 | P95): unchanged. Panel B — mu0 CrI width ratio (arm / C). mu0 is the population-level location parameter extracted via rstan::extract() from the stored Stan fits in ablation_fits. Its 95% CrI width measures pure estimation uncertainty, independent of tau. Values > 1 mean arm C has tighter knowledge of the population mean -- the information-gain force in isolation. Panel C — tau per arm with 95% CrI (three dodged points per pathogen). tau is parsed from the A_tau/B_tau/C_tau strings already in comparison_tbl. Shows whether arm C detects more between-study heterogeneity (the opposing force) or estimates it more precisely (narrower CrI despite same median). Panel D — predictive CrI ratio (arm / C): the combined confounded signal shown previously as Panel B. Now clearly labelled as the sum of both forces. All four panels share a single collected legend via unified ARM_COLOURS_FULL / ARM_SHAPES_FULL scales keyed by full label text so patchwork merges them correctly. Figure width scaled up by 1.35x to accommodate the extra panel. Co-Authored-By: Claude Sonnet 4.6 --- analysis/plot_data_format_ablation.R | 299 +++++++++++++++++---------- 1 file changed, 189 insertions(+), 110 deletions(-) diff --git a/analysis/plot_data_format_ablation.R b/analysis/plot_data_format_ablation.R index ada4d46..c472a5b 100644 --- a/analysis/plot_data_format_ablation.R +++ b/analysis/plot_data_format_ablation.R @@ -667,24 +667,8 @@ dist_label_lookup <- setNames( as.character(wide_best$pathogen) ) -# Long format for point + CrI layers — both P50 and P95 as facets. -.to_forest_long <- function(wb, med_col, lo_col, hi_col, metric_label) { - wb |> - select(pathogen, - A_med = {{ med_col }}, A_lo = {{ lo_col }}, A_hi = {{ hi_col }}, - B_med = B_med, B_lo = B_lo, B_hi = B_hi, - C_med = C_med, C_lo = C_lo, C_hi = C_hi) |> - pivot_longer( - cols = -pathogen, - names_to = c("arm", ".value"), - names_pattern = "^(.)_(.*)" - ) |> - filter(!is.na(med)) |> - mutate(arm = factor(arm, levels = c("A", "B", "C")), - metric = metric_label) -} +# ── 8c. Forest data — P50 and P95 facets ───────────────────────────────────── -# Manually build each half so column references are unambiguous. forest_p50 <- wide_best |> select(pathogen, A_med, A_lo, A_hi, B_med, B_lo, B_hi, C_med, C_lo, C_hi) |> pivot_longer(-pathogen, names_to = c("arm", ".value"), @@ -694,142 +678,237 @@ forest_p50 <- wide_best |> forest_p95 <- wide_best |> select(pathogen, - A_med = A_q95, A_lo = A_q95_lo, A_hi = A_q95_hi, - B_med = B_q95, B_lo = B_q95_lo, B_hi = B_q95_hi, - C_med = C_q95, C_lo = C_q95_lo, C_hi = C_q95_hi) |> + A_med = A_q95, A_lo = A_q95_lo, A_hi = A_q95_hi, + B_med = B_q95, B_lo = B_q95_lo, B_hi = B_q95_hi, + C_med = C_q95, C_lo = C_q95_lo, C_hi = C_q95_hi) |> pivot_longer(-pathogen, names_to = c("arm", ".value"), names_pattern = "^(.)_(.*)") |> filter(!is.na(med)) |> mutate(arm = factor(arm, levels = c("A", "B", "C")), metric = "95th percentile (P95)") -forest_both <- bind_rows(forest_p50, forest_p95) |> - mutate(metric = factor(metric, levels = c("Median (P50)", "95th percentile (P95)"))) - -# Dumbbell backbone segments for both metrics. +# Dumbbell backbone segments for both facets (A-to-B range, un-dodged). dumbbell_segs <- bind_rows( - wide_best |> - filter(!is.na(A_med), !is.na(B_med)) |> + wide_best |> filter(!is.na(A_med), !is.na(B_med)) |> transmute(pathogen, x = A_med, xend = B_med, metric = "Median (P50)"), - wide_best |> - filter(!is.na(A_q95), !is.na(B_q95)) |> + wide_best |> filter(!is.na(A_q95), !is.na(B_q95)) |> transmute(pathogen, x = A_q95, xend = B_q95, metric = "95th percentile (P95)") ) |> mutate(metric = factor(metric, levels = c("Median (P50)", "95th percentile (P95)"))) -# Dodging width and size scales. -COMB_DODGE <- 0.5 +# Dodging width and size scales (used by panels A and C). +COMB_DODGE <- 0.5 ARM_SIZES_C <- c("A" = 1.8, "B" = 1.8, "C" = 2.8) ARM_EBW_C <- c("A" = 0.35, "B" = 0.35, "C" = 0.75) -# ── 8c. Panel A: faceted dumbbell forest (P50 | P95) ───────────────────────── +# ── 8d. Extract mu0 CrI widths from Stan fits ───────────────────────────────── +# +# mu0 is the population-level location parameter (distribution's log scale). +# Its 95% CrI width captures pure estimation uncertainty, independent of τ. +# Ratio arm / C > 1 means arm C has tighter knowledge of the population mean. + +.extract_mu0_cri <- function(fit_slot) { + empty <- c(med = NA_real_, lo = NA_real_, hi = NA_real_, width = NA_real_) + if (is.null(fit_slot) || isTRUE(fit_slot$skipped) || is.null(fit_slot$fit)) + return(empty) + draws <- tryCatch( + rstan::extract(fit_slot$fit, pars = "mu0")$mu0, + error = function(e) NULL + ) + if (is.null(draws) || length(draws) == 0L) return(empty) + q <- quantile(draws, c(0.025, 0.5, 0.975), na.rm = TRUE) + c(med = q[2L], lo = q[1L], hi = q[3L], width = q[3L] - q[1L]) +} + +DISPLAY_TO_INTERNAL <- setNames(names(MAIN_DIST_TO_DISPLAY), MAIN_DIST_TO_DISPLAY) +ARM_FIT_KEYS <- c("A" = "individual_only", "B" = "summary_only", + "C" = "federated") + +mu0_tbl <- do.call(rbind, lapply(best_dist_tbl$pathogen, function(p) { + bd_int <- DISPLAY_TO_INTERNAL[[ best_dist_tbl$best_dist[best_dist_tbl$pathogen == p] ]] + if (is.na(bd_int)) return(NULL) + do.call(rbind, lapply(c("A", "B", "C"), function(arm) { + slot <- ablation_fits[[p]][[ ARM_FIT_KEYS[[arm]] ]][[ bd_int ]] + mt <- .extract_mu0_cri(slot) + data.frame(pathogen = p, arm = arm, + mu0_med = mt["med"], mu0_lo = mt["lo"], + mu0_hi = mt["hi"], mu0_width = mt["width"], + row.names = NULL, stringsAsFactors = FALSE) + })) +})) |> + filter(!is.na(mu0_width)) |> + mutate(arm = factor(arm, levels = c("A", "B", "C")), + pathogen = factor(pathogen, levels = levels(wide_best$pathogen))) + +# mu0 CrI width ratio: arm / C (> 1 → arm C is more precise about μ). +mu0_ratio <- mu0_tbl |> + select(pathogen, arm, mu0_width) |> + pivot_wider(names_from = arm, values_from = mu0_width, + names_prefix = "w_") |> + filter(!is.na(w_C)) |> + pivot_longer(cols = c(w_A, w_B), names_to = "arm_key", values_to = "ratio") |> + filter(!is.na(ratio)) |> + mutate( + ratio = ratio / w_C, + comparison = factor(arm_key, + levels = c("w_A", "w_B"), + labels = c("Individual-level only", + "Summary-statistics only")) + ) + +# ── 8e. Parse τ from comparison_tbl ────────────────────────────────────────── +# +# τ is the between-study heterogeneity SD. Only available when n_datasets ≥ 5 +# in that arm (stored as "— (n<5)" otherwise); missing rows are dropped. + +tau_long <- wide_best |> + rowwise() |> + mutate( + A_tau_med = .parse_cri(A_tau)[["med"]], + A_tau_lo = .parse_cri(A_tau)[["lo"]], + A_tau_hi = .parse_cri(A_tau)[["hi"]], + B_tau_med = .parse_cri(B_tau)[["med"]], + B_tau_lo = .parse_cri(B_tau)[["lo"]], + B_tau_hi = .parse_cri(B_tau)[["hi"]], + C_tau_med = .parse_cri(C_tau)[["med"]], + C_tau_lo = .parse_cri(C_tau)[["lo"]], + C_tau_hi = .parse_cri(C_tau)[["hi"]] + ) |> + ungroup() |> + select(pathogen, matches("^[ABC]_tau_(med|lo|hi)$")) |> + pivot_longer(-pathogen, + names_to = c("arm", ".value"), + names_pattern = "^(.)_tau_(.*)$") |> + filter(!is.na(med)) |> + mutate(arm = factor(arm, levels = c("A", "B", "C"))) + +# ── 8f. Unified arm labels for legend merging ───────────────────────────────── +# +# Convert the "A"/"B"/"C" factor to full label text in all forest/tau/mu0 data. +# All panels then use the same colour/shape scale keyed by label text, so +# patchwork's guides = "collect" produces a single merged legend. + +ARM_COLOURS_FULL <- setNames(ARM_COLOURS, ARM_LABELS[names(ARM_COLOURS)]) +ARM_SHAPES_FULL <- setNames(ARM_SHAPES, ARM_LABELS[names(ARM_SHAPES)]) +ARM_SIZES_C_FULL <- setNames(ARM_SIZES_C, ARM_LABELS[names(ARM_SIZES_C)]) +ARM_EBW_C_FULL <- setNames(ARM_EBW_C, ARM_LABELS[names(ARM_EBW_C)]) + +COMP_COLOURS <- ARM_COLOURS_FULL[c("Individual-level only", + "Summary-statistics only")] +COMP_SHAPES <- ARM_SHAPES_FULL[ c("Individual-level only", + "Summary-statistics only")] + +.relabel_arm <- function(arm_fac) { + factor(ARM_LABELS[as.character(arm_fac)], levels = unname(ARM_LABELS)) +} + +forest_both <- bind_rows(forest_p50, forest_p95) |> + mutate(metric = factor(metric, + levels = c("Median (P50)", "95th percentile (P95)")), + arm = .relabel_arm(arm)) + +tau_long <- tau_long |> mutate(arm = .relabel_arm(arm)) +mu0_tbl <- mu0_tbl |> mutate(arm = .relabel_arm(arm)) + +# ── 8g. Build all four panels ───────────────────────────────────────────────── + +# Shared y-axis theme for the narrow strip panels (labels suppressed). +y_shared <- list( + scale_y_discrete(drop = FALSE), + theme_ablation(), + theme(axis.text.y = element_blank(), + axis.ticks.y = element_blank(), + axis.title.y = element_blank()) +) +# Panel A: faceted dumbbell forest (P50 | P95). pA <- ggplot(forest_both, aes(x = med, y = pathogen, colour = arm, shape = arm)) + - # Grey dumbbell backbone (A-to-B range), drawn at the un-dodged y position. geom_segment( - data = dumbbell_segs, + data = dumbbell_segs, aes(x = x, xend = xend, y = pathogen, yend = pathogen), - colour = "grey78", - linewidth = 0.9, - lineend = "round", + colour = "grey78", linewidth = 0.9, lineend = "round", inherit.aes = FALSE ) + - # 95 % posterior CrI bars, dodged vertically. geom_errorbarh( aes(xmin = lo, xmax = hi, linewidth = arm), height = 0, position = position_dodge(width = COMB_DODGE) ) + - # Point estimates, same dodge. - geom_point( - aes(size = arm), - position = position_dodge(width = COMB_DODGE) - ) + + geom_point(aes(size = arm), position = position_dodge(width = COMB_DODGE)) + facet_wrap(~ metric, nrow = 1L, scales = "free_x") + - scale_x_log10( - breaks = c(1, 2, 5, 10, 20, 50, 100), - labels = c("1", "2", "5", "10", "20", "50", "100") - ) + + scale_x_log10(breaks = c(1, 2, 5, 10, 20, 50, 100), + labels = c("1", "2", "5", "10", "20", "50", "100")) + scale_y_discrete(labels = dist_label_lookup) + - scale_colour_manual(values = ARM_COLOURS, labels = ARM_LABELS) + - scale_shape_manual( values = ARM_SHAPES, labels = ARM_LABELS) + - scale_size_manual( values = ARM_SIZES_C, labels = ARM_LABELS) + - scale_linewidth_manual(values = ARM_EBW_C, labels = ARM_LABELS) + + scale_colour_manual(values = ARM_COLOURS_FULL) + + scale_shape_manual( values = ARM_SHAPES_FULL) + + scale_size_manual( values = ARM_SIZES_C_FULL) + + scale_linewidth_manual(values = ARM_EBW_C_FULL) + labs(x = "Days (log scale)", subtitle = "Best-fitting distribution per pathogen (main analysis model weights)") + theme_ablation() + - guides( - colour = guide_legend(override.aes = list(size = 2.5)), - size = "none", - linewidth = "none", - shape = "none" - ) - -# ── 8d. Panels B and C: gain metrics ───────────────────────────────────────── - -# Colours/shapes for the two comparison directions, using the same names as -# ARM_LABELS so the collected legend is unified across all three panels. -COMP_COLOURS <- c( - "Individual-level only" = ARM_COLOURS[["A"]], - "Summary-statistics only" = ARM_COLOURS[["B"]] -) -COMP_SHAPES <- c("Individual-level only" = 19, "Summary-statistics only" = 17) - -gain_best <- wide_best |> - select(pathogen, ratio_CA, ratio_CB, JS_CA, JS_CB) |> - pivot_longer( - cols = -pathogen, - names_to = c("metric", "comparison"), - names_pattern = "^(.*)_(C[AB])$" - ) |> - filter(!is.na(value)) |> - mutate( - comparison = factor( - comparison, - levels = c("CA", "CB"), - labels = c("Individual-level only", "Summary-statistics only") - ) - ) - -# Shared y scale for panels B and C — same factor, labels suppressed. -y_shared <- list( - scale_y_discrete(drop = FALSE), - theme_ablation(), - theme( - axis.text.y = element_blank(), - axis.ticks.y = element_blank(), - axis.title.y = element_blank() - ) -) - -pB <- ggplot( - filter(gain_best, metric == "ratio"), - aes(x = value, y = pathogen, colour = comparison, shape = comparison) - ) + + guides(colour = guide_legend(override.aes = list(size = 2.5)), + size = "none", + linewidth = "none", + shape = "none") + +# Panel B: μ₀ CrI width ratio — pure information gain. +# Values > 1 indicate arm C has tighter posterior for the population mean. +pB <- ggplot(mu0_ratio, + aes(x = ratio, y = pathogen, + colour = comparison, shape = comparison)) + geom_vline(xintercept = 1, linetype = "dashed", colour = "grey50", linewidth = 0.4) + geom_point(size = 2.2) + scale_colour_manual(values = COMP_COLOURS) + scale_shape_manual( values = COMP_SHAPES) + - labs(x = "Interval ratio\n(arm width / C width)") + + labs(x = expression(mu[0]~"CrI ratio (arm / C)"), + subtitle = "Information gain") + y_shared -pC <- ggplot( - filter(gain_best, metric == "JS"), - aes(x = value, y = pathogen, colour = comparison, shape = comparison) - ) + - geom_vline(xintercept = 0, linetype = "dashed", +# Panel C: τ per arm — between-study heterogeneity. +# Arm C may detect larger τ when cross-data-type contrast reveals more +# between-study variation; its τ posterior is also better estimated +# (narrower CrI) due to more studies. +pC <- ggplot(tau_long, + aes(x = med, y = pathogen, colour = arm, shape = arm)) + + geom_errorbarh(aes(xmin = lo, xmax = hi), + height = 0, linewidth = 0.35, + position = position_dodge(width = COMB_DODGE)) + + geom_point(size = 2.0, position = position_dodge(width = COMB_DODGE)) + + scale_colour_manual(values = ARM_COLOURS_FULL) + + scale_shape_manual( values = ARM_SHAPES_FULL) + + labs(x = expression(tau~"(heterogeneity SD)"), + subtitle = "Between-study heterogeneity") + + y_shared + +# Panel D: predictive CrI ratio — the combined (confounded) signal. +pred_ratio <- wide_best |> + select(pathogen, ratio_CA, ratio_CB) |> + pivot_longer(cols = c(ratio_CA, ratio_CB), + names_to = "comparison", + values_to = "ratio") |> + filter(!is.na(ratio)) |> + mutate(comparison = factor(comparison, + levels = c("ratio_CA", "ratio_CB"), + labels = c("Individual-level only", + "Summary-statistics only"))) + +pD <- ggplot(pred_ratio, + aes(x = ratio, y = pathogen, + colour = comparison, shape = comparison)) + + geom_vline(xintercept = 1, linetype = "dashed", colour = "grey50", linewidth = 0.4) + geom_point(size = 2.2) + scale_colour_manual(values = COMP_COLOURS) + scale_shape_manual( values = COMP_SHAPES) + - labs(x = "JS divergence (bits)") + + labs(x = "Predictive CrI ratio (arm / C)", + subtitle = "Combined effect") + y_shared -# ── 8d. Assemble and save ───────────────────────────────────────────────────── +# ── 8h. Assemble and save ───────────────────────────────────────────────────── -fig4 <- (pA | pB | pC) + - plot_layout(widths = c(2, 1, 1), guides = "collect") + +fig4 <- (pA | pB | pC | pD) + + plot_layout(widths = c(2, 1, 1, 1), guides = "collect") + plot_annotation(tag_levels = "A") & theme(legend.position = "bottom", plot.tag = element_text(face = "bold", size = 10)) @@ -837,9 +916,9 @@ fig4 <- (pA | pB | pC) + fig4_height <- max(4, n_pathogens * FIG_HEIGHT_ROW + 1.8) ggsave(file.path(OUTPUT_DIR, "fig_ablation_combined.pdf"), - fig4, width = FIG_WIDTH_WIDE, height = fig4_height, device = "pdf") + fig4, width = FIG_WIDTH_WIDE * 1.35, height = fig4_height, device = "pdf") ggsave(file.path(OUTPUT_DIR, "fig_ablation_combined.png"), - fig4, width = FIG_WIDTH_WIDE, height = fig4_height, dpi = FIG_DPI) + fig4, width = FIG_WIDTH_WIDE * 1.35, height = fig4_height, dpi = FIG_DPI) message("Figure 4 (combined) saved.") From 2b62f384d11f278e10a986e542d74cdd4821e6fe Mon Sep 17 00:00:00 2001 From: cm401 Date: Sun, 10 May 2026 21:09:44 +0800 Subject: [PATCH 09/18] Fix mu0 extraction: return data.frame with unname()d quantile values rstan::extract() returns draws whose quantile() output carries names like "2.5%" and "97.5%". Passing those directly as named-vector arguments to data.frame() caused the quantile names to become row names, which then propagated as corrupted column names when rows were stacked with do.call(rbind,...). Result: all mu0_width values were silently NA after filter(!is.na(mu0_width)). Fix: .extract_mu0_cri() now returns a one-row data.frame built from unname(quantile(...)) values. The per-arm rows are assembled with cbind() rather than named-vector data.frame() arguments. Confirmed: all 30 (pathogen x arm) combinations now have valid mu0_width values. Co-Authored-By: Claude Sonnet 4.6 --- analysis/plot_data_format_ablation.R | 42 ++++++++++++++++------------ 1 file changed, 24 insertions(+), 18 deletions(-) diff --git a/analysis/plot_data_format_ablation.R b/analysis/plot_data_format_ablation.R index c472a5b..cb8ef84 100644 --- a/analysis/plot_data_format_ablation.R +++ b/analysis/plot_data_format_ablation.R @@ -706,8 +706,11 @@ ARM_EBW_C <- c("A" = 0.35, "B" = 0.35, "C" = 0.75) # Its 95% CrI width captures pure estimation uncertainty, independent of τ. # Ratio arm / C > 1 means arm C has tighter knowledge of the population mean. +# Returns a one-row data.frame with unname()d values to avoid the +# "97.5%" row-name corruption that occurs when quantile output is passed +# directly to data.frame() as a named numeric vector. .extract_mu0_cri <- function(fit_slot) { - empty <- c(med = NA_real_, lo = NA_real_, hi = NA_real_, width = NA_real_) + empty <- data.frame(mu0_lo = NA_real_, mu0_hi = NA_real_, mu0_width = NA_real_) if (is.null(fit_slot) || isTRUE(fit_slot$skipped) || is.null(fit_slot$fit)) return(empty) draws <- tryCatch( @@ -715,8 +718,8 @@ ARM_EBW_C <- c("A" = 0.35, "B" = 0.35, "C" = 0.75) error = function(e) NULL ) if (is.null(draws) || length(draws) == 0L) return(empty) - q <- quantile(draws, c(0.025, 0.5, 0.975), na.rm = TRUE) - c(med = q[2L], lo = q[1L], hi = q[3L], width = q[3L] - q[1L]) + q <- unname(quantile(draws, c(0.025, 0.975), na.rm = TRUE)) + data.frame(mu0_lo = q[1L], mu0_hi = q[2L], mu0_width = q[2L] - q[1L]) } DISPLAY_TO_INTERNAL <- setNames(names(MAIN_DIST_TO_DISPLAY), MAIN_DIST_TO_DISPLAY) @@ -728,11 +731,10 @@ mu0_tbl <- do.call(rbind, lapply(best_dist_tbl$pathogen, function(p) { if (is.na(bd_int)) return(NULL) do.call(rbind, lapply(c("A", "B", "C"), function(arm) { slot <- ablation_fits[[p]][[ ARM_FIT_KEYS[[arm]] ]][[ bd_int ]] - mt <- .extract_mu0_cri(slot) - data.frame(pathogen = p, arm = arm, - mu0_med = mt["med"], mu0_lo = mt["lo"], - mu0_hi = mt["hi"], mu0_width = mt["width"], - row.names = NULL, stringsAsFactors = FALSE) + cbind( + data.frame(pathogen = p, arm = arm, stringsAsFactors = FALSE), + .extract_mu0_cri(slot) + ) })) })) |> filter(!is.na(mu0_width)) |> @@ -740,20 +742,24 @@ mu0_tbl <- do.call(rbind, lapply(best_dist_tbl$pathogen, function(p) { pathogen = factor(pathogen, levels = levels(wide_best$pathogen))) # mu0 CrI width ratio: arm / C (> 1 → arm C is more precise about μ). +# mu0 CrI width ratio: arm / C (> 1 → arm C is more precise about μ). +# Built via left_join to avoid pivot_wider column-existence issues. +mu0_c <- mu0_tbl |> + filter(as.character(arm) == "C") |> + select(pathogen, c_width = mu0_width) + mu0_ratio <- mu0_tbl |> - select(pathogen, arm, mu0_width) |> - pivot_wider(names_from = arm, values_from = mu0_width, - names_prefix = "w_") |> - filter(!is.na(w_C)) |> - pivot_longer(cols = c(w_A, w_B), names_to = "arm_key", values_to = "ratio") |> - filter(!is.na(ratio)) |> + filter(as.character(arm) %in% c("A", "B")) |> + left_join(mu0_c, by = "pathogen") |> + filter(!is.na(mu0_width), !is.na(c_width)) |> mutate( - ratio = ratio / w_C, - comparison = factor(arm_key, - levels = c("w_A", "w_B"), + ratio = mu0_width / c_width, + comparison = factor(as.character(arm), + levels = c("A", "B"), labels = c("Individual-level only", "Summary-statistics only")) - ) + ) |> + select(pathogen, comparison, ratio) # ── 8e. Parse τ from comparison_tbl ────────────────────────────────────────── # From aacb043880d33e182f7679c4040d3aae497298f4 Mon Sep 17 00:00:00 2001 From: cm401 Date: Sun, 10 May 2026 21:24:01 +0800 Subject: [PATCH 10/18] Rename ablation arm labels from A/B/C to I/S/F throughout plot script Replace "Individual-level only", "Summary-statistics only", "Federated" with "(I)", "(S)", "(F)" suffixes everywhere in plot_data_format_ablation.R to eliminate the clash with patchwork's A/B/C/D panel labels. Changes: - ARM_LABELS updated with (I)/(S)/(F) suffixes - COMP_COLOURS / COMP_SHAPES now derived dynamically via ARM_LABELS[c("A","B")] so they stay in sync with ARM_LABELS automatically - mu0_ratio and pred_ratio comparison factor labels now use ARM_LABELS[c("A","B")] instead of hardcoded strings - Figure 2 and supplementary gain figures updated: "vs A (individual-level)" -> "vs (I) individual-level" "vs B (summary-stats)" -> "vs (S) summary-stats" - Caption updated: "federated (C)" -> "federated (F)" Co-Authored-By: Claude Sonnet 4.6 --- analysis/plot_data_format_ablation.R | 42 +++++++++++++--------------- 1 file changed, 19 insertions(+), 23 deletions(-) diff --git a/analysis/plot_data_format_ablation.R b/analysis/plot_data_format_ablation.R index cb8ef84..3718b55 100644 --- a/analysis/plot_data_format_ablation.R +++ b/analysis/plot_data_format_ablation.R @@ -63,9 +63,9 @@ ARM_COLOURS <- c( "C" = "#009E73" # green ) ARM_LABELS <- c( - "A" = "Individual-level only", - "B" = "Summary-statistics only", - "C" = "Federated" + "A" = "Individual-level only (I)", + "B" = "Summary-statistics only (S)", + "C" = "Federated (F)" ) ARM_SHAPES <- c("A" = 19, "B" = 17, "C" = 18) # circle, triangle, diamond @@ -193,7 +193,7 @@ fig1 <- ggplot(forest_long, facet_wrap(~ dist, nrow = 1L, scales = "free_x") + labs( x = "Posterior predictive median (days, log scale)", - caption = "Horizontal bars: 95% credible intervals. Dashed line: federated (C) point estimate." + caption = "Horizontal bars: 95% credible intervals. Dashed line: federated (F) point estimate." ) + theme_ablation() + guides(colour = guide_legend(override.aes = list(size = 2.5))) @@ -231,7 +231,7 @@ gain_df <- comparison_tbl |> mutate( comparison = factor(comparison, levels = c("CA", "CB"), - labels = c("vs A (individual-level)", "vs B (summary-stats)")), + labels = c("vs (I) individual-level", "vs (S) summary-stats")), metric = factor(metric, levels = c("ratio", "JS", "OVL"), labels = c("Interval ratio\n(arm / federated)", @@ -252,12 +252,12 @@ y_scale <- scale_y_discrete(drop = FALSE) y_scale + labs(x = x_lab) + scale_colour_manual( - values = c("vs A (individual-level)" = ARM_COLOURS[["A"]], - "vs B (summary-stats)" = ARM_COLOURS[["B"]]) + values = c("vs (I) individual-level" = ARM_COLOURS[["A"]], + "vs (S) summary-stats" = ARM_COLOURS[["B"]]) ) + scale_shape_manual( - values = c("vs A (individual-level)" = 19, - "vs B (summary-stats)" = 17) + values = c("vs (I) individual-level" = 19, + "vs (S) summary-stats" = 17) ) + theme_ablation() + theme(legend.position = if (show_y) "none" else "bottom") @@ -495,8 +495,8 @@ gain_all_dists <- comparison_tbl |> mutate( comparison = factor(comparison, levels = c("CA", "CB"), - labels = c("vs A (individual-level)", - "vs B (summary-stats)")), + labels = c("vs (I) individual-level", + "vs (S) summary-stats")), metric = factor(metric, levels = c("ratio", "JS", "OVL"), labels = c("Interval ratio", @@ -522,12 +522,12 @@ fig_supp <- ggplot(gain_all_dists, linewidth = 0.35, ) + scale_colour_manual( - values = c("vs A (individual-level)" = ARM_COLOURS[["A"]], - "vs B (summary-stats)" = ARM_COLOURS[["B"]]) + values = c("vs (I) individual-level" = ARM_COLOURS[["A"]], + "vs (S) summary-stats" = ARM_COLOURS[["B"]]) ) + scale_shape_manual( - values = c("vs A (individual-level)" = 19, - "vs B (summary-stats)" = 17) + values = c("vs (I) individual-level" = 19, + "vs (S) summary-stats" = 17) ) + facet_grid(dist ~ metric, scales = "free_x") + labs(x = NULL) + @@ -756,8 +756,7 @@ mu0_ratio <- mu0_tbl |> ratio = mu0_width / c_width, comparison = factor(as.character(arm), levels = c("A", "B"), - labels = c("Individual-level only", - "Summary-statistics only")) + labels = ARM_LABELS[c("A", "B")]) ) |> select(pathogen, comparison, ratio) @@ -798,10 +797,8 @@ ARM_SHAPES_FULL <- setNames(ARM_SHAPES, ARM_LABELS[names(ARM_SHAPES)]) ARM_SIZES_C_FULL <- setNames(ARM_SIZES_C, ARM_LABELS[names(ARM_SIZES_C)]) ARM_EBW_C_FULL <- setNames(ARM_EBW_C, ARM_LABELS[names(ARM_EBW_C)]) -COMP_COLOURS <- ARM_COLOURS_FULL[c("Individual-level only", - "Summary-statistics only")] -COMP_SHAPES <- ARM_SHAPES_FULL[ c("Individual-level only", - "Summary-statistics only")] +COMP_COLOURS <- ARM_COLOURS_FULL[ARM_LABELS[c("A", "B")]] +COMP_SHAPES <- ARM_SHAPES_FULL[ ARM_LABELS[c("A", "B")]] .relabel_arm <- function(arm_fac) { factor(ARM_LABELS[as.character(arm_fac)], levels = unname(ARM_LABELS)) @@ -896,8 +893,7 @@ pred_ratio <- wide_best |> filter(!is.na(ratio)) |> mutate(comparison = factor(comparison, levels = c("ratio_CA", "ratio_CB"), - labels = c("Individual-level only", - "Summary-statistics only"))) + labels = ARM_LABELS[c("A", "B")])) pD <- ggplot(pred_ratio, aes(x = ratio, y = pathogen, From d05ff51f5ca565fb1154adbfd5134b00e1795c66 Mon Sep 17 00:00:00 2001 From: cm401 Date: Sun, 10 May 2026 21:35:41 +0800 Subject: [PATCH 11/18] Three panel improvements for Figure 4 (ablation combined plot) 1. Extract tau from Stan fits directly (Panel C) Replace the parsed-string tau_long (which silently dropped arms with n_datasets < 5) with Stan-extracted tau using .extract_tau_cri(), the same pattern used for mu0. Tau is always a parameter in the hierarchical model; arms with few datasets show wide prior-dominated CrIs, which is honest and ensures all three arms appear for every pathogen. 2. Make Panel A explicitly about posterior predictive CrI Change x-axis label to state "Posterior predictive estimate" and clarify that points are median/P95 and error bars are 95% credible intervals. 3. Annotate Panel C subtitle Add "(wide CrI: few datasets in arm)" so readers understand that width encodes estimation reliability, not just biological heterogeneity. Co-Authored-By: Claude Sonnet 4.6 --- analysis/plot_data_format_ablation.R | 65 +++++++++++++++++----------- 1 file changed, 40 insertions(+), 25 deletions(-) diff --git a/analysis/plot_data_format_ablation.R b/analysis/plot_data_format_ablation.R index 3718b55..0590bba 100644 --- a/analysis/plot_data_format_ablation.R +++ b/analysis/plot_data_format_ablation.R @@ -722,6 +722,23 @@ ARM_EBW_C <- c("A" = 0.35, "B" = 0.35, "C" = 0.75) data.frame(mu0_lo = q[1L], mu0_hi = q[2L], mu0_width = q[2L] - q[1L]) } +# Extracts τ posterior quantiles (2.5%, 50%, 97.5%) from a Stan fit slot. +# Returns a one-row data.frame with (med, lo, hi). τ is always a parameter +# in the hierarchical Stan model; arms with few datasets will have +# prior-dominated (wide) CrIs, which is honest and visually informative. +.extract_tau_cri <- function(fit_slot) { + empty <- data.frame(med = NA_real_, lo = NA_real_, hi = NA_real_) + if (is.null(fit_slot) || isTRUE(fit_slot$skipped) || is.null(fit_slot$fit)) + return(empty) + draws <- tryCatch( + rstan::extract(fit_slot$fit, pars = "tau")$tau, + error = function(e) NULL + ) + if (is.null(draws) || length(draws) == 0L) return(empty) + q <- unname(quantile(draws, c(0.025, 0.5, 0.975), na.rm = TRUE)) + data.frame(med = q[2L], lo = q[1L], hi = q[3L]) +} + DISPLAY_TO_INTERNAL <- setNames(names(MAIN_DIST_TO_DISPLAY), MAIN_DIST_TO_DISPLAY) ARM_FIT_KEYS <- c("A" = "individual_only", "B" = "summary_only", "C" = "federated") @@ -760,31 +777,29 @@ mu0_ratio <- mu0_tbl |> ) |> select(pathogen, comparison, ratio) -# ── 8e. Parse τ from comparison_tbl ────────────────────────────────────────── +# ── 8e. Extract τ from Stan fits (all arms) ─────────────────────────────────── # -# τ is the between-study heterogeneity SD. Only available when n_datasets ≥ 5 -# in that arm (stored as "— (n<5)" otherwise); missing rows are dropped. - -tau_long <- wide_best |> - rowwise() |> - mutate( - A_tau_med = .parse_cri(A_tau)[["med"]], - A_tau_lo = .parse_cri(A_tau)[["lo"]], - A_tau_hi = .parse_cri(A_tau)[["hi"]], - B_tau_med = .parse_cri(B_tau)[["med"]], - B_tau_lo = .parse_cri(B_tau)[["lo"]], - B_tau_hi = .parse_cri(B_tau)[["hi"]], - C_tau_med = .parse_cri(C_tau)[["med"]], - C_tau_lo = .parse_cri(C_tau)[["lo"]], - C_tau_hi = .parse_cri(C_tau)[["hi"]] - ) |> - ungroup() |> - select(pathogen, matches("^[ABC]_tau_(med|lo|hi)$")) |> - pivot_longer(-pathogen, - names_to = c("arm", ".value"), - names_pattern = "^(.)_tau_(.*)$") |> +# τ is the between-study heterogeneity SD. It is always sampled in the +# hierarchical Stan model, so we extract it directly from each fit rather than +# relying on the comparison_tbl formatted string (which suppresses τ when +# n_datasets < 5). Arms with few datasets show wide CrIs, reflecting +# prior-dominated estimates; this is honest and lets all three arms appear in +# Panel C for every pathogen. + +tau_long <- do.call(rbind, lapply(best_dist_tbl$pathogen, function(p) { + bd_int <- DISPLAY_TO_INTERNAL[[ best_dist_tbl$best_dist[best_dist_tbl$pathogen == p] ]] + if (is.na(bd_int)) return(NULL) + do.call(rbind, lapply(c("A", "B", "C"), function(arm) { + slot <- ablation_fits[[p]][[ ARM_FIT_KEYS[[arm]] ]][[ bd_int ]] + cbind( + data.frame(pathogen = p, arm = arm, stringsAsFactors = FALSE), + .extract_tau_cri(slot) + ) + })) +})) |> filter(!is.na(med)) |> - mutate(arm = factor(arm, levels = c("A", "B", "C"))) + mutate(arm = factor(arm, levels = c("A", "B", "C")), + pathogen = factor(pathogen, levels = levels(wide_best$pathogen))) # ── 8f. Unified arm labels for legend merging ───────────────────────────────── # @@ -846,7 +861,7 @@ pA <- ggplot(forest_both, scale_shape_manual( values = ARM_SHAPES_FULL) + scale_size_manual( values = ARM_SIZES_C_FULL) + scale_linewidth_manual(values = ARM_EBW_C_FULL) + - labs(x = "Days (log scale)", + labs(x = "Posterior predictive estimate (days, log scale)\nPoints: median or 95th percentile; error bars: 95% credible interval", subtitle = "Best-fitting distribution per pathogen (main analysis model weights)") + theme_ablation() + guides(colour = guide_legend(override.aes = list(size = 2.5)), @@ -881,7 +896,7 @@ pC <- ggplot(tau_long, scale_colour_manual(values = ARM_COLOURS_FULL) + scale_shape_manual( values = ARM_SHAPES_FULL) + labs(x = expression(tau~"(heterogeneity SD)"), - subtitle = "Between-study heterogeneity") + + subtitle = "Between-study heterogeneity\n(wide CrI: few datasets in arm)") + y_shared # Panel D: predictive CrI ratio — the combined (confounded) signal. From 6983a0fad826702954c6aa19299a36374aec5a10 Mon Sep 17 00:00:00 2001 From: cm401 Date: Sun, 10 May 2026 21:43:41 +0800 Subject: [PATCH 12/18] Revert Stan-based tau extraction; filter Figure 4 to tau-complete pathogens Reverts the previous Stan-based tau extraction (which included prior-dominated estimates for n < 5 datasets) and restores the correct parsed-string approach that only shows tau where n_datasets >= 5 in each arm. Adds a new filter step (8e-filter) that restricts all four Figure 4 panels to the subset of pathogens that have tau estimated in every arm (I, S, and F). This ensures Panel C is complete and coherent, and all panels share the same y-axis rows. Figure height is updated to reflect the filtered pathogen count. Also reverts the Panel A x-label and Panel C subtitle changes from the previous commit. Co-Authored-By: Claude Sonnet 4.6 --- analysis/plot_data_format_ablation.R | 99 ++++++++++++++++------------ 1 file changed, 58 insertions(+), 41 deletions(-) diff --git a/analysis/plot_data_format_ablation.R b/analysis/plot_data_format_ablation.R index 0590bba..d273fe2 100644 --- a/analysis/plot_data_format_ablation.R +++ b/analysis/plot_data_format_ablation.R @@ -722,23 +722,6 @@ ARM_EBW_C <- c("A" = 0.35, "B" = 0.35, "C" = 0.75) data.frame(mu0_lo = q[1L], mu0_hi = q[2L], mu0_width = q[2L] - q[1L]) } -# Extracts τ posterior quantiles (2.5%, 50%, 97.5%) from a Stan fit slot. -# Returns a one-row data.frame with (med, lo, hi). τ is always a parameter -# in the hierarchical Stan model; arms with few datasets will have -# prior-dominated (wide) CrIs, which is honest and visually informative. -.extract_tau_cri <- function(fit_slot) { - empty <- data.frame(med = NA_real_, lo = NA_real_, hi = NA_real_) - if (is.null(fit_slot) || isTRUE(fit_slot$skipped) || is.null(fit_slot$fit)) - return(empty) - draws <- tryCatch( - rstan::extract(fit_slot$fit, pars = "tau")$tau, - error = function(e) NULL - ) - if (is.null(draws) || length(draws) == 0L) return(empty) - q <- unname(quantile(draws, c(0.025, 0.5, 0.975), na.rm = TRUE)) - data.frame(med = q[2L], lo = q[1L], hi = q[3L]) -} - DISPLAY_TO_INTERNAL <- setNames(names(MAIN_DIST_TO_DISPLAY), MAIN_DIST_TO_DISPLAY) ARM_FIT_KEYS <- c("A" = "individual_only", "B" = "summary_only", "C" = "federated") @@ -777,29 +760,63 @@ mu0_ratio <- mu0_tbl |> ) |> select(pathogen, comparison, ratio) -# ── 8e. Extract τ from Stan fits (all arms) ─────────────────────────────────── +# ── 8e. Parse τ from comparison_tbl ────────────────────────────────────────── # -# τ is the between-study heterogeneity SD. It is always sampled in the -# hierarchical Stan model, so we extract it directly from each fit rather than -# relying on the comparison_tbl formatted string (which suppresses τ when -# n_datasets < 5). Arms with few datasets show wide CrIs, reflecting -# prior-dominated estimates; this is honest and lets all three arms appear in -# Panel C for every pathogen. - -tau_long <- do.call(rbind, lapply(best_dist_tbl$pathogen, function(p) { - bd_int <- DISPLAY_TO_INTERNAL[[ best_dist_tbl$best_dist[best_dist_tbl$pathogen == p] ]] - if (is.na(bd_int)) return(NULL) - do.call(rbind, lapply(c("A", "B", "C"), function(arm) { - slot <- ablation_fits[[p]][[ ARM_FIT_KEYS[[arm]] ]][[ bd_int ]] - cbind( - data.frame(pathogen = p, arm = arm, stringsAsFactors = FALSE), - .extract_tau_cri(slot) - ) - })) -})) |> +# τ is the between-study heterogeneity SD. Only available when n_datasets ≥ 5 +# in that arm (stored as "— (n<5)" otherwise); missing rows are dropped. + +tau_long <- wide_best |> + rowwise() |> + mutate( + A_tau_med = .parse_cri(A_tau)[["med"]], + A_tau_lo = .parse_cri(A_tau)[["lo"]], + A_tau_hi = .parse_cri(A_tau)[["hi"]], + B_tau_med = .parse_cri(B_tau)[["med"]], + B_tau_lo = .parse_cri(B_tau)[["lo"]], + B_tau_hi = .parse_cri(B_tau)[["hi"]], + C_tau_med = .parse_cri(C_tau)[["med"]], + C_tau_lo = .parse_cri(C_tau)[["lo"]], + C_tau_hi = .parse_cri(C_tau)[["hi"]] + ) |> + ungroup() |> + select(pathogen, matches("^[ABC]_tau_(med|lo|hi)$")) |> + pivot_longer(-pathogen, + names_to = c("arm", ".value"), + names_pattern = "^(.)_tau_(.*)$") |> filter(!is.na(med)) |> - mutate(arm = factor(arm, levels = c("A", "B", "C")), - pathogen = factor(pathogen, levels = levels(wide_best$pathogen))) + mutate(arm = factor(arm, levels = c("A", "B", "C"))) + +# ── 8e-filter. Restrict all panels to pathogens with τ for all three arms ───── +# +# Only pathogens where n_datasets ≥ 5 in every arm have τ estimated. Panels +# B, C, D require this; Panel A is also restricted so all four panels share the +# same y-axis rows. + +tau_complete <- tau_long |> + group_by(pathogen) |> + summarise(n_arms = n_distinct(as.character(arm)), .groups = "drop") |> + filter(n_arms == 3L) |> + pull(pathogen) |> + as.character() + +# Preserve original C-arm ordering (federated median, descending). +tau_levels <- intersect(levels(wide_best$pathogen), tau_complete) + +.restrict_pathogens <- function(df) { + df |> + filter(as.character(pathogen) %in% tau_complete) |> + mutate(pathogen = factor(as.character(pathogen), levels = tau_levels)) +} + +forest_p50 <- .restrict_pathogens(forest_p50) +forest_p95 <- .restrict_pathogens(forest_p95) +dumbbell_segs <- .restrict_pathogens(dumbbell_segs) +tau_long <- .restrict_pathogens(tau_long) +mu0_ratio <- .restrict_pathogens(mu0_ratio) +wide_best <- .restrict_pathogens(wide_best) # feeds pred_ratio in 8g + +dist_label_lookup <- dist_label_lookup[tau_levels] +n_tau_pathogens <- length(tau_levels) # ── 8f. Unified arm labels for legend merging ───────────────────────────────── # @@ -861,7 +878,7 @@ pA <- ggplot(forest_both, scale_shape_manual( values = ARM_SHAPES_FULL) + scale_size_manual( values = ARM_SIZES_C_FULL) + scale_linewidth_manual(values = ARM_EBW_C_FULL) + - labs(x = "Posterior predictive estimate (days, log scale)\nPoints: median or 95th percentile; error bars: 95% credible interval", + labs(x = "Days (log scale)", subtitle = "Best-fitting distribution per pathogen (main analysis model weights)") + theme_ablation() + guides(colour = guide_legend(override.aes = list(size = 2.5)), @@ -896,7 +913,7 @@ pC <- ggplot(tau_long, scale_colour_manual(values = ARM_COLOURS_FULL) + scale_shape_manual( values = ARM_SHAPES_FULL) + labs(x = expression(tau~"(heterogeneity SD)"), - subtitle = "Between-study heterogeneity\n(wide CrI: few datasets in arm)") + + subtitle = "Between-study heterogeneity") + y_shared # Panel D: predictive CrI ratio — the combined (confounded) signal. @@ -930,7 +947,7 @@ fig4 <- (pA | pB | pC | pD) + theme(legend.position = "bottom", plot.tag = element_text(face = "bold", size = 10)) -fig4_height <- max(4, n_pathogens * FIG_HEIGHT_ROW + 1.8) +fig4_height <- max(4, n_tau_pathogens * FIG_HEIGHT_ROW + 1.8) ggsave(file.path(OUTPUT_DIR, "fig_ablation_combined.pdf"), fig4, width = FIG_WIDTH_WIDE * 1.35, height = fig4_height, device = "pdf") From a0a7529cdefec349e533d7ee58a45ceb49f121b4 Mon Sep 17 00:00:00 2001 From: cm401 Date: Mon, 11 May 2026 17:14:42 +0800 Subject: [PATCH 13/18] add simulation study figure --- R/ploting_utils.R | 278 ++++++++++++++++++++++++++++++++++++++++++++++ 1 file changed, 278 insertions(+) diff --git a/R/ploting_utils.R b/R/ploting_utils.R index 60d2ce2..4b1a4f7 100644 --- a/R/ploting_utils.R +++ b/R/ploting_utils.R @@ -1609,6 +1609,284 @@ plot_simulation_study_figure <- function(summary_res, base_size = 9) { fig & scale_shape_manual(values = c("5" = 4, "10" = 3, "20" = 8, "25+" = 5)) } +#' Publication-ready simulation study figure (Figure 2) +#' +#' Two-panel figure (A: coverage, B: bias) showing predictive performance +#' across distribution families, summary types, and number of contributing +#' datasets. Colour is mapped to summary type via the AAAS palette, matching +#' the scheme used elsewhere in the package. +#' +#' @param summary_res A data frame returned by \code{create_results_summary()}. +#' @param base_size Base font size (default 9.5). +#' @return A \code{patchwork} object. +#' @export +plot_simulation_study_figure2 <- function(summary_res, base_size = 9.5) { + + dist_labels <- c( + burr12 = "Burr Type XII", gamma = "Gamma", + gengamma = "Gen. Gamma", lognormal = "Log-normal", + weibull = "Weibull" + ) + summary_labels <- c("1" = "Median+Range", "2" = "Median+IQR", + "3" = "Mean+SD", "5" = "Mixed") + + by_nd <- summary_res %>% + filter(summary_type != 4) %>% + mutate( + dist_label = factor(dist_labels[dist_type], + levels = c("Burr Type XII", "Gamma", "Gen. Gamma", + "Log-normal", "Weibull")), + summary_label = factor(summary_labels[as.character(summary_type)], + levels = c("Median+Range", "Median+IQR", + "Mean+SD", "Mixed")) + ) %>% + group_by(dist_label, summary_label, n_datasets) %>% + summarise( + coverage_p50 = mean(coverage_pred_median, na.rm = TRUE), + coverage_p95 = mean(coverage_pred_q95, na.rm = TRUE), + bias_p50 = median(bias_pred_median, na.rm = TRUE), + bias_p95 = median(bias_pred_q95, na.rm = TRUE), + mean_wis = mean(mean_wis, na.rm = TRUE), + .groups = "drop" + ) + + fig_theme <- theme_minimal(base_size = base_size, base_family = "Helvetica") + + theme( + panel.grid.minor = element_blank(), + panel.grid.major.x = element_blank(), + panel.grid.major.y = element_line(colour = "grey90", linewidth = 0.2), + panel.background = element_rect(fill = "white", colour = NA), + panel.border = element_rect(colour = "grey40", fill = NA, linewidth = 0.45), + panel.spacing.x = unit(5, "pt"), + panel.spacing.y = unit(6, "pt"), + strip.text = element_text(size = 8.5, face = "bold", + colour = "grey15", + margin = margin(b = 3, t = 3)), + strip.background = element_rect(fill = "grey96", colour = NA), + axis.text = element_text(size = 7, colour = "grey25"), + axis.title = element_text(size = 9, colour = "grey10"), + axis.ticks = element_line(colour = "grey40", linewidth = 0.25), + axis.ticks.length = unit(2, "pt"), + legend.title = element_text(size = 8.5, face = "bold"), + legend.text = element_text(size = 7.5), + legend.key.size = unit(12, "pt"), + legend.key.spacing.x = unit(2, "pt"), + plot.tag = element_text(face = "bold", size = 12), + plot.margin = margin(4, 6, 4, 4) + ) + + # ── Panel A: Coverage ────────────────────────────────────────────────────── + cov_long <- dplyr::bind_rows( + by_nd %>% transmute(dist_label, summary_label, n_datasets, + quantile = "Predictive median", value = coverage_p50), + by_nd %>% transmute(dist_label, summary_label, n_datasets, + quantile = "95th percentile", value = coverage_p95) + ) %>% + mutate(quantile = factor(quantile, + levels = c("Predictive median", "95th percentile"))) + + pA <- ggplot(cov_long) + + annotate("rect", xmin = -Inf, xmax = Inf, ymin = 0.90, ymax = 1.0, + fill = "#4CAF50", alpha = 0.12) + + geom_hline(yintercept = 0.95, linewidth = 0.4, colour = "grey40", + linetype = "dotted") + + geom_line(aes(x = n_datasets, y = value, colour = summary_label, + group = summary_label), + linewidth = 0.65, alpha = 0.85) + + geom_point(aes(x = n_datasets, y = value, colour = summary_label), + size = 1.5, alpha = 0.9) + + facet_grid(quantile ~ dist_label) + + ggsci::scale_colour_aaas(name = "Summary type") + + scale_x_continuous(breaks = c(3, 5, 10, 20, 30), limits = c(NA, 33)) + + scale_y_continuous(limits = c(NA, 1.02)) + + xlab(NULL) + ylab("Coverage") + + fig_theme + + theme(legend.position = "none") + + labs(tag = "A") + + # ── Panel B: Bias ────────────────────────────────────────────────────────── + bias_long <- dplyr::bind_rows( + by_nd %>% transmute(dist_label, summary_label, n_datasets, + quantile = "Predictive median", value = bias_p50), + by_nd %>% transmute(dist_label, summary_label, n_datasets, + quantile = "95th percentile", value = bias_p95) + ) %>% + mutate(quantile = factor(quantile, + levels = c("Predictive median", "95th percentile"))) + + pB <- ggplot(bias_long) + + annotate("rect", xmin = -Inf, xmax = Inf, ymin = -0.5, ymax = 0.5, + fill = "#4CAF50", alpha = 0.12) + + geom_hline(yintercept = 0, linewidth = 0.4, colour = "grey40", + linetype = "dotted") + + geom_line(aes(x = n_datasets, y = value, colour = summary_label, + group = summary_label), + linewidth = 0.65, alpha = 0.85) + + geom_point(aes(x = n_datasets, y = value, colour = summary_label), + size = 1.5, alpha = 0.9) + + facet_grid(quantile ~ dist_label, scales = "free_y") + + ggsci::scale_colour_aaas(name = "Summary type") + + scale_x_continuous(breaks = c(3, 5, 10, 20, 30), limits = c(NA, 33)) + + xlab("Number of contributing datasets") + ylab("Bias (days)") + + fig_theme + + theme(strip.text.x = element_blank()) + + labs(tag = "B") + + # ── Assemble ─────────────────────────────────────────────────────────────── + (pA / pB) + + patchwork::plot_layout(heights = c(1, 1), guides = "collect") & + theme(legend.position = "bottom", + legend.box = "horizontal", + legend.margin = margin(t = 3)) +} + + +#' Supplementary information simulation study figure +#' +#' Three-panel figure (A: coverage, B: bias, C: mean WIS) using jittered points +#' coloured by number of datasets and shaped by observations per study. All +#' summary types (including Freq Table) are shown on the x-axis. +#' +#' @param summary_res A data frame returned by \code{create_results_summary()}. +#' @param base_size Base font size (default 10). +#' @return A \code{patchwork} object. +#' @export +plot_simulation_study_si_figure <- function(summary_res, base_size = 10) { + + dist_labels <- c( + burr12 = "Burr Type XII", gamma = "Gamma", + gengamma = "Gen. Gamma", lognormal = "Log-normal", + weibull = "Weibull" + ) + summary_labels <- c( + "1" = "Median+Range", "2" = "Median+IQR", "3" = "Mean+SD", + "4" = "Freq Table", "5" = "Mixed" + ) + + ndataset_colours <- c("<10" = "#2ca02c", "<20" = "#ff7f0e", + "<30" = "#1f77b4", "30+" = "#d62728") + nobs_shapes <- c("5" = 4, "10" = 3, "20" = 8, "25+" = 5) + + # Keep one row per scenario (preserves sub-scenario spread, e.g. Burr XII κ clusters). + # Sensitivity scenarios (TauSens/Mu0Sens) are excluded — they belong to the + # sensitivity figure and would dilute the κ clusters and stretch the y-axis. + ss <- summary_res %>% + filter(!grepl("TauSens|Mu0Sens", scenario_name)) %>% + mutate( + dist_label = factor(dist_labels[dist_type], + levels = c("Burr Type XII", "Gamma", "Gen. Gamma", + "Log-normal", "Weibull")), + summary_label = factor(summary_labels[as.character(summary_type)], + levels = c("Median+Range", "Median+IQR", "Mean+SD", + "Freq Table", "Mixed")), + n_datasets_bin = cut(n_datasets, breaks = c(0, 9, 19, 29, Inf), + labels = c("<10", "<20", "<30", "30+")), + n_obs_bin = cut(n_obs_mean, breaks = c(0, 7.5, 15, 22.5, Inf), + labels = c("5", "10", "20", "25+")) + ) %>% + rename( + coverage_p50 = coverage_pred_median, + coverage_p95 = coverage_pred_q95, + bias_p50 = bias_pred_median, + bias_p95 = bias_pred_q95 + ) + + fig_theme <- theme_minimal(base_size = base_size, base_family = "Helvetica") + + theme( + panel.grid.minor = element_blank(), + panel.grid.major.x = element_blank(), + panel.grid.major.y = element_line(colour = "grey90", linewidth = 0.2), + panel.background = element_rect(fill = "white", colour = NA), + panel.border = element_rect(colour = "grey50", fill = NA, linewidth = 0.4), + panel.spacing.x = unit(6, "pt"), + panel.spacing.y = unit(6, "pt"), + strip.text = element_text(size = 9, face = "bold", colour = "grey15"), + strip.background = element_rect(fill = "grey96", colour = NA), + axis.text.x = element_text(size = 7, colour = "grey25", + angle = 45, hjust = 1), + axis.text.y = element_text(size = 7, colour = "grey25"), + axis.title = element_text(size = 9.5), + axis.ticks = element_line(colour = "grey50", linewidth = 0.25), + axis.ticks.length = unit(2, "pt"), + legend.title = element_text(size = 9, face = "bold"), + legend.text = element_text(size = 8), + legend.key.size = unit(12, "pt"), + plot.tag = element_text(face = "bold", size = 13), + plot.margin = margin(4, 6, 4, 4) + ) + + # ── Panel A: Coverage ────────────────────────────────────────────────────── + cov_long <- dplyr::bind_rows( + ss %>% transmute(dist_label, summary_label, n_datasets_bin, n_obs_bin, + quantile = "Median (P50)", value = coverage_p50), + ss %>% transmute(dist_label, summary_label, n_datasets_bin, n_obs_bin, + quantile = "95th percentile (P95)", value = coverage_p95) + ) %>% + mutate(quantile = factor(quantile, + levels = c("Median (P50)", "95th percentile (P95)"))) + + pA <- ggplot(cov_long, aes(x = summary_label, y = value)) + + geom_hline(yintercept = 0.95, linetype = "dashed", + colour = "#d62728", linewidth = 0.4) + + geom_point(aes(colour = n_datasets_bin, shape = n_obs_bin), + size = 2.2, alpha = 0.8, + position = position_jitter(width = 0.25, height = 0, seed = 42)) + + facet_grid(quantile ~ dist_label) + + scale_colour_manual("N datasets", values = ndataset_colours) + + scale_shape_manual( "N obs per study", values = nobs_shapes) + + scale_y_continuous(limits = c(0, 1.05)) + + xlab(NULL) + ylab("Coverage") + + fig_theme + + theme(legend.position = "none") + + labs(tag = "A") + + # ── Panel B: Bias ────────────────────────────────────────────────────────── + bias_long <- dplyr::bind_rows( + ss %>% transmute(dist_label, summary_label, n_datasets_bin, n_obs_bin, + quantile = "Median (P50)", value = bias_p50), + ss %>% transmute(dist_label, summary_label, n_datasets_bin, n_obs_bin, + quantile = "95th percentile (P95)", value = bias_p95) + ) %>% + mutate(quantile = factor(quantile, + levels = c("Median (P50)", "95th percentile (P95)"))) + + pB <- ggplot(bias_long, aes(x = summary_label, y = value)) + + geom_hline(yintercept = 0, linetype = "dashed", + colour = "#d62728", linewidth = 0.4) + + geom_point(aes(colour = n_datasets_bin, shape = n_obs_bin), + size = 2.2, alpha = 0.8, + position = position_jitter(width = 0.25, height = 0, seed = 42)) + + facet_grid(quantile ~ dist_label, scales = "free_y") + + scale_colour_manual("N datasets", values = ndataset_colours) + + scale_shape_manual( "N obs per study", values = nobs_shapes) + + xlab(NULL) + ylab("Bias (days)") + + fig_theme + + theme(legend.position = "none", + strip.text.x = element_blank()) + + labs(tag = "B") + + # ── Panel C: Mean WIS ────────────────────────────────────────────────────── + pC <- ggplot(ss, aes(x = summary_label, y = mean_wis)) + + geom_point(aes(colour = n_datasets_bin, shape = n_obs_bin), + size = 2.2, alpha = 0.8, + position = position_jitter(width = 0.25, height = 0, seed = 42)) + + facet_wrap(~ dist_label, nrow = 1) + + scale_colour_manual("N datasets", values = ndataset_colours) + + scale_shape_manual( "N obs per study", values = nobs_shapes) + + xlab(NULL) + ylab("Mean WIS") + + fig_theme + + theme(strip.text = element_blank()) + + labs(tag = "C") + + # ── Assemble ─────────────────────────────────────────────────────────────── + (pA / pB / pC) + + patchwork::plot_layout(heights = c(1, 1, 0.55), guides = "collect") & + theme(legend.position = "bottom", + legend.box = "horizontal", + legend.margin = margin(t = 4)) +} + + #' Supplementary simulation study figure — sensitivity scenarios #' #' Same layout as \code{plot_simulation_study_figure()} but restricted to From 00c030ad2c56b4b9f635298069dfc653471635ef Mon Sep 17 00:00:00 2001 From: cm401 Date: Mon, 11 May 2026 17:15:45 +0800 Subject: [PATCH 14/18] Update labels and layout for combined figure --- analysis/plot_data_format_ablation.R | 330 +++++++++++++++------------ 1 file changed, 182 insertions(+), 148 deletions(-) diff --git a/analysis/plot_data_format_ablation.R b/analysis/plot_data_format_ablation.R index d273fe2..1c7ebd4 100644 --- a/analysis/plot_data_format_ablation.R +++ b/analysis/plot_data_format_ablation.R @@ -8,25 +8,25 @@ # # Figure 1 — "Do the arms agree?" # Three-arm forest plot. Pathogens on the y-axis, pred_median on a -# log-scaled x-axis. Arms A, B, C dodged vertically per pathogen, with +# log-scaled x-axis. Arms I, S, F dodged vertically per pathogen, with # horizontal 95% CrI bars. A vertical reference line marks the federated -# (arm C) point estimate. Faceted by distribution family. +# (arm F) point estimate. Faceted by distribution family. # # Figure 2 — "How much does federation help?" # Side-by-side metrics strip using the same pathogen ordering. -# Left panel : interval ratio (A_width / C_width, B_width / C_width); +# Left panel : interval ratio (I_width / F_width, S_width / F_width); # reference line at 1 (equal precision). -# Middle panel: Jensen-Shannon divergence (JS_CA, JS_CB; bits). -# Right panel : overlap coefficient (OVL_CA, OVL_CB). +# Middle panel: Jensen-Shannon divergence (JS_FI, JS_FS; bits). +# Right panel : overlap coefficient (OVL_FI, OVL_FS). # All three panels share the y-axis (pathogens) and are assembled with # patchwork. Restricted to one distribution (PRIMARY_DIST) for clarity; # a supplementary version facets over all distributions. # # Figure 3 — "What does the gain look like in practice?" # Posterior predictive density overlays for the TOP_N_DENSITY pathogens -# with the largest JS_CA divergence (i.e. where individual-level data +# with the largest JS_FI divergence (i.e. where individual-level data # pulls the federated estimate furthest from the summary-stat arm). -# Three shaded ribbons (A, B, C) per pathogen panel, computed from +# Three shaded ribbons (I, S, F) per pathogen panel, computed from # compute_predictive_cdf(). # # Output @@ -41,13 +41,14 @@ library(ggplot2) library(dplyr) library(tidyr) library(patchwork) +library(ggsci) # ── 1. Settings ─────────────────────────────────────────────────────────────── # Distribution shown in Figure 2 (gain strip) and Figure 3 (densities). -PRIMARY_DIST <- "Log-normal" +PRIMARY_DIST <- "Burr XII" -# Number of pathogens shown in Figure 3, chosen by largest JS_CA divergence. +# Number of pathogens shown in Figure 3, chosen by largest JS_FI divergence. TOP_N_DENSITY <- 6L # Save dimensions (inches) and resolution. @@ -56,18 +57,14 @@ FIG_WIDTH_HALF <- 7 FIG_HEIGHT_ROW <- 0.55 # height per pathogen row in Figures 1 & 2 FIG_DPI <- 300 -# Arm colours (Wong colorblind-safe palette). -ARM_COLOURS <- c( - "A" = "#0072B2", # blue - "B" = "#D55E00", # vermilion - "C" = "#009E73" # green -) +# Arm colours from the AAAS palette (consistent with the rest of the repo). +ARM_COLOURS <- setNames(pal_aaas()(3), c("I", "S", "F")) ARM_LABELS <- c( - "A" = "Individual-level only (I)", - "B" = "Summary-statistics only (S)", - "C" = "Federated (F)" + "I" = "Individual-level only (I)", + "S" = "Summary-statistics only (S)", + "F" = "Federated (F)" ) -ARM_SHAPES <- c("A" = 19, "B" = 17, "C" = 18) # circle, triangle, diamond +ARM_SHAPES <- c("I" = 19, "S" = 17, "F" = 18) # circle, triangle, diamond OUTPUT_DIR <- here::here("results", "figures") dir.create(OUTPUT_DIR, showWarnings = FALSE, recursive = TRUE) @@ -87,6 +84,23 @@ ablation_fits <- stored$ablation_fits if (nrow(comparison_tbl) == 0L) stop("comparison_tbl is empty — check that data_format_ablation.R completed.") +# Back-compat: rename legacy A/B/C-prefixed and _CA/_CB-suffixed columns to +# the new I/S/F scheme so an older comparison_tbl still works without rerunning +# data_format_ablation.R. No-op once the RDS has been regenerated. +.rename_arm_cols <- function(nm) { + nm <- sub("^A_", "I_", nm) + nm <- sub("^B_", "S_", nm) + nm <- sub("^C_", "F_", nm) + nm <- sub("^ptail_A_", "ptail_I_", nm) + nm <- sub("^ptail_B_", "ptail_S_", nm) + nm <- sub("^ptail_C_", "ptail_F_", nm) + nm <- sub("_CA(_|$)", "_FI\\1", nm) + nm <- sub("_CB(_|$)", "_FS\\1", nm) + nm <- sub("_AB(_|$)", "_IS\\1", nm) + nm +} +names(comparison_tbl) <- .rename_arm_cols(names(comparison_tbl)) + # ── 3. Shared helpers ───────────────────────────────────────────────────────── @@ -119,9 +133,9 @@ theme_ablation <- function(base_size = 10) { tbl |> filter(.data$dist == !!dist) |> rowwise() |> - mutate(C_med = .parse_cri(C_pred_median)[["med"]]) |> + mutate(F_med = .parse_cri(F_pred_median)[["med"]]) |> ungroup() |> - arrange(C_med) |> + arrange(F_med) |> pull(pathogen) |> unique() } @@ -133,7 +147,7 @@ pathogen_order <- .pathogen_order(comparison_tbl) # # For each arm, parse the formatted CrI string into med/lo/hi and reshape to # long format. Dodge three arms vertically within each pathogen × distribution -# cell. Log-scale x-axis with a vertical reference line at arm C's median. +# cell. Log-scale x-axis with a vertical reference line at arm F's median. .parse_arm_col <- function(tbl, col, arm_label) { tbl |> @@ -151,19 +165,19 @@ pathogen_order <- .pathogen_order(comparison_tbl) } forest_long <- bind_rows( - .parse_arm_col(comparison_tbl, "A_pred_median", "A"), - .parse_arm_col(comparison_tbl, "B_pred_median", "B"), - .parse_arm_col(comparison_tbl, "C_pred_median", "C") + .parse_arm_col(comparison_tbl, "I_pred_median", "I"), + .parse_arm_col(comparison_tbl, "S_pred_median", "S"), + .parse_arm_col(comparison_tbl, "F_pred_median", "F") ) |> filter(!is.na(med)) |> mutate( pathogen = factor(pathogen, levels = pathogen_order), - arm = factor(arm, levels = c("A", "B", "C")) + arm = factor(arm, levels = c("I", "S", "F")) ) -# Federated reference line (arm C) per facet cell. +# Federated reference line (arm F) per facet cell. ref_lines <- forest_long |> - filter(arm == "C") |> + filter(arm == "F") |> select(pathogen, dist, ref_med = med) dodge_height <- 0.55 # vertical spread per pathogen row @@ -171,10 +185,10 @@ dodge_height <- 0.55 # vertical spread per pathogen row fig1 <- ggplot(forest_long, aes(x = med, y = pathogen, colour = arm, shape = arm)) + - # Reference line at arm C's point estimate + # Reference line at arm F's point estimate geom_vline(data = ref_lines, aes(xintercept = ref_med), - colour = ARM_COLOURS[["C"]], linewidth = 0.25, + colour = ARM_COLOURS[["F"]], linewidth = 0.25, linetype = "dashed") + # CrI bars geom_errorbarh(aes(xmin = lo, xmax = hi), @@ -188,7 +202,7 @@ fig1 <- ggplot(forest_long, breaks = c(1, 2, 5, 10, 20, 50), labels = c("1", "2", "5", "10", "20", "50") ) + - scale_colour_manual(values = ARM_COLOURS, labels = ARM_LABELS) + + scale_colour_aaas(labels = ARM_LABELS) + scale_shape_manual(values = ARM_SHAPES, labels = ARM_LABELS) + facet_wrap(~ dist, nrow = 1L, scales = "free_x") + labs( @@ -218,9 +232,9 @@ gain_df <- comparison_tbl |> filter(dist == PRIMARY_DIST) |> mutate(pathogen = factor(pathogen, levels = pathogen_order)) |> select(pathogen, - ratio_CA, ratio_CB, - JS_CA, JS_CB, - OVL_CA, OVL_CB) |> + ratio_FI, ratio_FS, + JS_FI, JS_FS, + OVL_FI, OVL_FS) |> pivot_longer( cols = -pathogen, names_to = c("metric", "comparison"), @@ -230,7 +244,7 @@ gain_df <- comparison_tbl |> filter(!is.na(value)) |> mutate( comparison = factor(comparison, - levels = c("CA", "CB"), + levels = c("FI", "FS"), labels = c("vs (I) individual-level", "vs (S) summary-stats")), metric = factor(metric, levels = c("ratio", "JS", "OVL"), @@ -251,10 +265,7 @@ y_scale <- scale_y_discrete(drop = FALSE) geom_point(size = 2.2, alpha = 0.85) + y_scale + labs(x = x_lab) + - scale_colour_manual( - values = c("vs (I) individual-level" = ARM_COLOURS[["A"]], - "vs (S) summary-stats" = ARM_COLOURS[["B"]]) - ) + + scale_colour_aaas() + scale_shape_manual( values = c("vs (I) individual-level" = 19, "vs (S) summary-stats" = 17) @@ -305,7 +316,7 @@ message("Figure 2 saved.") # ── 6. Figure 3: Posterior predictive density overlays ─────────────────────── # -# Selects the TOP_N_DENSITY pathogens with the largest JS_CA (arm A vs +# Selects the TOP_N_DENSITY pathogens with the largest JS_FI (arm I vs # federated) for PRIMARY_DIST. Calls compute_predictive_cdf() for each of the # three arms and overlays the posterior-median density with a ±95% CrI ribbon. # Falls back gracefully if a fit is absent. @@ -317,13 +328,13 @@ dist_code_map <- c( # Identify the top pathogens to plot. top_pathogens <- comparison_tbl |> - filter(dist == PRIMARY_DIST, !is.na(JS_CA)) |> - arrange(desc(JS_CA)) |> + filter(dist == PRIMARY_DIST, !is.na(JS_FI)) |> + arrange(desc(JS_FI)) |> slice_head(n = TOP_N_DENSITY) |> pull(pathogen) if (length(top_pathogens) == 0L) { - message("No JS_CA values available — skipping Figure 3.") + message("No JS_FI values available — skipping Figure 3.") } else { cdf_dist_name <- dist_code_map[[PRIMARY_DIST]] @@ -333,7 +344,7 @@ if (length(top_pathogens) == 0L) { .x_upper <- function(pathogen_name) { row <- filter(comparison_tbl, pathogen == pathogen_name & dist == PRIMARY_DIST) - vals <- sapply(c("A_pred_q95", "B_pred_q95", "C_pred_q95"), function(col) { + vals <- sapply(c("I_pred_q95", "S_pred_q95", "F_pred_q95"), function(col) { p <- tryCatch(.parse_cri(row[[col]])[["hi"]], error = function(e) NA_real_) p }) @@ -374,7 +385,7 @@ if (length(top_pathogens) == 0L) { for (pg in top_pathogens) { arm_names <- c("individual_only", "summary_only", "federated") - arm_labels <- c("A", "B", "C") + arm_labels <- c("I", "S", "F") # display labels matching ARM_LABELS keys primary_key <- names(dist_code_map[dist_code_map == cdf_dist_name & names(dist_code_map) == PRIMARY_DIST]) # Map PRIMARY_DIST label back to DIST_CODES key @@ -398,17 +409,17 @@ if (length(top_pathogens) == 0L) { density_df <- bind_rows(density_list) |> mutate( - arm = factor(arm, levels = c("A", "B", "C")), + arm = factor(arm, levels = c("I", "S", "F")), pathogen = factor(pathogen, levels = top_pathogens) ) - # Also annotate each panel with the JS_CA and JS_CB values. + # Also annotate each panel with the JS_FI and JS_FS values. js_labels <- comparison_tbl |> filter(pathogen %in% top_pathogens, dist == PRIMARY_DIST) |> mutate( label = sprintf( - "JS[CA] == %.3f~~~~~JS[CB] == %.3f", - round(JS_CA, 3L), round(JS_CB, 3L) + "JS[FI] == %.3f~~~~~JS[FS] == %.3f", + round(JS_FI, 3L), round(JS_FS, 3L) ), pathogen = factor(pathogen, levels = top_pathogens) ) |> @@ -425,8 +436,8 @@ if (length(top_pathogens) == 0L) { left_join(y_pos, by = "pathogen") # Ribbon alpha by arm: federated slightly more prominent. - ribbon_alpha <- c("A" = 0.20, "B" = 0.20, "C" = 0.25) - line_lwd <- c("A" = 0.6, "B" = 0.6, "C" = 0.9) + ribbon_alpha <- c("I" = 0.20, "S" = 0.20, "F" = 0.25) + line_lwd <- c("I" = 0.6, "S" = 0.6, "F" = 0.9) fig3 <- ggplot(density_df, aes(x = x, group = arm, fill = arm, colour = arm)) + @@ -437,15 +448,15 @@ if (length(top_pathogens) == 0L) { aes(x = x_pos, y = y_pos, label = label), inherit.aes = FALSE, size = 2.6, hjust = 0, parse = TRUE, colour = "grey30") + - scale_fill_manual(values = ARM_COLOURS, labels = ARM_LABELS) + - scale_colour_manual(values = ARM_COLOURS, labels = ARM_LABELS) + + scale_fill_aaas(labels = ARM_LABELS) + + scale_colour_aaas(labels = ARM_LABELS) + scale_linewidth_manual(values = line_lwd, labels = ARM_LABELS) + facet_wrap(~ pathogen, scales = "free", ncol = 2L) + labs( x = "Incubation period (days)", y = "Density", caption = paste0( - "Top ", TOP_N_DENSITY, " pathogens by JS divergence (arm A vs federated). ", + "Top ", TOP_N_DENSITY, " pathogens by JS divergence (arm I vs federated). ", PRIMARY_DIST, " distribution. ", "Shaded bands: 95% posterior credible intervals." ) @@ -482,11 +493,11 @@ if (length(top_pathogens) == 0L) { gain_all_dists <- comparison_tbl |> mutate(pathogen = factor(pathogen, levels = pathogen_order)) |> select(pathogen, dist, - ratio_CA, ratio_CB, - JS_CA, JS_CB, - OVL_CA, OVL_CB) |> + ratio_FI, ratio_FS, + JS_FI, JS_FS, + OVL_FI, OVL_FS) |> pivot_longer( - cols = c(ratio_CA, ratio_CB, JS_CA, JS_CB, OVL_CA, OVL_CB), + cols = c(ratio_FI, ratio_FS, JS_FI, JS_FS, OVL_FI, OVL_FS), names_to = c("metric", "comparison"), names_sep = "_", values_to = "value" @@ -494,7 +505,7 @@ gain_all_dists <- comparison_tbl |> filter(!is.na(value)) |> mutate( comparison = factor(comparison, - levels = c("CA", "CB"), + levels = c("FI", "FS"), labels = c("vs (I) individual-level", "vs (S) summary-stats")), metric = factor(metric, @@ -521,10 +532,7 @@ fig_supp <- ggplot(gain_all_dists, linetype = "dashed", linewidth = 0.35, ) + - scale_colour_manual( - values = c("vs (I) individual-level" = ARM_COLOURS[["A"]], - "vs (S) summary-stats" = ARM_COLOURS[["B"]]) - ) + + scale_colour_aaas() + scale_shape_manual( values = c("vs (I) individual-level" = 19, "vs (S) summary-stats" = 17) @@ -555,19 +563,19 @@ message("Supplementary gain figure (all distributions) saved.") # A single three-column figure suitable for the main paper body. # # Panel A (wide) — Dumbbell forest plot (PRIMARY_DIST only) -# A grey segment connects the arm A and arm B point estimates; its length +# A grey segment connects the arm I and arm S point estimates; its length # encodes how much the two single-format arms disagree. Three -# point + CrI layers (A, B, C) are drawn on top. Arm C is rendered more +# point + CrI layers (I, S, F) are drawn on top. Arm F is rendered more # prominently (larger point, heavier error bar) because it is the focus. # A log-scaled x-axis accommodates the range across pathogens. # -# Panel B (narrow) — Interval ratio (arm_width / C_width) -# One dot per arm × pathogen pair (CA in blue, CB in orange). +# Panel B (narrow) — Interval ratio (arm_width / F_width) +# One dot per arm × pathogen pair (FI in blue, FS in orange). # Reference line at 1: points to the right mean the single-format arm is # wider (less precise) than the federated model. # # Panel C (narrow) — Jensen-Shannon divergence (bits) -# JS_CA and JS_CB on the same scale. Larger values signal a bigger +# JS_FI and JS_FS on the same scale. Larger values signal a bigger # distributional shift between that arm and the federated result. # Reference line at 0 (identical distributions). # @@ -621,7 +629,7 @@ best_dist_tbl <- data.frame( ) # Fall back to PRIMARY_DIST where the main-analysis best dist is unavailable -# in the ablation results (e.g. Gen. gamma was not fitted in arm B/C). +# in the ablation results (e.g. Gen. gamma was not fitted in arm S/F). available_dists <- unique(comparison_tbl$dist) needs_fallback <- is.na(best_dist_tbl$best_dist) | !(best_dist_tbl$best_dist %in% available_dists) @@ -639,25 +647,25 @@ wide_best <- comparison_tbl |> rowwise() |> mutate( # Median (P50) posterior predictive estimates + 95 % CrI bounds - A_med = .parse_cri(A_pred_median)[["med"]], - A_lo = .parse_cri(A_pred_median)[["lo"]], - A_hi = .parse_cri(A_pred_median)[["hi"]], - B_med = .parse_cri(B_pred_median)[["med"]], - B_lo = .parse_cri(B_pred_median)[["lo"]], - B_hi = .parse_cri(B_pred_median)[["hi"]], - C_med = .parse_cri(C_pred_median)[["med"]], - C_lo = .parse_cri(C_pred_median)[["lo"]], - C_hi = .parse_cri(C_pred_median)[["hi"]], + I_med = .parse_cri(I_pred_median)[["med"]], + I_lo = .parse_cri(I_pred_median)[["lo"]], + I_hi = .parse_cri(I_pred_median)[["hi"]], + S_med = .parse_cri(S_pred_median)[["med"]], + S_lo = .parse_cri(S_pred_median)[["lo"]], + S_hi = .parse_cri(S_pred_median)[["hi"]], + F_med = .parse_cri(F_pred_median)[["med"]], + F_lo = .parse_cri(F_pred_median)[["lo"]], + F_hi = .parse_cri(F_pred_median)[["hi"]], # 95th percentile (P95) posterior predictive estimates + 95 % CrI bounds - A_q95 = .parse_cri(A_pred_q95)[["med"]], - A_q95_lo = .parse_cri(A_pred_q95)[["lo"]], - A_q95_hi = .parse_cri(A_pred_q95)[["hi"]], - B_q95 = .parse_cri(B_pred_q95)[["med"]], - B_q95_lo = .parse_cri(B_pred_q95)[["lo"]], - B_q95_hi = .parse_cri(B_pred_q95)[["hi"]], - C_q95 = .parse_cri(C_pred_q95)[["med"]], - C_q95_lo = .parse_cri(C_pred_q95)[["lo"]], - C_q95_hi = .parse_cri(C_pred_q95)[["hi"]] + I_q95 = .parse_cri(I_pred_q95)[["med"]], + I_q95_lo = .parse_cri(I_pred_q95)[["lo"]], + I_q95_hi = .parse_cri(I_pred_q95)[["hi"]], + S_q95 = .parse_cri(S_pred_q95)[["med"]], + S_q95_lo = .parse_cri(S_pred_q95)[["lo"]], + S_q95_hi = .parse_cri(S_pred_q95)[["hi"]], + F_q95 = .parse_cri(F_pred_q95)[["med"]], + F_q95_lo = .parse_cri(F_pred_q95)[["lo"]], + F_q95_hi = .parse_cri(F_pred_q95)[["hi"]] ) |> ungroup() @@ -670,41 +678,41 @@ dist_label_lookup <- setNames( # ── 8c. Forest data — P50 and P95 facets ───────────────────────────────────── forest_p50 <- wide_best |> - select(pathogen, A_med, A_lo, A_hi, B_med, B_lo, B_hi, C_med, C_lo, C_hi) |> + select(pathogen, I_med, I_lo, I_hi, S_med, S_lo, S_hi, F_med, F_lo, F_hi) |> pivot_longer(-pathogen, names_to = c("arm", ".value"), names_pattern = "^(.)_(.*)") |> filter(!is.na(med)) |> - mutate(arm = factor(arm, levels = c("A", "B", "C")), metric = "Median (P50)") + mutate(arm = factor(arm, levels = c("I", "S", "F")), metric = "Median (P50)") forest_p95 <- wide_best |> select(pathogen, - A_med = A_q95, A_lo = A_q95_lo, A_hi = A_q95_hi, - B_med = B_q95, B_lo = B_q95_lo, B_hi = B_q95_hi, - C_med = C_q95, C_lo = C_q95_lo, C_hi = C_q95_hi) |> + I_med = I_q95, I_lo = I_q95_lo, I_hi = I_q95_hi, + S_med = S_q95, S_lo = S_q95_lo, S_hi = S_q95_hi, + F_med = F_q95, F_lo = F_q95_lo, F_hi = F_q95_hi) |> pivot_longer(-pathogen, names_to = c("arm", ".value"), names_pattern = "^(.)_(.*)") |> filter(!is.na(med)) |> - mutate(arm = factor(arm, levels = c("A", "B", "C")), metric = "95th percentile (P95)") + mutate(arm = factor(arm, levels = c("I", "S", "F")), metric = "95th percentile (P95)") -# Dumbbell backbone segments for both facets (A-to-B range, un-dodged). +# Dumbbell backbone segments for both facets (I-to-S range, un-dodged). dumbbell_segs <- bind_rows( - wide_best |> filter(!is.na(A_med), !is.na(B_med)) |> - transmute(pathogen, x = A_med, xend = B_med, metric = "Median (P50)"), - wide_best |> filter(!is.na(A_q95), !is.na(B_q95)) |> - transmute(pathogen, x = A_q95, xend = B_q95, metric = "95th percentile (P95)") + wide_best |> filter(!is.na(I_med), !is.na(S_med)) |> + transmute(pathogen, x = I_med, xend = S_med, metric = "Median (P50)"), + wide_best |> filter(!is.na(I_q95), !is.na(S_q95)) |> + transmute(pathogen, x = I_q95, xend = S_q95, metric = "95th percentile (P95)") ) |> mutate(metric = factor(metric, levels = c("Median (P50)", "95th percentile (P95)"))) # Dodging width and size scales (used by panels A and C). COMB_DODGE <- 0.5 -ARM_SIZES_C <- c("A" = 1.8, "B" = 1.8, "C" = 2.8) -ARM_EBW_C <- c("A" = 0.35, "B" = 0.35, "C" = 0.75) +ARM_SIZES_C <- c("I" = 1.8, "S" = 1.8, "F" = 2.8) +ARM_EBW_C <- c("I" = 0.35, "S" = 0.35, "F" = 0.75) # ── 8d. Extract mu0 CrI widths from Stan fits ───────────────────────────────── # # mu0 is the population-level location parameter (distribution's log scale). # Its 95% CrI width captures pure estimation uncertainty, independent of τ. -# Ratio arm / C > 1 means arm C has tighter knowledge of the population mean. +# Ratio arm / F > 1 means arm F has tighter knowledge of the population mean. # Returns a one-row data.frame with unname()d values to avoid the # "97.5%" row-name corruption that occurs when quantile output is passed @@ -723,13 +731,13 @@ ARM_EBW_C <- c("A" = 0.35, "B" = 0.35, "C" = 0.75) } DISPLAY_TO_INTERNAL <- setNames(names(MAIN_DIST_TO_DISPLAY), MAIN_DIST_TO_DISPLAY) -ARM_FIT_KEYS <- c("A" = "individual_only", "B" = "summary_only", - "C" = "federated") +ARM_FIT_KEYS <- c("I" = "individual_only", "S" = "summary_only", + "F" = "federated") mu0_tbl <- do.call(rbind, lapply(best_dist_tbl$pathogen, function(p) { bd_int <- DISPLAY_TO_INTERNAL[[ best_dist_tbl$best_dist[best_dist_tbl$pathogen == p] ]] if (is.na(bd_int)) return(NULL) - do.call(rbind, lapply(c("A", "B", "C"), function(arm) { + do.call(rbind, lapply(c("I", "S", "F"), function(arm) { slot <- ablation_fits[[p]][[ ARM_FIT_KEYS[[arm]] ]][[ bd_int ]] cbind( data.frame(pathogen = p, arm = arm, stringsAsFactors = FALSE), @@ -738,25 +746,25 @@ mu0_tbl <- do.call(rbind, lapply(best_dist_tbl$pathogen, function(p) { })) })) |> filter(!is.na(mu0_width)) |> - mutate(arm = factor(arm, levels = c("A", "B", "C")), + mutate(arm = factor(arm, levels = c("I", "S", "F")), pathogen = factor(pathogen, levels = levels(wide_best$pathogen))) -# mu0 CrI width ratio: arm / C (> 1 → arm C is more precise about μ). -# mu0 CrI width ratio: arm / C (> 1 → arm C is more precise about μ). +# mu0 CrI width ratio: arm / F (> 1 → arm F is more precise about μ). +# mu0 CrI width ratio: arm / F (> 1 → arm F is more precise about μ). # Built via left_join to avoid pivot_wider column-existence issues. -mu0_c <- mu0_tbl |> - filter(as.character(arm) == "C") |> - select(pathogen, c_width = mu0_width) +mu0_f <- mu0_tbl |> + filter(as.character(arm) == "F") |> + select(pathogen, f_width = mu0_width) mu0_ratio <- mu0_tbl |> - filter(as.character(arm) %in% c("A", "B")) |> - left_join(mu0_c, by = "pathogen") |> - filter(!is.na(mu0_width), !is.na(c_width)) |> + filter(as.character(arm) %in% c("I", "S")) |> + left_join(mu0_f, by = "pathogen") |> + filter(!is.na(mu0_width), !is.na(f_width)) |> mutate( - ratio = mu0_width / c_width, + ratio = mu0_width / f_width, comparison = factor(as.character(arm), - levels = c("A", "B"), - labels = ARM_LABELS[c("A", "B")]) + levels = c("I", "S"), + labels = ARM_LABELS[c("I", "S")]) ) |> select(pathogen, comparison, ratio) @@ -768,23 +776,23 @@ mu0_ratio <- mu0_tbl |> tau_long <- wide_best |> rowwise() |> mutate( - A_tau_med = .parse_cri(A_tau)[["med"]], - A_tau_lo = .parse_cri(A_tau)[["lo"]], - A_tau_hi = .parse_cri(A_tau)[["hi"]], - B_tau_med = .parse_cri(B_tau)[["med"]], - B_tau_lo = .parse_cri(B_tau)[["lo"]], - B_tau_hi = .parse_cri(B_tau)[["hi"]], - C_tau_med = .parse_cri(C_tau)[["med"]], - C_tau_lo = .parse_cri(C_tau)[["lo"]], - C_tau_hi = .parse_cri(C_tau)[["hi"]] + I_tau_med = .parse_cri(I_tau)[["med"]], + I_tau_lo = .parse_cri(I_tau)[["lo"]], + I_tau_hi = .parse_cri(I_tau)[["hi"]], + S_tau_med = .parse_cri(S_tau)[["med"]], + S_tau_lo = .parse_cri(S_tau)[["lo"]], + S_tau_hi = .parse_cri(S_tau)[["hi"]], + F_tau_med = .parse_cri(F_tau)[["med"]], + F_tau_lo = .parse_cri(F_tau)[["lo"]], + F_tau_hi = .parse_cri(F_tau)[["hi"]] ) |> ungroup() |> - select(pathogen, matches("^[ABC]_tau_(med|lo|hi)$")) |> + select(pathogen, matches("^[ISF]_tau_(med|lo|hi)$")) |> pivot_longer(-pathogen, names_to = c("arm", ".value"), names_pattern = "^(.)_tau_(.*)$") |> filter(!is.na(med)) |> - mutate(arm = factor(arm, levels = c("A", "B", "C"))) + mutate(arm = factor(arm, levels = c("I", "S", "F"))) # ── 8e-filter. Restrict all panels to pathogens with τ for all three arms ───── # @@ -799,7 +807,7 @@ tau_complete <- tau_long |> pull(pathogen) |> as.character() -# Preserve original C-arm ordering (federated median, descending). +# Preserve original F-arm ordering (federated median, descending). tau_levels <- intersect(levels(wide_best$pathogen), tau_complete) .restrict_pathogens <- function(df) { @@ -812,6 +820,7 @@ forest_p50 <- .restrict_pathogens(forest_p50) forest_p95 <- .restrict_pathogens(forest_p95) dumbbell_segs <- .restrict_pathogens(dumbbell_segs) tau_long <- .restrict_pathogens(tau_long) +mu0_ratio_all <- mu0_ratio # unrestricted — used for the top-panel overview mu0_ratio <- .restrict_pathogens(mu0_ratio) wide_best <- .restrict_pathogens(wide_best) # feeds pred_ratio in 8g @@ -820,17 +829,15 @@ n_tau_pathogens <- length(tau_levels) # ── 8f. Unified arm labels for legend merging ───────────────────────────────── # -# Convert the "A"/"B"/"C" factor to full label text in all forest/tau/mu0 data. +# Convert the "I"/"S"/"F" factor to full label text in all forest/tau/mu0 data. # All panels then use the same colour/shape scale keyed by label text, so # patchwork's guides = "collect" produces a single merged legend. -ARM_COLOURS_FULL <- setNames(ARM_COLOURS, ARM_LABELS[names(ARM_COLOURS)]) ARM_SHAPES_FULL <- setNames(ARM_SHAPES, ARM_LABELS[names(ARM_SHAPES)]) ARM_SIZES_C_FULL <- setNames(ARM_SIZES_C, ARM_LABELS[names(ARM_SIZES_C)]) ARM_EBW_C_FULL <- setNames(ARM_EBW_C, ARM_LABELS[names(ARM_EBW_C)]) -COMP_COLOURS <- ARM_COLOURS_FULL[ARM_LABELS[c("A", "B")]] -COMP_SHAPES <- ARM_SHAPES_FULL[ ARM_LABELS[c("A", "B")]] +COMP_SHAPES <- ARM_SHAPES_FULL[ ARM_LABELS[c("I", "S")]] .relabel_arm <- function(arm_fac) { factor(ARM_LABELS[as.character(arm_fac)], levels = unname(ARM_LABELS)) @@ -874,7 +881,7 @@ pA <- ggplot(forest_both, scale_x_log10(breaks = c(1, 2, 5, 10, 20, 50, 100), labels = c("1", "2", "5", "10", "20", "50", "100")) + scale_y_discrete(labels = dist_label_lookup) + - scale_colour_manual(values = ARM_COLOURS_FULL) + + scale_colour_aaas() + scale_shape_manual( values = ARM_SHAPES_FULL) + scale_size_manual( values = ARM_SIZES_C_FULL) + scale_linewidth_manual(values = ARM_EBW_C_FULL) + @@ -887,21 +894,21 @@ pA <- ggplot(forest_both, shape = "none") # Panel B: μ₀ CrI width ratio — pure information gain. -# Values > 1 indicate arm C has tighter posterior for the population mean. +# Values > 1 indicate arm F has tighter posterior for the population mean. pB <- ggplot(mu0_ratio, aes(x = ratio, y = pathogen, colour = comparison, shape = comparison)) + geom_vline(xintercept = 1, linetype = "dashed", colour = "grey50", linewidth = 0.4) + geom_point(size = 2.2) + - scale_colour_manual(values = COMP_COLOURS) + + scale_colour_aaas() + scale_shape_manual( values = COMP_SHAPES) + - labs(x = expression(mu[0]~"CrI ratio (arm / C)"), + labs(x = expression(mu[0]~"CrI ratio (arm / F)"), subtitle = "Information gain") + y_shared # Panel C: τ per arm — between-study heterogeneity. -# Arm C may detect larger τ when cross-data-type contrast reveals more +# Arm F may detect larger τ when cross-data-type contrast reveals more # between-study variation; its τ posterior is also better estimated # (narrower CrI) due to more studies. pC <- ggplot(tau_long, @@ -910,7 +917,7 @@ pC <- ggplot(tau_long, height = 0, linewidth = 0.35, position = position_dodge(width = COMB_DODGE)) + geom_point(size = 2.0, position = position_dodge(width = COMB_DODGE)) + - scale_colour_manual(values = ARM_COLOURS_FULL) + + scale_colour_aaas() + scale_shape_manual( values = ARM_SHAPES_FULL) + labs(x = expression(tau~"(heterogeneity SD)"), subtitle = "Between-study heterogeneity") + @@ -918,14 +925,14 @@ pC <- ggplot(tau_long, # Panel D: predictive CrI ratio — the combined (confounded) signal. pred_ratio <- wide_best |> - select(pathogen, ratio_CA, ratio_CB) |> - pivot_longer(cols = c(ratio_CA, ratio_CB), + select(pathogen, ratio_FI, ratio_FS) |> + pivot_longer(cols = c(ratio_FI, ratio_FS), names_to = "comparison", values_to = "ratio") |> filter(!is.na(ratio)) |> mutate(comparison = factor(comparison, - levels = c("ratio_CA", "ratio_CB"), - labels = ARM_LABELS[c("A", "B")])) + levels = c("ratio_FI", "ratio_FS"), + labels = ARM_LABELS[c("I", "S")])) pD <- ggplot(pred_ratio, aes(x = ratio, y = pathogen, @@ -933,21 +940,48 @@ pD <- ggplot(pred_ratio, geom_vline(xintercept = 1, linetype = "dashed", colour = "grey50", linewidth = 0.4) + geom_point(size = 2.2) + - scale_colour_manual(values = COMP_COLOURS) + + scale_colour_aaas() + scale_shape_manual( values = COMP_SHAPES) + - labs(x = "Predictive CrI ratio (arm / C)", + labs(x = "Predictive CrI ratio (arm / F)", subtitle = "Combined effect") + y_shared -# ── 8h. Assemble and save ───────────────────────────────────────────────────── +# ── 8h. Top panel: information gain across all pathogens ───────────────────── +# +# Pathogens on the x-axis (ordered by federated pred_median), mu0 CrI width +# ratio on the y-axis. Each pathogen has two points: arm I vs F and arm S vs F. +# A horizontal reference line at 1 marks equal precision. + +pTop <- ggplot(mu0_ratio_all, + aes(x = pathogen, y = ratio, + colour = comparison, shape = comparison)) + + geom_hline(yintercept = 1, linetype = "dashed", + colour = "grey50", linewidth = 0.4) + + geom_point(size = 2.2) + + scale_x_discrete(drop = FALSE) + + scale_colour_aaas() + + scale_shape_manual(values = COMP_SHAPES) + + labs(y = expression(mu[0]~"CrI ratio (arm / F)"), + subtitle = "Information gain (values > 1 indicate federated arm is more precise about population mean)") + + theme_ablation() + + theme( + axis.title.y = element_text(), + axis.title.x = element_blank(), + axis.text.x = element_text(angle = 45, hjust = 1, size = 8) + ) + +# ── 8i. Assemble and save ───────────────────────────────────────────────────── + +bottom_row <- (pA | pB | pC | pD) + + plot_layout(widths = c(2, 1, 1, 1)) -fig4 <- (pA | pB | pC | pD) + - plot_layout(widths = c(2, 1, 1, 1), guides = "collect") + +fig4 <- (pTop / bottom_row) + + plot_layout(heights = c(1, 2), guides = "collect") + plot_annotation(tag_levels = "A") & theme(legend.position = "bottom", plot.tag = element_text(face = "bold", size = 10)) -fig4_height <- max(4, n_tau_pathogens * FIG_HEIGHT_ROW + 1.8) +fig4_height <- max(7, n_tau_pathogens * FIG_HEIGHT_ROW + 5) ggsave(file.path(OUTPUT_DIR, "fig_ablation_combined.pdf"), fig4, width = FIG_WIDTH_WIDE * 1.35, height = fig4_height, device = "pdf") From 4221d0834ea42b4f3a268b26ed6754cf4d4cc436 Mon Sep 17 00:00:00 2001 From: cm401 Date: Mon, 25 May 2026 16:56:28 +0100 Subject: [PATCH 15/18] small updates to code --- R/table_utils.R | 299 ++++++++++++++++++++++++++ analysis/make_supplementary_figures.R | 219 +++++++++++-------- analysis/plot_data_format_ablation.R | 98 ++++++--- analysis/run_meta_analysis.R | 23 +- 4 files changed, 506 insertions(+), 133 deletions(-) diff --git a/R/table_utils.R b/R/table_utils.R index 2f5c326..7a0b52c 100644 --- a/R/table_utils.R +++ b/R/table_utils.R @@ -740,3 +740,302 @@ generate_data_table <- function( paste(lines, collapse = "\n") } + + +# ── Simulation study table helpers ──────────────────────────────────────────── + +# Shared footnote for both simulation tables (no size command — the caller adds \footnotesize). +.sim_footnote <- paste0( + "\\textit{Note:} ", + "Coverage: empirical proportion of replicates where the true value fell inside the 95\\% posterior ", + "credible interval (nominal target: 95\\%). ", + "Predictive bias: median signed error (posterior median $-$ true value) across replicates. ", + "IQD: integrated quadratic distance between true and estimated predictive CDFs. ", + "WIS: weighted interval score (lower is better for IQD and WIS). ", + "Gen.\\,Gamma scenarios with fewer than 10 identifiable replicates excluded prior to aggregation." +) + +# Shared two-row column header builder. +# n_id_cols: number of identifying columns before the metric block. +# id_header: LaTeX for the id-column cells in the second header row. +.sim_col_headers <- function(n_id_cols, id_header, total_cols) { + span <- function(n, label) sprintf("\\multicolumn{%d}{c}{%s}", n, label) + + row1 <- paste0( + "\\rowcolor[HTML]{F0E68C}", + paste(rep("", n_id_cols), collapse = " & "), " & ", + span(3L, "Parameter coverage (\\%)"), " & ", + span(2L, "Predictive coverage (\\%)"), " & ", + span(2L, "Predictive bias (days)"), " & ", + span(2L, "Scoring rules"), + " \\\\" + ) + + m <- n_id_cols + 1L + cmidrule <- paste0( + sprintf("\\cmidrule(lr){%d-%d}", m, m + 2L), + sprintf("\\cmidrule(lr){%d-%d}", m + 3L, m + 4L), + sprintf("\\cmidrule(lr){%d-%d}", m + 5L, m + 6L), + sprintf("\\cmidrule(lr){%d-%d}", m + 7L, m + 8L) + ) + + row2 <- paste0( + "\\rowcolor[HTML]{F0E68C}", + id_header, + " & $\\mu_0$ & $\\tau$ & $\\phi$", + " & P50 & P95", + " & P50 & P95", + " & IQD & WIS", + " \\\\" + ) + + c(row1, cmidrule, row2) +} + +# Format the metric portion of one data row (N + 9 metric cells). +.sim_metric_cells <- function(n, cov_mu0, cov_tau, cov_phi, + cov_p50, cov_p95, bias_p50, bias_p95, iqd, wis) { + paste( + as.character(n), + sprintf("%.1f", cov_mu0 * 100), + sprintf("%.1f", cov_tau * 100), + sprintf("%.1f", cov_phi * 100), + sprintf("%.1f", cov_p50 * 100), + sprintf("%.1f", cov_p95 * 100), + sprintf("%+.2f", bias_p50), + sprintf("%+.2f", bias_p95), + sprintf("%.3f", iqd), + sprintf("%.2f", wis), + sep = " & " + ) +} + +# Group-header row spanning all columns. +.sim_group_header <- function(label, total_cols) { + paste0( + "\\rowcolor{gray!15}\\multicolumn{", total_cols, + "}{@{}l}{\\textbf{", label, "}} \\\\" + ) +} + +# Assemble the longtable from pre-built body lines and header components. +# The footnote is placed as a paragraph AFTER \end{longtable} rather than +# inside the table: a p{wide} multicolumn in \endlastfoot forces LaTeX to +# expand the last column to meet the stated width, pushing WIS far to the right. +.sim_longtable <- function(body_lines, col_spec, headers, caption, label, + total_cols, fn_text) { + cont_head <- c( + sprintf("\\multicolumn{%d}{l}{\\small\\textit{continued from previous page}} \\\\", + total_cols), + "\\toprule", + headers, + "\\midrule", + "\\endhead" + ) + + c( + "% Requires \\usepackage{booktabs}, \\usepackage{longtable},", + "% \\usepackage[table]{xcolor}, \\usepackage{pdflscape} in preamble.", + "\\begin{landscape}", + "\\small", + paste0("\\begin{longtable}{", col_spec, "}"), + paste0("\\caption{", caption, "}\\label{", label, "} \\\\"), + "\\toprule", + headers, + "\\midrule", + "\\endfirsthead", + cont_head, + sprintf("\\multicolumn{%d}{r}{\\small\\textit{continued on next page}} \\\\", + total_cols), + "\\endfoot", + "\\bottomrule", + "\\endlastfoot", + "%", + body_lines, + "\\end{longtable}", + "% Note placed outside longtable to avoid forcing column-width expansion.", + paste0("\\par\\smallskip\\footnotesize ", fn_text), + "\\end{landscape}" + ) +} + + +#' Generate a LaTeX simulation performance table by distribution and summary type +#' +#' Produces a landscape \pkg{longtable} summarising simulation study performance +#' aggregated over numbers of datasets and within-study sample sizes, with one +#' row per distribution-family/summary-type combination, grouped by distribution. +#' +#' @param summary_res Output of [create_results_summary()]. +#' @param caption LaTeX \code{\\caption\{\}} string. +#' @param label LaTeX \code{\\label\{\}} string. +#' @return A character string with the complete LaTeX \code{longtable} environment. +#' @export +generate_simulation_table1 <- function( + summary_res, + caption = paste0( + "Simulation study performance by distribution family and summary type. ", + "Coverage values are empirical 95\\% credible interval coverage (\\%; nominal target: 95\\%). ", + "Predictive bias is the median signed error (posterior median $-$ true value) across replicates. ", + "IQD: integrated quadratic distance; WIS: weighted interval score (lower is better). ", + "Metrics are averaged over numbers of datasets and within-study sample sizes." + ), + label = "tab:sim_dist_summary" +) { + dist_levels <- c("lognormal", "gamma", "weibull", "burr12", "gengamma") + dist_labels <- c("Log-normal", "Gamma", "Weibull", "Burr XII", "Gen.~Gamma") + st_levels <- 1:5 + st_labels <- c("Median + Range", "Median + IQR", "Mean + SD", + "Freq.~table", "Mixed") + + t1 <- summary_res %>% + mutate( + dist_label = factor(dist_type, levels = dist_levels, labels = dist_labels), + st_label = factor(summary_type, levels = st_levels, labels = st_labels) + ) %>% + group_by(dist_label, st_label) %>% + summarise( + n_scen = n(), + cov_mu0 = mean(coverage_mu0, na.rm = TRUE), + cov_tau = mean(coverage_tau, na.rm = TRUE), + cov_phi = mean(coverage_phi, na.rm = TRUE), + cov_p50 = mean(coverage_pred_median, na.rm = TRUE), + cov_p95 = mean(coverage_pred_q95, na.rm = TRUE), + bias_p50 = median(bias_pred_median, na.rm = TRUE), + bias_p95 = median(bias_pred_q95, na.rm = TRUE), + iqd = mean(mean_iqd, na.rm = TRUE), + wis = mean(mean_wis, na.rm = TRUE), + .groups = "drop" + ) %>% + arrange(dist_label, st_label) + + total_cols <- 11L + body_lines <- character(0L) + + for (i in seq_along(dist_labels)) { + d <- dist_labels[i] + rows <- t1[t1$dist_label == d, ] + if (nrow(rows) == 0L) next + + if (i > 1L) body_lines <- c(body_lines, "\\midrule") + body_lines <- c(body_lines, .sim_group_header(d, total_cols)) + + for (j in seq_len(nrow(rows))) { + r <- rows[j, ] + cells <- paste( + as.character(r$st_label), + .sim_metric_cells(r$n_scen, r$cov_mu0, r$cov_tau, r$cov_phi, + r$cov_p50, r$cov_p95, r$bias_p50, r$bias_p95, + r$iqd, r$wis), + sep = " & " + ) + body_lines <- c(body_lines, paste0("\\quad ", cells, " \\\\")) + } + } + + headers <- .sim_col_headers(2L, "Summary type & $N$", total_cols) + col_spec <- "@{} l r r r r r r r r r r @{}" + + paste(.sim_longtable(body_lines, col_spec, headers, caption, label, + total_cols, .sim_footnote), + collapse = "\n") +} + + +#' Generate a LaTeX simulation performance table by data availability +#' +#' Produces a landscape \pkg{longtable} summarising simulation study performance +#' by summary type, number of datasets, and within-study sample size, pooled +#' over distribution families, with rows grouped by summary type. +#' +#' @param summary_res Output of [create_results_summary()]. +#' @param caption LaTeX \code{\\caption\{\}} string. +#' @param label LaTeX \code{\\label\{\}} string. +#' @return A character string with the complete LaTeX \code{longtable} environment. +#' @export +generate_simulation_table2 <- function( + summary_res, + caption = paste0( + "Simulation study performance by summary type, number of datasets, and within-study sample size ", + "(pooled over distribution families). ", + "Coverage values are empirical 95\\% credible interval coverage (\\%; nominal target: 95\\%). ", + "Predictive bias is the median signed error (posterior median $-$ true value) across replicates. ", + "IQD: integrated quadratic distance; WIS: weighted interval score (lower is better)." + ), + label = "tab:sim_data_availability" +) { + st_levels <- 1:5 + st_labels <- c("Median + Range", "Median + IQR", "Mean + SD", + "Freq.~table", "Mixed") + + t2 <- summary_res %>% + mutate( + st_label = factor(summary_type, levels = st_levels, labels = st_labels), + n_datasets_bucket = factor( + case_when( + n_datasets < 10 ~ "$<$10", + n_datasets < 20 ~ "$<$20", + n_datasets < 30 ~ "$<$30", + n_datasets >= 30 ~ "$\\geq$30" + ), + levels = c("$<$10", "$<$20", "$<$30", "$\\geq$30") + ), + n_obs_bucket = factor( + case_when( + n_obs == 5 ~ "5", + n_obs == 10 ~ "10", + n_obs == 20 ~ "20", + n_obs > 25 ~ "$\\geq$25" + ), + levels = c("5", "10", "20", "$\\geq$25") + ) + ) %>% + group_by(st_label, n_datasets_bucket, n_obs_bucket) %>% + summarise( + n_scen = n(), + cov_mu0 = mean(coverage_mu0, na.rm = TRUE), + cov_tau = mean(coverage_tau, na.rm = TRUE), + cov_phi = mean(coverage_phi, na.rm = TRUE), + cov_p50 = mean(coverage_pred_median, na.rm = TRUE), + cov_p95 = mean(coverage_pred_q95, na.rm = TRUE), + bias_p50 = median(bias_pred_median, na.rm = TRUE), + bias_p95 = median(bias_pred_q95, na.rm = TRUE), + iqd = mean(mean_iqd, na.rm = TRUE), + wis = mean(mean_wis, na.rm = TRUE), + .groups = "drop" + ) %>% + arrange(st_label, n_datasets_bucket, n_obs_bucket) + + # summary type is a group header — data rows have 3 id cols (N datasets, n obs, N) + total_cols <- 12L + body_lines <- character(0L) + + for (i in seq_along(st_labels)) { + st <- st_labels[i] + rows <- t2[t2$st_label == st, ] + if (nrow(rows) == 0L) next + + if (i > 1L) body_lines <- c(body_lines, "\\midrule") + body_lines <- c(body_lines, .sim_group_header(st, total_cols)) + + for (j in seq_len(nrow(rows))) { + r <- rows[j, ] + cells <- paste( + as.character(r$n_datasets_bucket), + as.character(r$n_obs_bucket), + .sim_metric_cells(r$n_scen, r$cov_mu0, r$cov_tau, r$cov_phi, + r$cov_p50, r$cov_p95, r$bias_p50, r$bias_p95, + r$iqd, r$wis), + sep = " & " + ) + body_lines <- c(body_lines, paste0("\\quad ", cells, " \\\\")) + } + } + + headers <- .sim_col_headers(3L, "$N$ datasets & $n$ obs & $N$", total_cols) + col_spec <- "@{} l r r r r r r r r r r r @{}" + + paste(.sim_longtable(body_lines, col_spec, headers, caption, label, + total_cols, .sim_footnote), + collapse = "\n") +} diff --git a/analysis/make_supplementary_figures.R b/analysis/make_supplementary_figures.R index aaf87f8..d2b4391 100644 --- a/analysis/make_supplementary_figures.R +++ b/analysis/make_supplementary_figures.R @@ -466,6 +466,13 @@ library(dplyr) df$included <- factor(ifelse(df$included, "Included", "Excluded"), levels = c("Included", "Excluded")) + # Scale label font size down for crowded panels + n_rows <- nrow(df) + y_text_scale <- if (n_rows > 30) 0.45 + else if (n_rows > 20) 0.55 + else if (n_rows > 12) 0.65 + else 0.78 + ggplot2::ggplot(df, ggplot2::aes(y = label, colour = type, alpha = included)) + ggplot2::geom_segment( ggplot2::aes(x = lower, xend = upper, yend = label, linetype = type), @@ -502,7 +509,7 @@ library(dplyr) ggplot2::theme( panel.grid.minor = ggplot2::element_blank(), panel.grid.major.y = ggplot2::element_blank(), - axis.text.y = ggplot2::element_text(size = base_size * 0.75), + axis.text.y = ggplot2::element_text(size = base_size * y_text_scale), legend.position = "right", plot.title = ggplot2::element_text( size = base_size, face = "bold", @@ -813,13 +820,11 @@ build_supp_figure <- function(pathogen, error = function(e) { message(" [WARN] Panel F: ", conditionMessage(e)); NULL } ) - # ── Assemble with patchwork ────────────────────────────────────────────────── - # Layout: [A | B] on the top row (B omitted if NULL → A takes full width) - # [C] full-width (omitted if NULL) - # [D] full-width (omitted if NULL) - # [F] full-width (omitted if NULL) + # ── Assemble: two separate figures ────────────────────────────────────────── + # Figure 1 (CDF): panels A, B, C — model comparison and subgroup CDFs + # Figure 2 (Data): panels D, F — raw summary statistics and frequency tables # - # Heights are set relative to the number of datasets for D/F. + # Splitting keeps each PDF at a manageable height for \includegraphics in LaTeX. n_summary_ds <- sum(vapply(datasets_all, function(d) .detect_type(d) %in% c("A","B","C"), @@ -828,50 +833,100 @@ build_supp_figure <- function(pathogen, function(d) .detect_type(d) %in% c("D","E"), logical(1L))) - h_D <- if (!is.null(panel_D)) max(4, min(18, n_summary_ds * 0.70)) else 0 - h_F <- if (!is.null(panel_F)) max(4, min(16, n_freq_ds * 3.0)) else 0 - - # Row 1: Panel A (and optionally B) - if (!is.null(panel_A) && !is.null(panel_B)) { - row1 <- panel_A + panel_B + patchwork::plot_layout(ncol = 2L) - } else if (!is.null(panel_A)) { - row1 <- panel_A + # ── CDF figure ─────────────────────────────────────────────────────────────── + row1 <- if (!is.null(panel_A) && !is.null(panel_B)) { + panel_A + panel_B + patchwork::plot_layout(ncol = 2L) } else { - row1 <- NULL + panel_A } - # Collect non-NULL rows in order - rows <- Filter(Negate(is.null), list(row1, panel_C, panel_D, panel_F)) - heights <- c( + cdf_rows <- Filter(Negate(is.null), list(row1, panel_C)) + cdf_heights <- c( if (!is.null(row1)) 6 else NULL, - if (!is.null(panel_C)) 5.5 else NULL, - if (!is.null(panel_D)) h_D else NULL, - if (!is.null(panel_F)) h_F else NULL + if (!is.null(panel_C)) 5.5 else NULL ) - if (length(rows) == 0L) return(NULL) + fig_cdf <- if (length(cdf_rows) > 0L) { + patchwork::wrap_plots(cdf_rows, ncol = 1L, heights = cdf_heights) + } else NULL - fig <- patchwork::wrap_plots(rows, ncol = 1L, - heights = heights) + - patchwork::plot_annotation( - title = label, - theme = ggplot2::theme( - plot.title = ggplot2::element_text( - face = "bold", - size = base_size + 4L, - hjust = 0.5, - margin = ggplot2::margin(b = 6) - ) - ) - ) + # ── Data figure ────────────────────────────────────────────────────────────── + # Height scales with dataset count; no hard cap so crowded panels are legible. + h_D <- if (!is.null(panel_D)) max(4, n_summary_ds * 0.65) else 0 + h_F <- if (!is.null(panel_F)) max(4, n_freq_ds * 3.0) else 0 - fig + data_rows <- Filter(Negate(is.null), list(panel_D, panel_F)) + data_heights <- c( + if (!is.null(panel_D)) h_D else NULL, + if (!is.null(panel_F)) h_F else NULL + ) + + fig_data <- if (length(data_rows) > 0L) { + patchwork::wrap_plots(data_rows, ncol = 1L, heights = data_heights) + } else NULL + + list( + cdf = fig_cdf, + data = fig_data, + has_B = !is.null(panel_B), + has_C = !is.null(panel_C), + n_summary_ds = n_summary_ds, + n_freq_ds = n_freq_ds + ) } # ============================================================================= -# Main loop — generate and save one PDF per pathogen +# Main loop — generate and save two PDFs per pathogen # ============================================================================= +# +# Output files: +# _supplementary_cdf.pdf — panels A, B, C (model fits / CDFs) +# _supplementary_data.pdf — panels D, F (raw data) +# +# LaTeX inclusion (use two figure environments per pathogen): +# +# \begin{figure}[htbp] +# \centering +# \includegraphics[width=\linewidth]{figures/SI/_supplementary_cdf.pdf} +# \caption{ +# \textbf{Supplementary Figure~SX. Incubation period model fits for +# [PATHOGEN].} +# (\textbf{A})~Posterior predictive cumulative distribution functions +# (CDFs) for all converged parametric distributions fitted to the +# filtered dataset. Ribbons indicate 95\% credible intervals; the +# best-fitting distribution is annotated. +# (\textbf{B})~Comparison of CDFs fitted to all data vs.\ the filtered +# dataset for the best-fitting distribution; shown only when filtering +# removed at least one dataset. +# (\textbf{C})~Subgroup CDFs for the best-fitting distribution; the +# overall estimate (shaded ribbon) is shown alongside subgroup-specific +# estimates; shown only when subgroup analyses were performed. +# } +# \label{fig:_cdf} +# \end{figure} +# +# \begin{figure}[htbp] +# \centering +# \includegraphics[width=\linewidth]{figures/SI/_supplementary_data.pdf} +# \caption{ +# \textbf{Supplementary Figure~SX (continued). Raw incubation period +# data for [PATHOGEN].} +# (\textbf{D})~Raw summary statistics reported in source datasets. +# Points indicate the central estimate (circle: median; square: mean) +# and lines the reported uncertainty interval (solid: range; +# dashed: IQR; dotted: SD). Faded entries were excluded by the +# data-quality filter. +# (\textbf{F})~Empirical CDFs for datasets reporting frequency tables +# or interval-censored observations, overlaid with the posterior +# predictive CDF for the best-fitting distribution (coloured ribbon +# and line). For interval-censored data, the solid step line is the +# conservative ECDF at interval upper bounds; the dashed step line is +# the ECDF at interval lower bounds; the shaded band represents the +# uncertainty region. +# } +# \label{fig:_data} +# \end{figure} message("Loading main results (this may take a moment for large files)...") all_results <- readRDS(here::here("results", "main_results.rds")) @@ -895,7 +950,7 @@ for (pathogen in names(all_results)) { message("\n", strrep("-", 60)) message("Pathogen: ", pathogen) - fig <- tryCatch( + figs <- tryCatch( build_supp_figure( pathogen = pathogen, pathogen_results = all_results[[pathogen]], @@ -905,63 +960,47 @@ for (pathogen in names(all_results)) { ), error = function(e) { message(" [ERROR] ", conditionMessage(e)) - NULL + list(cdf = NULL, data = NULL, has_B = FALSE, has_C = FALSE, + n_summary_ds = 0L, n_freq_ds = 0L) } ) - if (is.null(fig)) { - message(" No figure produced — skipping.") - next + # ── CDF figure (panels A, B, C) ──────────────────────────────────────────── + if (!is.null(figs$cdf)) { + cdf_h <- 6.0 + (if (isTRUE(figs$has_C)) 5.5 else 0) + cdf_w <- if (isTRUE(figs$has_B)) 22 else 14 + cdf_path <- file.path(output_dir, paste0(pathogen, "_supplementary_cdf.pdf")) + ggplot2::ggsave( + filename = cdf_path, + plot = figs$cdf, + width = cdf_w, + height = max(cdf_h, 7), + units = "cm", + device = "pdf" + ) + message(" Saved: ", cdf_path) } - # Determine appropriate figure height from the assembled patchwork - # (each panel's height was set in cm; total ≈ sum of heights + title margin) - pdf_path <- file.path(output_dir, paste0(pathogen, "_supplementary.pdf")) - - # Collect panel info for this pathogen to set PDF dimensions - datasets_all <- tryCatch( - .get_datasets(all_results[[pathogen]][["all"]]), - error = function(e) list() - ) - if (is.null(datasets_all)) datasets_all <- list() - - # Width: wider when A and B are shown side-by-side (filtering removed data) - datasets_filt_names <- names(tryCatch( - .get_datasets(all_results[[pathogen]][["filtered"]]), - error = function(e) list() - )) - datasets_all_names <- names(datasets_all) - fig_width_cm <- if (!setequal(datasets_all_names, datasets_filt_names)) 22 else 18 - - has_C <- length(setdiff(names(all_results[[pathogen]]), - c("all", "filtered"))) > 0L - has_D <- any(vapply(datasets_all, - function(d) .detect_type(d) %in% c("A","B","C"), logical(1L))) - has_F <- any(vapply(datasets_all, - function(d) .detect_type(d) %in% c("D","E"), logical(1L))) + # ── Data figure (panels D, F) ─────────────────────────────────────────────── + if (!is.null(figs$data)) { + h_D <- if (figs$n_summary_ds > 0L) max(4, figs$n_summary_ds * 0.65) else 0 + h_F <- if (figs$n_freq_ds > 0L) max(4, figs$n_freq_ds * 3.0) else 0 + data_h <- max(h_D + h_F, 6) + data_path <- file.path(output_dir, paste0(pathogen, "_supplementary_data.pdf")) + ggplot2::ggsave( + filename = data_path, + plot = figs$data, + width = 18, + height = data_h, + units = "cm", + device = "pdf" + ) + message(" Saved: ", data_path) + } - n_summary_ds <- sum(vapply(datasets_all, - function(d) .detect_type(d) %in% c("A","B","C"), logical(1L))) - n_freq_ds <- sum(vapply(datasets_all, - function(d) .detect_type(d) %in% c("D","E"), logical(1L))) - - fig_height_cm <- 1.5 + # title margin - 6.0 + # row1 (A ± B) - (if (has_C) 5.5 else 0) + - (if (has_D) max(4, min(18, n_summary_ds * 0.70)) else 0) + - (if (has_F) max(4, min(16, n_freq_ds * 3.0)) else 0) - - fig_height_cm <- max(fig_height_cm, 10) - - ggplot2::ggsave( - filename = pdf_path, - plot = fig, - width = fig_width_cm, - height = fig_height_cm, - units = "cm", - device = "pdf" - ) - message(" Saved: ", pdf_path) + if (is.null(figs$cdf) && is.null(figs$data)) { + message(" No figures produced — skipping.") + } } message("\n", strrep("=", 60)) diff --git a/analysis/plot_data_format_ablation.R b/analysis/plot_data_format_ablation.R index 1c7ebd4..8b48e2d 100644 --- a/analysis/plot_data_format_ablation.R +++ b/analysis/plot_data_format_ablation.R @@ -705,7 +705,7 @@ dumbbell_segs <- bind_rows( # Dodging width and size scales (used by panels A and C). COMB_DODGE <- 0.5 -ARM_SIZES_C <- c("I" = 1.8, "S" = 1.8, "F" = 2.8) +ARM_SIZES_C <- c("I" = 3.0, "S" = 3.0, "F" = 4.2) ARM_EBW_C <- c("I" = 0.35, "S" = 0.35, "F" = 0.75) # ── 8d. Extract mu0 CrI widths from Stan fits ───────────────────────────────── @@ -821,6 +821,7 @@ forest_p95 <- .restrict_pathogens(forest_p95) dumbbell_segs <- .restrict_pathogens(dumbbell_segs) tau_long <- .restrict_pathogens(tau_long) mu0_ratio_all <- mu0_ratio # unrestricted — used for the top-panel overview +wide_best_all <- wide_best # unrestricted — used for P95 ratio top panel mu0_ratio <- .restrict_pathogens(mu0_ratio) wide_best <- .restrict_pathogens(wide_best) # feeds pred_ratio in 8g @@ -893,20 +894,6 @@ pA <- ggplot(forest_both, linewidth = "none", shape = "none") -# Panel B: μ₀ CrI width ratio — pure information gain. -# Values > 1 indicate arm F has tighter posterior for the population mean. -pB <- ggplot(mu0_ratio, - aes(x = ratio, y = pathogen, - colour = comparison, shape = comparison)) + - geom_vline(xintercept = 1, linetype = "dashed", - colour = "grey50", linewidth = 0.4) + - geom_point(size = 2.2) + - scale_colour_aaas() + - scale_shape_manual( values = COMP_SHAPES) + - labs(x = expression(mu[0]~"CrI ratio (arm / F)"), - subtitle = "Information gain") + - y_shared - # Panel C: τ per arm — between-study heterogeneity. # Arm F may detect larger τ when cross-data-type contrast reveals more # between-study variation; its τ posterior is also better estimated @@ -916,7 +903,7 @@ pC <- ggplot(tau_long, geom_errorbarh(aes(xmin = lo, xmax = hi), height = 0, linewidth = 0.35, position = position_dodge(width = COMB_DODGE)) + - geom_point(size = 2.0, position = position_dodge(width = COMB_DODGE)) + + geom_point(size = 3.5, position = position_dodge(width = COMB_DODGE)) + scale_colour_aaas() + scale_shape_manual( values = ARM_SHAPES_FULL) + labs(x = expression(tau~"(heterogeneity SD)"), @@ -934,35 +921,88 @@ pred_ratio <- wide_best |> levels = c("ratio_FI", "ratio_FS"), labels = ARM_LABELS[c("I", "S")])) +# Helper: predictive P95 CrI width ratio from a wide data frame that contains +# the parsed I/S/F q95 lo/hi columns. +.make_q95_ratio <- function(df) { + df |> + mutate( + ratio_FI = (I_q95_hi - I_q95_lo) / (F_q95_hi - F_q95_lo), + ratio_FS = (S_q95_hi - S_q95_lo) / (F_q95_hi - F_q95_lo) + ) |> + select(pathogen, ratio_FI, ratio_FS) |> + pivot_longer(cols = c(ratio_FI, ratio_FS), + names_to = "comparison", + values_to = "ratio") |> + filter(!is.na(ratio)) |> + mutate(comparison = factor(comparison, + levels = c("ratio_FI", "ratio_FS"), + labels = ARM_LABELS[c("I", "S")])) +} + +pred_q95_ratio_all <- .make_q95_ratio(wide_best_all) # all pathogens +pred_q95_ratio <- .make_q95_ratio(wide_best) # tau-complete pathogens + pD <- ggplot(pred_ratio, aes(x = ratio, y = pathogen, colour = comparison, shape = comparison)) + geom_vline(xintercept = 1, linetype = "dashed", colour = "grey50", linewidth = 0.4) + - geom_point(size = 2.2) + + geom_point(size = 3.5) + scale_colour_aaas() + scale_shape_manual( values = COMP_SHAPES) + - labs(x = "Predictive CrI ratio (arm / F)", - subtitle = "Combined effect") + + labs(x = "Predictive P50 CrI ratio (arm / F)", + subtitle = "Combined effect (P50)") + + y_shared + +pE2 <- ggplot(pred_q95_ratio, + aes(x = ratio, y = pathogen, + colour = comparison, shape = comparison)) + + geom_vline(xintercept = 1, linetype = "dashed", + colour = "grey50", linewidth = 0.4) + + geom_point(size = 3.5) + + scale_colour_aaas() + + scale_shape_manual(values = COMP_SHAPES) + + labs(x = "Predictive P95 CrI ratio (arm / F)", + subtitle = "Combined effect (P95)") + y_shared -# ── 8h. Top panel: information gain across all pathogens ───────────────────── +# ── 8h. Top panels: information gain across all pathogens ──────────────────── # -# Pathogens on the x-axis (ordered by federated pred_median), mu0 CrI width -# ratio on the y-axis. Each pathogen has two points: arm I vs F and arm S vs F. -# A horizontal reference line at 1 marks equal precision. +# Two full-width panels stacked vertically. Pathogens on the x-axis (ordered +# by federated pred_median), ratio on the y-axis. +# pTop — μ₀ CrI width ratio: pure information gain in the population mean. +# pTop2 — predictive P95 CrI width ratio: precision gain for the 95th %ile. pTop <- ggplot(mu0_ratio_all, aes(x = pathogen, y = ratio, colour = comparison, shape = comparison)) + geom_hline(yintercept = 1, linetype = "dashed", colour = "grey50", linewidth = 0.4) + - geom_point(size = 2.2) + + geom_point(size = 3.5) + scale_x_discrete(drop = FALSE) + scale_colour_aaas() + scale_shape_manual(values = COMP_SHAPES) + labs(y = expression(mu[0]~"CrI ratio (arm / F)"), - subtitle = "Information gain (values > 1 indicate federated arm is more precise about population mean)") + + subtitle = expression("Information gain: population mean "~mu[0]~ + " (all pathogens; values > 1 indicate federated arm is more precise)")) + + theme_ablation() + + theme( + axis.title.y = element_text(), + axis.title.x = element_blank(), + axis.text.x = element_text(angle = 45, hjust = 1, size = 8) + ) + +pTop2 <- ggplot(pred_q95_ratio_all, + aes(x = pathogen, y = ratio, + colour = comparison, shape = comparison)) + + geom_hline(yintercept = 1, linetype = "dashed", + colour = "grey50", linewidth = 0.4) + + geom_point(size = 3.5) + + scale_x_discrete(drop = FALSE) + + scale_colour_aaas() + + scale_shape_manual(values = COMP_SHAPES) + + labs(y = "Predictive P95 CrI ratio (arm / F)", + subtitle = "Precision gain for 95th percentile (all pathogens; values > 1 indicate federated arm is more precise)") + theme_ablation() + theme( axis.title.y = element_text(), @@ -972,16 +1012,16 @@ pTop <- ggplot(mu0_ratio_all, # ── 8i. Assemble and save ───────────────────────────────────────────────────── -bottom_row <- (pA | pB | pC | pD) + +bottom_row <- (pA | pC | pD | pE2) + plot_layout(widths = c(2, 1, 1, 1)) -fig4 <- (pTop / bottom_row) + - plot_layout(heights = c(1, 2), guides = "collect") + +fig4 <- (pTop / pTop2 / bottom_row) + + plot_layout(heights = c(2, 2, 5), guides = "collect") + plot_annotation(tag_levels = "A") & theme(legend.position = "bottom", plot.tag = element_text(face = "bold", size = 10)) -fig4_height <- max(7, n_tau_pathogens * FIG_HEIGHT_ROW + 5) +fig4_height <- max(12, n_tau_pathogens * FIG_HEIGHT_ROW + 9) ggsave(file.path(OUTPUT_DIR, "fig_ablation_combined.pdf"), fig4, width = FIG_WIDTH_WIDE * 1.35, height = fig4_height, device = "pdf") diff --git a/analysis/run_meta_analysis.R b/analysis/run_meta_analysis.R index 8efa316..8525e7b 100644 --- a/analysis/run_meta_analysis.R +++ b/analysis/run_meta_analysis.R @@ -196,17 +196,9 @@ dir.create(OUTDIR, showWarnings = FALSE, recursive = TRUE) } # Save a forest plot for one metamean result. +# meta 8.x uses grid graphics and auto-calculates the correct figure height +# when file= is supplied, so we let it open and close the PDF device itself. .save_forest <- function(m, df, label, path) { - n_studies <- nrow(df) - n_subgroups <- if (!is.null(m$subgroup)) - length(unique(na.omit(m$subgroup))) else 0L - # Allow extra lines for subgroup headers + within-subgroup pooled rows - n_rows <- n_studies + n_subgroups * 3L + 5L - fig_h <- max(7, 0.22 * n_rows + 3) - - pdf(path, width = 16, height = fig_h) - on.exit(dev.off(), add = TRUE) - forest( m, sortvar = TE, @@ -222,10 +214,13 @@ dir.create(OUTDIR, showWarnings = FALSE, recursive = TRUE) leftlabs = c("Study [dataset]", "N"), rightcols = c("effect", "ci"), rightlabs = c("Mean", "95% CI"), - print.subgroup.name = FALSE, - header.line = "both", - spacing = 0.8, - main = paste0(label, ": random-effects meta-analysis (MLN)") + print.subgroup.name = FALSE, + header.line = "both", + addrows.below.overall = 2L, + fontsize = 11, + colgap.forest.left = "4cm", + file = path, + width = 16 ) } From c97a897a46e692e5ca1ae4f3ca336a5c6b3416a2 Mon Sep 17 00:00:00 2001 From: cm401 Date: Mon, 25 May 2026 16:56:43 +0100 Subject: [PATCH 16/18] initial revision to create SI tables --- analysis/make_simulation_tables.R | 50 +++++++++++++++++++++++++++++++ 1 file changed, 50 insertions(+) create mode 100644 analysis/make_simulation_tables.R diff --git a/analysis/make_simulation_tables.R b/analysis/make_simulation_tables.R new file mode 100644 index 0000000..0c4694c --- /dev/null +++ b/analysis/make_simulation_tables.R @@ -0,0 +1,50 @@ +# analysis/make_simulation_tables.R +# ============================================================================= +# Generates LaTeX supplementary tables for the simulation study. +# +# Table S1: Performance by distribution family x summary type +# Table S2: Performance by data availability (N datasets x n obs x summary type) +# +# Usage: +# Rscript analysis/make_simulation_tables.R +# -- or -- +# source("analysis/make_simulation_tables.R") +# +# Outputs written to results/ +# simulation_table1.tex (Table Supplementary Information) +# simulation_table2.tex (Table Supplementary Information) +# ============================================================================= + +devtools::load_all(here::here(), quiet = TRUE) + +# ── Load simulation results ──────────────────────────────────────────────────── + +rds_path <- system.file("extdata", "simulation_results_all.rds", package = "ddsynth") +if (!nzchar(rds_path)) + rds_path <- file.path(here::here(), "vignettes", "simulation_results_all.rds") + +if (!file.exists(rds_path)) + stop("simulation_results_all.rds not found. Run analysis/new_simulation_study.R first.") + +res_out <- readRDS(rds_path) +cat(sprintf("Loaded %d rows across %d scenarios.\n", + nrow(res_out), length(unique(res_out$scenario_name)))) + +# ── Summarise ────────────────────────────────────────────────────────────────── + +summary_res <- create_results_summary(res_out) +cat(sprintf("Summary: %d scenario-level rows.\n", nrow(summary_res))) + +# ── Generate and write tables ───────────────────────────────────────────────── + +out_dir <- file.path(here::here(), "results") +if (!dir.exists(out_dir)) dir.create(out_dir, recursive = TRUE) + +tbl1 <- generate_simulation_table1(summary_res) +tbl2 <- generate_simulation_table2(summary_res) + +writeLines(tbl1, file.path(out_dir, "simulation_table1.tex")) +writeLines(tbl2, file.path(out_dir, "simulation_table2.tex")) + +cat("Written: results/simulation_table1.tex\n") +cat("Written: results/simulation_table2.tex\n") From fe6c0eb95bef5180ee8b06940a24f523face6e01 Mon Sep 17 00:00:00 2001 From: cm401 Date: Mon, 25 May 2026 17:16:47 +0100 Subject: [PATCH 17/18] minor fixes --- NAMESPACE | 4 ++++ R/ploting_utils.R | 4 ++-- analysis/plot_data_format_ablation.R | 2 +- 3 files changed, 7 insertions(+), 3 deletions(-) diff --git a/NAMESPACE b/NAMESPACE index 0c0e395..6ec00c8 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -42,11 +42,15 @@ export(generate_matched_moments_plot) export(generate_results_table) export(generate_results_table_split) export(generate_scenario_library) +export(generate_simulation_table1) +export(generate_simulation_table2) export(make_stan_init_fn) export(plot_main_figure) export(plot_main_figure_split) export(plot_simulation_study_figure) +export(plot_simulation_study_figure2) export(plot_simulation_study_sens_figure) +export(plot_simulation_study_si_figure) export(pre_inference_checks) export(prepare_stan_data_from_datasets) export(run_simulation_study_generalized) diff --git a/R/ploting_utils.R b/R/ploting_utils.R index 4b1a4f7..1d5602d 100644 --- a/R/ploting_utils.R +++ b/R/ploting_utils.R @@ -1733,7 +1733,7 @@ plot_simulation_study_figure2 <- function(summary_res, base_size = 9.5) { labs(tag = "B") # ── Assemble ─────────────────────────────────────────────────────────────── - (pA / pB) + + patchwork::wrap_plots(pA, pB, ncol = 1) + patchwork::plot_layout(heights = c(1, 1), guides = "collect") & theme(legend.position = "bottom", legend.box = "horizontal", @@ -1879,7 +1879,7 @@ plot_simulation_study_si_figure <- function(summary_res, base_size = 10) { labs(tag = "C") # ── Assemble ─────────────────────────────────────────────────────────────── - (pA / pB / pC) + + patchwork::wrap_plots(pA, pB, pC, ncol = 1) + patchwork::plot_layout(heights = c(1, 1, 0.55), guides = "collect") & theme(legend.position = "bottom", legend.box = "horizontal", diff --git a/analysis/plot_data_format_ablation.R b/analysis/plot_data_format_ablation.R index 8b48e2d..1a2ca03 100644 --- a/analysis/plot_data_format_ablation.R +++ b/analysis/plot_data_format_ablation.R @@ -530,7 +530,7 @@ fig_supp <- ggplot(gain_all_dists, aes(xintercept = ref_val), colour = "grey50", linetype = "dashed", - linewidth = 0.35, + linewidth = 0.35 ) + scale_colour_aaas() + scale_shape_manual( From 92ff13a4608fac89bae733a129632215f3407d32 Mon Sep 17 00:00:00 2001 From: cm401 Date: Mon, 25 May 2026 17:17:05 +0100 Subject: [PATCH 18/18] add documentation --- man/generate_simulation_table1.Rd | 32 ++++++++++++++++++++++++++ man/generate_simulation_table2.Rd | 32 ++++++++++++++++++++++++++ man/plot_simulation_study_figure2.Rd | 22 ++++++++++++++++++ man/plot_simulation_study_si_figure.Rd | 21 +++++++++++++++++ 4 files changed, 107 insertions(+) create mode 100644 man/generate_simulation_table1.Rd create mode 100644 man/generate_simulation_table2.Rd create mode 100644 man/plot_simulation_study_figure2.Rd create mode 100644 man/plot_simulation_study_si_figure.Rd diff --git a/man/generate_simulation_table1.Rd b/man/generate_simulation_table1.Rd new file mode 100644 index 0000000..d5962ba --- /dev/null +++ b/man/generate_simulation_table1.Rd @@ -0,0 +1,32 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/table_utils.R +\name{generate_simulation_table1} +\alias{generate_simulation_table1} +\title{Generate a LaTeX simulation performance table by distribution and summary type} +\usage{ +generate_simulation_table1( + summary_res, + caption = + paste0("Simulation study performance by distribution family and summary type. ", + "Coverage values are empirical 95\\\\\% credible interval coverage (\\\\\%; nominal target: 95\\\\\%). ", + "Predictive bias is the median signed error (posterior median $-$ true value) across replicates. ", + "IQD: integrated quadratic distance; WIS: weighted interval score (lower is better). ", + "Metrics are averaged over numbers of datasets and within-study sample sizes."), + label = "tab:sim_dist_summary" +) +} +\arguments{ +\item{summary_res}{Output of \code{\link[=create_results_summary]{create_results_summary()}}.} + +\item{caption}{LaTeX \code{\\caption\{\}} string.} + +\item{label}{LaTeX \code{\\label\{\}} string.} +} +\value{ +A character string with the complete LaTeX \code{longtable} environment. +} +\description{ +Produces a landscape \pkg{longtable} summarising simulation study performance +aggregated over numbers of datasets and within-study sample sizes, with one +row per distribution-family/summary-type combination, grouped by distribution. +} diff --git a/man/generate_simulation_table2.Rd b/man/generate_simulation_table2.Rd new file mode 100644 index 0000000..9e203e5 --- /dev/null +++ b/man/generate_simulation_table2.Rd @@ -0,0 +1,32 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/table_utils.R +\name{generate_simulation_table2} +\alias{generate_simulation_table2} +\title{Generate a LaTeX simulation performance table by data availability} +\usage{ +generate_simulation_table2( + summary_res, + caption = + paste0("Simulation study performance by summary type, number of datasets, and within-study sample size ", + "(pooled over distribution families). ", + "Coverage values are empirical 95\\\\\% credible interval coverage (\\\\\%; nominal target: 95\\\\\%). ", + "Predictive bias is the median signed error (posterior median $-$ true value) across replicates. ", + "IQD: integrated quadratic distance; WIS: weighted interval score (lower is better)."), + label = "tab:sim_data_availability" +) +} +\arguments{ +\item{summary_res}{Output of \code{\link[=create_results_summary]{create_results_summary()}}.} + +\item{caption}{LaTeX \code{\\caption\{\}} string.} + +\item{label}{LaTeX \code{\\label\{\}} string.} +} +\value{ +A character string with the complete LaTeX \code{longtable} environment. +} +\description{ +Produces a landscape \pkg{longtable} summarising simulation study performance +by summary type, number of datasets, and within-study sample size, pooled +over distribution families, with rows grouped by summary type. +} diff --git a/man/plot_simulation_study_figure2.Rd b/man/plot_simulation_study_figure2.Rd new file mode 100644 index 0000000..51aa025 --- /dev/null +++ b/man/plot_simulation_study_figure2.Rd @@ -0,0 +1,22 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/ploting_utils.R +\name{plot_simulation_study_figure2} +\alias{plot_simulation_study_figure2} +\title{Publication-ready simulation study figure (Figure 2)} +\usage{ +plot_simulation_study_figure2(summary_res, base_size = 9.5) +} +\arguments{ +\item{summary_res}{A data frame returned by \code{create_results_summary()}.} + +\item{base_size}{Base font size (default 9.5).} +} +\value{ +A \code{patchwork} object. +} +\description{ +Two-panel figure (A: coverage, B: bias) showing predictive performance +across distribution families, summary types, and number of contributing +datasets. Colour is mapped to summary type via the AAAS palette, matching +the scheme used elsewhere in the package. +} diff --git a/man/plot_simulation_study_si_figure.Rd b/man/plot_simulation_study_si_figure.Rd new file mode 100644 index 0000000..e78d81a --- /dev/null +++ b/man/plot_simulation_study_si_figure.Rd @@ -0,0 +1,21 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/ploting_utils.R +\name{plot_simulation_study_si_figure} +\alias{plot_simulation_study_si_figure} +\title{Supplementary information simulation study figure} +\usage{ +plot_simulation_study_si_figure(summary_res, base_size = 10) +} +\arguments{ +\item{summary_res}{A data frame returned by \code{create_results_summary()}.} + +\item{base_size}{Base font size (default 10).} +} +\value{ +A \code{patchwork} object. +} +\description{ +Three-panel figure (A: coverage, B: bias, C: mean WIS) using jittered points +coloured by number of datasets and shaped by observations per study. All +summary types (including Freq Table) are shown on the x-axis. +}