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/.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/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..4a8cee5 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) @@ -313,6 +317,16 @@ trials because their searches grow rapidly. These are exploratory tests without multiple-testing adjustment; passing a test does not prove equivalence. Compare only experiments with matching hosts, units, and endpoints. +Decide which experiments belong together (same isolate, outcome, host, and +route) before testing, and use `poolability_combinations()` and +`group_datasets()` to check that choice rather than to make it: with few doses +per study the test rarely rejects, different searches return different +groupings, and a group chosen from the results biases its pooled fit toward +pooling. Also check that experiments do not share subjects, since the test +assumes independent samples. The +[pooling article](https://seanthimons.github.io/singlehit/articles/pooling.html) +walks through these points with real Cryptosporidium challenge data. + Datasets are combined by **stacking** — each trial's dose groups are kept as separate binomial observations, so repeated doses across trials are preserved (matching the QMRA-wiki pooled-experiment convention), rather than summed. diff --git a/app/DESCRIPTION b/app/DESCRIPTION index fb6e4d4..412e0ae 100644 --- a/app/DESCRIPTION +++ b/app/DESCRIPTION @@ -22,16 +22,22 @@ Imports: generics, ggplot2, golem, + grDevices, mirai (>= 2.5.0), parallel, plotly, readr, + readxl, + scales, shiny, - singlehit (>= 0.3.0) + singlehit (>= 0.3.0), + stats, + tools Suggests: pkgload, testthat (>= 3.1.7), - withr + withr, + writexl Remotes: seanthimons/singlehit Config/testthat/edition: 3 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..2da52c1 100644 --- a/app/R/app_server.R +++ b/app/R/app_server.R @@ -5,12 +5,35 @@ #' @import shiny #' @noRd app_server <- function(input, output, session) { - data <- mod_upload_server("upload") + # Data loads both the workbook (explored on Library) and the dataset Fit + # analyses; Library sends an experiment back. The reactiveVal breaks that + # cycle: Data is built first and reads whatever Library last sent. + sent <- reactiveVal(NULL) + upload <- mod_upload_server("upload", incoming = sent) + library_tab <- mod_library_server("library", upload$workbook) + data <- upload$data + observeEvent(library_tab$sent(), { + sent(library_tab$sent()) + # A sent experiment replaces the current dataset; show it there first. + bslib::nav_select("main", "Data", session = session) + }) + observeEvent(upload$workbook_loaded(), bslib::nav_select("main", "Library", session = session), ignoreInit = TRUE) + pooling <- mod_pooling_server( + "pooling", upload$workbook, library_tab$experiments, library_tab$pool_set, + upload$workbook_name, upload$workbook_loaded + ) + observeEvent(library_tab$pool_set(), bslib::nav_select("main", "Pooling", session = session)) + observeEvent(pooling$sent(), { + sent(pooling$sent()) + bslib::nav_select("main", "Data", session = session) + }) analysis <- mod_fit_server("fit", data) mod_compare_server("compare", analysis) # Plots and Explore get the analysis with the bootstraps attached once a run # finishes; Compare does not use them. analysis_with_bands <- mod_bootstrap_server("bootstrap", analysis) - mod_figures_server("figures", analysis_with_bands) - mod_explore_server("explore", analysis_with_bands) + # Pooled data also carry each experiment's dose groups, for colouring and + # exploring the observed points by sub-study. + mod_figures_server("figures", analysis_with_bands, upload$components) + mod_explore_server("explore", analysis_with_bands, upload$components) } diff --git a/app/R/app_ui.R b/app/R/app_ui.R index fe6840f..f9a93c2 100644 --- a/app/R/app_ui.R +++ b/app/R/app_ui.R @@ -12,11 +12,14 @@ app_ui <- function(request) { # run starting can make the app unresponsive; this says why. shiny::useBusyIndicators(), bslib::page_navbar( + id = "main", title = "singlehit", bslib::nav_panel( "Data", bslib::layout_sidebar(sidebar = bslib::sidebar(mod_upload_ui("upload")), mod_upload_table_ui("upload")) ), + bslib::nav_panel("Library", mod_library_ui("library")), + bslib::nav_panel("Pooling", mod_pooling_ui("pooling")), bslib::nav_panel( "Fit", bslib::layout_sidebar(sidebar = bslib::sidebar(mod_fit_ui("fit")), mod_fit_results_ui("fit")) 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_library.R b/app/R/fct_library.R new file mode 100644 index 0000000..7e09e17 --- /dev/null +++ b/app/R/fct_library.R @@ -0,0 +1,401 @@ +#' 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", + "Strain or isolate" = "strain", + "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`, `strain`, `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) + } + # 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)), + 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"]]), + strain = library_strain(strain), + 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_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)) + 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, + strain = first$strain, + 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))) + + # 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))), + x = x_label, + y = "Probability of reported outcome", + colour = names(library_fields)[library_fields == colour][[1L]], + linetype = "Curve source", + 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", + labeller = ggplot2::label_wrap_gen(width = 20) + ) + } + plot +} diff --git a/app/R/fct_pooling.R b/app/R/fct_pooling.R new file mode 100644 index 0000000..b7daf55 --- /dev/null +++ b/app/R/fct_pooling.R @@ -0,0 +1,730 @@ +#' 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` (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) + } + 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 = " ") + 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) { + 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 = "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")) +) + +#' 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 short experiment label +#' +#' 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$short[[i]] + data$dose_unit <- selection$dose_unit[[1L]] + data + }) + stats::setNames(datasets, selection$short) +} + +#' 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 + ) +} + +#' 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_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)) NA_character_ else workbook, + experiments = paste(selection$label, collapse = "; "), + dose_unit = selection$dose_unit[[1L]], + set_origin = set_origin, + rationale = rationale, + 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, + 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 +} + +#' 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 +#' 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_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_library.R b/app/R/mod_library.R new file mode 100644 index 0000000..c9f8ca9 --- /dev/null +++ b/app/R/mod_library.R @@ -0,0 +1,223 @@ +#' Library module UI +#' +#' Explore the QMRA wiki workbook loaded on the Data tab: filter its +#' 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 +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")) + } + # 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, + 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("strain", "Strain or isolate"), + 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("empty")), + 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")) + ) +} + +#' 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 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) { + filters <- c("category", "pathogen", "outcome_group", "dose_unit", "host_group", "route_group", "strain", "model_label") + loaded <- reactive({ + req(rows()) + 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$empty <- renderUI({ + 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({ + 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, + Strain = shown$strain, + `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 + ) + }) + + 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, { + experiment <- chosen() + req(experiment$fittable) + count(count() + 1L) + sent(list( + data = library_experiment_data(loaded(), experiment$experiment_key), + label = library_experiment_label(experiment), + source = "Library tab", + sent = count() + )) + }) + 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..55f250a --- /dev/null +++ b/app/R/mod_pooling.R @@ -0,0 +1,466 @@ +#' 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")), + uiOutput(ns("explore_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 A list with `sent`: `NULL` until pooled data are sent to the Data +#' 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) { + 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) + # 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, 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({ + 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( + 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"), + 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, + Strain = shown$strain, + 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")) { + append_log(pool_log_entry(run, chosen, checks(), rationale(), workbook_name(), set_origin())) + } + }) + + output$log_status <- renderUI({ + entries <- log() + 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 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." + ) + }) + + 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 or explorations 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), + 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" + 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() + } + 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", + components = built$components, + 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)) + ) + }) + + 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) + + list(sent = sent) + }) +} diff --git a/app/R/mod_upload.R b/app/R/mod_upload.R index 98fcb6c..eb20d95 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", + "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." + ), 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,14 +41,17 @@ 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")), tableOutput(ns("table")) ) @@ -44,43 +59,101 @@ 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. +#' 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. -#' @return A reactive returning the standardized data, or `NULL` when nothing -#' valid has been loaded. +#' @param incoming Reactive returning `NULL` or a list with `data` (a data frame +#' 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; `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) { +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) + # 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 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) { + # A failed workbook read is stale once something else loads. + workbook_error(NULL) mapping(no_mapping) for (field in fields) { updateSelectInput(session, field, selected = "") } - source_path(path) + 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() + from <- if (is.null(sent$source)) "Library tab" else sent$source + 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 # names match what the real read sees. columns <- reactive({ - req(source_path()) + req(loaded()) + if (!is.null(loaded()$data)) { + 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()$path, n_max = 0L, show_col_types = FALSE, progress = FALSE, trim_ws = TRUE)), error = function(e) character() ) }) @@ -100,16 +173,33 @@ 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.null(loaded()$data)) { + function() singlehit::as_dose_response(loaded()$data, chosen$dose, chosen$positive, chosen$negative) + } else { + 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({ + 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() + ) + ) + } ) }) @@ -128,8 +218,12 @@ mod_upload_server <- function(id) { 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_name = reactive(workbook_name()), + components = reactive(loaded()$components), + workbook_loaded = reactive(workbook_loaded()) + ) }) } diff --git a/app/README.md b/app/README.md index 5c94388..4cec21f 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,7 +39,84 @@ 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:** explores the QMRA wiki workbook (`.xlsx`, "raw data" sheet), + 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, 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 + 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" 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 + 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. 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"). "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. + 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 + 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. + 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, 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. ## License diff --git a/app/dev/make_qmra_sample.R b/app/dev/make_qmra_sample.R new file mode 100644 index 0000000..a837b3f --- /dev/null +++ b/app/dev/make_qmra_sample.R @@ -0,0 +1,58 @@ +# 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, 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" +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", + # 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" +) +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 0000000..1d8cc1b Binary files /dev/null and b/app/inst/extdata/qmra-wiki-sample.xlsx differ diff --git a/app/tests/testthat/helper-data.R b/app/tests/testthat/helper-data.R index 9a3ca66..4d8287e 100644 --- a/app/tests/testthat/helper-data.R +++ b/app/tests/testthat/helper-data.R @@ -22,3 +22,31 @@ 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") +} + +# 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 new file mode 100644 index 0000000..1fb6c32 --- /dev/null +++ b/app/tests/testthat/test-app_server.R @@ -0,0 +1,37 @@ +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())) + }) +}) + +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) + }) +}) + +test_that("pooled data sent from the Pooling tab become the current dataset", { + shiny::testServer(app_server, { + session$setInputs(`upload-sample_workbook` = 1) + keys <- c( + "108|Cryptosporidium hominis/ Cryptosporidium parvum|infection|exponential", + "139|Cryptosporidium hominis/ Cryptosporidium parvum|infection|exponential" + ) + session$setInputs(`library-pool_set` = keys, `library-pool` = 1) + session$setInputs(`pooling-models` = "exponential", `pooling-alpha` = 0.05) + session$setInputs(`pooling-rationale` = "Both Iowa isolate.") + session$setInputs(`pooling-run` = 1) + 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_library.R b/app/tests/testthat/test-fct_library.R new file mode 100644 index 0000000..6cfdeb0 --- /dev/null +++ b/app/tests/testthat/test-fct_library.R @@ -0,0 +1,127 @@ +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)), 18L) + # 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), 18L) + 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") +}) + +test_that("strains are read, with unreported ones marked", { + expect_equal(library_strain(c(" Iowa isolate ", "NR", NA, "")), c("Iowa isolate", rep("Not reported", 3L))) + rows <- read_qmra_workbook(qmra_sample_path()) + experiments <- library_experiments(rows) + crypto <- experiments[grepl("Cryptosporidium", experiments$pathogen), ] + expect_setequal(crypto$strain, c("Iowa strain", "Iowa isolate", "TAMU isolate", "UCP isolate", "Moredun isolate")) + expect_no_error(ggplot2::ggplot_build(library_plot(library_curves(rows, points = 20L), experiments, facet = "strain"))) +}) diff --git a/app/tests/testthat/test-fct_pooling.R b/app/tests/testthat/test-fct_pooling.R new file mode 100644 index 0000000..f693a01 --- /dev/null +++ b/app/tests/testthat/test-fct_pooling.R @@ -0,0 +1,211 @@ +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 (different strains too); + # these differences only warn. + expect_equal(failed(checks), c("host", "strain")) + expect_equal(unique(checks$severity[!checks$passed]), "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$short) + expect_equal(unique(datasets[[2L]]$study_id), pool$selection$short[[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$short, "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; One strain or isolate") + 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))) +}) + +crypto_pool <- function() { + rows <- read_qmra_workbook(qmra_sample_path()) + experiments <- library_experiments(rows) + keys <- experiments$experiment_key[grepl("Cryptosporidium", experiments$pathogen)] + list(rows = rows, experiments = experiments, selection = pool_selection(rows, experiments, keys)) +} + +test_that("different isolates warn, and short labels carry the isolate", { + pool <- crypto_pool() + checks <- pool_check(pool$selection) + expect_false(pool_blocked(checks)) + expect_equal(failed(checks), "strain") + expect_true("#140 (TAMU isolate)" %in% pool$selection$short) +}) + +test_that("exploration finds the isolate groups and logs how much it searched", { + pool <- crypto_pool() + explore <- pool_explore(pool_datasets(pool$rows, pool$selection), c("exponential", "beta_poisson", "exact_beta_poisson")) + expect_equal(explore$models, c("exponential", "beta_poisson")) + expect_equal(explore$method, "exhaustive") + choices <- pool_group_choices(explore) + iowa <- choices[grepl("^Exponential", names(choices)) & grepl("#108", names(choices))] + keys <- pool_group_keys(explore, iowa, pool$selection) + expect_equal(sub("\\|.*", "", keys), c("108", "139")) + + entry <- pool_explore_log_entry(explore, pool$selection, "sample") + expect_equal(entry$kind, c("exploration", "exploration")) + expect_equal(unique(entry$n_pairwise_tests), 10) + expect_equal(unique(entry$partitions_considered), 52) + expect_match(entry$groups[[1L]], "\\{#108 \\(Iowa strain\\), #139 \\(Iowa isolate\\)\\}") +}) + +test_that("exploration refuses sets it cannot search", { + pool <- crypto_pool() + datasets <- pool_datasets(pool$rows, pool$selection) + expect_s3_class(pool_explore(datasets, "exact_beta_poisson"), "error") + expect_s3_class(pool_explore(datasets[1L], "exponential"), "error") +}) + +test_that("test and exploration log rows bind together", { + pool <- crypto_pool() + datasets <- pool_datasets(pool$rows, pool$selection) + explore <- pool_explore(datasets, "exponential") + run <- pool_run(datasets[1:2], "exponential") + rows <- rbind( + pool_explore_log_entry(explore, pool$selection, "sample"), + pool_log_entry(run, pool$selection[1:2, ], pool_check(pool$selection[1:2, ]), "Iowa isolates.", "sample") + ) + expect_equal(rows$kind, c("exploration", "test")) +}) + +test_that("Bell numbers count the ways to split a set", { + expect_equal(vapply(1:6, pool_bell, numeric(1L)), c(1, 2, 5, 15, 52, 203)) +}) + +test_that("a set within an explored set is data-informed", { + explored <- list(c("a", "b", "c")) + expect_true(pool_data_informed(c("a", "b"), explored)) + expect_false(pool_data_informed(c("a", "d"), explored)) + expect_false(pool_data_informed(c("a", "b"), list())) +}) + +test_that("pooled data for Fit sum shared doses under one set of metadata", { + pool <- crypto_pool() + chosen <- pool$selection[1:2, ] + data <- pool_fit_data(pool_datasets(pool$rows, chosen), chosen) + standardized <- singlehit::as_dose_response(data) + expect_equal(nrow(standardized), length(unique(data$dose))) + expect_equal(sum(standardized$total), sum(data$positive + data$negative)) + 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_library.R b/app/tests/testthat/test-mod_library.R new file mode 100644 index 0000000..8da0e8f --- /dev/null +++ b/app/tests/testthat/test-mod_library.R @@ -0,0 +1,69 @@ +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, args = list(rows = sample_rows()), { + expect_equal(nrow(experiments()), 18L) + expect_equal(nrow(filtered()), 18L) + + session$setInputs(filter_category = "Protozoa") + expect_equal(nrow(filtered()), 7L) + 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, args = list(rows = sample_rows()), { + key <- experiments()$experiment_key[experiments()$pathogen == "Rotavirus"] + session$setInputs(experiment = key, send = 1) + 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()$sent, 2L) + }) +}) + +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$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$sent()) + rows(read_qmra_workbook(qmra_sample_path())) + session$flushReact() + expect_equal(nrow(experiments()), 18L) + }) +}) + +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()), 18L) + }) +}) + +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..3e3eb67 --- /dev/null +++ b/app/tests/testthat/test-mod_pooling.R @@ -0,0 +1,236 @@ +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()) + }) +}) + +crypto_keys <- function(experiments) experiments$experiment_key[grepl("Cryptosporidium", experiments$pathogen)] + +test_that("a group taken from an exploration is marked data-informed and needs a rationale", { + args <- pooling_args() + shiny::testServer(mod_pooling_server, args = args, { + args$pool_set(list(keys = crypto_keys(experiments()), sent = 1L)) + session$flushReact() + session$setInputs(models = c("exponential", "beta_poisson"), alpha = 0.05) + session$setInputs(explore = 1) + expect_false(inherits(exploration()$explore, "error")) + expect_equal(unique(log()$kind), "exploration") + # The first exploration of a Library set is logged as a Library set. + expect_equal(unique(log()$set_origin), "Library") + + choices <- pool_group_choices(exploration()$explore) + session$setInputs(group = choices[[1L]], use_group = 1) + expect_equal(sub("\\|.*", "", keys()), c("108", "139")) + expect_equal(set_origin(), "exploratory grouping (Exponential)") + expect_false(ready()) + + session$setInputs(rationale = "Both Iowa isolate.") + session$setInputs(run = 1) + tests <- log()[log()$kind == "test", ] + expect_equal(unique(tests$set_origin), "exploratory grouping (Exponential)") + + session$setInputs(send = 1) + expect_match(unique(session$returned$sent()$data$study_id), "\\[data-informed set\\]$") + }) +}) + +test_that("re-picking an explored set on the Library tab stays data-informed", { + args <- pooling_args() + shiny::testServer(mod_pooling_server, args = args, { + all <- crypto_keys(experiments()) + args$pool_set(list(keys = all, sent = 1L)) + session$flushReact() + session$setInputs(models = "exponential", alpha = 0.05) + session$setInputs(explore = 1) + + args$pool_set(list(keys = all[1:2], sent = 2L)) + session$flushReact() + expect_equal(set_origin(), "Library, within a set explored this session") + }) +}) + +test_that("pooled data are sent on after a test, with their origin in the label", { + args <- pooling_args() + shiny::testServer(mod_pooling_server, args = args, { + args$pool_set(list(keys = crypto_keys(experiments())[1:2], sent = 1L)) + session$flushReact() + session$setInputs(models = "exponential", alpha = 0.05, rationale = "Both Iowa isolate.") + session$setInputs(send = 1) + expect_null(session$returned$sent()) + + session$setInputs(run = 1) + session$setInputs(send = 2) + sent <- session$returned$sent() + expect_equal(sent$source, "Pooling tab") + expect_match(sent$label, "^Pooled: ") + expect_no_error(singlehit::as_dose_response(sent$data)) + expect_true("sent to Fit" %in% log()$kind) + }) +}) + +test_that("reloading the workbook does not clear what was explored", { + args <- pooling_args() + shiny::testServer(mod_pooling_server, args = args, { + all <- crypto_keys(experiments()) + args$pool_set(list(keys = all, sent = 1L)) + session$flushReact() + session$setInputs(models = "exponential", alpha = 0.05) + session$setInputs(explore = 1) + + args$workbook_loaded(2L) + session$flushReact() + args$pool_set(list(keys = all[1:2], sent = 2L)) + session$flushReact() + expect_equal(set_origin(), "Library, within a set explored this session") + }) +}) + +test_that("the status counts an exploration's pairwise tests and says they are unadjusted", { + args <- pooling_args() + shiny::testServer(mod_pooling_server, args = args, { + args$pool_set(list(keys = crypto_keys(experiments()), sent = 1L)) + session$flushReact() + session$setInputs(models = c("exponential", "beta_poisson"), alpha = 0.05) + expect_false(grepl("adjusted", output$log_status$html)) + session$setInputs(explore = 1) + # Five experiments: 10 pairs under each of two models. + expect_match(output$log_status$html, "20 pairwise tests") + expect_match(output$log_status$html, "None of the p-values are adjusted") + }) +}) + +test_that("the exploration shows how to read the heatmap", { + args <- pooling_args() + shiny::testServer(mod_pooling_server, args = args, { + args$pool_set(list(keys = crypto_keys(experiments()), sent = 1L)) + session$flushReact() + session$setInputs(models = "exponential", alpha = 0.05) + session$setInputs(explore = 1) + html <- output$explore_results$html + expect_match(html, "How to read this heatmap") + expect_match(html, "articles/pooling.html") + }) +}) diff --git a/app/tests/testthat/test-mod_upload.R b/app/tests/testthat/test-mod_upload.R index 1118733..bed49c1 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,22 +37,104 @@ 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) }) }) + +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$data() + expect_equal(result$dose, c(1, 10, 100)) + expect_equal(unique(result$dose_unit), "CFU") + expect_equal(loaded()$label, "Test experiment (from the Library tab)") + + # Loading the example afterwards replaces the library data. + session$setInputs(example = 1) + 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)), 18L) + 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)), 18L) + 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)), 18L) + expect_equal(session$returned$workbook_loaded(), 1L) + + # The error is stale once a dataset loads. + session$setInputs(example = 1) + 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()) + }) +}) diff --git a/data-raw/cryptosporidium-isolates.R b/data-raw/cryptosporidium-isolates.R new file mode 100644 index 0000000..8fdb5c3 --- /dev/null +++ b/data-raw/cryptosporidium-isolates.R @@ -0,0 +1,38 @@ +# Build inst/extdata/cryptosporidium-isolates.csv, the Cryptosporidium human +# challenge experiments used in the pooling article +# (vignettes/articles/pooling.Rmd). They come from the QMRA +# wiki compilation (https://qmrawiki.org), whose source workbook is not in this +# repository; the counts are those reported in the cited studies. +# +# Needs readxl and the QMRA wiki workbook. Run from the package root: +# Rscript data-raw/cryptosporidium-isolates.R "QMRA wiki data_v4.xlsx" +args <- commandArgs(trailingOnly = TRUE) +if (length(args) != 1L) { + stop("Supply the path to the QMRA wiki workbook.") +} + +raw <- readxl::read_excel(args[[1L]], sheet = "raw data", .name_repair = "minimal") +raw <- raw[grepl("Cryptosporidium", raw$pathogen) & raw$experiment_id %in% c("108", "139", "140", "141", "181", "183"), ] + +# "DuPont, H. L., Chappell, ... (1995). The infectivity ..." -> "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/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/articles/pooling-results.rds b/vignettes/articles/pooling-results.rds new file mode 100644 index 0000000..a0611ad Binary files /dev/null and b/vignettes/articles/pooling-results.rds differ diff --git a/vignettes/articles/pooling.Rmd b/vignettes/articles/pooling.Rmd new file mode 100644 index 0000000..3cfb943 --- /dev/null +++ b/vignettes/articles/pooling.Rmd @@ -0,0 +1,358 @@ +--- +title: "Pooling experiments: tests, heatmaps, and forking paths" +--- + +```{r, include = FALSE} +knitr::opts_chunk$set(collapse = TRUE, comment = "#>", fig.width = 9, fig.height = 4.5) +``` + +When several challenge studies exist for one pathogen, it is tempting to +combine them into one dose-response curve: more subjects, narrower intervals, +one number to report. singlehit's pooling functions test whether that is +defensible. This article explains what the tests do, how to read their +results, and how to avoid the main trap, which is letting the test results +choose the groups. + +```{r setup} +library(singlehit) +``` + +## The data + +The bundled file holds six human Cryptosporidium challenge experiments +compiled by the [QMRA wiki](https://qmrawiki.org). Every experiment used +healthy adult volunteers, oral exposure, and doses counted in oocysts. Most +differ in the parasite isolate, and one records a different outcome. + +```{r data} +path <- system.file("extdata", "cryptosporidium-isolates.csv", package = "singlehit") +raw <- readr::read_csv(path, show_col_types = FALSE, col_types = readr::cols(study_id = "c")) +unique(raw[c("study_id", "isolate", "endpoint", "source")]) +trials <- split(raw, raw$study_id) +``` + +Each experiment has four to eight dose groups with a handful of volunteers in +each. That small size matters for everything below. + +## What a poolability test asks + +`poolability_test()` runs the likelihood-ratio test described by Haas, Rose, +and Gerba. For a set of experiments it: + +1. fits each experiment separately and adds up their deviances (the lack of + fit); +2. stacks every experiment's dose groups and fits one shared curve; +3. takes the difference in deviance as the test statistic. Sharing one curve + can only fit as well or worse, so the statistic is never negative; +4. compares it with a chi-squared distribution whose degrees of freedom are + the parameters saved by sharing: one per extra experiment for the + exponential model (*k*), two per extra experiment for beta-Poisson + (*alpha* and *N50*). + +The p-value answers one narrow question: *if these experiments really shared +one curve, how surprising would this much extra lack of fit be?* A small +p-value is evidence that they do not share a curve. A large one only means the +data could not tell them apart. + +## Decide the groups first + +The grouping should come from what the experiments are, decided before +looking at any test: same pathogen and isolate, same outcome, same host and +route, same dose unit. For these data the isolate is the obvious candidate, +because isolates of Cryptosporidium are known to differ in infectivity. + +singlehit enforces one of these conditions itself. Doses in different units +are not comparable, so datasets that record different `dose_unit` values are +refused: + +```{r unit-guard} +relabelled <- trials[c("139", "140")] +relabelled[["140"]]$dose_unit <- "genome copies" +tryCatch(poolability_test(relabelled), error = conditionMessage) +``` + +Hosts, routes, isolates, and outcomes are not checked. Comparing them is your +job. + +## A prespecified test + +Messner and colleagues (2001) analysed three isolates challenged in the same +programme: Iowa (139), TAMU (140), and UCP (141). A natural prespecified +question is whether the three isolates share one curve. + +```{r confirmatory} +isolates <- poolability_test(trials[c("139", "140", "141")]) +isolates[c("model", "lrt_statistic", "df", "p_value", "converged")] +``` + +```{r confirmatory-values, include = FALSE} +fp <- function(p) if (p < 0.001) "p < 0.001" else paste("p =", signif(p, 2)) +p_iso_exp <- isolates$p_value[isolates$model == "exponential"] +p_iso_bp <- isolates$p_value[isolates$model == "beta_poisson"] +``` + +Pooling is rejected under both models (exponential `r fp(p_iso_exp)`, +beta-Poisson `r fp(p_iso_bp)`). The +separate fits show why: the isolates differ several-fold in median +infectious dose. + +```{r separate-fits} +n50 <- vapply(c("139", "140", "141"), function(id) { + fit <- fit_dose_response(trials[[id]], "beta_poisson") + generics::tidy(fit)$estimate[generics::tidy(fit)$term == "n50"] +}, numeric(1)) +round(n50) +``` + +The defensible conclusion is to model the isolates separately, or to model +the variation between them explicitly (Messner and colleagues used a +hierarchical model for exactly this), rather than to report one pooled curve. +Had the test not rejected pooling, `fit_pooled_dose_response()` would fit the +shared curve to the stacked dose groups. + +## Exploring every pair + +Sometimes you want to see the structure first: which experiments look +alike? `group_datasets()` tests every pair and then proposes a grouping. The +pairwise p-values are attached to its result. + +The two `group_datasets()` searches in this article take about a minute +together, so their results were computed ahead of time by +`data-raw/pooling-article-results.R` and are read from a saved file; everything +else runs as the article is built. + +```{r frozen, include = FALSE} +frozen <- readRDS("pooling-results.rds") +if (!isTRUE(all.equal(frozen$data, as.data.frame(raw), check.attributes = FALSE))) { + stop("The example data changed: rerun data-raw/pooling-article-results.R.") +} +``` + +```{r explore, eval = FALSE} +# Collect the convergence warnings rather than printing one per refit. +warnings <- character() +explore <- withCallingHandlers( + group_datasets(trials, method = "agglomerative"), + warning = function(w) { + warnings <<- c(warnings, conditionMessage(w)) + invokeRestart("muffleWarning") + } +) +unique(warnings) +pairwise <- attr(explore, "pairwise") +``` + +```{r explore-frozen, echo = FALSE} +explore <- frozen$explore +warnings <- frozen$warnings +pairwise <- attr(explore, "pairwise") +# A cheap live test of one pair, so a change in the pooling code that moves +# the saved results stops the build instead of passing silently. +live <- suppressWarnings(poolability_test(trials[c("108", "139")])) +live_p <- stats::setNames(live$p_value, live$model) +# The tolerance allows for small optimizer differences between platforms. +if (!isTRUE(all.equal(live_p[["exponential"]], pairwise$exponential["108", "139"], tolerance = 1e-4)) || + !isTRUE(all.equal(live_p[["beta_poisson"]], pairwise$beta_poisson["108", "139"], tolerance = 1e-4))) { + stop("The saved results no longer match the pooling code: rerun data-raw/pooling-article-results.R.") +} +unique(warnings) +``` + +```{r heatmap, fig.height = 4.8} +cells <- do.call(rbind, lapply(names(pairwise), function(model) { + p <- pairwise[[model]] + grid <- expand.grid(a = rownames(p), b = colnames(p), stringsAsFactors = FALSE) + grid$p <- p[cbind(grid$a, grid$b)] + grid$model <- model + grid[grid$a != grid$b, ] +})) +cells$model <- factor(cells$model, levels = c("exponential", "beta_poisson"), labels = c("Exponential", "Beta-Poisson")) +cells$b <- factor(cells$b, levels = rev(sort(unique(cells$b)))) + +ggplot2::ggplot(cells, ggplot2::aes(a, b, fill = p)) + + ggplot2::geom_tile(colour = "white") + + ggplot2::geom_text(ggplot2::aes(label = ifelse(p < 0.001, "<0.001", signif(p, 2))), size = 3.2) + + ggplot2::scale_fill_gradient(low = "#F4E1D2", high = "#0066A1", limits = c(0, 1), name = "Pairwise p") + + ggplot2::facet_wrap(ggplot2::vars(model)) + + ggplot2::labs(x = NULL, y = NULL) + + ggplot2::theme_minimal() +``` + +```{r explore-values, include = FALSE} +pe <- pairwise$exponential +pb <- pairwise$beta_poisson +others_140 <- setdiff(rownames(pe), "140") +mean_off <- function(p) mean(p[upper.tri(p)]) +``` + +### Reading one cell + +Each square is one poolability test of two experiments. **Pale** means a +small p-value: one shared curve fits clearly worse, so there is evidence the +two differ. **Blue** means a large p-value: no evidence of a difference. The +matrix is symmetric, and the diagonal is blank because an experiment always +pools with itself. The grid reports no convergence flags, so read it +alongside the warnings above: any cell can come from a fit that did not +converge. + +In the exponential panel: + +- **140 (TAMU) is pale against everything** (largest p = + `r signif(max(pe["140", others_140]), 2)`). TAMU is far more infectious than + the other isolates, and no shared curve can fit it with any of them. +- **108 and 139 (both Iowa) give p = `r signif(pe["108", "139"], 2)`.** That + looks like two studies agreeing, but see "Overlapping subjects" below. +- **181 against 139 gives p = `r signif(pe["139", "181"], 2)`.** Experiment + 181 is a different isolate (TU502) *and records a different outcome* + (diarrhea rather than infection). It should never be pooled with infection + data, yet the test does not reject it. + +### Why the beta-Poisson panel is bluer + +The same pairs give larger p-values under beta-Poisson (mean +`r signif(mean_off(pb), 2)` against `r signif(mean_off(pe), 2)` for +exponential). The experiments did not become more alike; the test became +weaker. The second parameter lets each curve bend to fit its own data, the +test has two degrees of freedom instead of one (a higher bar to reject), and +with four to eight doses per experiment the approximate beta-Poisson model is +poorly determined. Beta-Poisson comparisons of small studies typically have +little power. + +## Three mistakes to avoid + +**"Not rejected" is not "the same".** A p-value of +`r signif(pe["139", "181"], 2)` does not show that experiments 139 and 181 +share a curve; it shows these data cannot tell them apart. With so few +volunteers, very different experiments often pass. + +**Pairwise results do not chain.** "A pools with B" and "B pools with C" does +not mean A, B, and C share one curve. A group needs one curve for all its +members, which is a stricter test than any pair. + +**Many tests, no correction.** Each panel holds `r choose(nrow(pe), 2)` +tests, `r 2 * choose(nrow(pe), 2)` in all. At alpha = 0.05 you would expect +a "rejection" or two by chance even if every experiment were identical, and +pairs can pass just as arbitrarily. None of these p-values are adjusted. + +## Groupings are not unique + +`group_datasets()` turns the pairwise results into groups. The agglomerative +method (above) repeatedly merges the most compatible pair of groups, testing +each merged group as a whole, until no merge passes; the exhaustive +method searches every way of splitting the experiments (203 ways for six) for +the fewest groups that all pass. + +```{r exhaustive-code, eval = FALSE} +exhaustive <- suppressWarnings(group_datasets(trials, method = "exhaustive")) +``` + +```{r exhaustive-frozen, include = FALSE} +exhaustive <- frozen$exhaustive +``` + +```{r exhaustive} +groups <- function(result, model) { + rows <- result[result$model == model, ] + vapply(split(rows$dataset, rows$group), paste, character(1), collapse = " + ") +} +list( + agglomerative = list(exponential = groups(explore, "exponential"), beta_poisson = groups(explore, "beta_poisson")), + exhaustive = list(exponential = groups(exhaustive, "exponential"), beta_poisson = groups(exhaustive, "beta_poisson")) +) +``` + +```{r groups-values, include = FALSE} +together <- function(result, model, a, b) { + rows <- result[result$model == model, ] + rows$group[rows$dataset == a] == rows$group[rows$dataset == b] +} +``` + +The same data give different groupings depending on the model and the +search method. Under beta-Poisson, the exhaustive search +`r if (together(exhaustive, "beta_poisson", "108", "139")) "keeps" else "separates"` +the two Iowa experiments, even though the heatmap shows no difference +between them (p = `r signif(pb["108", "139"], 2)`). Several groupings pass the tests, and the search reports one of +them. In both methods the TU502 diarrhea experiment (181) joins an infection +group. + +This is the forking-paths problem +([issue #20](https://github.com/seanthimons/singlehit/issues/20)): when the data choose the +groups, a pooled curve fitted to a chosen group, and its test, are biased +toward pooling. Use the exploration to sanity-check a grouping you decided +on, not to choose one. + +## Overlapping subjects + +Look at experiments 108 and 139 side by side: + +```{r overlap} +overlap <- merge( + trials[["108"]][c("dose", "positive", "negative")], + trials[["139"]][c("dose", "positive", "negative")], + by = "dose", suffixes = c("_108", "_139") +) +overlap +``` + +```{r overlap-values, include = FALSE} +same_rows <- sum(overlap$positive_108 == overlap$positive_139 & overlap$negative_108 == overlap$negative_139) +``` + +The counts are identical at `r same_rows` of `r nrow(overlap)` doses and +differ by one volunteer at the others. That strongly suggests both entries +report largely the same Iowa volunteers: DuPont and colleagues (1995) +reported an Iowa-isolate challenge, and Messner and colleagues (2001) +analysed Iowa-isolate data from the same investigators (Chappell and +Okhuysen are authors of both). The poolability test +assumes independent samples, so its p-value of +`r signif(pe["108", "139"], 2)` says nothing here, and pooling the two would +count the same people twice. Check where each experiment's subjects came +from before pooling. + +## Reporting + +A paragraph that reports this honestly might read: + +> 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 728e4e0..9648255 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 @@ -312,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.