Simple sim

.github/workflows/R-CMD-check.yaml

@@ -29,15 +29,27 @@ jobs:
 
     steps:
       - uses: actions/checkout@v4
-
-      - uses: r-lib/actions/setup-pandoc@v2
-
       - uses: r-lib/actions/setup-r@v2
         with:
           r-version: ${{ matrix.config.r }}
           http-user-agent: ${{ matrix.config.http-user-agent }}
           use-public-rspm: true
 
+      - uses: r-lib/actions/setup-pandoc@v2
+
+      # from https://github.com/r-quantities/workflows/blob/main/.github/workflows/build.yml
+      - name: Install macOS dependencies
+        if: ${{ runner.os == 'macOS' }}
+        run: |
+          brew install automake udunits
+          echo "UDUNITS2_INCLUDE=/opt/homebrew/Cellar/udunits/2.2.28/include" >> $GITHUB_ENV
+          echo "UDUNITS2_LIBS=/opt/homebrew/Cellar/udunits/2.2.28/lib" >> $GITHUB_ENV
+
+      # c.f. https://github.com/orgs/community/discussions/26947
+      - name: Install Windows dependencies
+        if: ${{ runner.os == 'Windows' }}
+        run: pacman -Syu --noconfirm mingw-w64-ucrt-x86_64-udunits
+
       - uses: r-lib/actions/setup-r-dependencies@v2
         with:
           extra-packages: any::rcmdcheck

DESCRIPTION

@@ -1,7 +1,7 @@
 Type: Package
 Package: serocalculator
 Title: Estimating Infection Rates from Serological Data
