Merge commit 'd7d1eba7c81d4e4a299d6b8642acd8b419f67add'

Author: d-morrison

.github/workflows/check-spelling.yaml

@@ -14,10 +14,18 @@ jobs:
     runs-on: ubuntu-latest
     name: Spellcheck
     container:
-      image: rocker/tidyverse:4.1.2
+      image: rocker/tidyverse:latest
     steps:
       - name: Checkout repo
         uses: actions/checkout@v3
 
-      - name: Run Spelling Check test
+      - uses: r-lib/actions/setup-r@v2
+        with:
+          use-public-rspm: true
+
+      - uses: r-lib/actions/setup-r-dependencies@v2
+        with:
+          packages: any::spelling
+
+      - name: Run Spelling Check Test
         uses: insightsengineering/r-spellcheck-action@v3.0.2

DESCRIPTION

@@ -1,12 +1,13 @@
 Type: Package
 Package: serocalculator
 Title: Estimating Infection Rates from Serological Data
-Version: 1.3.0.9048
+Version: 1.3.0.9055
 Authors@R: c(
     person("Peter", "Teunis", , "p.teunis@emory.edu", role = c("aut", "cph"),
            comment = "Author of the method and original code."),
     person("Kristina", "Lai", , "kwlai@ucdavis.edu", role = c("aut", "cre")),
     person("Chris", "Orwa", role = "aut"),
+    person("Kwan Ho", "Lee", , "ksjlee@ucdavis.edu", role = "ctb"),
     person("Kristen", "Aiemjoy", , "kaiemjoy@ucdavis.edu", role = "aut"),
     person("Douglas Ezra", "Morrison", , "demorrison@ucdavis.edu", role = "aut")
   )

NAMESPACE

@@ -5,6 +5,7 @@ S3method(autoplot,curve_params)
 S3method(autoplot,pop_data)
 S3method(autoplot,seroincidence)
 S3method(autoplot,seroincidence.by)
+S3method(autoplot,sim_results)
 S3method(autoplot,summary.seroincidence.by)
 S3method(print,seroincidence)
 S3method(print,seroincidence.by)
@@ -14,6 +15,7 @@ S3method(strata,default)
 S3method(summary,pop_data)
 S3method(summary,seroincidence)
 S3method(summary,seroincidence.by)
+export(analyze_sims)
 export(as_curve_params)
 export(as_noise_params)
 export(as_pop_data)
@@ -120,6 +122,7 @@ importFrom(tidyselect,contains)
 importFrom(tidyselect,ends_with)
 importFrom(utils,capture.output)
 importFrom(utils,download.file)
+importFrom(utils,head)
 importFrom(utils,tail)
 importFrom(utils,unzip)
 useDynLib(serocalculator, .registration = TRUE)

NEWS.md

@@ -2,9 +2,14 @@
 
 ## New features
 
