diff --git a/.Rbuildignore b/.Rbuildignore
index 3912071..e737dd8 100644
--- a/.Rbuildignore
+++ b/.Rbuildignore
@@ -1,3 +1,4 @@
 ^.*\.Rproj$
 ^\.Rproj\.user$
 ^\.github$
+^analysis$
diff --git a/NAMESPACE b/NAMESPACE
index 4ac80e2..937623b 100644
--- a/NAMESPACE
+++ b/NAMESPACE
@@ -44,6 +44,8 @@ export(generate_scenario_library)
 export(make_stan_init_fn)
 export(plot_main_figure)
 export(plot_main_figure_split)
+export(plot_simulation_study_figure)
+export(plot_simulation_study_sens_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 3b2068a..60d2ce2 100644
--- a/R/ploting_utils.R
+++ b/R/ploting_utils.R
@@ -1524,3 +1524,130 @@ plot_main_figure_split <- function(
       )
     )
 }
+
+
+# Internal core builder shared by the main and sensitivity simulation figures.
+# x_var:     column to put on the x-axis ("summary_type_label" or "dist_type").
+# facet_var: column to facet by, or NULL for no faceting.
+# color_var / shape_var: columns mapped to colour and shape aesthetics.
+.sim_figure_core <- function(summary_res, base_size = 9,
+                             x_var     = "summary_type_label",
+                             x_lab     = "Summary type",
+                             facet_var = "dist_type",
+                             color_var = "n_datasets_bucket",
+                             shape_var = "n_obs_bucket",
+                             color_lab = "N datasets",
+                             shape_lab = "N obs per study") {
+  dat <- .pred_buckets(summary_res)
+
+  .row <- function(data, y_col, ref_line = NULL, y_lab,
+                   hide_x = FALSE, hide_strip = FALSE) {
+    p <- ggplot(data,
+                aes(x     = .data[[x_var]],
+                    y     = .data[[y_col]],
+                    color = .data[[color_var]],
+                    shape = .data[[shape_var]]))
+
+    if (!is.null(ref_line))
+      p <- p + geom_hline(yintercept = ref_line, linetype = "dashed",
+                          color = "red", linewidth = 0.4)
+
+    p <- p +
+      geom_point(alpha    = 0.7,
+                 size     = 1,
+                 position = position_jitterdodge(jitter.width = 0.1,
+                                                 dodge.width  = 0.4)) +
+      scale_color_aaas() +
+      labs(x = if (hide_x) NULL else x_lab,
+           y = y_lab,
+           color = color_lab, shape = shape_lab) +
+      theme_minimal(base_size = base_size) +
+      theme(
+        axis.text.x  = if (hide_x) element_blank() else element_text(angle = 40, hjust = 1),
+        axis.ticks.x = if (hide_x) element_blank() else element_line(),
+        strip.text   = if (hide_strip) element_blank() else element_text()
+      )
+
+    if (!is.null(facet_var))
+      p <- p + facet_wrap(reformulate(facet_var), nrow = 1,
+                          labeller = label_value, scales = "free_x")
+
+    return(p)
+  }
+
+  p1 <- .row(dat, "mean_wis",             ref_line = NULL, y_lab = "Mean WIS",
+             hide_x = TRUE,  hide_strip = FALSE)
+  p2 <- .row(dat, "coverage_pred_median", ref_line = 0.95, y_lab = "Coverage - P50",
+             hide_x = TRUE,  hide_strip = TRUE)
+  p3 <- .row(dat, "bias_pred_median",     ref_line = 0,    y_lab = "Bias - P50 (days)",
+             hide_x = TRUE,  hide_strip = TRUE)
+  p4 <- .row(dat, "coverage_pred_q95",    ref_line = 0.95, y_lab = "Coverage - P95",
+             hide_x = TRUE,  hide_strip = TRUE)
+  p5 <- .row(dat, "bias_pred_q95",        ref_line = 0,    y_lab = "Bias - P95 (days)",
+             hide_x = FALSE, hide_strip = TRUE)
+
+  patchwork::wrap_plots(p1, p2, p3, p4, p5, ncol = 1) +
+    patchwork::plot_layout(guides = "collect") +
+    patchwork::plot_annotation(tag_levels = "A") &
+    theme(legend.position = "bottom",
+          plot.tag = element_text(face = "bold"))
+}
+
+#' Main-text simulation study figure
+#'
+#' Assembles a five-row composite figure from the simulation study summary,
+#' excluding sensitivity scenarios (\code{TauSens} / \code{Mu0Sens}).
+#'
+#' @param summary_res A data frame returned by \code{create_results_summary()}.
+#' @param base_size   Base font size passed to \code{theme_minimal()}.
+#' @return A \code{patchwork} object.
+#' @export
+plot_simulation_study_figure <- function(summary_res, base_size = 9) {
+  summary_res <- summary_res %>%
+    filter(!grepl("TauSens|Mu0Sens", scenario_name))
+  fig <- .sim_figure_core(summary_res, base_size = base_size)
+  fig & scale_shape_manual(values = c("5" = 4, "10" = 3, "20" = 8, "25+" = 5))
+}
+
+#' Supplementary simulation study figure — sensitivity scenarios
+#'
+#' Same layout as \code{plot_simulation_study_figure()} but restricted to
+#' \code{TauSens} and \code{Mu0Sens} scenarios.
+#'
+#' @param summary_res A data frame returned by \code{create_results_summary()}.
+#' @param base_size   Base font size passed to \code{theme_minimal()}.
+#' @return A \code{patchwork} object.
+#' @export
+plot_simulation_study_sens_figure <- function(summary_res, base_size = 9) {
+  summary_res <- summary_res %>%
+    filter(grepl("TauSens|Mu0Sens", scenario_name)) %>%
+    mutate(
+      sens_type = if_else(grepl("TauSens", scenario_name),
+                          "τ sensitivity", "μ₀ sensitivity"),
+      sens_value = case_when(
+        grepl("TauSens", scenario_name) ~
+          paste0("τ = ", stringr::str_extract(scenario_name, "(?<=tau)[0-9.]+")),
+        TRUE ~
+          stringr::str_replace(
+            scenario_name,
+            "^Mu0Sens_(?:lognormal|gamma|weibull|burr12|gengamma)_(.+)_D\\d+.*",
+            "μ₀ = \\1"
+          )
+      )
+    ) %>%
+    mutate(
+      sens_value = forcats::fct_reorder(sens_value, grepl("μ", sens_value))
+    )
+
+  fig <- .sim_figure_core(summary_res, base_size = base_size,
+                          x_var     = "dist_type",
+                          x_lab     = "Distribution",
+                          facet_var = NULL,
+                          color_var = "sens_value",
+                          shape_var = "sens_type",
+                          color_lab = "Parameter value",
+                          shape_lab = "Sensitivity")
+
+  fig & guides(color = guide_legend(nrow = 2, byrow = FALSE),
+               shape = guide_legend(nrow = 2, byrow = FALSE))
+}
diff --git a/R/table_utils.R b/R/table_utils.R
index e6fa70c..2f5c326 100644
--- a/R/table_utils.R
+++ b/R/table_utils.R
@@ -671,7 +671,7 @@ generate_data_table <- function(
 ) {
   stopifnot(is.list(datasets), length(datasets) > 0L, !is.null(names(datasets)))
 
-  col_spec <- "|C{4cm}|C{2.0cm}|C{2.5cm}|R{0.7cm}|C{5.5cm}|C{6cm}|"
+  col_spec <- "|C{4cm}|C{2.0cm}|C{2.5cm}|R{0.7cm}|C{7.3cm}|C{6cm}|"
 
   header_cells <- paste(
     "\\textbf{Reference}",
diff --git a/README.md b/README.md
index 155aa67..cbbcb62 100644
--- a/README.md
+++ b/README.md
@@ -6,10 +6,31 @@
 [![lifecycle](https://img.shields.io/badge/lifecycle-experimental-orange.svg)](https://lifecycle.r-lib.org/articles/stages.html#experimental)
 <!--[ ![DOI](https://zenodo.org/badge/DOI/XYZ)](https://doi.org/XYZ) -->
 
-Methods to compute delay distributions from summary statistics
+Bayesian hierarchical synthesis of incubation period distributions from individual-level data and published summary statistics.
 
 <br clear="right"/>
 
+## Overview
+
+The incubation period -- the interval between pathogen exposure and symptom onset -- is a critical parameter for quarantine policy and outbreak response, yet individual-level exposure data remain scarce in the published literature. For most pathogens, only summary statistics are available, and restricting inference to individual-level data alone would leave too few datasets for reliable between-study heterogeneity estimation.
+
+**ddsynth** introduces a Bayesian hierarchical framework that jointly models individual-level observations and published summary statistics under a unified federated analysis approach. Key features:
+
+- Accurately recovers incubation period distributions across a range of data availability scenarios
+- Outperforms approaches that use summary data alone
+- Quantifies between-study heterogeneity, including by outbreak country, pathogen variant, and exposure pathway
+- Covers 18 pathogens with outbreak potential, spanning six pathogen groups
+
+## Dataset coverage
+
+The package includes a curated dataset of incubation period studies spanning outbreaks worldwide:
+
+<img src="man/figures/map_figure.png" width="100%" alt="Geographic distribution of incubation period datasets across 18 pathogens" />
+
+Datasets cover six pathogen groups:
+
+<img src="man/figures/pathogen_grouping.png" width="100%" alt="Pathogen grouping diagram showing six disease categories with WHO R&D Blueprint priority pathogens highlighted" />
+
 ## Prerequisites
 
 ### R
diff --git a/analysis/run_meta_analysis.R b/analysis/run_meta_analysis.R
index 94cbe12..8efa316 100644
--- a/analysis/run_meta_analysis.R
+++ b/analysis/run_meta_analysis.R
@@ -159,18 +159,25 @@ dir.create(OUTDIR, showWarnings = FALSE, recursive = TRUE)
   )
 }
 
-# Extract tau posterior summary from the Stan lognormal fit.
-# Prefers the "all" slot (matching the unfiltered meta-analysis);
-# falls back to "filtered" if "all" is unavailable.
-.get_stan_tau <- function(main_res, pathogen_name) {
+# Extract tau posterior summary from the best-fitting Stan model.
+# Best distribution is taken from model_weights (highest stacking weight).
+# Prefers the "all" slot; falls back to "filtered" if "all" is unavailable.
+# Returns NULL when < 5 datasets were used (tau not reliably estimated).
+.get_stan_tau <- function(main_res, model_weights, pathogen_name) {
   pg <- main_res[[pathogen_name]]
   if (is.null(pg)) return(NULL)
 
-  slot_nm <- if (!is.null(pg[["all"]][["lognormal"]])) "all"
-             else if (!is.null(pg[["filtered"]][["lognormal"]])) "filtered"
+  # Identify best-fitting distribution for this pathogen
+  wts      <- model_weights[[pathogen_name]]
+  best_dist <- if (!is.null(wts) && length(wts) > 0L) names(which.max(wts)) else "lognormal"
+
+  slot_nm <- if (!is.null(pg[["all"]][[best_dist]])) "all"
+             else if (!is.null(pg[["filtered"]][[best_dist]])) "filtered"
              else return(NULL)
 
-  sfit <- pg[[slot_nm]][["lognormal"]][["fit"]]
+  fit_obj <- pg[[slot_nm]][[best_dist]]
+  if (length(fit_obj$datasets) < 5L) return(NULL)
+  sfit <- fit_obj$fit
   if (is.null(sfit)) return(NULL)
 
   s <- tryCatch(
@@ -180,10 +187,11 @@ dir.create(OUTDIR, showWarnings = FALSE, recursive = TRUE)
   if (is.null(s)) return(NULL)
 
   list(
-    median = s[1, "50%"],
-    lo     = s[1, "2.5%"],
-    hi     = s[1, "97.5%"],
-    slot   = slot_nm
+    median    = s[1, "50%"],
+    lo        = s[1, "2.5%"],
+    hi        = s[1, "97.5%"],
+    slot      = slot_nm,
+    best_dist = best_dist
   )
 }
 
@@ -288,13 +296,23 @@ message("\n", length(results), " pathogens analysed.")
 # Generate and save forest plots
 # =============================================================================
 
-# Load Stan results to annotate forest plots with Bayesian tau estimates.
+# Load Stan results and model weights for the tau comparison table.
 stan_res_path <- here("results", "main_results.rds")
-stan_res      <- if (file.exists(stan_res_path)) {
+weights_path  <- here("results", "model_weights.rds")
+
+stan_res <- if (file.exists(stan_res_path)) {
   message("Loading Stan results from ", stan_res_path)
   readRDS(stan_res_path)
 } else {
-  message("main_results.rds not found — forest plots will not include Stan tau")
+  message("main_results.rds not found — Stan tau will not be reported")
+  NULL
+}
+
+model_weights <- if (file.exists(weights_path)) {
+  message("Loading model weights from ", weights_path)
+  readRDS(weights_path)
+} else {
+  message("model_weights.rds not found — defaulting to lognormal for Stan tau")
   NULL
 }
 
@@ -319,7 +337,7 @@ summary_rows <- lapply(names(results), function(nm) {
   res <- results[[nm]]
   m   <- res$meta
 
-  stan_tau <- if (!is.null(stan_res)) .get_stan_tau(stan_res, nm) else NULL
+  stan_tau <- if (!is.null(stan_res)) .get_stan_tau(stan_res, model_weights, nm) else NULL
 
   data.frame(
     pathogen            = nm,
@@ -330,10 +348,11 @@ summary_rows <- lapply(names(results), function(nm) {
     ci_upper            = round(exp(m$upper.random), 2),
     I2_pct              = round(m$I2 * 100,          1),
     metamean_tau        = round(sqrt(m$tau2),         3),
-    stan_tau_median     = if (!is.null(stan_tau)) round(stan_tau$median, 3) else NA_real_,
-    stan_tau_lo         = if (!is.null(stan_tau)) round(stan_tau$lo,     3) else NA_real_,
-    stan_tau_hi         = if (!is.null(stan_tau)) round(stan_tau$hi,     3) else NA_real_,
-    stan_slot           = if (!is.null(stan_tau)) stan_tau$slot             else NA_character_,
+    stan_tau_median     = if (!is.null(stan_tau)) round(stan_tau$median,    3) else NA_real_,
+    stan_tau_lo         = if (!is.null(stan_tau)) round(stan_tau$lo,        3) else NA_real_,
+    stan_tau_hi         = if (!is.null(stan_tau)) round(stan_tau$hi,        3) else NA_real_,
+    stan_best_dist      = if (!is.null(stan_tau)) stan_tau$best_dist         else NA_character_,
+    stan_slot           = if (!is.null(stan_tau)) stan_tau$slot              else NA_character_,
     stringsAsFactors    = FALSE
   )
 })
@@ -345,9 +364,9 @@ write.csv(summary_tbl,
 
 saveRDS(results, file.path(OUTDIR, "meta_analysis_results.rds"))
 
-message("\nSummary table (metamean τ vs Stan log-normal τ):\n")
+message("\nSummary table (metamean τ vs Stan τ from best-fitting distribution):\n")
 print(summary_tbl[, c("pathogen", "k", "I2_pct", "metamean_tau",
-                       "stan_tau_median", "stan_tau_lo", "stan_tau_hi")],
+                       "stan_best_dist", "stan_tau_median", "stan_tau_lo", "stan_tau_hi")],
       row.names = FALSE)
 message("\nAll outputs written to: ", OUTDIR)
 
@@ -379,3 +398,340 @@ Discrepancies can arise from:
       towards zero relative to the REML estimator.
 ──────────────────────────────────────────────────────────────────────────────
 ")
+
+
+# =============================================================================
+# LaTeX table
+# =============================================================================
+
+generate_meta_analysis_table <- function(
+    summary_tbl,
+    caption = paste0(
+      "Classical meta-analysis of incubation period central estimates ",
+      "(random-effects model, log-transformed mean, \\texttt{sm = MLN}). ",
+      "$k$ = number of studies. ",
+      "Pooled mean and 95\\% confidence interval (CI) are back-transformed to days. ",
+      "$I^2$ measures the percentage of total variance due to between-study heterogeneity. ",
+      "$\\tau_{\\text{meta}}$ is the between-study SD on the log scale (REML). ",
+      "$\\tau_{\\text{Stan}}$ (median and 95\\% credible interval) is from the ",
+      "best-fitting Bayesian model (``all data'' slot), also on the log scale."
+    ),
+    label = "tab:meta_analysis"
+) {
+
+  dist_labels <- c(
+    lognormal = "Log-normal",
+    gamma     = "Gamma",
+    weibull   = "Weibull",
+    burr      = "Burr XII",
+    gengamma  = "Gen. Gamma"
+  )
+
+  pathogen_labels <- ddsynth:::.PATHOGEN_LABELS
+  group_map       <- ddsynth:::.RESULT_KEY_TO_GROUP
+  group_order     <- ddsynth:::.GROUP_ORDER
+
+  # Format a number: fixed decimal places, dash for NA
+  .fmt <- function(x, d = 2) ifelse(is.na(x), "--", formatC(x, digits = d, format = "f"))
+
+  # ── Build body rows grouped by epidemiological category ───────────────────
+  body_lines <- character(0L)
+  first_group <- TRUE
+
+  for (grp in group_order) {
+    keys_in_grp <- names(group_map)[group_map == grp]
+    rows_in_grp <- summary_tbl[summary_tbl$pathogen %in% keys_in_grp, , drop = FALSE]
+    if (nrow(rows_in_grp) == 0L) next
+
+    # Group header row
+    if (!first_group) body_lines <- c(body_lines, "\\midrule")
+    body_lines <- c(body_lines,
+      paste0("\\rowcolor{gray!15}",
+             "\\multicolumn{7}{l}{\\textit{", grp, "}} \\\\")
+    )
+    first_group <- FALSE
+
+    # Data rows (order matches group_order within group)
+    for (i in seq_len(nrow(rows_in_grp))) {
+      r        <- rows_in_grp[i, ]
+      plabel   <- pathogen_labels[r$pathogen]
+      if (is.na(plabel)) plabel <- r$pathogen
+
+      dist_lab <- dist_labels[r$stan_best_dist]
+      if (is.na(dist_lab)) dist_lab <- r$stan_best_dist
+
+      mean_ci  <- paste0(.fmt(r$pooled_mean), " [", .fmt(r$ci_lower), ", ", .fmt(r$ci_upper), "]")
+      stan_ci  <- if (!is.na(r$stan_tau_median))
+        paste0(.fmt(r$stan_tau_median, 3), " [", .fmt(r$stan_tau_lo, 3), ", ",
+               .fmt(r$stan_tau_hi, 3), "]")
+      else "--"
+
+      body_lines <- c(body_lines, paste0(
+        "\\quad ", plabel, " & ",
+        r$k, " & ",
+        mean_ci, " & ",
+        .fmt(r$I2_pct, 1), "\\% & ",
+        .fmt(r$metamean_tau, 3), " & ",
+        stan_ci, " & ",
+        dist_lab, " \\\\"
+      ))
+    }
+  }
+
+  # ── Column header ──────────────────────────────────────────────────────────
+  col_header <- paste0(
+    "\\rowcolor[HTML]{F0E68C}",
+    "Pathogen & $k$ & ",
+    "\\multicolumn{1}{c}{Pooled mean [95\\% CI] (days)} & ",
+    "\\multicolumn{1}{c}{$I^2$} & ",
+    "\\multicolumn{1}{c}{$\\tau_{\\text{meta}}$} & ",
+    "\\multicolumn{1}{c}{$\\tau_{\\text{Stan}}$ [95\\% CrI]} & ",
+    "\\multicolumn{1}{c}{Best model} \\\\"
+  )
+
+  # ── Assemble full table ────────────────────────────────────────────────────
+  lines <- c(
+    "% Requires \\usepackage{booktabs}, \\usepackage{longtable},",
+    "% \\usepackage[table]{xcolor} in preamble.",
+    "\\begin{longtable}{@{} l r l r r l l @{}}",
+    paste0("\\caption{", caption, "}\\label{", label, "} \\\\"),
+    "\\toprule",
+    col_header,
+    "\\midrule",
+    "\\endfirsthead",
+    "%",
+    "\\multicolumn{7}{l}{\\small\\textit{continued from previous page}} \\\\",
+    "\\toprule",
+    col_header,
+    "\\midrule",
+    "\\endhead",
+    "%",
+    "\\multicolumn{7}{r}{\\small\\textit{continued on next page}} \\\\",
+    "\\endfoot",
+    "%",
+    "\\bottomrule",
+    "\\endlastfoot",
+    "%",
+    body_lines,
+    "\\end{longtable}"
+  )
+
+  paste(lines, collapse = "\n")
+}
+
+tex <- generate_meta_analysis_table(summary_tbl)
+tex_path <- file.path(OUTDIR, "meta_analysis_table.tex")
+writeLines(tex, tex_path)
+message("LaTeX table written to: ", tex_path)
+cat(tex, "\n")
+
+
+# =============================================================================
+# Dataset composition table
+# =============================================================================
+
+# Classify a single dataset as "indiv" (D/E) or "sumstat" (A/B/C).
+.composition_type <- function(d) {
+  if (!is.null(d$freq_lower) && !is.null(d$freq_upper)) return("indiv")
+  if (!is.null(d$freq_value) && !is.null(d$freq_count)) return("indiv")
+  "sumstat"
+}
+
+# Aggregate composition stats from a list of datasets.
+# Returns list(k_total, N_total, k_indiv, N_indiv, k_sumstat, N_sumstat).
+.composition_stats <- function(datasets) {
+  k_total   <- 0L; N_total   <- 0L
+  k_indiv   <- 0L; N_indiv   <- 0L
+  k_sumstat <- 0L; N_sumstat <- 0L
+
+  for (d in datasets) {
+    n  <- if (!is.null(d$n)) d$n else if (!is.null(d$freq_count)) sum(d$freq_count) else NA_integer_
+    tp <- .composition_type(d)
+    k_total <- k_total + 1L
+    N_total <- N_total + if (!is.na(n)) n else 0L
+    if (tp == "indiv") {
+      k_indiv <- k_indiv + 1L; N_indiv <- N_indiv + if (!is.na(n)) n else 0L
+    } else {
+      k_sumstat <- k_sumstat + 1L; N_sumstat <- N_sumstat + if (!is.na(n)) n else 0L
+    }
+  }
+  list(k_total=k_total, N_total=N_total, k_indiv=k_indiv,
+       N_indiv=N_indiv, k_sumstat=k_sumstat, N_sumstat=N_sumstat)
+}
+
+# Format composition stats as a LaTeX data row string (6 cells, no pathogen label).
+.fmt_comp_cells <- function(s) {
+  paste(
+    s$k_total, format(s$N_total, big.mark = ","),
+    if (s$k_indiv > 0L) s$k_indiv else "--",
+    if (s$k_indiv > 0L) format(s$N_indiv, big.mark = ",") else "--",
+    if (s$k_sumstat > 0L) s$k_sumstat else "--",
+    if (s$k_sumstat > 0L) format(s$N_sumstat, big.mark = ",") else "--",
+    sep = " & "
+  )
+}
+
+generate_dataset_composition_table <- function(
+    full_registry,
+    caption = paste0(
+      "Composition of incubation period datasets included in the analysis. ",
+      "For each pathogen (and subgroup where applicable), $k$ is the number of datasets ",
+      "and $N$ is the total number of individual observations. ",
+      "\\textit{Individual-level data} denotes datasets provided as frequency tables ",
+      "(exact or interval-censored); \\textit{summary statistics only} denotes datasets ",
+      "reported as mean $\\pm$ SD, median $\\pm$ IQR, or median $\\pm$ range."
+    ),
+    label = "tab:dataset_composition"
+) {
+  pathogen_labels <- ddsynth:::.PATHOGEN_LABELS
+  subgroup_labels <- ddsynth:::.SUBGROUP_LABELS
+  group_map       <- ddsynth:::.RESULT_KEY_TO_GROUP
+  group_order     <- ddsynth:::.GROUP_ORDER
+
+  # Column header (two-row: group spans on top, k/N labels below)
+  hdr1 <- paste0(
+    "\\rowcolor[HTML]{F0E68C}",
+    " & \\multicolumn{2}{c}{\\textbf{Total}} &",
+    " \\multicolumn{2}{c}{\\textbf{Individual-level}} &",
+    " \\multicolumn{2}{c}{\\textbf{Summary statistics only}} \\\\"
+  )
+  hdr2 <- paste0(
+    "\\rowcolor[HTML]{F0E68C}",
+    "\\textbf{Pathogen} & $k$ & $N$ & $k$ & $N$ & $k$ & $N$ \\\\"
+  )
+
+  body_lines  <- character(0L)
+  first_group <- TRUE
+
+  # Index registry by pathogen name for quick lookup
+  reg_index <- stats::setNames(
+    lapply(full_registry, `[[`, "data"),
+    sapply(full_registry, `[[`, "name")
+  )
+
+  for (grp in group_order) {
+    keys_in_grp <- names(group_map)[group_map == grp]
+    keys_present <- intersect(keys_in_grp, names(reg_index))
+    if (length(keys_present) == 0L) next
+
+    if (!first_group) body_lines <- c(body_lines, "\\midrule")
+    body_lines <- c(body_lines,
+      paste0("\\rowcolor{gray!15}",
+             "\\multicolumn{7}{@{}l}{\\textbf{\\textit{", grp, "}}} \\\\")
+    )
+    first_group <- FALSE
+
+    first_in_group <- TRUE
+    for (nm in keys_in_grp) {
+      if (!nm %in% names(reg_index)) next
+      datasets <- reg_index[[nm]]
+      plabel   <- pathogen_labels[[nm]]
+      if (is.null(plabel) || is.na(plabel)) plabel <- nm
+
+      # Overall totals for this pathogen
+      stats_all <- .composition_stats(datasets)
+
+      # Detect subgroups (NULL entries → NA, then omit)
+      sgs <- unique(na.omit(sapply(datasets, function(d) {
+        sg <- d$subgroup; if (is.null(sg)) NA_character_ else sg
+      })))
+      unclass_datasets <- Filter(function(d) is.null(d$subgroup), datasets)
+
+      # Show subgroup rows only when ≥2 distinct subgroups OR ≥1 subgroup + unclassified.
+      # Skip when all datasets belong to a single subgroup with no unclassified entries
+      # (the parent row already captures everything).
+      has_subgroups <- (length(sgs) >= 2L) ||
+                       (length(sgs) == 1L && length(unclass_datasets) > 0L)
+
+      if (!first_in_group)
+        body_lines <- c(body_lines, "\\addlinespace[4pt]")
+
+      body_lines <- c(body_lines, paste0(
+        "\\textbf{", plabel, "} & ",
+        .fmt_comp_cells(stats_all), " \\\\"
+      ))
+      first_in_group <- FALSE
+
+      if (has_subgroups) {
+        for (sg in sgs) {
+          sg_datasets <- Filter(function(d) {
+            dsg <- d$subgroup; !is.null(dsg) && identical(dsg, sg)
+          }, datasets)
+          if (length(sg_datasets) == 0L) next
+          sg_stats <- .composition_stats(sg_datasets)
+          sg_lab <- if (!is.na(subgroup_labels[sg])) subgroup_labels[[sg]] else sg
+          body_lines <- c(body_lines, paste0(
+            "\\quad \\textit{", sg_lab, "} & ",
+            .fmt_comp_cells(sg_stats), " \\\\"
+          ))
+        }
+        if (length(unclass_datasets) > 0L) {
+          uc_stats <- .composition_stats(unclass_datasets)
+          body_lines <- c(body_lines, paste0(
+            "\\quad \\textit{Unclassified} & ",
+            .fmt_comp_cells(uc_stats), " \\\\"
+          ))
+        }
+      }
+    }
+  }
+
+  # Assemble full table
+  lines <- c(
+    "% Requires \\usepackage{booktabs}, \\usepackage{longtable},",
+    "% \\usepackage[table]{xcolor} in preamble.",
+    "\\begin{longtable}{@{} l r r r r r r @{}}",
+    paste0("\\caption{", caption, "}\\label{", label, "} \\\\"),
+    "\\toprule",
+    hdr1,
+    hdr2,
+    "\\midrule",
+    "\\endfirsthead",
+    "%",
+    "\\multicolumn{7}{l}{\\small\\textit{continued from previous page}} \\\\",
+    "\\toprule",
+    hdr1,
+    hdr2,
+    "\\midrule",
+    "\\endhead",
+    "%",
+    "\\multicolumn{7}{r}{\\small\\textit{continued on next page}} \\\\",
+    "\\endfoot",
+    "%",
+    "\\bottomrule",
+    "\\endlastfoot",
+    "%",
+    body_lines,
+    "\\end{longtable}"
+  )
+
+  paste(lines, collapse = "\n")
+}
+
+# Full registry covering all pathogens (for composition table)
+full_registry <- list(
+  list(name = "EVD",      data = datasets_EVD),
+  list(name = "MVD",      data = datasets_MVD),
+  list(name = "Lassa",    data = datasets_Lassa),
+  list(name = "CCHF",     data = datasets_CCHF_extended),
+  list(name = "COVID_19", data = datasets_COVID_19),
+  list(name = "SARS",     data = datasets_SARS),
+  list(name = "MERS",     data = datasets_MERS),
+  list(name = "Flu",      data = datasets_flu),
+  list(name = "Dengue",   data = datasets_Dengue),
+  list(name = "Zika",     data = datasets_Zika),
+  list(name = "RVF",      data = datasets_RVF),
+  list(name = "YFV",      data = datasets_YFV),
+  list(name = "Nipah",    data = datasets_Nipah),
+  list(name = "Mpox",     data = datasets_Mpox),
+  list(name = "Measles",  data = datasets_Measles),
+  list(name = "Smallpox", data = datasets_Smallpox),
+  list(name = "Cholera",  data = datasets_Cholera),
+  list(name = "Typhoid",  data = datasets_typhoid)
+)
+
+comp_tex      <- generate_dataset_composition_table(full_registry)
+comp_tex_path <- file.path(OUTDIR, "dataset_composition_table.tex")
+writeLines(comp_tex, comp_tex_path)
+message("Dataset composition table written to: ", comp_tex_path)
+cat(comp_tex, "\n")
diff --git a/man/figures/map_figure.png b/man/figures/map_figure.png
new file mode 100644
index 0000000..42280cf
Binary files /dev/null and b/man/figures/map_figure.png differ
diff --git a/man/figures/pathogen_grouping.png b/man/figures/pathogen_grouping.png
new file mode 100644
index 0000000..4bab29a
Binary files /dev/null and b/man/figures/pathogen_grouping.png differ
diff --git a/man/plot_simulation_study_figure.Rd b/man/plot_simulation_study_figure.Rd
new file mode 100644
index 0000000..dff400d
--- /dev/null
+++ b/man/plot_simulation_study_figure.Rd
@@ -0,0 +1,20 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/ploting_utils.R
+\name{plot_simulation_study_figure}
+\alias{plot_simulation_study_figure}
+\title{Main-text simulation study figure}
+\usage{
+plot_simulation_study_figure(summary_res, base_size = 9)
+}
+\arguments{
+\item{summary_res}{A data frame returned by \code{create_results_summary()}.}
+
+\item{base_size}{Base font size passed to \code{theme_minimal()}.}
+}
+\value{
+A \code{patchwork} object.
+}
+\description{
+Assembles a five-row composite figure from the simulation study summary,
+excluding sensitivity scenarios (\code{TauSens} / \code{Mu0Sens}).
+}
diff --git a/man/plot_simulation_study_sens_figure.Rd b/man/plot_simulation_study_sens_figure.Rd
new file mode 100644
index 0000000..cc8d993
--- /dev/null
+++ b/man/plot_simulation_study_sens_figure.Rd
@@ -0,0 +1,20 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/ploting_utils.R
+\name{plot_simulation_study_sens_figure}
+\alias{plot_simulation_study_sens_figure}
+\title{Supplementary simulation study figure — sensitivity scenarios}
+\usage{
+plot_simulation_study_sens_figure(summary_res, base_size = 9)
+}
+\arguments{
+\item{summary_res}{A data frame returned by \code{create_results_summary()}.}
+
+\item{base_size}{Base font size passed to \code{theme_minimal()}.}
+}
+\value{
+A \code{patchwork} object.
+}
+\description{
+Same layout as \code{plot_simulation_study_figure()} but restricted to
+\code{TauSens} and \code{Mu0Sens} scenarios.
+}
diff --git a/tests/testthat/test-datasets.R b/tests/testthat/test-datasets.R
index cca9ba3..f4cface 100644
--- a/tests/testthat/test-datasets.R
+++ b/tests/testthat/test-datasets.R
@@ -1,24 +1,35 @@
-# Smoke tests for built-in datasets: datasets_Nipah, datasets_MVD, datasets_EVD
-# Checks structure, recognised formats, and round-trip through prepare_stan_data_from_datasets().
+# Smoke tests for all built-in datasets.
+# Checks: correct type, expected entry count, all recognised summary formats,
+# and a round-trip through prepare_stan_data_from_datasets().
 
 # ---------------------------------------------------------------------------
-# Helper: verify every entry has at least one recognised summary-stat combination
+# Helpers
 # ---------------------------------------------------------------------------
+
 .check_recognised <- function(ds, ds_name) {
   for (nm in names(ds)) {
     d <- ds[[nm]]
-    has_range <- !is.null(d$median) && !is.null(d$min) && !is.null(d$max)
-    has_iqr   <- !is.null(d$median) && !is.null(d$Q1)  && !is.null(d$Q3)
+    has_range <- !is.null(d$median) && !is.null(d$min)  && !is.null(d$max)
+    has_iqr   <- !is.null(d$median) && !is.null(d$Q1)   && !is.null(d$Q3)
     has_mean  <- !is.null(d$mean)   && !is.null(d$sd)
     has_freq4 <- !is.null(d$freq_value) && !is.null(d$freq_count)
-    has_freq5 <- !is.null(d$freq_lower) && !is.null(d$freq_upper) && !is.null(d$freq_count)
+    has_freq5 <- !is.null(d$freq_lower) && !is.null(d$freq_upper) &&
+                 !is.null(d$freq_count)
     expect_true(has_range || has_iqr || has_mean || has_freq4 || has_freq5,
                 info = paste(ds_name, nm, "has no recognised format"))
   }
 }
 
+.round_trip <- function(ds, expected_n, ds_name) {
+  sd <- suppressWarnings(prepare_stan_data_from_datasets(ds, dist_type = 1))
+  expect_equal(sd$n_datasets, expected_n,
+               info = paste(ds_name, "n_datasets mismatch"))
+  expect_true(all(sd$n_obs > 0),
+              info = paste(ds_name, "has zero-n entry"))
+}
+
 # ---------------------------------------------------------------------------
-# datasets_Nipah
+# datasets_Nipah (11 entries)
 # ---------------------------------------------------------------------------
 
 test_that("datasets_Nipah is a named list with 11 entries", {
@@ -31,13 +42,11 @@ test_that("every datasets_Nipah entry has a recognised summary-stat format", {
 })
 
 test_that("datasets_Nipah round-trips through prepare_stan_data_from_datasets", {
-  sd <- prepare_stan_data_from_datasets(datasets_Nipah, dist_type = 1)
-  expect_equal(sd$n_datasets, 11L)
-  expect_true(all(sd$n_obs > 0))
+  .round_trip(datasets_Nipah, 11L, "datasets_Nipah")
 })
 
 # ---------------------------------------------------------------------------
-# datasets_MVD
+# datasets_MVD (2 entries)
 # ---------------------------------------------------------------------------
 
 test_that("datasets_MVD is a named list with 2 entries", {
@@ -50,13 +59,11 @@ test_that("every datasets_MVD entry has a recognised summary-stat format", {
 })
 
 test_that("datasets_MVD round-trips through prepare_stan_data_from_datasets", {
-  sd <- suppressWarnings(prepare_stan_data_from_datasets(datasets_MVD, dist_type = 1))
-  expect_equal(sd$n_datasets, 2L)
-  expect_true(all(sd$n_obs > 0))
+  .round_trip(datasets_MVD, 2L, "datasets_MVD")
 })
 
 # ---------------------------------------------------------------------------
-# datasets_EVD
+# datasets_EVD (11 entries)
 # ---------------------------------------------------------------------------
 
 test_that("datasets_EVD is a named list with 11 entries", {
@@ -69,23 +76,296 @@ test_that("every datasets_EVD entry has a recognised summary-stat format", {
 })
 
 test_that("datasets_EVD round-trips through prepare_stan_data_from_datasets", {
-  sd <- prepare_stan_data_from_datasets(datasets_EVD, dist_type = 1)
-  expect_equal(sd$n_datasets, 11L)
-  expect_true(all(sd$n_obs > 0))
+  .round_trip(datasets_EVD, 11L, "datasets_EVD")
+})
+
+# ---------------------------------------------------------------------------
+# datasets_SARS (22 entries)
+# ---------------------------------------------------------------------------
+
+test_that("datasets_SARS is a named list with 22 entries", {
+  expect_type(datasets_SARS, "list")
+  expect_equal(length(datasets_SARS), 22L)
+})
+
+test_that("every datasets_SARS entry has a recognised summary-stat format", {
+  .check_recognised(datasets_SARS, "datasets_SARS")
+})
+
+test_that("datasets_SARS round-trips through prepare_stan_data_from_datasets", {
+  .round_trip(datasets_SARS, 22L, "datasets_SARS")
+})
+
+# ---------------------------------------------------------------------------
+# datasets_MERS (11 entries)
+# ---------------------------------------------------------------------------
+
+test_that("datasets_MERS is a named list with 11 entries", {
+  expect_type(datasets_MERS, "list")
+  expect_equal(length(datasets_MERS), 11L)
+})
+
+test_that("every datasets_MERS entry has a recognised summary-stat format", {
+  .check_recognised(datasets_MERS, "datasets_MERS")
+})
+
+test_that("datasets_MERS round-trips through prepare_stan_data_from_datasets", {
+  .round_trip(datasets_MERS, 11L, "datasets_MERS")
+})
+
+# ---------------------------------------------------------------------------
+# datasets_Lassa (1 entry)
+# ---------------------------------------------------------------------------
+
+test_that("datasets_Lassa is a named list with 1 entry", {
+  expect_type(datasets_Lassa, "list")
+  expect_equal(length(datasets_Lassa), 1L)
+})
+
+test_that("every datasets_Lassa entry has a recognised summary-stat format", {
+  .check_recognised(datasets_Lassa, "datasets_Lassa")
+})
+
+test_that("datasets_Lassa round-trips through prepare_stan_data_from_datasets", {
+  .round_trip(datasets_Lassa, 1L, "datasets_Lassa")
+})
+
+# ---------------------------------------------------------------------------
+# datasets_Measles (12 entries)
+# ---------------------------------------------------------------------------
+
+test_that("datasets_Measles is a named list with 12 entries", {
+  expect_type(datasets_Measles, "list")
+  expect_equal(length(datasets_Measles), 12L)
+})
+
+test_that("every datasets_Measles entry has a recognised summary-stat format", {
+  .check_recognised(datasets_Measles, "datasets_Measles")
+})
+
+test_that("datasets_Measles round-trips through prepare_stan_data_from_datasets", {
+  .round_trip(datasets_Measles, 12L, "datasets_Measles")
+})
+
+# ---------------------------------------------------------------------------
+# datasets_Mpox (16 entries)
+# ---------------------------------------------------------------------------
+
+test_that("datasets_Mpox is a named list with 16 entries", {
+  expect_type(datasets_Mpox, "list")
+  expect_equal(length(datasets_Mpox), 16L)
+})
+
+test_that("every datasets_Mpox entry has a recognised summary-stat format", {
+  .check_recognised(datasets_Mpox, "datasets_Mpox")
+})
+
+test_that("datasets_Mpox round-trips through prepare_stan_data_from_datasets", {
+  .round_trip(datasets_Mpox, 16L, "datasets_Mpox")
+})
+
+# ---------------------------------------------------------------------------
+# datasets_Cholera (16 entries)
+# ---------------------------------------------------------------------------
+
+test_that("datasets_Cholera is a named list with 16 entries", {
+  expect_type(datasets_Cholera, "list")
+  expect_equal(length(datasets_Cholera), 16L)
+})
+
+test_that("every datasets_Cholera entry has a recognised summary-stat format", {
+  .check_recognised(datasets_Cholera, "datasets_Cholera")
+})
+
+test_that("datasets_Cholera round-trips through prepare_stan_data_from_datasets", {
+  .round_trip(datasets_Cholera, 16L, "datasets_Cholera")
+})
+
+# ---------------------------------------------------------------------------
+# datasets_RVF (2 entries)
+# ---------------------------------------------------------------------------
+
+test_that("datasets_RVF is a named list with 2 entries", {
+  expect_type(datasets_RVF, "list")
+  expect_equal(length(datasets_RVF), 2L)
+})
+
+test_that("every datasets_RVF entry has a recognised summary-stat format", {
+  .check_recognised(datasets_RVF, "datasets_RVF")
+})
+
+test_that("datasets_RVF round-trips through prepare_stan_data_from_datasets", {
+  .round_trip(datasets_RVF, 2L, "datasets_RVF")
+})
+
+# ---------------------------------------------------------------------------
+# datasets_CCHF (6 entries)
+# ---------------------------------------------------------------------------
+
+test_that("datasets_CCHF is a named list with 6 entries", {
+  expect_type(datasets_CCHF, "list")
+  expect_equal(length(datasets_CCHF), 6L)
+})
+
+test_that("every datasets_CCHF entry has a recognised summary-stat format", {
+  .check_recognised(datasets_CCHF, "datasets_CCHF")
+})
+
+test_that("datasets_CCHF round-trips through prepare_stan_data_from_datasets", {
+  .round_trip(datasets_CCHF, 6L, "datasets_CCHF")
+})
+
+# ---------------------------------------------------------------------------
+# datasets_COVID_19 (51 entries)
+# ---------------------------------------------------------------------------
+
+test_that("datasets_COVID_19 is a named list with 51 entries", {
+  expect_type(datasets_COVID_19, "list")
+  expect_equal(length(datasets_COVID_19), 51L)
+})
+
+test_that("every datasets_COVID_19 entry has a recognised summary-stat format", {
+  .check_recognised(datasets_COVID_19, "datasets_COVID_19")
+})
+
+test_that("datasets_COVID_19 round-trips through prepare_stan_data_from_datasets", {
+  .round_trip(datasets_COVID_19, 51L, "datasets_COVID_19")
+})
+
+# ---------------------------------------------------------------------------
+# datasets_Dengue (14 entries)
+# ---------------------------------------------------------------------------
+
+test_that("datasets_Dengue is a named list with 14 entries", {
+  expect_type(datasets_Dengue, "list")
+  expect_equal(length(datasets_Dengue), 14L)
+})
+
+test_that("every datasets_Dengue entry has a recognised summary-stat format", {
+  .check_recognised(datasets_Dengue, "datasets_Dengue")
+})
+
+test_that("datasets_Dengue round-trips through prepare_stan_data_from_datasets", {
+  .round_trip(datasets_Dengue, 14L, "datasets_Dengue")
+})
+
+# ---------------------------------------------------------------------------
+# datasets_YFV (3 entries)
+# ---------------------------------------------------------------------------
+
+test_that("datasets_YFV is a named list with 3 entries", {
+  expect_type(datasets_YFV, "list")
+  expect_equal(length(datasets_YFV), 3L)
+})
+
+test_that("every datasets_YFV entry has a recognised summary-stat format", {
+  .check_recognised(datasets_YFV, "datasets_YFV")
+})
+
+test_that("datasets_YFV round-trips through prepare_stan_data_from_datasets", {
+  .round_trip(datasets_YFV, 3L, "datasets_YFV")
+})
+
+# ---------------------------------------------------------------------------
+# datasets_flu (14 entries)
+# ---------------------------------------------------------------------------
+
+test_that("datasets_flu is a named list with 14 entries", {
+  expect_type(datasets_flu, "list")
+  expect_equal(length(datasets_flu), 14L)
+})
+
+test_that("every datasets_flu entry has a recognised summary-stat format", {
+  .check_recognised(datasets_flu, "datasets_flu")
+})
+
+test_that("datasets_flu round-trips through prepare_stan_data_from_datasets", {
+  .round_trip(datasets_flu, 14L, "datasets_flu")
+})
+
+# ---------------------------------------------------------------------------
+# datasets_typhoid (22 entries)
+# ---------------------------------------------------------------------------
+
+test_that("datasets_typhoid is a named list with 22 entries", {
+  expect_type(datasets_typhoid, "list")
+  expect_equal(length(datasets_typhoid), 22L)
+})
+
+test_that("every datasets_typhoid entry has a recognised summary-stat format", {
+  .check_recognised(datasets_typhoid, "datasets_typhoid")
+})
+
+test_that("datasets_typhoid round-trips through prepare_stan_data_from_datasets", {
+  .round_trip(datasets_typhoid, 22L, "datasets_typhoid")
+})
+
+# ---------------------------------------------------------------------------
+# datasets_Smallpox (4 entries)
+# ---------------------------------------------------------------------------
+
+test_that("datasets_Smallpox is a named list with 4 entries", {
+  expect_type(datasets_Smallpox, "list")
+  expect_equal(length(datasets_Smallpox), 4L)
+})
+
+test_that("every datasets_Smallpox entry has a recognised summary-stat format", {
+  .check_recognised(datasets_Smallpox, "datasets_Smallpox")
+})
+
+test_that("datasets_Smallpox round-trips through prepare_stan_data_from_datasets", {
+  .round_trip(datasets_Smallpox, 4L, "datasets_Smallpox")
+})
+
+# ---------------------------------------------------------------------------
+# datasets_Zika (2 entries)
+# ---------------------------------------------------------------------------
+
+test_that("datasets_Zika is a named list with 2 entries", {
+  expect_type(datasets_Zika, "list")
+  expect_equal(length(datasets_Zika), 2L)
+})
+
+test_that("every datasets_Zika entry has a recognised summary-stat format", {
+  .check_recognised(datasets_Zika, "datasets_Zika")
+})
+
+test_that("datasets_Zika round-trips through prepare_stan_data_from_datasets", {
+  .round_trip(datasets_Zika, 2L, "datasets_Zika")
 })
 
 # ---------------------------------------------------------------------------
-# source fields (all datasets)
+# Source fields: all datasets
 # ---------------------------------------------------------------------------
 
 test_that("every source field that is present is a non-empty character string", {
-  all_ds <- list(Nipah = datasets_Nipah, MVD = datasets_MVD, EVD = datasets_EVD)
+  all_ds <- list(
+    Nipah    = datasets_Nipah,
+    MVD      = datasets_MVD,
+    EVD      = datasets_EVD,
+    SARS     = datasets_SARS,
+    MERS     = datasets_MERS,
+    Lassa    = datasets_Lassa,
+    Measles  = datasets_Measles,
+    Mpox     = datasets_Mpox,
+    Cholera  = datasets_Cholera,
+    RVF      = datasets_RVF,
+    CCHF     = datasets_CCHF,
+    COVID_19 = datasets_COVID_19,
+    Dengue   = datasets_Dengue,
+    YFV      = datasets_YFV,
+    flu      = datasets_flu,
+    typhoid  = datasets_typhoid,
+    Smallpox = datasets_Smallpox,
+    Zika     = datasets_Zika
+  )
   for (ds_name in names(all_ds)) {
     for (nm in names(all_ds[[ds_name]])) {
       src <- all_ds[[ds_name]][[nm]]$source
       if (!is.null(src)) {
         expect_type(src, "character")
-        expect_gt(nchar(src), 0)
+        expect_gt(nchar(src), 0,
+                  label = paste(ds_name, nm, "source"))
       }
     }
   }
diff --git a/tests/testthat/test-should-attempt-gg.R b/tests/testthat/test-should-attempt-gg.R
new file mode 100644
index 0000000..d7068dd
--- /dev/null
+++ b/tests/testthat/test-should-attempt-gg.R
@@ -0,0 +1,153 @@
+# Tests for should_attempt_gg() -----------------------------------------------
+
+# Helpers ---------------------------------------------------------------------
+
+# Frequency-table dataset (type-4): CV ≈ 0.24, "rich"
+.gg_freq <- function(n = 30) {
+  list(freq_value = c(5, 7, 9, 11, 13),
+       freq_count = c(3, 8, 12, 8, 3), n = n)
+}
+
+# Interval-censored dataset (type-5): "rich"
+.gg_icens <- function(n = 30) {
+  list(freq_lower = c(4, 6, 8, 10, 12),
+       freq_upper = c(6, 8, 10, 12, 14),
+       freq_count = c(3, 8, 12, 8, 3), n = n)
+}
+
+# Range dataset (type-1): CV ≈ 0.22, not "rich"
+.gg_range <- function(n = 30) {
+  list(median = 9, min = 5, max = 13, n = n)
+}
+
+# IQR dataset (type-2): CV ≈ 0.30, not "rich"
+.gg_iqr <- function(n = 30) {
+  list(median = 9, Q1 = 6.3, Q3 = 11.7, n = n)
+}
+
+# Mean+SD dataset (type-3): not "rich"
+.gg_meansd <- function(cv = 0.30, mean = 9, n = 30) {
+  list(mean = mean, sd = mean * cv, n = n)
+}
+
+# ── 1. Returns TRUE when richness and spread are both acceptable ──────────────
+
+test_that("returns TRUE with sufficient rich datasets and low CV spread", {
+  ds <- list(d1 = .gg_freq(), d2 = .gg_freq(), d3 = .gg_range())
+  # rich_frac = 2/3 ≈ 0.67 >= 0.30; all CVs similar => spread << 2.5
+  expect_true(should_attempt_gg(ds, verbose = FALSE))
+})
+
+test_that("returns TRUE with only frequency-table datasets", {
+  ds <- list(d1 = .gg_freq(), d2 = .gg_freq())
+  expect_true(should_attempt_gg(ds, verbose = FALSE))
+})
+
+test_that("type-5 interval-censored datasets count as rich", {
+  ds <- list(d1 = .gg_icens(), d2 = .gg_range())
+  # rich_frac = 0.5 >= 0.30
+  expect_true(should_attempt_gg(ds, verbose = FALSE))
+})
+
+test_that("returns TRUE for a single dataset (no spread check possible)", {
+  ds <- list(d1 = .gg_freq())
+  expect_true(should_attempt_gg(ds, verbose = FALSE))
+})
+
+# ── 2. Returns FALSE when rich fraction is too low ───────────────────────────
+
+test_that("returns FALSE when no datasets are rich (all type-1/2/3)", {
+  ds <- list(d1 = .gg_range(), d2 = .gg_range(), d3 = .gg_iqr(), d4 = .gg_meansd())
+  # rich_frac = 0 < 0.30
+  expect_false(should_attempt_gg(ds, verbose = FALSE))
+})
+
+test_that("returns FALSE when rich fraction is below custom threshold", {
+  ds <- list(d1 = .gg_freq(), d2 = .gg_range(), d3 = .gg_range(), d4 = .gg_range())
+  # rich_frac = 0.25; default min_rich_fraction = 0.30 → FALSE
+  expect_false(should_attempt_gg(ds, verbose = FALSE))
+})
+
+test_that("returns TRUE when rich fraction exactly meets the threshold", {
+  # 1 of 3 datasets is rich: 1/3 ≈ 0.33 >= 0.30
+  ds <- list(d1 = .gg_freq(), d2 = .gg_range(), d3 = .gg_range())
+  expect_true(should_attempt_gg(ds, verbose = FALSE))
+})
+
+test_that("custom min_rich_fraction threshold is respected", {
+  ds <- list(d1 = .gg_freq(), d2 = .gg_range())  # rich_frac = 0.5
+  expect_true( should_attempt_gg(ds, min_rich_fraction = 0.40, verbose = FALSE))
+  expect_false(should_attempt_gg(ds, min_rich_fraction = 0.60, verbose = FALSE))
+})
+
+# ── 3. Returns FALSE when CV spread is too large ─────────────────────────────
+
+test_that("returns FALSE when CV spread exceeds threshold", {
+  # d1 CV ≈ 0.05 (concentrated), d2 CV ≈ 0.69 (dispersed) → spread ≈ 13.9
+  d_conc <- list(freq_value = c(8, 9, 10), freq_count = c(1, 8, 1), n = 10)
+  d_disp <- list(freq_value = c(2, 9, 20), freq_count = c(3, 4, 3), n = 10)
+  ds <- list(d1 = d_conc, d2 = d_disp)
+  expect_false(should_attempt_gg(ds, verbose = FALSE))
+})
+
+test_that("custom cv_spread_threshold is respected", {
+  # spread ≈ 2; just straddles default of 2.5
+  d1 <- list(freq_value = c(7, 9, 11), freq_count = c(2, 6, 2), n = 10)  # CV ≈ 0.19
+  d2 <- list(freq_value = c(4, 9, 14), freq_count = c(2, 6, 2), n = 10)  # CV ≈ 0.38
+  ds <- list(d1 = d1, d2 = d2)
+  expect_true( should_attempt_gg(ds, cv_spread_threshold = 3.0, verbose = FALSE))
+  expect_false(should_attempt_gg(ds, cv_spread_threshold = 1.5, verbose = FALSE))
+})
+
+# ── 4. CV spread check is skipped for a single valid CV ───────────────────────
+
+test_that("spread check is skipped when only one dataset has a valid CV", {
+  # Only one CV computable; no spread can be formed → proceed to richness check
+  ds <- list(
+    d1 = .gg_freq(),
+    d2 = list(n = 10)  # no stat fields → NA CV
+  )
+  expect_true(should_attempt_gg(ds, verbose = FALSE))
+})
+
+# ── 5. Verbose messaging ──────────────────────────────────────────────────────
+
+test_that("verbose=TRUE emits SKIP when richness is too low", {
+  ds <- list(d1 = .gg_range(), d2 = .gg_range())
+  expect_message(should_attempt_gg(ds, verbose = TRUE), "SKIP")
+})
+
+test_that("verbose=TRUE emits SKIP when CV spread is too large", {
+  d_conc <- list(freq_value = c(8, 9, 10), freq_count = c(1, 8, 1), n = 10)
+  d_disp <- list(freq_value = c(2, 9, 20), freq_count = c(3, 4, 3), n = 10)
+  ds <- list(d1 = d_conc, d2 = d_disp)
+  expect_message(should_attempt_gg(ds, verbose = TRUE), "SKIP")
+})
+
+test_that("verbose=TRUE emits OK when all checks pass", {
+  ds <- list(d1 = .gg_freq(), d2 = .gg_freq())
+  expect_message(should_attempt_gg(ds, verbose = TRUE), "OK")
+})
+
+test_that("verbose=FALSE suppresses all messages", {
+  ds <- list(d1 = .gg_range(), d2 = .gg_range())
+  expect_no_message(should_attempt_gg(ds, verbose = FALSE))
+})
+
+# ── 6. NA handling ────────────────────────────────────────────────────────────
+
+test_that("datasets with no recognisable format produce NA CV without error", {
+  ds <- list(
+    d1 = .gg_freq(),
+    d2 = .gg_freq(),
+    d3 = list(n = 10)  # no stat fields
+  )
+  expect_true(should_attempt_gg(ds, verbose = FALSE))
+})
+
+test_that("all-NA CVs bypass the spread check and proceed to richness", {
+  # All entries have unrecognised format → no valid CV → no spread check
+  ds <- list(d1 = list(n = 10), d2 = list(n = 15))
+  # rich_frac = 0 < 0.30 → FALSE (fails richness, not spread)
+  expect_false(should_attempt_gg(ds, verbose = FALSE))
+})
diff --git a/tests/testthat/test-utils.R b/tests/testthat/test-utils.R
index 929727a..6e968b5 100644
--- a/tests/testthat/test-utils.R
+++ b/tests/testthat/test-utils.R
@@ -12,3 +12,57 @@ test_that("check_scalar rejects non-scalars", {
   expect_error(ddsynth:::check_scalar("a"),     "single finite")
   expect_silent(ddsynth:::check_scalar(3.14))
 })
+
+# make_stan_init_fn -----------------------------------------------------------
+
+test_that("make_stan_init_fn returns a function", {
+  sd <- suppressWarnings(
+    prepare_stan_data_from_datasets(
+      list(d1 = list(median = 7, min = 3, max = 14, n = 20)),
+      dist_type = 1
+    )
+  )
+  init_fn <- make_stan_init_fn(sd)
+  expect_true(is.function(init_fn))
+})
+
+test_that("make_stan_init_fn closure returns a list with required names", {
+  sd <- suppressWarnings(
+    prepare_stan_data_from_datasets(
+      list(d1 = list(median = 7, min = 3, max = 14, n = 20),
+           d2 = list(median = 9, min = 5, max = 16, n = 30)),
+      dist_type = 1
+    )
+  )
+  init <- make_stan_init_fn(sd)()
+  expect_named(init, c("mu0", "log_tau", "log_phi", "log_kappa", "loc_d_raw"),
+               ignore.order = TRUE)
+})
+
+test_that("make_stan_init_fn initialises parameters from stan_data", {
+  sd <- suppressWarnings(
+    prepare_stan_data_from_datasets(
+      list(d1 = list(median = 7, min = 3, max = 14, n = 20),
+           d2 = list(median = 9, min = 5, max = 16, n = 30)),
+      dist_type = 1
+    )
+  )
+  init <- make_stan_init_fn(sd)()
+  expect_equal(init$mu0,       sd$mu0_mean)
+  expect_equal(init$log_tau,   sd$log_tau_mean)
+  expect_equal(init$log_phi,   sd$log_phi_mean)
+  expect_equal(init$log_kappa, sd$log_kappa_mean)
+})
+
+test_that("make_stan_init_fn sets loc_d_raw to zero vector of length n_datasets", {
+  sd <- suppressWarnings(
+    prepare_stan_data_from_datasets(
+      list(d1 = list(median = 7, min = 3, max = 14, n = 20),
+           d2 = list(median = 9, min = 5, max = 16, n = 30),
+           d3 = list(mean   = 8, sd  = 2,           n = 25)),
+      dist_type = 1
+    )
+  )
+  init <- make_stan_init_fn(sd)()
+  expect_equal(init$loc_d_raw, rep(0.0, 3L))
+})