Add plot_serocurve() for population-level serodynamic curve visualization

DESCRIPTION

@@ -1,6 +1,6 @@
 Package: serodynamics
 Title: What the Package Does (One Line, Title Case)
-Version: 0.0.0.9055
+Version: 0.0.0.9056
 Authors@R: c(
     person("Peter", "Teunis", , "p.teunis@emory.edu", role = c("aut", "cph"),
            comment = "Author of the method and original code."),

NAMESPACE

@@ -13,6 +13,7 @@ export(plot_jags_dens)
 export(plot_jags_effect)
 export(plot_jags_trace)
 export(plot_predicted_curve)
+export(plot_serocurve)
 export(post_summ)
 export(postprocess_jags_output)
 export(prep_data)

NEWS.md

@@ -1,5 +1,10 @@
 # serodynamics (development version)
 
+* Added `plot_serocurve()` for graphical visualization of population-level
+  serodynamic curves using posterior samples of the `mu.par` hyperparameter
+  (or optionally the "newperson" subject).  Supports 95% credible interval
+  ribbons, stratified curves with colour or faceting, and multiple
+  antigen-isotypes (#74).
 * Clarified Code Style Guidelines in `.github/copilot-instructions.md`:
   the UCD-SeRG Lab Manual takes precedence over the tidyverse style
   guide where they conflict, and functions should end with an explicit

R/Run_Mod.R

@@ -222,7 +222,7 @@ run_mod <- function(data,
   jags_out <- jags_out[, c("Iteration", "Chain", "Parameter", "Iso_type",
                            "Stratification", "Subject", "value")]
   jags_out <- rebuild_sr_model_attributes(jags_out, mod_atts)
-  
+
   # Adding population parameters optionally and priors in as attributes
   if (with_pop_params) {
     jags_out <- jags_out |>

R/nepal_sees_jags_output.R

@@ -27,7 +27,11 @@
 #'  `sees_npl_2`}
 #'  \item{value}{Estimated value of the parameter}
 #'  \item{attributes}{A [list] of `attributes` that summarize the jags inputs, 
-#'  priors, and optional jags_post mcmc object}
+#'  priors, optional jags_post mcmc object, and population-level parameters.
+#'  Includes a `population_params` attribute with posterior samples of the
+#'  population-level parameters (`mu.par`, `prec.par`, `prec.logy`), indexed
+#'  by `Iteration`, `Chain`, `Parameter`, `Iso_type`, `Stratification`, and
+#'  `Population_Parameter`.}
 #' }
 #' @source reference study: <https://doi.org/10.1016/S2666-5247(22)00114-8>
 "nepal_sees_jags_output"

R/plot_serocurve.R

@@ -0,0 +1,319 @@
+#' @title Plot Estimated Serodynamic Curves at the Population Level
+#' @description
+#' Plots the estimated antibody response curve derived from posterior samples
+#' of population-level (`mu.par`) or the predictive distribution from a fitted
+#' [run_mod()] model.  A median curve with an optional 95% credible interval
+#' ribbon is produced for each requested antigen-isotype and stratification
+#' combination.
+#'
+#' @param model An `sr_model` object returned by [run_mod()].
+#' @param antigen_iso A [character] vector of antigen-isotype combinations to
+#'   plot.  Defaults to all antigen-isotypes present in the subject-level
+#'   draws of `model` (`model$Iso_type`); in normal usage these match the
+#'   levels available in `attr(model, "population_params")`.
+#' @param strat A [character] vector of stratification levels to include.
+#'   Defaults to all stratification levels present in the subject-level
+#'   draws of `model` (`model$Stratification`); in normal usage these match
+#'   the levels available in `attr(model, "population_params")`.
+#' @param param_source [character]; which posterior samples to use for the
+#'   curve.  Options:
+#'   - `"population"` (default): uses population-level `mu.par` samples stored
+#'     in `attr(model, "population_params")`. Requires the model to have been
+#'     fitted with `run_mod(..., with_pop_params = TRUE)`.
+#'   - `"newperson"`: uses the predictive distribution for a new individual
+#'     drawn from the population-level prior.
+#' @param show_ci [logical]; if [TRUE] (default), draws a 95% credible
+#'   interval ribbon around the median curve.
+#' @param log_y [logical]; if [TRUE], applies a [log10] transformation to the
+#'   y-axis.  Defaults to [FALSE].
+#' @param log_x [logical]; if [TRUE], applies a pseudo-log10 transformation to
+#'   the x-axis.  Defaults to [FALSE].
+#' @param xlim (Optional) A numeric vector of length 2 giving custom x-axis
+#'   limits.
+#' @param facet_by_antigen_iso [logical]; if [TRUE], facets the plot by
+#'   antigen-isotype.  Defaults to [TRUE] when multiple antigen-isotypes are
+#'   requested.
+#' @param facet_by_strat [logical]; if [TRUE], facets the plot by
+#'   stratification level.  When [FALSE] (default), different stratification
+#'   levels are shown as different colours on the same panel.
+#' @param ncol [integer]; number of columns when faceting.  If [NULL]
+#'   (default), a sensible value is chosen automatically.
+#'
+#' @return A [ggplot2::ggplot] object.
+#' @export
+#'
+#' @example inst/examples/examples-plot_serocurve.R
+plot_serocurve <- function(
+    model,
+    antigen_iso = unique(model$Iso_type),
+    strat = unique(model$Stratification),
+    param_source = "population",
+    show_ci = TRUE,
+    log_y = FALSE,
+    log_x = FALSE,
+    xlim = NULL,
+    facet_by_antigen_iso = length(antigen_iso) > 1,
+    facet_by_strat = FALSE,
+    ncol = NULL) {
+
+  param_source <- match.arg(param_source, c("population", "newperson"))
+
+  antigen_iso_col <- "Iso_type"
+
+  # ---- Retrieve posterior samples of curve parameters --------------------
+  if (param_source == "population") {
+    pop_params <- attr(model, "population_params")
+    if (is.null(pop_params)) {
+      cli::cli_abort(
+        c(
+          "The {.arg model} object does not have a {.field population_params}",
+          " attribute.",
+          "i" = paste0(
+            "Re-fit the model with",
+            " {.code run_mod(..., with_pop_params = TRUE)}."
+          )
+        )
+      )
+    }
+    # The population_params tibble has columns:
+    # Iteration, Chain, Parameter, Iso_type, Stratification,
+    # Population_Parameter, value
+    # Filter to mu.par rows only, then pivot wider and transform from log scale.
+    param_samples <- pop_params |>
+      dplyr::filter(
+        .data$Population_Parameter == "mu.par",
+        .data$Iso_type %in% .env$antigen_iso,
+        .data$Stratification %in% .env$strat
+      ) |>
+      dplyr::select(
+        all_of(
+          c("Chain", "Iteration", "Parameter", "Iso_type", "Stratification",
+            "value")
+        )
+      ) |>
+      tidyr::pivot_wider(
+        names_from = "Parameter",
+        values_from = "value",
+        names_prefix = "log_"
+      ) |>
+      dplyr::mutate(
+        y0    = exp(.data$log_y0),
+        y1    = .data$y0 + exp(.data$log_y1),
+        t1    = exp(.data$log_t1),
+        alpha = exp(.data$log_alpha),
+        shape = exp(.data$log_shape) + 1
+      ) |>
+      dplyr::select(
+        -dplyr::starts_with("log_")
+      ) |>
+      dplyr::mutate(
+        Iso_type = factor(.data$Iso_type),
+        Stratification = factor(.data$Stratification)
+      )
+  } else {
+    # "newperson": predictive distribution for a new individual drawn from
+    # the population-level prior
+    newperson_rows <- model |>
+      dplyr::filter(
+        .data$Subject == "newperson",
+        .data$Iso_type %in% .env$antigen_iso,
+        .data$Stratification %in% .env$strat
+      )
+
+    if (nrow(newperson_rows) == 0) {
+      cli::cli_abort(
+        c(
+          paste0(
+            "No {.val newperson} subject found in {.arg model} for the ",
+            "requested {.arg antigen_iso}/{.arg strat}."
+          ),
+          "i" = paste0(
+            "Ensure the model was fit with a {.val newperson} subject ",
+            "included."
+          )
+        )
+      )
+    }
+
+    param_samples <- newperson_rows |>
+      dplyr::select(
+        all_of(
+          c("Chain", "Iteration", "Parameter", "Iso_type", "Stratification",
+            "value")
+        )
+      ) |>
+      tidyr::pivot_wider(
+        names_from = "Parameter",
+        values_from = "value"
+      ) |>
+      dplyr::mutate(
+        Iso_type = factor(.data$Iso_type),
+        Stratification = factor(.data$Stratification)
+      )
+  }
+
+  # ---- Compute predicted curves over a grid of time points ---------------
+  # Clamp the grid to `xlim` when supplied to avoid unnecessary computation
+  # outside the visible range.
+  if (!is.null(xlim)) {
+    tx <- seq(xlim[1], xlim[2], by = 5)
+  } else {
+    tx <- seq(0, 1200, by = 5)
+  }
+
+  serocourse_all <- param_samples |>
+    dplyr::reframe(
+      t = .env$tx,
+      res = ab(.data$t, .data$y0, .data$y1, .data$t1, .data$alpha,
+               .data$shape),
+      .by = all_of(
+        c("Chain", "Iteration", antigen_iso_col, "Stratification")
+      )
+    )
+
+  # ---- Summarise to median + 95 % CI -------------------------------------
+  curve_summary <- serocourse_all |>
+    dplyr::summarise(
+      .by = all_of(c(antigen_iso_col, "Stratification", "t")),
+      res_med  = stats::quantile(.data$res, probs = 0.50, na.rm = TRUE),
+      res_low  = stats::quantile(.data$res, probs = 0.025, na.rm = TRUE),
+      res_high = stats::quantile(.data$res, probs = 0.975, na.rm = TRUE)
+    )
+
+  # ---- Determine whether to colour by stratification ---------------------
+  n_strat <- length(unique(param_samples$Stratification))
+  multi_strat <- n_strat > 1 && !facet_by_strat
+
+  # ---- Build the ggplot --------------------------------------------------
+  p <- ggplot2::ggplot() +
+    ggplot2::theme_minimal() +
+    ggplot2::labs(x = "Time since onset", y = "Assay result") +
+    ggplot2::theme(legend.position = "bottom")
+
+  if (show_ci) {
+    if (multi_strat) {
+      p <- p +
+        ggplot2::geom_ribbon(
+          data = curve_summary,
+          ggplot2::aes(
+            x = .data$t,
+            ymin = .data$res_low,
+            ymax = .data$res_high,
+            fill = .data$Stratification
+          ),
+          alpha = 0.2,
+          inherit.aes = FALSE
+        )
+    } else {
+      p <- p +
+        ggplot2::geom_ribbon(
+          data = curve_summary,
+          ggplot2::aes(
+            x = .data$t,
+            ymin = .data$res_low,
+            ymax = .data$res_high,
+            fill = "ci"
+          ),
+          alpha = 0.2,
+          inherit.aes = FALSE
+        )
+    }
+  }
+
+  # Median line
+  if (multi_strat) {
+    p <- p +
+      ggplot2::geom_line(
+        data = curve_summary,
+        ggplot2::aes(
+          x = .data$t,
+          y = .data$res_med,
+          colour = .data$Stratification
+        ),
+        linewidth = 1,
+        inherit.aes = FALSE
+      )
+  } else {
+    p <- p +
+      ggplot2::geom_line(
+        data = curve_summary,
+        ggplot2::aes(
+          x = .data$t,
+          y = .data$res_med,
+          colour = "median"
+        ),
+        linewidth = 1,
+        inherit.aes = FALSE
+      )
+  }
+
+  # ---- Legend for single-stratification plots ----------------------------
+  if (!multi_strat) {
+    color_vals <- c("median" = "red")
+    color_labels <- c("median" = "Median")
+
+    p <- p +
+      ggplot2::scale_color_manual(
+        name = "",
+        values = color_vals,
+        labels = color_labels,
+        guide = ggplot2::guide_legend(override.aes = list(shape = NA))
+      )
+
+    if (show_ci) {
+      fill_vals <- c("ci" = "red")
+      fill_labels <- c("ci" = "95% credible interval")
+
+      p <- p +
+        ggplot2::scale_fill_manual(
+          name = "",
+          values = fill_vals,
+          labels = fill_labels,
+          guide = ggplot2::guide_legend(override.aes = list(color = NA))
+        )
+    }
+  }
+
+  # ---- Faceting ----------------------------------------------------------
+  facet_vars <- character(0)
+  if (facet_by_antigen_iso) facet_vars <- c(facet_vars, antigen_iso_col)
+  if (facet_by_strat)       facet_vars <- c(facet_vars, "Stratification")
+
+  if (length(facet_vars) > 0) {
+    facet_formula <- stats::as.formula(
+      paste("~", paste(facet_vars, collapse = " + "))
+    )
+    if (is.null(ncol)) {
+      n_panels <- length(unique(interaction(
+        curve_summary[[facet_vars[1]]],
+        if (length(facet_vars) > 1) curve_summary[[facet_vars[2]]] else NULL
+      )))
+      ncol <- if (n_panels == 1) {
+        1
+      } else if (n_panels <= 4) {
+        2
+      } else {
+        NULL
+      }
+    }
+    p <- p + ggplot2::facet_wrap(facet_formula, ncol = ncol)
+  }
+
+  # ---- Log scales --------------------------------------------------------
+  if (log_y) {
+    p <- p + ggplot2::scale_y_log10()
+  }
+  if (log_x) {
+    p <- p +
+      ggplot2::scale_x_continuous(
+        trans = scales::pseudo_log_trans(sigma = 1, base = 10)
+      )
+  }
+
+  # ---- Custom x-axis limits ----------------------------------------------
+  if (!is.null(xlim)) {
+    p <- p + ggplot2::coord_cartesian(xlim = xlim)
+  }
+
+  return(p)
+}

R/use_att_names.R

@@ -2,7 +2,7 @@
 #' @description
 #'  `use_att_names` takes prepared longitudinal data for antibody kinetic
 #'  modeling and names columns using attribute values to allow merging
-#'  with a modeled [run_mod()] output [tibble::tbl_df]. The column names include
+#'  with a modeled [run_mod] output [tibble::tbl_df]. The column names include
 #'  `Subject`, `Iso_type`, `t`, and `result`. 
 #' @param data A [data.frame] raw longitudinal data that has been
 #' prepared for antibody kinetic modeling using [as_case_data()].

inst/WORDLIST

@@ -70,4 +70,6 @@ unstratified
 vh
 wishdf
 SeRG
+Serodynamic
 UCD
+serodynamic