Clean up documentation

R/compare_models.R

@@ -1,27 +1,44 @@
-#' Compare models
+#' Generate table of metrics for model comparison
 #'
 #' @param data input data to fit into the models
 #' @param method method to compare models. Can be one of the built-in methods or a function to compute the returned metrics (see Details).
 #' @param ... models to be compared. Must be models created by serosv. If models' names are not provided, indices will be used instead for the `model` column in the returned data.frame.
 #'
 #'
 #' @return
-#' a data.frame of 4 columns
+#' a data.frame with the following columns
 #'   \item{label}{name or index of the model}
 #'   \item{type}{model type of the given model (a serosv model name)}
-#'   \item{AIC}{AIC value for the model (lower value indicates better fit)}
-#'   \item{BIC}{BIC value for the model (lower value indicates better fit)}
+#'   \item{metrics columns}{the columns for metrics of comparison, the number of which depends on the function that generate these metrics}
 #'
 #' @details
 #' Built-in comparison methods include:
-#' - computing AIC and BIC, which returns AIC, BIC values of the model if available
-#' - cross validation, which returns MSE and logloss (negative )
+#' \itemize{
+#' \item{computing AIC and BIC, which returns AIC, BIC values of the model if available}
+#' \item{cross validation (perform k-fold validation), which returns MSE and logloss (negative log Binomial likelihood)
+#'  for aggregated data, or AUC and logloss (negative log Bernoulli likelihood)
+#'  for linelisting data}
+#' }
 #'
 #' @importFrom magrittr %>%
 #' @importFrom purrr imap_dfr as_mapper
 #' @importFrom stringr str_detect
 #' @importFrom assertthat assert_that
 #'
