Unit testing and examples (#70)

* Updated calculate_response

* Added plot_hill tests

* Updated plot_exposure

* Updated resp_quantiles

* Added missing test case

* Fixed error due to zero-space char in example

* Updated plot_resp

* Avoid using chemical names in examples

* Changed all_of to any_of

DESCRIPTION

@@ -18,6 +18,7 @@ Imports:
     rlang,
     sf,
     stats,
+    stringr,
     tibble,
     tidyr,
     tidyselect,
@@ -33,7 +34,6 @@ Suggests:
     readxl,
     rmarkdown,
     scales,
-    stringr,
     testthat (>= 3.0.0)
 VignetteBuilder: 
     knitr

R/calculate_response.R

@@ -14,6 +14,25 @@
 #'
 #' @seealso [calc_internal_dose], [calc_invitro_concentration],
 #' [calc_concentration_response]
+#' 
+#' @examples
+#' # Create GeoTox object
+#' geoTox <- GeoTox()
+#'
+#' # The following fields are required inputs
+#' geoTox$IR <- 2
+#' geoTox$C_ext <- matrix(3)
+#' geoTox$C_ss <- 5
+#' geoTox$hill_params <- fit_hill(data.frame(logc = c(-1, 0, 1),
+#'                                           resp = c(10, 5, 0)))
+#'
+#' # Calculate response
+#' geoTox <- calculate_response(geoTox)
+#' 
+#' # The following fields will be computed
+#' geoTox$D_int
+#' geoTox$C_invitro
+#' geoTox$resp
 calculate_response <- function(x, ...) {
 
   # Update parameters
@@ -40,7 +59,7 @@ calculate_response <- function(x, ...) {
   x$C_invitro <- calc_invitro_concentration(x$D_int, x$C_ss)
 
   # Concentration response
-  if (is.null(x$C_ss)) {
+  if (is.null(x$hill_params)) {
     stop("GeoTox object must contain 'hill_params' field", call. = FALSE)
   }
   x$resp <- calc_concentration_response(x$C_invitro,

R/plot_exposure.R

@@ -1,6 +1,6 @@
 #' Plot exposure data.
 #'
-#' @param exposure exposure data.
+#' @param exposure list of exposure data named by region label.
 #' @param region_boundary "sf" data.frame mapping features to a "geometry"
 #' column. Used to color regions.
 #' @param group_boundary (optional) "sf" data.frame containing a "geometry"
@@ -10,6 +10,26 @@
 #'
 #' @return ggplot2 object.
 #' @export
+#' 
+#' @examples
+#' # Load package data
+#' exposure <- split(geo_tox_data$exposure, ~FIPS)
+#' region_boundary <- geo_tox_data$boundaries$county
+#' group_boundary <- geo_tox_data$boundaries$state
+#' 
+#' # Plot county exposure data
+#' # Use CASN as label to avoid long chemical names
+#' plot_exposure(exposure,
+#'               region_boundary,
+#'               chem_label = "casn",
+#'               ncol = 5)
+#'
+#' # Add state boundaries
+#' plot_exposure(exposure,
+#'               region_boundary,
+#'               group_boundary = group_boundary,
+#'               chem_label = "casn",
+#'               ncol = 5)
 plot_exposure <- function(exposure,
                           region_boundary,
                           group_boundary = NULL,
@@ -19,24 +39,46 @@ plot_exposure <- function(exposure,
   if (is.null(exposure)) {
     stop("No exposure data found.", call. = FALSE)
   }
+  if (is.null(region_boundary)) {
+    stop("No region_boundary data found.", call. = FALSE)
+  }
 
-  df <- tibble::tibble(id = names(exposure), data = exposure) |> 
+  df <- tibble::tibble("_temp_join_id_" = names(exposure), data = exposure) |> 
     tidyr::unnest(cols = "data") |> 
-    dplyr::inner_join(region_boundary |> dplyr::rename("id" = 1),
-                      by = dplyr::join_by("id"))
+    dplyr::inner_join(region_boundary |> dplyr::rename("_temp_join_id_" = 1),
+                      by = dplyr::join_by("_temp_join_id_")) |> 
+    dplyr::select(-"_temp_join_id_") |> 
+    # Fix for grid.Call error in examples due to
+    # 'mbcsToSbcs': for ​ (U+200B)
+    # Remove any zero-width space characters
+    dplyr::mutate(dplyr::across(tidyselect::any_of(chem_label),
+                                ~ stringr::str_remove_all(., "\u200b")))
 
-  ggplot2::ggplot(df, ggplot2::aes(fill = .data$norm)) +
-    ggplot2::geom_sf(ggplot2::aes(geometry = .data$geometry), color = NA) +
+  fig <- ggplot2::ggplot() +
+    ggplot2::geom_sf(
+      data = df,
+      ggplot2::aes(fill = .data$norm,
+                   geometry = .data$geometry),
+      color = NA) +
     ggplot2::facet_wrap(chem_label, ncol = ncol) +
     # Recolor subset as light grey
     ggplot2::geom_sf(
       data = df |> dplyr::filter(mean == 0),
       ggplot2::aes(geometry = .data$geometry),
       fill = "light grey",
       color = "light grey",
-      lwd = 0.01) +
-    # State borders
-    ggplot2::geom_sf(data = group_boundary, fill = NA, size = 0.15) +
+      lwd = 0.01)
+
+  # State borders
+  if (!is.null(group_boundary)) {
+    fig <- fig +
+      ggplot2::geom_sf(data = group_boundary,
+                       ggplot2::aes(geometry = .data$geometry),
+                       fill = NA,
+                       size = 0.15)
+  }
+
+  fig +
     ggplot2::scale_fill_viridis_c(
       name = "Normalized\nConcentration",
       direction = -1,
@@ -50,4 +92,4 @@ plot_exposure <- function(exposure,
       text = ggplot2::element_text(size = 12),
       panel.grid.major = ggplot2::element_blank(),
       panel.grid.minor = ggplot2::element_blank())
-}
+}

R/plot_resp.R

@@ -10,6 +10,44 @@
 #'
 #' @return ggplot2 object.
 #' @export
+#' 
+#' @examples
+#' # Use example boundary data from package
+#' region_boundary <- geo_tox_data$boundaries$county
+#' group_boundary <- geo_tox_data$boundaries$state
+#' n <- nrow(region_boundary)
+#' 
+#' # Single assay quantile
+#' df <- data.frame(id = region_boundary$FIPS,
+#'                  metric = "GCA.Eff",
+#'                  assay_quantile = 0.5,
+#'                  value = runif(n)^3)
+#' # Default plot
+#' plot_resp(df, region_boundary)
+#' # Add group boundary, a state border in this case
+#' plot_resp(df, region_boundary, group_boundary)
+#' # Change quantile label
+#' plot_resp(df, region_boundary, group_boundary,
+#'           assay_quantiles = c("Q50" = 0.5))
+#' 
+#' # Multiple assay quantiles
+#' df <- data.frame(id = rep(region_boundary$FIPS, 2),
+#'                  metric = "GCA.Eff",
+#'                  assay_quantile = rep(c(0.25, 0.75), each = n),
+#'                  value = c(runif(n)^3, runif(n)^3 + 0.15))
+#' plot_resp(df, region_boundary, group_boundary,
+#'           assay_quantiles = c("Q25" = 0.25, "Q75" = 0.75))
+#' 
+#' # Summary quantiles
+#' df <- data.frame(id = rep(region_boundary$FIPS, 4),
+#'                  assay_quantile = rep(rep(c(0.25, 0.75), each = n), 2),
+#'                  summary_quantile = rep(c(0.05, 0.95), each = n * 2),
+#'                  metric = "GCA.Eff",
+#'                  value = c(runif(n)^3, runif(n)^3 + 0.15,
+#'                            runif(n)^3 + 0.7, runif(n)^3 + 0.85))
+#' plot_resp(df, region_boundary, group_boundary,
+#'           assay_quantiles = c("A_Q25" = 0.25, "A_Q75" = 0.75),
+#'           summary_quantiles = c("S_Q05" = 0.05, "S_Q95" = 0.95))
 plot_resp <- function(
     df,
     region_boundary,
@@ -45,9 +83,12 @@ plot_resp <- function(
 
   metric <- df$metric[1]
 
-  fig <- ggplot2::ggplot(df, ggplot2::aes(fill = .data$value)) +
+  fig <- ggplot2::ggplot() +
     # Plot county data using fill, hide county borders by setting color = NA
-    ggplot2::geom_sf(ggplot2::aes(geometry = .data$geometry), color = NA) +
+    ggplot2::geom_sf(data = df,
+                     ggplot2::aes(fill = .data$value,
+                                  geometry = .data$geometry),
+                     color = NA) +
     # Add fill scale
     ggplot2::scale_fill_viridis_c(
       name = metric,
@@ -91,15 +132,18 @@ plot_resp <- function(
       # Create separate plots for each stat
       ggplot2::facet_wrap(
         ~assay_quantile,
-        ncol = length(assay_quantiles),
+        ncol = length(unique(df$assay_quantile)),
         labeller = ggplot2::labeller(
           assay_quantile = stats::setNames(names(assay_quantiles),
                                            assay_quantiles)))
   }
 
   if (!is.null(group_boundary)) {
-    fig <- fig + ggplot2::geom_sf(data = group_boundary, fill = NA,
-                                  size = 0.15)
+    fig <- fig + 
+      ggplot2::geom_sf(data = group_boundary,
+                       ggplot2::aes(geometry = .data$geometry),
+                       fill = NA,
+                       size = 0.15)
   }
 
   fig

R/resp_quantiles.R

@@ -1,6 +1,7 @@
 #' Get response quantiles
 #'
-#' @param resp calculated mixture response output from [calculate_response].
+#' @param resp calculated mixture response output from
+#' [calc_concentration_response].
 #' @param metric response metric, one of "GCA.Eff", "IA.Eff", "GCA.HQ.10"
 #' or "IA.HQ.10".
 #' @param assays assays to summarize. If NULL and multiple assays exist, then
@@ -18,6 +19,35 @@
 #'
 #' @return data frame with computed response quantiles.
 #' @export
+#' 
+#' @examples
+#' # Dummy response data
+#' resp <- list(
+#'   "r1" = data.frame(assay = c("a1", "a1", "a2", "a2"),
+#'                     sample = c(1, 2, 1, 2),
+#'                     GCA.Eff = c(1, 2, 3, 4),
+#'                     IA.Eff = c(5, 6, 7, 8),
+#'                     "GCA.HQ.10" = c(9, 10, 11, 12),
+#'                     "IA.HQ.10" = c(13, 14, 15, 16)))
+#'
+#' # Summarize single assay
+#' resp_quantiles(resp)
+#' # Specify assay
+#' resp_quantiles(resp, assays = "a1")
+#' # Specify quantiles
+#' resp_quantiles(resp, assays = "a1", assay_quantiles = c(0.25, 0.75))
+#' # Specify metric
+#' resp_quantiles(resp, assays = "a1", metric = "IA.HQ.10")
+#' 
+#' # Summarize across assays
+#' resp_quantiles(resp, assay_summary = TRUE)
+#' # Specify quantiles
+#' suppressWarnings(
+#'   resp_quantiles(resp,
+#'                  assay_summary = TRUE,
+#'                  assay_quantiles = c(0.25, 0.75),
+#'                  summary_quantiles = c(0.1, 0.9))
+#' )
 resp_quantiles <- function(
     resp,
     metric = c("GCA.Eff", "IA.Eff", "GCA.HQ.10", "IA.HQ.10"),