-Version: 1.3.0.9029
+Version: 1.3.0.9030
 Authors@R: c(
     person("Peter", "Teunis", , "p.teunis@emory.edu", role = c("aut", "cph"),
            comment = "Author of the method and original code."),
@@ -58,7 +58,9 @@ Suggests:
     ssdtools (>= 1.0.6.9016),
     testthat (>= 3.0.0),
     tidyverse,
+    units (>= 0.8.6),
     qrcode,
+    quantities,
     vdiffr,
     withr
 LinkingTo: 
@@ -67,10 +69,13 @@ Config/testthat/edition: 3
 Config/Needs/roxygen2: roxygen2, moodymudskipper/devtag
 Config/Needs/lint: r-lib/lintr
 Config/Needs/website: quarto, lewinfox/foodwebr
-Config/Needs/check: readr
+Config/Needs/check: readr, r-quantities/quantities, r-quantities/units
 Encoding: UTF-8
 Language: en-US
 LazyData: true
 NeedsCompilation: no
 Roxygen: list(markdown = TRUE, roclets = c("collate", "rd", "namespace", "devtag::dev_roclet"))
 RoxygenNote: 7.3.2
+Remotes:  
+    r-quantities/units, 
+    r-quantities/quantities

NAMESPACE

@@ -14,6 +14,7 @@ S3method(strata,seroincidence.by)
 S3method(summary,pop_data)
 S3method(summary,seroincidence)
 S3method(summary,seroincidence.by)
+export(ab_5p)
 export(as_curve_params)
 export(as_noise_params)
 export(as_pop_data)
@@ -43,6 +44,7 @@ export(set_id_var)
 export(sim.cs)
 export(sim.cs.multi)
 export(sim_pop_data)
+export(sim_pop_data_2)
 export(sim_pop_data_multi)
 export(strata)
 importFrom(Rcpp,sourceCpp)
@@ -62,10 +64,12 @@ importFrom(dplyr,if_any)
 importFrom(dplyr,if_else)
 importFrom(dplyr,inner_join)
 importFrom(dplyr,is.grouped_df)
+importFrom(dplyr,left_join)
 importFrom(dplyr,mutate)
 importFrom(dplyr,n)
 importFrom(dplyr,n_distinct)
 importFrom(dplyr,pull)
+importFrom(dplyr,reframe)
 importFrom(dplyr,relocate)
 importFrom(dplyr,rename)
 importFrom(dplyr,row_number)
@@ -98,6 +102,7 @@ importFrom(stats,pgamma)
 importFrom(stats,plnorm)
 importFrom(stats,qnorm)
 importFrom(stats,quantile)
+importFrom(stats,rexp)
 importFrom(stats,runif)
 importFrom(stats,xtabs)
 importFrom(tibble,as_tibble)

NEWS.md

@@ -2,6 +2,8 @@
 
 ## New features
 
+* Revised `methodology` vignette (#393)
+* Added support for units in `sim_pop_data()` (#393)
 * Improved warning messages for `get_biomarker_names_var()`
 * Added `get_*()` extractor functions to API (#380)
 * Added optional CI error bars to `autoplot.summary.seroincidence.by()` (#372)

R/ab_5p.R

@@ -0,0 +1,54 @@
+#' Calculate mean antibody response using 5-parameter model
+#'
+#' @param t time since last exposure
+#' @param y0 antibody response at exposure
+#' @param y1 antibody response at t1
+#' @param t1 time from last exposure to peak response
+#' @param alpha decay power function coefficient
+#' @param shape decay power function exponent
+#' [numeric] scalar; uses r > 1 scale for shape
+#'
+#' @returns a [numeric] [vector]
+#' @export
+#' @keywords internal
+#' @examples
+#' params <- typhoid_curves_nostrat_100[1, ]
+#' ab_5p(
+#'   t = units::as_units(50, "days"),
+#'   y0 = params$y0,
+#'   y1 = params$y1,
+#'   t1 = params$t1 |> units::as_units("days"),
+#'   alpha = params$alpha |> units::as_units("1/days"),
+#'   shape = params$r)
+ab_5p <- function(t, y0, y1, t1, alpha, shape) {
+  beta <- bt(y0, y1, t1)
+
+  yt <- if_else(
+    t < t1,
+    ab_5p_active_phase(t, y0, beta),
+    ab_5p_decay_phase(t, t1, y1, alpha, shape)
+  )
+  return(yt)
+}
+
+ab_5p_active_phase <- function(t, y0, beta) {
+  yt <- y0 * exp(units::set_units(beta * t, NULL))
+  return(yt)
+}
+
+ab_5p_decay_phase <- function(
+    t,
+    t1,
+    y1,
+    alpha,
+    shape) {
+  base_part_1 <- y1^(1 - shape)
+  base_part_2 <- units::set_units(
+    x = (1 - shape) * alpha * (t - t1),
+    value = NULL
+  )
+  base <- base_part_1 - base_part_2
+  expt <- (1 / (1 - shape))
+  ans <- base^expt
+  return(ans)
+}

R/f_dev.r

@@ -28,62 +28,7 @@ fdev <- function(lambda, csdata, lnpars, cond) {
 #' corresponding to input `lambda`
 #' @keywords internal
 #' @export
-#' @examples
-#' \donttest{
-#' library(dplyr)
-#' library(tibble)
-#'
-#' # load in longitudinal parameters
-#' curve_params <-
-#'   typhoid_curves_nostrat_100 %>%
-#'   filter(antigen_iso %in% c("HlyE_IgA", "HlyE_IgG"))
-#'
-#' # load in pop data
-#' xs_data <-
-#'   sees_pop_data_pk_100
-#'
-#' #Load noise params
-#' noise_params <- tibble(
-#'   antigen_iso = c("HlyE_IgG", "HlyE_IgA"),
-#'   nu = c(0.5, 0.5),                          # Biologic noise (nu)
-#'   eps = c(0, 0),                             # M noise (eps)
-#'   y.low = c(1, 1),                           # low cutoff (llod)
-#'   y.high = c(5e6, 5e6))                      # high cutoff (y.high)
-#'
-#' cur_antibody = "HlyE_IgA"
-#'
-#' cur_data <-
-#'   xs_data %>%
-#'   dplyr::filter(
-#'    .data$antigen_iso == cur_antibody) %>%
-#'   dplyr::slice_head(n = 100)
-#'
-#' cur_curve_params <-
-#'   curve_params %>%
-#'   dplyr::filter(.data$antigen_iso == cur_antibody) %>%
-#'   dplyr::slice_head(n = 100)
-#'
-#' cur_noise_params <-
-#'   noise_params %>%
-#'   dplyr::filter(.data$antigen_iso == cur_antibody)
-#'
-#' if (!is.element('d', names(cur_curve_params)))
-#' {
-#'   cur_curve_params <-
-#'     cur_curve_params %>%
-#'     dplyr::mutate(
-#'       alpha = .data$alpha * 365.25,
-#'       d = .data$r - 1)
-#' }
-#'
-#' lambda = 0.1
-#' f_dev0(
-#'     lambda = lambda,
-#'     csdata = cur_data,
-#'     lnpars = cur_curve_params,
-#'     cond = cur_noise_params
-#'   )
-#'}
+#' @example inst/examples/exm-f_dev.R
 f_dev0 <- function(lambda,
                    csdata,
                    lnpars,

R/graph.curve.params.R

@@ -8,6 +8,7 @@
 #' @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 log_x whether to log-transform the x-axis
 #'
 #' @returns a [ggplot2::ggplot()] object
 #' @export
@@ -30,7 +31,8 @@
   verbose = FALSE,
   show_quantiles = TRUE,
   show_all_curves = FALSE,
-  alpha_samples = 0.3
+  alpha_samples = 0.3,
+  log_x = FALSE
 ) {
   if (verbose) {
     message(
@@ -44,28 +46,10 @@
 
   tx2 <- 10^seq(-1, 3.1, 0.025)
 
-  bt <- function(y0, y1, t1) {
-    log(y1 / y0) / t1
+  if (inherits(curve_params$alpha, "units")) {
+    tx2 <- tx2 |> units::as_units("days")
   }
 
-  # 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
 
   dT <- # nolint: object_linter
@@ -92,14 +76,16 @@
     ) |>
     select(-"name") |>
     rowwise() |>
-    mutate(res = ab(
-      .data$t,
-      .data$y0,
-      .data$y1,
-      .data$t1,
-      .data$alpha,
-      .data$r
-    )) |>
+    mutate(
+      res = ab_5p(
+        .data$t,
+        .data$y0,
+        .data$y1,
+        .data$t1,
+        .data$alpha,
+        .data$r
+      )
+    ) |>
     ungroup()
 
   if (verbose) message("starting to compute quantiles")
@@ -206,6 +192,15 @@
       minor_breaks = NULL
     )
 
+  if (log_x) {
+    plot1 <-
+      plot1 +
+      ggplot2::scale_x_log10(
+        labels = scales::label_comma(),
+        minor_breaks = NULL
+      )
+  }
+
   if (show_quantiles) {
     plot1 <-
       plot1 +

R/graph.decay.curves.R

@@ -58,10 +58,18 @@
     # 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::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"))
+    ggplot2::theme(
+      plot.title = ggplot2::element_text(size = 20, face = "bold")
+    )
+
+  if (inherits(object$alpha, "units")) {
+    xlim <- xlim |> units::as_units("days")
+  }
 
   if (log_y) {
     plot1 <-
@@ -77,7 +85,8 @@
       fun = ab0,
       args = list(curve_params = cur_params),
       n = n_points,
-      ...)
+      ...
+    )
   }
 
   layers <-

R/serocalculator-package.R

@@ -18,10 +18,12 @@
 #' @importFrom dplyr if_else
 #' @importFrom dplyr inner_join
 #' @importFrom dplyr is.grouped_df
+#' @importFrom dplyr left_join
 #' @importFrom dplyr mutate
 #' @importFrom dplyr n
 #' @importFrom dplyr n_distinct
 #' @importFrom dplyr pull
+#' @importFrom dplyr reframe
 #' @importFrom dplyr relocate
 #' @importFrom dplyr rename
 #' @importFrom dplyr row_number
@@ -52,6 +54,7 @@
 #' @importFrom stats nlm
 #' @importFrom stats qnorm
 #' @importFrom stats quantile
+#' @importFrom stats rexp
 #' @importFrom stats runif
 #' @importFrom stats xtabs
 #' @importFrom tibble as_tibble

R/sim_age.R

@@ -0,0 +1,12 @@
+sim_age <- function(n_samples, age_range) {
+
+  age <- runif(n_samples,
+               min = age_range[1],
+               max = age_range[2])
+
+  if (inherits(age_range, "units")) {
+    age <- age |> units::as_units(units(age_range))
+  }
+
+  return(age)
+}

R/sim_pop_data.R

@@ -5,9 +5,9 @@
 #'  and adds noise, if desired.
 
 #' @param lambda a [numeric()] scalar indicating the incidence rate
-#' (in events per person-years)
+#' (in events per person-*year*s)
 #' @param n_samples number of samples to simulate
-#' @param age_range age range of sampled individuals, in years
+#' @param age_range age range of sampled individuals, in *years*
 #' @param age_fixed specify the curve parameters to use by age
 #' (does nothing at present?)
 #' @param antigen_isos Character vector with one or more antibody names.