+* Made `graph.curve.params()` the default sub-method for `autoplot.curve_params()` (#450)
+* Added `log_x` and `log_y` options to `graph.curve.params()` sub-method for 
+`autoplot.curve_params()` (#453)
+* Extended `sim_pop_data_multi()` to loop over multiple sample sizes (#444)
+* Added new functions `analyze_sims()` and `autoplot.sim_results()` (#444)
 * Rename `estimate_scr()` to `est_seroincidence_by()` (#439)
 * Rename `estimate_scr()` to `est_seroincidence()` (#432)
-* Rename argument `curve_params` to `sr_params` (#424)
+* Rename argument `curve_params` to `sr_params` for estimation functions (#424)
 * added documentation for `count_strata()` (#431)
 * Rename  `as_curve_params()` to `as_sr_params()` (#421)
 * Rename `load_curve_params()` to `load_sr_params()` (#421)

R/ab1.R

@@ -0,0 +1,16 @@
+# uses r > 1 scale for shape
+ab1 <- function(t, y0, y1, t1, alpha, shape) {
+  beta <- bt(y0, y1, t1)
+
+  yt <- 0
+
+  if (t <= t1) {
+    yt <- y0 * exp(beta * t)
+  }
+
+  if (t > t1) {
+    yt <- (y1^(1 - shape) - (1 - shape) * alpha * (t - t1))^(1 / (1 - shape))
+  }
+
+  return(yt)
+}

R/analyze_sims.R

@@ -0,0 +1,69 @@
+#' Analyze simulation results
+#'
+#' @param data a [tibble::tbl_df] with columns:
+#' * `lambda.sim`,
+#' * `incidence.rate`,
+#' * `SE`,
+#' * `CI.lwr`,
+#' * `CI.upr`
+#' for example, as produced by [summary.seroincidence.by()] with
+#' `lambda.sim` as a stratifying variable
+#'
+#' @returns a `sim_results` object (extends [tibble::tbl_df])
+#' @export
+#'
+#' @example inst/examples/exm-analyze_sims.R
+#'
+analyze_sims <- function(
+    data) {
+
+  to_return <-
+    data |>
+    split(
+      f = ~  sample_size + lambda.sim
+    ) |>
+    lapply(FUN = analyze_sims_one_stratum) |>
+    bind_rows()
+
+  class(to_return) <- union("sim_results", class(to_return))
+
+  return(to_return)
+}
+
+analyze_sims_one_stratum <- function(
+    data,
+    true_lambda = data$lambda.sim,
+    sample_size = data$sample_size) {
+
+  # Filter out rows where CI.lwr or CI.upr is Inf or NaN
+  data <- data |>
+    filter(is.finite(.data$CI.lwr) & is.finite(.data$CI.upr))
+
+  # Compute Bias
+  bias <- mean(data$incidence.rate - true_lambda, na.rm = TRUE)
+
+  # Standard Error (Mean of reported standard errors)
+  standard_error <- mean(data$SE, na.rm = TRUE)
+
+  # RMSE (Root Mean Square Error)
+  rmse <- mean((data$incidence.rate - true_lambda)^2, na.rm = TRUE) |> sqrt()
+
+  # Confidence Interval Width (Mean of Upper - Lower bounds, without Inf values)
+  ci_width <- mean(data$CI.upr - data$CI.lwr, na.rm = TRUE)
+
+  coverage_prop <-
+    mean(data$CI.lwr <= true_lambda & data$CI.upr >= true_lambda, na.rm = TRUE)
+
+  to_return <- tibble(
+    lambda.sim = mean(true_lambda),
+    sample_size = mean(sample_size),
+    Bias = bias,
+    Mean_Est_SE = standard_error,
+    Empirical_SE = stats::sd(data$incidence.rate, na.rm = TRUE),
+    RMSE = rmse,
+    Mean_CI_Width = ci_width,
+    CI_Coverage = coverage_prop
+  )
+  # Return computed statistics as a list
+  return(to_return)
+}

R/as_sr_params.R

@@ -67,6 +67,22 @@ as_sr_params <- function(data, antigen_isos = NULL) {
     ))
   }
 
+  # if `object` lacks an `iter` column, add it:
+  if (!is.element("iter", names(curve_data))) {
+    cli::cli_warn(
+      c(
+        "`data` is missing `iter` column",
+        "It will be inferred from row ordering."
+      )
+    )
+    curve_data <-
+      curve_data |>
+      mutate(
+        .by = any_of(c("antigen_iso", "chain")),
+        iter = row_number()
+      )
+  }
+
   # assign antigen attribute
   attr(curve_data, "antigen_isos") <- antigen_isos
 

R/autoplot.curve_params.R

@@ -1,9 +1,18 @@
-#' graph antibody decay curves by antigen isotype
+#' Graph antibody decay curves by antigen isotype
+#' @param object
+#' a `curve_params` object (constructed using [as_sr_params()]), which is
+#' a [data.frame()] containing MCMC samples of antibody decay curve parameters
+#' @param method a [character] string indicating whether to use
+#'  - [graph.curve.params()] (default) or
+#'  - [graph_seroresponse_model_1()] (previous default)
+#' as the graphing method.
 #'
-#' @inheritParams graph_seroresponse_model_1
-#' @inheritDotParams graph_seroresponse_model_1
-#' @param antigen_isos antigen isotypes to analyze
-#' (can subset `curve_params`)
+#' @param ... additional arguments passed to the sub-function
+#' indicated by the `method` argument.
+#' @details
+#' Currently, the backend for this method is [graph.curve.params()].
+#' Previously, the backend for this method was [graph_seroresponse_model_1()].
+#' That function is still available if preferred.
 #' @return a [ggplot2::ggplot()] object
 #' @export
 #' @examples
@@ -23,14 +32,12 @@
 #' }
 autoplot.curve_params <- function(
     object,
-    antigen_isos = unique(object$antigen_iso),
+    method = c("graph.curve.params", "graph_seroresponse_model_1"),
     ...) {
 
   # spaghettified in order to swap out implementations with minimal
   # disruption to API
-  object |>
-    graph_seroresponse_model_1(
-      antigen_isos = antigen_isos,
-      ...
-    )
+  method <- match.arg(method)
+  cur_function <- match.fun(method)
+  object |> cur_function(...)
 }

R/autoplot.seroincidence.by.R

@@ -24,7 +24,7 @@
 #' est2 <- est_seroincidence_by(
 #'   strata = c("catchment"),
 #'   pop_data = xs_data,
-#'   curve_params = curve,
+#'   sr_params = curve,
 #'   curve_strata_varnames= NULL,
 #'   noise_strata_varnames = NULL,
 #'   noise_params = noise,

R/autoplot.sim_results.R

@@ -0,0 +1,27 @@
+#' Plot simulation results
+#' `autoplot()` method for `sim_results` objects
+#'
+#' @param object a `sim_results` object (from [analyze_sims()])
+#' @param statistic which column of `object` should be the y-axis?
+#' @param ... unused
+#' @returns a [ggplot2::ggplot]
+#' @export
+#'
+#' @example inst/examples/exm-autoplot.sim_results.R
+autoplot.sim_results <- function(
+    object,
+    statistic = "Empirical_SE",
+    ...) {
+  object |>
+    dplyr::mutate(lambda.sim = factor(.data$lambda.sim)) |>
+    ggplot2::ggplot() +
+    ggplot2::aes(
+      x = .data$sample_size,
+      group = .data$lambda.sim,
+      col = .data$lambda.sim,
+      y = .data[[statistic]]
+    ) +
+    ggplot2::geom_point() +
+    ggplot2::geom_line() +
+    ggplot2::theme(legend.position = "bottom")
+}