From 87da56de25408bb3fd4c7a8cc88b63cc2b9d01f8 Mon Sep 17 00:00:00 2001 From: seanthimons Date: Wed, 28 Oct 2026 11:09:16 -0400 Subject: [PATCH 1/6] feat(app): add Library tab to browse and facet QMRA wiki experiments Upload the QMRA wiki workbook (or use a 14-experiment sample), filter experiments by normalised metadata, and plot preferred-model curves faceted and coloured by pathogen group, outcome, dose unit, host, route, or model. One fittable experiment can be sent to the Data and Fit tabs with its study id, host, unit, and endpoint. Groundwork for #35. --- .gitignore | 1 + app/DESCRIPTION | 9 +- app/NAMESPACE | 1 + app/R/app_server.R | 5 +- app/R/app_ui.R | 2 + app/R/fct_library.R | 383 +++++++++++++++++++++++++ app/R/mod_library.R | 194 +++++++++++++ app/R/mod_upload.R | 56 ++-- app/README.md | 18 +- app/dev/make_qmra_sample.R | 52 ++++ app/inst/extdata/qmra-wiki-sample.xlsx | Bin 0 -> 14413 bytes app/tests/testthat/helper-data.R | 4 + app/tests/testthat/test-fct_library.R | 118 ++++++++ app/tests/testthat/test-mod_library.R | 48 ++++ app/tests/testthat/test-mod_upload.R | 18 ++ 15 files changed, 885 insertions(+), 24 deletions(-) create mode 100644 app/R/fct_library.R create mode 100644 app/R/mod_library.R create mode 100644 app/dev/make_qmra_sample.R create mode 100644 app/inst/extdata/qmra-wiki-sample.xlsx create mode 100644 app/tests/testthat/test-fct_library.R create mode 100644 app/tests/testthat/test-mod_library.R diff --git a/.gitignore b/.gitignore index 7d8d312..67eb0f4 100644 --- a/.gitignore +++ b/.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/*~ diff --git a/app/DESCRIPTION b/app/DESCRIPTION index fb6e4d4..53d4394 100644 --- a/app/DESCRIPTION +++ b/app/DESCRIPTION @@ -22,16 +22,21 @@ Imports: generics, ggplot2, golem, + grDevices, mirai (>= 2.5.0), parallel, plotly, readr, + readxl, 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 diff --git a/app/NAMESPACE b/app/NAMESPACE index 61cd3a2..b9a2c76 100644 --- a/app/NAMESPACE +++ b/app/NAMESPACE @@ -2,6 +2,7 @@ export(run_app) import(shiny) +importFrom(ggplot2,.data) importFrom(golem, activate_js, add_resource_path, diff --git a/app/R/app_server.R b/app/R/app_server.R index 9dfdb65..960e3af 100644 --- a/app/R/app_server.R +++ b/app/R/app_server.R @@ -5,7 +5,10 @@ #' @import shiny #' @noRd app_server <- function(input, output, session) { - data <- mod_upload_server("upload") + library_sent <- mod_library_server("library") + data <- mod_upload_server("upload", incoming = library_sent) + # A sent experiment replaces the loaded data; show it there first. + observeEvent(library_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 diff --git a/app/R/app_ui.R b/app/R/app_ui.R index fe6840f..8b7201b 100644 --- a/app/R/app_ui.R +++ b/app/R/app_ui.R @@ -12,11 +12,13 @@ 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( "Fit", bslib::layout_sidebar(sidebar = bslib::sidebar(mod_fit_ui("fit")), mod_fit_results_ui("fit")) diff --git a/app/R/fct_library.R b/app/R/fct_library.R new file mode 100644 index 0000000..b9a83bb --- /dev/null +++ b/app/R/fct_library.R @@ -0,0 +1,383 @@ +#' Columns the Library tab needs from the workbook's "raw data" sheet +#' @noRd +library_required_columns <- c( + "pathogen_category", "pathogen", "experiment_id", "Exposure route", "Dose Units", + "Host type", "Preferred model", "Optimized parameters", "Response type", + "Dose", "Positive", "Negative" +) + +#' Fields offered as facets, colours, and filters on the Library tab +#' @noRd +library_fields <- c( + "Pathogen group" = "category", + "Outcome group" = "outcome_group", + "Dose unit" = "dose_unit", + "Host" = "host_group", + "Exposure route" = "route_group", + "Preferred model" = "model_label" +) + +#' Read the QMRA wiki workbook into one row per dose group +#' +#' Reads the "raw data" sheet and adds normalised metadata so that free-text +#' spellings ("mice", "C57BL/6 mice") fall into the same group. Experiments are +#' keyed as in `data-raw/plot-preferred-models.R`: the experiment id alone is +#' not unique, because the same id can be reported for several outcomes. +#' +#' @param path Path to an `.xlsx` workbook. +#' @return A data frame with the original metadata, numeric `dose`, `positive`, +#' and `negative` (`NA` in the single row kept for an experiment that reports +#' only parameters), and the normalised columns `experiment_key`, `category`, +#' `outcome_group`, `dose_unit`, `host_group`, `route_group`, `model`, +#' `model_label`, and `wiki_pooled`. +#' @noRd +read_qmra_workbook <- function(path) { + sheets <- readxl::excel_sheets(path) + if (!"raw data" %in% sheets) { + stop("The workbook has no \"raw data\" sheet.", call. = FALSE) + } + # "minimal": the sheet repeats one percentile column name, which is not used. + raw <- readxl::read_excel(path, sheet = "raw data", .name_repair = "minimal") + missing <- setdiff(library_required_columns, names(raw)) + if (length(missing) > 0L) { + stop(sprintf("The \"raw data\" sheet is missing: %s.", paste(missing, collapse = ", ")), call. = FALSE) + } + raw <- as.data.frame(raw[library_required_columns], stringsAsFactors = FALSE) + rows <- data.frame( + experiment_id = trimws(as.character(raw$experiment_id)), + pathogen = trimws(as.character(raw$pathogen)), + category = library_category(raw$pathogen_category), + response_type = trimws(as.character(raw[["Response type"]])), + outcome_group = library_outcome_group(raw[["Response type"]]), + dose_unit = library_dose_unit(raw[["Dose Units"]]), + host = trimws(as.character(raw[["Host type"]])), + host_group = library_host_group(raw[["Host type"]]), + route = trimws(as.character(raw[["Exposure route"]])), + route_group = library_route_group(raw[["Exposure route"]]), + preferred_model = trimws(as.character(raw[["Preferred model"]])), + model = library_model(raw[["Preferred model"]]), + parameters = as.character(raw[["Optimized parameters"]]), + dose = suppressWarnings(as.numeric(raw$Dose)), + positive = suppressWarnings(as.numeric(raw$Positive)), + negative = suppressWarnings(as.numeric(raw$Negative)), + stringsAsFactors = FALSE + ) + rows <- rows[!is.na(rows$experiment_id) & nzchar(rows$experiment_id), , drop = FALSE] + rows$experiment_key <- paste(rows$experiment_id, rows$pathogen, rows$response_type, rows$preferred_model, sep = "|") + # Some experiments record only the fitted parameters ("NR" counts). Keep one + # row for each of them so their workbook curve can still be drawn. + counted <- stats::complete.cases(rows[c("dose", "positive", "negative")]) + parameter_only <- !rows$experiment_key %in% rows$experiment_key[counted] & !duplicated(rows$experiment_key) + rows <- rows[counted | parameter_only, , drop = FALSE] + if (nrow(rows) == 0L) { + stop("The \"raw data\" sheet has no experiments.", call. = FALSE) + } + rows$model_label <- ifelse(is.na(rows$model), "Other or not reported", model_label_ui(rows$model)) + # Entries such as "142, 143, 144", "297 & 298", or "281,282 (pooled)" were + # already pooled by the wiki. + rows$wiki_pooled <- grepl(",|&|\\band\\b|pooled", rows$experiment_id, ignore.case = TRUE) + rownames(rows) <- NULL + rows +} + +#' @noRd +library_category <- function(x) { + x <- trimws(as.character(x)) + out <- tools::toTitleCase(tolower(x)) + out[is.na(x) | !nzchar(x)] <- "Not reported" + out +} + +#' Collapse free-text response types into three outcome groups +#' +#' Same rules as `data-raw/plot-preferred-models.R`. +#' @noRd +library_outcome_group <- function(x) { + x <- tolower(as.character(x)) + out <- ifelse( + grepl("death|mortality|stillbirth", x), + "Death or mortality", + ifelse(grepl("infection|culture|isolation|shedding|invasion", x), "Infection or detection", "Illness or clinical signs") + ) + out[is.na(x)] <- "Illness or clinical signs" + out +} + +#' @noRd +library_dose_unit <- function(x) { + x <- trimws(as.character(x)) + x[is.na(x) | !nzchar(x) | toupper(x) == "NR"] <- "Not reported" + x +} + +#' @noRd +library_host_group <- function(x) { + x <- tolower(as.character(x)) + out <- ifelse( + grepl("human|volunteer", x), + "Human", + ifelse(grepl("mice|mouse", x), "Mouse", "Other animal") + ) + out[is.na(x) | !nzchar(trimws(x))] <- "Not reported" + out +} + +#' @noRd +library_route_group <- function(x) { + x <- tolower(as.character(x)) + out <- ifelse( + grepl("oral|intragastric|ingest", x), + "Oral", + ifelse( + grepl("nasal|inhal|aerosol|trache|respir", x), + "Respiratory", + ifelse(grepl("intraperitoneal|intravenous|inject|subcutan|intramuscular|intradermal", x), "Injection", "Other") + ) + ) + out[is.na(x) | !nzchar(trimws(x))] <- "Not reported" + out +} + +#' Map the workbook's preferred model to a singlehit model name +#' @return `"exponential"`, `"beta_poisson"`, or `NA`. +#' @noRd +library_model <- function(x) { + x <- tolower(trimws(as.character(x))) + out <- rep(NA_character_, length(x)) + out[x %in% "exponential"] <- "exponential" + out[x %in% "beta-poisson"] <- "beta_poisson" + out +} + +#' Dose groups for one experiment, with the metadata singlehit carries +#' +#' @param rows Output of [read_qmra_workbook()]. +#' @param key One `experiment_key`. +#' @return A data frame accepted by [singlehit::as_dose_response()]. +#' @noRd +library_experiment_data <- function(rows, key) { + experiment <- rows[rows$experiment_key == key, , drop = FALSE] + if (nrow(experiment) == 0L) { + stop("That experiment is not in the workbook.", call. = FALSE) + } + if (anyNA(experiment$dose)) { + stop("The workbook reports only fitted parameters for this experiment, not dose-response counts.", call. = FALSE) + } + data.frame( + dose = experiment$dose, + positive = experiment$positive, + negative = experiment$negative, + study_id = experiment$experiment_id[[1L]], + host = library_or_unknown(experiment$host[[1L]]), + dose_unit = experiment$dose_unit[[1L]], + endpoint = library_or_unknown(experiment$response_type[[1L]]), + stringsAsFactors = FALSE + ) +} + +#' @noRd +library_or_unknown <- function(x) { + if (is.na(x) || !nzchar(x)) "Not reported" else x +} + +#' One row per experiment, with whether singlehit can fit it +#' +#' @param rows Output of [read_qmra_workbook()]. +#' @return A data frame keyed by `experiment_key`. `fittable` is `FALSE`, and +#' `problem` says why, when [singlehit::as_dose_response()] rejects the data. +#' @noRd +library_experiments <- function(rows) { + keys <- unique(rows$experiment_key) + first <- rows[match(keys, rows$experiment_key), , drop = FALSE] + problem <- vapply(keys, function(key) { + checked <- tryCatch(singlehit::as_dose_response(library_experiment_data(rows, key)), error = function(e) e) + if (inherits(checked, "error")) conditionMessage(checked) else NA_character_ + }, character(1L), USE.NAMES = FALSE) + out <- data.frame( + experiment_key = keys, + experiment_id = first$experiment_id, + pathogen = first$pathogen, + category = first$category, + response_type = first$response_type, + outcome_group = first$outcome_group, + dose_unit = first$dose_unit, + host = first$host, + host_group = first$host_group, + route = first$route, + route_group = first$route_group, + model = first$model, + model_label = first$model_label, + wiki_pooled = first$wiki_pooled, + doses = as.integer(tapply(!is.na(rows$dose), rows$experiment_key, sum)[keys]), + subjects = as.numeric(tapply(rows$positive + rows$negative, rows$experiment_key, sum, na.rm = TRUE)[keys]), + fittable = is.na(problem), + problem = problem, + stringsAsFactors = FALSE + ) + rownames(out) <- NULL + out +} + +#' Readable label for an experiment in a select input +#' @noRd +library_experiment_label <- function(experiments) { + paste0( + experiments$pathogen, " (", experiments$response_type, ", ", experiments$dose_unit, + "; experiment ", experiments$experiment_id, ")" + ) +} + +#' Pull a numeric parameter such as `k = 4.19E-03` out of free text +#' @noRd +library_parameter <- function(text, name) { + if (is.na(text)) { + return(NA_real_) + } + pattern <- paste0(name, "[^=]*=\\s*([0-9]+(?:\\.[0-9]+)?(?:[Ee][+-]?[0-9]+)?)") + match <- regmatches(text, regexec(pattern, text, ignore.case = TRUE, perl = TRUE))[[1L]] + if (length(match) < 2L) NA_real_ else as.numeric(match[[2L]]) +} + +#' Preferred-model curves for every experiment in the workbook +#' +#' Refits each experiment's preferred model from its counts and draws it over +#' the observed dose range. When the refit fails, the curve comes from the +#' workbook's own parameters over 0.01 to 100 times the ID50, as in +#' `data-raw/plot-preferred-models.R`. Experiments whose preferred model is not +#' exponential or beta-Poisson, or that have neither, are skipped. +#' +#' @param rows Output of [read_qmra_workbook()]. +#' @param points Points per curve. +#' @return A data frame of `experiment_key`, `dose`, `estimate`, and `source`, +#' with the skipped keys and reasons as the `"skipped"` attribute. +#' @noRd +library_curves <- function(rows, points = 100L) { + keys <- unique(rows$experiment_key) + skipped <- character() + curves <- lapply(keys, function(key) { + experiment <- rows[rows$experiment_key == key, , drop = FALSE] + model <- experiment$model[[1L]] + if (is.na(model)) { + skipped[[key]] <<- paste("unsupported model:", experiment$preferred_model[[1L]]) + return(NULL) + } + fit <- tryCatch( + singlehit::fit_dose_response(library_experiment_data(rows, key)[c("dose", "positive", "negative")], model), + error = function(e) e + ) + if (!inherits(fit, "error") && identical(fit$convergence, 0L)) { + curve <- singlehit::prediction_curve(fit, points = points, dose_range = range(fit$data$dose)) + return(data.frame(experiment_key = key, dose = curve$dose, estimate = curve$estimate, source = "Refit from counts")) + } + text <- experiment$parameters[[1L]] + k <- library_parameter(text, "k") + alpha <- library_parameter(text, "(?:a|\u03b1)") + n50 <- library_parameter(text, "N50") + median_dose <- if (model == "exponential") log(2) / k else n50 + if (!is.finite(median_dose) || median_dose <= 0 || (model == "beta_poisson" && (!is.finite(alpha) || alpha <= 0))) { + skipped[[key]] <<- "no usable fit or workbook parameters" + return(NULL) + } + dose <- exp(seq(log(median_dose / 100), log(median_dose * 100), length.out = points)) + estimate <- if (model == "exponential") { + singlehit::exponential_response(dose, k) + } else { + singlehit::beta_poisson_response(dose, alpha, n50) + } + data.frame(experiment_key = key, dose = dose, estimate = estimate, source = "Workbook parameters") + }) + out <- do.call(rbind, curves) + if (is.null(out)) { + out <- data.frame(experiment_key = character(), dose = numeric(), estimate = numeric(), source = character()) + } + attr(out, "skipped") <- skipped + out +} + +#' Colours for the pathogen groups, with a fallback for anything else +#' @noRd +library_palette <- function(values) { + known <- c(Bacteria = "#C65D00", Virus = "#0066A1", Protozoa = "#008C6A", Prion = "#6A1B9A") + values <- sort(unique(values)) + extra <- setdiff(values, names(known)) + fallback <- grDevices::hcl.colors(max(length(extra), 1L), "Dark 3")[seq_along(extra)] + c(known, stats::setNames(fallback, extra))[values] +} + +#' Overview plot of preferred-model curves +#' +#' @param curves Output of [library_curves()]. +#' @param experiments Output of [library_experiments()], already filtered. +#' @param facet One of [library_fields] or `"none"`. +#' @param colour One of [library_fields]. +#' @param free_x Whether each facet gets its own dose axis. +#' @return A ggplot object. +#' @importFrom ggplot2 .data +#' @noRd +library_plot <- function(curves, experiments, facet = "none", colour = "category", free_x = FALSE) { + if (!colour %in% library_fields) { + stop("Choose a field to colour by.", call. = FALSE) + } + if (!facet %in% c("none", library_fields)) { + stop("Choose a field to facet by.", call. = FALSE) + } + metadata <- experiments[c("experiment_key", unique(library_fields))] + data <- merge(curves, metadata, by = "experiment_key") + if (nrow(data) == 0L) { + stop("No curves to plot for these filters.", call. = FALSE) + } + data$colour <- data[[colour]] + counts <- table(data$colour[!duplicated(data$experiment_key)]) + groups <- names(counts) + labels <- stats::setNames(paste0(groups, " (", as.integer(counts), ")"), groups) + colours <- if (colour == "category") library_palette(groups) else NULL + colour_scale <- if (is.null(colours)) { + ggplot2::scale_colour_discrete(labels = labels) + } else { + ggplot2::scale_colour_manual(values = colours, breaks = groups, labels = labels) + } + + same_units <- facet == "dose_unit" || length(unique(data$dose_unit)) == 1L + x_label <- if (same_units && length(unique(data$dose_unit)) == 1L) { + sprintf("Dose (%s, log10 scale)", data$dose_unit[[1L]]) + } else if (same_units) { + "Dose (log10 scale; unit in panel title)" + } else { + "Recorded dose (log10 scale; study-specific units)" + } + caption <- paste( + "Curves use observed dose ranges when refit from counts; otherwise 0.01-100 x ID50 from workbook parameters (dashed).", + if (same_units) "" else "Dose units differ between curves, so horizontal positions are not comparable; facet by dose unit to compare." + ) + + plot <- ggplot2::ggplot( + data, + ggplot2::aes( + x = .data$dose, + y = .data$estimate, + group = .data$experiment_key, + colour = .data$colour, + linetype = .data$source + ) + ) + + ggplot2::geom_line(alpha = 0.6, linewidth = 0.55) + + colour_scale + + ggplot2::scale_linetype_manual(values = c("Refit from counts" = "solid", "Workbook parameters" = "dashed")) + + ggplot2::guides(colour = ggplot2::guide_legend(override.aes = list(alpha = 1))) + + ggplot2::scale_x_log10() + + ggplot2::scale_y_continuous(limits = c(0, 1)) + + ggplot2::labs( + subtitle = sprintf("%d experiments", length(unique(data$experiment_key))), + x = x_label, + y = "Probability of reported outcome", + colour = names(library_fields)[library_fields == colour][[1L]], + linetype = "Curve source", + caption = paste(strwrap(caption, width = 120), collapse = "\n") + ) + + ggplot2::theme_minimal(base_size = 12) + + ggplot2::theme(legend.position = "right", plot.caption = ggplot2::element_text(hjust = 0)) + if (facet != "none") { + plot <- plot + ggplot2::facet_wrap(ggplot2::vars(.data[[facet]]), scales = if (free_x) "free_x" else "fixed") + } + plot +} diff --git a/app/R/mod_library.R b/app/R/mod_library.R new file mode 100644 index 0000000..2e89ca1 --- /dev/null +++ b/app/R/mod_library.R @@ -0,0 +1,194 @@ +#' Library module UI +#' +#' Load the QMRA wiki workbook, filter its experiments, draw their preferred +#' models faceted by metadata, and send one experiment to the Data and Fit tabs. +#' +#' @param id Module id. +#' @noRd +mod_library_ui <- function(id) { + ns <- NS(id) + filter <- function(field, label) { + selectizeInput(ns(paste0("filter_", field)), label, choices = character(), multiple = TRUE, options = list(placeholder = "All")) + } + bslib::layout_sidebar( + sidebar = bslib::sidebar( + width = 320, + fileInput(ns("file"), "QMRA wiki workbook", accept = ".xlsx"), + actionButton(ns("sample"), "Use sample workbook"), + p( + class = "text-muted small mt-1", + "14 experiments cut down from the QMRA wiki workbook, for demos and testing. Data compiled from the ", + a("QMRA wiki", href = "https://qmrawiki.org", target = "_blank"), + "; see each row's reference." + ), + hr(), + h6("Filter experiments"), + filter("category", "Pathogen group"), + filter("pathogen", "Pathogen"), + filter("outcome_group", "Outcome group"), + filter("dose_unit", "Dose unit"), + filter("host_group", "Host"), + filter("route_group", "Exposure route"), + filter("model_label", "Preferred model"), + checkboxInput(ns("hide_wiki_pooled"), "Hide entries the wiki already pooled", value = FALSE), + hr(), + h6("Plot"), + selectInput(ns("facet"), "Facet by", choices = c("None" = "none", library_fields), selected = "outcome_group"), + selectInput(ns("colour"), "Colour by", choices = library_fields, selected = "category"), + checkboxInput(ns("free_x"), "Separate dose axis per facet", value = FALSE) + ), + uiOutput(ns("error")), + uiOutput(ns("empty")), + uiOutput(ns("status")), + plotOutput(ns("plot"), height = "620px"), + uiOutput(ns("send_panel")), + tableOutput(ns("table")) + ) +} + +#' Library module server +#' +#' @param id Module id. +#' @return A reactive that is `NULL` until an experiment is sent, then a list +#' of `data`, `label`, and `sent` (a counter that changes on every send). +#' @noRd +mod_library_server <- function(id) { + moduleServer(id, function(input, output, session) { + filters <- c("category", "pathogen", "outcome_group", "dose_unit", "host_group", "route_group", "model_label") + workbook <- reactiveVal(NULL) + + observeEvent(input$file, workbook(input$file$datapath)) + observeEvent(input$sample, workbook(app_sys("extdata", "qmra-wiki-sample.xlsx"))) + + rows <- reactive({ + req(workbook()) + tryCatch(read_qmra_workbook(workbook()), error = function(e) e) + }) + loaded <- reactive({ + req(!inherits(rows(), "error")) + rows() + }) + # Checking and refitting every experiment is the slow step (seconds for the + # full workbook), so it runs once per workbook, not on every filter change. + experiments <- reactive(library_experiments(loaded())) + curves <- reactive(library_curves(loaded())) + + observeEvent(experiments(), { + for (field in filters) { + updateSelectizeInput(session, paste0("filter_", field), choices = sort(unique(experiments()[[field]])), selected = character()) + } + }) + + filtered <- reactive({ + out <- experiments() + for (field in filters) { + chosen <- input[[paste0("filter_", field)]] + if (length(chosen) > 0L) { + out <- out[out[[field]] %in% chosen, , drop = FALSE] + } + } + if (isTRUE(input$hide_wiki_pooled)) { + out <- out[!out$wiki_pooled, , drop = FALSE] + } + out + }) + + output$error <- renderUI({ + if (inherits(rows(), "error")) { + div(class = "alert alert-danger", conditionMessage(rows())) + } + }) + output$empty <- renderUI({ + if (is.null(workbook())) { + p(class = "text-muted", "Upload the QMRA wiki workbook (.xlsx) or use the sample workbook to browse its experiments.") + } + }) + output$status <- renderUI({ + shown <- filtered() + drawn <- length(intersect(shown$experiment_key, curves()$experiment_key)) + p( + class = "text-muted", + sprintf( + "%d of %d experiments match the filters; %d have a curve to draw.", + nrow(shown), nrow(experiments()), drawn + ) + ) + }) + + output$plot <- renderPlot( + { + shown <- filtered() + validate(need(any(curves()$experiment_key %in% shown$experiment_key), "No curves to draw for these filters.")) + library_plot(curves(), shown, facet = input$facet, colour = input$colour, free_x = isTRUE(input$free_x)) + }, + res = 96 + ) + + output$send_panel <- renderUI({ + shown <- filtered() + req(nrow(shown) > 0L) + choices <- stats::setNames(shown$experiment_key, library_experiment_label(shown)) + keep <- intersect(isolate(input$experiment), choices) + div( + class = "d-flex flex-wrap align-items-end gap-3 mt-3", + selectInput( + session$ns("experiment"), + "Experiment to analyse", + choices = choices, + selected = if (length(keep) > 0L) keep else NULL, + width = "520px" + ), + uiOutput(session$ns("send_button")) + ) + }) + + chosen <- reactive({ + req(input$experiment) + match <- experiments()[experiments()$experiment_key == input$experiment, , drop = FALSE] + req(nrow(match) == 1L) + match + }) + + output$send_button <- renderUI({ + experiment <- chosen() + if (experiment$fittable) { + actionButton(session$ns("send"), "Send to Data and Fit", icon = icon("arrow-right"), class = "btn-primary mb-3") + } else { + div(class = "alert alert-warning small py-2 mb-3", "Cannot be fitted: ", experiment$problem) + } + }) + + output$table <- renderTable({ + shown <- filtered() + req(nrow(shown) > 0L) + data.frame( + Experiment = shown$experiment_id, + Pathogen = shown$pathogen, + Outcome = shown$response_type, + `Dose unit` = shown$dose_unit, + Host = shown$host, + Route = shown$route, + `Preferred model` = shown$model_label, + Doses = shown$doses, + Subjects = format(shown$subjects, big.mark = ","), + `Wiki-pooled` = ifelse(shown$wiki_pooled, "yes", ""), + Fittable = ifelse(shown$fittable, "yes", "no"), + check.names = FALSE + ) + }) + + sent <- reactiveVal(NULL) + count <- reactiveVal(0L) + observeEvent(input$send, { + experiment <- chosen() + req(experiment$fittable) + count(count() + 1L) + sent(list( + data = library_experiment_data(loaded(), experiment$experiment_key), + label = library_experiment_label(experiment), + sent = count() + )) + }) + sent + }) +} diff --git a/app/R/mod_upload.R b/app/R/mod_upload.R index 98fcb6c..7f1c7c0 100644 --- a/app/R/mod_upload.R +++ b/app/R/mod_upload.R @@ -37,6 +37,7 @@ mod_upload_table_ui <- function(id) { ns <- NS(id) tagList( uiOutput(ns("error")), + uiOutput(ns("origin")), uiOutput(ns("summary")), tableOutput(ns("table")) ) @@ -44,29 +45,34 @@ mod_upload_table_ui <- function(id) { #' Upload module server #' -#' Reads the chosen file with [singlehit::read_dose_response()], which -#' standardizes and validates the grouped counts. +#' Reads the chosen file with [singlehit::read_dose_response()], or a data +#' frame sent from the Library tab with [singlehit::as_dose_response()]; both +#' standardize and validate the grouped counts. #' #' @param id Module id. +#' @param incoming Reactive returning `NULL` or a list with `data` (a data frame +#' of dose groups), `label` (where it came from), and `sent` (a counter, so +#' sending the same experiment again reloads it). #' @return A reactive returning the standardized data, or `NULL` when nothing #' valid has been loaded. #' @noRd -mod_upload_server <- function(id) { +mod_upload_server <- function(id, incoming = reactive(NULL)) { moduleServer(id, function(input, output, session) { fields <- c("dose", "positive", "negative") no_mapping <- list(dose = NULL, positive = NULL, negative = NULL) - source_path <- reactiveVal(NULL) + # A path to read, or a list(data, label) sent from the Library tab. + loaded <- reactiveVal(NULL) mapping <- reactiveVal(no_mapping) # A new source starts from auto-detection. The mapping is reset here, in the - # same step as the path, so the new file is never parsed with the previous - # file's column choices. - load_source <- function(path) { + # same step as the source, so the new data are never parsed with the + # previous source's column choices. + load_source <- function(value) { mapping(no_mapping) for (field in fields) { updateSelectInput(session, field, selected = "") } - source_path(path) + loaded(value) } observeEvent(input$file, load_source(input$file$datapath)) @@ -74,13 +80,20 @@ mod_upload_server <- function(id) { input$example, load_source(system.file("extdata", "Ward_rotavirus.txt", package = "singlehit")) ) + observeEvent(incoming(), { + sent <- incoming() + load_source(list(data = sent$data, label = sent$label)) + }) # Same reader settings as singlehit::read_dose_response(), so the offered # names match what the real read sees. columns <- reactive({ - req(source_path()) + req(loaded()) + if (is.list(loaded())) { + return(names(loaded()$data)) + } tryCatch( - names(readr::read_delim(source_path(), n_max = 0L, show_col_types = FALSE, progress = FALSE, trim_ws = TRUE)), + names(readr::read_delim(loaded(), n_max = 0L, show_col_types = FALSE, progress = FALSE, trim_ws = TRUE)), error = function(e) character() ) }) @@ -100,17 +113,20 @@ mod_upload_server <- function(id) { }) parsed <- reactive({ - req(source_path()) + req(loaded()) chosen <- mapping() - tryCatch( - singlehit::read_dose_response( - source_path(), - dose = chosen$dose, - positive = chosen$positive, - negative = chosen$negative - ), - error = function(e) e - ) + read <- if (is.list(loaded())) { + function() singlehit::as_dose_response(loaded()$data, chosen$dose, chosen$positive, chosen$negative) + } else { + function() singlehit::read_dose_response(loaded(), dose = chosen$dose, positive = chosen$positive, negative = chosen$negative) + } + tryCatch(read(), error = function(e) e) + }) + + output$origin <- renderUI({ + if (is.list(loaded())) { + p(class = "text-muted", "From the Library tab: ", strong(loaded()$label)) + } }) output$error <- renderUI({ diff --git a/app/README.md b/app/README.md index 5c94388..0f3a5c4 100644 --- a/app/README.md +++ b/app/README.md @@ -36,7 +36,23 @@ depends on the app. "Use {plotly} to render plot" under Options for hover and zoom. The interface font is overridden in `inst/app/www/custom.css`. -Pooling and export are planned. +- **Library:** upload the QMRA wiki workbook (`.xlsx`, "raw data" sheet) or + use the bundled 14-experiment sample, filter experiments by pathogen group, + pathogen, outcome group, dose unit, host, exposure route, and preferred model, + and draw every experiment's preferred model on one log-dose plot, faceted and + coloured by any of those fields. Curves are refit from the counts; experiments + the workbook reports only as parameters are drawn from those (dashed). Dose + units differ between studies, so the plot says so unless it is faceted by + dose unit or filtered to one. "Send to Data and Fit" loads one experiment, + with its study id, host, unit, and endpoint, as the current data; experiments + singlehit cannot fit (fewer than three doses, parameters only) say why + instead. Free-text host, route, and outcome spellings are grouped by simple + rules in `R/fct_library.R`. The full workbook is not in the repository; the + sample is rebuilt from it with `dev/make_qmra_sample.R`. The data come from + the [QMRA wiki](https://qmrawiki.org) and the studies it cites; each sample + row keeps its reference, DOI, and wiki link. + +Pooling (#35) and export are planned. ## License diff --git a/app/dev/make_qmra_sample.R b/app/dev/make_qmra_sample.R new file mode 100644 index 0000000..9c3110e --- /dev/null +++ b/app/dev/make_qmra_sample.R @@ -0,0 +1,52 @@ +# Build the small QMRA wiki sample workbook used by the Library tab's "Use +# sample workbook" button and by the tests. The full workbook stays out of the +# repository; this keeps only the "raw data" sheet, only count and metadata +# columns, and a handful of experiments chosen to cover each pathogen group, +# several dose units, both preferred models, a wiki-pooled entry, a +# parameter-only entry, one with too few doses to fit, and one with an +# unsupported model. +# +# Run from the app/ directory: +# Rscript dev/make_qmra_sample.R "../QMRA wiki data_v4.xlsx" +args <- commandArgs(trailingOnly = TRUE) +if (length(args) != 1L) { + stop("Supply the path to the full QMRA wiki workbook.") +} + +keys <- c( + "18|Burkholderia mallei|infection|exponential", + "21|Burkholderia mallei|death|beta-Poisson", + "23|Burkholderia mallei|death|beta-Poisson", + "21,23|Burkholderia mallei|death|beta-Poisson", + "1|Yersinia pestis|death|exponential", + "Ojielo2003|Pseudomonas aeruginosa|death|exponential", + "71|Rotavirus|infection|beta-Poisson", + "31|Adenovirus|infection|exponential", + "63|Enterovirus|infection|exponential", + "1|Variola major|death|NR", + "108|Cryptosporidium hominis/ Cryptosporidium parvum|infection|exponential", + "46|Giardia|infection|exponential", + "48|Giardia|infection|exponential", + "252|scrapie|infection|exponential" +) +columns <- c( + "pathogen_category", "pathogen", "experiment_id", "Exposure route", "# of Doses", + "Agent Strain", "Dose Units", "Host type", "Preferred model", "LD50/ID50", + "Optimized parameters", "Response type", "Reference", "Dose", "Positive", + "Negative", "Total", "ref_doi", "wiki_link" +) + +raw <- readxl::read_excel(args[[1L]], sheet = "raw data", .name_repair = "minimal") +raw <- as.data.frame(raw[columns], check.names = FALSE) +key <- paste( + trimws(raw$experiment_id), trimws(raw$pathogen), trimws(raw[["Response type"]]), + trimws(raw[["Preferred model"]]), + sep = "|" +) +missing <- setdiff(keys, key) +if (length(missing) > 0L) { + stop("Not in the workbook: ", paste(missing, collapse = "; ")) +} +sample <- raw[key %in% keys, , drop = FALSE] +writexl::write_xlsx(list(`raw data` = sample), "inst/extdata/qmra-wiki-sample.xlsx") +cat(sprintf("Wrote %d rows from %d experiments.\n", nrow(sample), length(keys))) diff --git a/app/inst/extdata/qmra-wiki-sample.xlsx b/app/inst/extdata/qmra-wiki-sample.xlsx new file mode 100644 index 0000000000000000000000000000000000000000..5a70fc8eb6e58e34ec9419751e78374055de7392 GIT binary patch literal 14413 zcmaib1y~$Qvo?}IaCc8|C%9{H4FLj+6LfKRcMI+koZt?@34w&*?oM!bSmYma&UerG z^4$CUz*Kct*IQj(T{FEqO)Ja6KEsBBf_ed56i@~_GJJ~AP*5MAK|x^w5|Fr!wWG1M zqux6=TVn?uW>+iA@`PbU7&i3YBcYx`jatyLex$+bN_JjcD4O!TL1WGVM<#S~BXnZF7uY=R*QCG~pg zs3nPCO>c`xL2>GOTS7kgI2n=;Oi^sW_k3uI(J2W!D^uIDct z>K!u#IgxYF;yPo0kh&Z6(+7o*YL?F!#+2%NtDK_bNU!jrOmdAmt!Kt~ARG)X`zkE& zA-=7FtNAQ@mn)DUMmRAZ7@Hqn&9h%aLp`*m9w6Fwj^55!LoH8oW@bT45K^Tj zkKb~bjPtL%&pvRoH)J5Ua?BKNMXq;v%XcnER;aW0H0vzN{}ukzY?--jAnbB5P*6w! z*RwaabYNkAqUB>@F0_uNdh;lbq*s_?_z-MhF6*paw`QoAc6Ff9)j5DC(Vb=q3F zT#I8Bz7#nb7p;g?x=$P#v{NfmL(mZ%88n$Aj6AVRbq7&H!u+3vci(x;i!RM`7O%!^ z$61p~nQmtdfeX9gbL%`^;}Nf2v2f%k=E>8pi+qcXbH1N>|{s`CEgt zXx&=5*fJvN$S$e3`l%rrbbarU)51CQi35(_qNaV8{*oNP1;3c5o%pME;0(36lCP&- zE#}JVd2TA)os!EgpQ&x~0Gs{>5!NaHFIY_pxY6x}R}9w?(pgz0_yO>4M73I&-MFy2 zo#kgnMt?~0V~9r$h+>$_zC9iPrrKClEn6cpbK z6coUo;`0<42Qy=1M+cVQjO{5zQJD%3#NQL&T%9qX&wTVP`C0ik#__#B6-Gjkq@Bqn z>nt__H3A$I=&~c5+y7vb-$N+H&vuc)3^v z88=nZ?SR9+w;apC6AFPZ}XoK0<({t_lGa%j1v#X#_Az0@zdYAxqUnys}>(F$4?e3ZOF2+6mOPI zTbkUO+%7i`FB{`qSOS)Pnr~X)k2j9HyW2nh%GSLtheo)1@# z1ARBK9Q-0_e4=r1frV$x;(Yly_r;$nF5|OC*v2M!$6}*~@zUmg2=D!fj+>fu^NF3v zzF_dO-*1I4zxa5ZV$E0j34@9R|C>1>Iwyq+7Uo~%A}lSY7~xM6=jfEX(RacCRtU%~ zoKqOHKa^1L2ql!>F_+9Wp|p_9RhC!sVlNkm18RuB&I6sbq);UWrLhq$%26Z)6(bld zyzFslWG(LV+6x~v^~N_O1_5~JuDNI~E~SNNuClt47w8F&ONkhC5!m|#XOhO21q1M$ zhAA4FR!JE`-MnJw__}N3gdsPH2@zc-gfLm@cVfB9`bu5_PgqRMAWRRji1JtPPsNuY zFsw?=9rH8l*JJQcW(Je$uxaQIC@je9iXtV9Dd}IS(AF2VPfk(5st|;URSV|KPR0H% zekrf_qt$#kCQOVij+9V2pf*R)*mOfvf$bB-p)4U+c^)P<7X9$?`Fj(L_C#n!b+syy~$PfR`P;8c}tHO1mh_NTLE|r z6fw1bd+SFl)9r*EZdwbQ70RJV=D@6fsltF+e;K9~ zTZ{t`T8>{@M4Q{3e!R=XPv3TTM`qIL)v2^357W- zr>$~8#Z%I>g2pk|#Boz=Pr`gxVsLw)ekv@kj#q`9&p| z9Co00Ah}etX&pBK6+r`X#LO>%yVR6_@qD7ij-b@hR5`-5RE~>jTA=sL$=Bo-6F&ws z7q$n0!4FDLEYKik-Td@?5w=Wj7rXK_+@ z6PW-!|9D}>pqP+9`P;R)3dIQx(?jl3=~4GX2^l_KmF)pj$4$P|eF>Rn9-vas#jAp( zY!|d+K+WJ?UcbsIYl35q8{;_Q!HSk~lMz!KC_zO@PoeTz2}n60y?%$`uS^?ynmfaW zZZJ0A89^rRlrypyvff`y<0eX)-YoOo28OItXF@tro5A~Y9bfaeUL;5Uonya zva_q=s8Gahc3`lMPU>Bvqq1w2q!IQmH`kX_+~;QQsB`^*3HRV>{X)a2iLg6}>6gP^ z_auHruPtXV5>+z}{F9JyF&tJqY88`!bP18JEs5hRpwNI+;U}A|5ePPJQV&MrAKI!% zoBnBAjtr0-rN%X=C$0kUWa+cxLfV5B@4ZF&ulIyt?o9O$;KlB#*D<^Oywy}Y8}l5V zHq~~v6n3ovt;^xoMstRDX4T~zO71^x$-dx*@b__|Ud4fl52>XZt$T1EbE}ux?dw13 zBy)J`1?-s6G-!NG#vO)jg7fxDMXv~3O~Ys2m({z~+2Jzr1sgV1u{`;X)y>JKT|8UP z(^=qAeyvE4Db++B7%?`eNW7Bf+{LMU&>v|#h&(-DYY=cdW~1E52CK~u#9~(+XeX#h zT@_om?tx|wyCh75@ zWkjC-YRiO-cT@vNFd^S6DN(1j<0P9Ta^L}7@a*$xkxmxc6{>^IM0{pTklwA;9^f(q zn>Ot#;AbL!#W&lWkf#8CdyrX|T|2%%6T^Mf8oT}pG}QsYkI)Q11(3w;B*D^m8P{+k)5HiK}9w z9p{j_r*4g#s%oFYoVZ$ojECsXfNCqVKdFI4cO8(Q6DB$Sq_7Kg>+dZp%CfToQu1V} z7A`?z$Re)rj%ZPHvzhwpzf{u39|7M{Hv_3?kMf4iLTI4D3*i}+eAwKDNWE0u-q zU^Um^W2TsjI=yM|E*as)H)xjhTquW4YKy|$!S}&(x z;bhlkU;pj-aT-Bo>PZ>_i}}qJ*h70q{3VG`dA&82rIufX8=elVX14N_++0qCRufAZ z&q)3Gj#f%V&`#XU(W1xD)5y1f*R@iZ*$ykHFbS>#rfv)HrmVTjMW$;aDn=i+&Gx7c zM&9r{>%F)8#Bi}AE3aYLj5+?h@9n=URANG{oIhnN74f)+j1El35giM398ny-#B^B9 zLzw4AAz#B`b-pNRm^8m0PfTSdz8hN;`A@|>@M<`!d2F>lqCuxaaeL1S@Vhml!mmA8Kn2N8Ci9$L-JM{7z2Fm$BnziY|CpLoYSaYpG!W zOk~#d-6c?#d;)E?*9`yGjC!T}#3BCyTWv1ohs$~W*xsi?N&f#p@6YZXNbg<)j{I$6 z$>dPuDzy~?ogE4$;cMkLNxF1Pt{a9m&sLc;O}Zo~@G>#ahKvc-)(tZra7KGpKXGU(QJ4PEcdlgJP@R{e{RdKe%S1Q9R!J-5~{5q_5x?m&Y3m<;cA_E z;;H3cCQqXIfqS_edBcYonE^ecX5%TK--hA5Y{55H*_!jF6PKy=UWZ`Fz?@KR!!*xY z&l4XmJ79kj0ZuQQ$q}EtWcP{MNt&JC(peHhYxTlq+2eoiZCbs2YGVf_zup7AXX9t2 zwuEd|%`FZhhU(U_eeH#D@-_Gp^x(98EZ3s7-*^Wi-p(zzM%CKy7<3#wF{ySD2MKHk z55YQQGn}cVR{r!<+bq4jDH+f(Zl^yBYFYh&xxr$&-HCTO2eMm{Gkl~gWXOEqbUkMH zI5)axdr?q*q-lN*n z+>FoE27PVwv>SXoY0(3bQF#+Ep(oj4u z$4(@dZ`EtQz;ZxDqk)^vuRCKo+>_qB`kl4Ry;FlA1@h8xs?qRa$gl%v@b~rVmr}ve zVfoKBGA1qw@0uBMS7y7rRCgw>7rO5)(eKKPLY7OC4VTp(@FMO**nrhZFRfXt$FUY` zUKO~dyQo)<^$+{Mxhh0XKM{@P#s7KCxnwu4`2y0_ULGytld+ja4N-^6&mH`@r?=^u zzyt8$D;R(dxDG~6KaszoCaosDNR7Eq&;i%gslT8d`1Bi_XlZb*X6qTia(tZTvdW2m zWo*!j+Y5dorIm-k1G1h~7im2RG6u5i1f;eO#7xny=KD;gK?`@#<9lLLKDAiztw}%- zzj?_-Y8bX+t|#JGwl_DH`2Ma0rB47z1~%O$!Dw8mR*uWzXiTYATHuq+8lYE)RO#2m zl6vwpGWnZ<&B@o@c2Bjdo|7?AjqeQXx|ooqjbTnyi2VDD4Nr3ce!a58a(BB<4pa@j zbEZDgaa%hjl$-R-86WNLtdAJ_uugi2u4A4HXXH321~5TzzVQWDlJe%KOkfHRlUBPaFL)Rxjopazba*H`ccqeWcl5=>gR**I#fz8G;67IC!=%c{t2JH`hl;E0EFdM zTZVY)CFez>v;qlhdNQoK3&9;9>vbZEw{ZZ=$$rO-UCty%)sCK0EmmQkJN$d&Ae+hd zf-rTduXS_$os5^LtR9_yRbdF@Gs2SU3D;@2-+#!4`_xMPBEfb_{#FSarl+~4k#vpe z@?stCjt}icdMazfwK~a*brK+8=v?du3x1V$l=4Q8#B_Go#B}D@#F~AlHkZXP1}y^K zB(8TqNjqx3CpU@XH?B+@i_Clttp3rsD}rz1U^D4g6@;)(>KXOQfUn6XrYwIQSYf8C zzJGRMwRT<;nPE=UgQxGP`B_kX@-#N@ zsT)j4ql`veaC#FJiXNT}q2AZmunqm2(|&`MmS+s}o?-mH)q|DBQ11;o-){jk9qH7r z1`i}%FJDK^P~(j8dnziQ7>*{yvc=!M?%eT>+$1`A8)rWCI(1Wz`l3`GAM(4`;iq1W zzs_uFT;5=M44)?EpHM|!HgqI$vytN`fCqlaK73;`+We^qQ&3}YtlE-Hd71gfB z|2M^bDtH^^EzfLobKIgtd5h=S&^UjcW6oQKlkBJ-wXr-WBD(LD>_Cb9TZ`-Ve*jO# zc9g@TILwe~R+(-v{**P)Hg9}B-?`2GyYXA?1=Pl+x0b9gNjsZ@9{oLLnXcbigKnmGthcaRo&OPdyIiO?A3GPw+@cxPS>*k z#xVr@4UVgWQ&Qg0N47ngXJ26YT?=F}w!EV{Id8l|OrX#gA36ewoC)dYNSQ z{V#C6&DvE2v+?x{86>X}4sCy#T^FgjI!`^;2za9GI=?;}hVlB>I!!d*pwxD0yLGK)bou|Q7?IDIpuUe7ZOC7PJ?S!F9a`&PxDn>@9~Zm*nI+1 zIeAuRz0%opOJ?RT7_YDCp5_@T2;kRowkYa1ECL$`RZjXbhjrWNl$#ybPqDODIoF!= z>fnCdCd9vUGHj+UMB@IYM^gyi9p_}z*+1fpc3dswbKVu1{^nTSKD{fQjPL4WIuV5G zn9~N_J%u>iX)zVGkjC%p-93ZdgS9l-*4}~-H=BpkE@VeaOJvIe0)m2%g;`k(z?=1h zBOmXk`}1E)H|HBHjc90%WZsYS^Xr3q>x)|*Ohq&{);55of3}v!v}1AU{-WP4>t6W3 zcOIt2dm?0kokCh*X9n$`ornMI$wbA-DFNFtfk$`YlS}5L7O%7>->gx^pr0%MbPMB- z3L=X=h&f-*%a`>dD!Pa?5i)OWKa2`yxfJDPU`;Q0VJSzGDEPv2UQkmJjP+FFSgSA!j71vWYsRYS7inUx$@DeY zJW?YS>%ze{>FMe_qKsqDsngGXl*9Mo)bkx$aAE1_Ck>yq;>U?wR5l?A7l?{IR+Kqr zc6w&LB&JDZLT^5~ulP{*+~iHOfXBx&R442v#&AbpnvU)DM+Llno*`BExR(@>hw<1c z6xWGyF!%8TlD+R3Ws*T?-2?$Qp0P|ACF*;iv0tTL&KxkoV|B)7DE}Syw}G-G_$W|N zY;;giPusS@hL3}pfxYnuRY!YsYg31(&7K}D9eW~fOz&k?M{}_1TjS81YzmToy0Ba; z-lX!l1hcaZ7*i3WQ_diPwy$Ed@~4(ow4 zGT|wl1-61O*dZPi2$#NLVoBNqXjk6t?!5Sieq8So>f-!K4w+JyASXS`SzxObBOL^< zrbgC6w@pa%`>HgHmqgS;Q`>$Rzi!l+O&pi2ehWba)=hKh#sY+gKx{ zJVp>7K>fX$wBhS*c*C4NCkw4+#hO2OvB(fDX0%IwIV`uj9vZE$qt4r3RvUy*%cr2V zE?8%HTIO5hz*2#_8$ge!_HZ7*IH%-4dJA- z6}(5NXVk1;yt|}Cxf_q?n6Y+U#@9Vt)RD^73j45w1h7SrQ#Mwj~-7q~W+u7}Dm&QpxFz+;YmL)~SHphl1x0V}F zfSH>v4O((?b}PyoLqfW2L6nIt^Qvx2qF1`L2|aIHTw2Y*bEW}rSq++6#8L|8b3r_T zfq{8_-?U7~ zyV@%1%L_Ecq`4gjce0(V@N01kcM`eJEBd0wI>ZA4PCNekiH#_w2D){hy z4ReOyi@sXmr?}Q*=447Me(IqpNf)k4rBP%&PbF`!d2Dmr80ORcbW@`ZY-U* z$Eb~DSG-;sF`jid$yYay(vsWb(TFT5v{CPBQ{nNp!?%siYffue8WkB>!x42HqWlO2 z3+8F@g3Y@}JxE91&iH=%VvyFkL+6)>`YHFKZVj8!xhrFtjTr5@l_QAFaUd^^PFJ~l zC$=)#VW>>J@D*Ba=Akr_R!th7L<0)lhl6pVO;ux8m1^kBTuG$$qF2JF-71@CQg%Fd zD%Eie>iFvI?`nhugxz+QxDf(W_WM7wfBuEm~TC z6sR(@C#19#AAlMA3L(==N=dKi`+PgwSGFOWXxo_8^jDIC@5lD~^s;H(ZeQSaWG;eV zfwdSr>r_MPCd%8@CY(-ed3I-oo1NrqRzINgT(PR)&sYYwQF1+FNh$eYx7P)-mk@9j ze^VOK!*sENoF4IJo`W=l^|TOGsSJdqQ|ZseA=7aYwIeaUzR@w zyIAGVejw5di@se*G3<+^S%!&jQUqUb8Ov#j4E6as z0|`b-LDf!ut|6>xT{aviGV?^Nk1jbeg9RZUVLxQm4_+l}%je;wl9sISL$Y{s&6zPr z04w2b!X%}FG(IM=V3m;y+Wb06y3bnR{pH=V0-t>pzw9GKlLTJcL+2>0mQQ6bx^!wV zLYDf2joiLvM8xw0d7GB!>^QbnCl(9IFhP`&WI?)Hg(Fhzm>cn;EBgkF^M(f9%yCO> zoC5)+(PY^_oOM;Kh(;;&*r=Y@THXw%k%KDxYhYpG32F?^(bHwuK1vNykxJS3ICkT3 zzGYgCX3!GtHmS$X*sbCOjpzc?#N?nZ#F>hH(vNN!8|dXJKT_$lR`C?(jqH@hH|I9WUcL zhpiUKX7qD~K!cKHGxm^I+;3+mo(<95kV|I@T;BZ}9zqv#36QZF!;2T#k-P21Y1H=+ zTk#4xqtIF?6f8dzOJQEkB5bI^y!>UGv}hI6ocuP)H{%P@w%0=fB`I?4@v!g%=fIe5 zT%)YDwaMG0LHBLG;T%;CsJL$iyYP0nDE?;G!=i?xA5wyZace$=u+J7lzg0~p=XE&~ z=I8bNwQ}qubj7_0UZ-tgZ$&aUtHju1-feSh{AsPdq^_ol*eAYmlf&fw@F!S%S*vf7<( zCmiQT!C?WTRTs?aEe*vH!>g7DJg{#&BQ7yksm8M*K3K6TZ>zJy`rehXI;uD-=9Og5 zOQ8TP#^OosvYXMJpHo^^{rO++E+Lq+pHbJQBN@Px;akBsYdV;Vh`-u z7}(y(V8zAdc;tzn^1wx}^>wNo>t*hzUlu#3YKg$Whr_OgYWo@Syh^S*KRHq%tGx*_ zB+geJC;F^CWq{q2hfkLRiuQ0>r9`r_iwI&bcI$$393n^7UyYQ5#Uu{y4*bwPBd$ub zY9Qi2@r!tZ>zPWDYt;<|K{X+g&=kq^-Olr#xaH?=Z_2E^R~e->zpCJ3Hd#7#5q1}~ zYN;KPW}+)xSr>?SeE;~#A;6eX5ML7VV?I4NaEF}u-gwNLovJ%DuvDd~oF;w6A|hak z-?@@eago8n1I5Gr&3%2LS!12jHC`nBj;6qYp1>lNV{yo?!F027-KyZ|bk6lb38~3V z6IC7eeb3f%Zzn#oOZWr6H$gRC7J$15ovmCEGku3 zP*mJ=29WXR2^M^krcV|Kcx5pq62u?R?qA5c@FRXu;4jydkzo5J@s1P2Q*`rmh-_1a zOJ1@)eGz$Fs!1pVIRlJHl7D2^N{ebQ934D5QvS|!eBS+lmV9-DN$5+<{1@%V!-ZXJ zyX2c&Z0px-l$upHjaas~G>$M+b={&xg`Nu~^Fd`J+N=-na2#F-j?8AF$oN}u6W*WS z&aGNA(-9TjoPdM_gtC81TFYR_%!Yo5KhoiB{Z$VC{aK8%jJY{hO}U=AQR*}pQZw-u z&3n$VMGqcnKC#mq``IY8AG>F7U1gD{i5~vR4w(NOAHmd$g0p-jc8PFs?yd8jheX~~ zhm|dVIang8;Oa(|R9|gw)KRHkAIDGeqoNFHOlZF7Yw7AKjHbZzbZfT7nTt>JmCH&) z@mb~HDi{&mZ{MgW$=g{U@)7dl#JlxB&l`#NmHZs4+0XsiU1H4fv<+ft*ISA> zUb7c%k*!Lx4_XI-d5zv2HQn5=E;Cpj@R5}(z>yIZm21t}M?Bws22Mtk&_={(H)a** zeRw;L&r-ZgtVEna z%pOfzl;=Z7L4A=?*=}Z#1_FJLSY;RyKQ7bkp<*f6yj-kc#jZ%ObQQ(GjMtU^hRr(Q z8O3HB13&xyesr7sg}(qXkLabVR1E(Su`FVV-|NYtFEMjF7Z1U;vX`qRQ7+fmP1`b^ zLmwB4CnkGy&0gCwR)TuGCw$-dS0^+MNYboz5#YQjEl{h!hw%@`Ow~0Pg-?8JZ~nEP zh1W4pA)mzuaJm%3`e6=dqFT{Vgj`o@}2-)7LZ`RTTrq z!~@+>?g%6^lF|*vf$sZNXxzrHmNAqnHSr1FZW4K!a!be7Ss_Ir;Q?`D?(5TIC0|In%~HP7!@Ahl&mD>WgvRJlCjM=~wSTa$QJf^n5n6>^!?m zJI$5-&GOcl%6E0kgHF@|4UkSt4KgC^aPFCp>LMCskQpbW-s3Zo`xu3`yi#4MODrgm zMiV1X;k!MdYXxG4Q8Wb_Z7+)zQT!ECAvd7`ID@Q@UrB@^q({S6cx64Cvt$d$t zQ{8h~tJcxQbSUqOV@~H1^V^N$hM|tdK?Xg~RtDlq!P6fxqftZ0wb*$J$nM>GZHH8ZO#|0oYXpJ0wizaeA{2DzKq5c zilc|ltCYs^i1inzUJLYy(z5E2^F@e!MoHP}BKZmzljFy`wTsWEq%;T^lBQS`8Xg8i zW^x-dW*u4g+hK?HL8!j)*D+zzq^Lh)yFvrFbShx9`-Nvf9h-xJO`0#keHUj+bQ*{f zkSJ}%p6*orwq>ohX5R5-+H;mQ_cq-|jqII%X!Kx-goc?J^$wjQc_*z;YxXKGjhoY2 zYOMtwUwrir9QhtBo-b#tgOIMBtqST#i&sDntsQ1Qrg`%^pI7dAy^W!lJDQuvO4c0E zf64E&e5&;vdd(gXN!hUz*a)uI@@R(Iq#f4M32J5SpVERs%ionVINk6!fJv2+l|+w& za-krgha-S4!>BsMR?X2(`PmeMiX*4sE37E0k92z5gI1|0AP!uD-sgAMqr=}nl937D z_EETENmwGa-yWWIw7Z2qcwQcq(&#mx%=JEAt`~JYoCXn@S+4X~W<5OI?etd)L)IZQ zsMo>%IVfaRfo6{B#JM`DzP1=KIegjQd|^tZp%M{zd$O{E=$Fle~oiie9b@-UpuNpg*lVd%waMWQcqF^wvm{nL_YPEbg zESBk60ClGFAIy8ky^H-><%MUDlwTK{Zr>qCA(P3)E^4(e<*p}-Y9>Wkgi)hKG5g*g z&d|L1h|O?FW6DKSei4bXgP&~;#JKNbLp>I=Xtk)i+wM<3Kq>~8B(^PtAW6Hrp!hYCfom!`2CRi+Cy{6qa}(emN5WfqRY&@=!26+Jz3>;>=JSWqL(@;82Gz<0Nj)22;{2yh#i` zUzzl;hBVE?W69+D<;ONT(6T%1A;*Yjuhe-@b0r#oXHv2JFvcCRiLy`KO?9 z4H`phpKcptyC!^i;>B%Cb&mN9pHal44;Fo95@NLGlp|BRei?VnInRhiGqn;Zvc$C@e=YI>B;4(jV_wZY*F zr^ZvQ?KjMSPpmJFL&0BoeunaM7X@;|rv~gRr+(pM@($VXJr-4OV@Q3(A2sKugN=Er{d2ZZ-jTqopQT;$yB-n`shFD!oHgEaV~g%J=kYc$_g$Ohg%IS_8pIpN2%Yx@?)j@b>?&Wna$+J(E>i;{z?Bf4=%5euGLUDnJ9mphI& zWR}Yq;Ky@ij|i$=P54yLy5sNh#6sjbvy`v!9r7X?WuzoUU2c=RlDL+%9q-4M|NG%y zIfSY64R|~G7C6TMye<3vn#|G7()joLv9@ZGRYd0;>0O?y%@fg)U4re0d`;=c^_w!xS+BrAl* zzUmUsNj+QoahRK0vIt4{zme7=6}Kd>$Yk*LtCKL3F(a^=-irT9ZXZL-lVfb8fY;h3 z@Amn9WWF&t7j0#P>w%vKvDoI98jgoRP#k296<3sY> zm4w^{Kb-g1P!)|P>v=c`NMwy1FK(C`L>PKaCd4phDMit~3PeiUzoqA;?^fan5RpNI zV{8hnW0Z9_z>C`R0(5m^3svN9WzeNkF!JUOoz2PZ2d>^sI0;Za&XrNEZ%7DCm>c7T}vM-^YcH zALWfPAo!K1mgIV1>b&M~6ycZuJz3d=EKY~7-M7#0&}PN4ovL+6xA-WNEPdX3ui>i9 zY|98FI~=uMHafd;UC(n%3i`dz(|bkQ`fP)25Vhtmz?GLqh`pG=D0uk|3e3@qN}o2F zW#wdY)T!;yx0~x&!uhaPqQ`tA)48*aG1W_P>2N{_o8HnX-P_A`=pWzf2_6xx-;%M= zSSOdh(MZ zgGhe@|2(4e4=@Js@c%xr^QWCZkFxw@rwT~G|Lb7OpH}|dmiosEJ@Nlo`OhZRpWr_? zSN;L#ko>iV{JGWgr-eViRQzLMkN&rXKfhx9Y2eR=>mLI@SpRb2f0nU7p?^LE|A98X o{$C!6f5QK~-T#3jbNmbb_sy>?_Z;p?0}=T31?IVt^Qrd#0Gvon7ytkO literal 0 HcmV?d00001 diff --git a/app/tests/testthat/helper-data.R b/app/tests/testthat/helper-data.R index 9a3ca66..c62e540 100644 --- a/app/tests/testthat/helper-data.R +++ b/app/tests/testthat/helper-data.R @@ -22,3 +22,7 @@ example_analysis_with_bands <- function(models = c("exponential", "beta_poisson" boots <- singlehit::bootstrap_dose_response_models(analysis$fits, times = 20, seed = 1, backend = "sequential") singlehit::add_bootstrap(analysis, boots) } + +qmra_sample_path <- function() { + system.file("extdata", "qmra-wiki-sample.xlsx", package = "singlehitapp") +} diff --git a/app/tests/testthat/test-fct_library.R b/app/tests/testthat/test-fct_library.R new file mode 100644 index 0000000..911b515 --- /dev/null +++ b/app/tests/testthat/test-fct_library.R @@ -0,0 +1,118 @@ +test_that("free-text metadata falls into shared groups", { + expect_equal( + library_host_group(c("human", "Human, immune", "mice", "C57Bl/6J Mice", "guinea pig", NA)), + c("Human", "Human", "Mouse", "Mouse", "Other animal", "Not reported") + ) + expect_equal( + library_route_group(c("oral (in milk)", "intragastric", "intranasal", "aerosol", "intraperitoneal", "contact lens")), + c("Oral", "Oral", "Respiratory", "Respiratory", "Injection", "Other") + ) + expect_equal( + library_outcome_group(c("death", "stillbirths", "infection", "positive stool culture", "diarrhea")), + c("Death or mortality", "Death or mortality", "Infection or detection", "Infection or detection", "Illness or clinical signs") + ) + expect_equal(library_dose_unit(c(" CFU ", "NR", NA, "")), c("CFU", "Not reported", "Not reported", "Not reported")) + expect_equal(library_category(c("bacteria", "Bacteria", "virus", NA)), c("Bacteria", "Bacteria", "Virus", "Not reported")) + expect_equal(library_model(c("exponential", "beta-Poisson", "NR", NA)), c("exponential", "beta_poisson", NA, NA)) +}) + +test_that("workbook parameters are pulled out of free text", { + expect_equal(library_parameter("k = 4.19E-03", "k"), 4.19e-3) + expect_equal(library_parameter("a = 0.201 k = 4.87E-9 N50 = 1.42E+8", "N50"), 1.42e8) + expect_equal(library_parameter("a = 0.04 N50 = 8.26E+08", "(?:a|\u03b1)"), 0.04) + expect_true(is.na(library_parameter("NR", "k"))) + expect_true(is.na(library_parameter(NA, "k"))) +}) + +test_that("the sample workbook reads into one row per dose group", { + rows <- read_qmra_workbook(qmra_sample_path()) + expect_true(all(c("experiment_key", "dose", "positive", "negative", "dose_unit", "outcome_group") %in% names(rows))) + expect_equal(length(unique(rows$experiment_key)), 14L) + # The parameter-only experiment keeps a single row without counts. + parameter_only <- rows[rows$experiment_id == "Ojielo2003", ] + expect_equal(nrow(parameter_only), 1L) + expect_true(is.na(parameter_only$dose)) + expect_true(rows$wiki_pooled[rows$experiment_id == "21,23"][[1L]]) + expect_false(any(rows$wiki_pooled[rows$experiment_id == "21"])) +}) + +test_that("a workbook without the raw data sheet is rejected", { + skip_if_not_installed("writexl") + path <- withr::local_tempfile(fileext = ".xlsx") + writexl::write_xlsx(list(other = data.frame(x = 1)), path) + expect_error(read_qmra_workbook(path), "raw data") +}) + +test_that("experiments report whether singlehit can fit them", { + rows <- read_qmra_workbook(qmra_sample_path()) + experiments <- library_experiments(rows) + expect_equal(nrow(experiments), 14L) + expect_false(experiments$fittable[experiments$experiment_id == "Ojielo2003"]) + expect_match(experiments$problem[experiments$experiment_id == "Ojielo2003"], "only fitted parameters") + expect_false(experiments$fittable[experiments$pathogen == "Enterovirus"]) + expect_match(experiments$problem[experiments$pathogen == "Enterovirus"], "three distinct dose groups") + expect_true(experiments$fittable[experiments$pathogen == "Rotavirus"]) + expect_equal(experiments$doses[experiments$experiment_id == "Ojielo2003"], 0L) +}) + +test_that("experiment data carry singlehit's metadata columns", { + rows <- read_qmra_workbook(qmra_sample_path()) + key <- rows$experiment_key[rows$pathogen == "Rotavirus"][[1L]] + data <- library_experiment_data(rows, key) + standardized <- singlehit::as_dose_response(data) + expect_equal(unique(standardized$dose_unit), "FFU") + expect_equal(unique(standardized$study_id), "71") + expect_error(library_experiment_data(rows, "nope"), "not in the workbook") +}) + +test_that("curves are refit from counts or drawn from workbook parameters", { + rows <- read_qmra_workbook(qmra_sample_path()) + curves <- library_curves(rows, points = 20L) + sources <- tapply(curves$source, curves$experiment_key, unique) + expect_equal(unname(sources[grepl("^Ojielo2003", names(sources))]), "Workbook parameters") + expect_equal(unname(sources[grepl("^71\\|", names(sources))]), "Refit from counts") + expect_true(all(curves$estimate >= 0 & curves$estimate <= 1)) + expect_equal(names(attr(curves, "skipped")), "1|Variola major|death|NR") +}) + +test_that("the overview plot builds for every facet and colour", { + rows <- read_qmra_workbook(qmra_sample_path()) + experiments <- library_experiments(rows) + curves <- library_curves(rows, points = 20L) + for (facet in c("none", library_fields)) { + plot <- library_plot(curves, experiments, facet = facet, colour = "category") + expect_s3_class(plot, "ggplot") + expect_no_error(ggplot2::ggplot_build(plot)) + } + for (colour in library_fields) { + expect_no_error(ggplot2::ggplot_build(library_plot(curves, experiments, colour = colour, free_x = TRUE))) + } +}) + +test_that("the dose axis says when units differ", { + rows <- read_qmra_workbook(qmra_sample_path()) + experiments <- library_experiments(rows) + curves <- library_curves(rows, points = 20L) + mixed <- library_plot(curves, experiments) + expect_match(mixed$labels$x, "study-specific units") + expect_match(mixed$labels$caption, "facet by dose unit") + cfu <- library_plot(curves, experiments[experiments$dose_unit == "CFU", ]) + expect_equal(cfu$labels$x, "Dose (CFU, log10 scale)") + by_unit <- library_plot(curves, experiments, facet = "dose_unit") + expect_match(by_unit$labels$x, "unit in panel title") +}) + +test_that("an unknown pathogen group still gets a colour", { + colours <- library_palette(c("Virus", "Fungus")) + expect_equal(names(colours), c("Fungus", "Virus")) + expect_false(anyNA(colours)) +}) + +test_that("plots with nothing to draw or bad fields fail clearly", { + rows <- read_qmra_workbook(qmra_sample_path()) + experiments <- library_experiments(rows) + curves <- library_curves(rows, points = 20L) + expect_error(library_plot(curves, experiments[0, ]), "No curves") + expect_error(library_plot(curves, experiments, facet = "pathogen"), "facet") + expect_error(library_plot(curves, experiments, colour = "nope"), "colour") +}) diff --git a/app/tests/testthat/test-mod_library.R b/app/tests/testthat/test-mod_library.R new file mode 100644 index 0000000..5737ae2 --- /dev/null +++ b/app/tests/testthat/test-mod_library.R @@ -0,0 +1,48 @@ +test_that("the sample workbook loads and the filters narrow the experiments", { + shiny::testServer(mod_library_server, { + session$setInputs(sample = 1) + expect_equal(nrow(experiments()), 14L) + expect_equal(nrow(filtered()), 14L) + + session$setInputs(filter_category = "Protozoa") + expect_equal(nrow(filtered()), 3L) + session$setInputs(filter_category = character(), filter_dose_unit = c("CFU", "PFU")) + expect_setequal(unique(filtered()$dose_unit), c("CFU", "PFU")) + session$setInputs(filter_dose_unit = character(), hide_wiki_pooled = TRUE) + expect_false(any(filtered()$wiki_pooled)) + }) +}) + +test_that("sending an experiment returns its data with a changing counter", { + shiny::testServer(mod_library_server, { + session$setInputs(sample = 1) + key <- experiments()$experiment_key[experiments()$pathogen == "Rotavirus"] + session$setInputs(experiment = key, send = 1) + first <- session$returned() + expect_equal(first$sent, 1L) + expect_match(first$label, "Rotavirus") + expect_equal(unique(first$data$dose_unit), "FFU") + + session$setInputs(send = 2) + expect_equal(session$returned()$sent, 2L) + }) +}) + +test_that("an experiment that cannot be fitted is not sent", { + shiny::testServer(mod_library_server, { + session$setInputs(sample = 1) + key <- experiments()$experiment_key[experiments()$experiment_id == "Ojielo2003"] + session$setInputs(experiment = key, send = 1) + expect_null(session$returned()) + }) +}) + +test_that("a file that is not a QMRA workbook shows an error", { + skip_if_not_installed("writexl") + path <- withr::local_tempfile(fileext = ".xlsx") + writexl::write_xlsx(list(other = data.frame(x = 1)), path) + shiny::testServer(mod_library_server, { + session$setInputs(file = list(datapath = path)) + expect_s3_class(rows(), "error") + }) +}) diff --git a/app/tests/testthat/test-mod_upload.R b/app/tests/testthat/test-mod_upload.R index 1118733..efa0158 100644 --- a/app/tests/testthat/test-mod_upload.R +++ b/app/tests/testthat/test-mod_upload.R @@ -56,3 +56,21 @@ test_that("the example button loads the Ward rotavirus data", { expect_equal(sum(data$total), 59) }) }) + +test_that("data sent from the Library tab load without a file", { + sent <- shiny::reactiveVal(NULL) + data <- data.frame(dose = c(1, 10, 100), positive = c(1, 5, 9), negative = c(9, 5, 1), study_id = "x", dose_unit = "CFU") + shiny::testServer(mod_upload_server, args = list(incoming = sent), { + sent(list(data = data, label = "Test experiment", sent = 1L)) + session$flushReact() + result <- session$returned() + expect_equal(result$dose, c(1, 10, 100)) + expect_equal(unique(result$dose_unit), "CFU") + expect_true(is.list(loaded())) + + # Loading the example afterwards replaces the library data. + session$setInputs(example = 1) + expect_false(is.list(loaded())) + expect_equal(nrow(session$returned()), 8L) + }) +}) From 84c036c8b8739a77c238127eece76eb64a3b85f1 Mon Sep 17 00:00:00 2001 From: seanthimons Date: Wed, 28 Oct 2026 13:22:13 -0400 Subject: [PATCH 2/6] feat(app): load workbooks on the Data tab and fix Library plot layout The Data tab is now the single place to load data: delimited files and the Ward example become the dataset Fit analyses, while an .xlsx QMRA wiki workbook (or the bundled sample) is loaded for the Library tab, which keeps exploration and send-to-Fit. A bad workbook reports its error and keeps the previous one. The Library plot no longer collapses in the fillable layout ("figure margins too large"), and narrow facets get power-of-ten dose labels and wrapped panel titles. --- app/DESCRIPTION | 1 + app/R/app_server.R | 17 +++- app/R/fct_library.R | 11 ++- app/R/mod_library.R | 45 ++++------ app/R/mod_upload.R | 116 ++++++++++++++++++++------ app/README.md | 16 ++-- app/tests/testthat/test-app_server.R | 10 +++ app/tests/testthat/test-mod_library.R | 28 ++++--- app/tests/testthat/test-mod_upload.R | 71 +++++++++++++--- 9 files changed, 222 insertions(+), 93 deletions(-) create mode 100644 app/tests/testthat/test-app_server.R diff --git a/app/DESCRIPTION b/app/DESCRIPTION index 53d4394..412e0ae 100644 --- a/app/DESCRIPTION +++ b/app/DESCRIPTION @@ -28,6 +28,7 @@ Imports: plotly, readr, readxl, + scales, shiny, singlehit (>= 0.3.0), stats, diff --git a/app/R/app_server.R b/app/R/app_server.R index 960e3af..49d37cf 100644 --- a/app/R/app_server.R +++ b/app/R/app_server.R @@ -5,10 +5,19 @@ #' @import shiny #' @noRd app_server <- function(input, output, session) { - library_sent <- mod_library_server("library") - data <- mod_upload_server("upload", incoming = library_sent) - # A sent experiment replaces the loaded data; show it there first. - observeEvent(library_sent(), bslib::nav_select("main", "Data", session = session)) + # 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_sent <- mod_library_server("library", upload$workbook) + data <- upload$data + observeEvent(library_sent(), { + sent(library_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) analysis <- mod_fit_server("fit", data) mod_compare_server("compare", analysis) # Plots and Explore get the analysis with the bootstraps attached once a run diff --git a/app/R/fct_library.R b/app/R/fct_library.R index b9a83bb..276808f 100644 --- a/app/R/fct_library.R +++ b/app/R/fct_library.R @@ -364,7 +364,8 @@ library_plot <- function(curves, experiments, facet = "none", colour = "category colour_scale + ggplot2::scale_linetype_manual(values = c("Refit from counts" = "solid", "Workbook parameters" = "dashed")) + ggplot2::guides(colour = ggplot2::guide_legend(override.aes = list(alpha = 1))) + - ggplot2::scale_x_log10() + + # Powers of ten ("10^5") stay short enough not to collide in narrow facets. + ggplot2::scale_x_log10(breaks = scales::breaks_log(n = 4), labels = scales::label_log()) + ggplot2::scale_y_continuous(limits = c(0, 1)) + ggplot2::labs( subtitle = sprintf("%d experiments", length(unique(data$experiment_key))), @@ -372,12 +373,16 @@ library_plot <- function(curves, experiments, facet = "none", colour = "category y = "Probability of reported outcome", colour = names(library_fields)[library_fields == colour][[1L]], linetype = "Curve source", - caption = paste(strwrap(caption, width = 120), collapse = "\n") + caption = paste(strwrap(caption, width = 90), collapse = "\n") ) + ggplot2::theme_minimal(base_size = 12) + ggplot2::theme(legend.position = "right", plot.caption = ggplot2::element_text(hjust = 0)) if (facet != "none") { - plot <- plot + ggplot2::facet_wrap(ggplot2::vars(.data[[facet]]), scales = if (free_x) "free_x" else "fixed") + plot <- plot + ggplot2::facet_wrap( + ggplot2::vars(.data[[facet]]), + scales = if (free_x) "free_x" else "fixed", + labeller = ggplot2::label_wrap_gen(width = 20) + ) } plot } diff --git a/app/R/mod_library.R b/app/R/mod_library.R index 2e89ca1..82b55ec 100644 --- a/app/R/mod_library.R +++ b/app/R/mod_library.R @@ -1,7 +1,8 @@ #' Library module UI #' -#' Load the QMRA wiki workbook, filter its experiments, draw their preferred -#' models faceted by metadata, and send one experiment to the Data and Fit tabs. +#' Explore the QMRA wiki workbook loaded on the Data tab: filter its +#' experiments, draw their preferred models faceted by metadata, and send one +#' experiment back to the Data and Fit tabs. #' #' @param id Module id. #' @noRd @@ -10,18 +11,13 @@ mod_library_ui <- function(id) { filter <- function(field, label) { selectizeInput(ns(paste0("filter_", field)), label, choices = character(), multiple = TRUE, options = list(placeholder = "All")) } + # Not fillable: the plot sits above a long table, and in a fillable layout + # the page squeezes the plot to fit the window until base graphics fails + # with "figure margins too large". Here the main panel scrolls instead. bslib::layout_sidebar( + fillable = FALSE, sidebar = bslib::sidebar( width = 320, - fileInput(ns("file"), "QMRA wiki workbook", accept = ".xlsx"), - actionButton(ns("sample"), "Use sample workbook"), - p( - class = "text-muted small mt-1", - "14 experiments cut down from the QMRA wiki workbook, for demos and testing. Data compiled from the ", - a("QMRA wiki", href = "https://qmrawiki.org", target = "_blank"), - "; see each row's reference." - ), - hr(), h6("Filter experiments"), filter("category", "Pathogen group"), filter("pathogen", "Pathogen"), @@ -37,10 +33,9 @@ mod_library_ui <- function(id) { selectInput(ns("colour"), "Colour by", choices = library_fields, selected = "category"), checkboxInput(ns("free_x"), "Separate dose axis per facet", value = FALSE) ), - uiOutput(ns("error")), uiOutput(ns("empty")), uiOutput(ns("status")), - plotOutput(ns("plot"), height = "620px"), + plotOutput(ns("plot"), height = "620px", fill = FALSE), uiOutput(ns("send_panel")), tableOutput(ns("table")) ) @@ -49,23 +44,16 @@ mod_library_ui <- function(id) { #' Library module server #' #' @param id Module id. +#' @param rows Reactive returning the workbook rows from [read_qmra_workbook()], +#' or `NULL` before a workbook is loaded on the Data tab. #' @return A reactive that is `NULL` until an experiment is sent, then a list #' of `data`, `label`, and `sent` (a counter that changes on every send). #' @noRd -mod_library_server <- function(id) { +mod_library_server <- function(id, rows) { moduleServer(id, function(input, output, session) { filters <- c("category", "pathogen", "outcome_group", "dose_unit", "host_group", "route_group", "model_label") - workbook <- reactiveVal(NULL) - - observeEvent(input$file, workbook(input$file$datapath)) - observeEvent(input$sample, workbook(app_sys("extdata", "qmra-wiki-sample.xlsx"))) - - rows <- reactive({ - req(workbook()) - tryCatch(read_qmra_workbook(workbook()), error = function(e) e) - }) loaded <- reactive({ - req(!inherits(rows(), "error")) + req(rows()) rows() }) # Checking and refitting every experiment is the slow step (seconds for the @@ -93,14 +81,9 @@ mod_library_server <- function(id) { out }) - output$error <- renderUI({ - if (inherits(rows(), "error")) { - div(class = "alert alert-danger", conditionMessage(rows())) - } - }) output$empty <- renderUI({ - if (is.null(workbook())) { - p(class = "text-muted", "Upload the QMRA wiki workbook (.xlsx) or use the sample workbook to browse its experiments.") + if (is.null(rows())) { + p(class = "text-muted", "Load the QMRA wiki workbook (.xlsx), or the sample workbook, on the Data tab to browse its experiments.") } }) output$status <- renderUI({ diff --git a/app/R/mod_upload.R b/app/R/mod_upload.R index 7f1c7c0..3413d6d 100644 --- a/app/R/mod_upload.R +++ b/app/R/mod_upload.R @@ -1,7 +1,7 @@ #' Upload module UI #' -#' Sidebar inputs for choosing a delimited file or the bundled example, with an -#' optional column-mapping panel. +#' Sidebar inputs for choosing a delimited file, the QMRA wiki workbook, or the +#' bundled examples, with an optional column-mapping panel. #' #' @param id Module id. #' @noRd @@ -9,16 +9,28 @@ mod_upload_ui <- function(id) { ns <- NS(id) auto <- c("Auto-detect" = "") tagList( - fileInput(ns("file"), "Dose-response file", accept = c(".csv", ".tsv", ".txt")), + fileInput(ns("file"), "Dose-response file or QMRA wiki workbook", accept = c(".csv", ".tsv", ".txt", ".xlsx")), + p( + class = "text-muted small", + "A delimited file is one experiment and goes straight to Fit. An .xlsx workbook is browsed on the Library tab,", + "where you choose an experiment to send here." + ), actionButton(ns("example"), "Use Ward rotavirus example data"), p(class = "text-muted small mt-1", "Ward et al. (1986) human rotavirus challenge study."), + actionButton(ns("sample_workbook"), "Use QMRA wiki sample workbook"), + p( + class = "text-muted small mt-1", + "14 experiments cut down from the QMRA wiki workbook, for demos and testing. Data compiled from the ", + a("QMRA wiki", href = "https://qmrawiki.org", target = "_blank"), + "; see each row's reference." + ), bslib::accordion( id = ns("mapping_panel"), open = FALSE, class = "mt-3", bslib::accordion_panel( "Column mapping (optional)", - p(class = "text-muted small", "Only needed when column names are not recognized automatically."), + p(class = "text-muted small", "Only needed when a delimited file's column names are not recognized automatically."), selectInput(ns("dose"), "Dose column", choices = auto), selectInput(ns("positive"), "Positive-response column", choices = auto), selectInput(ns("negative"), "Negative-response column", choices = auto) @@ -29,13 +41,15 @@ mod_upload_ui <- function(id) { #' Upload module table UI #' -#' Main-panel output showing the standardized data or the import error. +#' Main-panel output showing the loaded workbook, then the standardized data +#' or the import error. #' #' @param id Module id. #' @noRd mod_upload_table_ui <- function(id) { ns <- NS(id) tagList( + uiOutput(ns("workbook_status")), uiOutput(ns("error")), uiOutput(ns("origin")), uiOutput(ns("summary")), @@ -45,29 +59,52 @@ mod_upload_table_ui <- function(id) { #' Upload module server #' -#' Reads the chosen file with [singlehit::read_dose_response()], or a data -#' frame sent from the Library tab with [singlehit::as_dose_response()]; both -#' standardize and validate the grouped counts. +#' Holds two independent things: the current dataset, which Fit analyses, and +#' the QMRA wiki workbook, which the Library tab explores. A delimited file is +#' read with [singlehit::read_dose_response()] and a data frame sent from the +#' Library tab with [singlehit::as_dose_response()]; both standardize and +#' validate the grouped counts. An `.xlsx` file is read as a workbook with +#' [read_qmra_workbook()] and leaves the current dataset alone. #' #' @param id Module id. #' @param incoming Reactive returning `NULL` or a list with `data` (a data frame #' of dose groups), `label` (where it came from), and `sent` (a counter, so #' sending the same experiment again reloads it). -#' @return A reactive returning the standardized data, or `NULL` when nothing -#' valid has been loaded. +#' @return A list of reactives: `data`, the standardized data or `NULL`; +#' `workbook`, the workbook rows or `NULL`; and `workbook_loaded`, a counter +#' that changes each time a workbook loads successfully. #' @noRd mod_upload_server <- function(id, incoming = reactive(NULL)) { moduleServer(id, function(input, output, session) { fields <- c("dose", "positive", "negative") no_mapping <- list(dose = NULL, positive = NULL, negative = NULL) - # A path to read, or a list(data, label) sent from the Library tab. + # list(path, label) for a file, or list(data, label) sent from the Library tab. loaded <- reactiveVal(NULL) mapping <- reactiveVal(no_mapping) + # A failed workbook read keeps the previous workbook and reports the error. + workbook <- reactiveVal(NULL) + workbook_name <- reactiveVal(NULL) + workbook_error <- reactiveVal(NULL) + workbook_loaded <- reactiveVal(0L) + + load_workbook <- function(path, name) { + rows <- tryCatch(read_qmra_workbook(path), error = function(e) e) + if (inherits(rows, "error")) { + workbook_error(sprintf("Could not read %s: %s", name, conditionMessage(rows))) + return(invisible()) + } + workbook_error(NULL) + workbook(rows) + workbook_name(name) + workbook_loaded(workbook_loaded() + 1L) + } # A new source starts from auto-detection. The mapping is reset here, in the # same step as the source, so the new data are never parsed with the # previous source's column choices. load_source <- function(value) { + # A failed workbook read is stale once something else loads. + workbook_error(NULL) mapping(no_mapping) for (field in fields) { updateSelectInput(session, field, selected = "") @@ -75,25 +112,40 @@ mod_upload_server <- function(id, incoming = reactive(NULL)) { loaded(value) } - observeEvent(input$file, load_source(input$file$datapath)) + observeEvent(input$file, { + name <- input$file$name + if (is.null(name)) name <- basename(input$file$datapath) + if (tolower(tools::file_ext(name)) == "xlsx") { + load_workbook(input$file$datapath, name) + } else { + load_source(list(path = input$file$datapath, label = name)) + } + }) observeEvent( input$example, - load_source(system.file("extdata", "Ward_rotavirus.txt", package = "singlehit")) + load_source(list( + path = system.file("extdata", "Ward_rotavirus.txt", package = "singlehit"), + label = "Ward rotavirus example data" + )) + ) + observeEvent( + input$sample_workbook, + load_workbook(app_sys("extdata", "qmra-wiki-sample.xlsx"), "the QMRA wiki sample workbook") ) observeEvent(incoming(), { sent <- incoming() - load_source(list(data = sent$data, label = sent$label)) + load_source(list(data = sent$data, label = paste(sent$label, "from the Library tab"))) }) # Same reader settings as singlehit::read_dose_response(), so the offered # names match what the real read sees. columns <- reactive({ req(loaded()) - if (is.list(loaded())) { + if (!is.null(loaded()$data)) { return(names(loaded()$data)) } tryCatch( - names(readr::read_delim(loaded(), n_max = 0L, show_col_types = FALSE, progress = FALSE, trim_ws = TRUE)), + names(readr::read_delim(loaded()$path, n_max = 0L, show_col_types = FALSE, progress = FALSE, trim_ws = TRUE)), error = function(e) character() ) }) @@ -115,18 +167,32 @@ mod_upload_server <- function(id, incoming = reactive(NULL)) { parsed <- reactive({ req(loaded()) chosen <- mapping() - read <- if (is.list(loaded())) { + read <- if (!is.null(loaded()$data)) { function() singlehit::as_dose_response(loaded()$data, chosen$dose, chosen$positive, chosen$negative) } else { - function() singlehit::read_dose_response(loaded(), dose = chosen$dose, positive = chosen$positive, negative = chosen$negative) + function() singlehit::read_dose_response(loaded()$path, dose = chosen$dose, positive = chosen$positive, negative = chosen$negative) } tryCatch(read(), error = function(e) e) }) output$origin <- renderUI({ - if (is.list(loaded())) { - p(class = "text-muted", "From the Library tab: ", strong(loaded()$label)) - } + req(loaded()) + h6(class = "mt-2", "Current dataset: ", loaded()$label) + }) + + output$workbook_status <- renderUI({ + tagList( + if (!is.null(workbook_error())) div(class = "alert alert-danger", workbook_error()), + if (!is.null(workbook())) { + div( + class = "alert alert-info py-2", + sprintf( + "Workbook loaded: %d experiments from %s. Explore them and pick one to analyse on the Library tab.", + length(unique(workbook()$experiment_key)), workbook_name() + ) + ) + } + ) }) output$error <- renderUI({ @@ -144,8 +210,10 @@ mod_upload_server <- function(id, incoming = reactive(NULL)) { data_preview(parsed()) }) - reactive({ - if (inherits(parsed(), "error")) NULL else parsed() - }) + list( + data = reactive(if (inherits(parsed(), "error")) NULL else parsed()), + workbook = reactive(workbook()), + workbook_loaded = reactive(workbook_loaded()) + ) }) } diff --git a/app/README.md b/app/README.md index 0f3a5c4..6a8a9cb 100644 --- a/app/README.md +++ b/app/README.md @@ -8,8 +8,11 @@ depends on the app. ## What it does today -- **Data:** upload a delimited file or use the bundled Ward rotavirus example, with optional - column mapping when names are not recognized automatically. +- **Data:** the one place to load data. A delimited file, or the bundled Ward + rotavirus example, becomes the dataset Fit analyses, with optional column + mapping when names are not recognized automatically. An `.xlsx` QMRA wiki + workbook, or the bundled sample workbook, is loaded for the Library tab + instead and leaves the current dataset alone. - **Fit:** choose at least two models (exponential, beta-Poisson, exact beta-Poisson), see the dose-trend test, parameter estimates, and goodness of fit. The analysis runs without bootstrapping. @@ -36,15 +39,16 @@ depends on the app. "Use {plotly} to render plot" under Options for hover and zoom. The interface font is overridden in `inst/app/www/custom.css`. -- **Library:** upload the QMRA wiki workbook (`.xlsx`, "raw data" sheet) or - use the bundled 14-experiment sample, filter experiments by pathogen group, +- **Library:** explores the QMRA wiki workbook (`.xlsx`, "raw data" sheet) + loaded on the Data tab, or the bundled 14-experiment sample, and will be + the starting point for pooling. Filter experiments by pathogen group, pathogen, outcome group, dose unit, host, exposure route, and preferred model, and draw every experiment's preferred model on one log-dose plot, faceted and coloured by any of those fields. Curves are refit from the counts; experiments the workbook reports only as parameters are drawn from those (dashed). Dose units differ between studies, so the plot says so unless it is faceted by - dose unit or filtered to one. "Send to Data and Fit" loads one experiment, - with its study id, host, unit, and endpoint, as the current data; experiments + dose unit or filtered to one. "Send to Data and Fit" makes one experiment the current dataset, + with its study id, host, unit, and endpoint; experiments singlehit cannot fit (fewer than three doses, parameters only) say why instead. Free-text host, route, and outcome spellings are grouped by simple rules in `R/fct_library.R`. The full workbook is not in the repository; the diff --git a/app/tests/testthat/test-app_server.R b/app/tests/testthat/test-app_server.R new file mode 100644 index 0000000..fbf4d69 --- /dev/null +++ b/app/tests/testthat/test-app_server.R @@ -0,0 +1,10 @@ +test_that("an experiment sent from the Library becomes the current dataset", { + shiny::testServer(app_server, { + session$setInputs(`upload-sample_workbook` = 1) + key <- "71|Rotavirus|infection|beta-Poisson" + session$setInputs(`library-experiment` = key, `library-send` = 1) + expect_equal(nrow(data()), 8L) + expect_equal(unique(data()$dose_unit), "FFU") + expect_false(is.null(upload$workbook())) + }) +}) diff --git a/app/tests/testthat/test-mod_library.R b/app/tests/testthat/test-mod_library.R index 5737ae2..1c4a4fb 100644 --- a/app/tests/testthat/test-mod_library.R +++ b/app/tests/testthat/test-mod_library.R @@ -1,6 +1,10 @@ +sample_rows <- function() { + rows <- read_qmra_workbook(qmra_sample_path()) + shiny::reactive(rows) +} + test_that("the sample workbook loads and the filters narrow the experiments", { - shiny::testServer(mod_library_server, { - session$setInputs(sample = 1) + shiny::testServer(mod_library_server, args = list(rows = sample_rows()), { expect_equal(nrow(experiments()), 14L) expect_equal(nrow(filtered()), 14L) @@ -14,8 +18,7 @@ test_that("the sample workbook loads and the filters narrow the experiments", { }) test_that("sending an experiment returns its data with a changing counter", { - shiny::testServer(mod_library_server, { - session$setInputs(sample = 1) + shiny::testServer(mod_library_server, args = list(rows = sample_rows()), { key <- experiments()$experiment_key[experiments()$pathogen == "Rotavirus"] session$setInputs(experiment = key, send = 1) first <- session$returned() @@ -29,20 +32,19 @@ test_that("sending an experiment returns its data with a changing counter", { }) test_that("an experiment that cannot be fitted is not sent", { - shiny::testServer(mod_library_server, { - session$setInputs(sample = 1) + shiny::testServer(mod_library_server, args = list(rows = sample_rows()), { key <- experiments()$experiment_key[experiments()$experiment_id == "Ojielo2003"] session$setInputs(experiment = key, send = 1) expect_null(session$returned()) }) }) -test_that("a file that is not a QMRA workbook shows an error", { - skip_if_not_installed("writexl") - path <- withr::local_tempfile(fileext = ".xlsx") - writexl::write_xlsx(list(other = data.frame(x = 1)), path) - shiny::testServer(mod_library_server, { - session$setInputs(file = list(datapath = path)) - expect_s3_class(rows(), "error") +test_that("the Library is empty until a workbook is loaded", { + rows <- shiny::reactiveVal(NULL) + shiny::testServer(mod_library_server, args = list(rows = rows), { + expect_null(session$returned()) + rows(read_qmra_workbook(qmra_sample_path())) + session$flushReact() + expect_equal(nrow(experiments()), 14L) }) }) diff --git a/app/tests/testthat/test-mod_upload.R b/app/tests/testthat/test-mod_upload.R index efa0158..35eab64 100644 --- a/app/tests/testthat/test-mod_upload.R +++ b/app/tests/testthat/test-mod_upload.R @@ -1,7 +1,7 @@ test_that("the example data loads through singlehit's public reader", { shiny::testServer(mod_upload_server, { session$setInputs(example = 1) - data <- session$returned() + data <- session$returned$data() expect_s3_class(data, "tbl_df") expect_true(all(c("dose", "positive", "negative", "total", "response") %in% names(data))) }) @@ -12,7 +12,7 @@ test_that("invalid uploads return NULL instead of failing the session", { writeLines(c("dose,positive", "1,2"), bad) shiny::testServer(mod_upload_server, { session$setInputs(file = list(datapath = bad)) - expect_null(session$returned()) + expect_null(session$returned$data()) }) }) @@ -22,10 +22,10 @@ test_that("column mapping reads files with unrecognized column names", { shiny::testServer(mod_upload_server, { session$setInputs(file = list(datapath = path)) expect_equal(columns(), c("conc", "infected", "healthy")) - expect_null(session$returned()) + expect_null(session$returned$data()) session$setInputs(dose = "conc", positive = "infected", negative = "healthy") - data <- session$returned() + data <- session$returned$data() expect_equal(data$dose, c(1, 10, 100)) expect_equal(data$positive, c(1, 5, 9)) }) @@ -37,20 +37,20 @@ test_that("loading a new source discards the previous column mapping", { shiny::testServer(mod_upload_server, { session$setInputs(file = list(datapath = path)) session$setInputs(dose = "conc", positive = "infected", negative = "healthy") - expect_false(is.null(session$returned())) + expect_false(is.null(session$returned$data())) # The stale selections are still in `input`, as they would be for a moment # in a real session; the example must still parse with auto-detection. session$setInputs(example = 1) - expect_true(all(c("dose", "positive", "negative") %in% names(session$returned()))) - expect_equal(nrow(session$returned()), 8L) + expect_true(all(c("dose", "positive", "negative") %in% names(session$returned$data()))) + expect_equal(nrow(session$returned$data()), 8L) }) }) test_that("the example button loads the Ward rotavirus data", { shiny::testServer(mod_upload_server, { session$setInputs(example = 1) - data <- session$returned() + data <- session$returned$data() expect_equal(nrow(data), 8L) expect_equal(range(data$dose), c(0.009, 90000)) expect_equal(sum(data$total), 59) @@ -63,14 +63,61 @@ test_that("data sent from the Library tab load without a file", { shiny::testServer(mod_upload_server, args = list(incoming = sent), { sent(list(data = data, label = "Test experiment", sent = 1L)) session$flushReact() - result <- session$returned() + result <- session$returned$data() expect_equal(result$dose, c(1, 10, 100)) expect_equal(unique(result$dose_unit), "CFU") - expect_true(is.list(loaded())) + expect_match(loaded()$label, "Test experiment from the Library tab") # Loading the example afterwards replaces the library data. session$setInputs(example = 1) - expect_false(is.list(loaded())) - expect_equal(nrow(session$returned()), 8L) + expect_equal(loaded()$label, "Ward rotavirus example data") + expect_equal(nrow(session$returned$data()), 8L) + }) +}) + +test_that("the sample workbook loads as a workbook and leaves the dataset alone", { + shiny::testServer(mod_upload_server, { + session$setInputs(example = 1) + session$setInputs(sample_workbook = 1) + expect_equal(length(unique(session$returned$workbook()$experiment_key)), 14L) + expect_equal(session$returned$workbook_loaded(), 1L) + expect_equal(nrow(session$returned$data()), 8L) + + # Loading the same workbook again still counts as a new load. + session$setInputs(sample_workbook = 2) + expect_equal(session$returned$workbook_loaded(), 2L) + }) +}) + +test_that("an uploaded .xlsx is read as a workbook, whatever the case of its extension", { + shiny::testServer(mod_upload_server, { + session$setInputs(file = list(datapath = qmra_sample_path(), name = "QMRA.XLSX")) + expect_equal(length(unique(session$returned$workbook()$experiment_key)), 14L) + expect_null(loaded()) + }) +}) + +test_that("a delimited file after a workbook keeps the workbook", { + path <- withr::local_tempfile(fileext = ".csv") + writeLines(c("dose,positive,negative", "1,1,9", "10,5,5", "100,9,1"), path) + shiny::testServer(mod_upload_server, { + session$setInputs(sample_workbook = 1) + session$setInputs(file = list(datapath = path, name = "trial.csv")) + expect_equal(nrow(session$returned$data()), 3L) + expect_false(is.null(session$returned$workbook())) + expect_equal(loaded()$label, "trial.csv") + }) +}) + +test_that("a bad workbook reports an error and keeps the previous workbook", { + skip_if_not_installed("writexl") + bad <- withr::local_tempfile(fileext = ".xlsx") + writexl::write_xlsx(list(other = data.frame(x = 1)), bad) + shiny::testServer(mod_upload_server, { + session$setInputs(sample_workbook = 1) + session$setInputs(file = list(datapath = bad, name = "bad.xlsx")) + expect_match(workbook_error(), "raw data") + expect_equal(length(unique(session$returned$workbook()$experiment_key)), 14L) + expect_equal(session$returned$workbook_loaded(), 1L) }) }) From ba62ac22f9855b8045814f59a99db10b3463003d Mon Sep 17 00:00:00 2001 From: seanthimons Date: Wed, 28 Oct 2026 14:01:46 -0400 Subject: [PATCH 3/6] feat(pooling)!: refuse to pool datasets in different dose units poolability_test(), fit_pooled_dose_response(), poolability_combinations(), and group_datasets() now check dose_unit metadata before fitting. When any dataset records a dose_unit, every dataset must record one and all must name the same unit (ignoring case and surrounding spaces); otherwise the error lists each dataset's unit. Datasets without unit metadata are not checked. Hosts and endpoints are still not enforced. Refs #21. BREAKING CHANGE: calls that pooled or tested datasets recording different dose units, or a unit on only some datasets, now error instead of running. --- R/pooling.R | 54 +++++++++++++++++++++++++++----- README.md | 10 ++++-- man/fit_pooled_dose_response.Rd | 5 +-- man/group_datasets.Rd | 3 +- man/poolability_combinations.Rd | 3 +- man/poolability_test.Rd | 10 ++++-- tests/testthat/test-pooled-fit.R | 11 +++++++ tests/testthat/test-pooling.R | 32 +++++++++++++++++++ vignettes/getting-started.Rmd | 5 +-- 9 files changed, 115 insertions(+), 18 deletions(-) diff --git a/R/pooling.R b/R/pooling.R index 4f70472..1d43dae 100644 --- a/R/pooling.R +++ b/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,7 +268,7 @@ 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) { @@ -267,6 +276,37 @@ standardize_datasets <- function(datasets) { } ) }) + 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) { diff --git a/README.md b/README.md index c20961d..52cf8ca 100644 --- a/README.md +++ b/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) diff --git a/man/fit_pooled_dose_response.Rd b/man/fit_pooled_dose_response.Rd index 511fbf2..abcfba7 100644 --- a/man/fit_pooled_dose_response.Rd +++ b/man/fit_pooled_dose_response.Rd @@ -30,8 +30,9 @@ two binomial observations. } \details{ The datasets are not tested for poolability here; run -\code{\link[=poolability_test]{poolability_test()}} first. Check metadata compatibility as well, since it -is not enforced. Single-dataset input is allowed and fits that dataset alone. +\code{\link[=poolability_test]{poolability_test()}} first. Datasets that carry \code{dose_unit} metadata must +share one unit (see \code{\link[=poolability_test]{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") diff --git a/man/group_datasets.Rd b/man/group_datasets.Rd index 3d518ee..7b559e5 100644 --- a/man/group_datasets.Rd +++ b/man/group_datasets.Rd @@ -48,7 +48,8 @@ models. Group numbers are labels within a model, not rankings. This function returns assignments, not pooled datasets or fitted models. Use \code{\link[=poolability_combinations]{poolability_combinations()}} to inspect overlapping candidate combinations. -Check metadata compatibility before grouping; it is not enforced here. +Datasets must share one \code{dose_unit} when they carry one, as in +\code{\link[=poolability_test]{poolability_test()}}; hosts and endpoints are not checked. } \examples{ path <- system.file("extdata", "pooling-example.csv", package = "singlehit") diff --git a/man/poolability_combinations.Rd b/man/poolability_combinations.Rd index 4c87b88..dcdd5a6 100644 --- a/man/poolability_combinations.Rd +++ b/man/poolability_combinations.Rd @@ -40,7 +40,8 @@ equivalence, and results with \code{converged = FALSE} may be unreliable. \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 \code{dose_unit} when they carry one, as in +\code{\link[=poolability_test]{poolability_test()}}; hosts and endpoints are not checked. Use \code{\link[=group_datasets]{group_datasets()}} to choose non-overlapping groups from these trials. } \examples{ diff --git a/man/poolability_test.Rd b/man/poolability_test.Rd index 2b46fd0..0c5a3af 100644 --- a/man/poolability_test.Rd +++ b/man/poolability_test.Rd @@ -38,8 +38,14 @@ Datasets are combined by stacking their rows (each dataset's dose groups are kept as separate binomial observations); they are not aggregated, so repeated doses across datasets are preserved. -Check that hosts, dose units, and endpoints match before testing. -The current pooling functions do not enforce metadata compatibility. +Doses are only comparable within one unit, so \code{dose_unit} is +enforced: when any dataset carries a \code{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, +\code{\link[=as_dose_response]{as_dose_response()}} requires a single exact value. Datasets without a +\code{dose_unit} column are not checked. Hosts and endpoints are not enforced; +check that they match before testing. \code{poolable = TRUE} means the test did not detect a significantly worse pooled fit at \code{alpha}; it does not prove equivalence or adequate absolute fit. Inspect \code{converged} before interpreting the result. diff --git a/tests/testthat/test-pooled-fit.R b/tests/testthat/test-pooled-fit.R index 0607aac..df7b227 100644 --- a/tests/testthat/test-pooled-fit.R +++ b/tests/testthat/test-pooled-fit.R @@ -84,3 +84,14 @@ test_that("fit_pooled_dose_response validates its inputs", { bad <- data.frame(dose = c(1, 10, 100), pos = c(0, 0, 1), neg = c(5, 5, 4)) expect_error(fit_pooled_dose_response(list(good = ward_fixture(), bad = bad)), "bad") }) + +test_that("fit_pooled_dose_response refuses mixed dose units", { + a <- ward_fixture() + b <- ward_fixture() + a$dose_unit <- "FFU" + b$dose_unit <- "PFU" + expect_error(fit_pooled_dose_response(list(a = a, b = b), models = "exponential"), "different dose units") + b$dose_unit <- "FFU" + pooled <- fit_pooled_dose_response(list(a = a, b = b), models = "exponential") + expect_equal(unique(pooled$exponential$data$dose_unit), "FFU") +}) diff --git a/tests/testthat/test-pooling.R b/tests/testthat/test-pooling.R index b4e5406..343384e 100644 --- a/tests/testthat/test-pooling.R +++ b/tests/testthat/test-pooling.R @@ -124,3 +124,35 @@ test_that("bundled synthetic collection demonstrates passing and failing combina expect_false(g$group[g$dataset == "A"] == g$group[g$dataset == "C"]) } }) + +with_unit <- function(data, unit) { + data$dose_unit <- unit + data +} + +test_that("pooling functions refuse datasets in different dose units", { + a <- with_unit(ward_fixture(), "FFU") + b <- with_unit(ward_fixture(), "PFU") + expect_error(poolability_test(list(a = a, b = b)), "different dose units.*'a' \\(FFU\\).*'b' \\(PFU\\)") + expect_error(poolability_combinations(list(a = a, b = b)), "different dose units") + expect_error(group_datasets(list(a = a, b = b)), "different dose units") +}) + +test_that("dose units match regardless of case and surrounding spaces", { + a <- with_unit(ward_fixture(), "FFU") + b <- with_unit(ward_fixture(), " ffu ") + result <- poolability_test(list(a = a, b = b)) + expect_true(all(result$poolable)) +}) + +test_that("a dose unit on only some datasets is an error", { + a <- with_unit(ward_fixture(), "FFU") + expect_error( + poolability_test(list(a = a, b = ward_fixture())), + "Only some datasets record `dose_unit` \\(missing from 'b'\\)" + ) +}) + +test_that("datasets without unit metadata are not checked", { + expect_no_error(poolability_test(list(a = ward_fixture(), b = ward_fixture()))) +}) diff --git a/vignettes/getting-started.Rmd b/vignettes/getting-started.Rmd index 728e4e0..ee76d53 100644 --- a/vignettes/getting-started.Rmd +++ b/vignettes/getting-started.Rmd @@ -269,8 +269,9 @@ names(collection$datasets) The output names A, B, and C. Keep the pathogen name once on the collection. Split by `study_id` before coercion because it identifies an individual experiment. `read_dose_response()` handles one experiment at a time and rejects -mixed identity values. Check that hosts, units, and endpoints match before -pooling; the current pooling functions do not enforce that comparison. +mixed identity values. The pooling functions refuse datasets that record +different dose units; they do not compare hosts or endpoints, so check those +match before pooling. ### Inspect every combination From 01d3f7794e2c5fca8f6b3ea7ea9ab9a369864f5d Mon Sep 17 00:00:00 2001 From: seanthimons Date: Wed, 28 Oct 2026 14:41:14 -0400 Subject: [PATCH 4/6] feat(app): add Pooling tab for poolability testing Choose a pooling set on the Library tab and test it on the new Pooling tab with singlehit's poolability_test() and fit_pooled_dose_response(). A written rationale is required before each run, and every run goes into a session log (rationale, acknowledged warnings, p-values, workbook) that counts unadjusted tests and downloads as CSV. Results show the LRT with a low-power caveat, separate versus pooled fits, the spread of parameter estimates, and the wiki's own pooled entry when one matches (21 + 23 reproduces the wiki's "21,23"). Data-driven grouping is not offered. The checks live in one pool_rules table and are a deliberately loose, tunable first cut: only data problems and dose units block; outcome, host, route, and pathogen differences warn. Refs #35, #20. --- app/R/app_server.R | 11 +- app/R/app_ui.R | 1 + app/R/fct_pooling.R | 416 ++++++++++++++++++++++++++ app/R/mod_library.R | 53 +++- app/R/mod_pooling.R | 261 ++++++++++++++++ app/R/mod_upload.R | 4 +- app/README.md | 39 ++- app/tests/testthat/test-app_server.R | 9 + app/tests/testthat/test-fct_pooling.R | 128 ++++++++ app/tests/testthat/test-mod_library.R | 27 +- app/tests/testthat/test-mod_pooling.R | 128 ++++++++ app/tests/testthat/test-mod_upload.R | 4 + 12 files changed, 1062 insertions(+), 19 deletions(-) create mode 100644 app/R/fct_pooling.R create mode 100644 app/R/mod_pooling.R create mode 100644 app/tests/testthat/test-fct_pooling.R create mode 100644 app/tests/testthat/test-mod_pooling.R diff --git a/app/R/app_server.R b/app/R/app_server.R index 49d37cf..f2d2d73 100644 --- a/app/R/app_server.R +++ b/app/R/app_server.R @@ -10,14 +10,19 @@ app_server <- function(input, output, session) { # cycle: Data is built first and reads whatever Library last sent. sent <- reactiveVal(NULL) upload <- mod_upload_server("upload", incoming = sent) - library_sent <- mod_library_server("library", upload$workbook) + library_tab <- mod_library_server("library", upload$workbook) data <- upload$data - observeEvent(library_sent(), { - sent(library_sent()) + 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) + 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)) analysis <- mod_fit_server("fit", data) mod_compare_server("compare", analysis) # Plots and Explore get the analysis with the bootstraps attached once a run diff --git a/app/R/app_ui.R b/app/R/app_ui.R index 8b7201b..f9a93c2 100644 --- a/app/R/app_ui.R +++ b/app/R/app_ui.R @@ -19,6 +19,7 @@ app_ui <- function(request) { 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")) diff --git a/app/R/fct_pooling.R b/app/R/fct_pooling.R new file mode 100644 index 0000000..788a1ab --- /dev/null +++ b/app/R/fct_pooling.R @@ -0,0 +1,416 @@ +#' Experiments chosen for a poolability test, with stable labels +#' +#' Experiment ids repeat across pathogens ("1" is both Yersinia pestis and +#' Variola major), so each experiment is labelled by pathogen and id, adding +#' the outcome when that is still ambiguous. +#' +#' @param rows Output of [read_qmra_workbook()]. +#' @param experiments Output of [library_experiments()]. +#' @param keys `experiment_key` values, in the order chosen. +#' @return The matching rows of `experiments` with `label` and `signature` +#' (the experiment's dose groups, to spot the same counts entered twice). +#' @noRd +pool_selection <- function(rows, experiments, keys) { + selection <- experiments[match(unique(keys), experiments$experiment_key), , drop = FALSE] + selection <- selection[!is.na(selection$experiment_key), , drop = FALSE] + if (nrow(selection) == 0L) { + selection$label <- character() + selection$signature <- character() + return(selection) + } + label <- paste0(selection$pathogen, " #", selection$experiment_id) + ambiguous <- label %in% label[duplicated(label)] + label[ambiguous] <- paste0(label[ambiguous], " (", selection$response_type[ambiguous], ")") + selection$label <- make.unique(label, sep = " ") + selection$signature <- vapply(selection$experiment_key, function(key) { + counts <- rows[rows$experiment_key == key & !is.na(rows$dose), c("dose", "positive", "negative"), drop = FALSE] + if (nrow(counts) == 0L) { + return(NA_character_) + } + counts <- counts[order(counts$dose, counts$positive, counts$negative), , drop = FALSE] + paste(counts$dose, counts$positive, counts$negative, sep = ":", collapse = ";") + }, character(1L), USE.NAMES = FALSE) + rownames(selection) <- NULL + selection +} + +#' @noRd +pool_normalise <- function(x) tolower(trimws(as.character(x))) + +#' @noRd +pool_differs <- function(values, what) { + values <- unique(values) + if (length(unique(pool_normalise(values))) > 1L) { + sprintf("The experiments differ in %s: %s.", what, paste(values, collapse = "; ")) + } else { + character() + } +} + +# PINCH POINT: these rules, and above all their severities, are a first, +# deliberately loose cut and need revisiting (see #20, #21, #35). "block" stops +# the test from running; "warn" is shown and logged, but the test still runs. +# Only data problems and dose units block for now. Outcome, host, route, and +# pathogen differences warn, because the QMRA wiki itself pools across hosts +# and routes (for example "21,23" is guinea pig plus diabetic rat). +# Each check takes a pool_selection() and returns the problems it finds. +pool_rules <- list( + list( + id = "size", + severity = "block", + label = "At least two experiments", + check = function(selection) { + if (nrow(selection) < 2L) "Choose at least two experiments to test." else character() + } + ), + list( + id = "fittable", + severity = "block", + label = "Every experiment can be fitted", + check = function(selection) { + bad <- selection[!selection$fittable, , drop = FALSE] + if (nrow(bad) == 0L) character() else paste0(bad$label, ": ", bad$problem) + } + ), + list( + id = "wiki_pooled", + severity = "block", + label = "No entry the wiki already pooled", + check = function(selection) { + bad <- selection$label[selection$wiki_pooled] + if (length(bad) == 0L) { + character() + } else { + sprintf("%s already pools other experiments, so its subjects would be counted twice.", bad) + } + } + ), + list( + id = "duplicate", + severity = "block", + label = "No experiment chosen twice", + check = function(selection) { + id <- paste(selection$pathogen, selection$experiment_id) + repeated <- unique(id[duplicated(id)]) + signature <- selection$signature + copies <- unique(signature[!is.na(signature) & duplicated(signature)]) + c( + if (length(repeated) > 0L) { + sprintf("%s appears more than once (under different outcomes or models).", repeated) + }, + vapply(copies, function(value) { + sprintf("%s have identical dose-response counts.", paste(selection$label[signature %in% value], collapse = " and ")) + }, character(1L), USE.NAMES = FALSE) + ) + } + ), + list( + id = "dose_unit", + severity = "block", + label = "One reported dose unit", + check = function(selection) { + c( + if (any(selection$dose_unit == "Not reported")) { + sprintf("No dose unit is reported for %s.", paste(selection$label[selection$dose_unit == "Not reported"], collapse = ", ")) + }, + pool_differs(selection$dose_unit, "dose unit") + ) + } + ), + list(id = "outcome_group", severity = "warn", label = "One outcome group", check = function(selection) pool_differs(selection$outcome_group, "outcome group")), + list(id = "response_type", severity = "warn", label = "One outcome description", check = function(selection) pool_differs(selection$response_type, "reported outcome")), + list(id = "host", severity = "warn", label = "One host", check = function(selection) pool_differs(selection$host, "host")), + list(id = "route", severity = "warn", label = "One exposure route", check = function(selection) pool_differs(selection$route, "exposure route")), + list(id = "pathogen", severity = "warn", label = "One pathogen", check = function(selection) pool_differs(selection$pathogen, "pathogen")) +) + +#' Run every pooling rule against a selection +#' +#' @param selection Output of [pool_selection()]. +#' @return One row per rule: `id`, `severity`, `label`, `passed`, and +#' `message` (the problems found, joined). +#' @noRd +pool_check <- function(selection) { + out <- lapply(pool_rules, function(rule) { + problems <- rule$check(selection) + data.frame( + id = rule$id, + severity = rule$severity, + label = rule$label, + passed = length(problems) == 0L, + message = paste(problems, collapse = " "), + stringsAsFactors = FALSE + ) + }) + do.call(rbind, out) +} + +#' Whether any blocking rule failed +#' @noRd +pool_blocked <- function(checks) { + any(!checks$passed & checks$severity == "block") +} + +#' Datasets for the pooling functions, named by experiment label +#' +#' Each dataset's `study_id` is its label, so plots keep experiments apart, and +#' every dataset gets the first experiment's spelling of the shared dose unit +#' so the pooled rows carry one unit. +#' @noRd +pool_datasets <- function(rows, selection) { + datasets <- lapply(seq_len(nrow(selection)), function(i) { + data <- library_experiment_data(rows, selection$experiment_key[[i]]) + data$study_id <- selection$label[[i]] + data$dose_unit <- selection$dose_unit[[1L]] + data + }) + stats::setNames(datasets, selection$label) +} + +#' Run the poolability test, the pooled fit, and the separate fits +#' +#' Warnings from fits that do not converge are collected rather than stopping +#' the run. An error, for example a model with no workable starting value, +#' is returned as the condition. +#' +#' @param datasets Output of [pool_datasets()]. +#' @param models Models to test. +#' @param alpha Significance level. +#' @return A list of `test`, `pooled`, `separate` (fits by model, then by +#' experiment, `NULL` where a fit failed), `datasets`, `models`, `alpha`, and +#' `warnings`; or an error condition. +#' @noRd +pool_run <- function(datasets, models, alpha = 0.05) { + warnings <- character() + collect <- function(expr) { + withCallingHandlers(expr, warning = function(w) { + warnings <<- c(warnings, conditionMessage(w)) + invokeRestart("muffleWarning") + }) + } + tryCatch( + { + test <- collect(singlehit::poolability_test(datasets, models = models, alpha = alpha)) + pooled <- collect(singlehit::fit_pooled_dose_response(datasets, models = models)) + separate <- lapply(stats::setNames(models, models), function(model) { + lapply(datasets, function(data) { + tryCatch(collect(singlehit::fit_dose_response(data, model)), error = function(e) { + warnings <<- c(warnings, conditionMessage(e)) + NULL + }) + }) + }) + list( + test = test, + pooled = pooled, + separate = separate, + datasets = datasets, + models = models, + alpha = alpha, + warnings = unique(warnings) + ) + }, + error = function(e) e + ) +} + +#' The test results as a display table +#' +#' Deliberately avoids saying the experiments "are poolable": with few doses +#' per experiment the test has little power, so not rejecting pooling is not +#' evidence that the experiments share one curve (#20). +#' @noRd +pool_test_table <- function(run) { + test <- run$test + data.frame( + Model = model_label_ui(test$model), + `Pooled deviance` = format_number(test$deviance_pooled), + `Separate deviance` = format_number(test$deviance_unpooled), + `LRT statistic` = format_number(test$lrt_statistic), + df = test$df, + `Critical value` = format_number(test$chi_square_critical), + `p-value` = format_p(test$p_value), + Result = ifelse( + test$poolable, + sprintf("Pooling not rejected at alpha = %s", run$alpha), + sprintf("Pooling rejected at alpha = %s", run$alpha) + ), + Converged = ifelse(test$converged, "yes", "no"), + check.names = FALSE + ) +} + +#' Observed points, separate curves, and pooled curves for plotting +#' @noRd +pool_curves <- function(run, points = 100L) { + observed <- do.call(rbind, lapply(names(run$datasets), function(label) { + data <- singlehit::as_dose_response(run$datasets[[label]]) + data.frame(experiment = label, dose = data$dose, response = data$response) + })) + curve <- function(fit, label, kind, model) { + if (is.null(fit)) { + return(NULL) + } + out <- singlehit::prediction_curve(fit, points = points, dose_range = range(fit$data$dose)) + data.frame(experiment = label, kind = kind, model = model_label_ui(model), dose = out$dose, estimate = out$estimate) + } + curves <- do.call(rbind, lapply(run$models, function(model) { + separate <- lapply(names(run$separate[[model]]), function(label) { + curve(run$separate[[model]][[label]], label, "Separate fit", model) + }) + rbind(do.call(rbind, separate), curve(run$pooled[[model]], "Pooled", "Pooled fit", model)) + })) + list(observed = observed, curves = curves) +} + +#' Plot each experiment's own fit against the pooled fit, one panel per model +#' @noRd +pool_plot <- function(run) { + parts <- pool_curves(run) + separate <- parts$curves[parts$curves$kind == "Separate fit", , drop = FALSE] + pooled <- parts$curves[parts$curves$kind == "Pooled fit", , drop = FALSE] + observed <- merge(parts$observed, data.frame(model = unique(parts$curves$model)), by = NULL) + # Panels in the order the models were chosen, not alphabetically. + levels <- model_label_ui(run$models) + separate$model <- factor(separate$model, levels = levels) + pooled$model <- factor(pooled$model, levels = levels) + observed$model <- factor(observed$model, levels = levels) + unit <- run$datasets[[1L]]$dose_unit[[1L]] + ggplot2::ggplot(mapping = ggplot2::aes(x = .data$dose)) + + ggplot2::geom_line( + data = separate, + ggplot2::aes(y = .data$estimate, colour = .data$experiment), + linewidth = 0.5, + alpha = 0.8 + ) + + ggplot2::geom_line(data = pooled, ggplot2::aes(y = .data$estimate), colour = "black", linewidth = 1.1) + + ggplot2::geom_point(data = observed, ggplot2::aes(y = .data$response, colour = .data$experiment), size = 2) + + ggplot2::facet_wrap(ggplot2::vars(.data$model)) + + ggplot2::scale_x_log10(breaks = scales::breaks_log(n = 5), labels = scales::label_log()) + + ggplot2::scale_y_continuous(limits = c(0, 1)) + + ggplot2::labs( + x = sprintf("Dose (%s, log10 scale)", unit), + y = "Proportion responding", + colour = "Experiment", + caption = "Coloured: each experiment's own fit and observed proportions. Black: the pooled fit to the stacked dose groups." + ) + + ggplot2::theme_minimal(base_size = 12) + + ggplot2::theme(legend.position = "bottom", legend.direction = "vertical", plot.caption = ggplot2::element_text(hjust = 0)) +} + +#' Parameter estimates for each experiment and for the pooled fit +#' +#' Wald intervals from [generics::tidy()], which are `NA` when the fit's +#' Hessian is singular (common for beta-Poisson with few doses). +#' @noRd +pool_parameters <- function(run) { + one <- function(fit, label, model) { + if (is.null(fit)) { + return(NULL) + } + tidy <- as.data.frame(generics::tidy(fit, conf.int = TRUE)) + data.frame( + model = model_label_ui(model), + experiment = label, + term = tidy$term, + estimate = tidy$estimate, + conf.low = tidy$conf.low, + conf.high = tidy$conf.high + ) + } + do.call(rbind, lapply(run$models, function(model) { + separate <- lapply(names(run$separate[[model]]), function(label) one(run$separate[[model]][[label]], label, model)) + rbind(do.call(rbind, separate), one(run$pooled[[model]], "Pooled", model)) + })) +} + +#' Forest plot of parameter estimates, to show spread between experiments +#' +#' A direct view of heterogeneity that does not depend on the test's power: +#' if the experiments' intervals barely overlap, a single pooled curve hides +#' real differences whatever the test says. +#' @noRd +pool_forest <- function(run) { + parameters <- pool_parameters(run) + panels <- unique(paste(parameters$model, parameters$term, sep = ": ")) + parameters$panel <- factor(paste(parameters$model, parameters$term, sep = ": "), levels = panels) + parameters$experiment <- factor(parameters$experiment, levels = rev(c(names(run$datasets), "Pooled"))) + parameters$pooled <- parameters$experiment == "Pooled" + intervals <- parameters[is.finite(parameters$conf.low) & is.finite(parameters$conf.high), , drop = FALSE] + reference <- parameters[parameters$pooled, , drop = FALSE] + ggplot2::ggplot(parameters, ggplot2::aes(x = .data$estimate, y = .data$experiment)) + + ggplot2::geom_vline(data = reference, ggplot2::aes(xintercept = .data$estimate), linetype = "dashed", colour = "grey40") + + ggplot2::geom_linerange( + data = intervals, + ggplot2::aes(xmin = .data$conf.low, xmax = .data$conf.high, colour = .data$pooled), + linewidth = 0.8 + ) + + ggplot2::geom_point(ggplot2::aes(colour = .data$pooled, shape = .data$pooled), size = 2.6) + + ggplot2::scale_colour_manual(values = c(`FALSE` = "#0066A1", `TRUE` = "black"), guide = "none") + + ggplot2::scale_shape_manual(values = c(`FALSE` = 16, `TRUE` = 18), guide = "none") + + # Plain numbers: power-of-ten labels read badly over less than a decade. + ggplot2::scale_x_log10(breaks = scales::breaks_log(n = 4), labels = function(x) format_number(x, digits = 3)) + + ggplot2::facet_wrap(ggplot2::vars(.data$panel), scales = "free_x") + + ggplot2::labs( + x = "Estimate (log10 scale) with 95% Wald interval", + y = NULL, + caption = "Dashed line: pooled estimate. A point without a bar has no interval (the fit's information matrix is singular)." + ) + + ggplot2::theme_minimal(base_size = 12) + + ggplot2::theme(plot.caption = ggplot2::element_text(hjust = 0)) +} + +#' Experiment ids listed in a wiki-pooled entry such as "281,282 (pooled)" +#' @noRd +pool_wiki_ids <- function(id) { + id <- gsub("\\(pooled\\)", "", id, ignore.case = TRUE) + parts <- trimws(unlist(strsplit(id, ",|&|\\band\\b"))) + parts[nzchar(parts)] +} + +#' The wiki's own pooled entry for exactly these experiments, if any +#' +#' Lets a pooled fit be checked against what the wiki reports, for example +#' 21 + 23 for Burkholderia mallei against the wiki's "21,23". +#' @return A display table, or `NULL` when no entry matches. +#' @noRd +pool_wiki_reference <- function(rows, experiments, selection) { + if (length(unique(selection$pathogen)) != 1L) { + return(NULL) + } + candidates <- experiments[experiments$wiki_pooled & experiments$pathogen == selection$pathogen[[1L]], , drop = FALSE] + match <- vapply(candidates$experiment_id, function(id) setequal(pool_wiki_ids(id), selection$experiment_id), logical(1L)) + candidates <- candidates[match, , drop = FALSE] + if (nrow(candidates) == 0L) { + return(NULL) + } + data.frame( + `Wiki entry` = candidates$experiment_id, + Outcome = candidates$response_type, + `Preferred model` = candidates$model_label, + `Reported parameters` = rows$parameters[match(candidates$experiment_key, rows$experiment_key)], + check.names = FALSE + ) +} + +#' Log rows for one run: one row per model +#' @noRd +pool_log_entry <- function(run, selection, checks, rationale, workbook, time = Sys.time()) { + warned <- checks$label[!checks$passed & checks$severity == "warn"] + data.frame( + time = format(time, "%Y-%m-%d %H:%M:%S"), + workbook = if (is.null(workbook)) "" else workbook, + experiments = paste(selection$label, collapse = "; "), + dose_unit = selection$dose_unit[[1L]], + rationale = rationale, + warnings_acknowledged = paste(warned, collapse = "; "), + model = run$test$model, + alpha = run$alpha, + lrt_statistic = run$test$lrt_statistic, + df = run$test$df, + p_value = run$test$p_value, + pooling_rejected = !run$test$poolable, + converged = run$test$converged, + stringsAsFactors = FALSE + ) +} diff --git a/app/R/mod_library.R b/app/R/mod_library.R index 82b55ec..3ed90cb 100644 --- a/app/R/mod_library.R +++ b/app/R/mod_library.R @@ -1,8 +1,9 @@ #' Library module UI #' #' Explore the QMRA wiki workbook loaded on the Data tab: filter its -#' experiments, draw their preferred models faceted by metadata, and send one -#' experiment back to the Data and Fit tabs. +#' experiments, draw their preferred models faceted by metadata, send one +#' experiment back to the Data and Fit tabs, or choose several to test for +#' poolability on the Pooling tab. #' #' @param id Module id. #' @noRd @@ -37,6 +38,27 @@ mod_library_ui <- function(id) { uiOutput(ns("status")), plotOutput(ns("plot"), height = "620px", fill = FALSE), uiOutput(ns("send_panel")), + div( + class = "mt-2 mb-3 p-3 border rounded", + h6("Pooling set"), + p( + class = "text-muted small", + "Experiments to test for poolability on the Pooling tab. Only experiments that can be fitted and that the wiki", + "has not already pooled are offered. The list follows the filters but keeps what you have already chosen." + ), + div( + class = "d-flex flex-wrap align-items-end gap-3", + selectizeInput( + ns("pool_set"), + NULL, + choices = character(), + multiple = TRUE, + width = "520px", + options = list(placeholder = "Choose two or more experiments") + ), + actionButton(ns("pool"), "Test poolability", icon = icon("layer-group"), class = "btn-outline-primary mb-3") + ) + ), tableOutput(ns("table")) ) } @@ -46,8 +68,11 @@ mod_library_ui <- function(id) { #' @param id Module id. #' @param rows Reactive returning the workbook rows from [read_qmra_workbook()], #' or `NULL` before a workbook is loaded on the Data tab. -#' @return A reactive that is `NULL` until an experiment is sent, then a list -#' of `data`, `label`, and `sent` (a counter that changes on every send). +#' @return A list of reactives: `sent`, `NULL` until an experiment is sent and +#' then a list of `data`, `label`, and `sent` (a counter that changes on every +#' send); `pool_set`, `NULL` until "Test poolability" is pressed and then a +#' list of `keys` and a `sent` counter; and `experiments`, the +#' [library_experiments()] table, so other tabs need not rebuild it. #' @noRd mod_library_server <- function(id, rows) { moduleServer(id, function(input, output, session) { @@ -160,6 +185,24 @@ mod_library_server <- function(id, rows) { ) }) + observeEvent(filtered(), { + all <- experiments() + selected <- intersect(isolate(input$pool_set), all$experiment_key) + shown <- filtered() + eligible <- shown$experiment_key[shown$fittable & !shown$wiki_pooled] + offered <- all[all$experiment_key %in% union(eligible, selected), , drop = FALSE] + choices <- stats::setNames(offered$experiment_key, library_experiment_label(offered)) + updateSelectizeInput(session, "pool_set", choices = choices, selected = selected) + }) + + pool_sent <- reactiveVal(NULL) + pool_count <- reactiveVal(0L) + observeEvent(input$pool, { + req(length(input$pool_set) > 0L) + pool_count(pool_count() + 1L) + pool_sent(list(keys = input$pool_set, sent = pool_count())) + }) + sent <- reactiveVal(NULL) count <- reactiveVal(0L) observeEvent(input$send, { @@ -172,6 +215,6 @@ mod_library_server <- function(id, rows) { sent = count() )) }) - sent + list(sent = sent, pool_set = pool_sent, experiments = experiments) }) } diff --git a/app/R/mod_pooling.R b/app/R/mod_pooling.R new file mode 100644 index 0000000..fca387c --- /dev/null +++ b/app/R/mod_pooling.R @@ -0,0 +1,261 @@ +#' Pooling module UI +#' +#' Tests whether experiments chosen on the Library tab can share one +#' dose-response curve, with the checks and the session log that guard against +#' choosing groups after seeing the fits (#20). +#' +#' @param id Module id. +#' @noRd +mod_pooling_ui <- function(id) { + ns <- NS(id) + bslib::layout_sidebar( + fillable = FALSE, + sidebar = bslib::sidebar( + width = 340, + checkboxGroupInput( + ns("models"), + "Models", + choices = c("Exponential" = "exponential", "Beta-Poisson" = "beta_poisson", "Exact beta-Poisson (very slow)" = "exact_beta_poisson"), + selected = c("exponential", "beta_poisson") + ), + numericInput(ns("alpha"), "Significance level (alpha)", value = 0.05, min = 0.001, max = 0.5, step = 0.01), + textAreaInput( + ns("rationale"), + "Why do these experiments belong together?", + rows = 4, + placeholder = "Decide this before running: for example, same pathogen, host, and route, so one curve is plausible." + ), + p( + class = "text-muted small", + "Required. It is recorded with the result, so groups are justified before the test, not chosen from it." + ), + uiOutput(ns("run_button")), + hr(), + h6("Session log"), + uiOutput(ns("log_status")), + downloadButton(ns("download_log"), "Download log (CSV)", class = "btn-sm") + ), + uiOutput(ns("empty")), + uiOutput(ns("selection_ui")), + uiOutput(ns("results_ui")) + ) +} + +#' Pooling module server +#' +#' @param id Module id. +#' @param rows Reactive returning the workbook rows, or `NULL`. +#' @param experiments Reactive returning the [library_experiments()] table. +#' @param pool_set Reactive returning `NULL` or a list of `keys` and `sent`. +#' @param workbook_name Reactive returning where the workbook came from. +#' @param workbook_loaded Reactive counter that changes on every workbook load, +#' including reloading the same file, which leaves `rows()` unchanged. +#' @return Nothing useful; the log is kept in the session. +#' @noRd +mod_pooling_server <- function(id, rows, experiments, pool_set, workbook_name = reactive(NULL), workbook_loaded = reactive(0L)) { + moduleServer(id, function(input, output, session) { + keys <- reactiveVal(character()) + result <- reactiveVal(NULL) + # The log lasts the whole session, across workbooks (each row records its + # workbook), so reloading cannot erase the record of tests already run. + log <- reactiveVal(NULL) + # A rationale is written for one set of experiments. After the set or the + # workbook changes, the old text no longer counts until it is edited. + rationale_stale <- reactiveVal(FALSE) + + start_over <- function(new_keys) { + keys(new_keys) + result(NULL) + rationale_stale(TRUE) + updateTextAreaInput(session, "rationale", value = "") + } + # Ahead of the other observers, so nothing reads the new workbook with the + # previous workbook's keys. + observeEvent(workbook_loaded(), start_over(character()), priority = 1) + observeEvent(pool_set(), start_over(pool_set()$keys), priority = 1) + observeEvent(input$rationale, rationale_stale(FALSE), ignoreInit = TRUE) + + selection <- reactive({ + req(rows(), length(keys()) > 0L) + pool_selection(rows(), experiments(), keys()) + }) + checks <- reactive(pool_check(selection())) + rationale <- reactive({ + if (rationale_stale() || is.null(input$rationale)) "" else trimws(input$rationale) + }) + ready <- reactive({ + !pool_blocked(checks()) && nzchar(rationale()) && length(input$models) > 0L && + is.numeric(input$alpha) && isTRUE(input$alpha > 0 && input$alpha < 1) + }) + + output$empty <- renderUI({ + if (is.null(rows())) { + p(class = "text-muted", "Load the QMRA wiki workbook on the Data tab, then choose a pooling set on the Library tab.") + } else if (length(keys()) == 0L) { + p(class = "text-muted", "Choose experiments in the Pooling set on the Library tab and press Test poolability.") + } + }) + + output$selection_ui <- renderUI({ + req(nrow(selection()) > 0L) + tagList( + h5("Experiments"), + tableOutput(session$ns("selection_table")), + h5("Checks"), + p( + class = "text-muted small", + "Blocking checks stop the test; warnings are shown and logged. These rules are a first, loose cut and will be", + "revised (see pool_rules in R/fct_pooling.R)." + ), + tableOutput(session$ns("checks_table")) + ) + }) + + output$selection_table <- renderTable({ + shown <- selection() + data.frame( + Experiment = shown$label, + Outcome = shown$response_type, + `Dose unit` = shown$dose_unit, + Host = shown$host, + Route = shown$route, + Doses = shown$doses, + Subjects = format(shown$subjects, big.mark = ","), + check.names = FALSE + ) + }) + + output$checks_table <- renderTable({ + shown <- checks() + data.frame( + Check = shown$label, + Severity = ifelse(shown$severity == "block", "Blocks", "Warns"), + Result = ifelse(shown$passed, "OK", ifelse(shown$severity == "block", "BLOCKED", "Warning")), + Detail = shown$message, + check.names = FALSE + ) + }) + + # Why Run is disabled, or "" when it is not. Held in a reactiveVal, which + # ignores a set to the same value, so the button is redrawn only when its + # state changes and not on every keystroke in the rationale. + run_blocker <- reactiveVal("Choose a pooling set on the Library tab first.") + observe(run_blocker( + if (is.null(rows()) || length(keys()) == 0L) { + "Choose a pooling set on the Library tab first." + } else if (pool_blocked(checks())) { + "A blocking check failed; see Checks." + } else if (!nzchar(rationale())) { + "Write why these experiments belong together first." + } else if (length(input$models) == 0L) { + "Choose at least one model." + } else if (!isTRUE(ready())) { + "Choose a significance level between 0 and 1." + } else { + "" + } + )) + + output$run_button <- renderUI({ + reason <- run_blocker() + tagList( + { + button <- actionButton(session$ns("run"), "Run poolability test", icon = icon("play"), class = "btn-primary w-100") + # Set the attribute directly: actionButton()'s own `disabled` argument + # is recent. The server also refuses a run that is not ready. + if (nzchar(reason)) tagAppendAttributes(button, disabled = NA) else button + }, + if (nzchar(reason)) p(class = "text-muted small mt-1", reason) + ) + }) + + observeEvent(input$run, { + req(isTRUE(ready())) + chosen <- selection() + run <- withProgress(message = "Fitting separate and pooled models", { + tryCatch( + pool_run(pool_datasets(rows(), chosen), models = input$models, alpha = input$alpha), + error = function(e) e + ) + }) + result(list(run = run, selection = chosen)) + if (!inherits(run, "error")) { + entry <- pool_log_entry(run, chosen, checks(), rationale(), workbook_name()) + previous <- log() + entry <- cbind(run = if (is.null(previous)) 1L else max(previous$run) + 1L, entry) + log(rbind(previous, entry)) + } + }) + + output$log_status <- renderUI({ + entries <- log() + tests <- if (is.null(entries)) 0L else length(unique(entries$run)) + p( + class = "small", + sprintf("%d poolability %s run this session.", tests, if (tests == 1L) "test" else "tests"), + if (tests > 1L) " The p-values are not adjusted for running several tests." + ) + }) + + output$download_log <- downloadHandler( + filename = function() paste0("singlehit-pooling-log-", format(Sys.Date()), ".csv"), + content = function(file) { + entries <- log() + if (is.null(entries)) { + entries <- data.frame(note = "No poolability tests have been run in this session.") + } + utils::write.csv(entries, file, row.names = FALSE) + } + ) + + output$results_ui <- renderUI({ + current <- result() + req(current) + if (inherits(current$run, "error")) { + return(div(class = "alert alert-danger mt-3", "The test could not run: ", conditionMessage(current$run))) + } + run <- current$run + reference <- pool_wiki_reference(rows(), experiments(), current$selection) + tagList( + hr(), + h5("Poolability test"), + div( + class = "alert alert-secondary small", + "The likelihood-ratio test compares one pooled curve with a separate curve per experiment. With a few doses per", + "experiment it has little power, so \"not rejected\" is not evidence that the experiments share one curve;", + "read it together with the spread of the estimates below." + ), + if (length(run$warnings) > 0L) { + div(class = "alert alert-warning small", strong("Fit warnings: "), paste(run$warnings, collapse = " ")) + }, + tableOutput(session$ns("test_table")), + if (!is.null(reference)) { + tagList( + h6("The wiki's own pooled entry for these experiments"), + tableOutput(session$ns("wiki_table")) + ) + }, + h5("Separate and pooled fits"), + plotOutput(session$ns("fit_plot"), height = "480px", fill = FALSE), + h5("Spread of the estimates"), + plotOutput(session$ns("forest_plot"), height = "380px", fill = FALSE) + ) + }) + + run_ok <- reactive({ + current <- result() + req(current, !inherits(current$run, "error")) + current$run + }) + output$test_table <- renderTable(pool_test_table(run_ok())) + output$wiki_table <- renderTable({ + current <- result() + req(current) + pool_wiki_reference(rows(), experiments(), current$selection) + }) + output$fit_plot <- renderPlot(pool_plot(run_ok()), res = 96) + output$forest_plot <- renderPlot(pool_forest(run_ok()), res = 96) + + invisible(NULL) + }) +} diff --git a/app/R/mod_upload.R b/app/R/mod_upload.R index 3413d6d..9bb8df1 100644 --- a/app/R/mod_upload.R +++ b/app/R/mod_upload.R @@ -71,7 +71,8 @@ mod_upload_table_ui <- function(id) { #' of dose groups), `label` (where it came from), and `sent` (a counter, so #' sending the same experiment again reloads it). #' @return A list of reactives: `data`, the standardized data or `NULL`; -#' `workbook`, the workbook rows or `NULL`; and `workbook_loaded`, a counter +#' `workbook`, the workbook rows or `NULL`; `workbook_name`, where it came +#' from; and `workbook_loaded`, a counter #' that changes each time a workbook loads successfully. #' @noRd mod_upload_server <- function(id, incoming = reactive(NULL)) { @@ -213,6 +214,7 @@ mod_upload_server <- function(id, incoming = reactive(NULL)) { list( data = reactive(if (inherits(parsed(), "error")) NULL else parsed()), workbook = reactive(workbook()), + workbook_name = reactive(workbook_name()), workbook_loaded = reactive(workbook_loaded()) ) }) diff --git a/app/README.md b/app/README.md index 6a8a9cb..deb2b05 100644 --- a/app/README.md +++ b/app/README.md @@ -39,9 +39,9 @@ depends on the app. "Use {plotly} to render plot" under Options for hover and zoom. The interface font is overridden in `inst/app/www/custom.css`. -- **Library:** explores the QMRA wiki workbook (`.xlsx`, "raw data" sheet) - loaded on the Data tab, or the bundled 14-experiment sample, and will be - the starting point for pooling. Filter experiments by pathogen group, +- **Library:** explores the QMRA wiki workbook (`.xlsx`, "raw data" sheet), + or the bundled 14-experiment sample, loaded on the Data tab, and is the + starting point for pooling. Filter experiments by pathogen group, pathogen, outcome group, dose unit, host, exposure route, and preferred model, and draw every experiment's preferred model on one log-dose plot, faceted and coloured by any of those fields. Curves are refit from the counts; experiments @@ -54,9 +54,36 @@ depends on the app. rules in `R/fct_library.R`. The full workbook is not in the repository; the sample is rebuilt from it with `dev/make_qmra_sample.R`. The data come from the [QMRA wiki](https://qmrawiki.org) and the studies it cites; each sample - row keeps its reference, DOI, and wiki link. - -Pooling (#35) and export are planned. + row keeps its reference, DOI, and wiki link. Choose several experiments in + the "Pooling set" and press "Test poolability" to open them on the Pooling + tab. + +- **Pooling:** tests whether the chosen experiments can share one curve, with + singlehit's likelihood-ratio `poolability_test()` and a pooled fit from + `fit_pooled_dose_response()`. Before running you must write why the + experiments belong together, and every run is kept in a session log (with + the rationale, the warnings you ran past, and the p-values) that downloads + as CSV and counts how many unadjusted tests have been run. The log lasts the whole + session, across workbooks, and the rationale clears whenever the pooling set + or workbook changes, so each set gets its own justification. Results show the test with the caveat that "not rejected" + is not evidence of equivalence, each experiment's own fit against the pooled + fit, the spread of the parameter estimates, and the wiki's own pooled entry + when the chosen experiments match one (21 + 23 for Burkholderia mallei + reproduces the wiki's "21,23"). The data-driven `group_datasets()` and + `poolability_combinations()` are deliberately not offered: picking groups + from the fits is the forking path described in #20. + + **The checks are a tunable pinch point.** They live in one table, + `pool_rules` in `R/fct_pooling.R`, and are a first, loose cut to be revised. + Only data problems block a test: fewer than two experiments, an experiment + that cannot be fitted, double counting (an entry the wiki already pooled, + the same experiment twice, identical counts), and dose units that differ or + are not reported (singlehit itself also refuses mixed units). Different + outcomes, hosts, routes, or pathogens only warn, because the wiki itself + pools across hosts and routes. + +Hierarchical (family-level) models, bootstrap intervals for pooled fits, and +export are planned. ## License diff --git a/app/tests/testthat/test-app_server.R b/app/tests/testthat/test-app_server.R index fbf4d69..cfac324 100644 --- a/app/tests/testthat/test-app_server.R +++ b/app/tests/testthat/test-app_server.R @@ -8,3 +8,12 @@ test_that("an experiment sent from the Library becomes the current dataset", { expect_false(is.null(upload$workbook())) }) }) + +test_that("a pooling set chosen on the Library reaches the Pooling tab", { + shiny::testServer(app_server, { + session$setInputs(`upload-sample_workbook` = 1) + keys <- c("21|Burkholderia mallei|death|beta-Poisson", "23|Burkholderia mallei|death|beta-Poisson") + session$setInputs(`library-pool_set` = keys, `library-pool` = 1) + expect_equal(library_tab$pool_set()$keys, keys) + }) +}) diff --git a/app/tests/testthat/test-fct_pooling.R b/app/tests/testthat/test-fct_pooling.R new file mode 100644 index 0000000..8614396 --- /dev/null +++ b/app/tests/testthat/test-fct_pooling.R @@ -0,0 +1,128 @@ +sample_pool <- function(ids) { + rows <- read_qmra_workbook(qmra_sample_path()) + experiments <- library_experiments(rows) + keys <- experiments$experiment_key[experiments$experiment_id %in% ids] + list(rows = rows, experiments = experiments, selection = pool_selection(rows, experiments, keys)) +} + +failed <- function(checks) checks$id[!checks$passed] + +test_that("two experiments in one unit pass every blocking check", { + pool <- sample_pool(c("21", "23")) + checks <- pool_check(pool$selection) + expect_false(pool_blocked(checks)) + # The wiki pooled guinea pig with diabetic rat; differing hosts only warn. + expect_equal(failed(checks), "host") + expect_equal(checks$severity[checks$id == "host"], "warn") +}) + +test_that("data problems and dose units block the test", { + pool <- sample_pool(c("21", "21,23", "71", "Ojielo2003")) + checks <- pool_check(pool$selection) + expect_true(pool_blocked(checks)) + expect_true(all(c("fittable", "wiki_pooled", "dose_unit") %in% failed(checks))) + expect_true(all(c("outcome_group", "response_type", "route", "pathogen") %in% failed(checks))) + expect_false("size" %in% failed(checks)) +}) + +test_that("a single experiment, an unreported unit, or repeated counts block", { + expect_true("size" %in% failed(pool_check(sample_pool("21")$selection))) + + pool <- sample_pool(c("21", "23")) + unreported <- pool$selection + unreported$dose_unit[[2L]] <- "Not reported" + expect_match(pool_check(unreported)$message[pool_check(unreported)$id == "dose_unit"], "No dose unit is reported") + + copied <- pool$selection + copied$signature[[2L]] <- copied$signature[[1L]] + expect_true("duplicate" %in% failed(pool_check(copied))) + + repeated <- pool$selection + repeated$experiment_id[[2L]] <- repeated$experiment_id[[1L]] + expect_match(pool_check(repeated)$message[pool_check(repeated)$id == "duplicate"], "appears more than once") +}) + +test_that("labels keep experiments with the same id apart", { + pool <- sample_pool("1") + expect_setequal(pool$selection$label, c("Yersinia pestis #1", "Variola major #1")) +}) + +test_that("every rule has a known severity", { + expect_true(all(vapply(pool_rules, function(rule) rule$severity %in% c("block", "warn"), logical(1L)))) + expect_false(anyDuplicated(vapply(pool_rules, `[[`, character(1L), "id")) > 0L) +}) + +test_that("pooling 21 and 23 reproduces the wiki's own pooled entry", { + pool <- sample_pool(c("21", "23")) + run <- pool_run(pool_datasets(pool$rows, pool$selection), c("exponential", "beta_poisson")) + expect_false(inherits(run, "error")) + beta <- generics::tidy(run$pooled$beta_poisson) + # The wiki reports a = 2.13E-01, N50 = 4.77E+02 for "21,23". + expect_equal(signif(beta$estimate[beta$term == "alpha"], 3), 0.213) + expect_equal(signif(beta$estimate[beta$term == "n50"], 3), 477) + expect_true(run$test$poolable[run$test$model == "beta_poisson"]) + expect_false(run$test$poolable[run$test$model == "exponential"]) + + reference <- pool_wiki_reference(pool$rows, pool$experiments, pool$selection) + expect_equal(reference$`Wiki entry`, "21,23") + expect_match(reference$`Reported parameters`, "4.77E+02", fixed = TRUE) +}) + +test_that("the results table never says the experiments are poolable", { + pool <- sample_pool(c("21", "23")) + run <- pool_run(pool_datasets(pool$rows, pool$selection), c("exponential", "beta_poisson")) + table <- pool_test_table(run) + expect_equal(table$Result, c("Pooling rejected at alpha = 0.05", "Pooling not rejected at alpha = 0.05")) + expect_false(any(grepl("are poolable", unlist(table)))) +}) + +test_that("pooled datasets carry one label per experiment and one unit spelling", { + pool <- sample_pool(c("21", "23")) + datasets <- pool_datasets(pool$rows, pool$selection) + expect_equal(names(datasets), pool$selection$label) + expect_equal(unique(datasets[[2L]]$study_id), pool$selection$label[[2L]]) + expect_equal(unique(unlist(lapply(datasets, `[[`, "dose_unit"))), "CFU") +}) + +test_that("errors are returned, not thrown", { + datasets <- list( + a = data.frame(dose = c(1, 10, 100), positive = c(1, 5, 9), negative = c(9, 5, 1), dose_unit = "CFU"), + b = data.frame(dose = c(1, 10, 100), positive = c(1, 5, 9), negative = c(9, 5, 1), dose_unit = "PFU") + ) + expect_s3_class(pool_run(datasets, "exponential"), "error") +}) + +test_that("the plots build", { + pool <- sample_pool(c("21", "23")) + run <- pool_run(pool_datasets(pool$rows, pool$selection), c("exponential", "beta_poisson")) + expect_no_error(ggplot2::ggplot_build(pool_plot(run))) + expect_no_error(ggplot2::ggplot_build(pool_forest(run))) + parameters <- pool_parameters(run) + expect_setequal(unique(parameters$experiment), c(pool$selection$label, "Pooled")) +}) + +test_that("wiki-pooled ids are split on commas, ampersands, and 'and'", { + expect_equal(pool_wiki_ids("281,282 (pooled)"), c("281", "282")) + expect_equal(pool_wiki_ids("297 & 298"), c("297", "298")) + expect_equal(pool_wiki_ids("300 and 301"), c("300", "301")) + expect_equal(pool_wiki_ids("142, 143, 144"), c("142", "143", "144")) +}) + +test_that("log entries record the run, the rationale, and acknowledged warnings", { + pool <- sample_pool(c("21", "23")) + checks <- pool_check(pool$selection) + run <- pool_run(pool_datasets(pool$rows, pool$selection), "beta_poisson") + entry <- pool_log_entry(run, pool$selection, checks, "Same pathogen, route, and outcome.", "sample") + expect_equal(nrow(entry), 1L) + expect_equal(entry$rationale, "Same pathogen, route, and outcome.") + expect_equal(entry$warnings_acknowledged, "One host") + expect_equal(entry$dose_unit, "CFU") + expect_false(entry$pooling_rejected) +}) + +test_that("a selection that matches nothing is empty, not an error", { + rows <- read_qmra_workbook(qmra_sample_path()) + selection <- pool_selection(rows, library_experiments(rows), "not|a|real|key") + expect_equal(nrow(selection), 0L) + expect_true("size" %in% failed(pool_check(selection))) +}) diff --git a/app/tests/testthat/test-mod_library.R b/app/tests/testthat/test-mod_library.R index 1c4a4fb..68a3b84 100644 --- a/app/tests/testthat/test-mod_library.R +++ b/app/tests/testthat/test-mod_library.R @@ -21,13 +21,13 @@ test_that("sending an experiment returns its data with a changing counter", { shiny::testServer(mod_library_server, args = list(rows = sample_rows()), { key <- experiments()$experiment_key[experiments()$pathogen == "Rotavirus"] session$setInputs(experiment = key, send = 1) - first <- session$returned() + first <- session$returned$sent() expect_equal(first$sent, 1L) expect_match(first$label, "Rotavirus") expect_equal(unique(first$data$dose_unit), "FFU") session$setInputs(send = 2) - expect_equal(session$returned()$sent, 2L) + expect_equal(session$returned$sent()$sent, 2L) }) }) @@ -35,16 +35,35 @@ test_that("an experiment that cannot be fitted is not sent", { shiny::testServer(mod_library_server, args = list(rows = sample_rows()), { key <- experiments()$experiment_key[experiments()$experiment_id == "Ojielo2003"] session$setInputs(experiment = key, send = 1) - expect_null(session$returned()) + expect_null(session$returned$sent()) }) }) test_that("the Library is empty until a workbook is loaded", { rows <- shiny::reactiveVal(NULL) shiny::testServer(mod_library_server, args = list(rows = rows), { - expect_null(session$returned()) + expect_null(session$returned$sent()) rows(read_qmra_workbook(qmra_sample_path())) session$flushReact() expect_equal(nrow(experiments()), 14L) }) }) + +test_that("the pooling set offers only eligible experiments and is handed on", { + shiny::testServer(mod_library_server, args = list(rows = sample_rows()), { + session$flushReact() + expect_null(session$returned$pool_set()) + keys <- experiments()$experiment_key[experiments()$experiment_id %in% c("21", "23")] + session$setInputs(pool_set = keys, pool = 1) + expect_equal(session$returned$pool_set()$keys, keys) + expect_equal(session$returned$pool_set()$sent, 1L) + expect_equal(nrow(session$returned$experiments()), 14L) + }) +}) + +test_that("pressing Test poolability with nothing chosen hands nothing on", { + shiny::testServer(mod_library_server, args = list(rows = sample_rows()), { + session$setInputs(pool = 1) + expect_null(session$returned$pool_set()) + }) +}) diff --git a/app/tests/testthat/test-mod_pooling.R b/app/tests/testthat/test-mod_pooling.R new file mode 100644 index 0000000..2797b70 --- /dev/null +++ b/app/tests/testthat/test-mod_pooling.R @@ -0,0 +1,128 @@ +pooling_args <- function() { + rows <- shiny::reactiveVal(read_qmra_workbook(qmra_sample_path())) + list( + rows = rows, + # Built from the rows, as on the Library tab, so a new workbook brings new + # experiments. + experiments = shiny::reactive(library_experiments(rows())), + pool_set = shiny::reactiveVal(NULL), + workbook_name = shiny::reactive("the sample workbook"), + workbook_loaded = shiny::reactiveVal(1L) + ) +} + +keys_for <- function(experiments, ids) experiments$experiment_key[experiments$experiment_id %in% ids] + +test_that("the test runs only with a rationale and no blocking check", { + args <- pooling_args() + shiny::testServer(mod_pooling_server, args = args, { + args$pool_set(list(keys = keys_for(experiments(), c("21", "23")), sent = 1L)) + session$flushReact() + session$setInputs(models = c("exponential", "beta_poisson"), alpha = 0.05, rationale = "") + expect_false(ready()) + + # Pressing Run without a rationale does nothing. + session$setInputs(run = 1) + expect_null(result()) + expect_null(log()) + + session$setInputs(rationale = "Same pathogen, route, and outcome.") + expect_true(ready()) + session$setInputs(run = 2) + expect_false(inherits(result()$run, "error")) + expect_equal(nrow(log()), 2L) + expect_equal(unique(log()$run), 1L) + expect_equal(unique(log()$workbook), "the sample workbook") + + session$setInputs(run = 3) + expect_equal(unique(log()$run), c(1L, 2L)) + }) +}) + +test_that("a blocking check keeps the test from running", { + args <- pooling_args() + shiny::testServer(mod_pooling_server, args = args, { + args$pool_set(list(keys = keys_for(experiments(), c("21", "71")), sent = 1L)) + session$setInputs(models = "exponential", alpha = 0.05, rationale = "Testing.") + expect_true(pool_blocked(checks())) + expect_false(ready()) + session$setInputs(run = 1) + expect_null(result()) + }) +}) + +test_that("a new workbook clears the selection and results but keeps the log", { + args <- pooling_args() + shiny::testServer(mod_pooling_server, args = args, { + args$pool_set(list(keys = keys_for(experiments(), c("21", "23")), sent = 1L)) + session$flushReact() + session$setInputs(models = "beta_poisson", alpha = 0.05, rationale = "Same pathogen.") + session$setInputs(run = 1) + expect_equal(nrow(log()), 1L) + + # Reloading the same file leaves the rows identical; the load counter still + # starts a new analysis. + args$rows(read_qmra_workbook(qmra_sample_path())) + args$workbook_loaded(2L) + session$flushReact() + expect_equal(keys(), character()) + expect_null(result()) + expect_equal(nrow(log()), 1L) + }) +}) + +test_that("a new pooling set needs its own rationale", { + args <- pooling_args() + shiny::testServer(mod_pooling_server, args = args, { + args$pool_set(list(keys = keys_for(experiments(), c("21", "23")), sent = 1L)) + session$flushReact() + session$setInputs(models = "beta_poisson", alpha = 0.05, rationale = "Same pathogen.") + expect_true(ready()) + + args$pool_set(list(keys = keys_for(experiments(), c("23", "21")), sent = 2L)) + session$flushReact() + # The old text is still in the input here (testServer does not apply the + # update that clears it), but it no longer counts. + expect_equal(rationale(), "") + expect_false(ready()) + session$setInputs(run = 1) + expect_null(result()) + + session$setInputs(rationale = "Same pathogen, checked again.") + expect_true(ready()) + }) +}) + +test_that("the log downloads as CSV", { + args <- pooling_args() + shiny::testServer(mod_pooling_server, args = args, { + empty <- read.csv(output$download_log) + expect_equal(names(empty), "note") + + args$pool_set(list(keys = keys_for(experiments(), c("21", "23")), sent = 1L)) + session$flushReact() + session$setInputs(models = "beta_poisson", alpha = 0.05, rationale = "Same pathogen.") + session$setInputs(run = 1) + logged <- read.csv(output$download_log) + expect_equal(nrow(logged), 1L) + expect_true(all(c("run", "workbook", "experiments", "rationale", "p_value") %in% names(logged))) + }) +}) + +test_that("loading a different workbook does not trip over the old selection", { + args <- pooling_args() + shiny::testServer(mod_pooling_server, args = args, { + args$pool_set(list(keys = keys_for(experiments(), c("21", "23")), sent = 1L)) + session$flushReact() + session$setInputs(models = "beta_poisson", alpha = 0.05, rationale = "Same pathogen.") + expect_true(ready()) + + # A workbook without 21 and 23, set before the load counter as the Data + # tab does. + rows <- read_qmra_workbook(qmra_sample_path()) + args$rows(rows[!rows$experiment_id %in% c("21", "23"), ]) + args$workbook_loaded(2L) + expect_no_error(session$flushReact()) + expect_equal(keys(), character()) + }) +}) diff --git a/app/tests/testthat/test-mod_upload.R b/app/tests/testthat/test-mod_upload.R index 35eab64..b12e918 100644 --- a/app/tests/testthat/test-mod_upload.R +++ b/app/tests/testthat/test-mod_upload.R @@ -119,5 +119,9 @@ test_that("a bad workbook reports an error and keeps the previous workbook", { expect_match(workbook_error(), "raw data") expect_equal(length(unique(session$returned$workbook()$experiment_key)), 14L) expect_equal(session$returned$workbook_loaded(), 1L) + + # The error is stale once a dataset loads. + session$setInputs(example = 1) + expect_null(workbook_error()) }) }) From b596cb68627e793c83519249000885e560fd6c52 Mon Sep 17 00:00:00 2001 From: seanthimons Date: Wed, 28 Oct 2026 17:26:16 -0400 Subject: [PATCH 5/6] feat(pooling): explore groups and send pooled data to Fit - Add strain/isolate metadata and filters to the experiment library - Add exploratory poolability grouping, logging, and data-informed set notices - Send pooled experiments to Data and Fit; add Cryptosporidium sample data and pooling article --- .Rbuildignore | 1 + README.md | 10 + app/R/app_server.R | 6 +- app/R/fct_library.R | 15 +- app/R/fct_pooling.R | 307 ++++++++++++++++++- app/R/mod_library.R | 5 +- app/R/mod_pooling.R | 216 ++++++++++++- app/R/mod_upload.R | 5 +- app/README.md | 46 ++- app/dev/make_qmra_sample.R | 12 +- app/inst/extdata/qmra-wiki-sample.xlsx | Bin 14413 -> 16145 bytes app/tests/testthat/test-app_server.R | 17 + app/tests/testthat/test-fct_library.R | 13 +- app/tests/testthat/test-fct_pooling.R | 88 +++++- app/tests/testthat/test-mod_library.R | 10 +- app/tests/testthat/test-mod_pooling.R | 108 +++++++ app/tests/testthat/test-mod_upload.R | 8 +- data-raw/cryptosporidium-isolates.R | 38 +++ data-raw/pooling-article-results.R | 40 +++ inst/extdata/cryptosporidium-isolates.csv | 33 ++ vignettes/articles/pooling-results.rds | Bin 0 -> 1276 bytes vignettes/articles/pooling.Rmd | 358 ++++++++++++++++++++++ vignettes/getting-started.Rmd | 20 +- 23 files changed, 1290 insertions(+), 66 deletions(-) create mode 100644 data-raw/cryptosporidium-isolates.R create mode 100644 data-raw/pooling-article-results.R create mode 100644 inst/extdata/cryptosporidium-isolates.csv create mode 100644 vignettes/articles/pooling-results.rds create mode 100644 vignettes/articles/pooling.Rmd diff --git a/.Rbuildignore b/.Rbuildignore index 92416cb..735abe3 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -13,3 +13,4 @@ ^app$ ^AGENTS\.md$ ^CLAUDE\.md$ +^vignettes/articles$ diff --git a/README.md b/README.md index 52cf8ca..4a8cee5 100644 --- a/README.md +++ b/README.md @@ -317,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. diff --git a/app/R/app_server.R b/app/R/app_server.R index f2d2d73..cfa8117 100644 --- a/app/R/app_server.R +++ b/app/R/app_server.R @@ -18,11 +18,15 @@ app_server <- function(input, output, session) { bslib::nav_select("main", "Data", session = session) }) observeEvent(upload$workbook_loaded(), bslib::nav_select("main", "Library", session = session), ignoreInit = TRUE) - mod_pooling_server( + 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 diff --git a/app/R/fct_library.R b/app/R/fct_library.R index 276808f..7e09e17 100644 --- a/app/R/fct_library.R +++ b/app/R/fct_library.R @@ -14,6 +14,7 @@ library_fields <- c( "Dose unit" = "dose_unit", "Host" = "host_group", "Exposure route" = "route_group", + "Strain or isolate" = "strain", "Preferred model" = "model_label" ) @@ -28,7 +29,7 @@ library_fields <- c( #' @return A data frame with the original metadata, numeric `dose`, `positive`, #' and `negative` (`NA` in the single row kept for an experiment that reports #' only parameters), and the normalised columns `experiment_key`, `category`, -#' `outcome_group`, `dose_unit`, `host_group`, `route_group`, `model`, +#' `outcome_group`, `dose_unit`, `host_group`, `route_group`, `strain`, `model`, #' `model_label`, and `wiki_pooled`. #' @noRd read_qmra_workbook <- function(path) { @@ -42,6 +43,9 @@ read_qmra_workbook <- function(path) { if (length(missing) > 0L) { stop(sprintf("The \"raw data\" sheet is missing: %s.", paste(missing, collapse = ", ")), call. = FALSE) } + # Optional: the agent strain or isolate, often what separates experiments + # on the same pathogen (Cryptosporidium Iowa, TAMU, and UCP isolates). + strain <- if ("Agent Strain" %in% names(raw)) raw[["Agent Strain"]] else rep(NA_character_, nrow(raw)) raw <- as.data.frame(raw[library_required_columns], stringsAsFactors = FALSE) rows <- data.frame( experiment_id = trimws(as.character(raw$experiment_id)), @@ -54,6 +58,7 @@ read_qmra_workbook <- function(path) { host_group = library_host_group(raw[["Host type"]]), route = trimws(as.character(raw[["Exposure route"]])), route_group = library_route_group(raw[["Exposure route"]]), + strain = library_strain(strain), preferred_model = trimws(as.character(raw[["Preferred model"]])), model = library_model(raw[["Preferred model"]]), parameters = as.character(raw[["Optimized parameters"]]), @@ -110,6 +115,13 @@ library_dose_unit <- function(x) { x } +#' @noRd +library_strain <- function(x) { + x <- trimws(as.character(x)) + x[is.na(x) | !nzchar(x) | toupper(x) == "NR"] <- "Not reported" + x +} + #' @noRd library_host_group <- function(x) { x <- tolower(as.character(x)) @@ -205,6 +217,7 @@ library_experiments <- function(rows) { host_group = first$host_group, route = first$route, route_group = first$route_group, + strain = first$strain, model = first$model, model_label = first$model_label, wiki_pooled = first$wiki_pooled, diff --git a/app/R/fct_pooling.R b/app/R/fct_pooling.R index 788a1ab..45b613e 100644 --- a/app/R/fct_pooling.R +++ b/app/R/fct_pooling.R @@ -7,14 +7,18 @@ #' @param rows Output of [read_qmra_workbook()]. #' @param experiments Output of [library_experiments()]. #' @param keys `experiment_key` values, in the order chosen. -#' @return The matching rows of `experiments` with `label` and `signature` -#' (the experiment's dose groups, to spot the same counts entered twice). +#' @return The matching rows of `experiments` with `label` (pathogen and id, +#' for tables and the log), `short` (for plots: just the id, with the strain +#' when reported, when every experiment is the same pathogen), and +#' `signature` (the experiment's dose groups, to spot the same counts +#' entered twice). #' @noRd pool_selection <- function(rows, experiments, keys) { selection <- experiments[match(unique(keys), experiments$experiment_key), , drop = FALSE] selection <- selection[!is.na(selection$experiment_key), , drop = FALSE] if (nrow(selection) == 0L) { selection$label <- character() + selection$short <- character() selection$signature <- character() return(selection) } @@ -22,6 +26,11 @@ pool_selection <- function(rows, experiments, keys) { ambiguous <- label %in% label[duplicated(label)] label[ambiguous] <- paste0(label[ambiguous], " (", selection$response_type[ambiguous], ")") selection$label <- make.unique(label, sep = " ") + short <- if (length(unique(selection$pathogen)) == 1L) paste0("#", selection$experiment_id) else label + strain <- selection$strain != "Not reported" + short[strain] <- paste0(short[strain], " (", selection$strain[strain], ")") + short[ambiguous] <- paste0(short[ambiguous], " (", selection$response_type[ambiguous], ")") + selection$short <- make.unique(short, sep = " ") selection$signature <- vapply(selection$experiment_key, function(key) { counts <- rows[rows$experiment_key == key & !is.na(rows$dose), c("dose", "positive", "negative"), drop = FALSE] if (nrow(counts) == 0L) { @@ -121,6 +130,7 @@ pool_rules <- list( list(id = "response_type", severity = "warn", label = "One outcome description", check = function(selection) pool_differs(selection$response_type, "reported outcome")), list(id = "host", severity = "warn", label = "One host", check = function(selection) pool_differs(selection$host, "host")), list(id = "route", severity = "warn", label = "One exposure route", check = function(selection) pool_differs(selection$route, "exposure route")), + list(id = "strain", severity = "warn", label = "One strain or isolate", check = function(selection) pool_differs(selection$strain, "strain or isolate")), list(id = "pathogen", severity = "warn", label = "One pathogen", check = function(selection) pool_differs(selection$pathogen, "pathogen")) ) @@ -151,20 +161,20 @@ pool_blocked <- function(checks) { any(!checks$passed & checks$severity == "block") } -#' Datasets for the pooling functions, named by experiment label +#' Datasets for the pooling functions, named by short experiment label #' -#' Each dataset's `study_id` is its label, so plots keep experiments apart, and +#' Each dataset's `study_id` is its short label, so plots keep experiments apart, and #' every dataset gets the first experiment's spelling of the shared dose unit #' so the pooled rows carry one unit. #' @noRd pool_datasets <- function(rows, selection) { datasets <- lapply(seq_len(nrow(selection)), function(i) { data <- library_experiment_data(rows, selection$experiment_key[[i]]) - data$study_id <- selection$label[[i]] + data$study_id <- selection$short[[i]] data$dose_unit <- selection$dose_unit[[1L]] data }) - stats::setNames(datasets, selection$label) + stats::setNames(datasets, selection$short) } #' Run the poolability test, the pooled fit, and the separate fits @@ -393,24 +403,295 @@ pool_wiki_reference <- function(rows, experiments, selection) { ) } -#' Log rows for one run: one row per model +#' One log row per model, in the schema shared by tests and explorations +#' +#' Columns that do not apply to a kind of entry are `NA`, so rows of both +#' kinds bind together. #' @noRd -pool_log_entry <- function(run, selection, checks, rationale, workbook, time = Sys.time()) { - warned <- checks$label[!checks$passed & checks$severity == "warn"] +pool_log_rows <- function(kind, selection, models, alpha, workbook, set_origin, rationale = NA_character_, + warnings_acknowledged = NA_character_, method = NA_character_, groups = NA_character_, + n_pairwise_tests = NA_integer_, partitions_considered = NA_real_, lrt_statistic = NA_real_, + df = NA_real_, p_value = NA_real_, pooling_rejected = NA, converged = NA, time = Sys.time()) { data.frame( + kind = kind, time = format(time, "%Y-%m-%d %H:%M:%S"), - workbook = if (is.null(workbook)) "" else workbook, + workbook = if (is.null(workbook)) NA_character_ else workbook, experiments = paste(selection$label, collapse = "; "), dose_unit = selection$dose_unit[[1L]], + set_origin = set_origin, rationale = rationale, - warnings_acknowledged = paste(warned, collapse = "; "), - model = run$test$model, + warnings_acknowledged = warnings_acknowledged, + model = models, + alpha = alpha, + method = method, + groups = groups, + n_pairwise_tests = n_pairwise_tests, + partitions_considered = partitions_considered, + lrt_statistic = lrt_statistic, + df = df, + p_value = p_value, + pooling_rejected = pooling_rejected, + converged = converged, + stringsAsFactors = FALSE + ) +} + +#' Log rows for a confirmatory test: one row per model +#' @noRd +pool_log_entry <- function(run, selection, checks, rationale, workbook, set_origin = "Library", time = Sys.time()) { + warned <- checks$label[!checks$passed & checks$severity == "warn"] + pool_log_rows( + kind = "test", + selection = selection, + models = run$test$model, alpha = run$alpha, + workbook = workbook, + set_origin = set_origin, + rationale = rationale, + warnings_acknowledged = paste(warned, collapse = "; "), lrt_statistic = run$test$lrt_statistic, df = run$test$df, p_value = run$test$p_value, pooling_rejected = !run$test$poolable, converged = run$test$converged, - stringsAsFactors = FALSE + time = time + ) +} + +#' Log rows for an exploration: one row per model, with the groups found +#' @noRd +pool_explore_log_entry <- function(explore, selection, workbook, set_origin = "Library", time = Sys.time()) { + groups <- vapply(explore$models, function(model) pool_groups_text(explore, model), character(1L)) + pool_log_rows( + kind = "exploration", + selection = selection, + models = explore$models, + alpha = explore$alpha, + workbook = workbook, + set_origin = set_origin, + method = explore$method, + groups = unname(groups), + n_pairwise_tests = explore$n * (explore$n - 1L) / 2L, + partitions_considered = if (explore$method == "exhaustive") pool_bell(explore$n) else NA_real_, + time = time + ) +} + +#' Number of ways to partition n items (the Bell number) +#' @noRd +pool_bell <- function(n) { + row <- 1 + for (i in seq_len(n)) { + nxt <- numeric(length(row) + 1L) + nxt[[1L]] <- row[[length(row)]] + for (j in seq_along(row)) { + nxt[[j + 1L]] <- nxt[[j]] + row[[j]] + } + row <- nxt + } + row[[1L]] +} + +#' Most experiments an exploration will search, to keep it to seconds +#' @noRd +pool_explore_limit <- 12L + +#' Exploratory search for mutually poolable groups +#' +#' Runs [singlehit::group_datasets()]: every pairwise poolability test, then a +#' grouping (all partitions for up to six experiments, a greedy merge above +#' that). This is the data-driven grouping #20 warns about, so the app labels +#' it exploratory and logs it, and a set chosen from it is marked as such. +#' The exact beta-Poisson model is left out because it is too slow to search. +#' +#' @return A list of `groups` (dataset, model, group), `pairwise` (p-value +#' matrices by model), `method`, `models`, `alpha`, `n`, and `warnings`; or +#' an error condition. +#' @noRd +pool_explore <- function(datasets, models, alpha = 0.05) { + models <- intersect(models, c("exponential", "beta_poisson")) + n <- length(datasets) + if (length(models) == 0L) { + return(simpleError("Exploration uses the exponential and beta-Poisson models; choose at least one of them.")) + } + if (n < 2L || n > pool_explore_limit) { + return(simpleError(sprintf("Exploration needs between 2 and %d experiments (this set has %d).", pool_explore_limit, n))) + } + method <- if (n <= 6L) "exhaustive" else "agglomerative" + warnings <- character() + tryCatch( + { + groups <- withCallingHandlers( + singlehit::group_datasets(datasets, models = models, alpha = alpha, method = method), + warning = function(w) { + warnings <<- c(warnings, conditionMessage(w)) + invokeRestart("muffleWarning") + } + ) + list( + groups = as.data.frame(groups), + pairwise = attr(groups, "pairwise"), + method = method, + models = models, + alpha = alpha, + n = n, + warnings = unique(warnings) + ) + }, + error = function(e) e + ) +} + +#' The groups found for one model, as text such as "{A, B} {C}" +#' @noRd +pool_groups_text <- function(explore, model) { + groups <- explore$groups[explore$groups$model == model, , drop = FALSE] + members <- split(groups$dataset, groups$group) + paste0("{", vapply(members, paste, character(1L), collapse = ", "), "}", collapse = " ") +} + +#' Groups of two or more that can become the pooling set +#' +#' @return A named character vector for a select input: values are +#' "model::group", names describe the members. +#' @noRd +pool_group_choices <- function(explore) { + choices <- character() + for (model in explore$models) { + groups <- explore$groups[explore$groups$model == model, , drop = FALSE] + for (group in sort(unique(groups$group))) { + members <- groups$dataset[groups$group == group] + if (length(members) >= 2L) { + label <- sprintf("%s, group %d: %s", model_label_ui(model), group, paste(members, collapse = ", ")) + choices[[label]] <- paste(model, group, sep = "::") + } + } + } + choices +} + +#' Experiment keys for a chosen group +#' @noRd +pool_group_keys <- function(explore, choice, selection) { + parts <- strsplit(choice, "::", fixed = TRUE)[[1L]] + groups <- explore$groups + members <- groups$dataset[groups$model == parts[[1L]] & groups$group == as.integer(parts[[2L]])] + selection$experiment_key[match(members, selection$short)] +} + +#' Heatmap of pairwise poolability p-values, one panel per model +#' @noRd +pool_explore_plot <- function(explore) { + tiles <- do.call(rbind, lapply(explore$models, function(model) { + matrix <- explore$pairwise[[model]] + cells <- expand.grid(a = rownames(matrix), b = colnames(matrix), stringsAsFactors = FALSE) + cells$p <- matrix[cbind(cells$a, cells$b)] + cells$model <- model_label_ui(model) + cells[cells$a != cells$b, , drop = FALSE] + })) + labels <- rownames(explore$pairwise[[explore$models[[1L]]]]) + tiles$a <- factor(tiles$a, levels = labels) + tiles$b <- factor(tiles$b, levels = rev(labels)) + tiles$model <- factor(tiles$model, levels = model_label_ui(explore$models)) + ggplot2::ggplot(tiles, ggplot2::aes(x = .data$a, y = .data$b, fill = .data$p)) + + ggplot2::geom_tile(colour = "white") + + ggplot2::geom_text(ggplot2::aes(label = format_p(.data$p, digits = 2)), size = 3.4) + + ggplot2::scale_fill_gradient(low = "#F4E1D2", high = "#0066A1", limits = c(0, 1), name = "Pairwise p") + + ggplot2::facet_wrap(ggplot2::vars(.data$model)) + + ggplot2::labs( + x = NULL, + y = NULL, + caption = paste(strwrap(paste( + "Unadjusted pairwise likelihood-ratio tests. A high p-value is not evidence two experiments share a curve;", + "with few doses the test rarely rejects." + ), width = 90), collapse = "\n") + ) + + ggplot2::theme_minimal(base_size = 12) + + ggplot2::theme(axis.text.x = ggplot2::element_text(angle = 30, hjust = 1), plot.caption = ggplot2::element_text(hjust = 0)) +} + +#' One dataset of the stacked experiments, for the Data and Fit tabs +#' +#' singlehit's reader needs one value per metadata column, so the pooled set +#' gets combined values. It then adds up counts at doses shared between +#' experiments: the fitted parameters are unchanged, but goodness of fit and +#' the bootstrap use the merged dose groups, while the Pooling test keeps each +#' experiment's rows separate. A data-informed set says so in its `study_id`, +#' which the Fit, Compare, and Plots tabs show. +#' @noRd +pool_fit_data <- function(datasets, selection, data_informed = FALSE) { + counts <- do.call(rbind, lapply(datasets, function(data) data[c("dose", "positive", "negative")])) + rownames(counts) <- NULL + counts$study_id <- paste("Pooled:", paste(selection$label, collapse = " + ")) + if (data_informed) { + counts$study_id <- paste(counts$study_id, "[data-informed set]") + } + counts$host <- paste(unique(selection$host), collapse = " + ") + counts$dose_unit <- selection$dose_unit[[1L]] + counts$endpoint <- paste(unique(selection$response_type), collapse = " + ") + counts +} + +#' Whether a set of keys lies within a set that was explored +#' +#' A set chosen after seeing an exploration of the same experiments is +#' data-informed however it reached the Pooling tab, so its test p-values are +#' biased toward pooling. +#' @noRd +pool_data_informed <- function(keys, explored) { + length(keys) > 0L && any(vapply(explored, function(set) all(keys %in% set), logical(1L))) +} + +#' "How to read this heatmap", a collapsed note under the exploration plot +#' +#' A summary of the pooling article on the package website. +#' @noRd +pool_heatmap_help <- function() { + bslib::accordion( + open = FALSE, + class = "mb-3", + bslib::accordion_panel( + "How to read this heatmap", + tags$ul( + class = "small mb-1", + tags$li( + strong("Each square is one test of two experiments: "), + "does one shared curve fit clearly worse than a curve for each? Pale (small p) is evidence they differ;", + "blue (large p) means no evidence of a difference." + ), + tags$li( + strong("\"Not rejected\" is not \"the same\". "), + "With a few doses per experiment the test rarely rejects, so very different experiments often pass,", + "including ones with different outcomes." + ), + tags$li( + strong("Beta-Poisson is usually bluer. "), + "Its second parameter and two degrees of freedom make the test weaker; it does not mean the experiments", + "agree more." + ), + tags$li( + strong("Pairs do not chain. "), + "A with B and B with C passing does not mean A, B, and C share one curve." + ), + tags$li( + strong("Groupings are not unique. "), + "The model and the search method can give different groups from the same data. Choosing a group from", + "them biases its test toward pooling." + ), + tags$li( + strong("Check where the subjects came from. "), + "The test assumes independent experiments; entries that report the same volunteers look alike for that", + "reason alone." + ) + ), + p( + class = "small mb-0", + a( + "Read the full explanation with worked examples", + href = "https://seanthimons.github.io/singlehit/articles/pooling.html", + target = "_blank" + ) + ) + ) ) } diff --git a/app/R/mod_library.R b/app/R/mod_library.R index 3ed90cb..c9f8ca9 100644 --- a/app/R/mod_library.R +++ b/app/R/mod_library.R @@ -26,6 +26,7 @@ mod_library_ui <- function(id) { filter("dose_unit", "Dose unit"), filter("host_group", "Host"), filter("route_group", "Exposure route"), + filter("strain", "Strain or isolate"), filter("model_label", "Preferred model"), checkboxInput(ns("hide_wiki_pooled"), "Hide entries the wiki already pooled", value = FALSE), hr(), @@ -76,7 +77,7 @@ mod_library_ui <- function(id) { #' @noRd mod_library_server <- function(id, rows) { moduleServer(id, function(input, output, session) { - filters <- c("category", "pathogen", "outcome_group", "dose_unit", "host_group", "route_group", "model_label") + filters <- c("category", "pathogen", "outcome_group", "dose_unit", "host_group", "route_group", "strain", "model_label") loaded <- reactive({ req(rows()) rows() @@ -176,6 +177,7 @@ mod_library_server <- function(id, rows) { `Dose unit` = shown$dose_unit, Host = shown$host, Route = shown$route, + Strain = shown$strain, `Preferred model` = shown$model_label, Doses = shown$doses, Subjects = format(shown$subjects, big.mark = ","), @@ -212,6 +214,7 @@ mod_library_server <- function(id, rows) { sent(list( data = library_experiment_data(loaded(), experiment$experiment_key), label = library_experiment_label(experiment), + source = "Library tab", sent = count() )) }) diff --git a/app/R/mod_pooling.R b/app/R/mod_pooling.R index fca387c..e76ca86 100644 --- a/app/R/mod_pooling.R +++ b/app/R/mod_pooling.R @@ -37,7 +37,8 @@ mod_pooling_ui <- function(id) { ), uiOutput(ns("empty")), uiOutput(ns("selection_ui")), - uiOutput(ns("results_ui")) + uiOutput(ns("results_ui")), + uiOutput(ns("explore_ui")) ) } @@ -50,7 +51,8 @@ mod_pooling_ui <- function(id) { #' @param workbook_name Reactive returning where the workbook came from. #' @param workbook_loaded Reactive counter that changes on every workbook load, #' including reloading the same file, which leaves `rows()` unchanged. -#' @return Nothing useful; the log is kept in the session. +#' @return A list with `sent`: `NULL` until pooled data are sent to the Data +#' tab, then a list of `data`, `label`, `source`, and `sent` (a counter). #' @noRd mod_pooling_server <- function(id, rows, experiments, pool_set, workbook_name = reactive(NULL), workbook_loaded = reactive(0L)) { moduleServer(id, function(input, output, session) { @@ -62,17 +64,47 @@ mod_pooling_server <- function(id, rows, experiments, pool_set, workbook_name = # A rationale is written for one set of experiments. After the set or the # workbook changes, the old text no longer counts until it is edited. rationale_stale <- reactiveVal(FALSE) + # Where the current set came from, every set of experiments explored this + # session, and the exploration of the current set. + origin <- reactiveVal("Library") + explored <- reactiveVal(list()) + exploration <- reactiveVal(NULL) + sent <- reactiveVal(NULL) - start_over <- function(new_keys) { + start_over <- function(new_keys, from = "Library") { keys(new_keys) + origin(from) result(NULL) + exploration(NULL) rationale_stale(TRUE) updateTextAreaInput(session, "rationale", value = "") } # Ahead of the other observers, so nothing reads the new workbook with the # previous workbook's keys. + # Explored sets are kept across workbook loads, like the log: reloading + # must not make a data-informed set look prespecified. observeEvent(workbook_loaded(), start_over(character()), priority = 1) observeEvent(pool_set(), start_over(pool_set()$keys), priority = 1) + + # A set chosen from an exploration, or lying within a set already explored + # (however it was chosen), is data-informed: its test p-values are biased + # toward pooling. Re-picking the same experiments on the Library tab does + # not reset this. + set_origin <- reactive({ + if (origin() != "Library") { + origin() + } else if (pool_data_informed(keys(), explored())) { + "Library, within a set explored this session" + } else { + "Library" + } + }) + + append_log <- function(entry) { + previous <- log() + entry <- cbind(run = if (is.null(previous)) 1L else max(previous$run) + 1L, entry) + log(rbind(previous, entry)) + } observeEvent(input$rationale, rationale_stale(FALSE), ignoreInit = TRUE) selection <- reactive({ @@ -99,6 +131,15 @@ mod_pooling_server <- function(id, rows, experiments, pool_set, workbook_name = output$selection_ui <- renderUI({ req(nrow(selection()) > 0L) tagList( + if (set_origin() != "Library") { + div( + class = "alert alert-warning", + strong("Data-informed set. "), + sprintf("Source: %s. ", set_origin()), + "These experiments were chosen after seeing poolability results for them, so the test below is biased", + "toward pooling and its p-values are not confirmatory. Say so if you report a pooled fit." + ) + }, h5("Experiments"), tableOutput(session$ns("selection_table")), h5("Checks"), @@ -119,6 +160,7 @@ mod_pooling_server <- function(id, rows, experiments, pool_set, workbook_name = `Dose unit` = shown$dose_unit, Host = shown$host, Route = shown$route, + Strain = shown$strain, Doses = shown$doses, Subjects = format(shown$subjects, big.mark = ","), check.names = FALSE @@ -180,20 +222,26 @@ mod_pooling_server <- function(id, rows, experiments, pool_set, workbook_name = }) result(list(run = run, selection = chosen)) if (!inherits(run, "error")) { - entry <- pool_log_entry(run, chosen, checks(), rationale(), workbook_name()) - previous <- log() - entry <- cbind(run = if (is.null(previous)) 1L else max(previous$run) + 1L, entry) - log(rbind(previous, entry)) + append_log(pool_log_entry(run, chosen, checks(), rationale(), workbook_name(), set_origin())) } }) output$log_status <- renderUI({ entries <- log() - tests <- if (is.null(entries)) 0L else length(unique(entries$run)) + count <- function(kind) if (is.null(entries)) 0L else length(unique(entries$run[entries$kind == kind])) + tests <- count("test") + explorations <- count("exploration") + # One exploration row per model, each with that model's pairwise tests. + pairwise <- if (is.null(entries)) 0 else sum(entries$n_pairwise_tests[entries$kind == "exploration"]) p( class = "small", - sprintf("%d poolability %s run this session.", tests, if (tests == 1L) "test" else "tests"), - if (tests > 1L) " The p-values are not adjusted for running several tests." + sprintf( + "%d poolability %s and %d %s run this session.", + tests, if (tests == 1L) "test" else "tests", + explorations, if (explorations == 1L) "exploration" else "explorations" + ), + if (explorations > 0L) sprintf(" The explorations made %d pairwise tests.", as.integer(pairwise)), + if (tests > 1L || explorations > 0L) " None of the p-values are adjusted for running several tests." ) }) @@ -202,7 +250,7 @@ mod_pooling_server <- function(id, rows, experiments, pool_set, workbook_name = content = function(file) { entries <- log() if (is.null(entries)) { - entries <- data.frame(note = "No poolability tests have been run in this session.") + entries <- data.frame(note = "No poolability tests or explorations have been run in this session.") } utils::write.csv(entries, file, row.names = FALSE) } @@ -238,7 +286,149 @@ mod_pooling_server <- function(id, rows, experiments, pool_set, workbook_name = h5("Separate and pooled fits"), plotOutput(session$ns("fit_plot"), height = "480px", fill = FALSE), h5("Spread of the estimates"), - plotOutput(session$ns("forest_plot"), height = "380px", fill = FALSE) + plotOutput(session$ns("forest_plot"), height = "380px", fill = FALSE), + h5("Analyse the pooled data"), + if (all(!run$test$poolable)) { + div(class = "alert alert-warning small", "The test rejected pooling under every model chosen.") + }, + p( + class = "text-muted small", + "Sends the stacked experiments to the Data tab as one dataset for Fit, Compare, Bootstrap, and Plots. singlehit", + "adds up counts at doses the experiments share: the fitted parameters are the same as the pooled fit here,", + "but goodness of fit and the bootstrap use the merged dose groups, while the test above keeps each", + "experiment's rows separate." + ), + actionButton(session$ns("send"), "Send pooled data to Data and Fit", icon = icon("arrow-right"), class = "btn-primary") + ) + }) + + send_count <- reactiveVal(0L) + observeEvent(input$send, { + current <- result() + req(current, !inherits(current$run, "error"), identical(current$selection$experiment_key, selection()$experiment_key)) + chosen <- current$selection + informed <- set_origin() != "Library" + data <- tryCatch(pool_fit_data(pool_datasets(rows(), chosen), chosen, informed), error = function(e) e) + if (inherits(data, "error")) { + showNotification(conditionMessage(data), type = "error") + return() + } + send_count(send_count() + 1L) + label <- paste("Pooled:", paste(chosen$label, collapse = " + ")) + if (informed) { + label <- paste0(label, " [data-informed set]") + } + sent(list(data = data, label = label, source = "Pooling tab", sent = send_count())) + append_log(pool_log_rows( + kind = "sent to Fit", + selection = chosen, + models = NA_character_, + alpha = current$run$alpha, + workbook = workbook_name(), + set_origin = set_origin() + )) + }) + + # Exploration ----------------------------------------------------------- + + output$explore_ui <- renderUI({ + req(nrow(selection()) > 0L) + blocked <- pool_blocked(checks()) + too_many <- nrow(selection()) > pool_explore_limit + button <- actionButton(session$ns("explore"), "Explore pairwise poolability", icon = icon("magnifying-glass"), class = "btn-outline-secondary") + tagList( + hr(), + h5("Exploratory: which experiments pool with which"), + div( + class = "alert alert-warning small", + strong("Not for choosing groups to report as pooled. "), + "This searches the set for mutually poolable groups (singlehit's group_datasets(): every pairwise test, then", + "a grouping). Groups found this way are chosen by the data, so a test of them is biased toward pooling (#20).", + "Each exploration is logged, and a set taken from it is marked data-informed. The exact beta-Poisson model", + sprintf("is left out, and at most %d experiments are searched.", pool_explore_limit) + ), + if (blocked || too_many) tagAppendAttributes(button, disabled = NA) else button, + if (blocked) p(class = "text-muted small mt-1", "A blocking check failed; see Checks."), + if (too_many) p(class = "text-muted small mt-1", sprintf("Choose at most %d experiments to explore.", pool_explore_limit)), + uiOutput(session$ns("explore_results")) + ) + }) + + observeEvent(input$explore, { + req(!pool_blocked(checks()), nrow(selection()) <= pool_explore_limit) + chosen <- selection() + models <- if (length(input$models) > 0L) input$models else c("exponential", "beta_poisson") + explore <- withProgress(message = "Testing every pair and grouping", { + tryCatch(pool_explore(pool_datasets(rows(), chosen), models, alpha = input$alpha), error = function(e) e) + }) + exploration(list(explore = explore, selection = chosen)) + if (!inherits(explore, "error")) { + # Log first: set_origin() describes the set as it was before this + # exploration. + append_log(pool_explore_log_entry(explore, chosen, workbook_name(), set_origin())) + explored(c(explored(), list(chosen$experiment_key))) + } + }) + + explore_ok <- reactive({ + current <- exploration() + req(current, !inherits(current$explore, "error")) + current + }) + + output$explore_results <- renderUI({ + current <- exploration() + req(current) + if (inherits(current$explore, "error")) { + return(div(class = "alert alert-danger mt-3", "The exploration could not run: ", conditionMessage(current$explore))) + } + explore <- current$explore + choices <- pool_group_choices(explore) + tagList( + p( + class = "text-muted small mt-3", + sprintf( + "%d pairwise tests; groups by the %s method%s.", + explore$n * (explore$n - 1L) / 2L, + explore$method, + if (explore$method == "exhaustive") sprintf(" over all %d ways to split the set", pool_bell(explore$n)) else " (greedy merges)" + ) + ), + if (length(explore$warnings) > 0L) { + div(class = "alert alert-warning small", strong("Fit warnings: "), paste(explore$warnings, collapse = " ")) + }, + plotOutput(session$ns("explore_plot"), height = "420px", fill = FALSE), + pool_heatmap_help(), + tableOutput(session$ns("explore_groups")), + if (length(choices) > 0L) { + div( + class = "d-flex flex-wrap align-items-end gap-3", + selectInput(session$ns("group"), "Group to test", choices = choices, width = "560px"), + actionButton(session$ns("use_group"), "Use this group as the pooling set", class = "btn-outline-warning mb-3") + ) + } else { + p(class = "text-muted", "No two experiments were grouped together.") + } + ) + }) + + output$explore_plot <- renderPlot(pool_explore_plot(explore_ok()$explore), res = 96) + output$explore_groups <- renderTable({ + explore <- explore_ok()$explore + data.frame( + Model = model_label_ui(explore$models), + Groups = vapply(explore$models, function(model) pool_groups_text(explore, model), character(1L)), + check.names = FALSE + ) + }) + + observeEvent(input$use_group, { + current <- explore_ok() + req(input$group %in% pool_group_choices(current$explore)) + model <- strsplit(input$group, "::", fixed = TRUE)[[1L]][[1L]] + start_over( + pool_group_keys(current$explore, input$group, current$selection), + from = sprintf("exploratory grouping (%s)", model_label_ui(model)) ) }) @@ -256,6 +446,6 @@ mod_pooling_server <- function(id, rows, experiments, pool_set, workbook_name = output$fit_plot <- renderPlot(pool_plot(run_ok()), res = 96) output$forest_plot <- renderPlot(pool_forest(run_ok()), res = 96) - invisible(NULL) + list(sent = sent) }) } diff --git a/app/R/mod_upload.R b/app/R/mod_upload.R index 9bb8df1..f3b312c 100644 --- a/app/R/mod_upload.R +++ b/app/R/mod_upload.R @@ -20,7 +20,7 @@ mod_upload_ui <- function(id) { actionButton(ns("sample_workbook"), "Use QMRA wiki sample workbook"), p( class = "text-muted small mt-1", - "14 experiments cut down from the QMRA wiki workbook, for demos and testing. Data compiled from the ", + "18 experiments cut down from the QMRA wiki workbook, for demos and testing. Data compiled from the ", a("QMRA wiki", href = "https://qmrawiki.org", target = "_blank"), "; see each row's reference." ), @@ -135,7 +135,8 @@ mod_upload_server <- function(id, incoming = reactive(NULL)) { ) observeEvent(incoming(), { sent <- incoming() - load_source(list(data = sent$data, label = paste(sent$label, "from the Library tab"))) + from <- if (is.null(sent$source)) "Library tab" else sent$source + load_source(list(data = sent$data, label = paste0(sent$label, " (from the ", from, ")"))) }) # Same reader settings as singlehit::read_dose_response(), so the offered diff --git a/app/README.md b/app/README.md index deb2b05..5f9247e 100644 --- a/app/README.md +++ b/app/README.md @@ -40,9 +40,10 @@ depends on the app. The interface font is overridden in `inst/app/www/custom.css`. - **Library:** explores the QMRA wiki workbook (`.xlsx`, "raw data" sheet), - or the bundled 14-experiment sample, loaded on the Data tab, and is the + or the bundled 18-experiment sample, loaded on the Data tab, and is the starting point for pooling. Filter experiments by pathogen group, - pathogen, outcome group, dose unit, host, exposure route, and preferred model, + pathogen, outcome group, dose unit, host, exposure route, strain or isolate, + and preferred model, and draw every experiment's preferred model on one log-dose plot, faceted and coloured by any of those fields. Curves are refit from the counts; experiments the workbook reports only as parameters are drawn from those (dashed). Dose @@ -63,15 +64,40 @@ depends on the app. `fit_pooled_dose_response()`. Before running you must write why the experiments belong together, and every run is kept in a session log (with the rationale, the warnings you ran past, and the p-values) that downloads - as CSV and counts how many unadjusted tests have been run. The log lasts the whole - session, across workbooks, and the rationale clears whenever the pooling set - or workbook changes, so each set gets its own justification. Results show the test with the caveat that "not rejected" + as CSV and counts how many unadjusted tests have been run. The log lasts + the whole session, across workbooks, and the rationale clears whenever the + pooling set or workbook changes, so each set gets its own justification. + Results show the test with the caveat that "not rejected" is not evidence of equivalence, each experiment's own fit against the pooled fit, the spread of the parameter estimates, and the wiki's own pooled entry when the chosen experiments match one (21 + 23 for Burkholderia mallei - reproduces the wiki's "21,23"). The data-driven `group_datasets()` and - `poolability_combinations()` are deliberately not offered: picking groups - from the fits is the forking path described in #20. + reproduces the wiki's "21,23"). "Send pooled data to Data and Fit" then + makes the tested set one dataset for Fit, Compare, Bootstrap, and Plots. + singlehit adds up counts at doses the experiments share, so the fitted + parameters match the pooled fit but goodness of fit and the bootstrap use + the merged dose groups, while the test keeps each experiment's rows apart. + + An **exploratory** section runs singlehit's `group_datasets()` on the set + (every pairwise test as a heatmap, then a grouping over all ways to split + up to six experiments, or greedy merges for up to 12) and lets you take a + group it finds as the pooling set. This is the data-driven grouping #20 + warns about, so it is labelled exploratory and every exploration is logged + with how many tests it made. A set taken from it, or any set lying within + one already explored this session however it was chosen, is marked + data-informed in the tab, the log, and the label of pooled data sent to + Fit: its test is biased toward pooling. On the sample's Cryptosporidium + experiments the exponential grouping roughly follows the isolate (Iowa + together, TAMU alone, UCP with Moredun), while the beta-Poisson grouping + does not, which is why the isolate is the kind of grouping to decide in + advance rather than read off the exploration. A "How to read this + heatmap" note under the plot summarises the interpretation; the + [pooling article](https://seanthimons.github.io/singlehit/articles/pooling.html) + works through it in full. + + Neither the checks nor the test can tell when two entries report the same + subjects. Experiments 108 and 139 (both Iowa isolate) have identical counts + at six of eight doses, so pooling them, or sending them to Fit, where counts + at shared doses are added, would count the same volunteers twice. **The checks are a tunable pinch point.** They live in one table, `pool_rules` in `R/fct_pooling.R`, and are a first, loose cut to be revised. @@ -79,8 +105,8 @@ depends on the app. that cannot be fitted, double counting (an entry the wiki already pooled, the same experiment twice, identical counts), and dose units that differ or are not reported (singlehit itself also refuses mixed units). Different - outcomes, hosts, routes, or pathogens only warn, because the wiki itself - pools across hosts and routes. + outcomes, hosts, routes, strains, or pathogens only warn, because the wiki + itself pools across hosts and routes. Hierarchical (family-level) models, bootstrap intervals for pooled fits, and export are planned. diff --git a/app/dev/make_qmra_sample.R b/app/dev/make_qmra_sample.R index 9c3110e..a837b3f 100644 --- a/app/dev/make_qmra_sample.R +++ b/app/dev/make_qmra_sample.R @@ -2,9 +2,9 @@ # sample workbook" button and by the tests. The full workbook stays out of the # repository; this keeps only the "raw data" sheet, only count and metadata # columns, and a handful of experiments chosen to cover each pathogen group, -# several dose units, both preferred models, a wiki-pooled entry, a -# parameter-only entry, one with too few doses to fit, and one with an -# unsupported model. +# several dose units, both preferred models, several isolates of one pathogen, +# a wiki-pooled entry, a parameter-only entry, one with too few doses to fit, +# and one with an unsupported model. # # Run from the app/ directory: # Rscript dev/make_qmra_sample.R "../QMRA wiki data_v4.xlsx" @@ -25,6 +25,12 @@ keys <- c( "63|Enterovirus|infection|exponential", "1|Variola major|death|NR", "108|Cryptosporidium hominis/ Cryptosporidium parvum|infection|exponential", + # Iowa, TAMU, UCP, and Moredun isolates: one pathogen, host, route, and unit, + # where the isolate explains which experiments pool. + "139|Cryptosporidium hominis/ Cryptosporidium parvum|infection|exponential", + "140|Cryptosporidium hominis/ Cryptosporidium parvum|infection|exponential", + "141|Cryptosporidium hominis/ Cryptosporidium parvum|infection|beta-Poisson", + "183|Cryptosporidium hominis/ Cryptosporidium parvum|infection|beta-Poisson", "46|Giardia|infection|exponential", "48|Giardia|infection|exponential", "252|scrapie|infection|exponential" diff --git a/app/inst/extdata/qmra-wiki-sample.xlsx b/app/inst/extdata/qmra-wiki-sample.xlsx index 5a70fc8eb6e58e34ec9419751e78374055de7392..1d8cc1b7c8ff8b14da9bf2e7db65773b76722bbf 100644 GIT binary patch delta 11816 zcmZX41yoc~*C?SNjkGiZf)1s0BOL-FDcvb0ofo8A5kx>hLO^Qh8enJ?kcL5|bLj4R zcl7($_ul{3nzhe9`|LXV?6ddGz4OcELtuFeCZ4pVg{D3o8rm5)`p2J`41kf>w7@61 zxiwwy^#JqWJLnJlJfUiQ=kG@~$i;bg?G2)ur)n!%JPOw41{1~BhlMF87Y(*L@N$qD zPbRuK_qwuf&exZ!0K=K-!1|l31C{^}--E45Q)2QH6)%{k4P7|n=Zt`#I?Ih$lQ!C{x%n`?v4b1m=-cc07g zV`ts0ThIzOLNJSK*9qR(a(;Zi?JTI9vP&N=5#TS8k@;yYi%Zg9;`~`VLeQ_f{$?eG z$Qe3NHg|J9I=%+u4w*I#XllB2xjsL+j+k>&H*WSmT>w`u5y0Wr>PAWm8=z(aS9+Ww z1{^OgOy6EFY@Gg4xV_%mIea`PwPYzdaaR&e@lP-feU zXjEic%+Kq%Z)Hw$))4`}bva%@oO*Gqj~qa&GK^-8T6|sIt}m<7Bift+SQA`-eRE7L zb4KE7{0|(Mn>7M1o=snyH2LnI+{{ZBND|ij`MkYtyxHG4>N*5K^40b&-M*^W!}a@f z-OP`~7dhykZr;DWpYX%C(`!FcF5ygqUgIc|Ve!6p!r4>$ijDhg`;qwyO*@89VTF<& zgcZdXDFaNNuYny}bHs6uTcpvX%o;JQvIivE0d*jBoO35Uzm2)PTD>?0o{By009r-F zxoV1Uel59FZeEv$dVU*Zsdl(x&?I`do7Y)(eo8Vme)C)!HXx~Nt9xcB4XcH+R{y?X z$u_EkixJwuoFr9j5!=$zFehxBZCpCmYDl$`tu8Vy^qczNbb-wyU+bKIVCV4&FR3{1 z@c#nL*%^n@(xPi-SW@vMQ=zq~$W*V7Ff3M}pTbztmhzafV@gikdfur|R(LC9?k;Eb&|#RcT## zm7&BOS$!wyj;bZZdd~>!@y|c5noQh{r^%%BQxcX#0aYh{Dx83%aw`9%-WDl z_o<$>o86u(MQXG}K7;Ml2@&b_>c016oSyBw9*PBunmvk|Yl?#ox4cagW5&JeB9*y$ zH7D?{!9=XQ%cVej-EFg4!d7?Ijh&*uIxW#B))bj;9Zv`;8d!_>$a(qhvZPn6GjsZ-EW z`U%KUB)ixMZ5PR|m^j4Bv1Y;jRx^Q>vrN(`qMVmcl!ldKX<4TrQYi12aKh=Y z&$}5|i=x@*YK1K*pRd#cE>S)o&Eeejmc}dteFn5ya9QjBIUSt{I`@#);h$MNb<78u zL*Z;=3yEAYq+!xx0Z$jRS=GWh)W}5kLtlwylOqlrC&$KqSxY}Hs|jXTN+)H#5S>=| za4TU?Epe>Fvf-+)hAm5Q$46Qw^W_H>>L_WYO#Q6FaM^4&b@pD;w725Uj70#G274cB z-<~2=1L3kCNG0v5BkdKn{yZP!b6F)7ImF3m$O|sS$vAhPVpWcC&w1qjA8MvxwMr(5 z#(!AaT+}SPh!Vrq#Y<*DmPEoRY9UbeJaoWtM3GG7@bNL2SK#=MG)NQqpBd1Q*QTIE z?VkjOpwKgYC{;y-#eVF^uZ-#G-fJPw?3DEMh5`wVXUog7u9AX+?tMS1NLlSXb6AYN zusdzY+i=;jL^qw9JUC8_C2RmVi}Wr2Qx@v@Ovk{ipQ6Tk6;&MnOf>E2u6v4c2OtIoApUT1CJzZz+c*?j2SzJg2nzh%9Ve1!)HTq zGK43bE}}624Zknv|BbYj2a4;PdJ_KcFYKsDyU_^~9}AThs3qn71RW}?xw%wPTJjR# zBJDXU1QlEARFVGwZjVc7kpFkb4W1X@LX-kBd10+ zB?!KNhYW5LJPF*4^D(bg5lco80T~nsW%*f(m_SuEwL%)G+&Gwx0P`t-K@-ZP;djx$ z7R(D|F>DUof4q^`gz3#GH%>`*{Ei*)-#Ja)e1GD6+MF-2OvEd4_#mSp{*CTgO}_kr zIxo1B<5PB}jOduulC$Ry9qL^bgN8OMUaM&n=o;YhvHdG!;OU)S_Iv!00?XS_>L9Y; zWBDKI;K-;vT>~Xkq29^02HaQmz#5^$_GFvMe!hPVCo1iIAr2~b|Afav!NWxG5%V*l zD%*e|c+C^Bayp}7;_s_7L{K-3(7s%ihjlnW?C0b?%n|V z-;8S#?v0h*KVe7Z+OCEcOlj&$LU~OLc!;tZ3)ZeDsNk;!W;GB#h(x8Cn<^%lSbwon z0@)~XuMXjl|4B9!>JbAkL>XN2yyyFg4eTwQR6NU^#qmD4RbbSGJh-CfZ?c=)^k)%`@7kM zQHt4Ia1wwDlTUzybjD$1A(-h9673I#(Wpef3_zWD(L>OCegau-mcJNSQ9u#!sK8C2 z_S){%+z7S~lnda+^8F%kyYepH-*-b12>7rX!(8+573ASJn->FQo@oaA%P;N3M75E> zN+(mySqMwp66GJW*8Hl-abL1)uJC>k#Q$=Qw+KvDJ%ToSe4WjZapID2Q+_Nwnt@WG zpm8dWu{9mq*2J>6CQGAyLEmI3+3=R)*iO0^+x&Fsg^pqa)*jF`=Am=a34cvi2Vk~V zNyrDRjTytky|%-plf#vh!@bCVh$J3wn!9QMK=48^6G0R6C#WP?WpEP69R8+4{Zsjc z&c?^EXX3dp>Ugm#hw=?v86i(!BS%Yx-hcB;j@Z^4E>xiwv`psP`VoETYjLQIS&5}4 zZo%LyUjSYe2lL8&iP!Ao0I#x#s0AZlX(Ksy<)cS<&8U&(M0|&vg)qP;4Xd)8CS&cm zdzGK)2X*~1$2GEL?M9G-d*D-d*TP7>h~agnd$E27S6up1Sh#qu86vxxZKEkUBS?i(0%rpe&blp%HFU$=<{!Kr$q&lj*3otO(Iz}be1$R{r- zMb2Ug0po~xck+BZ&L(o~MwbXT)ecokd+(B&H~LPhob+w{#QP1;3`s7@FLf;D^iDZ) zV+KtOueRQ-e_n_;UUE@^?r*td$Jd$PX2f;f$fUrLyl(%WiRWxIrcvCww5W4eFxUR>LWo9Ve?VQ>i=CUunsil2>-b4) z4gU1i+h`ZNEftbEeOh|yP|h6{c-Kt#kx7>Ro1?3*1^9yxfy7^GS2m;q)m2s%RW)BT z>&$KX%z-uimxB(jdKU}cRhpv46bXCO_bynU8b|Fb$5NiaxUU1fXn5z7skpbiQ}`B> zsi2+S*|~L#K`9a)qjb%k_rbYY)9-!^V>?P?OcnTxSVX6{YAmGyocnqeCCUHPc=CqJ zuG7K133$|vC)Z4l-NbSd$HGnVY9CvvV=;%2ZawXDnH{AOrY5}nJ3?rLpVw&e9{h-A z!LeTo(eIvJn8Jsm@s5V&pyEJsDW6A4(xY%P+`XYdH{~oKn6CHLVO9(heiP3gtS-4C8`B7Hv`d}Z=rdPTKyedZRLE}oQQ`Mf z+^HrE!#=KhtfD|OZK!K#F^}OS1b*bV;@FR-YgRurrb4~&1phB);2NX?n9oTcN-ute z@e#QE9NI4}<}tDYx2+2=F9x@*oyz5YtO9qQ?Z$W4&c_p7V~qpH;%7#M8(_|`R_x~c zR=kguRHSA|jva*5_}s1E_#ev^&MOzCoi~*OCEgzO7j#EDE3i-wFcw?|YLtY5PUAOs z&| zb0<0v#Iscyo;6!U%M*ubN{!v>z2{c!*~$TafWHGv{NTG3$oA<;+J&Hk-|#@nFNoOv z8+S?*Oe+(;-ww3rIpI6J5JLZB*hcA-z2wi(PXE8=&y@Rh&09zg!xCi2$DDUr31<6@ z-RB`sP3%BX8#ksEpmVeK-IhnI4(2T`pmjoj)qX^$zbeq5^5*U?4Pri7 zg6!x(fzwmi>=4&q6?R6gciBE+y!3%7CxAwRDP?;#j-uWM!^Fo40()bq4d z5RjbiND4!cvqbt^f~JN|zdJEZ54l>g7F){e5keTB?w9umtw@-rqW9tJS-1YR9b|>Y z!3$DIg&6<|#ua73H9Xi>g=dxB-KGPsCal$#@&<&^Vhy7)-crrczCIN+3JP7 z!BUZVZ7OIxNWr7oK32Z0^kk|7#`pjf6(!ol&9MuLeh`qTS`5E)K?tJmiqQm?y-sNTcP0Lp}_tvh;+#Kb} zY1nmhcgyx!g4mgl8DKK`A5z!CCNtjEJB+NDx^@is7Rec4#khxucLh$qwsLdqshG@D zQk;d|ZR!5NgBkLBV#p1QZtTc_3X#bzUzZjBw$EWPw&;Ca)KCwa${s z1Ggbpy}0nocEltD;{li*#^8buYQZYZnfH|YLo(kL`YluZNzpZFWn~8xMjJ>cbk$1; zuj~Mz$UTVY)s)v+GF_lm1nR1xMw8zQx}G4ZsROGjP?)!#@YyMuHE@P?e^Y8v_tP1O%SD;I zUPJ&boluT_h>Cw4vk7Xv>Lr0!qDZWS>Uf{{7Uj}zBG+|5Mk5k9Axj{+!sejaWj9#r zyKf3bW6}@D=O7Bj{_ulcPBpNOd^`+GU(wP}@^+NDBn%U&K2=23xTU>SxjzfEC6fm1 zLbDtv&7V1pvHb}+0!#aKt$>y%huW#on?vK8CkG(X>%R~voNRdRKqImjCVr+vdp0sj zSFFY{9k7n+cfIh6;*$P^LVZAYrVG)!h)6!yPlf&kH=`z?f|_S0ki#(1RW1Ew&;{%k zmkwI6YOIgmx#~#(gE{aAQP1jeAs|}iTw#Amzrg?U45)c(9Jg;Dt83A8k??N3C=x`d zm%I8k1S+t98b4;W` zOaNdTRZuQ@1SX=+iovtjTXm-tm)p(V11aQ}&SM-WbA@(WL9Z)r<$syOVJSk>;IqGt zK=RRFYFluPD(^j-8jty1uvfBnLr~3j1JA!0_ur(x2Gt`;z-zy%=h~{WSYC4-ePOJ5 z>u#vD<17FT(wz!Jw9Z0Fm%!Hg7X$=2YokM>_S2}ir{};H2~z(@cDXnKvVI7pmNqa>Y~be6i>n-yiYzXwsWD z@NjE!|CsG^DCfX}t=iEn&cEjDxw!1l3D)Y*mo+>lVdqc!DNpC1-v)KfT5u1wKXXru zftEq~;)x$v@Yf||7MAXN4SwL`qqKLNFXl0^qBH?O?{!v-pP3||FJfKQoSYYdNeSj; z*q=2$*LjdW4LFYnbI;hZ`%C1&OzeIx6xyN49$cjER0Jjs|#2^dW?a#kj z{M3Zt+`eDeqNVD$8_$#(E?(%Eqt8Ef1!7X@77FaNgPa>1M$f%v-kAO zTiq9L?p9_uBgWuT0&g$Vo1Zvm*uu5~XeKYaPg*LPQ|d9P1n27P&?Y;BDq89`(5QBN zpMLAlf4?1`MTg% zXO2A@JpU+Bf7jiBKFpS9hh12p0+aeozFbAXK|b@s?(f$%%s43 zK0RFU=ZN;*R?pX00_EUJphe8w6HK6+0zQsznv8A{=L}Z&u`_*I^kRG}hpTRC6Z2i> z=AnuuQ7T7wAc?OfL>uc(O1Cp++??K$;KZ4uRjb}ood?=vZ%{=KD&W5Uobf@W_*8yZ z-Sl&hb~0&s;a*E=Oo zO*sx|V*F{urdK56uZ@sHFPHzOUC=hk`MSI50Rx6qznjgA!+WXY*u6W4T-FlhZq`0% zzcw@Yeqq|{(CZyen^hs`Y^yGU($(iApYao!KsZl!lPNC!(%?ci1JG;Y#^71t9#N1^u1{#ZHhkRc>*)TqU(o7;k=oLpK+kx3mY8&1 z-&w(E@BvEv(@mQCjn)xw@vanNck0wt#WcPPj<{S+OY{0Xdce<9;Dgs+ab4#Nj7$!u zVT%{L?p0A^KQDI07S2fLHoy4B@K8|ij`uxcVPmYJCLn^A`_A-Q!15~F`=xJa(1afg z=0Pb6f2YT2#k-7b$IlCWcIFho?I`l^NC`@!&yO09A+>(+o9X(?JPA4afu1 zq5X>>gaXnfthJBX{Tv7B{h@V3bVz(41fya>wBHO0@1|QlAno7QhbR_=biIo}OI-q1 zFn;;<@4sscn6+{8c_(cAZ13U+Ef=-7kLBN~VEF~1n@9^_H{RqXsc?Bg+u(cnQ*PC> zls<2??VwF2sDB^h-KzMiH_n2@F0XE;X*}iqCW&VeefNy%Qz7a+JrEq7IWx7MggGW$ z-iTrx!{+nlhi)(mvXOODxdkR~fbWLv^Qb>zzuiu0LMC@`19u-T7($N^!p6FL(}09c z!)>72Qn}l2hb&`){{6*w2xzdd*lic-`9jpynQY7Xm42%lgR#FzyHokV9xKc(QAvl-Xbtv-CseBJq=emuv=z|1KQA^|4zx+BE=%oTL6b=2} z2^qiJyHNa6!rYlKfu%@yHn(i+%2Fo&-0+W_AvLFLp5`O9Lq)@fO99p^AL`ww*ap;S zB30BmWE5wa)HX3)Jud*MFDt?IN1e|?$FUb#h}3s0W9YZO=tQIom(2P*1q<~>bjIm0 zz4Come{{#Sz6YM{6{&)`I6~SfUV+8rlOQzB{zIK1N=2Rcw39k?Tw&4jK9hgZhFn7D zKFJ*&XNfvRL+Qm#@U>%OqbRHi=6}3c*JlE~dFyky88KKlC}Avn{r2==`KTCt#ktx6 z-Mu~Ag>BB!-#7^SoFFz~uuXN~YW;eA1y-DXDfx?SFtDMe<#LH=(4^V>2KAnke&BEe zK(Bg28-@;0!fv^`3DM9#KSM(!LqkLJe#7l%YxdU4Qq%pdos*3lm$&1aZ3C6J(?Ycw zr%#ACkhqX&9Fi**rIZ!ifLFC|`g9cEes@@wO8oIck);(&>gT#f=~uR9^tf59S;8aF zR}df9HmxVgruv>o&abGI&hEEfXZUd)Ev^lEgq5X^0EW>SobN(7J7t+vxbj5}7zT;3 ze^ZQvMFDy*}`RAl26?i9yRE}|Dh;f-r{?H zixXMyx6G&38@)s!rbt9m45!u)ozEGZiEHN+_IGa$`c7c^JU?rVTZXPy@D3}ptsfP< z@PAThrU`ss<@P^32qVf!z%Hq)8YIwcvTpZZT=88)x-okcQ@yI>f75(AsHip+JoW;5 zx}}!NgY0MT%^+&ObEM(2)bGn#^W1t?h(Sc-O?#3#YmwE@qT?MeY^8YyI+ixyo`~+f zG{B72zA%TM_k}QnxK-l^>8e?xtah(d_Ai7&p;ds&_l5<8+NGbsw}9=nk(X;6(W(;3 ziBQ40J;=@X#l_zH#fkEPAKaS*On%WYdoYu8`xxY=mNd*yDQwaaXd@h(B{J^(_Wh`} zsI1gz)Pvprp*0E=YV|AFG2alU9!6FCI5Y+75>}W1wmbN22gdI1BokfEfYGSZFN%q>=WmO2{@<{mJ zyAXBYB4v9cFdIG|{myDKjGi|2B&H`tydDuKodXZ`18`0a1H1?~ke@8>e-M4_&eZ|X z@ot^h^ZsZ){c4lsI6br>px^_&e#;5j@bw}QjzI5IaS__8YO=_$G3}MK?R9I7soY!~=cY9cP^b2k^#Fign7~;`{ zg>8{bqThIeY0usQYyKMk*w(82&S!OCk?Gl0(J}H}$Srq+d01Om$+cXJdh&u?*HyUN zEyumZNLeM?_l%4zUa`w_&z%=^+v8ams!HbSzMi)!fFXQk(8v>!9CV%(z zTz^#kL(eON&6en}ec_%cf~z`A_)thqJ|yVtsgjEq9#C+OR5%SiU3+ZhMdl&}u=M+@ z2@~Fn`oZY(f`7fLxSt>EyWa0;)re@)9tBStw*8H6i!UqzG+r{K+}89W9mx8@J>y{i z>WEqM>KJkNwm#~uBl3275liBOo9Is`J~v?@S7PqrB1`ynPepubdy)C3o?~=~aMAp- z22*6nq9#{~Fx3aqBc{pu<+8dkx7MK{}@@c4kmdF~o^N}6@N=$jwCAVGBBUtX~ zs-6}Vp7zIB7dqkvKrH$ zMSW2&`Us8QHY_$r8exdphe0{#n1t6+y}%fIsiiyl#@MZcPhLTRigO1|yP972%foc< zoX2o$9tQI2VaG1U;;GWqMc^)Z@%q^2c@l5BDr-YTyK~CpX2XZI=yV-%6y54?C6AOO zj;h2S{vNbPKh*Onv}|^J>w;Y5c=1l`U}Y^zYsIRLJ8qq|Bk8?ARRD0pyS-!avbNIl zolak8L};GisgF4}p&m{K#p2coL4kT;0ozkw^YGkg^#PZ}dNvy&Dqdjfy?a0a9c8Jc zaEEG5#*72AgMR;eyB@;c5`%zmg$_e3eSRlI^0W(P0vpepneWnu9$(o^R1J5_o=* zT$S2kTk)%zzql9v8&i~W$1I>F_i(Pk6fSY=7PDIRHX`8rh*#MLKqWQe;5*k)M4DUH zU~hyAvy?tLksYY5|2(XAJmD9lB_HoBHn(e0jPqzXS87uc_lmn`a+I>!e+e z`i^FYnER|oDt2Xtj@BPFo*%ks{-5Psoj;-V>(98<#|_(JL}GJvR^IIOF)xuMEY+Gf zE<{s*Uqz-t^s~VyIZDKS3__|7Vq0k`t?qh;@(t7~tdqGsAM|7t^+JP z@+e0*f7JIyC0SK{Y0S%8ry@h^xvzA^NF5%+*ncAuSt|H5e0fd!(;J(IH2i3c8yGAK z0=JbPbD^}f*D1i@b)HWttrhPB;0eLD_Tt#Fuk$ST+XqJ_-e zm@J0T9a~TRs0wl!j8t2`7FE-e?+wvwO4D?QCGXxp8G=@Z;d_m-_6DUm<#9M(-!nWg zU~+!K4BIbIVe{AN7s^wHGJdO7#vp9Stxl7fwj*ddE&$kI zx#Q@1*W-m2aQKb@nZ!i$HI6c3-sNE~_Yp$1qEGg``-c{-s zPYN4EOBW<~^fPFA4!@tfM+oof>;89+0d=@swL-|`3II~c?kcw$-|H(McKr?R;^tto zB*?BkeYl@@yQE}dk?T8D$wlVqtd$&7PLywETr*7(%0 zzf>@Ihorx)x%@H!`F+>(M*w#oM*&4ZR_TPe*DrEPYTB041TGU?S{z*bwGo^22Qn6R zk6*Y9Xaa=Jq4N#*9jvPp_o8iP>F!CrGtaTaWf^6EU0$h8L4A$!AqcN%6`4f-DtW#t z_vrP#&Ocbb#dbaXkQot$Pu(3H0~w20ZLbf0MA+m88_2ph`6he_(|CWO@i{EDPfX12 zb+a#-X<~zAA`BLnwKL_-A(9+pXEKe*}3+ z8Og8RMcb+VW#;!DbD<~m;2#7fA%)xXiNM62=EI1>!Tj(BmHQZIefus)xqbw@-kaIt z_~apVuk<7D{MgW)(n~s?(qYuM?;{|xKHW66j*u85i(j~#8c^TOcH&CeB%)8mE1$De z`<|K}cu+EHko(%8cjtuudH%|-kxzD6TUCEE><=HB+V6(-%q`%)W4vdA3-vrP+jB9+ zAlN;_z4N@+vAUk($dk>}knMUK$cH}^*y2m=l7rcGG)*wA2i|M^PUF)fY!ezi zt#Y#=PLIEKTQV+9`pvVS*l2l#7Z7lBG<)OMac%mv9ktJvIF$%Mi$*P`ukSTX&pg_D z*)^QwZDIEnyPX9`imH`ofw&`edfYI(r_4XQ6z(eHkPz!)J()*G4$#I&-{eMo`rZl-9yE@Ejy-n{*&*NcZ$5&dDr!*y{{lPUbV;9k_ zT)2Dpl=8GuQgQpo5Q2CY%Z}&jies=Ygb`5u*@{!>M}sk?WFz+KWgNTab`swA3=IY@ z;bO}yqM8bVXwRz5&J9aMoD?0DJUllz*J?)h%?eH8ZNOvb=Q!C(nI0{&(oZKv+cKI@ zl%+hrdGNTLGZaAvoV@qD)GBI@lxwZK+XORyzXltVBvMbg>ka<&W*(`dg?nowi7+6; z#w(KCH66rQ^8N%OPky&nTWyFa`8?{DhhOs8gvR}{mO<|>TK2B~Q@M!BqGf{lKJmv? zBKgWcY1|3@9W^z~INogAHx%cg+x9$MfB4wvlOyI+x$qDn!lgL0qii&J(L>nsj0jQK zVbWPD5#(GWF+U0}c4d0Jhxw$p?O+eU6n;L~h(N2B7JOxtolmE0IOyx0Lr+AE?b(PSkDg56jKSRKAM zG4W~-gP%O37C+N@irU~b65coU=(O|O^=YV@m6X@6#@6C%scdh|CEGCM;`7LBr#7l9 z%J7!KFSDavg+>M?b^dLyPwU?bW8EAFC7VX%I~a=L?h)h6^LwOOcc3ecetE-KKk~!* z&bX8ZKCgG???can^l({_Po}W%kkgs4e*F226JQtEDjT9u zC&#U+KB7p}kU5llv|V=FQoo_!rkU_lY?qgp={TqE(PePt6f8XHF7pGS^<>XK;Z3C} z=U*&8hibNaV-x(wOqBD}?Q8v(utyRQ^qJ3e!eHrr-_Uy_kcU zEFAvRBaF)_lNhB`=^A)S>>q$T;`a`DsDCPbUoOO-8^*_AeN0t8Pcl-;z8k`N-nslw+A!1#auyqzSVjUmPKaZ`Ya_J2Q103Royp`!p=-NTRU uf>IdJryuPFB{5@W)j_af~66kE}uuvHZkP68{HSV|Gpe delta 10076 zcmZ{K1yqz>*EYx_2nZC;$O!?&#W6GD`otx>#_q*M=au`wVhWG18SB;v#>BcHAl!ow+n zE~`5~D%@rRPCG(R*Z`TiS!CVi>5=RC>c!E@>A~pEmInji=Y5X6K5PNt!|~|Q)c2*1 zbMyX^(o&ae2UkZ#{a0I^j&Q&OXaIa9rPD5!ThQ-MbadvWC4J9P+uK*8)61c$Xh#Pb zI+-ibCj$3iaG>8ses4e@Z(j}19;_ZGWfSf7uO1${pcYqWJBCY{eSP4sqNYV0=NiP# zW*X0@Lzg>3OM_n=U7QL}mI;a56lNVL<%>DfUV=m8|9E|-sXy4}R3d{8Uh z->&y&58!AOK-vp&*4@2~UM>t>-pb1x-soIYK_yL3BmD_8k^bJ^SKp7gMy~d)^h4|7 z$2&ws0H5o!*(=oW!EC7`Lq>-7+5D?IgcrgKwX%zX$Jg-$%>%V(^@hXnVQ+7j>tC5> zmpfOBsFrbHQTn>2lC4yl;i6;X1dBc+!{6T!0Nk%uw)gro(i+ak@6Wd}Ka^TRd$ty- zn6xsJOD*Xfgnd=HxcC02ckBCuuxynJ()an^t%sh9!_3^v_Y3I{vMViW@?|^-%k?hb ztlfi*LB&HSYR%jT%R{7G3Cj~3c4;j%>Oi&HZxRq&k%{zi_dri0sWK85=grd@lBNoP z6VWvL24vEJ{<*w7WLc(H>W+UbG31bewQBTAO&+qgMm&WFRkwkm<>ga%`iZA%FVCD2 z6SbEQvWlk@gKHdcaP31%SgXjUSgOdp+08eGD~~pe!ZJeE|BBipF2#QN8)5iI{wq{% z*6QO1`n(J=Vqv)U^3$MMVqMM8n|LrjEJN{#C7aD{LrqcV!G9U4%^E?2ryJu6$?G4??ayj7k5=f85~QU<*4Y(xn7r`@pVc0b zX6N1eFGE;{{E;MKyi@`w?q#)#JTfM?tN5Oj6^>N3jt;9g$KbQ@4E3WCZmGoeZ;JnG zc%Gd{EcL$SB}r|ZI&$}$B9Wd07$5yPT!J)>dt1v6l<`#>dX!zer>as98TJ1SD_SO0$n$dSOOnplc+P$6*|^h&HJ_n}XH zuTHLeI?p}`O30>&%f$iM!MTUr3gTbMF7v-Jq2e0A&l-c?j6E`;H z1$-FEJm&Os*{;VaK7POd$w9pE^LKJNc3;#}RbS2_w0F1w0aK4H|V8@9^X)&R?EUQb0KF_e|`C-G55*{I=rIT3*y5zZRN7bSzgO!Khf9mgWSN$S@1%RgWGT5 zk5x(4VIVxKW{i!kURjOW>`T6TeATH{LZ6r7h^(0oZn(PYbNOuG$qOC7pc^nTaR{-G zd}J{*&du}z1cs$sal!r1L;fM0n{|3WR*|xC>@dA%tjdd0v|{04)_qi+*El-Hgr!Rr zE?+K{H97X-&-8n7b>JhTt9|Ok5BX7wx;qbTaSc|zh;*?!rN~$c7RYI>O-<42Fw}h z`R^2kSk$!|(?|IkUyuCkO`los0k_|)Ix@qCmwpCP(V+}mirJ+KPx z#eZ)rs?~~uO0}?2L}_PfvH8J@tSlT%y5Drb zz5`qAPa%N6?bdg5eik(?FjDbDs&uh22XRi3w3L`IYSBY6X0X4yRG%wnyJUj>WCiZ| zW*cCS-%M}Hz~aUr^hF?ZFxhe^cg?Nk+eTP^pg<-*v0iNcTSc(tj>0)qpUY(rB=@iP zxUIP1L^o0g(ShLj_yiX$U$3kHH`z7Bi313X1%P(<$lpbZ=u>SMRq&eNx3oYzQeQE#lP@eqSc<8hX z^i_fNkAx1iX-oxj3Ybj%Xgrdbuzf^iEXoUk0MJJE4u7W}JbV_E&vq$Gj6^IM5WM9~ z4gHSSq6e1i-Vr|Dc@Dox;-b79@(4Y3lh+DSR(VtD#KnTU6nNg8gidPPqXziwr$Jw{ z!j*?_0ye?CZ~f)ucs5obI=(#RG6ndY8MGz-k#%}rj$_Sl=L#AmB5{ya6X;4V6!zoD z$mJYow|XZE%Ezc{KnUltdvsS@s#cL0l$t0dJ}K}EbYM_pYB5rn>3i8ylL2yZw%`_Y zd%akEKO_ihOnYjaT$*@hoXo?wvzYB~)%axi&4D!uzHhP5sXCTU=D@H9^;w+gzExYF z!6qCDNj%JPUg$1J=b^tCvRLc>&$q_X{^Qt^Dkv}x+hm|VZpGwDOOw5ku|1B?e#5tka&G$NnokY`b}oN*pg<~L}}BI4Gv8M^IE5UjYzEjC}mk@22;om z?XQQWLIkd4xrN)=Fve|w1kXoJ!u11tV;PubXIL+pdpCbvbm%YT0!;&}T8?2>7}Xq= zZ4OafJ5hASFHOWQXrEOYZ<=5eM1@ED$y;VlOyo6|>)36R(D+1OM#is6o?2PlSlNWZ zfQH_gy+Ey^Z4wS2?m1#kbd|v&uA_WYb4ct2F%?!xQ~1 zZcNiX4KjKb0D-tFwF9m?xWpgyT3|(g{08LKxo6Ffgl%0I_90ksmZ$N8}Kt=zb&}e3}n-TpVwEgPb*9Z;2 zpX=MK1WErN?0=w_hnH{YFZe+w|1L3Sy=!oiT#v@g41>{$*GpJ`M0q7u^<&v)>Q24- zh@#5zv$oAd4@p;6^#fcxLNUJORsCE$0vhO^8{0VBgL<*2B$slTI0hu%OJNicbba^K zT3~g*&5i&B9W^FhS>5jk%bc8gRRb2Re&S|KF7~s&F-`YeYGi3zT*a#OKrL&n#=yii z{a-SrlB^dM(WSO20G#v}?Fii1K-%**8ng!w1c2y`kCW4D zrgPFn^?u|$TN1bZUzN}Qc67$jTyDeMatyT8o4adYE}{UT#b_~g;`h6M(EJVsxpG_I z*QhhzhF?Hv*RzWq2#vRqy{^3n*5w`w5Xtp-eOPZYIS=)cOS^paRtr&Q1w96URU^k? zNZsN~;uRjd^(G3`6vTN!!{VAfmowcEaXMsiJvF%GbevPZXJmU2RvK8>AAR<^$H20Y z`f&8zS7^164;p>eOyJD#M+H(K^j-89)u%m;n*hdEI4Wz$y>J@o^R8$tQ`5!jrD{w) z33~j@s7ezLjQwMd1)0zfl-QA?(_O@6lSm6ad~0G;9XsnG$9bcKvOs z?Ta0}Mcz#YRG-ifEXLp#SA7;pvfe*G%)fMU_I8WAmKYujEBaSVFnV%3UnF}!7QSO< zXm{_T$TD=kAjx80?}{SwGL#Rzs-P*g#%->ZB-78pCmZPpw;OO?c*3v+wI4NAu|`4!&~xWKqU!I{YS8lt$u1z1{U?i%ZZewjkWNCfImt| z>IjJBC<6w-lY_og}8z3x{n<`y`oBej#n*=0Y+Egmj%}G zBw`lcHMD9%>@@--(l7qiJo9J^2?O-7qtLs>;uyKko{@@v&~->CObhEKqz;4TCwhpo ztUW*oI|{O=v4u7ES1ICq5seAEO9`L5NoeVq@rY%eC=MGLQ4Jmv#+W;%@W6?cRg*rp z9hmXBvB>Z5^pCC&{I$*%V7;KO`;jDi>)kB|n23ZkhTSZdA>Apu>#DH=IW7v)!6?e` zz^s&`E9<&4TgLd*j>%i--o)kJw0$7{LVqVfO%ly>ZlXpz_dw_vuBt`Ho0^2B??H7T z&U>0j;BOVgbFht^aLyXNU$z0&sU*q$;*Dcy72*htHs*xut29?lNi=bx2zh;)0?Wd2 zhbLr|^%G80&O3gpM*x*7zvxKal3Gi#!l6b>1|Ls}J?<=HUx*XkNloUhIn}4TvrGq; z?4GcI`AlG`GmEC>H7&dIDJ{G0DXmfGZ%0)AE%Q1_e>%^LE_!4|M-oDT6235QB|Fg^ zTn_w%Z^}|SzHuB4EDOP1ruU6zHWx<#u|;2(!RuzctmB&puY>!N>;!ir6o&&@(N$XX zn|1eIyXmPA!5s}v$~xad!r%k56OaU`!53k)xWl6NJ_ck3)vj5pS%z+ zJ4>|pw@T%44&|*m6C%3Ptx35BVu9fDK9GT-b+h4`!Pprw8m?v#_bwH|N&ux$qlPx~ zW7Tx8jtBOZ+8uIE`-V#dl=qfeVHlb>8Lok+9p%=jhXZzAyEw9<4?baK$jS=D3L@}| zX8)w!xDX#WOLX(MO20`vc2-Szr%;m;eM7mt|0b(dGcd7cfI1`g>HnSh^?)@BRfGH} z*chybZI{$nij@vE`a{YC!!jK@jRCga4GiPY3O8No=WmtiSiUtDA&t!L@0jYP*JHFx z>MzC0giiih)#Hw3CUhLroRg$oe!c+9oa+cYr5b*4zFO@`|bH}|HI!n(WhqU-u< zK@N$6w{j2h0{)9Ubcyf<*AWECH;+y5-z9r|_iyPsGsG7k|8}O9R#gW;Av|&>T--Fy zoj#8(jx6tUEIMS}-ER{|9!nhB4>VCsgXiU2-V5JsGO%czScSMPY&p{^hJ5|NW7nG` z_&f1asVP5_;HAH1^1R`y#>`@*mH*aWfn^YO3*I1?_#04BJT=dt;bwHKQKwt5&*dW~)V zdnttUr_BV1(21VA z$Dh}sXFn7p%VW|j#Lc&V3lU0Wt}g`gN>$I(*_O}KO&b1!RXZ-7L~>i5-ciH%8xS;k ztG4N(GF9aZ<&DIl$*c^pmaKjQM}m(JPI3J;u!2+rVt zJw3hLyR|&KhJ2jI=IG!ETDm7I*Uda*< z-Q~sPjZ4m(G3QyxlL~tu7jGx@E}PxLMx`23?VU}19lcwMbm%!`kQpPJ;mi3>8tTJ@ ziwX#r|7g-fbmHIWEk?N;DEvI3D()@a?&IV+^g-Kv2F!lJ5-`Q}L%?}+y;jk+hGi+Q zAJ9CM9vjsxNTJsHnu4y?jr!>QfIAxEiz4E-j>wER`?E5VxAksCxDS z#X@i~94|QG;5koeb~tZ@OQ2U}{f_8)|AY$xjddwtwFtMBH;etXzy zqgM=uo~&_2^hzhhIY#cCJ-&Y8?7U@C7^i&2y;0*^^f4OvQ2QYwslM2p3C4{;nb&1y zZI%?nLqmT!L0A&&rk{D~$eox~CbT`R^QgChO}z?w%4?3OlmC=6oedEU4i5g($ruKo zMtmx7cG@sD_2CXVSP%+?w+0=&6(%U;AZkxv7~<7)=7Yb*U;2TbR;nA(^mNjGp(%-o zP4U`)*Te@LEJzrO$TZP-Piy~bX#7QQLc&znC|QT}v8hrkX4F`4r}*q=f7Evhbf`r(6O3-*53({jNwmKH!*ZF1)2`o2273 z&wWX7N8~(tZ+?EB4tD<5KnY*f#Uz4vLfZwZkMQB@{?{(->;&j$bS#;zz z+Xo7?e=j3q2x1gH+lOCIme4;uFV-iOYSBKr{j&-xsgGLtao)mvrQ&+=etn|aW9@T{ zy%_!dx}XOtR1TR%ublJZ3r{}hC6VpMiZ30~JvoSYrstxPCW^i=n!)qpX*%C0gfpCo z`K87BRsm;8Cugpz%Im{1cs!pTML(}h#-Uh`ShupR_+7yR*XHDnz%zXl9Y`&eGyRC- zoo&mY5Bn<3nlE zo4z82TxO!|^j%eaAdv=;-V<&si0IC z+4lH&;cjZA(zGCbTCBk=K-!s%!_!&4Q%@Bk)%3}$=XINNBkXs3cDs|O?2DqWV_36g zZv}Q_+p+Gg=Q8yN;2Twu2^F@pfBA$x(aAuH37pJu>3bevh|Lw`;Z)mw5VE_A@$EDX z|CWl@GiUv&nsCIjdIVWi`hk2Wds1Q=5AIjo?a;+Q9NmmH;Oe0C5kv0FB!#4kLo2}` zUh4CNQ5G##N@7~6GD}^e>1Dc9U@7<#buq6c?h-AbevLMw!%_7y-3zZ2*WHRK9P5o# z2Ohq3)Y!I*jJ&<0=~#DLfXu1v;Pp(>tq_(dh7hy0+yNB<;+1&0g>7@LX$$if?zlNV zp`M_^7>3NB?q<67G=oe~KGxfnc4xgQjF8gq3anf4R2Ak&_fplD-mCPn(yO?%xweoA zJ$<|w!)YwnDfETf_842Yrmtn!GhB?i=PPilYxFtC4|*}-vx|*Sd>n(x3z-J3+8PELYYOLPLBN3$9Bf3nl6dD_|@NcPue3$&2fk#UUEa@yqydV zyzr4<@C!X;GM>qmDn69|#J!k7T~k4f`sMU-);_j2>FLLSwC^_# zdWOv6;OY(z)=xk7dasN3XFU_dh-)?9#BnAkcx!XoFK03M@>7TmdBw|6fysQ#r_YiY z#XJsWB*c7wE$jo*CnB@3jiCNaUIO*sohK!U&}VwaRZ1nUJqZKvoXz*H5 zW942*5<*D_k#OaYniQQ4ZghFz_*w^LF1?#afow{Rk`mcmQ|Ek&7=%8m~9 zbQvnoLM&*%J{hL@W-?|@+V+K#{SyR;@$pz>!MAsg3=zop`ho4*$9VU}Vw46JG3~pS z;Fm2E3eVUU&1K(?{GuHZ{zvztXW1DV_gO*|^(*{W&o^$j#m(RLQYy0dU*u9XYStwu zM%cOipl->lH`d#wPrs*i;*cZj)A9cEn;KRm-Z~@8y^B6Z>7uXTpM_m&M8F{$A*#Q4pefM)ipnK(|o-4N!6m% z;CR;QP65638G^M+vwWMDJgL~3}dxx{OQLx z-D5m;o`YGJAHM&Rx4b8+@UfBGAQQlgNZb^PNxdEPt$z6}<}Fm(vI1G7SGlzuXSj}9 z_pYoe2Pp~bGfzl#+-**X)wdBIN;<^n*SHi#u?32>?+-8UXn08Az0{H@Hd0gM`=$6? z=)SMq>hC^=RW)Hv<;K+6yTd97>9o6PFg&`id(QQb2wh}i-VF{EKldG;_TG5}FfI-} zmi}Hh{oUkxcV?5+IqB@2)PbLm#i;BIPU3XV=6Y+aszoj@*LS91I;3d8g!k$>*&F`g zfys0NwYRTDs4tJsrxqQ!*=h354j?i?(wSY#4r=$+Cc|FF@0kkK|0>4m_$T&>nyoEK zMKRRYGI<<^t{8br5HfuPV*+wT4n2Y``RrC~;reGYJ(h ze3+x|oqFm%<)ipz3~BEqQ4CZ3m~(RWjQ)k*)S#dY>!l>S#JU?21#1k;Qv;-~KWXAE5j#J^;wpQ)KBD z`R1)Sm8g1BJJZR?357IWvI7>Yx6VzEv~tr8%Ci^AL%zQshfY4iRM0#6(NHobIU+ zqtMQc1S+$DZ0&Vz7o2Q~Y20P4X*Jc1lgU1ks&Rpg>MDsd9BB+V#F=HZz^l1+O(n(fFr)G!*aaro)Nuq?PR+ z;$W}gizl@WEhtTBP#Z*XM36zmg7!ELv928){MQX=gwZv0g+!{F^EDQPJXU(0J}YS* zPkBYZsCW7xno}H~>6&vT?wEy%MB>}f6|Qjgv|KJ?lEa(rVp#y)iueS7FU2qE8gu&& z8KHR)nH~jTeOc2#ZPjbDlt^Ym@**I@rA+!lw(rGtgWuh3vKb9GeDHZgm>rFzct2sn zjT`jlWa&}hXF-JIRgUvSF!<`$XTFS~xsI4`!u`s4`gu1&q3n3Nw6n|<+& zo>XNvglPQ0x8V?XoD>k23|dkMXJcQ=-aJzA2o4sMi!j(F&TI|S z-j}|HOGYWXJmnDMXweZ2l2yaQ=0XHlajAQolPE@d zU@82Hs1o17no;U+`>?C~{>8YPk_gp6M0uQliP_vcy*E3AI;??J)o~`u_c}Xja>sF~ z5cfann{7m6A0r%L-fO#6KlyMtwsa4d5o=u?D;r)7EcR^$MrIZDDW%=)nUvaILf7I| z^_|dQ)X2(wd*y=Y$qt*%HAdxtxZkb+>;KyfrSxub(8iZaG130-hfoZRn@^`VUzmSn z7G3 "DuPont et al. (1995)" +short_citation <- function(reference) { + author <- sub(",.*", "", trimws(reference)) + year <- sub(".*?\\((\\d{4})\\).*", "\\1", reference, perl = TRUE) + paste0(author, " et al. (", year, ")") +} + +isolates <- data.frame( + study_id = raw$experiment_id, + isolate = trimws(raw[["Agent Strain"]]), + host = trimws(raw[["Host type"]]), + route = trimws(raw[["Exposure route"]]), + dose_unit = trimws(raw[["Dose Units"]]), + endpoint = trimws(raw[["Response type"]]), + dose = as.numeric(raw$Dose), + positive = as.integer(raw$Positive), + negative = as.integer(raw$Negative), + source = short_citation(raw$Reference) +) +stopifnot(nrow(isolates) == 32L, !anyNA(isolates)) +utils::write.csv(isolates, "inst/extdata/cryptosporidium-isolates.csv", row.names = FALSE) +cat(sprintf("Wrote %d rows from %d experiments.\n", nrow(isolates), length(unique(isolates$study_id)))) diff --git a/data-raw/pooling-article-results.R b/data-raw/pooling-article-results.R new file mode 100644 index 0000000..55adfcf --- /dev/null +++ b/data-raw/pooling-article-results.R @@ -0,0 +1,40 @@ +# Precompute the slow results in the pooling article +# (vignettes/articles/pooling.Rmd): the two group_datasets() searches take +# about a minute together, so the article reads them from +# vignettes/articles/pooling-results.rds instead of rerunning them on every +# render (see #39). Rerun this after changing the pooling code or the example +# data. The article stops with a message if the data no longer match, or if a +# live test of one pair disagrees with the saved pairwise p-values. +# +# Run from the package root: +# Rscript data-raw/pooling-article-results.R +pkgload::load_all(".", quiet = TRUE) + +path <- "inst/extdata/cryptosporidium-isolates.csv" +raw <- readr::read_csv(path, show_col_types = FALSE, col_types = readr::cols(study_id = "c")) +trials <- split(raw, raw$study_id) + +warnings <- character() +explore <- withCallingHandlers( + group_datasets(trials, method = "agglomerative"), + warning = function(w) { + warnings <<- c(warnings, conditionMessage(w)) + invokeRestart("muffleWarning") + } +) +exhaustive <- suppressWarnings(group_datasets(trials, method = "exhaustive")) + +saveRDS( + list( + explore = explore, + exhaustive = exhaustive, + warnings = unique(warnings), + # The data themselves, not a file checksum, so line-ending conversion on + # checkout does not look like a change. + data = as.data.frame(raw), + singlehit_version = as.character(utils::packageVersion("singlehit")), + created = format(Sys.Date()) + ), + "vignettes/articles/pooling-results.rds" +) +cat("Wrote vignettes/articles/pooling-results.rds\n") diff --git a/inst/extdata/cryptosporidium-isolates.csv b/inst/extdata/cryptosporidium-isolates.csv new file mode 100644 index 0000000..62584e0 --- /dev/null +++ b/inst/extdata/cryptosporidium-isolates.csv @@ -0,0 +1,33 @@ +"study_id","isolate","host","route","dose_unit","endpoint","dose","positive","negative","source" +"108","Iowa strain","human","oral","oocysts","infection",30,1,4,"DuPont et al. (1995)" +"108","Iowa strain","human","oral","oocysts","infection",100,3,5,"DuPont et al. (1995)" +"108","Iowa strain","human","oral","oocysts","infection",300,2,1,"DuPont et al. (1995)" +"108","Iowa strain","human","oral","oocysts","infection",500,5,1,"DuPont et al. (1995)" +"108","Iowa strain","human","oral","oocysts","infection",1000,2,0,"DuPont et al. (1995)" +"108","Iowa strain","human","oral","oocysts","infection",10000,3,0,"DuPont et al. (1995)" +"108","Iowa strain","human","oral","oocysts","infection",1e+05,1,0,"DuPont et al. (1995)" +"108","Iowa strain","human","oral","oocysts","infection",1e+06,1,0,"DuPont et al. (1995)" +"139","Iowa isolate","human","oral","oocysts","infection",30,2,3,"Messner et al. (2001)" +"139","Iowa isolate","human","oral","oocysts","infection",100,4,4,"Messner et al. (2001)" +"139","Iowa isolate","human","oral","oocysts","infection",300,2,1,"Messner et al. (2001)" +"139","Iowa isolate","human","oral","oocysts","infection",500,5,1,"Messner et al. (2001)" +"139","Iowa isolate","human","oral","oocysts","infection",1000,2,0,"Messner et al. (2001)" +"139","Iowa isolate","human","oral","oocysts","infection",10000,3,0,"Messner et al. (2001)" +"139","Iowa isolate","human","oral","oocysts","infection",1e+05,1,0,"Messner et al. (2001)" +"139","Iowa isolate","human","oral","oocysts","infection",1e+06,1,0,"Messner et al. (2001)" +"140","TAMU isolate","human","oral","oocysts","infection",10,2,1,"Messner et al. (2001)" +"140","TAMU isolate","human","oral","oocysts","infection",30,2,1,"Messner et al. (2001)" +"140","TAMU isolate","human","oral","oocysts","infection",100,3,0,"Messner et al. (2001)" +"140","TAMU isolate","human","oral","oocysts","infection",500,5,0,"Messner et al. (2001)" +"141","UCP isolate","human","oral","oocysts","infection",500,3,2,"Messner et al. (2001)" +"141","UCP isolate","human","oral","oocysts","infection",1000,2,1,"Messner et al. (2001)" +"141","UCP isolate","human","oral","oocysts","infection",5000,2,3,"Messner et al. (2001)" +"141","UCP isolate","human","oral","oocysts","infection",10000,4,0,"Messner et al. (2001)" +"181","TU502","human","oral","oocysts","diarrhea",10,2,3,"Chappell et al. (2006)" +"181","TU502","human","oral","oocysts","diarrhea",30,3,2,"Chappell et al. (2006)" +"181","TU502","human","oral","oocysts","diarrhea",100,5,2,"Chappell et al. (2006)" +"181","TU502","human","oral","oocysts","diarrhea",500,3,1,"Chappell et al. (2006)" +"183","Moredun isolate","human","oral","oocysts","infection",100,2,2,"Okhuysen et al. (2002)" +"183","Moredun isolate","human","oral","oocysts","infection",300,2,3,"Okhuysen et al. (2002)" +"183","Moredun isolate","human","oral","oocysts","infection",1000,1,2,"Okhuysen et al. (2002)" +"183","Moredun isolate","human","oral","oocysts","infection",3000,3,1,"Okhuysen et al. (2002)" diff --git a/vignettes/articles/pooling-results.rds b/vignettes/articles/pooling-results.rds new file mode 100644 index 0000000000000000000000000000000000000000..a0611ade890994ed171eb83dfb65d5649eddc28c GIT binary patch literal 1276 zcmVk$<9so*4??&xwCmxOG;X3)E2D@{uCdmKMXdu5~@EchDv?>5K$`qA(~=Z6*Qt2 z@sX9CopUoYP0VWzN!f+rp67SYz31+nx%Yumgpe2c0!7H{M+hN*NBf-(QG`5I2%!K% z-YU}op-O~MIr5>BEFTu5oNR33_K**rqsx}3r;oS-HX zy^k<#ZAuS9o(hCejWmVKpLAY3{ru5gFG|5nLmO}KJQaIVIdJ5qBOgek74^!wM{bl_ z-x~eAeeX241HUcnGeb3>ejjVBZQChso|HbUiY8ut_^5Q<$ge|xANmvUzz+PHOU@qK zC%#iB?fLkJIl|tj$VbNpk%rKcJ^^DwL(%*T5GAz>IL%lugwqA^krM% zJN4xx{^A*FkLUN3r&_O;#x6X*^Q)WdrP0dnd%xaSBb|KV%rzsIE^s^WdyN$DUU%(? zbh_C$v1MXFI`?g~`dM;7IyE#In|{L2@xTuJj>rBPKmX7%>60_w=;0H;NY8!wL)pdP z1dj`NV8^%WiRk^G{Ya#r!r|eM1`hK4ATHp69r&HyH1YG@sd_#S$PeNI9@sgyfIl;P z-2<7`U*~L|MHvg#v(T5a>6(2C?Q6}nFKAmuGoz#Ui|tsR%jqqSPx567J!mOC zXf3~h3XP`FXco|DOgH$Bx~)Bdh1;R3suE~sUSoqNrnNT+1Az^p#T zStlRW`Rl(xuO#iPLoucC0cNlQd!@0?I+qZ&8^;a8kCO$#st^PfBI|kr%S+^!mjgoC zlbFyF%bhn`9<)aN{oFX%k{h44#_{3Y_{VN@+#Wi(-{hC;Ut{^Xalqe=b9~r|H=n;d zFTiTsO}X@y%k+2IaesGSbH!VrpL@Q1KU6yX-FeRU zs%QT z#c)BnW|h@D60)Y@WYVtPwgTGw9!;bN7*_4t$p5Kz-Y}Iyk7Pz#8R#Smd=fcNkiGTo60G&n|F-R}xRoo-n_yD76J&xzNH~+?UZw|Ly zj_dKu?!f$Fs4*nET_R#+TSHiEh(ulM0rz;txL-?B9WpD${RuhE?78>Jx=Pd@#^)&; m@PV Pairwise likelihood-ratio tests (exponential model, +> `r choose(nrow(pe), 2)` unadjusted comparisons) showed that the TAMU +> isolate differed from every other experiment (p < 0.01), and a test of the +> three isolates studied by Messner et al. (2001) rejected a common curve +> (p < 0.001). We therefore modelled isolates separately. Beta-Poisson +> comparisons had low power (two degrees of freedom, four to eight doses per +> study) and some fits did not converge, so we did not use them to group +> experiments. Experiments 108 and 139 appear to share subjects and were not +> pooled. + +## In the Shiny app + +The app's Pooling tab runs the same functions. It requires a written reason +before each test, keeps a session log of every test and exploration, and +marks any set chosen from an exploration, or lying within one, as +*data-informed* in the tab, the log, and the label of pooled data sent to +Fit. Its Library tab can filter experiments by strain or isolate, so a +grouping can be decided before anything is tested. + +```{r premises, include = FALSE} +# The text above interprets these results; stop the build if they change. +stopifnot( + p_iso_exp < 0.001, + p_iso_bp < 0.05, + all(pe["140", others_140] < 0.01), + pe["108", "139"] > 0.05, + pe["139", "181"] > 0.05, + mean_off(pb) > mean_off(pe), + !identical(groups(explore, "beta_poisson"), groups(exhaustive, "beta_poisson")), + !together(exhaustive, "beta_poisson", "108", "139"), + # 181 (diarrhea) shares a group with an infection experiment in every result. + all(vapply(list(explore, exhaustive), function(result) { + all(vapply(c("exponential", "beta_poisson"), function(model) { + any(vapply(c("108", "139", "140", "141", "183"), function(id) together(result, model, id, "181"), logical(1))) + }, logical(1))) + }, logical(1))), + same_rows >= 6, + pb["108", "139"] > 0.05, + length(warnings) > 0, + max(n50) / min(n50) > 3 +) +``` diff --git a/vignettes/getting-started.Rmd b/vignettes/getting-started.Rmd index ee76d53..9648255 100644 --- a/vignettes/getting-started.Rmd +++ b/vignettes/getting-started.Rmd @@ -313,7 +313,19 @@ pair[c("model", "p_value", "poolable", "converged", "conclusion")] For this example, A and B are candidates for a shared fit; keep C separate. The report and grouping functions test poolability and return assignments. -They do not create a pooled analysis or establish goodness of fit. A later -pooled fit must preserve separate trial observations at repeated doses; -simply coercing stacked trials would aggregate those observations or reject -mixed identities. +They do not create a pooled analysis or establish goodness of fit. +`fit_pooled_dose_response()` fits the shared curve, stacking the trials so +that observations at repeated doses stay separate: + +```{r pooling-fit} +pooled <- fit_pooled_dose_response(collection$datasets[c("A", "B")], models = "exponential") +pooled$exponential +``` + +These synthetic trials were built so that the answer is clear. With real +studies it rarely is: each study has few doses, the test has little power, +and choosing groups from the test results biases the pooled fit. The +[pooling article](https://seanthimons.github.io/singlehit/articles/pooling.html) +works through real Cryptosporidium challenge data: how to read the pairwise +results, why "not rejected" is not "the same", why groupings are not unique, +and how to report a pooled fit. From 8a6a2a5b51a9dc8ec29a3f6a2b8531cbf694b933 Mon Sep 17 00:00:00 2001 From: seanthimons Date: Wed, 28 Oct 2026 18:43:12 -0400 Subject: [PATCH 6/6] feat(pooling): expose pooled sub-study observations - Carry per-experiment dose groups into Plots and Explore - Colour pooled observations by sub-study metadata and inspect fitted values --- app/R/app_server.R | 6 +- app/R/fct_explore.R | 26 +++++++- app/R/fct_figures.R | 94 ++++++++++++++++++++++++++- app/R/fct_pooling.R | 33 ++++++++++ app/R/mod_explore.R | 20 +++++- app/R/mod_figures.R | 27 +++++++- app/R/mod_pooling.R | 21 +++++- app/R/mod_upload.R | 15 +++-- app/README.md | 7 ++ app/tests/testthat/helper-data.R | 24 +++++++ app/tests/testthat/test-app_server.R | 1 + app/tests/testthat/test-fct_explore.R | 12 ++++ app/tests/testthat/test-fct_figures.R | 31 +++++++++ app/tests/testthat/test-fct_pooling.R | 9 +++ app/tests/testthat/test-mod_explore.R | 16 +++++ app/tests/testthat/test-mod_figures.R | 18 +++++ app/tests/testthat/test-mod_upload.R | 13 ++++ 17 files changed, 356 insertions(+), 17 deletions(-) diff --git a/app/R/app_server.R b/app/R/app_server.R index cfa8117..2da52c1 100644 --- a/app/R/app_server.R +++ b/app/R/app_server.R @@ -32,6 +32,8 @@ app_server <- function(input, output, session) { # 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) } diff --git a/app/R/fct_explore.R b/app/R/fct_explore.R index f1d0feb..edb0c21 100644 --- a/app/R/fct_explore.R +++ b/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 diff --git a/app/R/fct_figures.R b/app/R/fct_figures.R index 61e51d5..6cc687f 100644 --- a/app/R/fct_figures.R +++ b/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 diff --git a/app/R/fct_pooling.R b/app/R/fct_pooling.R index 45b613e..b7daf55 100644 --- a/app/R/fct_pooling.R +++ b/app/R/fct_pooling.R @@ -632,6 +632,39 @@ pool_fit_data <- function(datasets, selection, data_informed = FALSE) { counts } +#' Each experiment's dose groups with its metadata, sent alongside pooled data +#' +#' The pooled dataset that Fit analyses has one row per dose, so it cannot say +#' which experiment a count came from. This table keeps that, for colouring +#' and exploring the observed points on the Plots and Explore tabs. +#' @return One row per experiment and dose: `experiment` (short label), +#' `experiment_id`, `pathogen`, `strain`, `host`, `route`, `outcome`, +#' `dose_unit`, `dose`, `positive`, `negative`, `total`, and `response`. +#' @noRd +pool_components <- function(datasets, selection) { + out <- do.call(rbind, lapply(seq_len(nrow(selection)), function(i) { + data <- singlehit::as_dose_response(datasets[[i]][c("dose", "positive", "negative")]) + data.frame( + experiment = selection$short[[i]], + experiment_id = selection$experiment_id[[i]], + pathogen = selection$pathogen[[i]], + strain = selection$strain[[i]], + host = selection$host[[i]], + route = selection$route[[i]], + outcome = selection$response_type[[i]], + dose_unit = selection$dose_unit[[1L]], + dose = data$dose, + positive = data$positive, + negative = data$negative, + total = data$total, + response = data$response, + stringsAsFactors = FALSE + ) + })) + rownames(out) <- NULL + out +} + #' Whether a set of keys lies within a set that was explored #' #' A set chosen after seeing an exploration of the same experiments is diff --git a/app/R/mod_explore.R b/app/R/mod_explore.R index 7de397d..f8a58eb 100644 --- a/app/R/mod_explore.R +++ b/app/R/mod_explore.R @@ -78,14 +78,30 @@ mod_explore_ui <- function(id) { #' #' @param id Module id. #' @param analysis Reactive returning a `qdr_analysis`, or `NULL`. +#' @param components Reactive returning the per-experiment dose groups of +#' pooled data, or `NULL`. Used only when they add up to the analysed data. #' @return A reactive returning the data frame being explored. #' @noRd -mod_explore_server <- function(id, analysis) { +mod_explore_server <- function(id, analysis, components = reactive(NULL)) { moduleServer(id, function(input, output, session) { output$has_analysis <- renderText(if (is.null(analysis())) "no" else "yes") outputOptions(output, "has_analysis", suspendWhenHidden = FALSE) - all_datasets <- reactive(explore_datasets(req(analysis()))) + matched_components <- reactive({ + if (!is.null(analysis()) && components_match(components(), analysis()$data)) components() + }) + all_datasets <- reactive(explore_datasets(req(analysis()), matched_components())) + # Offer the sub-study observations only while they exist, and leave them + # for the observed groups if they disappear while selected. + # A reactiveVal ignores a set to the same value, so re-running the analysis + # does not touch the dataset input (which would upset a pending preset). + has_substudies <- reactiveVal(FALSE) + observe(has_substudies(!is.null(matched_components()))) + observeEvent(has_substudies(), { + selected <- isolate(input$dataset) + if (!has_substudies() && identical(selected, "substudies")) selected <- "observed" + updateSelectInput(session, "dataset", choices = explore_choices(has_substudies()), selected = selected) + }, ignoreInit = TRUE) # Judged on the data, not on whether a bootstrap is attached: a bootstrap # with too few converged replicates gives bare curves with no bands to draw. output$has_bands <- renderText({ diff --git a/app/R/mod_figures.R b/app/R/mod_figures.R index 0754c01..e36a7d1 100644 --- a/app/R/mod_figures.R +++ b/app/R/mod_figures.R @@ -8,6 +8,7 @@ mod_figures_ui <- function(id) { sidebar = bslib::sidebar( radioButtons(ns("type"), "View", choices = figure_types), uiOutput(ns("band_hint")), + uiOutput(ns("colour_by_ui")), conditionalPanel( "input.type == 'single'", ns = ns, @@ -29,9 +30,25 @@ mod_figures_ui <- function(id) { #' #' @param id Module id. #' @param analysis Reactive returning a `qdr_analysis`, or `NULL`. +#' @param components Reactive returning the per-experiment dose groups of +#' pooled data (see [pool_components()]), or `NULL`. #' @noRd -mod_figures_server <- function(id, analysis) { +mod_figures_server <- function(id, analysis, components = reactive(NULL)) { moduleServer(id, function(input, output, session) { + # Only components that add up to the analysed data are used, so a stale + # table can never colour another dataset's points. + matched <- reactive(!is.null(analysis()) && components_match(components(), analysis()$data)) + output$colour_by_ui <- renderUI({ + req(matched()) + keep <- isolate(input$colour_by) + selectInput( + session$ns("colour_by"), + "Colour observed points by", + choices = c("None (pooled dose groups)" = "", substudy_fields), + selected = if (is.null(keep)) "" else keep + ) + }) + # Refresh the model list whenever the analysis changes so a model from an # earlier run cannot stay selected. observeEvent( @@ -55,7 +72,13 @@ mod_figures_server <- function(id, analysis) { if (input$type == "single") { req(input$model %in% names(analysis()$fits)) } - build_figure(analysis(), input$type, input$model) + build_figure( + analysis(), + input$type, + input$model, + components = if (matched()) components(), + colour_by = if (matched() && isTRUE(nzchar(input$colour_by))) input$colour_by + ) }) output$empty <- renderUI({ diff --git a/app/R/mod_pooling.R b/app/R/mod_pooling.R index e76ca86..55f250a 100644 --- a/app/R/mod_pooling.R +++ b/app/R/mod_pooling.R @@ -52,7 +52,9 @@ mod_pooling_ui <- function(id) { #' @param workbook_loaded Reactive counter that changes on every workbook load, #' including reloading the same file, which leaves `rows()` unchanged. #' @return A list with `sent`: `NULL` until pooled data are sent to the Data -#' tab, then a list of `data`, `label`, `source`, and `sent` (a counter). +#' tab, then a list of `data`, `label`, `source`, `components` (each +#' experiment's dose groups with its metadata, from [pool_components()]), +#' and `sent` (a counter). #' @noRd mod_pooling_server <- function(id, rows, experiments, pool_set, workbook_name = reactive(NULL), workbook_loaded = reactive(0L)) { moduleServer(id, function(input, output, session) { @@ -308,7 +310,14 @@ mod_pooling_server <- function(id, rows, experiments, pool_set, workbook_name = req(current, !inherits(current$run, "error"), identical(current$selection$experiment_key, selection()$experiment_key)) chosen <- current$selection informed <- set_origin() != "Library" - data <- tryCatch(pool_fit_data(pool_datasets(rows(), chosen), chosen, informed), error = function(e) e) + built <- tryCatch( + { + datasets <- pool_datasets(rows(), chosen) + list(data = pool_fit_data(datasets, chosen, informed), components = pool_components(datasets, chosen)) + }, + error = function(e) e + ) + data <- if (inherits(built, "error")) built else built$data if (inherits(data, "error")) { showNotification(conditionMessage(data), type = "error") return() @@ -318,7 +327,13 @@ mod_pooling_server <- function(id, rows, experiments, pool_set, workbook_name = if (informed) { label <- paste0(label, " [data-informed set]") } - sent(list(data = data, label = label, source = "Pooling tab", sent = send_count())) + sent(list( + data = data, + label = label, + source = "Pooling tab", + components = built$components, + sent = send_count() + )) append_log(pool_log_rows( kind = "sent to Fit", selection = chosen, diff --git a/app/R/mod_upload.R b/app/R/mod_upload.R index f3b312c..eb20d95 100644 --- a/app/R/mod_upload.R +++ b/app/R/mod_upload.R @@ -68,11 +68,13 @@ mod_upload_table_ui <- function(id) { #' #' @param id Module id. #' @param incoming Reactive returning `NULL` or a list with `data` (a data frame -#' of dose groups), `label` (where it came from), and `sent` (a counter, so -#' sending the same experiment again reloads it). +#' of dose groups), `label` (where it came from), optionally `source` and +#' `components` (per-experiment dose groups for pooled data), and `sent` (a +#' counter, so sending the same experiment again reloads it). #' @return A list of reactives: `data`, the standardized data or `NULL`; #' `workbook`, the workbook rows or `NULL`; `workbook_name`, where it came -#' from; and `workbook_loaded`, a counter +#' from; `components`, the per-experiment dose groups that came with pooled +#' data from the Pooling tab, or `NULL`; and `workbook_loaded`, a counter #' that changes each time a workbook loads successfully. #' @noRd mod_upload_server <- function(id, incoming = reactive(NULL)) { @@ -136,7 +138,11 @@ mod_upload_server <- function(id, incoming = reactive(NULL)) { observeEvent(incoming(), { sent <- incoming() from <- if (is.null(sent$source)) "Library tab" else sent$source - load_source(list(data = sent$data, label = paste0(sent$label, " (from the ", from, ")"))) + load_source(list( + data = sent$data, + label = paste0(sent$label, " (from the ", from, ")"), + components = sent$components + )) }) # Same reader settings as singlehit::read_dose_response(), so the offered @@ -216,6 +222,7 @@ mod_upload_server <- function(id, incoming = reactive(NULL)) { data = reactive(if (inherits(parsed(), "error")) NULL else parsed()), workbook = reactive(workbook()), workbook_name = reactive(workbook_name()), + components = reactive(loaded()$components), workbook_loaded = reactive(workbook_loaded()) ) }) diff --git a/app/README.md b/app/README.md index 5f9247e..4cec21f 100644 --- a/app/README.md +++ b/app/README.md @@ -76,6 +76,13 @@ depends on the app. singlehit adds up counts at doses the experiments share, so the fitted parameters match the pooled fit but goodness of fit and the bootstrap use the merged dose groups, while the test keeps each experiment's rows apart. + Each experiment's dose groups and metadata (experiment, strain or isolate, + host, route, outcome) travel with the pooled data: the Plots tab can colour + the observed points by any of them, one point per experiment and dose sized + by subjects, and the Explore tab adds a "Sub-study observations" dataset + with the pooled model's fitted value at each dose (one row per experiment, + dose, and model). Both are used only while they add up to the data being + analysed. An **exploratory** section runs singlehit's `group_datasets()` on the set (every pairwise test as a heatmap, then a grouping over all ways to split diff --git a/app/tests/testthat/helper-data.R b/app/tests/testthat/helper-data.R index c62e540..4d8287e 100644 --- a/app/tests/testthat/helper-data.R +++ b/app/tests/testthat/helper-data.R @@ -26,3 +26,27 @@ example_analysis_with_bands <- function(models = c("exponential", "beta_poisson" qmra_sample_path <- function() { system.file("extdata", "qmra-wiki-sample.xlsx", package = "singlehitapp") } + +# Three Cryptosporidium isolates from the sample workbook, pooled as the +# Pooling tab sends them: the pooled data, the per-experiment components, and +# an analysis of the pooled data. +pooled_crypto <- function(bootstrap = FALSE) { + rows <- read_qmra_workbook(qmra_sample_path()) + experiments <- library_experiments(rows) + keys <- experiments$experiment_key[experiments$experiment_id %in% c("140", "141", "183")] + selection <- pool_selection(rows, experiments, keys) + datasets <- pool_datasets(rows, selection) + data <- singlehit::as_dose_response(pool_fit_data(datasets, selection)) + analysis <- suppressWarnings(singlehit::analyze_dose_response( + data, + models = c("exponential", "beta_poisson"), + bootstrap_times = 0L, + exact_bootstrap_times = 0L, + check_trend = FALSE + )) + if (bootstrap) { + boots <- suppressWarnings(singlehit::bootstrap_dose_response_models(analysis$fits, times = 20, seed = 1, backend = "sequential")) + analysis <- singlehit::add_bootstrap(analysis, boots) + } + list(data = data, components = pool_components(datasets, selection), analysis = analysis) +} diff --git a/app/tests/testthat/test-app_server.R b/app/tests/testthat/test-app_server.R index e67829a..1fb6c32 100644 --- a/app/tests/testthat/test-app_server.R +++ b/app/tests/testthat/test-app_server.R @@ -32,5 +32,6 @@ test_that("pooled data sent from the Pooling tab become the current dataset", { session$setInputs(`pooling-send` = 1) expect_match(unique(data()$study_id), "^Pooled: ") expect_equal(unique(data()$dose_unit), "oocysts") + expect_setequal(unique(upload$components()$experiment), c("#108 (Iowa strain)", "#139 (Iowa isolate)")) }) }) diff --git a/app/tests/testthat/test-fct_explore.R b/app/tests/testthat/test-fct_explore.R index 9c6d21e..494b668 100644 --- a/app/tests/testthat/test-fct_explore.R +++ b/app/tests/testthat/test-fct_explore.R @@ -64,3 +64,15 @@ test_that("mapping_matches() compares esquisse's slot mapping with the wanted va expect_false(mapping_matches(list(xvar = "dose"), wanted)) expect_false(mapping_matches(NULL, wanted)) }) + +test_that("pooled data add the sub-study observations with fitted values", { + pooled <- pooled_crypto() + datasets <- explore_datasets(pooled$analysis, pooled$components) + substudies <- datasets$substudies + expect_equal(nrow(substudies), 2L * nrow(pooled$components)) + expect_true(all(c("experiment", "strain", "model", "fitted", "response") %in% names(substudies))) + expect_true(all(substudies$fitted >= 0 & substudies$fitted <= 1)) + expect_null(explore_datasets(pooled$analysis)$substudies) + expect_true("substudies" %in% explore_choices(TRUE)) + expect_false("substudies" %in% explore_choices(FALSE)) +}) diff --git a/app/tests/testthat/test-fct_figures.R b/app/tests/testthat/test-fct_figures.R index 7b5f826..63f4521 100644 --- a/app/tests/testthat/test-fct_figures.R +++ b/app/tests/testthat/test-fct_figures.R @@ -56,3 +56,34 @@ test_that("new bands move someone off the overlay, which never draws them, and n expect_null(view_for_bands("single", TRUE)) expect_null(view_for_bands(NULL, TRUE)) }) + +test_that("components must add up to the analysed data", { + pooled <- pooled_crypto() + expect_true(components_match(pooled$components, pooled$data)) + changed <- pooled$components + changed$positive[[1L]] <- changed$positive[[1L]] + 1L + expect_false(components_match(changed, pooled$data)) + expect_false(components_match(NULL, pooled$data)) + expect_false(components_match(pooled$components, example_data())) +}) + +test_that("observed points can be coloured by a sub-study field in every view", { + pooled <- pooled_crypto(bootstrap = TRUE) + for (type in c("overlay", "panels", "single")) { + plain <- build_figure(pooled$analysis, type, "exponential") + coloured <- build_figure(pooled$analysis, type, "exponential", pooled$components, "strain") + # The merged points are replaced, not added to. + expect_equal(length(coloured$layers), length(plain$layers)) + built <- ggplot2::ggplot_build(coloured) + points <- built$data[[length(coloured$layers)]] + expect_equal(nrow(points) %% nrow(pooled$components), 0L) + expect_equal(coloured$labels[[if (type == "overlay") "fill" else "colour"]], "Strain or isolate") + } +}) + +test_that("figures are unchanged without components or a field", { + pooled <- pooled_crypto() + plain <- build_figure(pooled$analysis, "panels") + expect_equal(length(build_figure(pooled$analysis, "panels", components = pooled$components)$layers), length(plain$layers)) + expect_identical(recolour_observed(plain, example_data(), pooled$components, "strain"), plain) +}) diff --git a/app/tests/testthat/test-fct_pooling.R b/app/tests/testthat/test-fct_pooling.R index 824d752..f693a01 100644 --- a/app/tests/testthat/test-fct_pooling.R +++ b/app/tests/testthat/test-fct_pooling.R @@ -200,3 +200,12 @@ test_that("pooled data for Fit sum shared doses under one set of metadata", { expect_match(unique(standardized$study_id), "^Pooled: .*#108 \\+ .*#139$") expect_equal(unique(standardized$dose_unit), "oocysts") }) + +test_that("components keep each experiment's dose groups and add up to the pooled data", { + pooled <- pooled_crypto() + components <- pooled$components + expect_setequal(unique(components$experiment), c("#140 (TAMU isolate)", "#141 (UCP isolate)", "#183 (Moredun isolate)")) + expect_true(all(c("strain", "host", "route", "outcome", "total", "response") %in% names(components))) + expect_equal(sum(components$total), sum(pooled$data$total)) + expect_true(components_match(components, pooled$data)) +}) diff --git a/app/tests/testthat/test-mod_explore.R b/app/tests/testthat/test-mod_explore.R index 9d25f86..4c77607 100644 --- a/app/tests/testthat/test-mod_explore.R +++ b/app/tests/testthat/test-mod_explore.R @@ -192,3 +192,19 @@ test_that("the bands button is hidden when the bootstrap gives no bands", { expect_equal(output$has_bands, "no") }) }) + +test_that("sub-study observations are offered only while they match", { + pooled <- pooled_crypto() + components <- shiny::reactiveVal(pooled$components) + shiny::testServer(mod_explore_server, args = list(analysis = shiny::reactive(pooled$analysis), components = components), { + session$flushReact() + expect_true(has_substudies()) + session$setInputs(dataset = "substudies") + expect_true("experiment" %in% names(dataset())) + + components(NULL) + session$flushReact() + expect_false(has_substudies()) + expect_null(all_datasets()$substudies) + }) +}) diff --git a/app/tests/testthat/test-mod_figures.R b/app/tests/testthat/test-mod_figures.R index c8e011c..abffd64 100644 --- a/app/tests/testthat/test-mod_figures.R +++ b/app/tests/testthat/test-mod_figures.R @@ -112,3 +112,21 @@ test_that("bands arriving while on the panel or single view leave the choice alo expect_equal(updates, 0L) } }) + +test_that("the colour-by choice appears only for matching sub-study data", { + pooled <- pooled_crypto() + components <- shiny::reactiveVal(pooled$components) + shiny::testServer(mod_figures_server, args = list(analysis = shiny::reactive(pooled$analysis), components = components), { + session$setInputs(type = "panels") + expect_true(matched()) + expect_match(output$colour_by_ui$html, "Colour observed points by") + session$setInputs(colour_by = "strain") + expect_equal(figure()$labels$colour, "Strain or isolate") + + # Stale components are ignored even though the input keeps its value. + components(example_data()) + session$flushReact() + expect_false(matched()) + expect_null(figure()$labels$colour) + }) +}) diff --git a/app/tests/testthat/test-mod_upload.R b/app/tests/testthat/test-mod_upload.R index 087585c..bed49c1 100644 --- a/app/tests/testthat/test-mod_upload.R +++ b/app/tests/testthat/test-mod_upload.R @@ -125,3 +125,16 @@ test_that("a bad workbook reports an error and keeps the previous workbook", { expect_null(workbook_error()) }) }) + +test_that("components sent with pooled data are kept until something else loads", { + sent <- shiny::reactiveVal(NULL) + pooled <- pooled_crypto() + shiny::testServer(mod_upload_server, args = list(incoming = sent), { + sent(list(data = pooled$data, label = "Pooled", source = "Pooling tab", components = pooled$components, sent = 1L)) + session$flushReact() + expect_equal(session$returned$components(), pooled$components) + + session$setInputs(example = 1) + expect_null(session$returned$components()) + }) +})