From 05a941c722f7b1625897a645f3141e3c7c691c6b Mon Sep 17 00:00:00 2001
From: SkylarMarvel <57526164+SkylarMarvel@users.noreply.github.com>
Date: Mon, 4 Nov 2024 16:45:28 -0700
Subject: [PATCH] Unit tests (#71)

* Updated plot_sensitivity

* Added example for plot_sensitivity

* Added unit test for GeoTox S3 object
---
 R/GeoTox.R                                |  42 +++++++-
 R/plot_sensitivity.R                      |  22 +++++
 _pkgdown.yml                              |   3 +-
 man/GeoTox.Rd                             |  40 ++++++++
 man/plot_sensitivity.Rd                   |  19 ++++
 tests/testthat/test-GeoTox.R              | 112 ++++++++++++++++++++++
 tests/testthat/test-plot_sensitivity.R    |  89 +++++++++++++++++
 tests/testthat/test-simulate_population.R |  19 ----
 8 files changed, 324 insertions(+), 22 deletions(-)
 create mode 100644 tests/testthat/test-GeoTox.R
 create mode 100644 tests/testthat/test-plot_sensitivity.R

diff --git a/R/GeoTox.R b/R/GeoTox.R
index 5effaf5..1c61f64 100644
--- a/R/GeoTox.R
+++ b/R/GeoTox.R
@@ -5,6 +5,46 @@
 #' GeoTox analysis.
 #'
 #' @export
+#' 
+#' @examples
+#' # Use a subset of the package data for demonstration purposes
+#' set.seed(2357)
+#' n <- 10 # Population size
+#' m <- 5 # Number of regions
+#' idx <- if (m < 100) sample(1:100, m) else 1:100
+#' 
+#' geoTox <- GeoTox() |> 
+#'   # Set region and group boundaries (for plotting)
+#'   set_boundaries(region = geo_tox_data$boundaries$county,
+#'                  group  = geo_tox_data$boundaries$state) |> 
+#'   # Simulate populations for each region
+#'   simulate_population(age           = split(geo_tox_data$age, ~FIPS)[idx],
+#'                       obesity       = geo_tox_data$obesity[idx, ],
+#'                       exposure      = split(geo_tox_data$exposure, ~FIPS)[idx],
+#'                       simulated_css = geo_tox_data$simulated_css,
+#'                       n             = n) |> 
+#'   # Estimated Hill parameters
+#'   set_hill_params(geo_tox_data$dose_response |>
+#'                     fit_hill(assay = "endp", chem = "casn") |> 
+#'                     dplyr::filter(!tp.sd.imputed, !logAC50.sd.imputed)) |>
+#'   # Calculate response
+#'   calculate_response() |>
+#'   # Perform sensitivity analysis
+#'   sensitivity_analysis()
+#' 
+#' # Print GeoTox object
+#' geoTox
+#' 
+#' # Plot hill fits
+#' plot(geoTox, type = "hill")
+#' # Plot exposure data
+#' plot(geoTox, type = "exposure")
+#' # Plot response data
+#' plot(geoTox)
+#' plot(geoTox, assays = "TOX21_H2AX_HTRF_CHO_Agonist_ratio")
+#' # Plot sensitivity data
+#' plot(geoTox, type = "sensitivity")
+#' plot(geoTox, type = "sensitivity", assay = "TOX21_H2AX_HTRF_CHO_Agonist_ratio")
 GeoTox <- function() {
   structure(
     list(
@@ -169,7 +209,5 @@ plot.GeoTox <- function(x,
                   ncol            = ncol)
   } else if (type == "sensitivity") {
     plot_sensitivity(x, ...)
-  } else {
-    stop("Invalid type.")
   }
 }
diff --git a/R/plot_sensitivity.R b/R/plot_sensitivity.R
index c286791..e40f018 100644
--- a/R/plot_sensitivity.R
+++ b/R/plot_sensitivity.R
@@ -10,6 +10,25 @@
 #'
 #' @return ggplot2 object.
 #' @export
+#' 
+#' @examples
+#' # Required GeoTox fields are generated by first running [calculate_response]
+#' # and [sensitivity_analysis] on a GeoTox object. This will create the fields
+#' # `resp` and `sensitivity`. For this example, dummy data will be used.
+#' make_data <- function(n = 5, metric = "GCA.Eff") {
+#'   list(stats::setNames(data.frame(1:n, runif(n)),
+#'                        c("sample", metric)))
+#' }
+#' 
+#' x <- GeoTox()
+#' x$resp <- make_data()
+#' x$sensitivity <- list(age = make_data(),
+#'                       obesity = make_data(),
+#'                       css_params = make_data(),
+#'                       fit_params = make_data(),
+#'                       C_ext = make_data())
+#' 
+#' plot_sensitivity(x)
 plot_sensitivity <- function(x,
                              metric = "GCA.Eff",
                              assay = NULL,
@@ -20,6 +39,9 @@ plot_sensitivity <- function(x,
   if (is.null(x$sensitivity)) {
     stop("No sensitivity data found.", call. = FALSE)
   }
+  if (is.null(x$resp)) {
+    stop("No baseline response data found.", call. = FALSE)
+  }
   
   df <- get_sensitivity_df(x, metric = metric, assay = assay)
   fig <- plot_sensitivity_df(df, y = y, xlab = xlab, ylab = ylab)
diff --git a/_pkgdown.yml b/_pkgdown.yml
index 9dd1e2e..75651b6 100644
--- a/_pkgdown.yml
+++ b/_pkgdown.yml
@@ -18,6 +18,8 @@ reference:
   - compute_sensitivity
   - plot_sensitivity
   - sensitivity_analysis
+  - set_boundaries
+  - set_hill_params
   - simulate_population
 - subtitle: Simulate
   desc: Functions that perform Monte Carlo simulation
@@ -45,7 +47,6 @@ reference:
   - starts_with("get_")
   - starts_with("hill_")
   - starts_with("obj_")
-  - starts_with("set_")
 - subtitle: Datasets
   desc: Package datasets
 - contents:
diff --git a/man/GeoTox.Rd b/man/GeoTox.Rd
index 55dad7c..0e43924 100644
--- a/man/GeoTox.Rd
+++ b/man/GeoTox.Rd
@@ -20,6 +20,46 @@ GeoTox()
 An S3 object that can be used to help organize the data and results of a
 GeoTox analysis.
 }
+\examples{
+# Use a subset of the package data for demonstration purposes
+set.seed(2357)
+n <- 10 # Population size
+m <- 5 # Number of regions
+idx <- if (m < 100) sample(1:100, m) else 1:100
+
+geoTox <- GeoTox() |> 
+  # Set region and group boundaries (for plotting)
+  set_boundaries(region = geo_tox_data$boundaries$county,
+                 group  = geo_tox_data$boundaries$state) |> 
+  # Simulate populations for each region
+  simulate_population(age           = split(geo_tox_data$age, ~FIPS)[idx],
+                      obesity       = geo_tox_data$obesity[idx, ],
+                      exposure      = split(geo_tox_data$exposure, ~FIPS)[idx],
+                      simulated_css = geo_tox_data$simulated_css,
+                      n             = n) |> 
+  # Estimated Hill parameters
+  set_hill_params(geo_tox_data$dose_response |>
+                    fit_hill(assay = "endp", chem = "casn") |> 
+                    dplyr::filter(!tp.sd.imputed, !logAC50.sd.imputed)) |>
+  # Calculate response
+  calculate_response() |>
+  # Perform sensitivity analysis
+  sensitivity_analysis()
+
+# Print GeoTox object
+geoTox
+
+# Plot hill fits
+plot(geoTox, type = "hill")
+# Plot exposure data
+plot(geoTox, type = "exposure")
+# Plot response data
+plot(geoTox)
+plot(geoTox, assays = "TOX21_H2AX_HTRF_CHO_Agonist_ratio")
+# Plot sensitivity data
+plot(geoTox, type = "sensitivity")
+plot(geoTox, type = "sensitivity", assay = "TOX21_H2AX_HTRF_CHO_Agonist_ratio")
+}
 \seealso{
 \link{plot_resp}, \link{plot_hill}, \link{plot_exposure}, \link{plot_sensitivity}
 }
diff --git a/man/plot_sensitivity.Rd b/man/plot_sensitivity.Rd
index 227363f..ebf3295 100644
--- a/man/plot_sensitivity.Rd
+++ b/man/plot_sensitivity.Rd
@@ -33,3 +33,22 @@ ggplot2 object.
 \description{
 Plot results of sensitivity analysis.
 }
+\examples{
+# Required GeoTox fields are generated by first running [calculate_response]
+# and [sensitivity_analysis] on a GeoTox object. This will create the fields
+# `resp` and `sensitivity`. For this example, dummy data will be used.
+make_data <- function(n = 5, metric = "GCA.Eff") {
+  list(stats::setNames(data.frame(1:n, runif(n)),
+                       c("sample", metric)))
+}
+
+x <- GeoTox()
+x$resp <- make_data()
+x$sensitivity <- list(age = make_data(),
+                      obesity = make_data(),
+                      css_params = make_data(),
+                      fit_params = make_data(),
+                      C_ext = make_data())
+
+plot_sensitivity(x)
+}
diff --git a/tests/testthat/test-GeoTox.R b/tests/testthat/test-GeoTox.R
new file mode 100644
index 0000000..a8f7d7c
--- /dev/null
+++ b/tests/testthat/test-GeoTox.R
@@ -0,0 +1,112 @@
+test_that("errors", {
+  
+  geoTox <- GeoTox()
+  
+  expect_error(plot(geoTox),
+               "No response data found")
+  
+  geoTox$resp <- list()
+  
+  expect_error(plot(geoTox),
+               "No region boundary data found")
+})
+
+test_that("default params", {
+  
+  geoTox <- GeoTox()
+  
+  expect_equal(geoTox$par,
+               list(n = 1e3,
+                    IR_params = NULL,
+                    obesity = list(obes_prev  = "OBESITY_CrudePrev",
+                                   obes_sd    = "OBESITY_SD",
+                                   obes_label = "FIPS"),
+                    exposure = list(expos_mean  = "mean",
+                                    expos_sd    = "sd",
+                                    expos_label = "casn"),
+                    internal_dose = list(time    = 1,
+                                         BW      = 1,
+                                         scaling = 1),
+                    resp = list(max_mult = 1.5)))
+})
+
+test_that("package data subset", {
+  
+  set.seed(2357)
+  n <- 10 # Population size
+  m <- 5 # Number of regions
+  idx <- if (m < 100) sample(1:100, m) else 1:100
+
+  expect_no_error(
+    geoTox <- GeoTox() |>
+      # Set region and group boundaries (for plotting)
+      set_boundaries(region = geo_tox_data$boundaries$county,
+                     group  = geo_tox_data$boundaries$state) |>
+      # Simulate populations for each region
+      simulate_population(age           = split(geo_tox_data$age, ~FIPS)[idx],
+                          obesity       = geo_tox_data$obesity[idx, ],
+                          exposure      = split(geo_tox_data$exposure, ~FIPS)[idx],
+                          simulated_css = geo_tox_data$simulated_css,
+                          n             = n) |> 
+      # Estimated Hill parameters
+      set_hill_params(geo_tox_data$dose_response |>
+                        fit_hill(assay = "endp", chem = "casn") |>
+                        dplyr::filter(!tp.sd.imputed, !logAC50.sd.imputed)) |>
+      # Calculate response
+      calculate_response() |>
+      # Perform sensitivity analysis
+      sensitivity_analysis()
+  )
+  
+  # Print GeoTox object
+  expect_no_error(capture_output(print(geoTox)))
+  
+  # Plot hill fits
+  expect_no_error(plot(geoTox, type = "hill"))
+  # Plot exposure data
+  expect_no_error(plot(geoTox, type = "exposure"))
+  # Plot response data
+  # Warning from resp_quantiles()
+  expect_warning(plot(geoTox),
+                 "Multiple assays found, using first assay")
+  expect_no_error(plot(geoTox, assays = "TOX21_H2AX_HTRF_CHO_Agonist_ratio"))
+  # Plot sensitivity data
+  expect_warning(expect_warning(
+    plot(geoTox, type = "sensitivity"),
+    "Multiple assays found, using first assay"),
+    "Removed \\d+ NA value")
+  expect_no_error(plot(geoTox,
+                       type = "sensitivity",
+                       assay = "TOX21_H2AX_HTRF_CHO_Agonist_ratio"))
+})
+
+test_that("print corner cases", {
+  
+  geoTox <- GeoTox()
+  
+  expect_no_error(capture_output(print(geoTox)))
+  
+  geoTox$age <- list()
+  geoTox$obesity <- c(0)
+  geoTox$par <- NULL
+  
+  expect_no_error(capture_output(print(geoTox)))
+  
+  geoTox <- GeoTox() |>
+    set_hill_params(geo_tox_data$dose_response |>
+                      dplyr::filter(endp == "TOX21_H2AX_HTRF_CHO_Agonist_ratio",
+                                    casn == "510-15-6") |>
+                      fit_hill())
+
+  expect_no_error(capture_output(print(geoTox)))
+  
+  geoTox <- geoTox |>
+    set_hill_params(geo_tox_data$dose_response |>
+                      fit_hill(assay = "endp", chem = "casn"))
+  
+  expect_no_error(capture_output(print(geoTox)))
+  
+  geoTox$resp <- c(1:5)
+  
+  expect_no_error(capture_output(print(geoTox)))
+})
diff --git a/tests/testthat/test-plot_sensitivity.R b/tests/testthat/test-plot_sensitivity.R
new file mode 100644
index 0000000..ebf92c4
--- /dev/null
+++ b/tests/testthat/test-plot_sensitivity.R
@@ -0,0 +1,89 @@
+test_that("errors", {
+  
+  geoTox <- GeoTox()
+  
+  expect_error(plot_sensitivity(geoTox),
+               "No sensitivity data found")
+  
+  geoTox$sensitivity <- list()
+  
+  expect_error(plot_sensitivity(geoTox),
+               "No baseline response data found")
+})
+
+test_that("no assay", {
+  
+  make_data <- function(n = 5, metric = "GCA.Eff") {
+    list(stats::setNames(data.frame(1:n, runif(n)),
+                         c("sample", metric)))
+  }
+  
+  x <- GeoTox()
+  n <- 5
+  x$resp <- make_data(n)
+  x$sensitivity <- list(age = make_data(n),
+                        obesity = make_data(n),
+                        css_params = make_data(n),
+                        fit_params = make_data(n),
+                        C_ext = make_data(n))
+  
+  expect_no_error(plot_sensitivity(x))
+})
+
+test_that("with assay", {
+  
+  make_data <- function(n = 5, metric = "GCA.Eff", assay = c("a1", "a2")) {
+    list(stats::setNames(data.frame(rep(assay, each = n),
+                                    rep(1:n, 2),
+                                    runif(2 * n)),
+                         c("assay", "sample", metric)))
+  }
+  
+  x <- GeoTox()
+  n <- 5
+  x$resp <- make_data(n)
+  x$sensitivity <- list(age = make_data(n),
+                        obesity = make_data(n),
+                        css_params = make_data(n),
+                        fit_params = make_data(n),
+                        C_ext = make_data(n))
+  
+  expect_warning(plot_sensitivity(x),
+                 "Multiple assays found, using first assay")
+  expect_no_error(plot_sensitivity(x, assay = "a1"))
+  expect_no_error(plot_sensitivity(x, assay = "a2"))
+})
+
+test_that("data NAs and missing", {
+  
+  make_data <- function(n = 5, metric = "GCA.Eff", p_NA = 0.5) {
+    df <- stats::setNames(data.frame(1:n, runif(n)),
+                          c("sample", metric))
+    df[sample(n, n * p_NA), metric] <- NA
+    list(df)
+  }
+  
+  x <- GeoTox()
+  n <- 5
+  x$resp <- make_data(n)
+  x$sensitivity <- list(age = make_data(n),
+                        obesity = make_data(n),
+                        css_params = make_data(n),
+                        fit_params = make_data(n),
+                        C_ext = make_data(n))
+  
+  expect_warning(plot_sensitivity(x),
+                 "Removed \\d+ NA value")
+  
+  p_NA <- 1
+  x$resp <- make_data(n, p_NA = p_NA)
+  x$sensitivity <- list(age = make_data(n, p_NA = p_NA),
+                        obesity = make_data(n, p_NA = p_NA),
+                        css_params = make_data(n, p_NA = p_NA),
+                        fit_params = make_data(n, p_NA = p_NA),
+                        C_ext = make_data(n, p_NA = p_NA))
+  
+  expect_error(expect_warning(plot_sensitivity(x),
+                              "Removed \\d+ NA value"),
+               "No data to plot")
+})
diff --git a/tests/testthat/test-simulate_population.R b/tests/testthat/test-simulate_population.R
index 29af5d6..67e17ad 100644
--- a/tests/testthat/test-simulate_population.R
+++ b/tests/testthat/test-simulate_population.R
@@ -127,25 +127,6 @@ test_that("clear downstream - obesity", {
   
 })
 
-test_that("default params", {
-  
-  geoTox <- GeoTox()
-  
-  expect_equal(geoTox$par,
-               list(n = 1e3,
-                    IR_params = NULL,
-                    obesity = list(obes_prev  = "OBESITY_CrudePrev",
-                                   obes_sd    = "OBESITY_SD",
-                                   obes_label = "FIPS"),
-                    exposure = list(expos_mean  = "mean",
-                                    expos_sd    = "sd",
-                                    expos_label = "casn"),
-                    internal_dose = list(time    = 1,
-                                         BW      = 1,
-                                         scaling = 1),
-                    resp = list(max_mult = 1.5)))
-})
-
 test_that("update params", {
   
   # New params