+#' @examples
+#' comparison_table <- suppressWarnings(
+#'   compare_models(
+#'     data = hav_bg_1964,
+#'     method = "CV",
+#'     polynomial_mod = ~polynomial_model(.x, k=1),
+#'     penalized_spline = penalized_spline_model,
+#'     farrington = ~farrington_model(.x, start=list(alpha=0.3,beta=0.1,gamma=0.03))
+#'   )
+#' )
+#' # view table of metrics
+#' comparison_table
+#' # view the model fitted with the whole dataset
+#' comparison_table$plots
 #' @export
 compare_models <- function(data, method="AIC/BIC",...){
   list(...) %>%
@@ -118,15 +135,11 @@ cv <- function(dat, mod_func, k=4){
 
     }else{
       # if data is linelisting
-      # make sure pred is slightly higher than 0 and lower than 1 to avoid log(0)
-      eps <- .Machine$double.eps  # smallest positive floating point number ~ 2.2e-16
-      pred <- pmax(eps, pmin(1 - eps, pred))
-
       # compute logloss (negative bernoulli loglikelihood)
-      curr_metric$logloss <- -sum(test_dat$status*log(pred) + (1 - test_dat$status)*log(1-pred),
+      curr_metric$logloss <- -sum(dbinom(test_dat$status, 1, prob=pred, log=TRUE),
                                   na.rm = TRUE)
       # and estimate auc
-      curr_metric$auc <- as.numeric(pROC::auc(test_dat$status, pred))
+      curr_metric$auc <- as.numeric(pROC::auc(test_dat$status, pred, quiet=TRUE))
     }
 
     curr_metric

R/mixture_model.R

@@ -152,9 +152,14 @@ estimate_from_mixture <- function(age, antibody_level, threshold_status = NULL,
   model <- list()
 
   # sort, just in case
-  age <- age[order(age)]
   antibody_level <- antibody_level[order(age)]
   log_antibody <- log(antibody_level + 1)
+  threshold_status <- if(!is.null(threshold_status)){
+    threshold_status[order(age)]
+  }else{
+    NULL
+  }
+  age <- age[order(age)]
 
   # get mu_s and mu_i from mixture model
   mu_s <- mixture_model$info$parameters$mu[1]

R/plots.R

@@ -1,11 +1,11 @@
 #' Helper to adjust styling of a plot
 #'
-#' @param sero - color for seroprevalence line
-#' @param ci - color for confidence interval
-#' @param foi - color for force of infection line
-#' @param sero_line - linetype for seroprevalence line
-#' @param foi_line - linetype for force of infection line
-#' @param xlabel - x label
+#' @param sero color for seroprevalence line
+#' @param ci color for confidence interval
+#' @param foi color for force of infection line
+#' @param sero_line linetype for seroprevalence line
+#' @param foi_line linetype for force of infection line
+#' @param xlabel x label
 #'
 #' @return list of updated aesthetic values
 #' @export

R/semiparametric_models.R

@@ -128,13 +128,23 @@ penalized_spline_model <- function(data,
   model$datatype <- data$type
 
   # s <- mgcv:::s
-  spos <- pos/tot
+  neg <- tot - pos
 
   if (framework == "pl"){
-    model$info <- mgcv::gam(spos ~ s(age, bs = s, sp=sp), family = binomial(link = link))
+    model$info <- if(data$type == "aggregated"){
+        mgcv::gam(cbind(pos, neg) ~ s(age, bs = s, sp=sp), family = binomial(link = link))
+      }else{
+        mgcv::gam(pos ~ s(age, bs = s, sp=sp), family = binomial(link = link))
+      }
+
     model$sp <- model$info$fitted.values
   }else if(framework == "glmm"){
-    model$info <- mgcv::gamm(spos ~ s(age, bs = s, sp=sp), family = binomial(link = link))
+    model$info <- if(data$type == "aggregated"){
+        mgcv::gamm(cbind(pos, neg) ~ s(age, bs = s, sp=sp), family = binomial(link = link))
+      }else{
+        mgcv::gamm(pos ~ s(age, bs = s, sp=sp), family = binomial(link = link))
+      }
+
     model$sp <- model$info$gam$fitted.values
   }else{
     stop(paste0('Invalid value for framework. Expected "pl" or "glmm", got ', framework))

R/to_titer.R

@@ -81,7 +81,9 @@ to_titer <- function(df, model="4PL", positive_threshold=NULL, ci = .95,
   out
 }
 
-#' Preprocess data
+#' Standardize raw serological test data for titer conversion
+#'
+#' Validate and prepare raw serological test results for use with `to_titer()`
 #'
 #' @param df data.frame with columns for plate id, sample id, result, dilution factor, and (optionally) negative controls
 #' @param plate_id_col name of the column storing plates id
@@ -168,16 +170,21 @@ data2function <- function(df) {
 
 # Function to convert samples' OD to LC
 #' @importFrom purrr map_dfr
+#' @import dplyr
 process_samples <- function(plate, std_crv, midpoint=2, positive_threshold = NULL) {
+  # convert from OD to concentration (with dilution taken into account)
   out <- plate %>%
-    bind_cols(purrr::map_dfr(.$result, std_crv))
+    bind_cols(purrr::map_dfr(.$result, std_crv)) %>%
+    # after mapping from OD to LC
+    # multiply by the dilution factors to get the sample actual concentration
+    mutate(across(c(lower, median, upper), ~ 10^.x * dilution_factors))
 
   if(is.numeric(positive_threshold)){
     out <- out %>% mutate(
       positive = if_else(
         # check if upper bound for titer is available
         !is.na(upper),
-        upper >= positive_threshold, # positive if above threshold
+        upper >= (positive_threshold), # positive if above threshold
         # if upper is NA, then OD (or result) is either too high or too low
         # label positive when result is higher than midpoint
         result > midpoint
@@ -238,8 +245,8 @@ get_negative_controls <- function(plate, std_crv){
   plate %>%
     select(starts_with("negative")) %>%
     unique() %>%
-    tidyr::pivot_longer(everything(), names_to = "dilution", values_to = "result") %>%
-    mutate(across(dilution, ~ stringr::str_remove(.x, "negative_") %>%  as.integer()))
+    tidyr::pivot_longer(everything(), names_to = "dilution_factors", values_to = "result") %>%
+    mutate(across(dilution_factors, ~ stringr::str_remove(.x, "negative_") %>%  as.integer()))
 }
 
 # ======== Plot functions ==========
@@ -339,10 +346,10 @@ add_thresholds <- function(dilution_factors, positive_threshold = 0.1,
 #' The figure below demonstrates the interpretation of the plot.
 #' \figure{interpret_titer_qc.png}{options: width="70\%"}
 #'
-#' @param x - output of `to_titer()`
-#' @param n_plates - maximum number of plates to plot
-#' @param n_samples - maximum number of samples per plate to plot
-#' @param n_dilutions - number of dilutions used for testing
+#' @param x output of `to_titer()`
+#' @param n_plates maximum number of plates to plot
+#' @param n_samples maximum number of samples per plate to plot
+#' @param n_dilutions number of dilutions used for testing
 #'
 #' @importFrom magrittr %>%
 #' @importFrom purrr walk

R/utils.R

@@ -1,8 +1,8 @@
 
 #' Estimate force of infection
 #'
-#' @param t - time (in this case age) vector
-#' @param sp - seroprevalence vector
+#' @param t time (in this case age) vector
+#' @param sp seroprevalence vector
 #'
 #' @importFrom stats approx
 #'
@@ -55,7 +55,6 @@ pava<- function(pos=pos,tot=rep(1,length(pos)))
   return(list(pai1=pai1,pai2=pai2))
 }
 
-# TODO: update aggregate func here to be "pipe"-able
 #' Aggregate data
 #'
 #' Generate a dataframe with `t`, `pos` and `tot` columns from

docs/articles/model_selection.md

@@ -91,7 +91,7 @@ cv_out
     ##   <chr>              <dbl> <dbl> <chr>                  <list>     <list>
     ## 1 griffith            277. 0.594 polynomial_model       <plynml_m> <gg>  
     ## 2 penalized_spline    197. 0.601 penalized_spline_model <pnlzd_s_> <gg>  
-    ## 3 farrington          356. 0.600 farrington_model       <frrngtn_> <gg>  
+    ## 3 farrington          185. 0.600 farrington_model       <frrngtn_> <gg>  
     ## 4 local_polynomial    237. 0.588 lp_model               <lp_model> <gg>
 
 With aggregated data
@@ -117,7 +117,6 @@ cv_out <- compare_models(
     local_polynomial = lp_model 
   ) %>% suppressWarnings()
 
-
 aic_bic_out
 ```
 

docs/llms.txt

@@ -297,9 +297,9 @@ parvob19 |>
 - [`transform_data()`](https://oucru-modelling.github.io/serosv/reference/transform_data.md)
   : Aggregate data
 - [`standardize_data()`](https://oucru-modelling.github.io/serosv/reference/standardize_data.md)
-  : Preprocess data
+  : Standardize raw serological test data for titer conversion
 - [`compare_models()`](https://oucru-modelling.github.io/serosv/reference/compare_models.md)
-  : Compare models
+  : Generate table of metrics for model comparison
 - [`compute_ci(`*`<default>`*`)`](https://oucru-modelling.github.io/serosv/reference/compute_ci.default.md)
   : Compute confidence interval for a model of serosv
 - [`compute_ci(`*`<fp_model>`*`)`](https://oucru-modelling.github.io/serosv/reference/compute_ci.fp_model.md)