feat(app): add QMRA wiki library, pooling, and pooled-data exploration

.Rbuildignore

@@ -13,3 +13,4 @@
 ^app$
 ^AGENTS\.md$
 ^CLAUDE\.md$
+^vignettes/articles$

.gitignore

@@ -40,6 +40,7 @@ docs/
 # Local-only development artifacts and source data (kept out of the public repo)
 /dev/
 /data/raw/QMRA wiki data_v4.xlsx
+/QMRA wiki data_v4.xlsx
 
 # Translation temp files
 po/*~

R/pooling.R

@@ -14,8 +14,9 @@
 #'   [plot_dose_response()] with a fit and its bootstrap to draw bands, because a
 #'   pooled fit is not a `qdr_analysis`.
 #' @details The datasets are not tested for poolability here; run
-#'   [poolability_test()] first. Check metadata compatibility as well, since it
-#'   is not enforced. Single-dataset input is allowed and fits that dataset alone.
+#'   [poolability_test()] first. Datasets that carry `dose_unit` metadata must
+#'   share one unit (see [poolability_test()]); host and endpoint are not
+#'   checked. Single-dataset input is allowed and fits that dataset alone.
 #' @examples
 #' path <- system.file("extdata", "pooling-example.csv", package = "singlehit")
 #' raw <- readr::read_csv(path, show_col_types = FALSE)
@@ -55,8 +56,14 @@ fit_pooled_dose_response <- function(datasets, models = c("exponential", "beta_p
 #'   `deviance_pooled`, `deviance_unpooled`, `lrt_statistic`, `df`,
 #'   `chi_square_critical`, `p_value`, `poolable`, `converged`, and a
 #'   human-readable `conclusion`.
-#' @details Check that hosts, dose units, and endpoints match before testing.
-#' The current pooling functions do not enforce metadata compatibility.
+#' @details Doses are only comparable within one unit, so `dose_unit` is
+#' enforced: when any dataset carries a `dose_unit` column, every dataset must
+#' carry one and all must name the same unit, compared without regard to case
+#' or surrounding spaces ("CFU" and "cfu" match; "CFU" and "log10 CFU" do
+#' not). Otherwise an error names each dataset's unit. Within one dataset,
+#' [as_dose_response()] requires a single exact value. Datasets without a
+#' `dose_unit` column are not checked. Hosts and endpoints are not enforced;
+#' check that they match before testing.
 #' `poolable = TRUE` means the test did not detect a significantly worse pooled
 #' fit at `alpha`; it does not prove equivalence or adequate absolute fit.
 #' Inspect `converged` before interpreting the result.
@@ -132,7 +139,8 @@ poolability_test <- function(datasets, models = c("exponential", "beta_poisson")
 #'   Each row describes one combination under one model.
 #' @details A passing pair need not remain poolable after another dataset is
 #' added. Inspect the row for the complete combination you intend to pool.
-#' Check metadata compatibility before testing; it is not enforced here.
+#' Datasets must share one `dose_unit` when they carry one, as in
+#' [poolability_test()]; hosts and endpoints are not checked.
 #' Use [group_datasets()] to choose non-overlapping groups from these trials.
 #' @examples
 #' path <- system.file("extdata", "pooling-example.csv", package = "singlehit")
@@ -197,7 +205,8 @@ poolability_combinations <- function(
 #' @details Group numbers are labels within a model, not rankings. This function
 #' returns assignments, not pooled datasets or fitted models. Use
 #' [poolability_combinations()] to inspect overlapping candidate combinations.
-#' Check metadata compatibility before grouping; it is not enforced here.
+#' Datasets must share one `dose_unit` when they carry one, as in
+#' [poolability_test()]; hosts and endpoints are not checked.
 #' @examples
 #' path <- system.file("extdata", "pooling-example.csv", package = "singlehit")
 #' raw <- readr::read_csv(path, show_col_types = FALSE)
@@ -259,14 +268,45 @@ standardize_datasets <- function(datasets) {
     stop("`datasets` must be a non-empty list of dose-response datasets.", call. = FALSE)
   }
   names(datasets) <- dataset_names(datasets)
-  purrr::imap(datasets, function(dataset, name) {
+  standardized <- purrr::imap(datasets, function(dataset, name) {
     tryCatch(
       as_dose_response(dataset),
       error = function(cnd) {
         stop(sprintf("Dataset '%s': %s", name, conditionMessage(cnd)), call. = FALSE)
       }
     )
   })
+  check_shared_dose_unit(standardized)
+  standardized
+}
+
+# Doses in different units are not comparable, so pooled or jointly tested
+# datasets must share one. Datasets without unit metadata are not checked.
+check_shared_dose_unit <- function(standardized) {
+  has_unit <- vapply(standardized, function(dataset) "dose_unit" %in% names(dataset), logical(1))
+  if (!any(has_unit)) {
+    return(invisible(standardized))
+  }
+  if (!all(has_unit)) {
+    stop(
+      sprintf(
+        "Only some datasets record `dose_unit` (missing from %s). Add `dose_unit` to every dataset or remove it from all of them.",
+        paste0("'", names(standardized)[!has_unit], "'", collapse = ", ")
+      ),
+      call. = FALSE
+    )
+  }
+  units <- vapply(standardized, function(dataset) as.character(dataset$dose_unit[[1L]]), character(1))
+  if (length(unique(tolower(trimws(units)))) > 1L) {
+    stop(
+      sprintf(
+        "Datasets record different dose units, so their doses are not comparable: %s. Pool only datasets that share a unit.",
+        paste0("'", names(standardized), "' (", units, ")", collapse = ", ")
+      ),
+      call. = FALSE
+    )
+  }
+  invisible(standardized)
 }
 
 dataset_names <- function(datasets) {

README.md

@@ -211,9 +211,11 @@ poolability_test(collection$datasets)
 ```
 
 This list is an organizational convention, not a validated collection object.
-Check that hosts, dose units, and endpoints match before testing poolability;
-the current pooling functions do not enforce those checks. Recording matching
-metadata does not by itself establish that experiments should be pooled.
+The pooling functions refuse datasets whose `dose_unit` metadata differ, since
+doses in different units are not comparable. They do not check hosts or
+endpoints, so confirm those match before testing poolability. Recording
+matching metadata does not by itself establish that experiments should be
+pooled.
 
 ## Example
 
@@ -288,6 +290,8 @@ dataset separately and the stacked combination, then compares the deviance
 difference to a chi-squared distribution (per model). Datasets that pool
 significantly worse than when fit separately are kept distinct.
 `group_datasets()` extends this to find which trials are mutually poolable.
+Datasets that record `dose_unit` must all record the same unit; mixed units are
+an error.
 
 ```r
 trials <- list(trial_a = data_a, trial_b = data_b, trial_c = data_c)
@@ -313,6 +317,16 @@ trials because their searches grow rapidly. These are exploratory tests without
 multiple-testing adjustment; passing a test does not prove equivalence.
 Compare only experiments with matching hosts, units, and endpoints.
 
+Decide which experiments belong together (same isolate, outcome, host, and
+route) before testing, and use `poolability_combinations()` and
+`group_datasets()` to check that choice rather than to make it: with few doses
+per study the test rarely rejects, different searches return different
+groupings, and a group chosen from the results biases its pooled fit toward
+pooling. Also check that experiments do not share subjects, since the test
+assumes independent samples. The
+[pooling article](https://seanthimons.github.io/singlehit/articles/pooling.html)
+walks through these points with real Cryptosporidium challenge data.
+
 Datasets are combined by **stacking** — each trial's dose groups are kept as
 separate binomial observations, so repeated doses across trials are preserved
 (matching the QMRA-wiki pooled-experiment convention), rather than summed.

app/DESCRIPTION

@@ -22,16 +22,22 @@ Imports:
     generics,
     ggplot2,
     golem,
+    grDevices,
     mirai (>= 2.5.0),
     parallel,
     plotly,
     readr,
+    readxl,
+    scales,
     shiny,
-    singlehit (>= 0.3.0)
+    singlehit (>= 0.3.0),
+    stats,
+    tools
 Suggests:
     pkgload,
     testthat (>= 3.1.7),
-    withr
+    withr,
+    writexl
 Remotes:
     seanthimons/singlehit
 Config/testthat/edition: 3

app/NAMESPACE

@@ -2,6 +2,7 @@
 
 export(run_app)
 import(shiny)
+importFrom(ggplot2,.data)
 importFrom(golem,
   activate_js,
   add_resource_path,

app/R/app_server.R

@@ -5,12 +5,35 @@
 #' @import shiny
 #' @noRd
 app_server <- function(input, output, session) {
-  data <- mod_upload_server("upload")
+  # Data loads both the workbook (explored on Library) and the dataset Fit
+  # analyses; Library sends an experiment back. The reactiveVal breaks that
+  # cycle: Data is built first and reads whatever Library last sent.
+  sent <- reactiveVal(NULL)
+  upload <- mod_upload_server("upload", incoming = sent)
+  library_tab <- mod_library_server("library", upload$workbook)
+  data <- upload$data
+  observeEvent(library_tab$sent(), {
+    sent(library_tab$sent())
+    # A sent experiment replaces the current dataset; show it there first.
+    bslib::nav_select("main", "Data", session = session)
+  })
+  observeEvent(upload$workbook_loaded(), bslib::nav_select("main", "Library", session = session), ignoreInit = TRUE)
+  pooling <- mod_pooling_server(
+    "pooling", upload$workbook, library_tab$experiments, library_tab$pool_set,
+    upload$workbook_name, upload$workbook_loaded
+  )
+  observeEvent(library_tab$pool_set(), bslib::nav_select("main", "Pooling", session = session))
+  observeEvent(pooling$sent(), {
+    sent(pooling$sent())
+    bslib::nav_select("main", "Data", session = session)
+  })
   analysis <- mod_fit_server("fit", data)
   mod_compare_server("compare", analysis)
   # Plots and Explore get the analysis with the bootstraps attached once a run
   # finishes; Compare does not use them.
   analysis_with_bands <- mod_bootstrap_server("bootstrap", analysis)
-  mod_figures_server("figures", analysis_with_bands)
-  mod_explore_server("explore", analysis_with_bands)
+  # Pooled data also carry each experiment's dose groups, for colouring and
+  # exploring the observed points by sub-study.
+  mod_figures_server("figures", analysis_with_bands, upload$components)
+  mod_explore_server("explore", analysis_with_bands, upload$components)
 }

app/R/app_ui.R

@@ -12,11 +12,14 @@ app_ui <- function(request) {
     # run starting can make the app unresponsive; this says why.
     shiny::useBusyIndicators(),
     bslib::page_navbar(
+      id = "main",
       title = "singlehit",
       bslib::nav_panel(
         "Data",
         bslib::layout_sidebar(sidebar = bslib::sidebar(mod_upload_ui("upload")), mod_upload_table_ui("upload"))
       ),
+      bslib::nav_panel("Library", mod_library_ui("library")),
+      bslib::nav_panel("Pooling", mod_pooling_ui("pooling")),
       bslib::nav_panel(
         "Fit",
         bslib::layout_sidebar(sidebar = bslib::sidebar(mod_fit_ui("fit")), mod_fit_results_ui("fit"))

app/R/fct_explore.R

@@ -17,11 +17,16 @@ explore_dataset_choices <- c(
 #'   `lower`/`upper` bounds and `level` as a factor (widest interval first); otherwise they are `NA`.
 #' - `model_residuals`: observed groups with fitted values and deviance
 #'   residuals per model.
+#' - `substudies` (pooled data only): each experiment's dose groups with its
+#'   metadata and the pooled model's fitted probability at that dose, one row
+#'   per experiment, dose, and model.
 #'
 #' @param analysis A `qdr_analysis`.
+#' @param components Per-experiment dose groups of pooled data that match the
+#'   analysis (see [pool_components()]), or `NULL`.
 #' @return A named list of data frames.
 #' @noRd
-explore_datasets <- function(analysis) {
+explore_datasets <- function(analysis, components = NULL) {
   fits <- analysis$fits
   models <- names(fits)
 
@@ -51,11 +56,28 @@ explore_datasets <- function(analysis) {
     )
   }))
 
-  list(
+  out <- list(
     observed = as.data.frame(analysis$data),
     curves = curves,
     model_residuals = model_residuals
   )
+  if (!is.null(components)) {
+    out$substudies <- do.call(rbind, lapply(models, function(name) {
+      cbind(
+        components,
+        model = model_label_ui(name),
+        fitted = stats::predict(fits[[name]], newdata = components$dose),
+        stringsAsFactors = FALSE
+      )
+    }))
+  }
+  out
+}
+
+#' Explore dataset choices, with the sub-study observations when available
+#' @noRd
+explore_choices <- function(has_substudies) {
+  if (has_substudies) c(explore_dataset_choices, "Sub-study observations (pooled data)" = "substudies") else explore_dataset_choices
 }
 
 #' Input ids inside esquisse's module, relative to the Explore module

app/R/fct_figures.R

@@ -17,9 +17,10 @@ figure_types <- c(
 #' @param model Model name for `type = "single"`.
 #' @return A ggplot object.
 #' @noRd
-build_figure <- function(analysis, type = c("overlay", "panels", "single"), model = NULL) {
+build_figure <- function(analysis, type = c("overlay", "panels", "single"), model = NULL,
+                         components = NULL, colour_by = NULL) {
   type <- match.arg(type)
-  switch(
+  plot <- switch(
     type,
     overlay = singlehit::plot_model_overlay(analysis),
     panels = ggplot2::autoplot(analysis),
@@ -30,6 +31,95 @@ build_figure <- function(analysis, type = c("overlay", "panels", "single"), mode
       singlehit::plot_dose_response(analysis$fits[[model]], bootstrap = analysis$bootstraps[[model]])
     }
   )
+  if (!is.null(components) && !is.null(colour_by) && colour_by %in% substudy_fields) {
+    plot <- recolour_observed(plot, analysis$data, components, colour_by, aesthetic = if (type == "overlay") "fill" else "colour")
+  }
+  plot
+}
+
+#' Sub-study fields the observed points of pooled data can be coloured by
+#' @noRd
+substudy_fields <- c(
+  "Experiment" = "experiment",
+  "Pathogen" = "pathogen",
+  "Strain or isolate" = "strain",
+  "Host" = "host",
+  "Exposure route" = "route",
+  "Outcome" = "outcome"
+)
+
+#' Whether per-experiment dose groups add up to the analysed data
+#'
+#' Pooled data reach Fit with counts at shared doses added together; the
+#' components are the same counts per experiment. They belong to an analysis
+#' only when their per-dose sums equal its data, which guards against pairing
+#' an analysis with another dataset's components.
+#' @noRd
+components_match <- function(components, data) {
+  if (is.null(components) || is.null(data) || nrow(components) == 0L) {
+    return(FALSE)
+  }
+  sums <- stats::aggregate(cbind(positive, negative) ~ dose, data = components, FUN = sum)
+  observed <- as.data.frame(data)[c("dose", "positive", "negative")]
+  observed <- observed[order(observed$dose), , drop = FALSE]
+  sums <- sums[order(sums$dose), , drop = FALSE]
+  nrow(sums) == nrow(observed) &&
+    isTRUE(all.equal(sums$dose, observed$dose)) &&
+    all(sums$positive == observed$positive) &&
+    all(sums$negative == observed$negative)
+}
+
+#' Replace a singlehit plot's observed points with per-experiment points
+#'
+#' singlehit draws the observed dose groups as one point layer. For pooled
+#' data those points merge experiments at shared doses, so they are swapped
+#' for one point per experiment and dose, coloured by a sub-study field and
+#' sized by the number of subjects. The layer is found by its geometry and
+#' its data; if it is not there the plot is returned unchanged.
+#'
+#' @param aesthetic `"colour"`, or `"fill"` for the overlay, where colour
+#'   already distinguishes the models.
+#' @noRd
+recolour_observed <- function(plot, data, components, field, aesthetic = c("colour", "fill")) {
+  aesthetic <- match.arg(aesthetic)
+  is_observed <- vapply(plot$layers, function(layer) {
+    inherits(layer$geom, "GeomPoint") && is.data.frame(layer$data) &&
+      "dose" %in% names(layer$data) && all(layer$data$dose %in% data$dose)
+  }, logical(1L))
+  if (!any(is_observed)) {
+    return(plot)
+  }
+  plot$layers <- plot$layers[-which(is_observed)[[1L]]]
+  components$group_value <- components[[field]]
+  label <- names(substudy_fields)[substudy_fields == field][[1L]]
+  # No facet column: a layer without `model` is drawn in every model's panel.
+  # Each experiment is nudged by a fixed fraction of a log10 unit, the same in
+  # every panel, so experiments that share a dose and response do not hide
+  # each other.
+  experiments <- unique(components$experiment)
+  offset <- (match(components$experiment, experiments) - (length(experiments) + 1) / 2) * 0.02
+  components$plot_dose <- components$dose * 10^offset
+  points <- if (aesthetic == "colour") {
+    ggplot2::geom_point(
+      data = components,
+      ggplot2::aes(x = .data$plot_dose, y = .data$response, colour = .data$group_value, size = .data$total),
+      alpha = 0.8
+    )
+  } else {
+    ggplot2::geom_point(
+      data = components,
+      ggplot2::aes(x = .data$plot_dose, y = .data$response, fill = .data$group_value, size = .data$total),
+      shape = 21,
+      colour = "grey20",
+      alpha = 0.8
+    )
+  }
+  # In the overlay, colour already shows the models; a separate palette keeps
+  # the experiments' fills from matching a model's line.
+  fill_scale <- if (aesthetic == "fill") ggplot2::scale_fill_viridis_d(end = 0.9)
+  plot + points + fill_scale +
+    ggplot2::scale_size_area(max_size = 5, name = "Subjects", breaks = function(x) unique(round(pretty(x)))) +
+    ggplot2::labs(!!aesthetic := label)
 }
 
 #' Named model choices for a select input