Merge pull request #284 from UCD-SERG/new-methodology

New methodology

.Rbuildignore

@@ -23,5 +23,9 @@ allpopsamples_hlye.csv$
 ^CRAN-SUBMISSION$
 ^README\.qmd$
 ^codecov\.yml$
+^_quarto\.yml$
 ^\.lintr$
+^vignettes$
 ^man/df_to_array\.Rd$
+^vignettes/\.quarto$
+^vignettes/methodology\.qmd$

.gitignore

@@ -15,3 +15,5 @@ serocalculator*.tar.gz
 serocalculator*.tgz
 README.html
 README_files
+
+/.quarto/

DESCRIPTION

@@ -1,7 +1,7 @@
 Type: Package
 Package: serocalculator
 Title: Estimating Infection Rates from Serological Data
-Version: 1.2.0.9011
+Version: 1.2.0.9012
 Authors@R: c(
     person("Peter", "Teunis", , "p.teunis@emory.edu", role = c("aut", "cph"),
            comment = "Author of the method and original code."),
@@ -49,11 +49,13 @@ Suggests:
     pak,
     parallel,
     readr,
+    quarto,
     rmarkdown,
     spelling,
     ssdtools (>= 1.0.6.9016),
     testthat (>= 3.0.0),
-    tidyverse
+    tidyverse,
+    qrcode
 LinkingTo: 
     Rcpp
 Config/testthat/edition: 3
@@ -67,3 +69,5 @@ Remotes:
     moodymudskipper/devtag
 Roxygen: list(markdown = TRUE, roclets = c("collate", "rd", "namespace", "devtag::dev_roclet"))
 RoxygenNote: 7.3.2
+VignetteBuilder: 
+    quarto

NAMESPACE

@@ -16,6 +16,7 @@ S3method(summary,seroincidence)
 S3method(summary,seroincidence.by)
 export(as_curve_params)
 export(as_pop_data)
+export(autoplot)
 export(check_pop_data)
 export(est.incidence)
 export(est.incidence.by)
@@ -25,7 +26,7 @@ export(fdev)
 export(getAdditionalData)
 export(get_additional_data)
 export(graph.curve.params)
-export(graph.loglik)
+export(graph_loglik)
 export(llik)
 export(load_curve_params)
 export(load_noise_params)

NEWS.md

@@ -6,6 +6,8 @@
 
 * Updated enteric fever example article with upgraded code and visualizations (#290)
 
+* Added `Methodology` vignette (#284)
+
 * Added template for reporting Issues 
 (from `usethis::use_tidy_issue_template()`) (#270)
 
@@ -30,6 +32,7 @@
 * Added `.gitattributes` file (<https://git-scm.com/docs/gitattributes>)
 copied from <https://github.com/tidyverse/ggplot2>
 
+* Added QR code to `README.qmd`
 * Added additional automated checks through 
 [GitHub actions](https://docs.github.com/en/actions), 
 including:

R/est.incidence.R

@@ -119,7 +119,7 @@ est.incidence <- function(
 
   if (build_graph) {
     if (verbose) message("building likelihood graph")
-    graph <- graph.loglik(
+    graph <- graph_loglik(
       highlight_points = lambda_start,
       highlight_point_names = "lambda_start",
       pop_data = pop_data,

R/graph.decay.curves.R

@@ -4,109 +4,88 @@
 #' @param verbose verbose output
 #' @param xlim range of x values to graph
 #' @param n_curves how many curves to plot (see details).
-#' @param n_points Number of points to interpolate along the x axis (passed to [ggplot2::geom_function()])
-#' @param rows_to_graph which rows of `curve_params` to plot (overrides `n_curves`).
-#' @param alpha (passed to [ggplot2::geom_function()]) how transparent the curves should be:
+#' @param n_points Number of points to interpolate along the x axis
+#' (passed to [ggplot2::geom_function()])
+#' @param rows_to_graph which rows of `curve_params` to plot
+#' (overrides `n_curves`).
+#' @param alpha (passed to [ggplot2::geom_function()])
+#' how transparent the curves should be:
 #' * 0 = fully transparent (invisible)
 #' * 1 = fully opaque
-#' @param log_x should the x-axis be on a logarithmic scale (`TRUE`) or linear scale (`FALSE`, default)?
-#' @param log_y should the Y-axis be on a logarithmic scale (default, `TRUE`) or linear scale (`FALSE`)?
+#' @param log_x should the x-axis be on a logarithmic scale (`TRUE`)
+#' or linear scale (`FALSE`, default)?
+#' @param log_y should the Y-axis be on a logarithmic scale
+#' (default, `TRUE`) or linear scale (`FALSE`)?
 #' @inheritParams ggplot2::geom_function
 #' @inheritDotParams ggplot2::geom_function
 #' @returns a [ggplot2::ggplot()] object
 #' @details
 #' ## `n_curves` and `rows_to_graph`
-#' In most cases, `curve_params` will contain too many rows of MCMC samples for all of these samples to be plotted at once.
-#' * Setting the  `n_curves` argument to a value smaller than the number of rows in `curve_params` will cause this function to select the first `n_curves` rows to graph.
-#' * Setting `n_curves` larger than the number of rows in ` will result all curves being plotted.
-#' * If the user directly specifies the `rows_to_graph` argument, then `n_curves` has no effect.
+#' In most cases, `curve_params` will contain too many rows of MCMC
+#' samples for all of these samples to be plotted at once.
+#' * Setting the  `n_curves` argument to a value smaller than the
+#' number of rows in `curve_params` will cause this function to select
+#' the first `n_curves` rows to graph.
+#' * Setting `n_curves` larger than the number of rows in ` will
+#' result all curves being plotted.
+#' * If the user directly specifies the `rows_to_graph` argument,
+#' then `n_curves` has no effect.
 #' @examples
 #' library(dplyr) # loads the `%>%` operator and `dplyr::filter()`
 #'
 #' load_curve_params("https://osf.io/download/rtw5k/") %>%
 #'   dplyr::filter(antigen_iso == "HlyE_IgG") %>%
 #'   serocalculator:::plot_curve_params_one_ab()
 #'
-plot_curve_params_one_ab = function(
+plot_curve_params_one_ab <- function(
     object,
     verbose = FALSE,
     alpha = .4,
     n_curves = 100,
     n_points = 1000,
     log_x = FALSE,
     log_y = TRUE,
-    rows_to_graph = 1:min(n_curves, nrow(object)),
-    xlim = c(10^-1, 10^3.1),
+    rows_to_graph = seq_len(min(n_curves, nrow(object))),
+    xlim = c(10 ^ -1, 10 ^ 3.1),
     ...) {
   plot1 <-
     ggplot2::ggplot() +
     # ggplot2::scale_x_log10() +
     ggplot2::theme_linedraw() +
-    ggplot2::theme(
-      axis.line = ggplot2::element_line()
-    ) +
-    ggplot2::labs(
-      x = "Days since fever onset",
-      y = "Antibody Concentration"
-    ) +
-    ggplot2::ggtitle("Decay Curve") +
-    ggplot2::theme(
-      plot.title =
-        ggplot2::element_text(
-          size = 20,
-          face = "bold"
-        )
-    )
+    ggplot2::theme(axis.line = ggplot2::element_line()) +
+    ggplot2::labs(x = "Days since fever onset", y = "Antibody concentration") +
+    ggplot2::ggtitle("Antibody Response Curve") +
+    ggplot2::theme(plot.title =
+                   ggplot2::element_text(size = 20, face = "bold"))
 
   if (log_y) {
-    plot1 <- plot1 +
-      ggplot2::scale_y_log10(
-        # limits = c(0.9, 2000),
-        labels = scales::label_comma(),
-        # breaks = c(0.1, 1, 10, 100, 1000),
-        minor_breaks = NULL
-      )
+    plot1 <-
+      plot1 +
+      ggplot2::scale_y_log10(labels = scales::label_comma(),
+                             minor_breaks = NULL)
   }
 
   layer_function <- function(cur_row) {
     cur_params <- object[cur_row, ]
     ggplot2::geom_function(
       alpha = alpha,
-      # aes(color = cur_row),
       fun = ab0,
       args = list(curve_params = cur_params),
       n = n_points
     )
   }
 
   layers <-
-    lapply(
-      X = rows_to_graph,
-      FUN = layer_function
-    )
+    lapply(X = rows_to_graph, FUN = layer_function)
 
   plot1 <- plot1 + layers
 
   if (log_x) {
     plot1 <- plot1 +
-      ggplot2::scale_x_log10(
-        limits = xlim,
-        labels = scales::label_comma()
-      )
+      ggplot2::scale_x_log10(limits = xlim, labels = scales::label_comma())
   } else {
     plot1 <- plot1 + ggplot2::xlim(xlim)
   }
 
   return(plot1)
 }
-
-
-# ggplot() +
-#   geom_line(data = serocourse.all, aes(x= t, y = res, group = iter)) +
-#   facet_wrap(~antigen_iso, ncol=2) +
-#   scale_y_log10(limits = c(0.9, 2000), breaks = c(1, 10, 100, 1000), minor_breaks = NULL) +
-#   theme_minimal()  +
-#   theme(axis.line=element_line()) +
-#   labs(x="Days since fever onset", y="ELISA units")
-
-# mcmc %>% ungroup() %>% slice_head(by = antigen_iso, n = 10) %>% droplevels() %>% plot_curve_params_one_ab(alpha  = .4) %>% print()

R/graph.loglik.R

@@ -3,8 +3,10 @@
 #' @param highlight_points a possible highlighted value
 #' @param x sequence of lambda values to graph
 #' @param highlight_point_names labels for highlighted points
-#' @param log_x should the x-axis be on a logarithmic scale (`TRUE`) or linear scale (`FALSE`, default)?
-#' @param previous_plot if not NULL, the current data is added to the existing graph
+#' @param log_x should the x-axis be on a logarithmic scale (`TRUE`)
+#' or linear scale (`FALSE`, default)?
+#' @param previous_plot if not NULL, the current data is added to the
+#' existing graph
 #' @param curve_label if not NULL, add a label for the curve
 #' @inheritParams log_likelihood
 #' @inheritDotParams log_likelihood -lambda
@@ -33,30 +35,34 @@
 #'   y.high = c(5e6, 5e6))                      # High cutoff (y.high)
 #'
 #' # Graph the log likelihood
-#' lik_HlyE_IgA <- graph.loglik(
-#'   pop_data = xs_data,
-#'   curve_params = dmcmc,
-#'   noise_params = cond,
-#'   antigen_isos = "HlyE_IgA",
-#'   log_x = TRUE
+#' lik_HlyE_IgA <- # nolint: object_name_linter
+#'   graph_loglik(
+#'     pop_data = xs_data,
+#'     curve_params = dmcmc,
+#'     noise_params = cond,
+#'     antigen_isos = "HlyE_IgA",
+#'     log_x = TRUE
 #' )
 #'
-#' lik_HlyE_IgA
+#' lik_HlyE_IgA # nolint: object_name_linter
 #'
-graph.loglik = function(
+graph_loglik <- function(
     pop_data,
     curve_params,
     noise_params,
-    antigen_isos,
+    antigen_isos = pop_data %>% get_biomarker_levels(),
     x = 10^seq(-3, 0, by = .1),
     highlight_points = NULL,
     highlight_point_names = "highlight_points",
     log_x = FALSE,
     previous_plot = NULL,
     curve_label = paste(antigen_isos, collapse = " + "),
     ...) {
-  if (!is.list(curve_params) &&
-    !is.element("d", names(curve_params))) {
+
+  needs_rescale <-
+    !is.list(curve_params) &&
+    !is.element("d", names(curve_params))
+  if (needs_rescale) {
     curve_params <-
       curve_params %>%
       dplyr::mutate(

R/reexports.R

@@ -0,0 +1,3 @@
+#' @export
+#' @importFrom ggplot2 autoplot
+ggplot2::autoplot