Add quantiles parameter to graph.curve.params()

DESCRIPTION

@@ -1,7 +1,7 @@
 Type: Package
 Package: serocalculator
 Title: Estimating Infection Rates from Serological Data
-Version: 1.3.0.9047
+Version: 1.3.0.9048
 Authors@R: c(
     person("Peter", "Teunis", , "p.teunis@emory.edu", role = c("aut", "cph"),
            comment = "Author of the method and original code."),

NEWS.md

@@ -21,6 +21,7 @@
 * Improved y-limit calculation in `graph.curve.params()` (#368)
 * Added option for `graph.curve.params()` to show all curves (#368)
 * Added color-coding for `graph.curve.params()` (#383)
+* Added `quantiles` parameter to `graph.curve.params()` and corresponding test in `test-graph.curve.params.R` (#434)
 * Removed `warn.missing.strata()` from API (#366)
 
 * Added more details about contributing PRs in `Contributing.md` (#280)

R/graph.curve.params.R

@@ -7,30 +7,24 @@
 #' @param antigen_isos antigen isotypes
 #' @param alpha_samples `alpha` parameter passed to [ggplot2::geom_line]
 #' (has no effect if `show_all_curves = FALSE`)
-#' @param show_quantiles whether to show point-wise (over time) quantiles
+#' @param quantiles Optional [numeric] [vector] of quantiles to plot
+#' (e.g., 10%, 50%, and 90% = `c(0.1, 0.5, 0.9)`). If `NULL`, no quantile
+#' lines are shown.
 #'
-#' @returns a [ggplot2::ggplot()] object
+#' @returns a [ggplot2::ggplot()] object showing the antibody dynamic
+#' kinetics of selected antigen/isotype combinations, with optional posterior
+#' distribution quantile curves.
 #' @export
 #'
-#' @examples
-#' curve <-
-#'   typhoid_curves_nostrat_100 |>
-#'   dplyr::filter(antigen_iso %in% c("HlyE_IgA", "HlyE_IgG"))
+#' @examples inst/examples/exm-graph.curve.params.R
 #'
-#' plot1 <- graph.curve.params(curve)
-#'
-#' print(plot1)
-#'
-#' plot2 <- graph.curve.params(curve, show_all_curves = TRUE)
-#' show(plot2)
-#'
-graph.curve.params <- function( # nolint: object_name_linter
+graph.curve.params <- function(# nolint: object_name_linter
   curve_params,
   antigen_isos = unique(curve_params$antigen_iso),
   verbose = FALSE,
-  show_quantiles = TRUE,
   show_all_curves = FALSE,
-  alpha_samples = 0.3
+  alpha_samples = 0.3,
+  quantiles = c(0.1, 0.5, 0.9)  # numeric, flexible
 ) {
   if (verbose) {
     message(
@@ -39,34 +33,12 @@
     )
   }
 
-  curve_params <- curve_params |>
+  object <- object |>
     dplyr::filter(.data$antigen_iso %in% antigen_isos)
 
   tx2 <- 10^seq(-1, 3.1, 0.025)
 
-  bt <- function(y0, y1, t1) {
-    log(y1 / y0) / t1
-  }
-
-  # uses r > 1 scale for shape
-  ab <- 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)
-  }
-
-
-  d <- curve_params
+  d <- object
 
   dT <- # nolint: object_linter
     data.frame(t = tx2) |>
@@ -83,149 +55,125 @@
       )
     )
 
-
   serocourse_all <-
     cbind(d, dT) |>
     tidyr::pivot_longer(
       cols = dplyr::starts_with("time"),
       values_to = "t"
     ) |>
-    select(-"name") |>
-    rowwise() |>
-    mutate(res = ab(
-      .data$t,
-      .data$y0,
-      .data$y1,
-      .data$t1,
-      .data$alpha,
-      .data$r
-    )) |>
-    ungroup()
-
-  if (verbose) message("starting to compute quantiles")
-  serocourse_sum <- serocourse_all |>
-    summarise(
-      .by = c("antigen_iso", "t"),
-      res.med = quantile(.data$res, 0.5),
-      res.low = quantile(.data$res, 0.025),
-      res.high = quantile(.data$res, 0.975),
-      res.p75 = quantile(.data$res, 0.75),
-      res.p25 = quantile(.data$res, 0.25),
-      res.p10 = quantile(.data$res, 0.10),
-      res.p90 = quantile(.data$res, 0.90)
+    dplyr::select(-c("name")) |>
+    dplyr::rowwise() |>
+    dplyr::mutate(
+      res = ab1(
+        .data$t,
+        .data$y0,
+        .data$y1,
+        .data$t1,
+        .data$alpha,
+        .data$r
+      )
     ) |>
-    pivot_longer(
-      names_to = "quantile",
-      cols = c(
-        "res.med",
-        "res.low",
-        "res.high",
-        "res.p25",
-        "res.p75",
-        "res.p10",
-        "res.p90"
-      ),
-      names_prefix = "res.",
-      values_to = "res"
-    )
-
+    dplyr::ungroup()
+
+  if (!is.null(quantiles)) {
+    serocourse_sum <- serocourse_all |>
+      dplyr::group_by(.data$antigen_iso, .data$t) |>
+      dplyr::summarise(
+        res_vals = list(.data$res),
+        .groups = "drop"
+      ) |>
+      dplyr::mutate(
+        quantiles_df = purrr::map(
+          .data$res_vals,
+          ~ tibble::tibble(
+            quantile = quantiles,
+            res = stats::quantile(.x, probs = quantiles, na.rm = TRUE)
+          )
+        )
+      ) |>
+      tidyr::unnest(c("quantiles_df"))
+  }
 
   range <-
-    serocourse_sum |>
-    dplyr::filter(.data$quantile %in% c("med", "p10", "p90")) |>
+    serocourse_all |>
     dplyr::summarize(
-      min = min(.data$res, 0.9),
-      max = max(.data$res, 2000)
+      min = min(.data$res, na.rm = TRUE),
+      max = max(.data$res, na.rm = TRUE)
     )
 
-
-  plot1 <-
-    serocourse_sum |>
-    ggplot2::ggplot() +
-    ggplot2::aes(x = .data$t,
-                 y = .data$res) +
-    ggplot2::facet_wrap(
-      ~ .data$antigen_iso,
-      ncol = 2
-    ) +
+  plot1 <- ggplot2::ggplot() +
+    ggplot2::aes(x = .data$t, y = .data$res) +
+    ggplot2::facet_wrap(~ antigen_iso, ncol = 2) +
     ggplot2::theme_minimal() +
     ggplot2::theme(axis.line = ggplot2::element_line()) +
-    ggplot2::labs(x = "Days since fever onset",
-                  y = "ELISA units",
-                  col = if_else(show_all_curves, "MCMC chain", "")) +
+    ggplot2::labs(
+      x   = "Days since fever onset",
+      y   = "ELISA units",
+      col = if (show_all_curves) "MCMC chain" else NULL
+    ) +
     ggplot2::theme(legend.position = "bottom")
 
-  if (show_all_curves) {
-
-    range <-
-      serocourse_all |>
-      dplyr::summarize(
-        min = min(.data$res, 0.9),
-        max = max(.data$res, 2000)
-      )
 
-    group_vars <-
-      c("iter", "chain") |>
+  if (show_all_curves) {
+    group_vars <- c("iter", "chain") |>
       intersect(names(serocourse_all))
 
     if (length(group_vars) > 1) {
       serocourse_all <-
         serocourse_all |>
-        mutate(
-          iter = interaction(across(all_of(group_vars)))
+        dplyr::mutate(
+          iter = interaction(
+            dplyr::across(dplyr::all_of(group_vars))
+          )
         )
+
       plot1 <-
         plot1 +
-        geom_line(
+        ggplot2::geom_line(
           data = serocourse_all,
           alpha = alpha_samples,
           aes(
-            color = .data$chain |> factor(),
+            color = factor(.data$chain),
             group = .data$iter
           )
-        ) +
-        ggplot2::expand_limits(y = range)
+        )
     } else {
-
       plot1 <-
         plot1 +
-        geom_line(data = serocourse_all,
-                  alpha = alpha_samples,
-                  aes(group = .data$iter)) +
-        ggplot2::expand_limits(y = range)
-
+        ggplot2::geom_line(
+          data = serocourse_all,
+          alpha = alpha_samples,
+          aes(group = .data$iter)
+        )
     }
-
+    plot1 <-
+      plot1 + ggplot2::expand_limits(y = unlist(range))
   }
 
-  plot1 <-
-    plot1 +
+  plot1 <- plot1 +
     ggplot2::scale_y_log10(
       limits = unlist(range),
       labels = scales::label_comma(),
       minor_breaks = NULL
     )
 
-  if (show_quantiles) {
-    plot1 <-
-      plot1 +
-      ggplot2::geom_line(
-        ggplot2::aes(col = "median"),
-        data = serocourse_sum |> filter(.data$quantile == "med"),
-        linewidth = 1
-      ) +
+  if (!is.null(quantiles)) {
+    plot1 <- plot1 +
       ggplot2::geom_line(
-        ggplot2::aes(col = "10% quantile"),
-        data = serocourse_sum |> filter(quantile == "p10"),
-        linewidth = .5
-      ) +
-      ggplot2::geom_line(
-        ggplot2::aes(col = "90% quantile"),
-        data = serocourse_sum |> filter(quantile == "p90"),
-        linewidth = .5
+        data = serocourse_sum,
+        aes(
+          color = paste0(.data$quantile * 100, "% quantile"),
+          group = .data$quantile
+        ),
+        linewidth = 0.75
       )
+
+    if (show_all_curves) {
+      plot1 <- plot1 +
+        ggplot2::labs(col = "MCMC chain")
+    }
   }
 
-  return(plot1)
 
+  return(plot1)
 }

inst/examples/exm-graph.curve.params.R

@@ -0,0 +1,15 @@
+# Load example dataset
+curve <- typhoid_curves_nostrat_100 |>
+  dplyr::filter(antigen_iso %in% c("HlyE_IgA", "HlyE_IgG"))
+
+# Plot without showing all curves
+plot1 <- graph.curve.params(curve)
+print(plot1)
+
+# Plot with additional quantiles and show all curves
+plot2 <- graph.curve.params(
+  curve,
+  show_all_curves = TRUE,
+  quantiles = c(0.1, 0.5, 0.9)
+)
+print(plot2)