From ee4c8e933c61f527999516e3dcea5b9b44832b51 Mon Sep 17 00:00:00 2001 From: Kwan-Jenny <68584166+Kwan-Jenny@users.noreply.github.com> Date: Sat, 31 Jan 2026 13:32:28 -0800 Subject: [PATCH 01/30] Delete R/create_incidence_table.R --- R/create_incidence_table.R | 17 ----------------- 1 file changed, 17 deletions(-) delete mode 100644 R/create_incidence_table.R diff --git a/R/create_incidence_table.R b/R/create_incidence_table.R deleted file mode 100644 index 6a965f4..0000000 --- a/R/create_incidence_table.R +++ /dev/null @@ -1,17 +0,0 @@ -# create table of incidence.rate of each region -create_incidence_table <- function(...) { - # Capture input objects and their names - est_list <- list(...) - country_names <- names(est_list) - - # Extract incidence.rate from the summary() of each estimate object - incidence_rates <- sapply(est_list, function(x) summary(x)$incidence.rate) - - # Create a tidy tibble - incidence_table <- tibble( - Country = country_names, - Incidence_Rate = incidence_rates - ) - - return(incidence_table) -} \ No newline at end of file From 4207ef5da238a59553371ede4ba06ae89d754690 Mon Sep 17 00:00:00 2001 From: Kwan-Jenny <68584166+Kwan-Jenny@users.noreply.github.com> Date: Sat, 31 Jan 2026 13:32:40 -0800 Subject: [PATCH 02/30] Delete R/create_noise_df.R --- R/create_noise_df.R | 12 ------------ 1 file changed, 12 deletions(-) delete mode 100644 R/create_noise_df.R diff --git a/R/create_noise_df.R b/R/create_noise_df.R deleted file mode 100644 index f4640d6..0000000 --- a/R/create_noise_df.R +++ /dev/null @@ -1,12 +0,0 @@ -create_noise_df <- function(country) { - noise_df <- tibble( - antigen_iso = c("IgG"), - Country = factor(c(country)), # Use input argument for Country - y.low = c(25), - eps = c(0.25), - nu = c(0.5), - y.high = c(200000) - ) - - return(noise_df) -} \ No newline at end of file From a05bd51d5a5f7508ed9621043fd2be32394a9178 Mon Sep 17 00:00:00 2001 From: Kwan-Jenny <68584166+Kwan-Jenny@users.noreply.github.com> Date: Sat, 31 Jan 2026 13:34:07 -0800 Subject: [PATCH 03/30] Delete R/plot_summary_metrics.R --- R/plot_summary_metrics.R | 14 -------------- 1 file changed, 14 deletions(-) delete mode 100644 R/plot_summary_metrics.R diff --git a/R/plot_summary_metrics.R b/R/plot_summary_metrics.R deleted file mode 100644 index ce58750..0000000 --- a/R/plot_summary_metrics.R +++ /dev/null @@ -1,14 +0,0 @@ -plot_summary_metrics <- function(summary_metrics) { - summary_metrics |> - ggplot() + - aes(x = sample_size, y = empirical_se, color = Age_Group) + - geom_line() + - geom_point() + - labs( - title = "Empirical Standard Error vs. Sample Size", - x = "Sample Size", - y = "Empirical Standard Error", - color = "Age Group" - ) + - theme_minimal() -} \ No newline at end of file From 1b39d19ffde0e9ebb397b5c248c8f6b7813fd057 Mon Sep 17 00:00:00 2001 From: Kwan-Jenny <68584166+Kwan-Jenny@users.noreply.github.com> Date: Sat, 31 Jan 2026 13:34:17 -0800 Subject: [PATCH 04/30] Delete R/simulate_seroincidence.R --- R/simulate_seroincidence.R | 98 -------------------------------------- 1 file changed, 98 deletions(-) delete mode 100644 R/simulate_seroincidence.R diff --git a/R/simulate_seroincidence.R b/R/simulate_seroincidence.R deleted file mode 100644 index 1c7908b..0000000 --- a/R/simulate_seroincidence.R +++ /dev/null @@ -1,98 +0,0 @@ -# Define the simulation function -simulate_seroincidence <- function( - dmcmc, # Curve parameters - nrep, - n_sim, - observed, - range = NULL, - batch_size = 40, parallel = TRUE, - antibodies = c("IgG"), # Antigen-isotypes - dlims = rbind("IgG" = c(min = 0, max = 0.5)), # Biologic noise distribution - # Noise parameters - cond = tibble( - antigen_iso = c("IgG"), - nu = c(0.5), - eps = c(0.25), - y.low = c(25), - y.high = c(200000) - ) -) { - - lambda <- observed # Simulated incidence rate per person-year - - - - # Calculate number of batches - n_batches <- ceiling(n_sim / batch_size) - - # Parallel processing: Use future_map for efficiency - if (parallel) { - plan(multisession, workers = max(1, future::availableCores() - 1)) - - results <- future_map(1:n_batches, function(batch) { - # Run simulations within each batch - replicate(batch_size, expr = { - csdata <- sim_pop_data( - curve_params = dmcmc, - lambda = lambda, - n_samples = nrep, - age_range = range, - antigen_isos = antibodies, - n_mcmc_samples = 0, - renew_params = FALSE, # Keep the same params for speed - add_noise = TRUE, - noise_limits = dlims, - format = "long" - ) - - # Estimate seroincidence - est <- est.incidence( - pop_data = csdata, - curve_params = dmcmc, - noise_params = cond, - lambda_start = 0.1, - build_graph = FALSE, # Disable visualization for speed - verbose = FALSE, - print_graph = FALSE, - antigen_isos = antibodies - ) - - list(csdata = csdata, est1 = est) - }, simplify = FALSE) # Keep list structure - }, .options = furrr_options(seed = TRUE)) # Set seed for reproducibility - - # Flatten nested list - results <- unlist(results, recursive = FALSE) - } else { - # Sequential execution (for debugging or single-core use) - results <- lapply(1:n_sim, function(i) { - csdata <- sim.cs( - curve_params = dmcmc, - lambda = lambda, - n.smpl = nrep, - age_range = range, - antigen_isos = antibodies, - n.mc = 0, - renew_params = FALSE, - add.noise = TRUE, - noise_limits = dlims, - format = "long" - ) - - est <- est.incidence( - pop_data = csdata, - curve_params = dmcmc, - noise_params = cond, - lambda_start = 0.1, - build_graph = FALSE, - verbose = FALSE, - print_graph = FALSE, - antigen_isos = antibodies - ) - - list(csdata = csdata, est1 = est) - }) - } - - return(results) -} \ No newline at end of file From 992fdb0e5d764f2919fb934466cd5e01ada81b7b Mon Sep 17 00:00:00 2001 From: Kwan-Jenny <68584166+Kwan-Jenny@users.noreply.github.com> Date: Sat, 31 Jan 2026 13:34:27 -0800 Subject: [PATCH 05/30] Delete R/prepare_df_for_serocalculator.R --- R/prepare_df_for_serocalculator.R | 22 ---------------------- 1 file changed, 22 deletions(-) delete mode 100644 R/prepare_df_for_serocalculator.R diff --git a/R/prepare_df_for_serocalculator.R b/R/prepare_df_for_serocalculator.R deleted file mode 100644 index 00f4c29..0000000 --- a/R/prepare_df_for_serocalculator.R +++ /dev/null @@ -1,22 +0,0 @@ -# Create function to clearly define cross-sectional data -prepare_df_for_serocalculator <- function(df, age_col = "age", value_col = "value") { - # Ensure correct column names - df <- df %>% - rename(age = all_of(age_col)) - - # Assign attributes for serocalculator - attr(df, "age_var") <- "age" - attr(df, "value_var") <- value_col - - # Check if serocalculator recognizes attributes - get_value_var <- serocalculator:::get_value_var - detected_value_var <- get_value_var(df) - - if (is.null(detected_value_var)) { - warning("serocalculator did not detect the 'value' column. Check column naming.") - } else { - message("serocalculator recognized 'value' column: ", detected_value_var) - } - - return(df) -} \ No newline at end of file From 94f3bd94ea6fcb6ea0fcca35a07a22e3e56f76d4 Mon Sep 17 00:00:00 2001 From: Kwan-Jenny <68584166+Kwan-Jenny@users.noreply.github.com> Date: Sat, 31 Jan 2026 13:34:37 -0800 Subject: [PATCH 06/30] Delete R/postprocess_jags_output.R --- R/postprocess_jags_output.R | 47 ------------------------------------- 1 file changed, 47 deletions(-) delete mode 100644 R/postprocess_jags_output.R diff --git a/R/postprocess_jags_output.R b/R/postprocess_jags_output.R deleted file mode 100644 index 32bfd22..0000000 --- a/R/postprocess_jags_output.R +++ /dev/null @@ -1,47 +0,0 @@ -#' Postprocess JAGS output -#' -#' @param jags_output output from [runjags::run.jags()] -#' -#' @returns a [tibble::tbl_df] -#' @export -#' -postprocess_jags_output <- function(jags_output) { - mcmc_list <- coda::as.mcmc.list(jags_output) - - mcmc_df <- ggmcmc::ggs(mcmc_list) - - wide_predpar_df <- mcmc_df |> - mutate( - parameter = sub("^(\\w+)\\[.*", "\\1", .data$Parameter), - index_id = as.numeric(sub("^\\w+\\[(\\d+),.*", "\\1", .data$Parameter)), - antigen_iso = as.numeric(sub("^\\w+\\[\\d+,(\\d+).*", "\\1", .data$Parameter)) - ) |> - mutate( - index_id = factor(index_id, labels = c(unique(dL_clean$index_id), "newperson")), - antigen_iso = factor(antigen_iso, labels = unique(dL_clean$antigen_iso)) - ) |> - filter(index_id == "newperson") |> - select(-Parameter) |> - pivot_wider(names_from = "parameter", values_from = "value") |> - rowwise() |> - droplevels() |> - ungroup() |> - rename(r = shape) - - # Assuming wide_predpar_df is your data frame - curve_params <- wide_predpar_df - - # Set class and attributes for serocalculator - class(curve_params) <- c("curve_params", class(curve_params)) - antigen_isos <- unique(curve_params$antigen_iso) - attr(curve_params, "antigen_isos") <- antigen_isos - - to_return <- curve_params |> - mutate( - iter = Iteration, - chain = Chain, - ) |> - select(antigen_iso, iter, chain, y0, y1, t1, alpha, r) - - return(to_return) -} From 19eb4d202b9309af2f5e426d0c5d80b4c95ce29a Mon Sep 17 00:00:00 2001 From: Kwan-Jenny <68584166+Kwan-Jenny@users.noreply.github.com> Date: Sat, 31 Jan 2026 13:34:46 -0800 Subject: [PATCH 07/30] Delete R/generate_final_table.R --- R/generate_final_table.R | 28 ---------------------------- 1 file changed, 28 deletions(-) delete mode 100644 R/generate_final_table.R diff --git a/R/generate_final_table.R b/R/generate_final_table.R deleted file mode 100644 index d84a3f0..0000000 --- a/R/generate_final_table.R +++ /dev/null @@ -1,28 +0,0 @@ -# Define a function to generate final tables -generate_final_table <- function(results_list, sample_size) { - # Initialize an empty list to store the results - summary_results <- list() - - # Loop through each of the 100 results and extract the required columns - for (i in 1:200) { - # Extract the summary for each result - result_summary <- summary(results_list[[i]]$est1) - - # Select the required columns - extracted_columns <- result_summary %>% - select(incidence.rate, SE, CI.lwr, CI.upr) - - # Add a column for the index (optional, for tracking) - extracted_columns <- extracted_columns %>% - mutate(index = i) - - # Append to the list - summary_results[[i]] <- extracted_columns - } - - # Combine all results into a single data frame - final_table <- bind_rows(summary_results) %>% - mutate(sample_size = sample_size) # Add sample size column for clarity - - return(final_table) -} From fb6f507154e42704ecf908baf4e47966462919a2 Mon Sep 17 00:00:00 2001 From: Kwan-Jenny <68584166+Kwan-Jenny@users.noreply.github.com> Date: Sat, 31 Jan 2026 13:34:56 -0800 Subject: [PATCH 08/30] Delete R/create_xs_data.R --- R/create_xs_data.R | 29 ----------------------------- 1 file changed, 29 deletions(-) delete mode 100644 R/create_xs_data.R diff --git a/R/create_xs_data.R b/R/create_xs_data.R deleted file mode 100644 index 9a0f0a3..0000000 --- a/R/create_xs_data.R +++ /dev/null @@ -1,29 +0,0 @@ -# create function for generating specific region data -create_xs_data <- function(df, filter_countries, filter_antigen_iso, value_col) { - df %>% - # First, create/rename columns - mutate( - id = sid, - Country = site_name, - study = study_name, - age = age, - antigen_iso = factor(isotype_name), - value = {{ value_col } - } # value_col is specified by the user, e.g. n_ipab_MFI - ) %>% - # Then filter by the desired catchment(s) and antigen iso(s) - filter( - Country %in% filter_countries, - antigen_iso %in% filter_antigen_iso - ) %>% - # Create age categories - mutate( - ageCat = factor(case_when( - age < 5 ~ "<5", - age >= 5 & age <= 15 ~ "5-15", - age > 15 ~ "16+" - )) - ) %>% - # Optionally select only the needed columns - select(id, Country, study, age, antigen_iso, value, ageCat) -} \ No newline at end of file From 40ddcd8c662197b90f1820211af2c13d4af243e1 Mon Sep 17 00:00:00 2001 From: Kwan-Jenny <68584166+Kwan-Jenny@users.noreply.github.com> Date: Sat, 31 Jan 2026 13:37:11 -0800 Subject: [PATCH 09/30] Delete R/process_shigella_data.R --- R/process_shigella_data.R | 27 --------------------------- 1 file changed, 27 deletions(-) delete mode 100644 R/process_shigella_data.R diff --git a/R/process_shigella_data.R b/R/process_shigella_data.R deleted file mode 100644 index 8bb1d01..0000000 --- a/R/process_shigella_data.R +++ /dev/null @@ -1,27 +0,0 @@ -process_shigella_data <- function(data, study_filter, antigen) { - # Filter the data for the specific study - filtered_data <- data %>% - filter(study_name == study_filter) - - # Capture the column name of the antigen - antigen_col <- ensym(antigen) - - # Manipulate and restructure the data - processed_data <- filtered_data %>% - select(isotype_name, sid, timepoint, `Actual day`, !!antigen_col) %>% - mutate( - index_id = sid, - antigen_iso = isotype_name, - visit = timepoint, - timeindays = `Actual day`, - result = !!antigen_col - ) %>% - group_by(index_id, antigen_iso) %>% - arrange(visit) %>% - mutate(visit_num = rank(visit, ties.method = "first")) %>% - ungroup() %>% - # Remove rows with NA in timeindays - filter(!is.na(timeindays)) - - return(processed_data) -} \ No newline at end of file From 9d1b8acf3d5e67ccdeb98c85ac2c7813a293412a Mon Sep 17 00:00:00 2001 From: Kwan-Jenny <68584166+Kwan-Jenny@users.noreply.github.com> Date: Sat, 31 Jan 2026 13:38:12 -0800 Subject: [PATCH 10/30] Add process_shigella_data function for data reshaping This function filters and reshapes longitudinal Shigella data for serodynamics workflows, standardizing column names and preparing the dataset for conversion. --- R/process_shigella_data.R | 49 +++++++++++++++++++++++++++++++++++++++ 1 file changed, 49 insertions(+) create mode 100644 R/process_shigella_data.R diff --git a/R/process_shigella_data.R b/R/process_shigella_data.R new file mode 100644 index 0000000..215222f --- /dev/null +++ b/R/process_shigella_data.R @@ -0,0 +1,49 @@ +#' Reshape Shigella longitudinal data for serodynamics workflows +#' +#' Filters a dataset to a given study and antigen column, standardizes column +#' names, and returns a visit-ordered long dataset suitable for conversion to +#' `serodynamics::as_case_data()`. +#' +#' @param data A data frame containing longitudinal measurements. +#' @param study_filter Character scalar. Value of `study_name` to keep (e.g. "SOSAR"). +#' @param antigen Unquoted column name for the antigen measurement (e.g. n_ipab_MFI). +#' +#' @return A tibble with standardized columns: +#' \describe{ +#' \item{index_id}{Participant ID (copied from `sid`).} +#' \item{antigen_iso}{Isotype label (copied from `isotype_name`).} +#' \item{visit}{Visit label (copied from `timepoint`).} +#' \item{timeindays}{Time since infection (copied from `Actual day`).} +#' \item{result}{Antibody measurement (from `antigen`).} +#' } +#' +#' @examples +#' \dontrun{ +#' dat_long <- process_shigella_data(df, "SOSAR", n_ipab_MFI) +#' dL <- serodynamics::as_case_data(dat_long, +#' id_var = "index_id", biomarker_var = "antigen_iso", +#' time_in_days = "timeindays", value_var = "result" +#' ) +#' } +#' +#' @export +process_shigella_data <- function(data, study_filter, antigen) { + antigen_col <- rlang::ensym(antigen) + + data |> + dplyr::filter(.data$study_name == study_filter) |> + dplyr::select(.data$isotype_name, .data$sid, .data$timepoint, .data$`Actual day`, + .data$cohort_name, !!antigen_col) |> + dplyr::mutate( + index_id = .data$sid, + antigen_iso = .data$isotype_name, + visit = .data$timepoint, + timeindays = .data$`Actual day`, + result = !!antigen_col + ) |> + dplyr::group_by(.data$index_id, .data$antigen_iso) |> + dplyr::arrange(.data$visit, .by_group = TRUE) |> + dplyr::mutate(visit_num = rank(.data$visit, ties.method = "first")) |> + dplyr::ungroup() |> + dplyr::filter(!is.na(.data$timeindays)) +} From 58db1c78c9e6d497402f1f0ec31d11c795e2f061 Mon Sep 17 00:00:00 2001 From: Kwan-Jenny <68584166+Kwan-Jenny@users.noreply.github.com> Date: Sat, 31 Jan 2026 13:39:01 -0800 Subject: [PATCH 11/30] Add function to summarize antibody trajectories for newperson This function summarizes population-level antibody trajectories for a hypothetical new individual, computing median and credible intervals. It allows customization of various parameters for plotting and data return options. --- R/fig2_overall_newperson.R | 106 +++++++++++++++++++++++++++++++++++++ 1 file changed, 106 insertions(+) create mode 100644 R/fig2_overall_newperson.R diff --git a/R/fig2_overall_newperson.R b/R/fig2_overall_newperson.R new file mode 100644 index 0000000..eccba0b --- /dev/null +++ b/R/fig2_overall_newperson.R @@ -0,0 +1,106 @@ +#' Summarize population-level ("newperson") antibody trajectories from overall models +#' +#' Computes median and credible interval bands of the antibody trajectory for a +#' hypothetical new individual drawn from the population distribution ("newperson"). +#' +#' @param overall_models Named list of posterior draws in long format (one per antigen), +#' with columns at least: Subject, Iso_type, Chain, Iteration, Parameter, value. +#' @param osps Character vector of antigen names (must match names in `overall_models`). +#' @param ids Character vector of Subject IDs to include (default: "newperson"). +#' @param isotypes Character vector of isotypes (default: c("IgG","IgA")). +#' @param t_grid Numeric vector of time points (days) to evaluate. +#' @param cred Credible level (default 0.95). +#' @param log_y Logical; if TRUE, applies log10 scale to y when returning plot. +#' @param xlim Optional numeric length-2 vector for x-axis limits when returning plot. +#' @param ylab Y-axis label for plot. +#' @param line_color Line color for plot. +#' @param ribbon_alpha Alpha for credible ribbon in plot. +#' @param facet_scales Passed to ggplot facet scales. +#' @param return_data If TRUE, returns a list with `plot` and `data`. +#' +#' @return By default, a ggplot. If `return_data = TRUE`, returns list(plot = p, data = df). +#' +#' @details +#' Requires that each model draw can be pivoted to wide parameters including: +#' y0, y1, t1, alpha, shape. +#' +#' @export +fig2_overall_newperson <- function( + overall_models, + osps = c("IpaB","Sf2a","Sf3a","Sf6","Sonnei"), + ids = "newperson", + isotypes = c("IgG","IgA"), + t_grid = seq(0, 210, by = 5), + cred = 0.95, + log_y = TRUE, + xlim = c(0, 210), + ylab = "Normalized MFI", + line_color = "#1f77b4", + ribbon_alpha = 0.20, + facet_scales = "fixed", + return_data = FALSE +) { + q_lo <- (1 - cred) / 2 + q_hi <- 1 - q_lo + + get_sum_overall <- function(model_df, osp, iso) { + model_df |> + dplyr::filter(.data$Subject %in% ids, .data$Iso_type == iso) |> + dplyr::select(.data$Chain, .data$Iteration, .data$Iso_type, + .data$Parameter, .data$value, .data$Subject) |> + tidyr::pivot_wider(names_from = .data$Parameter, values_from = .data$value) |> + dplyr::mutate( + antigen = osp, + iso = factor(.data$Iso_type, levels = c("IgG","IgA")) + ) |> + tidyr::crossing(t = t_grid) |> + dplyr::mutate( + res = serodynamics:::ab( + t = .data$t, + y0 = .data$y0, + y1 = .data$y1, + t1 = .data$t1, + alpha = .data$alpha, + shape = .data$shape + ) + ) |> + dplyr::filter(is.finite(.data$res), .data$res > 0) |> + dplyr::summarise( + res.med = stats::quantile(.data$res, 0.50, na.rm = TRUE), + res.low = stats::quantile(.data$res, q_lo, na.rm = TRUE), + res.high = stats::quantile(.data$res, q_hi, na.rm = TRUE), + .by = c("antigen","iso","t") + ) + } + + sum_all <- purrr::map_dfr(osps, function(osp) { + purrr::map_dfr(isotypes, function(iso) { + get_sum_overall(overall_models[[osp]], osp, iso) + }) + }) |> + dplyr::mutate(antigen = factor(.data$antigen, levels = osps)) + + p <- ggplot2::ggplot(sum_all, ggplot2::aes(x = .data$t)) + + ggplot2::geom_ribbon( + ggplot2::aes(ymin = .data$res.low, ymax = .data$res.high), + alpha = ribbon_alpha + ) + + ggplot2::geom_line( + ggplot2::aes(y = .data$res.med), + linewidth = 1.1, + color = line_color, + lineend = "round" + ) + + ggplot2::labs(x = "Days since fever onset", y = ylab) + + ggplot2::facet_grid(rows = ggplot2::vars(antigen), + cols = ggplot2::vars(iso), + scales = facet_scales) + + ggplot2::theme_minimal() + + ggplot2::theme(legend.position = "none") + + if (log_y) p <- p + ggplot2::scale_y_log10() + if (!is.null(xlim)) p <- p + ggplot2::coord_cartesian(xlim = xlim) + + if (return_data) return(list(plot = p, data = sum_all)) + p +} From 137a59858d3ac6fc929bce622d8792acccf3ad8c Mon Sep 17 00:00:00 2001 From: Kwan-Jenny <68584166+Kwan-Jenny@users.noreply.github.com> Date: Sat, 31 Jan 2026 22:32:43 -0800 Subject: [PATCH 12/30] Add function to predict posterior at specified times This function generates posterior predictions of antibody trajectories at specified time points for selected subjects and antigen/isotype. It processes posterior draws and evaluates the antibody curve using specified parameters. --- R/predict_posterior_at_times.R | 79 ++++++++++++++++++++++++++++++++++ 1 file changed, 79 insertions(+) create mode 100644 R/predict_posterior_at_times.R diff --git a/R/predict_posterior_at_times.R b/R/predict_posterior_at_times.R new file mode 100644 index 0000000..5b59a3b --- /dev/null +++ b/R/predict_posterior_at_times.R @@ -0,0 +1,79 @@ +#' Posterior predictions at specified times for given subjects and antigen/isotype +#' +#' Generates draw-level posterior predictions of the antibody trajectory at user-specified +#' time points, for one or more subjects and a selected antigen/isotype. This is a low-level +#' helper used by residual-based posterior predictive diagnostics. +#' +#' @param model A data frame of posterior draws in long format with columns: +#' \code{Subject}, \code{Iso_type}, \code{Chain}, \code{Iteration}, \code{Parameter}, \code{value}. +#' @param ids Character vector of subject IDs to include (matched against \code{Subject}). +#' @param antigen_iso Character scalar specifying the antigen/isotype to include +#' (matched against \code{Iso_type}). +#' @param times Numeric vector of time points (days) at which to evaluate predictions. +#' +#' @return A tibble with one row per (posterior draw \eqn{\times} time \eqn{\times} subject), +#' including the evaluated prediction \code{res}. Output includes at least: +#' \describe{ +#' \item{id}{Subject ID (character).} +#' \item{t}{Time (days) at which prediction was evaluated.} +#' \item{Chain}{MCMC chain index (if present in \code{model}).} +#' \item{Iteration}{MCMC iteration index (if present in \code{model}).} +#' \item{sample_id}{Row index for the draw (added if missing).} +#' \item{y0, y1, t1, alpha, shape}{Model parameters (wide).} +#' \item{res}{Predicted antibody level at time \code{t}.} +#' } +#' +#' @details +#' This function pivots posterior draws to wide format (parameters as columns), +#' expands them over \code{times}, and evaluates the antibody curve via +#' \code{serodynamics:::ab()} using parameters \code{y0}, \code{y1}, \code{t1}, +#' \code{alpha}, and \code{shape}. +#' +#' @seealso \code{\link{compute_residual_metrics}} +#' +#' @examples +#' \dontrun{ +#' preds <- predict_posterior_at_times( +#' model = overall_sf2a, +#' ids = "newperson", +#' antigen_iso = "IgG", +#' times = c(0, 30, 90, 180) +#' ) +#' } +#' +#' @keywords internal +predict_posterior_at_times <- function(model, ids, antigen_iso, times) { + + sr_model_sub <- model |> + dplyr::filter(.data$Subject %in% ids, .data$Iso_type == antigen_iso) + + param_wide <- sr_model_sub |> + dplyr::select(.data$Chain, .data$Iteration, .data$Iso_type, + .data$Parameter, .data$value, .data$Subject) |> + tidyr::pivot_wider(names_from = .data$Parameter, values_from = .data$value) |> + dplyr::arrange(.data$Chain, .data$Iteration) |> + dplyr::mutate( + antigen_iso = factor(.data$Iso_type), + id = as.character(.data$Subject) + ) |> + dplyr::select(-c("Iso_type", "Subject")) + + if (!"sample_id" %in% names(param_wide)) { + param_wide <- param_wide |> + dplyr::mutate(sample_id = dplyr::row_number()) + } + + dt <- tibble::tibble(t = times) |> + dplyr::mutate(idx = dplyr::row_number()) |> + tidyr::pivot_wider(names_from = "idx", values_from = "t", names_prefix = "time") |> + dplyr::slice(rep(seq_len(dplyr::n()), each = nrow(param_wide))) + + predictions <- cbind(param_wide, dt) |> + tidyr::pivot_longer(cols = dplyr::starts_with("time"), values_to = "t") |> + dplyr::select(-"name") |> + dplyr::mutate( + res = serodynamics:::ab(.data$t, .data$y0, .data$y1, .data$t1, .data$alpha, .data$shape) + ) + + predictions +} From f054420d81ea1d2a1ce2ac26cb851769a913a606 Mon Sep 17 00:00:00 2001 From: Kwan-Jenny <68584166+Kwan-Jenny@users.noreply.github.com> Date: Sat, 31 Jan 2026 22:33:42 -0800 Subject: [PATCH 13/30] Add compute_residual_metrics function for antibody data This function computes residuals between observed antibody measurements and posterior median predictions, returning pointwise residuals or aggregated error metrics (MAE, RMSE, SSE) at various summary levels. --- R/compute_residual_metrics.R | 153 +++++++++++++++++++++++++++++++++++ 1 file changed, 153 insertions(+) create mode 100644 R/compute_residual_metrics.R diff --git a/R/compute_residual_metrics.R b/R/compute_residual_metrics.R new file mode 100644 index 0000000..0a83b9d --- /dev/null +++ b/R/compute_residual_metrics.R @@ -0,0 +1,153 @@ +#' Residual-based posterior predictive metrics for longitudinal antibody curves +#' +#' Computes residuals between observed antibody measurements and posterior median +#' predictions evaluated at the observed time points. Returns pointwise residuals +#' or aggregated error metrics (MAE, RMSE, SSE) at multiple summary levels. +#' +#' @param model A data frame of posterior draws in long format with columns: +#' \code{Subject}, \code{Iso_type}, \code{Chain}, \code{Iteration}, \code{Parameter}, \code{value}. +#' @param dataset A \code{serodynamics} case dataset produced by +#' \code{serodynamics::as_case_data()} (must contain \code{id} and \code{antigen_iso} +#' columns, and time/value attributes). +#' @param ids Character vector of subject IDs to include (matched against \code{dataset$id} +#' and \code{model$Subject}). +#' @param antigen_iso Character scalar specifying the antigen/isotype to analyze +#' (matched against \code{dataset$antigen_iso} and \code{model$Iso_type}). +#' @param scale Scale on which to compute residuals. One of \code{"original"} or \code{"log"}. +#' If \code{"log"}, residuals are computed on the natural log scale and observations/predictions +#' \eqn{\le 0} are removed. +#' @param summary_level Level at which to summarize metrics. One of: +#' \describe{ +#' \item{\code{"pointwise"}}{Return pointwise residuals at each observed time point.} +#' \item{\code{"id_antigen"}}{Summarize by \code{id} and \code{antigen_iso}.} +#' \item{\code{"antigen"}}{Summarize by \code{antigen_iso} only.} +#' \item{\code{"overall"}}{Single summary across all included observations.} +#' } +#' +#' @return A tibble. For \code{summary_level = "pointwise"}, returns per-observation residuals. +#' Otherwise returns MAE, RMSE, SSE, and \code{n_obs} at the requested summary level. +#' +#' @details +#' The posterior predictive summary uses the posterior median at each observed time point, +#' with a 95\% credible interval computed from draw-level predictions. Predictions are +#' generated by \code{\link{predict_posterior_at_times}}. +#' +#' @examples +#' \dontrun{ +#' # Per-ID error metrics (original scale) +#' m_id <- compute_residual_metrics( +#' model = overall_sf2a, +#' dataset = dL_clean_sf2a, +#' ids = unique(dL_clean_sf2a$id), +#' antigen_iso = "IgG", +#' scale = "original", +#' summary_level = "id_antigen" +#' ) +#' +#' # Pointwise residuals on log scale +#' r_pw <- compute_residual_metrics( +#' model = overall_sf2a, +#' dataset = dL_clean_sf2a, +#' ids = "SOSAR-22008", +#' antigen_iso = "IgA", +#' scale = "log", +#' summary_level = "pointwise" +#' ) +#' } +#' +#' @export +compute_residual_metrics <- function(model, + dataset, + ids, + antigen_iso, + scale = c("original", "log"), + summary_level = c("id_antigen", "pointwise", "antigen", "overall")) { + + scale <- match.arg(scale) + summary_level <- match.arg(summary_level) + + time_var <- serodynamics:::get_timeindays_var(dataset) + value_var <- serocalculator::get_values_var(dataset) + + observed_data <- dataset |> + dplyr::rename(t = {{ time_var }}, obs = {{ value_var }}) |> + dplyr::select(dplyr::all_of(c("id", "t", "obs", "antigen_iso"))) |> + dplyr::mutate(id = as.character(.data$id)) |> + dplyr::filter(.data$id %in% ids, .data$antigen_iso == antigen_iso) + + if (nrow(observed_data) == 0) { + rlang::abort("No observed data found for the specified IDs and antigen_iso.") + } + + obs_times <- sort(unique(observed_data$t)) + + predictions_all <- predict_posterior_at_times( + model = model, ids = ids, antigen_iso = antigen_iso, times = obs_times + ) + + pred_summary <- predictions_all |> + dplyr::summarise( + .by = c("id", "t"), + pred_med = stats::median(.data$res, na.rm = TRUE), + pred_lower = stats::quantile(.data$res, probs = 0.025, na.rm = TRUE), + pred_upper = stats::quantile(.data$res, probs = 0.975, na.rm = TRUE) + ) |> + dplyr::mutate(id = as.character(.data$id)) + + residual_data <- observed_data |> + dplyr::inner_join(pred_summary, by = c("id", "t")) + + if (scale == "log") { + residual_data <- residual_data |> + dplyr::filter(.data$obs > 0, .data$pred_med > 0) |> + dplyr::mutate( + residual = log(.data$obs) - log(.data$pred_med), + abs_residual = abs(.data$residual), + sq_residual = .data$residual^2 + ) + } else { + residual_data <- residual_data |> + dplyr::mutate( + residual = .data$obs - .data$pred_med, + abs_residual = abs(.data$residual), + sq_residual = .data$residual^2 + ) + } + + if (summary_level == "pointwise") { + return(residual_data |> + dplyr::select("id", "antigen_iso", "t", "obs", "pred_med", + "pred_lower", "pred_upper", + "residual", "abs_residual", "sq_residual")) + } + + if (summary_level == "id_antigen") { + return(residual_data |> + dplyr::summarise( + .by = c("id", "antigen_iso"), + MAE = mean(.data$abs_residual, na.rm = TRUE), + RMSE = sqrt(mean(.data$sq_residual, na.rm = TRUE)), + SSE = sum(.data$sq_residual, na.rm = TRUE), + n_obs = dplyr::n() + )) + } + + if (summary_level == "antigen") { + return(residual_data |> + dplyr::summarise( + .by = "antigen_iso", + MAE = mean(.data$abs_residual, na.rm = TRUE), + RMSE = sqrt(mean(.data$sq_residual, na.rm = TRUE)), + SSE = sum(.data$sq_residual, na.rm = TRUE), + n_obs = dplyr::n() + )) + } + + residual_data |> + dplyr::summarise( + MAE = mean(.data$abs_residual, na.rm = TRUE), + RMSE = sqrt(mean(.data$sq_residual, na.rm = TRUE)), + SSE = sum(.data$sq_residual, na.rm = TRUE), + n_obs = dplyr::n() + ) +} From 9efd9b13a84ce8f3e076e90cc6b9151927539dcc Mon Sep 17 00:00:00 2001 From: Kwan-Jenny <68584166+Kwan-Jenny@users.noreply.github.com> Date: Sat, 31 Jan 2026 22:35:00 -0800 Subject: [PATCH 14/30] Add model comparison function with residual metrics This function computes and compares residual metrics for serotype-specific and overall models, providing insights on their performance differences. --- R/model_comparison_table.R | 82 ++++++++++++++++++++++++++++++++++++++ 1 file changed, 82 insertions(+) create mode 100644 R/model_comparison_table.R diff --git a/R/model_comparison_table.R b/R/model_comparison_table.R new file mode 100644 index 0000000..01a12c1 --- /dev/null +++ b/R/model_comparison_table.R @@ -0,0 +1,82 @@ +#' Compare serotype-specific vs overall models using residual metrics +#' +#' Computes per-ID residual metrics for two models on the intersection of IDs +#' present in both datasets, then reports absolute and percent differences. +#' +#' @param model_serospecific Posterior draws (long format) for the serotype-specific model. +#' @param data_serospecific Case dataset used for the serotype-specific model. +#' @param model_overall Posterior draws (long format) for the overall model. +#' @param data_overall Case dataset used for the overall model. +#' @param antigen_iso Character scalar antigen/isotype label. +#' @param scale "original" or "log". +#' @param tie_tol Numeric tolerance to declare a tie. +#' +#' @return A tibble with per-ID MAE/RMSE for each model and deltas/winner labels. +#' +#' @export +make_model_comparison_table <- function(model_serospecific, data_serospecific, + model_overall, data_overall, + antigen_iso, + scale = c("original", "log"), + tie_tol = 1e-8) { + scale <- match.arg(scale) + + ids_common <- intersect(unique(data_serospecific$id), unique(data_overall$id)) + + m_sero <- compute_residual_metrics( + model = model_serospecific, + dataset = data_serospecific, + ids = ids_common, + antigen_iso = antigen_iso, + scale = scale, + summary_level = "id_antigen" + ) |> + dplyr::select(.data$id, .data$antigen_iso, .data$MAE, .data$RMSE, .data$n_obs) |> + dplyr::rename( + MAE_serospecific = .data$MAE, + RMSE_serospecific = .data$RMSE, + n_obs_serospecific = .data$n_obs + ) + + m_over <- compute_residual_metrics( + model = model_overall, + dataset = data_overall, + ids = ids_common, + antigen_iso = antigen_iso, + scale = scale, + summary_level = "id_antigen" + ) |> + dplyr::select(.data$id, .data$antigen_iso, .data$MAE, .data$RMSE, .data$n_obs) |> + dplyr::rename( + MAE_overall = .data$MAE, + RMSE_overall = .data$RMSE, + n_obs_overall = .data$n_obs + ) + + dplyr::full_join(m_sero, m_over, by = c("id", "antigen_iso")) |> + dplyr::mutate( + delta_MAE = .data$MAE_overall - .data$MAE_serospecific, + delta_RMSE = .data$RMSE_overall - .data$RMSE_serospecific, + pct_improve_MAE = 100 * .data$delta_MAE / .data$MAE_overall, + pct_improve_RMSE = 100 * .data$delta_RMSE / .data$RMSE_overall, + best_MAE = dplyr::case_when( + is.na(.data$MAE_overall) | is.na(.data$MAE_serospecific) ~ NA_character_, + abs(.data$delta_MAE) <= tie_tol ~ "tie", + .data$delta_MAE > 0 ~ "serospecific", + TRUE ~ "overall" + ), + best_RMSE = dplyr::case_when( + is.na(.data$RMSE_overall) | is.na(.data$RMSE_serospecific) ~ NA_character_, + abs(.data$delta_RMSE) <= tie_tol ~ "tie", + .data$delta_RMSE > 0 ~ "serospecific", + TRUE ~ "overall" + ), + best_overall = dplyr::case_when( + .data$best_MAE == "serospecific" & .data$best_RMSE == "serospecific" ~ "serospecific", + .data$best_MAE == "overall" & .data$best_RMSE == "overall" ~ "overall", + .data$best_MAE == "tie" & .data$best_RMSE == "tie" ~ "tie", + TRUE ~ "mixed" + ) + ) |> + dplyr::arrange(.data$id) +} From 4655d41e2b46a2eaa57470f2b057e5c335d9d6a8 Mon Sep 17 00:00:00 2001 From: Kwan-Jenny <68584166+Kwan-Jenny@users.noreply.github.com> Date: Sat, 31 Jan 2026 22:35:49 -0800 Subject: [PATCH 15/30] Add fmt_mci function to format median and intervals --- R/fmt_mci.R | 17 +++++++++++++++++ 1 file changed, 17 insertions(+) create mode 100644 R/fmt_mci.R diff --git a/R/fmt_mci.R b/R/fmt_mci.R new file mode 100644 index 0000000..9982358 --- /dev/null +++ b/R/fmt_mci.R @@ -0,0 +1,17 @@ +#' Format median and credible interval as a single string +#' +#' @param med Median value. +#' @param lo Lower bound. +#' @param hi Upper bound. +#' @param digits Number of digits. +#' @param sci Logical; if TRUE uses scientific notation. +#' +#' @return A string like "1.23 (0.50ā€“2.00)". +#' +#' @export +fmt_mci <- function(med, lo, hi, digits = 2, sci = FALSE) { + f <- function(x) { + if (sci) formatC(x, format = "e", digits = digits) else formatC(x, format = "f", digits = digits) + } + sprintf("%s (%sā€“%s)", f(med), f(lo), f(hi)) +} From ae38681c4c9f838b156e8f53b20e44119b5a0560 Mon Sep 17 00:00:00 2001 From: Kwan-Jenny <68584166+Kwan-Jenny@users.noreply.github.com> Date: Sat, 31 Jan 2026 22:36:20 -0800 Subject: [PATCH 16/30] Add function to extract and summarize newperson parameters This function processes posterior draws for the 'newperson' subject, filtering and reshaping the data into a wide format suitable for analysis. --- R/prep_newperson_params.R | 22 ++++++++++++++++++++++ 1 file changed, 22 insertions(+) create mode 100644 R/prep_newperson_params.R diff --git a/R/prep_newperson_params.R b/R/prep_newperson_params.R new file mode 100644 index 0000000..2c4b5c1 --- /dev/null +++ b/R/prep_newperson_params.R @@ -0,0 +1,22 @@ +#' Extract "newperson" parameter draws and summarize for Table 2 +#' +#' @param draws_long Posterior draws in long format with columns: +#' Subject, Iso_type, Chain, Iteration, Parameter, value. +#' @param antigen_label Character scalar antigen name to attach. +#' +#' @return A tibble of newperson draws in wide format with columns: +#' Iteration, Chain, antigen, Iso_type, y0, y1, t1, alpha, rho. +#' +#' @export +prep_newperson_params <- function(draws_long, antigen_label) { + draws_long |> + dplyr::filter(.data$Subject == "newperson") |> + dplyr::mutate(antigen = antigen_label) |> + dplyr::select(.data$Iteration, .data$Chain, .data$antigen, + .data$Iso_type, .data$Parameter, .data$value) |> + tidyr::pivot_wider(names_from = .data$Parameter, values_from = .data$value) |> + # In your model output, decay-shape parameter is named "shape" + dplyr::rename(rho = .data$shape) |> + dplyr::select(.data$Iteration, .data$Chain, .data$antigen, .data$Iso_type, + .data$y0, .data$y1, .data$t1, .data$alpha, .data$rho) +} From 87c45a27e4f6e549dfab0e20fb071ceca9306acc Mon Sep 17 00:00:00 2001 From: Kwan-Jenny <68584166+Kwan-Jenny@users.noreply.github.com> Date: Sat, 31 Jan 2026 22:58:51 -0800 Subject: [PATCH 17/30] Update R/fig2_overall_newperson.R Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- R/fig2_overall_newperson.R | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/R/fig2_overall_newperson.R b/R/fig2_overall_newperson.R index eccba0b..81af6c8 100644 --- a/R/fig2_overall_newperson.R +++ b/R/fig2_overall_newperson.R @@ -73,11 +73,12 @@ fig2_overall_newperson <- function( ) } - sum_all <- purrr::map_dfr(osps, function(osp) { - purrr::map_dfr(isotypes, function(iso) { + sum_all <- lapply(osps, function(osp) { + dplyr::bind_rows(lapply(isotypes, function(iso) { get_sum_overall(overall_models[[osp]], osp, iso) - }) + })) }) |> + dplyr::bind_rows() |> dplyr::mutate(antigen = factor(.data$antigen, levels = osps)) p <- ggplot2::ggplot(sum_all, ggplot2::aes(x = .data$t)) + From e7648cdab6a7c0e514026e52681eac99669a82e3 Mon Sep 17 00:00:00 2001 From: Kwan-Jenny <68584166+Kwan-Jenny@users.noreply.github.com> Date: Sat, 31 Jan 2026 22:59:04 -0800 Subject: [PATCH 18/30] Update R/prep_newperson_params.R Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- R/prep_newperson_params.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/prep_newperson_params.R b/R/prep_newperson_params.R index 2c4b5c1..7371d58 100644 --- a/R/prep_newperson_params.R +++ b/R/prep_newperson_params.R @@ -15,7 +15,7 @@ prep_newperson_params <- function(draws_long, antigen_label) { dplyr::select(.data$Iteration, .data$Chain, .data$antigen, .data$Iso_type, .data$Parameter, .data$value) |> tidyr::pivot_wider(names_from = .data$Parameter, values_from = .data$value) |> - # In your model output, decay-shape parameter is named "shape" + # The model output uses "shape" for the decay-shape parameter, which is renamed to "rho" for consistency dplyr::rename(rho = .data$shape) |> dplyr::select(.data$Iteration, .data$Chain, .data$antigen, .data$Iso_type, .data$y0, .data$y1, .data$t1, .data$alpha, .data$rho) From 3702aa10b8119e810fe2861053491f6ba74c4d62 Mon Sep 17 00:00:00 2001 From: Kwan-Jenny <68584166+Kwan-Jenny@users.noreply.github.com> Date: Sat, 31 Jan 2026 22:59:21 -0800 Subject: [PATCH 19/30] Update R/model_comparison_table.R Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- R/model_comparison_table.R | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/R/model_comparison_table.R b/R/model_comparison_table.R index 01a12c1..df88e8a 100644 --- a/R/model_comparison_table.R +++ b/R/model_comparison_table.R @@ -57,8 +57,16 @@ make_model_comparison_table <- function(model_serospecific, data_serospecific, dplyr::mutate( delta_MAE = .data$MAE_overall - .data$MAE_serospecific, delta_RMSE = .data$RMSE_overall - .data$RMSE_serospecific, - pct_improve_MAE = 100 * .data$delta_MAE / .data$MAE_overall, - pct_improve_RMSE = 100 * .data$delta_RMSE / .data$RMSE_overall, + pct_improve_MAE = dplyr::case_when( + is.na(.data$MAE_overall) ~ NA_real_, + abs(.data$MAE_overall) <= .Machine$double.eps ~ NA_real_, + TRUE ~ 100 * .data$delta_MAE / .data$MAE_overall + ), + pct_improve_RMSE = dplyr::case_when( + is.na(.data$RMSE_overall) ~ NA_real_, + abs(.data$RMSE_overall) <= .Machine$double.eps ~ NA_real_, + TRUE ~ 100 * .data$delta_RMSE / .data$RMSE_overall + ), best_MAE = dplyr::case_when( is.na(.data$MAE_overall) | is.na(.data$MAE_serospecific) ~ NA_character_, abs(.data$delta_MAE) <= tie_tol ~ "tie", From be6012c894c3b9f97424b8770f05704c19f52fb9 Mon Sep 17 00:00:00 2001 From: Kwan-Jenny <68584166+Kwan-Jenny@users.noreply.github.com> Date: Sat, 31 Jan 2026 22:59:33 -0800 Subject: [PATCH 20/30] Update R/compute_residual_metrics.R Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- R/compute_residual_metrics.R | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/R/compute_residual_metrics.R b/R/compute_residual_metrics.R index 0a83b9d..550286f 100644 --- a/R/compute_residual_metrics.R +++ b/R/compute_residual_metrics.R @@ -76,7 +76,7 @@ compute_residual_metrics <- function(model, dplyr::filter(.data$id %in% ids, .data$antigen_iso == antigen_iso) if (nrow(observed_data) == 0) { - rlang::abort("No observed data found for the specified IDs and antigen_iso.") + stop("No observed data found for the specified IDs and antigen_iso.") } obs_times <- sort(unique(observed_data$t)) From cca53928d6dcb9c2750cfabab97eaf4cd4a6b494 Mon Sep 17 00:00:00 2001 From: Kwan-Jenny <68584166+Kwan-Jenny@users.noreply.github.com> Date: Sat, 31 Jan 2026 22:59:39 -0800 Subject: [PATCH 21/30] Update R/process_shigella_data.R Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- R/process_shigella_data.R | 1 + 1 file changed, 1 insertion(+) diff --git a/R/process_shigella_data.R b/R/process_shigella_data.R index 215222f..a600277 100644 --- a/R/process_shigella_data.R +++ b/R/process_shigella_data.R @@ -26,6 +26,7 @@ #' ) #' } #' +#' @importFrom rlang ensym #' @export process_shigella_data <- function(data, study_filter, antigen) { antigen_col <- rlang::ensym(antigen) From 9b608fe4dc8f3b91f1beabcd6cd2b01c3d4f845f Mon Sep 17 00:00:00 2001 From: Copilot <198982749+Copilot@users.noreply.github.com> Date: Mon, 30 Mar 2026 17:37:18 -0700 Subject: [PATCH 22/30] Fix package dependencies, namespace exports, and package infrastructure (#12) * Initial plan * Fix DESCRIPTION dependencies and remove ::: usage - Move tibble, tidyr, purrr from Suggests to Imports - Add serocalculator and rlang to Imports - Remove unused cohort_name selection in process_shigella_data - Replace serodynamics::: calls with internal helper functions - Add utils_internal.R with get_timeindays_var() and ab() implementations Co-authored-by: Kwan-Jenny <68584166+Kwan-Jenny@users.noreply.github.com> * Add package infrastructure: tests, workflows, vignettes, and mock data - Add testthat tests for main functions - Add R CMD check and pkgdown GitHub Actions workflows - Add data-raw/ with mock data generation script - Add R/data.R with mock data documentation - Add getting-started vignette - Update README.Rmd with package overview and usage - Update .Rbuildignore to exclude non-package files Co-authored-by: Kwan-Jenny <68584166+Kwan-Jenny@users.noreply.github.com> * Fix rlang::rnorm() to rnorm() in mock data generation Use base R rnorm() instead of non-existent rlang::rnorm() Co-authored-by: Kwan-Jenny <68584166+Kwan-Jenny@users.noreply.github.com> * Add permissions block to R-CMD-check workflow Set explicit GITHUB_TOKEN permissions (contents: read) for security Co-authored-by: Kwan-Jenny <68584166+Kwan-Jenny@users.noreply.github.com> * Update documentation to remove reference to serodynamics:::ab() Update comment to reflect internal implementation instead of external ::: call Co-authored-by: Kwan-Jenny <68584166+Kwan-Jenny@users.noreply.github.com> * Align internal utils with serodynamics and update tests Co-authored-by: Kwan-Jenny <68584166+Kwan-Jenny@users.noreply.github.com> Agent-Logs-Url: https://github.com/UCD-SERG/shigella/sessions/d36469e9-e5db-41bb-ad8f-cc452230b453 * Fix namespace export errors and add mock data .rda files Agent-Logs-Url: https://github.com/UCD-SERG/shigella/sessions/1327221c-702c-4c06-907f-1b381d75ecd9 Co-authored-by: Kwan-Jenny <68584166+Kwan-Jenny@users.noreply.github.com> * Fix failing tests by adding model_comparison_table export and sync roxygen note Agent-Logs-Url: https://github.com/UCD-SERG/shigella/sessions/47d8cc5c-e58f-4f67-a13f-b5f56651a62e Co-authored-by: Kwan-Jenny <68584166+Kwan-Jenny@users.noreply.github.com> * Fix remaining docs-check and lint CI issues Agent-Logs-Url: https://github.com/UCD-SERG/shigella/sessions/d635bda4-41b7-4feb-8fc3-330a60fec1ac Co-authored-by: Kwan-Jenny <68584166+Kwan-Jenny@users.noreply.github.com> * Fix compute_residual_metrics indentation at line 137 Agent-Logs-Url: https://github.com/UCD-SERG/shigella/sessions/6af850e8-229a-4ffd-a552-ff9467c5bbbb Co-authored-by: Kwan-Jenny <68584166+Kwan-Jenny@users.noreply.github.com> --------- Co-authored-by: copilot-swe-agent[bot] <198982749+Copilot@users.noreply.github.com> Co-authored-by: Kwan-Jenny <68584166+Kwan-Jenny@users.noreply.github.com> --- .Rbuildignore | 9 ++ .github/workflows/R-CMD-check.yaml | 48 +++++++ .github/workflows/pkgdown.yaml | 46 ++++++ .gitignore | 5 +- DESCRIPTION | 13 +- NAMESPACE | 10 +- R/compute_residual_metrics.R | 77 ++++++---- R/data.R | 65 +++++++++ R/fig2_overall_newperson.R | 41 ++++-- R/model_comparison_table.R | 98 +++++++++++-- R/predict_posterior_at_times.R | 41 ++++-- R/process_shigella_data.R | 15 +- R/shigella-package.R | 1 + R/utils_internal.R | 51 +++++++ README.Rmd | 59 +++++++- data-raw/mock_data.R | 85 +++++++++++ data/.gitignore | 2 + data/mock_case_data.rda | Bin 0 -> 984 bytes data/mock_posterior_draws.rda | Bin 0 -> 56194 bytes man/compute_residual_metrics.Rd | 87 ++++++++++++ man/fig2_overall_newperson.Rd | 69 +++++++++ man/fmt_mci.Rd | 25 ++++ man/make_model_comparison_table.Rd | 40 ++++++ man/mock_case_data.Rd | 39 +++++ man/mock_posterior_draws.Rd | 43 ++++++ man/model_comparison_table.Rd | 35 +++++ man/postprocess_jags_output.Rd | 17 --- man/predict_posterior_at_times.Rd | 64 +++++++++ man/prep_newperson_params.Rd | 21 +++ man/process_shigella_data.Rd | 42 ++++++ man/utils_internal.Rd | 9 ++ tests/testthat.R | 12 ++ .../testthat/test-compute_residual_metrics.R | 102 +++++++++++++ tests/testthat/test-model_comparison_table.R | 68 +++++++++ .../test-predict_posterior_at_times.R | 68 +++++++++ tests/testthat/test-process_shigella_data.R | 46 ++++++ tests/testthat/test-utils_internal.R | 45 ++++++ vignettes/getting-started.Rmd | 134 ++++++++++++++++++ 38 files changed, 1541 insertions(+), 91 deletions(-) create mode 100644 .github/workflows/R-CMD-check.yaml create mode 100644 .github/workflows/pkgdown.yaml create mode 100644 R/data.R create mode 100644 R/utils_internal.R create mode 100644 data-raw/mock_data.R create mode 100644 data/mock_case_data.rda create mode 100644 data/mock_posterior_draws.rda create mode 100644 man/compute_residual_metrics.Rd create mode 100644 man/fig2_overall_newperson.Rd create mode 100644 man/fmt_mci.Rd create mode 100644 man/make_model_comparison_table.Rd create mode 100644 man/mock_case_data.Rd create mode 100644 man/mock_posterior_draws.Rd create mode 100644 man/model_comparison_table.Rd delete mode 100644 man/postprocess_jags_output.Rd create mode 100644 man/predict_posterior_at_times.Rd create mode 100644 man/prep_newperson_params.Rd create mode 100644 man/process_shigella_data.Rd create mode 100644 man/utils_internal.Rd create mode 100644 tests/testthat.R create mode 100644 tests/testthat/test-compute_residual_metrics.R create mode 100644 tests/testthat/test-model_comparison_table.R create mode 100644 tests/testthat/test-predict_posterior_at_times.R create mode 100644 tests/testthat/test-process_shigella_data.R create mode 100644 tests/testthat/test-utils_internal.R create mode 100644 vignettes/getting-started.Rmd diff --git a/.Rbuildignore b/.Rbuildignore index a0f6a52..d92b3e9 100644 --- a/.Rbuildignore +++ b/.Rbuildignore @@ -8,3 +8,12 @@ ^vignettes/articles/*_files$ ^vignettes/articles$ ^codecov\.yml$ +^data-raw$ +^README\.Rmd$ +^\.lintr\.R$ +^_quarto\.yml$ +^.*\.qmd$ +^.*\.png$ +^Vary.*\.r$ +^sim_.*\.r$ +^simulation.*\.qmd$ diff --git a/.github/workflows/R-CMD-check.yaml b/.github/workflows/R-CMD-check.yaml new file mode 100644 index 0000000..f01eb71 --- /dev/null +++ b/.github/workflows/R-CMD-check.yaml @@ -0,0 +1,48 @@ +on: + push: + branches: [main, master] + pull_request: + branches: [main, master] + +name: R-CMD-check + +permissions: + contents: read + +jobs: + R-CMD-check: + runs-on: ${{ matrix.config.os }} + + name: ${{ matrix.config.os }} (${{ matrix.config.r }}) + + strategy: + fail-fast: false + matrix: + config: + - {os: ubuntu-latest, r: 'release'} + - {os: ubuntu-latest, r: 'devel', http-user-agent: 'release'} + + env: + GITHUB_PAT: ${{ github.token }} + R_KEEP_PKG_SOURCE: yes + + 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-r-dependencies@v2 + with: + extra-packages: any::rcmdcheck + needs: check + + - uses: r-lib/actions/check-r-package@v2 + with: + upload-snapshots: true + build_args: 'c("--no-manual","--compact-vignettes=gs+qpdf")' diff --git a/.github/workflows/pkgdown.yaml b/.github/workflows/pkgdown.yaml new file mode 100644 index 0000000..f238786 --- /dev/null +++ b/.github/workflows/pkgdown.yaml @@ -0,0 +1,46 @@ +on: + push: + branches: [main, master] + release: + types: [published] + workflow_dispatch: + +name: pkgdown + +permissions: + contents: write + pages: write + +jobs: + pkgdown: + runs-on: ubuntu-latest + # Only restrict concurrency for non-PR jobs + concurrency: + group: pkgdown-${{ github.event_name != 'pull_request' || github.run_id }} + env: + GITHUB_PAT: ${{ github.token }} + steps: + - uses: actions/checkout@v4 + + - uses: r-lib/actions/setup-pandoc@v2 + + - uses: r-lib/actions/setup-r@v2 + with: + use-public-rspm: true + + - uses: r-lib/actions/setup-r-dependencies@v2 + with: + extra-packages: any::pkgdown, local::. + needs: website + + - name: Build site + run: pkgdown::build_site_github_pages(new_process = FALSE, install = FALSE) + shell: Rscript {0} + + - name: Deploy to GitHub pages šŸš€ + if: github.event_name != 'pull_request' + uses: JamesIves/github-pages-deploy-action@v4 + with: + clean: false + branch: gh-pages + folder: docs diff --git a/.gitignore b/.gitignore index cbe14e1..5b09519 100644 --- a/.gitignore +++ b/.gitignore @@ -8,7 +8,10 @@ simulation1_edited_files simulation1_edited.html *.rmarkdown **/.quarto/ -data/ +data/* +!data/.gitignore +!data/mock_posterior_draws.rda +!data/mock_case_data.rda /.quarto/ *.pdf README_files diff --git a/DESCRIPTION b/DESCRIPTION index 6b43288..50e8fa8 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -9,14 +9,19 @@ Description: What the package does (one paragraph). License: MIT + file LICENSE Encoding: UTF-8 Roxygen: list(markdown = TRUE) -RoxygenNote: 7.3.2 +RoxygenNote: 7.3.3 Imports: + cli, dplyr, future, future.apply, ggplot2, gridExtra, - serodynamics (>= 0.0.0.9011) + rlang, + serocalculator, + serodynamics (>= 0.0.0.9011), + tibble, + tidyr URL: https://ucd-serg.github.io/shigella/ Remotes: UCD-SERG/serodynamics @@ -46,11 +51,11 @@ Suggests: scales, spelling, table1, - tibble, - tidyr, + testthat (>= 3.0.0), tidyverse Depends: R (>= 3.5) LazyData: true +Config/testthat/edition: 3 Config/Needs/website: quarto Language: en-US diff --git a/NAMESPACE b/NAMESPACE index 56e7b54..0d108b7 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,6 +1,12 @@ # Generated by roxygen2: do not edit by hand -export(postprocess_jags_output) +export(compute_residual_metrics) +export(fig2_overall_newperson) +export(fmt_mci) +export(make_model_comparison_table) +export(model_comparison_table) +export(prep_newperson_params) +export(process_shigella_data) importFrom(dplyr,any_of) importFrom(dplyr,filter) importFrom(dplyr,mutate) @@ -13,3 +19,5 @@ importFrom(ggplot2,geom_point) importFrom(ggplot2,ggplot) importFrom(ggplot2,labs) importFrom(ggplot2,theme_minimal) +importFrom(rlang,.data) +importFrom(rlang,ensym) diff --git a/R/compute_residual_metrics.R b/R/compute_residual_metrics.R index 550286f..815721e 100644 --- a/R/compute_residual_metrics.R +++ b/R/compute_residual_metrics.R @@ -1,36 +1,44 @@ #' Residual-based posterior predictive metrics for longitudinal antibody curves #' -#' Computes residuals between observed antibody measurements and posterior median -#' predictions evaluated at the observed time points. Returns pointwise residuals -#' or aggregated error metrics (MAE, RMSE, SSE) at multiple summary levels. +#' Computes residuals between observed antibody measurements and posterior +#' median predictions evaluated at the observed time points. Returns pointwise +#' residuals or aggregated error metrics (MAE, RMSE, SSE) at multiple +#' summary levels. #' #' @param model A data frame of posterior draws in long format with columns: -#' \code{Subject}, \code{Iso_type}, \code{Chain}, \code{Iteration}, \code{Parameter}, \code{value}. +#' \code{Subject}, \code{Iso_type}, \code{Chain}, \code{Iteration}, +#' \code{Parameter}, \code{value}. #' @param dataset A \code{serodynamics} case dataset produced by -#' \code{serodynamics::as_case_data()} (must contain \code{id} and \code{antigen_iso} -#' columns, and time/value attributes). -#' @param ids Character vector of subject IDs to include (matched against \code{dataset$id} -#' and \code{model$Subject}). +#' \code{serodynamics::as_case_data()} (must contain \code{id} and +#' \code{antigen_iso} columns, and time/value attributes). +#' @param ids Character vector of subject IDs to include (matched against +#' \code{dataset$id} and \code{model$Subject}). #' @param antigen_iso Character scalar specifying the antigen/isotype to analyze #' (matched against \code{dataset$antigen_iso} and \code{model$Iso_type}). -#' @param scale Scale on which to compute residuals. One of \code{"original"} or \code{"log"}. -#' If \code{"log"}, residuals are computed on the natural log scale and observations/predictions -#' \eqn{\le 0} are removed. +#' @param scale Scale on which to compute residuals. One of \code{"original"} +#' or \code{"log"}. If \code{"log"}, residuals are computed on the natural +#' log scale and observations/predictions \eqn{\le 0} are removed. #' @param summary_level Level at which to summarize metrics. One of: #' \describe{ -#' \item{\code{"pointwise"}}{Return pointwise residuals at each observed time point.} -#' \item{\code{"id_antigen"}}{Summarize by \code{id} and \code{antigen_iso}.} +#' \item{\code{"pointwise"}}{ +#' Return pointwise residuals at each observed time point. +#' } +#' \item{\code{"id_antigen"}}{ +#' Summarize by \code{id} and \code{antigen_iso}. +#' } #' \item{\code{"antigen"}}{Summarize by \code{antigen_iso} only.} #' \item{\code{"overall"}}{Single summary across all included observations.} #' } #' -#' @return A tibble. For \code{summary_level = "pointwise"}, returns per-observation residuals. -#' Otherwise returns MAE, RMSE, SSE, and \code{n_obs} at the requested summary level. +#' @return A tibble. For \code{summary_level = "pointwise"}, returns +#' per-observation residuals. Otherwise returns MAE, RMSE, SSE, and +#' \code{n_obs} at the requested summary level. #' #' @details -#' The posterior predictive summary uses the posterior median at each observed time point, -#' with a 95\% credible interval computed from draw-level predictions. Predictions are -#' generated by \code{\link{predict_posterior_at_times}}. +#' The posterior predictive summary uses the posterior median at each observed +#' time point, with a 95\% credible interval computed from draw-level +#' predictions. Predictions are generated by +#' \code{\link{predict_posterior_at_times}}. #' #' @examples #' \dontrun{ @@ -61,12 +69,17 @@ compute_residual_metrics <- function(model, ids, antigen_iso, scale = c("original", "log"), - summary_level = c("id_antigen", "pointwise", "antigen", "overall")) { + summary_level = c( + "id_antigen", + "pointwise", + "antigen", + "overall" + )) { scale <- match.arg(scale) summary_level <- match.arg(summary_level) - time_var <- serodynamics:::get_timeindays_var(dataset) + time_var <- get_timeindays_var(dataset) value_var <- serocalculator::get_values_var(dataset) observed_data <- dataset |> @@ -76,13 +89,18 @@ compute_residual_metrics <- function(model, dplyr::filter(.data$id %in% ids, .data$antigen_iso == antigen_iso) if (nrow(observed_data) == 0) { - stop("No observed data found for the specified IDs and antigen_iso.") + cli::cli_abort( + "No observed data found for the specified IDs and antigen_iso." + ) } obs_times <- sort(unique(observed_data$t)) predictions_all <- predict_posterior_at_times( - model = model, ids = ids, antigen_iso = antigen_iso, times = obs_times + model = model, + ids = ids, + antigen_iso = antigen_iso, + times = obs_times ) pred_summary <- predictions_all |> @@ -116,9 +134,18 @@ compute_residual_metrics <- function(model, if (summary_level == "pointwise") { return(residual_data |> - dplyr::select("id", "antigen_iso", "t", "obs", "pred_med", - "pred_lower", "pred_upper", - "residual", "abs_residual", "sq_residual")) + dplyr::select( + "id", + "antigen_iso", + "t", + "obs", + "pred_med", + "pred_lower", + "pred_upper", + "residual", + "abs_residual", + "sq_residual" + )) } if (summary_level == "id_antigen") { diff --git a/R/data.R b/R/data.R new file mode 100644 index 0000000..dd553d6 --- /dev/null +++ b/R/data.R @@ -0,0 +1,65 @@ +# Mock data documentation + +#' Mock posterior draws for testing +#' +#' A mock dataset of posterior parameter draws in long format, mimicking the +#' structure expected by functions like \code{\link{compute_residual_metrics}} +#' and \code{\link{predict_posterior_at_times}}. +#' +#' @format A data frame with columns: +#' \describe{ +#' \item{Subject}{Character. Subject ID (e.g., "newperson", "SOSAR-22008")} +#' \item{Iso_type}{Character. Isotype (e.g., "IgG", "IgA")} +#' \item{Chain}{Integer. MCMC chain number} +#' \item{Iteration}{Integer. MCMC iteration number} +#' \item{Parameter}{Character. Parameter name (y0, y1, t1, alpha, shape)} +#' \item{value}{Numeric. Parameter value} +#' } +#' +#' @details +#' This is synthetic data generated for testing and examples. Real Shigella +#' posterior draws will be added to the package separately. +#' +#' Parameters represent: +#' \itemize{ +#' \item \code{y0}: Baseline antibody level +#' \item \code{y1}: Peak antibody level +#' \item \code{t1}: Time to peak (days) +#' \item \code{alpha}: Decay rate parameter +#' \item \code{shape}: Decay shape parameter (rho) +#' } +#' +#' @examples +#' head(mock_posterior_draws) +#' table(mock_posterior_draws$Parameter) +"mock_posterior_draws" + +#' Mock case data for testing +#' +#' A mock longitudinal antibody dataset compatible with +#' \code{serodynamics::as_case_data()}, for testing functions like +#' \code{\link{compute_residual_metrics}}. +#' +#' @format A data frame with class \code{c("case_data", "data.frame")} +#' and columns: +#' \describe{ +#' \item{id}{Character. Subject ID} +#' \item{antigen_iso}{Character. Isotype (e.g., "IgG", "IgA")} +#' \item{timepoint}{Numeric. Time in days since infection} +#' \item{value}{Numeric. Antibody measurement} +#' } +#' +#' @details +#' This is synthetic data generated for testing and examples. Real Shigella +#' case data will be added separately. +#' +#' The dataset has attributes: +#' \itemize{ +#' \item \code{attr(mock_case_data, "timeindays") = "timepoint"} +#' \item \code{attr(mock_case_data, "value_var") = "value"} +#' } +#' +#' @examples +#' head(mock_case_data) +#' table(mock_case_data$id) +"mock_case_data" diff --git a/R/fig2_overall_newperson.R b/R/fig2_overall_newperson.R index 81af6c8..2c2e35b 100644 --- a/R/fig2_overall_newperson.R +++ b/R/fig2_overall_newperson.R @@ -1,24 +1,34 @@ -#' Summarize population-level ("newperson") antibody trajectories from overall models +utils::globalVariables(c("antigen", "iso")) + +#' Summarize population-level ("newperson") antibody trajectories +#' from overall models #' #' Computes median and credible interval bands of the antibody trajectory for a -#' hypothetical new individual drawn from the population distribution ("newperson"). +#' hypothetical new individual drawn from the population distribution +#' ("newperson"). #' -#' @param overall_models Named list of posterior draws in long format (one per antigen), -#' with columns at least: Subject, Iso_type, Chain, Iteration, Parameter, value. -#' @param osps Character vector of antigen names (must match names in `overall_models`). +#' @param overall_models Named list of posterior draws in long format +#' (one per antigen), +#' with columns at least: Subject, Iso_type, Chain, Iteration, +#' Parameter, value. +#' @param osps Character vector of antigen names +#' (must match names in `overall_models`). #' @param ids Character vector of Subject IDs to include (default: "newperson"). -#' @param isotypes Character vector of isotypes (default: c("IgG","IgA")). +#' @param isotypes Character vector of isotypes +#' (default: c("IgG", "IgA")). #' @param t_grid Numeric vector of time points (days) to evaluate. #' @param cred Credible level (default 0.95). #' @param log_y Logical; if TRUE, applies log10 scale to y when returning plot. -#' @param xlim Optional numeric length-2 vector for x-axis limits when returning plot. +#' @param xlim Optional numeric length-2 vector for x-axis limits when +#' returning plot. #' @param ylab Y-axis label for plot. #' @param line_color Line color for plot. #' @param ribbon_alpha Alpha for credible ribbon in plot. #' @param facet_scales Passed to ggplot facet scales. #' @param return_data If TRUE, returns a list with `plot` and `data`. #' -#' @return By default, a ggplot. If `return_data = TRUE`, returns list(plot = p, data = df). +#' @return By default, a ggplot. If `return_data = TRUE`, returns +#' list(plot = p, data = df). #' #' @details #' Requires that each model draw can be pivoted to wide parameters including: @@ -27,9 +37,9 @@ #' @export fig2_overall_newperson <- function( overall_models, - osps = c("IpaB","Sf2a","Sf3a","Sf6","Sonnei"), + osps = c("IpaB", "Sf2a", "Sf3a", "Sf6", "Sonnei"), ids = "newperson", - isotypes = c("IgG","IgA"), + isotypes = c("IgG", "IgA"), t_grid = seq(0, 210, by = 5), cred = 0.95, log_y = TRUE, @@ -48,14 +58,17 @@ fig2_overall_newperson <- function( dplyr::filter(.data$Subject %in% ids, .data$Iso_type == iso) |> dplyr::select(.data$Chain, .data$Iteration, .data$Iso_type, .data$Parameter, .data$value, .data$Subject) |> - tidyr::pivot_wider(names_from = .data$Parameter, values_from = .data$value) |> + tidyr::pivot_wider( + names_from = .data$Parameter, + values_from = .data$value + ) |> dplyr::mutate( antigen = osp, - iso = factor(.data$Iso_type, levels = c("IgG","IgA")) + iso = factor(.data$Iso_type, levels = c("IgG", "IgA")) ) |> tidyr::crossing(t = t_grid) |> dplyr::mutate( - res = serodynamics:::ab( + res = ab( t = .data$t, y0 = .data$y0, y1 = .data$y1, @@ -69,7 +82,7 @@ fig2_overall_newperson <- function( res.med = stats::quantile(.data$res, 0.50, na.rm = TRUE), res.low = stats::quantile(.data$res, q_lo, na.rm = TRUE), res.high = stats::quantile(.data$res, q_hi, na.rm = TRUE), - .by = c("antigen","iso","t") + .by = c("antigen", "iso", "t") ) } diff --git a/R/model_comparison_table.R b/R/model_comparison_table.R index df88e8a..61a6cfd 100644 --- a/R/model_comparison_table.R +++ b/R/model_comparison_table.R @@ -1,9 +1,71 @@ +#' Compare summary residual metrics between two model outputs +#' +#' Builds a compact comparison table from pre-computed summary metrics for an +#' overall model and a pointwise model. +#' +#' @param metrics_overall Data frame with one row containing at least +#' \code{MAE}, \code{RMSE}, \code{SSE}, and \code{n_obs} +#' for the overall model. +#' @param metrics_pointwise Data frame with one row containing at least +#' \code{MAE}, \code{RMSE}, \code{SSE}, and \code{n_obs} +#' for the pointwise model. +#' @param model_overall_label Label used for the overall model row. +#' @param model_pointwise_label Label used for the pointwise model row. +#' +#' @return A tibble with one row per model and comparison columns: +#' \code{delta_MAE}, \code{delta_RMSE}, \code{pct_improve_MAE}, +#' \code{pct_improve_RMSE}. +#' +#' @export +model_comparison_table <- function(metrics_overall, + metrics_pointwise, + model_overall_label = "Overall Model", + model_pointwise_label = "Pointwise Model") { + out <- dplyr::bind_rows( + dplyr::mutate( + metrics_overall, + Model = model_overall_label + ), + dplyr::mutate( + metrics_pointwise, + Model = model_pointwise_label + ) + ) |> + dplyr::select( + "Model", + "MAE", + "RMSE", + "SSE", + "n_obs" + ) + + base_mae <- metrics_overall$MAE[[1]] + base_rmse <- metrics_overall$RMSE[[1]] + + out |> + dplyr::mutate( + delta_MAE = .data$MAE - base_mae, + delta_RMSE = .data$RMSE - base_rmse, + pct_improve_MAE = dplyr::case_when( + is.na(base_mae) ~ NA_real_, + abs(base_mae) <= .Machine$double.eps ~ NA_real_, + TRUE ~ 100 * .data$delta_MAE / base_mae + ), + pct_improve_RMSE = dplyr::case_when( + is.na(base_rmse) ~ NA_real_, + abs(base_rmse) <= .Machine$double.eps ~ NA_real_, + TRUE ~ 100 * .data$delta_RMSE / base_rmse + ) + ) +} + #' Compare serotype-specific vs overall models using residual metrics #' #' Computes per-ID residual metrics for two models on the intersection of IDs #' present in both datasets, then reports absolute and percent differences. #' -#' @param model_serospecific Posterior draws (long format) for the serotype-specific model. +#' @param model_serospecific Posterior draws (long format) +#' for the serotype-specific model. #' @param data_serospecific Case dataset used for the serotype-specific model. #' @param model_overall Posterior draws (long format) for the overall model. #' @param data_overall Case dataset used for the overall model. @@ -11,7 +73,8 @@ #' @param scale "original" or "log". #' @param tie_tol Numeric tolerance to declare a tie. #' -#' @return A tibble with per-ID MAE/RMSE for each model and deltas/winner labels. +#' @return A tibble with per-ID MAE/RMSE for each model and +#' deltas/winner labels. #' #' @export make_model_comparison_table <- function(model_serospecific, data_serospecific, @@ -31,7 +94,13 @@ make_model_comparison_table <- function(model_serospecific, data_serospecific, scale = scale, summary_level = "id_antigen" ) |> - dplyr::select(.data$id, .data$antigen_iso, .data$MAE, .data$RMSE, .data$n_obs) |> + dplyr::select( + .data$id, + .data$antigen_iso, + .data$MAE, + .data$RMSE, + .data$n_obs + ) |> dplyr::rename( MAE_serospecific = .data$MAE, RMSE_serospecific = .data$RMSE, @@ -46,7 +115,13 @@ make_model_comparison_table <- function(model_serospecific, data_serospecific, scale = scale, summary_level = "id_antigen" ) |> - dplyr::select(.data$id, .data$antigen_iso, .data$MAE, .data$RMSE, .data$n_obs) |> + dplyr::select( + .data$id, + .data$antigen_iso, + .data$MAE, + .data$RMSE, + .data$n_obs + ) |> dplyr::rename( MAE_overall = .data$MAE, RMSE_overall = .data$RMSE, @@ -68,21 +143,26 @@ make_model_comparison_table <- function(model_serospecific, data_serospecific, TRUE ~ 100 * .data$delta_RMSE / .data$RMSE_overall ), best_MAE = dplyr::case_when( - is.na(.data$MAE_overall) | is.na(.data$MAE_serospecific) ~ NA_character_, + is.na(.data$MAE_overall) | is.na(.data$MAE_serospecific) ~ + NA_character_, abs(.data$delta_MAE) <= tie_tol ~ "tie", .data$delta_MAE > 0 ~ "serospecific", TRUE ~ "overall" ), best_RMSE = dplyr::case_when( - is.na(.data$RMSE_overall) | is.na(.data$RMSE_serospecific) ~ NA_character_, + is.na(.data$RMSE_overall) | is.na(.data$RMSE_serospecific) ~ + NA_character_, abs(.data$delta_RMSE) <= tie_tol ~ "tie", .data$delta_RMSE > 0 ~ "serospecific", TRUE ~ "overall" ), best_overall = dplyr::case_when( - .data$best_MAE == "serospecific" & .data$best_RMSE == "serospecific" ~ "serospecific", - .data$best_MAE == "overall" & .data$best_RMSE == "overall" ~ "overall", - .data$best_MAE == "tie" & .data$best_RMSE == "tie" ~ "tie", + .data$best_MAE == "serospecific" & + .data$best_RMSE == "serospecific" ~ "serospecific", + .data$best_MAE == "overall" & + .data$best_RMSE == "overall" ~ "overall", + .data$best_MAE == "tie" & + .data$best_RMSE == "tie" ~ "tie", TRUE ~ "mixed" ) ) |> diff --git a/R/predict_posterior_at_times.R b/R/predict_posterior_at_times.R index 5b59a3b..b7a505c 100644 --- a/R/predict_posterior_at_times.R +++ b/R/predict_posterior_at_times.R @@ -1,18 +1,24 @@ -#' Posterior predictions at specified times for given subjects and antigen/isotype +#' Posterior predictions at specified times for given subjects +#' and antigen/isotype #' -#' Generates draw-level posterior predictions of the antibody trajectory at user-specified -#' time points, for one or more subjects and a selected antigen/isotype. This is a low-level -#' helper used by residual-based posterior predictive diagnostics. +#' Generates draw-level posterior predictions of the antibody trajectory at +#' user-specified time points, for one or more subjects and a selected +#' antigen/isotype. This is a low-level helper used by +#' residual-based posterior predictive diagnostics. #' #' @param model A data frame of posterior draws in long format with columns: -#' \code{Subject}, \code{Iso_type}, \code{Chain}, \code{Iteration}, \code{Parameter}, \code{value}. -#' @param ids Character vector of subject IDs to include (matched against \code{Subject}). +#' \code{Subject}, \code{Iso_type}, \code{Chain}, \code{Iteration}, +#' \code{Parameter}, \code{value}. +#' @param ids Character vector of subject IDs to include +#' (matched against \code{Subject}). #' @param antigen_iso Character scalar specifying the antigen/isotype to include #' (matched against \code{Iso_type}). -#' @param times Numeric vector of time points (days) at which to evaluate predictions. +#' @param times Numeric vector of time points (days) at which to evaluate +#' predictions. #' -#' @return A tibble with one row per (posterior draw \eqn{\times} time \eqn{\times} subject), -#' including the evaluated prediction \code{res}. Output includes at least: +#' @return A tibble with one row per +#' (posterior draw \eqn{\times} time \eqn{\times} subject), including the +#' evaluated prediction \code{res}. Output includes at least: #' \describe{ #' \item{id}{Subject ID (character).} #' \item{t}{Time (days) at which prediction was evaluated.} @@ -26,8 +32,8 @@ #' @details #' This function pivots posterior draws to wide format (parameters as columns), #' expands them over \code{times}, and evaluates the antibody curve via -#' \code{serodynamics:::ab()} using parameters \code{y0}, \code{y1}, \code{t1}, -#' \code{alpha}, and \code{shape}. +#' an internal implementation of the antibody kinetics model using parameters +#' \code{y0}, \code{y1}, \code{t1}, \code{alpha}, and \code{shape}. #' #' @seealso \code{\link{compute_residual_metrics}} #' @@ -50,7 +56,10 @@ predict_posterior_at_times <- function(model, ids, antigen_iso, times) { param_wide <- sr_model_sub |> dplyr::select(.data$Chain, .data$Iteration, .data$Iso_type, .data$Parameter, .data$value, .data$Subject) |> - tidyr::pivot_wider(names_from = .data$Parameter, values_from = .data$value) |> + tidyr::pivot_wider( + names_from = .data$Parameter, + values_from = .data$value + ) |> dplyr::arrange(.data$Chain, .data$Iteration) |> dplyr::mutate( antigen_iso = factor(.data$Iso_type), @@ -65,14 +74,18 @@ predict_posterior_at_times <- function(model, ids, antigen_iso, times) { dt <- tibble::tibble(t = times) |> dplyr::mutate(idx = dplyr::row_number()) |> - tidyr::pivot_wider(names_from = "idx", values_from = "t", names_prefix = "time") |> + tidyr::pivot_wider( + names_from = "idx", + values_from = "t", + names_prefix = "time" + ) |> dplyr::slice(rep(seq_len(dplyr::n()), each = nrow(param_wide))) predictions <- cbind(param_wide, dt) |> tidyr::pivot_longer(cols = dplyr::starts_with("time"), values_to = "t") |> dplyr::select(-"name") |> dplyr::mutate( - res = serodynamics:::ab(.data$t, .data$y0, .data$y1, .data$t1, .data$alpha, .data$shape) + res = ab(.data$t, .data$y0, .data$y1, .data$t1, .data$alpha, .data$shape) ) predictions diff --git a/R/process_shigella_data.R b/R/process_shigella_data.R index a600277..d1f56db 100644 --- a/R/process_shigella_data.R +++ b/R/process_shigella_data.R @@ -5,8 +5,10 @@ #' `serodynamics::as_case_data()`. #' #' @param data A data frame containing longitudinal measurements. -#' @param study_filter Character scalar. Value of `study_name` to keep (e.g. "SOSAR"). -#' @param antigen Unquoted column name for the antigen measurement (e.g. n_ipab_MFI). +#' @param study_filter Character scalar. Value of `study_name` to keep +#' (e.g. "SOSAR"). +#' @param antigen Unquoted column name for the antigen measurement +#' (e.g. n_ipab_MFI). #' #' @return A tibble with standardized columns: #' \describe{ @@ -33,8 +35,13 @@ process_shigella_data <- function(data, study_filter, antigen) { data |> dplyr::filter(.data$study_name == study_filter) |> - dplyr::select(.data$isotype_name, .data$sid, .data$timepoint, .data$`Actual day`, - .data$cohort_name, !!antigen_col) |> + dplyr::select( + .data$isotype_name, + .data$sid, + .data$timepoint, + .data$`Actual day`, + !!antigen_col + ) |> dplyr::mutate( index_id = .data$sid, antigen_iso = .data$isotype_name, diff --git a/R/shigella-package.R b/R/shigella-package.R index d809e4e..a0576ad 100644 --- a/R/shigella-package.R +++ b/R/shigella-package.R @@ -14,5 +14,6 @@ #' @importFrom ggplot2 ggplot #' @importFrom ggplot2 labs #' @importFrom ggplot2 theme_minimal +#' @importFrom rlang .data ## usethis namespace: end NULL diff --git a/R/utils_internal.R b/R/utils_internal.R new file mode 100644 index 0000000..91ccd9f --- /dev/null +++ b/R/utils_internal.R @@ -0,0 +1,51 @@ +#' Internal utility functions +#' +#' @keywords internal +#' @name utils_internal +NULL + +# NOTE: ab() and get_timeindays_var() are internal re-implementations of +# functions from the serodynamics package (ucdavis/serodynamics). +# They are duplicated here to avoid using ::: (non-exported access), +# which is not allowed in R packages passing R CMD check. +# If serodynamics updates these functions, this file must be updated to match. + +#' Get time variable name from a case_data object +#' +#' @param dataset A serodynamics case_data object +#' @return Character scalar with the time variable name +#' @keywords internal +#' @noRd +get_timeindays_var <- function(dataset) { + var <- attr(dataset, "timeindays") + if (is.null(var)) var <- "timeindays" + return(var) +} + +#' Antibody kinetics model function +#' +#' Evaluates the antibody trajectory model at specified time points. +#' +#' @param t Numeric vector of time points (days) +#' @param y0 Baseline antibody level +#' @param y1 Peak antibody level +#' @param t1 Time to peak (days) +#' @param alpha Decay rate parameter +#' @param shape Decay shape parameter (rho) +#' @return Numeric vector of predicted antibody levels +#' @keywords internal +#' @noRd +bt <- function(y0, y1, t1) log(y1 / y0) / t1 + +# Mirrors serodynamics:::ab() ā€” see ucdavis/serodynamics/R/ab.R +ab <- function(t, y0, y1, t1, alpha, shape) { + # NOTE: this preserves serodynamics behavior; when shape == 1 the decay branch + # uses 1 / (1 - shape), matching upstream implementation. + beta <- bt(y0, y1, t1) + yt <- ifelse( + t <= t1, + y0 * exp(beta * t), + (y1^(1 - shape) - (1 - shape) * alpha * (t - t1))^(1 / (1 - shape)) + ) + return(yt) +} diff --git a/README.Rmd b/README.Rmd index d17fbb6..73daa41 100644 --- a/README.Rmd +++ b/README.Rmd @@ -12,11 +12,66 @@ knitr::opts_chunk$set( ``` # shigella -This repo is for the Shigella analyses + +Analysis tools for longitudinal antibody kinetics data from Shigella infection studies. [![Codecov test coverage](https://codecov.io/gh/UCD-SERG/shigella/graph/badge.svg)](https://app.codecov.io/gh/UCD-SERG/shigella) [![R-CMD-check](https://github.com/UCD-SERG/shigella/actions/workflows/R-CMD-check.yaml/badge.svg)](https://github.com/UCD-SERG/shigella/actions/workflows/R-CMD-check.yaml) -The goal of shigella is to ... +## Overview + +The `shigella` package provides tools for analyzing longitudinal antibody kinetics following Shigella infection. Key features include: + +- Data processing and reshaping for `serodynamics` workflows +- Residual-based model evaluation and comparison +- Visualization of population-level antibody trajectories +- Mock data for testing and development + +## Installation + +You can install the development version from GitHub: + +```r +# install.packages("remotes") +remotes::install_github("UCD-SERG/shigella") +``` + +## Usage + +```r +library(shigella) + +# Process raw data +dat_long <- process_shigella_data( + data = raw_data, + study_filter = "SOSAR", + antigen = n_ipab_MFI +) + +# Compute residual metrics +metrics <- compute_residual_metrics( + model = posterior_draws, + dataset = case_data, + ids = unique(case_data$id), + antigen_iso = "IgG", + scale = "original", + summary_level = "id_antigen" +) + +# Compare models +comparison <- model_comparison_table( + metrics_overall = metrics_model1, + metrics_pointwise = metrics_model2 +) +``` + +## Documentation + +- [Package website](https://ucd-serg.github.io/shigella/) +- [Getting started vignette](https://ucd-serg.github.io/shigella/articles/getting-started.html) + +## Note on Data + +This package currently uses mock datasets for testing and examples. Real Shigella antibody data will be added in a future release. diff --git a/data-raw/mock_data.R b/data-raw/mock_data.R new file mode 100644 index 0000000..ce32604 --- /dev/null +++ b/data-raw/mock_data.R @@ -0,0 +1,85 @@ +## Code to prepare mock datasets for testing and examples +## This script generates mock data that mimics the structure expected by shigella functions + +set.seed(2024) + +# Mock posterior draws dataset ------------------------------------------------- +# Structure: Subject, Iso_type, Chain, Iteration, Parameter, value +# Parameters: y0, y1, t1, alpha, shape (rho) + +n_chains <- 3 +n_iter <- 100 +subjects <- c("newperson", "SOSAR-22008", "SOSAR-22015") +isotypes <- c("IgG", "IgA") +parameters <- c("y0", "y1", "t1", "alpha", "shape") + +mock_posterior_draws <- expand.grid( + Subject = subjects, + Iso_type = isotypes, + Chain = 1:n_chains, + Iteration = 1:n_iter, + Parameter = parameters, + stringsAsFactors = FALSE +) + +# Generate plausible parameter values +mock_posterior_draws$value <- NA +for (i in seq_len(nrow(mock_posterior_draws))) { + param <- mock_posterior_draws$Parameter[i] + mock_posterior_draws$value[i] <- switch(param, + y0 = runif(1, 100, 500), # baseline antibody level + y1 = runif(1, 1000, 5000), # peak antibody level + t1 = runif(1, 5, 20), # time to peak (days) + alpha = runif(1, 0.001, 0.05), # decay rate + shape = runif(1, 0.5, 2.0) # decay shape (rho) + ) +} + +# Mock case dataset ------------------------------------------------------------- +# Structure compatible with serodynamics::as_case_data() +# Required: id, antigen_iso, time (in days), value columns + +n_timepoints <- 6 +timepoints <- c(0, 7, 14, 30, 60, 90) +ids_case <- c("SOSAR-22008", "SOSAR-22015", "SOSAR-22020", "SOSAR-22025") + +mock_case_data_raw <- expand.grid( + id = ids_case, + antigen_iso = c("IgG", "IgA"), + timepoint = timepoints, + stringsAsFactors = FALSE +) + +# Generate mock antibody values using a simple kinetics model +mock_case_data_raw$value <- apply(mock_case_data_raw, 1, function(row) { + t <- as.numeric(row["timepoint"]) + # Simple rise-decay pattern with noise + y0 <- runif(1, 100, 300) + y1 <- runif(1, 1000, 3000) + t1 <- 10 + alpha <- 0.02 + shape <- 1.2 + + if (t <= t1) { + pred <- y0 + (y1 - y0) * (t / t1) + } else { + pred <- y0 + (y1 - y0) * exp(-alpha * ((t - t1)^shape)) + } + + # Add measurement noise (20% CV) + pred * rnorm(1, 1, 0.2) +}) + +# Convert to serodynamics case_data format +mock_case_data <- mock_case_data_raw +attr(mock_case_data, "timeindays") <- "timepoint" +attr(mock_case_data, "value_var") <- "value" +class(mock_case_data) <- c("case_data", "data.frame") + +# Save datasets ----------------------------------------------------------------- +usethis::use_data(mock_posterior_draws, overwrite = TRUE) +usethis::use_data(mock_case_data, overwrite = TRUE) + +message("Mock datasets created successfully!") +message("- mock_posterior_draws: ", nrow(mock_posterior_draws), " rows") +message("- mock_case_data: ", nrow(mock_case_data), " rows") diff --git a/data/.gitignore b/data/.gitignore index 53c221b..069e8b7 100644 --- a/data/.gitignore +++ b/data/.gitignore @@ -1,2 +1,4 @@ *.rda *.rds +!mock_posterior_draws.rda +!mock_case_data.rda diff --git a/data/mock_case_data.rda b/data/mock_case_data.rda new file mode 100644 index 0000000000000000000000000000000000000000..3c404e5f1bab0986c12034adfd36ad765b08622a GIT binary patch literal 984 zcmV;}11J1KT4*^jL0KkKS(>;sx&Q==|NsC0{rB^~|L6Co|Mle2(?HPCqW~ea zLMW!lfB*n!41?5Y20#D+28MtD00000XaE7A0002c01W^H05TdG0Avj`XaS&T$)F7w zhJl~}!~kf}0079)Xvi3iG-xpcKn(*R009(`njoGco`imq2c!&-Q%w!1JwVz(4KxOZ zLFx}Fr>JNKL((z;0BNHj05l$^jT#yN4K!^-)a=&K90(jx8XzaXs%7AaYOTKSJ9Rh# zblwzDOn}##VpL8QXOJLbanf8#vZDe-t5lxm)nYC5-F0ZnyP#6N6U+!qz-gfo`SPn& zttuwXW+p%a%mNG*?7YCjuptRW&S|&BnqyEzMTopooFHTjfx~OdMTH{0G|g7M#gQTa zr9cQ!U9^@&FaXULP)!03x)cz`fZM^=c)FZx>N)(0C;=z~OE*{adpb&no&;OIf>&^y z1Rrbz_x{wGp3CXi?Q|DqPp;LRlK=+$Jl}ax(xS}5Dn3*(=3JB-k}048P5>YX=;U@t zyBMID0O9r>jL43YkVd9nNaXQ6e8|!W;{q~A&!;wOda6;9H8wRrn=_-%&zP|Xc>$v| zVR{CmWB7L%D!&ACi^oKfh&e@R1jSQU44yQ8=0rkoLt}t@4ff2oQ)alK#Ixv3B@<#b zAn2O)TLnRboCI zmcQ70xd)kj=95i;kdA)E*5k<#Y$v_C&TgvhaUV!G)V+lxH+m(c+vM+;yNP$L+(MAW z$0POGxlGOVs;X8_pL$901mE+v=1stoHjr$f*2TrF{j{~Tb+*`cBxS9Oq>69DJB=Zo zgl0(q(*X)46EpDCA!2Iah*)Cg%ygmaBy5rpkM2+aF=#k8@1G?fP-UHPe(xs(8_}F) z1DI2ykb(p3U=B|?qCo^G8W!u(2!?~l1_1RIz|z~XvEI&w?0G%%TERO8iZ%Pl@`_AC ziXbdDSK{nAM#Th(C>~W6CfS^RLiIA0Ql*fpCIE6L&LA8q7I4h1Qfbxy0l4{4p-yRl zKVuA`Oh#t`zyK903P^PcdcAS@C!hu6BYjeUZ0WUTggKoA3UDDHabSY~i@744C`e6Q G8eIS;QM+jX literal 0 HcmV?d00001 diff --git a/data/mock_posterior_draws.rda b/data/mock_posterior_draws.rda new file mode 100644 index 0000000000000000000000000000000000000000..af6922e195ef69b3aadcd0281564e44b09586f65 GIT binary patch literal 56194 zcma%ibyOWs%;*I!?(T5W;>F$F-QC^Y9g4fVODXQ|?oiy_-3pZQ_`Uc3|Gu0(XJFTsICWjwm|6##u{a)N^PTVi-iM;H(XX0KOJAH&ac-l9Kj-5BWbuX&%(>>aY z&XQxX4B^amY~>nWx?Td8)I40mu6Hp2h}G9GP0trg*~97dYAsjONK@5W399k<+t0j! z%!!rGmU8uz?pzI@neXN=?H=Vk-HqFwJRXG0^W8jM7qwnLeOAJ0t`0kDIout1T75b@ z8&_IAx;81UJHzj4Vqw`Mv^B0Bmu$Xw`gGGV+p`ODJB3gNfDS5L?H|Up$L)#_Q^$2M14e^P+>q3V}Hu1 z$WZ_khmvEU2%;gWuSZ~_Qj(x>Fq9x*Ab=ShhJ=)S2t@<{Mvf*-DlY_N!l5J;6aq@x zpa29JNL7V^dWmQQFeLlm8}hNF_ZWgn(3D|FbE#sIyCQ>ts6-JN!LVpWFeLq`%%EUg zR2Co_Wni%AZxqRFa9DIlAyP*e@-L#1bD%^XAlkPt)FX;y{uk<2h%^v}Tu}%Zg#t0TyR*7IU&-VFM>h9C{Iu# zHxP{m5G*E&B3TIzi^(WNS^-0@_C>G`l&Aqjqy0je2E%anM(s0*aYAB)s=VNaJ|(Z( zJhxJMp8s~oB3sDl7-_C`!!Acuh*!dw>mrSWl;gZRBO_pT@I8|Hv$cogDruOE7e8wV z`MlvQ>2irvM?fn9KtRZ9bs8bCOs--2u(HOdEsn__dR_9f8gUiUk^|tdD8aBeb3$YU zFyz88XsuN=yQxF&hR;48f558U~C0 zMUk8Uhs9MCBBO&LH|qz+^nnsHfoLp(!Q!?kk~3eZFVeg)_yiQm*V9&qwk$~!1V9RhZVUX95(34q)rTBQ000x9kBbB+IZx1}K!H-TFpGJP zfCMRsq8N31mC}9z-0mj27<~fz>(#nmm$#{p&F3QPH;`|DTvSfaLI8|bW3rWRx+-XB zM?*<$IZWm=qCO3rP^4_i5wfg3wKn*>&3=1lBZID)f?#?&L17h#@}jl99R5#2L;X72 zVmtbj2;aqq7$!r?WbIWVKhR*36qh5Z7MKIDQf3g+(|F^D5;JX33ZH9+#c3 z9{f;&8&$f8gjN; z4b{Vbo7M~jS%5KAZEjbU4sGd1tx$VY~V>;IA-GmuVv6XdIqu`{~eFbl>fyg0rr__LTNTUC6 zN3$ZF&oknZ5g}wz44cwdD?e?DCM>S(8)^H@m|TmWs>>_JZj1R+&rYcxr(v>;m5x5A zeHCd1oyxue+o6!l>E;4kmn72YgB-LkLo&uk9Ft0V&x#q=qL7UfN0Cqaitw0#Tsh&w8NyG{MFsvzlj5Wld^xHd;gmuZZ zEGdI2sEt}U{J2G?EYmK^jqCQXjKCyMa1;T$tJ

{X18X-r+H?1m`c8; z=g&=Q({o3H{}CA+0i-@NYxNcn*h=_Fy z_5>At^}}5qUa3?FDf9N8s@QcF4c$~tzc(Z}dXx3NiUMtjyA7cpmNm4DbE8O<%_fyP zJhOKvy?WSm9g76ibylmK4YfR$`dGeKViNJ>a0RX==Mg(s1#9}$UYR9%m1dO=HY*Vt z=Jm|1%ZNlskvSWnSFWDoY@HGy>B|8+F*x=M$rAPR4oKLjq9i@|MlCCrNbba8gDC+E z|5I;$xHfUGHe{XHjco&OylPjtK>XV1SaZ5 zKmJ;8Ej7Jus??5$?!9=?|F1R=_If0of5t0$mJF1X3iYg|)ufDrH3c+lhm@SeG9Gz_iqHeTLzM-P zX#>&h&jTDrxOe;W`hyDb=tEKB~>WHslQ4?#)Y4Iqu3o(=F(f4AXJk&1?gdq1@ z=Re>Z?l!T8++1XFftkZ*_2(eUa-=1&&H@H!l9K@J`xDiclGUc#iNz~eBYsAFPiF1y z`aSUN`ZvP%Xezf)h!b_&Aulu2zbEA1%FOX!?452AB-V^BF0PS1xg4k(dn#4wdnZai zVHe6vLaew?Pb11UxAc8n750<5nW*uUSxc86&ZK)Y88@k$w_?ew{e;%^u?}Y$p!-f0 zo+yV$5wOpU#xV2t)z?~aHn7Gv7 z=@6)t0rNKNA)KM0Yg@*>#3X(Lj#Zp8&{YgN`gAm9+OLwRln3!lQm@R8qQf7^o?glH2$dR ze8K`PB!yKuCCRV8?aYc@m!_tcaC~>iAxe3U!DxD{^s2b>XQhB{2$2?C6GNgAO(PiE z7$mw0mb*sXx-IciT_V)SjnE%`qS#~YWqe)lPo|KnVbEn~>MzdzNrJhd^OygmM1Z+zA2|bDeU3MEf-yeQtep_cxB=*f} z`5WrO))5>XY6FhBeLQWaSF}-vq=Rmc$@fI>GJTtXYb7LRUkJw}sX#aM`=1DHyLC(u zvjWY+p3!h*w&l67ovDjNo9aL%ZNIDA$F{QA8n*jNxm-$!zICyA^84>ccz_EUT#1y9 zEaq~W;LMQxF=9q9(DlTg*qt0eR?c(Ovhf>7D&yf7yrQl>GjwEFHa>3`ys4+_UDMp3 zTX}`1vZB`Yk*#4_k`D;McKheF)hyYUs4{W!wI^(Cj&y8JcW$0eg%4II6U<26lav&B zkzZZgpMhk_KsKQ&oKKyaz17zyUc37#l`R+%k-s7@w`z+oX-;Yk)(B~Vo zC7;BoPvNZP=OII%iHI;dr)Z10^NxT43RQwDT(-5Nw6^qwT=)k~d@|6j)P1A%ZjS9* z%blvx;*XmI(W)`t-180}xIg%ne#<@hBZC_;`-);InZL8sPmY+QIQ-~Apc&`wFNAAg zAWmI6kb%dS^>G|ZkCovxuN3xo&ePoZ#7(>l6%X8tGfs_-(X&;cHF}q!;$6v51BZNS zv7aBGnKM#eJ<0lerF8>>&&*5o-)m#5$Bl&ym_h9wVa9;V{OB69$5ND z+Zcz`h)KU9+3O>`+@^IcnsqY^^`UN(vpM6`uS}VW)dpc;)j3ygE7j# zJ?il+X-ySIf7uO7-+FpC8yt``Nq23nHnry2<*v(CJ-6^+i+rDRz=kgYGXvxBM3g?hcDa%0f%9hd8$>i{3O=N?Is2N9E2{uz!{@D*VL-hgK z<8HBxs~aRryfmBpu$kc9RdKCdGJ2OBjP*FA;ZKsy=Lm8aE3`N22Rv?CJ1=-Gui(B> z)dw^tiz8k4==;)%f6H8&LrBN|(a(qo*bm)!NVeW`OllR;(|Vo+w3Y7q+IzDTb}}aw ztz?XTpR7jj-fg4A!^CnfMs>vurOaHF&u>H3&cpRJ?i?u9&RSZtxKc3$|;#-PJF#I_z|Nn^AmUr61~a7|8b zF*M!AxTUh=ArT)hFYY4Nm+JR=#${*QB=>~zc5esi+UTn7{M7Kn*~twEEyy*Xv8tq#<2{kl)4?aamETl!n`=ivC$=W6{%FMdp#j^9 zNQjbD<{P?s&s)ck_)5SLL|?miM&o zN&1H+M*f!7CZ(!yNzPVr!&BAx)%sq3{aR5W>4((bC_e>%7*xRq`$gM+SH}yPlZS~~ z>xEX({$*_qxg*cpaFSt{&5l7Qg0EGlUcDQ`!QJ7Ea2B!F2R{ajJd=Nh?-M2y{fkQH zvtzX_)y1}XXPj4WCjfmZHn5_+*8@;s6;MfaWtx9SMr=dT2e3P(r$RSpeWUUhxA zQ~m}SDn6Y57`b2XL#i5N7Qf8Y0N4uzjj^ON|8(k#jzwFsWF?vh^M+qFq4~281lho%nKVK!=kc}NBG_g zEt@rf*>aIhzfnP?ufQoDg6>M~A7VC81%qgP;2|zsdSH%gW1CKFv=*(fL`?bjhaPik zi&i&;=`ZSDX?)!&FeU%U;Y)H;qnvv>WKo7vCj57|1}Wln)U!ovCz{xzU+DK_Wyjs+ z)D#uX$R`H5w5qK%;ZoDU9Os4u6^!b%Hx&;P^jN}4Tk%$05z}JIq{)%FmV-{-??ZS+ zZAq-HX>Gq=m+gDMD>-0Q(=1|ov3CR96xf=AYB8@YZalB9VNGGS2h( zAQ1O4mLY>x%$Puy{zV;!)zBx&ZLqi24DI%ik1GMXO^;Gm;J)}^bati{ z4#@P=&(3JjN%r2jmQrjo+;w-%T`osCG>E%UH+t?MKt8tN%B@avaxITGAElRfOO$w< zHox>aAq-}nYATujXJhKDwz+o>Cq4R)jkbW!W<7?<7BlC0V2P1Zd2H3>^56k8jG(Xf z#S_b8Y0jQ#dZ>O%;B?LS!$87`_(f#XjwVgYoftMvi!ogn#{t8>$qw2puUG%Cv_H&3 z71r>JbH>Sc483Ux>~Q|HLivB%ArpBbG&J#;h=TRmJ1#m^**0wOULQve5A-bYdIf9+ z-_g0CC|)QiFp_@uvZi;gKb@;zatWj{s3TmFccI7hz}Iv`)!4Si?1Us=WbZM)@F?|le$#iQAc*1Dn*}}A>R?%q{HiVM;=J513)WlcfE468j}dKgF9? zlW%f_GnPh7;?ItxYd5l~G7UE$kLE%B&jnN6k9@KEu{i&lQjCW;iPUIYK|5~FMqttP zZ{}PoQ6cdUwRKix6Ta|ET3MCp1oz?iT(cJo3VQ-7ThIqg_|v75DpQMpo=6{l7LZ#@ zd9tP&e>4}hHzV?E60B|4QeMtbgW)c99>VR7Z3ra#E zM}s8+d^J!(2~x<56H`3kUjv+khQ6Yq&1u{S!@Ba1-|Cy7Qlx;vk_lt*fnV`;w7RcA zxG2%^XccD9P0HxQy3 z2>?KzB>a^JAeUAG0)zllbL0uYRB!-!ek3>m_#1^hSV9jVM4Au;j0Xpz2r}bWm-d0f zV!6Tq(C-R_NZ*2hC@}aX06?xZFhLMhs%EYeUev}c#c;BnvwL`k&ou|ye-7&3%VY`R z0-Ufwe-qB+3ihG2#NnD6pG7z}7rGo1>`xHdyX9o^wnPacNApEuEQjvq#{nWja0xyp zMC_iF4|%y{W{MgfCb?umhkMJr%vsb)!X$8VMIoRGW#-M3vXW4dI9e*3Td*WZeFM$7 zp`SLxl^#PGKMtTtcv9 zNi2Z6=*hsGyct+7-~e#eOBfr=(S?(nECU8b_63IMq=Xn1x_ckH{U$7P>Q4mcw+5ko z=)OD@6w>}=s)WWl!@OM7y%|rJOtiUj!6LyI7GV;3FD3L)iMD}EEZIWYBZ%ZA1uw&p z4Fm7ADI%;hRCoWb{VJ zA`=Wr-W6i`hRw+Z&6Nv>(gDpFbeZ6@opY)SLH_J)F0bIj=n@jbl|z(y6VZQ14DPcp ze3>g+S%EIUCObR&=CXo;d=E-wj6kfby>gWtY=mf|AEJ)7|fMJ0MEGtGkLY~lx)8hEbx>!?9YTD zm_Zim-zVc54}~Rx*nAw%;2UFhMo2=KC%~c71~0JSIvohv#LdDMK8~0P+Q;U~yZw6I z3@ZBO@E*ojvVx0~6*lkLC4>bk?8eb&^Kh2#ar_bt7;XKP)Pg~gJU-C`2$9JZ8N!0V z^aa6TiIV^*k{g#9D3Ynunkay8{(+#xj!Q;R;;*knPB7b7PfRdDQk()U_4TFqPu8&) z?WWA79iLg*&;5leK0hCvsQ4h@;hAq2&+M-TWwmKv3zRTK_+Fo6(Z-(9xj$6c@HIb$ zNIH+tS-I0>=GO{VD_FW?1>=~>!2y7C%3EKpzx3Kwo1>DlJk)^TL&v>rGG%+GtyRQ; zuHQbdK2ILqm@+@+v+G-JpH{yzZ&}Ue$y*^m<$Lx}j(pH>4b>MqoW_h=xp^6^reA9v zn2)qHAxW1?bnVZpe{w`;=(IV$uUjg;)IPCu(slD*d_R!p^_ZbVs(2ye2Tbn4Ji*5) z#|~*)$l<8dZph({Ir44ns>`gcOlgBkds3GY{Q$uLK(Xf>|7WIq*+W9f3-w=rV6lGF z>xG)KJam8k&m-7hKGa6Hn}BFsKsZ>PQ6LHe<5I2fw6kIW_Jvj}#cigk+qryudxc(6 z@nLC3ydE|=;I{Ii%6wxSRJg{bI~mYvF#6isSny>OTr2@!y7TAoL$3eH%)wvAI7U9~ zmxWPGX;*GdyU-X<9&=hg6E&L%wGhwXCq>S`Z8pSQR!IC&vIM!?DV=youc2HtJTA#VwrK zv+aS#jX@<~Q_k1vy6Ck|Y!AZx9^#e#po8lf!Wmhz+XWWL*!P@Mc}KdWzuMBt1HAvM zMnoU{X0MjJVg&62mpa+kCxipQ`I-#zBKHbz%jl2fL$H|@SAlqI3%&b?W|eo*fJDMY zW>^pn9`dp2KcBY6nuhrq(?2U6saY%q1~@PPKgA+t%oy6Yj~F{&U?$wb@j5!G4VFzGFT zZ`jaW9ld{v#isDcho;W{EoXgy_iU7EchTJc5(97TiKfKo_0ntUTYw+{fTBFLulx_S zb-TUn%xL(Z9J;S{kW~2cgf7NMSC&&}()R1C#ivocOO@X84OL%v=*zE!z9-8{Y@0=S zXX2)5j5do5&pzL=>9tI_d`M7wPLhVZGk&B&&(uH1|5e-#hSMMiO7p%Xme@6GA zaWq!gb6hMhr(xFe2oJqkd&9$|IXk#&Hb<*Wr{7am)$IFEx5|ue?70^Q=>9dJV?AU3 z?V70o0)oL1>hd9r@qBJ~t&hYIn*}Gb3vZ(}-X9&Ar%$ywv#6f7@d7RnTfYmjM{kOy zbbcE(4B&Q0N<2b-jjYALZDmWDguhLhY@^*7q|k}oD=AAd2rsL2dGtSzJF{}eZ+gKX zWfV@wnQhY_kA%ZkCU^V^4+Db6mA?tUWPv(toj_l9sraM!ZRb=);b6RrjqRtwIahJI zhOfsHL1`Oct0u4Ec05JeE? zzpaZ!>`c7lb$=7lA~EP@X=IfN8He{*h&lBDJ^dotH?ZhF!#C*$pD~Rh{<+Ieb!j&i9O5b>e4~n5Qg&(M#H(cn2X;(%Z#KV;@DyjEY zdN)1Ei_rw5uiNe{4>31gH$GD<4coW$7kZTKxILlixa%S1Mt7bd90&yP+Md*%M1#R& zC^Uqh&<0Y9H_p`Z3QTvAU1oJLOgA2L`#k$G-3Z-fOfz6|rvL!L*XA+`YSp+ljcDQv zbEUFsW~uIKLjCNvXkA|06K$7mRe8uPEIgN1B<&Piv3ViQInBMyz*Im{7ae;26DpDV zPn`MfJdre~(Yz6cp+isA5~EtB^dZUeN`cjv+yZT#4cg z20Z>ChqCyP;+LhEXoW+UD@vlhq8`}syem)OyVIN)LG!(ZP4xP-Xr3lf>=_EljE~xJ zrG7yW%Xf}Mu0g9YIiHdtAxmkBT8!7vF))U-_Kr#WchM0aG*{JC0?H=6f5QQw2!Ye{ zY{eJz*_q#LIx_c~T~RUAsaY8f^5o7>m}N;91RCX(;$bEVf(3&8bk4Bk*@_6@o8%~% zNj_R<=Y9M2@cJW_QhK&LIZU#l-G1NOOaDoL??Ea~o=@up4q8utV>MrQV_3J2EG1{L zxG$V>s~F*dBJgd4)QI3eVSja=kNK!Az^w~-97+A_XG=={nK>mD*0RZ-lf0JLTNHq~ znGmSzYTEc`5`)^|03H7a9Vie662{kaOl^G0w+3U=&F>Iy%sp3xC18r@+fFT^ME7dt zx=;E3KWH`THg*N?w?Bz^$v&Ah_*VmvWxW#V11&Bl@Ns$r*C(}|s8Ev+Tllfw>8{{k z2d1MW2NbPL(H}!aH^V}g&5RAZ#ces{K?0vq3AlF}KM$ytZc5ihQHY|bc;zLK{wn=p zG09l?b9s77)I#$_ET3yei5N56j)(KAw0WbK!ZlDI6$CT%ld?5aBxh>v0`u3;D&ZS2 zdL!Q8Y_C}*F_xFE;S1})fj0S4VU8N(#X~nB%GWWjAHb`Ld7IZ|)dmOjn|PvgQ$-KJ zax?3_K^{giNiuTGC(?d-RAdwLRJ^P}<|va;pyJTiS{j;6zvkFewBf{SS}plTyWhBV zH$kC$daYXhxJVqcWb-QJb<5Hqd10Udbx6=jv;!^L<%ihq_p#`NmHb+&Rz$t&rvPJG zouqc3ep}0_W8>t97(iNtlgEFD!c(Wa55pz%L!#6bZG&h_L!oGT1i(h1PIZ^>=r3=! zIUnr`&Up|@f+B`_{FD;0+?4)Y*I&UGQOS0UC^%wNf_(2_=2!erpN8KBv+pQ!ZM9Pe zN7Pxxnnka4y*BFoBm4)&pCjBg|G=Z1iZ&ZT#15Byu_^r@06`d#4$tamEmAm_4W8TI zBF0U53XyGvks|C?Fgzh{hNE|!jr0=}lkYHrsGnfNZ$3A(+gc}kPTxLuR_JHIj8Yy) zsnxo1*BPo*%OBA&s5>*GSGby&pfnsqkk;0Va56YG6?`^zQC67y;CGwXgapt zZH#pky?fwIhcv)b537Uvn;Cz67uj zPqJ`q5GtGWA4HvWD;-j?<(Jh1OuO522bx|{9J`QNPV(n)Wxwi`jkAl0m9=ZC8Kw^b zB61E=Mz~dTvLce}u%581u;pt@Y4HiotX6Dwp{^hKn%}3Tp0<81NtW>-zxqQzRs^)^ zS(rDSGeeFMveKs8!^7*bjYPV?irp`W?Eo7$FD=~sPF7yMFoNS)N% zv$_w=e-CH9nMR8nLRnp2(KksK;>g?ngj87> zW@du5i~VIi>eoC&uES1YB_6TddZ zQ(cNxo-uTM{UC2P&M9`4D|}aafIGqk0CJ8!z8faN@acN~ZAGlHSuu4oeJ8|d`jrYL zykPtAP#@!TP9g!TnLri$Lxxakhx=WNQH4hpg|FfGpsq{4`Z3}O;%XJ!bhz^jGkZ*D zlQa}_I=YLCX;*%-Uo%lJCyMxQIXXSQZjyGyvf{s>s|MYlTAnDh5HnwY2cq!c-{mMR za;)>V{;KJA;wc6tZ$D4*xT-u+u986Ry91pOfE#ne)7-~MC#OsY&j;`i5>W_*48^E< zG2b{8C9@6Cw$(^0HKEtXc#tKNVFq1$+I-s8^Y*+98} zbZvC!ejldO7QL*z#AF)Fildlj4qlVW!rf%|mU+!DdRRE^rns;~}$u zT=mIBI>7Q(HROIet`<7qc9t1YrL_a-Yq-qf0WjRn$IyxLsm|=fNbqkc{Cn4w{ufQA-xx>AEf`5->1N)@E~B37-Ihj3S>fLNLOp4oFzf96o>kA~UJc znXtfm-tl&KN>%M(C?gQ;u8H&+qqz-W*bfY{mqB7u!y5N6+lrb;li7CCTrC;GfB}(n zcxPESOoHeuYHJLL{2^3B?hd@09nXe_31!`3plz8+qG34wn6Sqpnt!eO;{X7ZMwXW* z*(iN!7}zn1ls6&$kZrVd<*LSh*ELQ z6#VDLlf~C$y$zU^^y@bXDEfZI28zqA)IyDeW*scdObggXmGyaOg>{*zehZOEQx*}x1 ziKhNOC0}ijTd>i20tOV3Ve~LN{u?^6r|_N3wzDqoF4ZH0~B8nbXo@HrNKPiOvb!v0h{oH9x1)_6RgB zOfl$RC}!JLMm_i}ZJMy_s-Ne)?N z1o`kRM3g`PwnCk^f1vY*MReE;vz#qg=l%LGPUTYln{WBkU2~KxkmnCpqF0}?_rfA5 z>IL-@Y@+%)&j$DceCg!k%)K^0nSL3GDU}p?ezctD8Q)JFAkZ@TSONM1U5T!n^Vl^1 z1eJnA6a;6t{HIO~e$47UPGQj4@CR>F-pQyHyvfKgn@>w?K#*ID4JHv_GDmiLH*R z`og-l4HZJpFP^z4r$2F|iJ!8PTyVi{5NFQ53?!B6pvqNZnY>Jhla4;jdv+FVyhzex zra&^*lu`sZuGpr=UUr&}w|>xR-`xbh`GmE?9yI3CL*Rh3*BMG+;STNRxj`DMqjD!n zY|amBP98Li8rE+2cXSKn{pU!LXvhfC9=%stoub_cP@P9ZK;a~N(f#`Tc4mRPi-)7$ z&4fT@;WwrV=zHzBTa7G?9VVpMnBYVx2<7l9Jamq{p+QEG9LGLdPaaKZ&{e3n=} z$U8Ls$28skwC?tkS3-DsqvrfW)cTeArOoEW-q9ZJg-&v$vA6hEK4g6coJT4rKI7c_-6P@=2iG^LTyBU$oF|8T6Q*iU399({Q zdvUM3_qu+b7BcjrS=eSA3*rQyU7-JMt4tF+WW<`aNEmkMt8}GI8i1PMxi!d@u5InR z<(ndRBmoHv$|taKo^*sqAOvAj8VU=BhCTss!9i8}_`ChmH2RO|*4vd!-w*!84!T<1 z@#u0Zs6Dgkh_L?LreP^>WaqHInjaL87jj)eMG%OpzvcCnZhk$wu0<|U^}(QNvey)K zDp$&qP&SE@n{;sg*-I`6uo8wu-q|&t9V^HDy}s(VZeF~AAq*6ZND@w)Um}kGeSJ+4 z=M(VHdPir}?WcAdO-3w4&!^Xq=H)&a1D9+o0y5TuVDSZm^n=YC)RwCFWSOdjT?(R` zMi%v^p`V3nX&-z%F~rOeI><&#qRba{jgoFk87g4R1mCmaT}Orw-|7h^vS=2U2jQ$w zUHnFE{h#qSfeaT5_G{GDt?fzDzUUv)9Z!3&?&&VH@3#rCzM3%4VPXm+;Dn*<}fmG~%jhnw585#l6>dvk|OBurZh3S=8MAw!KR zTR32Z8miGVwW0(6i#h&EH$LcJGQ6Z{(JXSTblMIG-(4j<915)qS|ghN`vW;4dnZ=! zu$WC-Vex`jo2l*9XPR4kG^PdyS-Pl}Ei4|)G-0`TMk24WGp&5rikqehTCMVz#2I(O z3WHoWEM3%G^?tXjk9;**c!?W}0)cA4I6|@0gEQxk<+vB$C~;N<_RJqI=gtyn4CN+J*V%bN33)*e>jXpT0?8D`7KdM+|^+r=6 zcXo`b+UiP4?xp+3U+Hgscb;^HeM*BNto$WAQuun4&)u)=?4M)e&0&s&wRf(1WB%ce zw$fzr%r*Os%>Rqwgm3Ff717Ny(6|jqv(R;ykY2W31gp`_>V!KE>bZf^07GLz-;7?B zs6@+Msk&iz`c8|lY#)DO>Quf-zVI!Rrk=!AdAP2I!{L+TO4~{rpz&1_k~bjtJ^pmsfeZ|OUHC-QN3?Rc+k7}2O$zJ!clF5mV> zkN6ujX9-rWf}TqvQp1>Ub)kVnRlnrlG;4HaAlmojTJ^?Z{ojUa=zAl7*wpCAV^n#j zE#eJtPOzC!3cH*wE&YWDn|m*ImVPRsRrx~~!2rN}JBov4xN!@He(1{>7-{6B6;+jW z6e*D`uj8Jw_Y1)(dI(}gNna-Z9#`7zy8O-|%#@KKwj12Ds~fn!KA2T5PLSrjAclro z1|3OY|CuD~9mQYDvj^fH;Kx2x zuM7MNVIfvAFiol)5}3q2hb_K}V9up?B~>iBzhkamfvQw%d~Axs1o{q5kM2?&fMF>0 zBnAP*d=n#&-_{bl-dUX&C4ve-@4y_;OPdyUXMZ$!@JK#+I;S!hTh#Wt9`f2XBVGxE z0lf4_;_po$yh5u-RNT%iKeiM`YEm9&YKg(?o%*(Bwa{;a7a)XGE)H{Tmmim9L5Vf0MyM6yHX-bBZ*TZF7X_@c2^FX^z zI3)10x`oN5zyGV|=b~gE^NK!iHCLN|&05t0M{g&>jX6nc%bWYtLEVp5D5J$@iHFvi zc3Qen7%8(wpMxo}8bK5Km6M*qTa6=!`rXx$nguo~Yqny2+;2Wr5(q)7*|)@r5Ckzl zarW=zM*!Eo7%eoD9GMoQk%e!)vVxJ5nN+hG@5lcI5W8)ml3dLkekvzP5WLxbjKV= z1|pP*2jCflS(&?rz`8t2eh)^Hj>QlqzbT%+cbQit?yA8?*-_q|zEj)p2b}tO*y_e6 z3NdL@;W{Mj|679gdHd9ZMveuTDuN%GSIW0`6NwIfRe)){KD^t{T*W0%QzuQVjB+Sx zuZL?>1*kGMOlEehuv~1Z#JeL$)ntv)(G|>ukDv$9M+e`hL1up-Vj6Y4)nT>uG1lu?e`2WP2^$5!vYCulR77g}Rp7+Ou*u zR`Ti7Z;+LT>Z{sUw_1M$ylr~Dy!j2keAVf#g!f#tsEiYxGoCRj(QEPZzGG6SA8~Xl zKY~kPf-*%v#UMR! zTv2gg&n<;NX#7)j0zjcQL3QI41uinRXcPJwPMe^R2bv1nmCL z0S)ryPG$;?^xyvVr?fPJGO(z|EvdH*LRmPbH;`*dlxW3*e`ZmZp+GPQkK_<&67Eji zE4>0~OiAytq9-4aYw-kH9z+Yc_EFsbC?QW!xCRA_*M<~mUAhZC*0<|0 zlID9_pPEXrK=IvDy>O+`4Z00Hg@@l-1kY*eYzU$pM1T_vG{`XA`(b0Rnm+EZ472E8 zsHr2KL)-zuCCS6Y^lTPkg~)OSPB;Vk$OU*E^}Lef#_;aQMYN3I17jEvpm}%*5_W<6-$h#RjdY+r=AlEI9Urb{`Xm}9QkI6 z_l^v)mn6|eCx_|ZRUam?v<&Ri63fOGSFnPbUb}zUT_i`wtkVo&jJw&XhdPy5jbycQrEqiA_te{0Et6JBIT&-I=q8B5wpx=VMkNg918$zCNO zG2Jb2@Z7jki#Oo?^vTOKbt3zlHYxd))xIE zgr%?Y7#zf48bBB@cv&FP;D?zsGvYxCYBx!?7jQ+KVM&C9{T~42KpVf>eEbMV^PzlJ z>h|Z*@O8G^MEW`bBllQ+ipl=9bq&%QZ)YLe^{_*iqHSE|x5##;{zn~B8q_}>s}7ns zC3vjloPD)QLPDCJ+Jpyh3fyF6$eQXCz(FE2!(p$|?sHj0ijj%CN-e~Y00g2)00LD# zvsOzaeOB#uyIro)^273>_iETyl(UsnchNi8yAH9+L6ej^ScJ(6`ae~mNWI!*7uwI@ zJ@i=w&v2$ee0*Y@PFFhI(^wfrW|xo&Ad(r9Mw3E!2}nSw!5a9>U&=6o0l36pUJ$=v z)`yb2x9zbXTGHG$p0Z;cPe_wgD}fLhX?8dX#^)0hPb9$Qp(I0yhWnPm^ajh(I~0vH zvJ+iH9bUhGX$5EWglw1UhSd!E z2fl~R6{PBZg|z`aLSByZX4)zDA=As~Ad!pOS1;N3|If_o`&BXWRe0FEwY!;!n%r(Q z8nx!9OTpmlI60j;7^{7hW239Jy7embJY^MiY!;58aB=J5{q4-hoeb|N*L6w9eL1as z77Qk@>Yqr}=;IKr>s+-O2QHH29xp9rVSn(^Z|SL`xKIQMpKs#x=d94>Op;)1_-V zvMODJ;I>buaLv^r3VgDE@4O0c{Xw9-N|1bt`P;q0?=6xs+WJP(=G`s+94mgkwvFHq zm0ImB&qi5)1QIViU*3Ue$JFk69+irHW~g58ce~y0lPS^GzIO|K3c%?5oQ1%9HLbQl zMyrB7G$i-#Ob^R7*IM)e^#5s*o8>k?W^j*a28~MXjdWxi+m-Oq|jL2v8F;@L> z*qS<}7hI)N*ka)$r6n5WB{1`fxN^i||1!?!5X-hGP^5I@gHGhugcBhKm_&aI&ndWI z7;0@BY7~-y00ee1Z%yafzrWr>g>j`jqNSOaiBRES%=Ni@pwp^PaqRbNzi;jvseVd- z6Cq6lkfoX^McibhS{=1r=^G30H~s!zMvtQ`KY4;2_Jr5XNd6C z{C~&O8c76#NUw(4=4{!sX3Wc_!Wqfv!2JBDWBx6@axW?Xj|?CY46tAsV(o{kznGfe z=}{Ya?RT!!-_~=703ZNB3o}n-1ZT>sEd~?P-qzH88h^&GjDJ^^NQ}C#`EJH(jd-_} zlWEZML?z*)Q^zah;w)PK+*a1cPbEVh(qkKpO*=N{%gckmQmD|TEM0_a-OofX<8_Oa z5`Z8C5Ib7DSBgeuE#t)HDLELbBWw(^CZuF!PDsevo{Htp*KSd=PjmfQCBNkoA_$}r zG~P9jgb)u^dn|oOl?_z~>*`u#T9cnFw~?cySF>1vdkggU=2ATB=98&C7+a0AaX!{? z`AjO^)~5(hS?u5ajo`V%l)R3;InH{X(TMK8=wzAVA{#@cge59Px%b@N7B?Kj$(s)qo;7+`Yh(Ef3zeTtd2y$a-H#Kko01U!q0B{m0Z6)T>~Vo%@6mt+XX z{9X6w!`L|NntjL9NK-&F-)43l5=S>pY1Hg>iyr_05J3PSf&v=`n*_C63k7=wT!Q^s zDzqf+IUEvyE!W4-L0xjX>(I7X1E@eINPFP{jCWK50KXjBt| z(sG8od8hwOKN7huNxaujlK5WTtI7A#`q~`t9-KU>$1b(oUhaOvDvl|F;D3(f&*Iq4 ziYrDhf|7?*wdrj=>I2z~um%9b1Hsl%d3d2AsXe9r{r&y@{KftL6$I5!wsj!0s(cUzhX_0&7!&EUm6(N_c&yAY9td`gCW*q4nXi4eGco|66F?V;ARc^i1M$^d z&G|O!8Yq6ZJBfP?9+@T6#E*}y!G+gSIbQXCBH6-ddQO<%_3HuyZ2wju(=TmacBuj5 zwUKjcfQ|eC03d+RI8U*1<@1GA-yR|}@~EPU)62s%Fn5k_msZvnB{nyD_IIh8vvCOU z-lKY?Bj!8r3-!fq^0YLl?Gp?y#;n#u|(oM zVVej*9!Fxmk@D^>$(Q*`Su$@Oc$3sS#^0T%DE50^WbQNsg0}gwTL1-OlADocsAQG< z22%b%Ey3P+TO~gg@#pzNm%k%CfM5(T5)Segs=b8;ro zVCyjWhIW!@%XD|Xb~W6YsR*Ro^%bqzCO^MzOkx;ug;M;jN>b$5kHogb)A#oEq)FrbA9L0vCES z2Ja;H8ta-#L2IUsl*)xP4@sd%qhZ@cda52Cbj8M`3p^exs&z4{!Dzq?V7eyYjhqlX zq5OR^oHwnX5z+II6~vo(6hpik?c+S&YMF9Ts~f8ACp-8!weAE*L`oqf=W!WwIa_ytU zO+GaIbwqr*_>{>hm5x$3l~)tz3+48_Pb_-xd&en>=b=KNZ8C=VHrJky6u{A-Z>2p+c^sVcI(B|#f zV|IjvH$x=v1CbDY#*-(@n~MDyw|^U4FdOR6|C zkKMajn##L0`JlCWPUxTF+YM4PPE*=2DsKD_8=h7D0 zOuxU*{&JylRNu65jk9K|BLYMj6+~v<;~3AUUcT*yqO9f4HFcAf^dMjZ_7H&pH{Q{j zY4tWbxsJoBYUFUfZm{c(yC}*0>ArAf-C5VxddHF(7dGw_@%kp!DdqTa2sYF`FCz*h z=I+n-(>9wjK_hu(`TE#enn${u$_iI_0$XNm=yhEeb4YL?9cpK%3u`R$sTm8FtzCUz z7PO@}u~Q2%1e8f_fciF{=AiA2A*ZQt(Xuy>N7*!rqeSj6+^Xe&?JYkisYa{V8bbv3 zl0}W>nu{4|^j8N-`}McG<14gu#>>Gt!xd~J1#;K78trGWzVNZLI!f*p~G_ipeH(sB*99E!bus#wBCZZ>r|XO6Qy zS4x%tVPaWM8cSK#92v4_i$|1lRT@@JC(lpb8*EdVQ|)KDMmp%1@732G{js`1+OXB7IDaciC5u8I z6RnEMP+OE#(0hCMZoC#!)pRM8>6`ve(J8JKlAD8rOgAolgFz6O;*$#`D%&6_*oG^u zgIj{V%`#Rx0*Sk-gN6~($Sd8t(XdZx@i}j4)9Typ+I~uHh61XI9X6RC3D<3W=Ss9$ zu_@Tz+#haKSW$5Eugf|bRNbwOzE_sy#+vi;8C$^rI8p-^!Gu*Vu|^5bgYAO9t1ULg_2 z&<&m`&-mcajj*#0+P>KWz-ahqSaP&E+P#~R=(3;1p7)g;zf_dq=YT_m3zn2c*woF$k zR4LXn8g@~It$%xz_~?{&K-|Oblx)}$O)?aQO!eW=x>uHXA!Ri<&rMqJ(=oK`(?Po8 zvUO%~W{z;t7k!^Zx6FQ^L^9@#}ggjM=CN+)>o~#S3C51Klu-eWQ|mTAwt^yni2z zIbk?f=Ca|Xe(Sr5vjv0-G&X=*!eV@_i8< zCNjxC7E3WsVRxNngHI2D-GYF^tDu~YIcw)LP>`L? z^Ss>>2$GQrY(}|PBe;fds@uJ>^5hY*cb_xY38cEFC=8|FIocVvpnS1FwQih@TyGx3 z4&!y0nM}nT?Ok%luP(6kw81QCQr((88!v)%ap;SCL})no&Lg`$QKk}Dwupk0#DaH1 z?`+vvPGK}sg8?bp4VogVvrOyZM1A5~oW-Y=qg!=|IaKpJHNs`u=n^3)&v^{1YX^3)cf<=JHq~|$ISZkoJ_Dvj~z$l>iAKc2OA=M1=wf0CW%p6J11*Meza%01b)&0syq4Nh7~uzj>^je026t z(VTPlRcU7@jftj1Y{xftZ7{E!vpA4(0~|%)^V!@Ym2;<6l3J0iW9l3S7NHtSR{!A( zpKo|eS683Y^6))$UK84XK``R|4EHZbt-pab+B{fWeFL95$YoeL+#8eb|I->=z3uY8 zYFoipcK$xiXZ&VYB3|_-%MSi*huicX&V@8u`{EsZoHTW{w+e2gq{RBaW;fNtl0G#%s=XK9s>0r0!Y_ruX(Oy-Y;f`$If{#qLfcF=uDbzZDplUR6 zee7e2-#>Ql{VAP(-cGygzDAPqV1Dgg?yBXFDRS8krZQ_Nru9%-3~ zYG>{#YdSwZ&-J@YH5e&i03L$dqCf&WQosN-NP-C*1OPo0Nd$FK0#<+|j-mtrMOhF) z0L_sEjshf-1GG(>I@B+I?lAd^Z!fPK^e z03Q+{04Ye2NzEbvAP-_75&57&B<7F+o)RWTfh3U?00UzVTh`<^Ahf<#Vtj`x<){F48Z zHylU=0ze+#wjZuzowDw3x;X#7vix*|-u=EC?II*E8_#nj9 z3BFNtLZ3VB^?&DKc8>je@~l($^=}WuD^S$-Phz1^~kV zm|m?xq4oEeeTLdO$9vAjVQK$V$(uuU*DkJg&@VgJ0Jse#sLZ7S^vL-!AG|i1XnhBsh1`bGo6DUCi z1|Wg__mv0WrBdYE|IB;)z0l+z)aE=~8>!Zeze!15B8M~_dUje^HCg%keIHA`jU)cV zbAhKF1OotJfb+f9Iji^k2{e&(K3ne}ScNR>k5fh@lec0~>0Ymto}9FL&!=G{k85LL z-xewnu7t@uQ5371tO165ux@UbFG+%9cIudwrGeL(fXh5gS<|F=q;Te}CNnEDer-1c z>9OB}9(_}IiF0c|<1lB}*Q>Lf(WfyhMChUXfMq$g-5&i&EV9rFKD(SspkM+3K>(H4 z$diwmoEgU7vDGx{i7h#6@OTVdDDG`tAtK{o*M47@B}Goxb@)oe^BB!KBg?zsbHlh# zx?MUtr|V)uj?KgI+w_!NmgR+ck@OY8aj$`fWuxwNl&&a}0T2NIL8EX~ zQKqF)8%f)k|1NokuWKi31F?fk$&y2vKj`QQr>pdSFdv3aqd52i06_rf6VKQE{;d$e zd`rK2dPG`t_DY4R!`A8DX&ZO7?l)MFnVu^y(Ww@_T)AZSyl^ zPbgU5fR=gX*%@c_fCvEu3SV#5`j;o@SRe!dfB=F30Yrma#Gye)gQQx^ zLJIZB{wWVHQJ>Gyw|;^d%crAL+)xkWdD~+I=OtkDnD7HJ!}~n(RfV}S3luZTpB~8V zL(4@fyZ&LjYgSpA{t0gGWAyZB0GNCs!Y~R1SO}|me%ttlyCKQbzXgOfuLMTbEvM)p zHl>z_Gb;PRW{?2@Ac1>_Fr4Rto!@Zy+r3qpPUv#fDqhPCBz4kz#G){JPW0_ek+`M% z5qB3DU6Qr%eyYly_2+lv=I~=@H7J{;y(wFx;-qI-)SQI$yAa%x##zR_(B6~txXqh_ z9fExq%)MZu-*p@i01bC5!@PGciKTsQS|gUfS2yPgH`?`ijGJ*^b0}#d;B?}%$awuI z+8WI!fU3|<-N-ci74m&4x^;&A77@(wx!Qm8DNAlGu3I3o(;CcAI4U1(I1D88!KjAa=d<;(aiCv@fI?^U?1kh7+^5XCwLsJiB%wNu73sD^*m=U z_MiX&-Iv55PJ+1H1flKnb#3_BuU0lIBu%~qk{#?w`&kT^ZH*qjY}CB2H^o-3V57#`$d2&R^!%jxm+^>2*&mzPCA?l9O(Z z;t&AMeGNi#73^n-kwT&Z(AzjU~#C$sm1ZaHw1-CJ*}+Jiu;^`YcLLc^Bw=m&mc z=9_{c0E7WJ2my$3%dzlR`HEy=v=GnBSU?~I5C{iJIOF+`)0Ypo5jigQJ_i(Xqz(z` z|FqQPka*S0zFQ34b*}!lf*vQkxVv8kx85^bY3U%TqBm?rpMd=foasn_q~4f6R; z%AR3tr>3{!a*$|2k{0_c<*v&u|6}rtXaN8L1P}rUAT8CSxkz~j#9UwOuKH+x6fVCz zdf#Gvh8}Tu`>tT-S7wP{?y{qs#(zhZ`ThF42F)W>le~w>(x<;%eK^NK)?Bol@@{_5 zs^4&*{`UeSW-afoU}j{IV7@FpQr+mMcx%>sUT5aJf(a6;q9j$*WIAUL^t_|=@LPRc zQl)MrcGCLA|0_zJ&09VO4}~Vnn!HP)T#i{=3m7=aEX)@aqxNRPueLRawJOSa~#|epQScQSXm$ zxqoW-8bu5nni#!m49BqmCxf^<{c%Dz))9Yd5`C%PPunLnYBKN<#zCyoGBh%V$T&1Je4|iccnq)*=m%Tgs0rfg5mry;b%{#jBhIBPx#{H5 zb&L^`LrH@=xRKD(4-JR{06_we?IL@hxLQ9`HE-L6(@^Vk+@U+q%D=|2VPsNZ`hMv~ zTC6r@JS0Z6LA003b^H%=NVLn{p|)O+Cs8F{$yPCfz>x?j$zhS_7I8 zfrBxxTPfwb0P^%W0suh)T>lc@VyEmiXh=7gEpawvP|$D0DX_BlnO9ay3*`&14I@z9 zNxd&zT$ze)4e$s62q1tEK>(eZ*``7BmhbdL005ycFOa;o5WgqfZ=0RWYEf?FUkkQ` zYxBkQGxO*SAM_+j$}bmpbk5GdS#OdtC^uJ}^(4C%`ThjrOS#dS@#lD#3@0zlqm#dV zQ(H<69&Y=k?eU@JJpY}B{^=kD0D=YYESz5wT<(tb2hghJJ##UY+H?^+^Q+ewH*sj^ z;X9nny_`d91OiATk8jtZv0(ck0!So@C_n;8B#36te~L7yuyBXTbqr%q1NeAyu#q4M zx9S9W2R@~@JQ+=~`lnrHHDh$(({Lwr)%k5K?-JkbuBTXj0Bll-yME_SC!-qxhvrtj zc@(yD1$Hg5Cj$|AOt!Lv!PyiTVgNwG0F$vgk-!K61OZMGX*dKccWe~1R8imbQ5gwt zm5*9*)_iP3O6v?YWoMt~wWN0Gb0NNi9X0L{A`tNe2lx?WSP6<$r!fhWZ;> ztkxSsaYrWjn}28fH)qYp4m?WY@uON_BI&Dl^i@_)qU1mj0tgnDet%(#^3HS+2X5z~ zTNrm9`o>iq?~hw%8Xm1>8V8fa!*v3?z1=ZqdhSJl`yq_6qaozc2}$dy(mF-Db~SB1 z3fiLYgPBu(-{9tv`5=CkbluER;bE=D*7>#knG%)Ps-jg*`C56J23rVnb7anJAP@ow z3xA8ne)8dsj<0WRxbfFPpAH?xH~SWw)ib`1eCKz4aP-xRM0TP4yh z)^WP;Baza5lKsJvcq9OUTKOT|ZO3>S)cbQ_%RF@0X5amE$hIT#z7A_~W#2tK`E-pg z{?8_|8Dp16%U|co7`?01eGc5ZKUR_0uLt(j4#z=H{Yg@XLrCfXTkfJ zpb@zt%|$lT72}xnK(8C0`#jW0~fRtwThW6c+;Y2{%odwD&p#kECWP7HA) zTzFjuU>E+jb-fUJEG2kf0dL_WY6c2G#Q-pq*((aODV^w7=+698;rm<0rk=^6m8rASx)>tOW6fA*tv1su2-(K zuddA*fmVRdZBf1hkAcH#xUlWLb06*C(*Pg@5G!@f`bF+c*PYs{%$8FAr2Xm4)WmC69KN9@W^n`*}7%14m<5q08dQE zaPXZ}RK_p%;gv2yxb(pubq^OdX4|0;%)IirCO~JnMkZUa^U_$HC#PUZvA(V+Puvs7 zmIaw7zQE?iNH!_g+)uDv+4uk;1P}=qlK<8x><@JYQHAWM77zWU`fW+l$zh>k{AF0KQFZHRW@rx zatVtIw6AET=#i|N>!B1c_)HMIm~_s9-Pd}TbCx0p4?-(>oB>^EFOaO)z3CBPe~5wH z5ron$ib;fr2DDFwaCf{;|I>Mg*rWc&A7Zqjk;z^G1daWBJ+YG5J$v0H+g%|X|9a@* z`r0@}n26j)UdL_!wEk1k=_zrb28J-S=jb_ccAZ@9#UYkHotaP1cizdqdP&)MsfZe; zq8(KdA)WSvHjBK3NbpRDYR#Y77lNmgyz9V?HsD39iAVJ98J_<}0WSvvYIQBf5WGX- zJht(M0K)^!#Zr~O?|G*g+iPgowyr2LYc6@c( zT=5p9s!Mi8dsZ?wJL0jmb(5Dg#*02j76fME_h>50C&>0J%fHr0Ul-%G?qZl42p3u` z?=^0zZ%1s^({hJrDOnKX1C-w6vgP@sEr6e|tv&1+c%_Nm9*NqapkzH*cvpfe$_ z@T-p(*?RRdsH(8)FVQllcnrY$&w_=uiXzf(J^=*_t!wgG#Ag<$c(I=D6=X zIuI*APvWSG+ksh>FdXd90ZhXk$=@ymm$>yJ|9K6@4`b1}n8oHyoT`X;@Lu-+2_^BN z>_dsrwn;7ZskcTe1RlVH4ObSAPxqynNI1o{7`rC|D6w=ajgv`PhD3)Rxfp~de@)pa zF%X<8)1)TtZygDp5^2ViYim(eNpv(qrZgwA3ImD&1OS2qpx-z3R%>4I{>mYZlp&9m zHoKFUYT2LCgR$3~%4lHF7W)rp4V(TG{j-F(#@8DP1Hpg*#onRUM*ko_s|Rsr@-pBy zU^y@VZ>O|CHNo!5k{hL^SU4(Lv6z`@khv%VX7(knSHJB8NyPchC)wU}$N0)p@u^wv z;yylH6a{&!WH-)FNY6cRsQTGe!?k!QY~GaO|}YqabymyCrApFK~jpXPaY zW}EQ;E${So{c7G0INOsyPL?N)uWDzY$Ae80|C`mda*~n9(!#NJlh@v`dtTP5D8_x} zF19eobYLk&M3LF1|m^De9`rc}W;N^6D_G=8~^|jB}cI*9ARQwSe(@XV-ZwgGKZ7eaf}IWoV4? z=$l2PZ3^fi)u%M-vhRRKe{dO_t- zkq18yGy?!(fZ?h~ArJ3)zJJ}8`!t^G-$z*?o*e;v(@8*xZ=E9-ju{4;c{so6aT(G` zug2`zc_Z!2cbh^3SWZKR*+ zWZxWBlmIJOEIq^w;SZ!ZA46mKn$t092ZHxMeIu5~Y2%T*e_ryt7^pN^qp#|ZK?D*T zdp^loW%$c?e(RAYMmH1^!9h?^!lCO+={hwNK2HjVsGEFuSFA(f8Ym`uS{Sa#IV$y3Iiwu1wAR_w$<}zELoO!R<9b%&xah$Fq%R}(|K&{`@B0xLqCdfGc?^#u=4 zpvM7in3M_%x>C@1H$Jl3xhFd-6I!-*N(h%5^9Yts1CA8^Qec6`i=+Z z7v9QCe4)}#ACW49j9B>9u;fmOmEP=pySuTjXq<8{q)XTxE?YY+7C)ygi`qVWrgqkw z#^~;1f(QT;Ba`EKbS&*!tHu03Jaqkdbt6~Q%EQ8}U~Y#WKb?__APr4t9t;2i6-xFI z-c#ZZIw7tphe@4rx#7gi-`)4$CR*n#>I01f0FntJs7S@mbp!%P z0tg^6r|{z(39SYL5gqvbQ}FX+fyxXnJ}8|TF(!VC;=4r;-A|0%cBmh+5tJPxNAndT z+i?!Lwkaq&Y2J%+j zLV4C@PRG0^TA+YSqYw-l7oJS)VYuBlpqYa{j7TQS@0VTpPrAmTDU} zJZna&gz7~PrtWOp*w`=`HCIA{G$U zsSY`=^WN6y=q_z7nv1TfBP=_54{#DIMI{r&6D;|aW_~y+IO91TgwbF{+^@ZH( zfz43-O|CT&9Mq34G4xVg6C?uXxtN<455}Q>`*B}n;$M}F7ds`tZT;>w5iMHip1j@T zxiQHdPDEyE$zz)EuY+IzqZ8#7 z$n+@Jx__zG3T6z`hE`w~MAd8O4heSJAuz$ zShreuqla?^ExG|2Y7;@})PaQJZ^t8d$_pc+*ZHuMQ*OaG}kA^IQ1Pie2;h(J# z*EK}}L{8uM*q(*VS+hzX-OW)Jfq^?EfL>E5;t)`W+P1?x zA&wJ@cgJ(+dK^3sP&QN#8R_G=aT-Chm{tH4nV+?U!cwl781gKVo@2b}ok)}kAUQHf z2DcN-W7ggwyhd7*^-ds$03$*QGTLXG1o72&IMxpoF9s5-?2nw{@B#=YS5C(RJ#vX{ zJW?=tsU^lhiJVmIkflO#R7B8LH3h=$c1HoDF(bi+5qg%0>E2qI&t;-9(;GfD@5Bc+ z6?E@1mlbq$I$PU{-}wyX@JB+HCGty-|1UwH~&2YL>?OJr6@`PCs@Z%#N-e zn7}hclxHdcfeZR74pYw%5uC^3a?7JaAQ&v+1%j+sN6t!Q0ZNL$x-p#jatZf|KaNW{ zakPua*>Guh;S)lNv^Wg-8zYQwWcH5iURUS+BrB7Gj42rU6COis;sS&P0hsORR_5Ij zY6PU4UABGe<1*S(o^^C;g&id5IfnvWSoLqfeKyd`Pwy#kw{9P=OL^Rg8PUV36ofD} z>c#i;y-Ep%LsOQRM8LtriMw5OnIG40`>0Q0@1t(8U-7>d!og?GRjLy1U*j3t{&@rf z00hacl4pD>F)`MTj@Xd}@)3K}!W7a`&9*Giy?>qwXuvKALI9kI28c5TV7Dv3u{0;i z359Mi&L3C7O7u8p#pGo3XTUM4YI>814n%Vo(*M#AwF^-YLmSu*f`DwX(~PL*rqD* zWR{o=KmeQss~=(>v@2{_k@Ddqt|*~vX^QI(0kA>Z2pk*IGaM&MeP_icYAc_R^@*9q z0kv?+Yv6BEvAd)YVQPI6OuJF{o+>1MRf?-e7{4W-4F*Gc z!d^%K%n~2pK+UVOseaZ$YL3q@&m}0CKjxXraQyra;S*@MfC5M)iYsAMY=WiJf0=WC zX|t;9Lp=64SM97ndbZzL8%DzVZxcEBt{{q8fAyS$uw4)L_1*2|UsfV*E7w56;o)nC zo)GB$VF9tJP0)qA_hHv9Y6s#37dG%M#Pmg-VyZ!8gpa=at$4M_jHItD}wu}AjUpI_Q9Yt$)F zVq$a#q3J$I;QrHupXELhJ9*4jSwzWXu|KWonrZQ_+pAW;_OqF3bZI zY!Y4wzX4~`7cN=8_>7>f!Y~p0%2TlPh?Os77F%0aE;lDgcuYPKUB6Br-H!CQ4mz=n zdkx;DqDe8<2aSNl6iEO8N*U4tBoasOXu6OH022&uwp8aLk^UA-i)uJibprv^I|d= z*l>Fmv8??2xN;tRGD4j&c0XU`nR1$3NBf30PwQTpe!9= z|I2X})6Zk}Y#5|Rg3V$mQ)Iw_Iu8>u1+&CB%242gE*fP!Zx&;!!bpXeC2??97>+naDSXM9JNX zy@0YQEr9Vt4=8+rPzpG951`AUe{=|gp(C7br~ZkV5j@fx71Fl_d#>Mdk*-QJI={;9 zvIDV|YbFE%CM1A<_lp09zX#O^`*~r^sM&qVZ;H3$k0263Bt~PE6%S$pW#82OZ*UK{ zi8bQ4BovJFqwCv?ITvdRAxNM%;WRzKxD>aHqaa8K0h7rW*QF0axHmlxD+WO+6^efv zZnwr8zpsae-0G@XNmOCVVDbGaRlQrY(M@4he7v_v5mOCJFudx5KbB1Jb>+s^s*P1N z>?|+}K;dfWukm^w=FI!EdPbueTrD8Vhc^MYaopz^%aUEREGrbJ_NJ&rVTAOx%G-d2 z{i|))Hv5GVw14t?2bJ>2ZMb%gNBbEKMzMB%v1};!*BN%Lq?Uw9002t~op_3Q0FntI z*-o8~w|>^pa4~AiSK?qLa&BHCKa~yD3~$wq_cxFav+9CI_26m*hEt7runC-#D8&8S zf|1B!$Ozw-9=7}k?7)F7pG01rSD0kfzc|wu=p!)%8fhYcslZSiF^B`Z2XG8| zeM5sHt+SC<=C1=XW8lF2oI?W2VUM-=SG;~^*n0oV+%cvGkz<(e>ac&Ru}aXmWAE)u zY|nD30I;L(E>2V4-qT<>8&eFb_5cU~6LOHwzKAQU8<4a(daV<4=NJu-p8MXo@*{^Y z+M{%}aCXzOerIB-4M9|IqTd7HkM!pN9|0r-W`zF#+(JkPrEJ^dB0(f! z{~OKvCbpd%L2W~yB(r@B!Cl70>MRNw8OLeUleb{})@{a^2ZJ$kkNfKpA`c*NQbh$L zvPbQZh@1AJz~Y{ob4S1jyq0$~O2kR7c9`&)e@a1Qwoqj@CPCwX&=^{Q*^m$3CbK-m>~O4w{2i-gclABy|6Ua|K+})TAFA zWZPY-f#cMzbat=!+XTO~BB5Xl=kSEOf^fs?Ew$3(EtKj79YgPKMS)Rkv05eD5;uG) zoR4+y(3VyW=-?VOHfd2dirtUhnTJfl5WLvgceSrF%B4BF<;mXMGoT)YXd_GpKCz7Y!y)0yhuC;-%fgAwN_Fe(e3#7jva zKOq^)W0rMG zF2y@{R`Xkg_`|CpD!`!as&7X{J}GH-)IDeE%yg<~!`)PxZ-RL75omX{snha`M2AjG zaKXGsTpLuSk3fo3fc7aqyI|v>py(`^$ov#cYCdf2lI=h2`;i6$uG$X&DcKSMBoaib zQt1^rm%;~yQBR**3td(2$O7ctU9D0h3Hwo6uSlPor6ucW;9?F3hw{}zU_3y70w`q4 z#HtotE4lF)&K03+XVe!+zIome4_n zC4?^>l_yiL*}mlX+Dra}{-*u-h}Nh}QFdjI$#zY;y?u!~vyIG_d;6~chldb9k?KAm zrH#}bANaO?e)z05cw0j6lN$g)-JpX|CE=i4=5_451wJdexN1VUWY|zFS1*rs0K(ew z&K<5@TMIUT{Rl}`Syz>|n%Xrj=4eM+O~YO%VXsD0hxS5gg6wWJ$hs z)onqCxdHp%UPFVEMo{W_l3iZ2c)58Q`s9kcW)Ed z9;Vj}Z-Nd5v+0v;U?VwXz&rnBCHie8emm>wP=5`DazUe9MPJ;o`fh*%00jNDN4Nk$ z0GC>%cbH-Nb9MYf#)73g!7v5nSy!H((@e_=mNLQ8ez{51`Dt|wCpQfXuiVDTQkXIm zk}NC-WqjD1T=P2z_1&)fZ!pXn3J&qq1eFKt90xh&2b9%C9W3|3%_DpG5CEe>A(0?D z4*r$WgH%#Jk4zX%L0{=D4tK7=ev62`%_oAxZO`j4m;MK?c^m#wXxiKk0!mK>eN|hd zg?Smu)uGh%)_n4401XqL0j!mgCtB7pVqjFW1N1n-6)6_)wBhh&XM!gNm_R2mEL`bo zhkMDA0&ywK`;Q)#TIV;Z4fJKpJ7s|6mM5UYtJi7@t=Khr*Pp3n(8E!#&Z8@8foYWY zZo`KYWgyr=^6mZzn0N)*3O}x9;UhxHK>ubJ9hBo0=hM*F*36SI1P~YwE-_PH0YhJ7 zh=wsulKsuzJ6PK{YonNx`4QvVp!(sbuT8|-h$uBTK0a&xS?IN{KQ17TaN&wCB-MMvLei+NiGm!ULrW< z6&iET2#5kj4^Y6^3>Zk@y^Wq-%|mEc?TCz>Ja>=|RyAAy%sC&c9g4c4GQH~K-|4?( z95#`ER|!Z>=|lw6!4MLm)eEa13Fv@97S9h-l9h@^@*1)^gqF(+&$NoK=7*0#E#5CQ zMJ?tABNUpZCMotwge5yTI%AtZvY;~bC7&LfR9?|YKxgU+`DY2C{0VH#AP5Fwb5Qwx zGmCZ~Ne`029jLHhVICxBJ9wDbyQ=N7blzO|B6%Rt?IaFg3sw;4;P9W=QTHK-2o{5} znp~m#le8de8wL?XR}cGUYv#)oNJrSLGHEEHS4Ajj4glmCEu-|6$8I#R3GIQS5X9gb zsUW4fSQt+bh_|%N&7t(2Q~*W3bL!v;N2>b=Y^{;RZ)_Gh5nDBIE4smtJ3nu%3DIg4Xr< z_T1gBj!$fvF#S-A^e}$yXPgGTTF~*Kf0s7B3Nj+nqfa%X#uuBaD*)j^0|YMruw40a z11Qk6=?>oRQBF8iqBQH|sFTcwiqdFlX(fF7A@>2<-<#-1o%OaK0?}y9G>Xu_1-Gj- z=yUz|Ta>9|q}^@&x??(tI78d(Hv8|SIEr9WK2ADN^Tw^T%I&Y<{2tCspYrtx-|vpG z;B}db2|h_WY(6_T`0;dcVgpODN3hg^zwS2M{+YwPyGDF-jU}ES zmM&^*T-=BbWSei>Srg#5_tEID)xTbGl)J%PK1k}&NAcRRNy#ec_X>dY<+#EPlnT%fUlg( zL5TMPs#~$7j#EW}e&m?|NkF#0k$f|Ch7N7n4g+)LMrhjkf!>!CxrsPe-_<6O#K$z* z?b-*4&YO|lgX1+EP<9&d9&R+4_w;k+Tze8f>E2%F`S)$I`lviAU`Or$A{z|)_*x9F z`-8Io&AYdgjTZ>@2WV}_f%SDYUFK%Az`5Eg#kfta>ghpx{{5;+{63drhjo~C8xotA zP5h;dyM^_eenII}KRYU#?*shZ=H4ql(OOYaee~#85&g6E%^qL;u&(OPef!wN!b;D{ zb-MqjX$Y4T#FbK~ly_0fofLdEFEPj(_xuxgx{#KvPA_v+3{<%3spce9TCTwjqO?Hu zRu&8MmAPHehuAV>Kh%&R}?;Pt4& zNoZk-_i}1#yFv7G^`!d;3vJMO>-z5R>(Wxxn^3a@z^6*gN@Jn;Y2s5q7_vHkDv69T z>zA!o>&d@rhj!JEyPU+mO+EBOe`m&lhxA@?27GEOPA%h&fz0|RnCJ``eqUSxj~z$y zp#S?c{9WIv0Gt|c9a2pBFF4r&2H;~uuz)B6Wt1_)->MC~d|oDK%MD_*6aMI%hq=;H zG|^+-V?2?<)1{m=AO9f#24ZIy49HYTSB^2;GlH`yprBczAmDboW%)cF{>!zOde0dL0(jv+Vm6H&j-Cz=WQ#piC* zzuf4O`kxI16MD0)VBX-nu18QSK8^L;SqdrW=j;n9V^(rgwGdJNBKMNbtz!21WnGVF z{V|Y^xl*7djIdW>P}&oBy`80kiPW$L!5IaSG974;T{^2;-4f~)5U|Tla~iABK0C6h z*+x7N>BSb;`=n{vdz#_>xtB6-tH`r>ozLD=sY&UlVj-m5(9>Gy&tw<#upK=7Dv2K) z#1Jn-l-rK#27AzWFguBRJOkT9=br)qDEi%BH$bMGk_PpaW}8L7yhy|__rY79`zMh2d~V$yK|E3_81>@`HW>C4;OL15CHy`oVuO%IvP(pap-HDKq?d< z07HTt3LY&($M6$Ng6O@j#y_ht(zV76KBj~MsJPa>KWE)ZyigKBhhw?J-qYE>P_!b7AvVGc_3Z*-1EP2JU%;P2H7wU&~~}iyuV`uUV6E{v=+YbM2R7ZjFn5 zTai(t(KpSpk`Ztjtfu(Mz=Ti9qvrP7&OO&);d^p7WYhn&Jq9wnY{vw6N(^bf-iy7r zXD6s>lL-Vw@a(l8rA=<3F+Ml9e<;hJMgg$A9W^53Mw?$WMP9jr6{J=6j6o+M1`?U+ zs^t3J7iv`VRpP>D6uY58QdC&^vhaaGO(M47SWxck$!YUh&+ zGW`b9*UG6mN5#yWKV>D+ct8=U!UpZ8T|PLIStbt?E3QXb)cqBRCaCS^kFOaTz{xFA zE_cVl6N{;S2lDQ*(^339A^|)^D7ih-WgSwbxT!sUB77?k4grw;t3QCsX4aNcywXHt z@qVr6ZWI>XMwCnLXyScxuM??|@O|ZW;>Q)rWKyYEWBbb|PP1Qsa3yXqUX_@VsrEnd z8?SGU|IPm?kfojblvi)HoOY0ZkD01^b)}u^w?s9fQ(lG`9}2be;+QgCVez{%$?w>*AB#An?A=wYaFhRM(Y2r;&6hBjs?C9B?sq}rXpN9=suKsDk zhxsbqNNSBf(VksH=0V5TV4wK8(=Am~gqFkryW$9-2t*7dI+rOUFG!%m!hE1e&xrP4 zNb$Zq^Kdh2cO(GcN7t?Fw>!g45B)=bh7oE;2>Actw7;BqHgax-O46`rj2JnES&?<& z?hAzLwRh3Cu-nn8t%X&G(+-@?Lb#$t&3`#lOPK26(T>aKeYu?fL261;semRkXfY0n z2`WP^v5sD=k2%0ylZ?S*)<69#O{09Ro`^iOciN@73d;ZLV(&HUI+>C!G@>c|$zO4m ztCyI!5uYnI@APL9T7HG?z4Ye33TFQwm^2ceoBqhR#gQ{_rZOm*9n-65sQ1*003Y4f zyL{AWC_Bm3vZi!|Ho2lB~R;X)}8otVASRINhc=whPa8rMy zJba(QQW{w`l}bX9K~(M1-dyJ6?%lU%$Z!K-4PmkUz;c+H4 zB>$yL=Q^n%9WGhcD%E;9iB`sF-3w2;p4H|Z+bp?2hcO5cBgq2Srou5sH`UuF>8Am3 zlk0|Te|k6=z!6l-vZq|74j;BOBa-)38BJ@|ZvDnihl2#cMF{O~RlD{; z4|ICy)6gczFNb~R=fA#Xa3>>GB6KeyxiO21X8RohGuT@7Gczg**1CMW%PD%h&QP(9 z(Ho(4jIOM$SSWD6@iugIbV4=?1WkkxM*ZH^sWIlbd~54}<&rM^>uL4gX5M!@p_9#1 z`24KUhQP8G2lc^zRAxCWNWpzyS?Q?Fu+dR`^r3j_57_a|3^^YcnNM`VOE#@r54d;^ zC3o7J-~7fLojtQWn5|W-BG>xo-|@7BtQV~#q)+z67Lt*7ikF*p8uuhU-}GInrGQuj zk8*-y1i5fNc2!J=S@pj}em&D0kmB%5XvPUY6bjWbBmT+*+xaB)p!?O0}2OR_z<|+iFkZ@u_Hp z{AS-dHBT(PI=#xD85g$p+$7#z6R_JM(W?fytxTThU40Nywk{=MKf8u zFfl+p4?#8a?pby{UHkbWbg-8{ms{pjTL}we;z@d*b8rx%sY?})SN~<&dU1_w9ehdj zMbk@M`f!d6>L0_b)ND70QB#t@o}_pufu$BfHH_ zsubRvX49^-ndenw4n5Jem|p4)%v7?~CK9Dvu{W2N`CkvwJFh(B>$M!D1L1@QFa{U* zgXG&|82Lf!#x?n2R~Pi;4;%OK2CLjc6Bkvy!m6ATx=h>zVwZ5T4M&4l z?%PoxGBx>c9sft-QK3F{P*WHA%(G>`gC4z5T zKHOrd{ZR!q2G%P+znwcxu=Z-w+TmN8YObS5ym+;^urb8nMHK~X+FCeMw% z_A8k9->QEi{AO&qt!TXK;gu#}&57G~z+$q^+W52LXcaQ;?Hf3D`X7*jn0ssEEj0NVkrSLu&ZMiwxHC&GWIYgc?jEGK;UxB{a1sF4H{ zKf`(1@|2N`Nf@8p%C8bj?ji9Gbujkujss|`n3XqWD!QG@PuJ3Kxn+O|vT2n3bR?&mx$>CX6sNzJh~ApTx_eU^VD*ZA|79@VOt-{$IZK45Y|PVP509Yw&~ti2>mn z*)uR@*;Mi+Z2+VwsCtjT-zi{jCKONvB;&4px+qWP6ojO`37%g?ebF`#v`s=1dy!8w zkH!O)F=j2Z2}rgt{8&VeFo2#CK3!9R(}HeL4D4)#n9b1^ZiN{bvd1i+^o_m>FyPJ! zsR3BX5gShlRA(@fv5NRGw9P7q)k!n@!{d_PR1xC(`L77{`B9tK@~69>Tl5{yxw<2r zSebE!Ippb_5C}ZTu5qY7K5{YxnAMe2i0QTBV-c$oh}KNW~q@G_$yn zU?86?oKdk(LrTrV)o%E7tXI4Dc4H?b_oHQu504XTuaK?vS!O3h`nGsVfDP>up@IK)j)xjNm0Z#nUT={ImXqlzYkRIQ7P;*Aei?A7&Tu^Goe3C6@lSC0?<^p;#c<;z(@-vTaoZSL2&*9g05T{H zn@@k}crgZv?{PuRj2Z#Ar&Il z>Eni*I~R{K+rmO-k0(}^ZZ5%*Mw?8w;HU!oT5*~*m`4uT`#E%vgHy+3hK`@QDU-({ zt}uaBf*6AOp_LCueqwxBArPX&AhiKj($ z)tf-JOEB@E)4Xu=9Dn1pf8)#I{YpO^+g%ODgAUtC9T|zacd%QZ`h6cCsdqh-F=?9? zvp^*Drj5JfWL^KZ5ATB(HbZ}Kn7??p5A0~V|Mwete*NQPJ&h^hc%Z#uEw2A9F^)wG z2jtfcGV93c_*09b|<0QtPzWL#wE--DQ&9fz!?Sg;Sc^-2)bmuo3Cd1b`%<>mw zDF<gELrOt6JAX2@Kus=p=(A1SiZ*Xb?x5TQ` z+m?WitEaSE0#^=et@8cDgLk5H>=}&JO8t9{ijqpwNTKI)|9laP8SgQSQkBiGcFG`Y zILH~U!3|ur@9Z1|P0?z(&lIa|1?FUdWg;!+tg~@PgpdYAp!aUZqV6P;cZ;tb_1YH+ z#@G>dHbiL|F^K1}1;qm_YwM2xcvU5)?YA;T#wcN3ln@J$gbYj7>;b^I4@P%onWvOZ zP+tDb%b$xQxHJDEBCfBk!i5K-?OPZ?F=*fihKLeW%JA>aB#yN%xI^2RhCT^W-syOft@m_^>NS*Hf(ICFAl9oi4Eyl)g(XU0M z$fS{v2zOd&=})}(3j2C7kNbXVwwG$!n8OhTr%aIlDH4d+%^DV127DmeZr{dfm%W5( zAjs@`%m8j5bg`)Ekzb<%0?H%ZY2*Zl_|oOmR`UVHqS}NHYa*o23(}(^QT1$Cq#Y3p zDoT9c@G3bA2`!M^lO3&CCY-CRa^7dIl%xh1B_jULv;SO|;>{t>A>r09f$yjGf8Yi8 z|IY&b(s|PQ-cRfG1a@EMME)U>anI`vfB>jizFxYiNqAe@db0I#28<&LxerLDSAg*2 z^8r83dn_ymI(>ZX)SBuk+k1@hp!LyJw~ffzUXRTU-3hdnahY*Z)PB=zH#3dBipJFI z4jxohg>-o@+h=r--AHHwZ)Du!rH)7i)&6e5aT5QjRdSN~K;Y`Ch#&^uw1OXY;7A5I z-=NDGp3N>|J0EFXm4Z z`Iq;|4my6Ne@Z-pwB?Z|07IY=dU2QoJ2aCGY)YkqDC+R~9ZnMLFn$>gKk?}~c|fQ> z&-$)IfTTrL7+UdF!*(4)Bj$Re+=mNx>F@5{X-RjBNQM`@ z>vM}T#9uCLL9%8zs&sn0;NsXEQGm*(EQ5(gDB(MI7UmrTzR_VXg1k^3TATb*zmcNl zYcoy$W|2~gVy-<_X$|8^5+LLF8pkJqYOggMeS>V!jIxas7A(K=<-6^6Jzr=rI*&?g z1&^zk*Jp_7Vz8D?ORPJrgi)LbLA-L@cOopd|2O3+8fk;F^N=%3JKiD9n z-N!vFC8|GggdmeW6&8)qMb>vR&GdLs865WYru-(2Bzw|(HyG)=>G`B?`d4cX^@gSI zUYpx%rF!`TEB%4Ke~Dw`pHo)9@^;WT=K40aB`9(e?&WpwDOT!jJLmHHH!?IG$?py+ z+H`VQ^6na4uFA)ld@qvRh1Irx(>`jxk1nyMM$`GM2CijwFJU?GN9abI30)z(d}}Ma zVUYRvM5X_Mz6nU%3%kv}@!#Hx3Wrn%~k%|*Eqat<}FFv)+aS@-e)2oZ4zGni~})3 z;PRs5-w6NMUhWm}1@=X6=iqI{b)pI!zO?awUMOu@ZF93qA`v?XvhXdHk@vZEq8HM6 z&0LmYAwK%SWg-azc%*AQ9jOH3)As1dz;=@)41U5v_0O^uXWk+i*C+pPp@1w0N!Q1X zQLsrQhVYYwBe(Qwc&_F!2B`0RgaT!*~#4 zL<^J=uG1<+s z2!yzp7ON+L>&x3+NJp0rMP8mh|8K(s?R|L;E$@JEXWs=X>vzt4F`_e6+Rsz}O^t@d zhvmi?fl?K^m-;(cp-MRhA9wRu`LersVrx3`2@d}CVTr;~hg!^#&8FM>GSQf*A=b2? z|1^b0Ldj^S{HN031LAkh;d_ONquj!!n{$xTFzGw7_P)(7fH)=nEh6VUu;MV+I@p7s zfv3sdcxz+CZERr3J6PhD$`5nQ_db z@%dxf$V_g059ZebD>=mTKd~5P9%8!#4Jqbbf6m7>2jh|=rz&7zX_XMz94KA`=?Vbh zXbXPf$>Fo}ja5S5`ZNMWy8rL}0D4Jjd$KT26K8LU$ltb+7}XV+YExZW*W6%Urz{z8 zwf;v?mZ@VF>P;D*pusd#Wgvx?u{2flr>&Z$%MIi|(stW8-mDu7O7=u$8%A)+)AVyVf&`bm!_I{9 z=HQ(zY?hhx6ehmgS&-KAnp+!raBi*bygj%&xc2Eo)LiMkoA;$iBL$i+{Wu}zTFS#p zU>^!{b-A4=RGA(kM+k$(c#pcv=9jpmFBmM8rl`YfaHSy?1XO`SojYdktb%Zyb@IKa zGB=-YRmHSquzswh+Z#GWp(T*&klrb??>lGt(4MU1_8|7v@d(d`C#LP;{kwF`{>wvi z@c!NG1?M%il1b8@u+*v-^MK~OLlf49zM%uQ(bE?Y@Dhr-7J751tCK+lG*l@z=5Xak zW#06?#Yj#RMjHMb2h{n_aHTh^MAm2+VlfOhZ^g-(1uy)bbc^y{_HJH-sqs?wG>Z-4i6MJJ^U4(@TJS&csi5p2ISfNtBg673>R8}r8~{`jG@eT^!+ z%q~p%7}qtDj7C5~-mw$$K?GojXASUh1Mf}Srh8hA9y<-~(__CoM3L8~lk%SLeW=OO zMl#eT*(f@wb~~*Jz0@72e4)a#M{g7quYTKAgvOaZ_Zbt1_+F_1% z+lTUQc;CeoRjq}HkkQvsM}Uwxx<77lLY`)}l&$n*tGHZ6e|0zc^rLdt`K;A=pDM-JJdaCjtm=~> z=0Oq`jdN0ybGDp0q*!MPn-nshp@~jQ-NUtbno1V;f4KJjxC|qH&l?B4NnPOn+j~&E z|FJ{^W51)ix27!&dfF37B0@|JI;D3D>5zbtfQ|#^F4hQCd)q< zu*vFq^s}gKQ{USl%+g_xd1?YIH0g@5uTz-{8L}>$TTzcE~ z(_`Uq?BZQ2eosS+w^$-*JWUjCFYRi3QZ^fB;^5~wstS;^*tRr8Zzx!;_%|G#R?{Hm zs0Iz1rzF0)^?CE0Kjg`BR&6-ua8w?1;7f})!#VX9x4DN|fWL10=0+No3>A(2?4Yd0~X2rNNDYDVV)rsYHuUH#> zZx>$RRja+eCMqRQm*uiasBpRl8)vtzl24h00{uVi+lQ+ zEUc{-@q&zqAWUnD+kKpV!Qv7gS8Y~X>W9U${TztLA?dwrN<2$!!= z{$QJ!$&LzDp|FE#C;ncu+xs4HXS{Fc+F_H*zW z2+h(N?n4cp35R;wEj=)OO2kjaO<0XGHr5q48NzLjk~7ClTw3k8KhPjF({Ak$PPAa; zfiJM7D_}*c)bSm!z_X1x%`P6Qk2BF7Rf?gh!XqG7L-nh27(!VSqz(xhsNzt0f+jiX zEP&aQ#IC z*mVK|@|ziZ=VC*L$aoo&nR;u7D&O@mU9)a^b_;QL=f)rcU}bMn1~d?Q9=d0b(_Mko zl*^rgR6=n56?(!ca0LiJDG9R+s4L0@A>hI6pT`c29b7hmUg$oKRKijT{q*=q&y(m)M-3XV+duy7?Pf-=q^tl6(qRM;8<*etFO?zcXRe#6xKa8-1SF>2w~YK;%2`)Od=s$#QCD$cv0 zGnbnqs*R%L1Ck|Fg@NqDcST4)lOF>Oi__kbe4=)>(i&=!YB@p&m8xp6l=PRCXVS53b6aq zM-E!>qA<61I^?@~hf|Q5Pgo!M8M^K*A{XQYMT6<*_zxqvg?{HNvEX0^k@HGwH7E^N zbk#1@miYhm0dZWLH|hQT7!U7vS+T}}^;@O&d8F>?6T-cv&(5*=P}WnM;Y0CN4z=TD z&600oMqLftbwX`7gh}l`6uo7uP~=gXi=T&&z7#rxJXPS7RX8bNUi|{UTq`Y`McGvzDFSWzwuc zyXw@-yUNR&($GG~*%i#y|6oj_P4R|Ql`ZT1w|5gS;b0Qd))Xl@$z8t366(T)XA3D) zG10884!#5O%}4V5dX37?%CFnmmQ4$}qA$W<{1HSHHh~bxz^q^)m@0Ga??WcDPki`} zSrCY7KOy7$T1903RekF^O;jH3qh8lAJ;Pg8KjMV1ij_Jkns*y567hu*<7poLgsak$ z%8f!7jjwJbQA2sTJ_;^14S#t4k^8e~bCvW|**cEkAYd?6bzE~enJ`jK6SAj0TGz9w zziI~iHfZ-TN4~TLaEt&^>&NF?uk9$ps8J)m@ZjM)JtEDZ5_Y?rZml@oHY$p&&b!^g z2~XOSVrR!0BbZuO%ki*R@L+r;VTf8(MQvMo6_~o3D-o8Gr+-tTUm23UO^esodQRiJ zlV$efWueHYOxjtOnHky$4;D|zp@ig&)EIi1YaHh8WG&?b z`l{s{K2LRA(QDhSBiKr}xwE^eh~r((_m)vGogQAnKKal1?wPB%&NX)SvkefzdqW5# z4DsknJb0n0hO9w9_9ZA}2G~3?wTc&&i~hxFZ9ASAj>KViLXwivkh;CBb8JR~ZcVc| z)cdr1>&eqpvTDGQUWtDGRn~GajgDsJFvLs{?KnyIFJ-DrL+Ru6m;Sp3P5vwNv-jx8 z$u5wu?Z({xLnXP64AOLQ+GwPl|4aQuURLM3efDK*Q2)J>@9q>bl+styxpi6&eUE+7 zroJ}XF)c2{UeETiUT_n?hjV#2-jD_z=^bx!?E`YZCr9+2QVo!RP%(s20uYHryjX-? zRP}7f=CiG9IXaYlzxS%sXJcmCcO@LuIWbo;WoL)wVy3|3NZEbQ4B+cOPM%2hqp*xR}1YQD04Dvgn^9{$Z!$(sh4nGcIp#f~L!$hO0AhR#9!Wyk$o@r&grt z@O05@W3|hT%Iwqi>+$QB;W|JEsZQEWk@TN#J_aft^z8$Humc1DqP7U2JCXPP9<1GW zuSzhIe&{RARF(~h*oiDO*Ezyt^VXB4IMI;i@jqN+! z(O&2JsUIr-l3%ws^Gx4bsynJ<*Oj#FI{UHKdpEc4ckDhKX{5Vcl3Kl$u^cc05>hKQ zHNC98GOgVgom5ooe*qtMYkqZZxq{`Jndz@pX4&ExgbE^n0=H)sTjO!{iY;*QoKF~$ zPIT^vm)zXDp;6ZQ8vn|xQ9X$oEEP?lx5B5X0yFqJ^Wg#2@7^Fypuy_ys>C$Jv^<)0 zRi|M2TfKAH+F(A)zDamCaQMC4j-#3jdf6#&-gm0QyVhLP&3(%{{O(1JS)I~{q^Szv zJ7EBxyxRqZD2n4XV~T3B_T)C&SKr0=lJDHA{koylOr&=$v&8RWz@b@wz{=#DqcTpT zp9Os(!A|SYX0$e0DqxR`_bQp}o60g`fBIb=1-<1*%18-GheHF{cE`A;GGm*Ll_fLKEAFlHs|h$0sc43jh`mr7E|u(G!_bO3ICNSsx*oeA&tC$ zvied2u?OA*|Cmdt15BqyY3Q>k{ zJjprg<}Mz3s}BmufD(OXBlc@sFv9$WwB-YAs%d!|Km2MsROY_{XLemO|)?zfraw?J#*SQU?T)f#k zrm)K+Ai-g69U(JzJIszX`gv^>*16#&8RciZHTqv!C^a+9 z(^B!*DbQ)!c*?w%emoIsu3=Xl^V>(Jvx_#vITond5v7u~C)a^{O=0Y-8dSkuoPu;O zUS<5J9$j}b`K1Vw*7l|q4Jh6cTV6%`VIB#xzX~G{_dn>naE?C4WO_T|9KEFJt$)P7H zvnD+!Cj(}mlRNLGHCi)i)IlXqU;aKh32f#U&}IJ*FZ{FT!90?P3ZGFfqL(x8m~SiE z%g*;J=z+1Q{pwmSH6+W>oj$_tqtMf8)w~{2%*yKyrDQpql==H!u`{@`_xXl$?#$W@ zKSJx2l%^_uc6ZZr|L*fr=F(3`uNNVITZHScFP~EL?bQ7>ImUdtnh(C&#m?lk);3`M zPO`4V>J3T=`c!FhGTnxHm1lZSC80yx&EjWMNfB`K(!EpqcJbHYGf&? z`BcQyx7A8IFLxl`?zyQfmYphc#YV^su`VHcy3vR?4m{$BGhnw zachRtJ{!mU7|&Aeo|0Xv?gpio*%NB3UDj`{u%fQK!{uXXzClG<@aSG2{BoaT-*fZ) zzXh}r$4(Uy3@TE}Id<)yUuW5BGwR8> zy1k!YiD%>2|Gg&mC;HPJ;~BM|U8(5NK^Pzi*Td_uK%`bhIsJ`Omu-A6*RxDCo3^Ui zy+q1f?DYTo##L^yy&c8wdu@2R{SuUT=pSQxt}ZV%K}uL2wwun(r@!Sl9G9(Py(gEx zrP{a1_L?==Vd(XoY}aqs9y*AC4f^750halKKU_VlpdOz1aMf|p&UU!xe>Uq&z?!>Q*q>~C)zLOUg zd%{^;_v3VRCs|rJYqm^G`a5a-4aoC;b`-=tcPWDrM^HG0-yg7S|LR8DVSP=dfdI?3nHXp8gd`4}5+5BvA!_jt%L?inOsXrMdx=@qDLILm+2($t@aLA33d|6SYBU2XqS zY^(ye8Dn2#E?ZV!l2F@)$Jx8tsQ&AZ@aFf}7I2R+bFgK;!RVSs)uzCn_s|gCbt!*u z^i0|AL}|24llTW#-0y%1Sd4Xt6n~X7B9C)Q_NRLdv&L ze_`^@rG~+xgJgjP;hUx(*c-RIz$b+$YO^E&6jJTh8`RlV?}g6IGDIG8V~I z{eD(gt`U&K4|uZbN(uK9J%a+JPMIWvKflq}hwpqra#a|hm{qq=5MhQGQw->TQVWBaX&nnT zrNZ3j5Tfu7Oc{e&C{ggMd6i0KW;t^d4d=ExhUn-@28!yxtSU~I_0PG=^Y-t#Ja+rN zAo+`uh5!KxrL{W+`pT=Rm}zmHj&am}uJ2{$rP#w^_xYMuzA6;Uf3>A3^Q`z`~R=OiyNn57dGWIc^-Ic4NsMdKf#j%i@!z6tfVuA#+R4pi?h+5+Ryl* z9`-(&vlFJ>OE>a4iQkq}y&6?J_BLh-xz>exYy!{F?7TY)yB}I)g4M|X`MT-z#|hud zZB%{kIxwG^7DB%+Ms7<)`QAieGWI8R?&Ixnh>8#a%{f*CteF6*y#ckHuCB#*U8c?E z5;KVEwV}N)Us}NP=+yR+_PjAD|1Pkt*Z0=idZ`&KcvtbMQHbkzb@KbmEh&_Tv-oBS z;1EL|qMC2L`8(gmVJjlCfm-+GZsDtJwT#xTg?5L5?4>iK8dtV`Zgj=>pG0@tUGh7p z`Atf<_yl`!IQCa~X8(3XtWa5#pIuer{{=ipDfw$M<{Hl}Bh=|UkcyqQol{HcPd5Gu z;yQIHfv1z@79`ak&rTmmej$SPH#*5%9&8s0!@}r^*PGj z{3$(e8A*H6E4a7<0MCp407e*IZ9DSUPE$m?Z>Ogs=ylT|*Z$F*Gp1ovx~nNsdMiQ2 zGp_yPEwm&UWpZulm0>sdJrT9?3z}Ac;il`3rD4 zb1dm^{NGiXf1M?t>%+4*Nf0?pp3PCEVqGPIuI^kzdA+ZmWfEED2#X!S3J?Gx#} z(eK)8Qyn{{H!6%U0D8yjWcbp5&fD92v-C&yYK^jS)~`_g9UEcqgTppuy4LJ!mwv|+ zv_EF?MOO+sX@*e=49<&X{eh~7V$(*Pe^WMH2P%?{T8AA{9v1@5rkYj71`mJf^t7w> zF#oKDj=#S|nVe@U`ZaH`OeO2BFiOB*KTw zTS=&?&)Gbvr=R^``YQM+UB~?$uIMrm(mcsZrt-!}vrT_SPl`ffd=Ni#qX1!pA=%K^ z+-B*iwoCuGK5l-poC^JhRj#z(tKea(FWDUx;eMO5sVKb&vj&pDwEV&Vf|d5eGial! zE@nG{DV2YZ_U22+xOZVC45tt=yinVEEADttYSvH6F1S=)nU5yNci|T2>t8ot&>rfH zx4M_ZM`VjEZjH@9fSSU1%v$icX}A_EQb2<#B}{f%@5yNyjrZb-v4tY2@%DkE$lapB zXd z+`M;NAg{A?{@F@bR(xbJH1K__*VIhLKCXc^iLr^@O3_ktF=XOVCfobIb$rexRjWEv zey!sSnuWXyt;znSii5LEm>$iIK zI>&^A^YO)RU<{MUHN{*MM@}Oid+9c zQh0CLAc0!igjv;#jElO{5f69R&drRf2PEr!BllAMT#wXoj5<-TZo$K4v!c()W$Gzx ziV{O2Nhqx*Y>m!3qpD31@K!c?ANnf2y&j;a7yz6#iou2ZJ7}6$zfADG=ovSfw|I*Y zrnX-EUT3_##VOwF@oHhqw@It6FUxv`moF;%9)}MAn;wg>&i38v7mD#q`$)cvH-G5a z=cl6_Zl*bO+6vf9=(6mJj3uqT22A8}m1ZfR%B)EsSD&F5gZX2KKfU!DYF1MrWDU&B5C8+DhUxoF%WXKuTfKb`_iFzBPWe zetXuR@o4uMkfB>256x3MW8g=rZQ|XrE3-oPWo(U{P_wLgCwkR}g1`^u8vL?kN$O## z+va;d8uE_f5tlKdDbFq!r~y3mf(r&36T-TjEg&jYulld2e?tETo(2VtvhPCbGyW|p|D0KQ-QtAk^#$6urs4;J28s+|_M;WVX zI@hWHoV3i_G_?;5-^I!^77dy1_WSF(f%+^*EjrRWBwAOOxU$d|v+tMgHcZuaKQ6Fs zQkF)aH4JW|H&A><6Uh1K|2Y3;_1o$Yg76vo`udSf^2C#&WXSF*V14D-aSPbMVW=Ll zPb{+`u~TAM{(ru;hyJr4&^o4Z1|7vyMSX1qXR_) zO1^*Ao{0~q_TrA9*)=L^K3&>BG9^2`lcA%V!j9PcAX?q4sVTYe@$MPm-TjOkBJ>tC z>j&$Scm=oDGCY37k3Gg!V=U2)sxb_{O8($?rjY0W1!?md%<$E1!L<1QRj@f~Gt8&? z>9#%%4=Y_#G1t6G`0mpiJwJrQTsl<2R5km};*n=u@rC*4_}7oV@btfuBiYuPt$!JU zaOfW|EqJ$1{gU6f#Fx~rAnF(OI#d}RdJe_DIU1jExTG}LZ1Xi_H?dCdigviUeZ5L-)1`z;bFX8h_D+n6 z+Rpux-~s_fY6=n840<+&?YTQ>EEd0S(`m|;z|Cd{jbq&Vv)i@Lds&YAEwvH8LmOa} zp4+5y1M@bxk9mv!d{i+*qpsWPoTtGtCk#<_ftaRCp7J7{pMY%vN!WX-j6c^Pp4pwX{;j?C<)wb`^p#4ihh zqOi@!T4Y4_)^wOMYll)hyw&3(c<^8+SPiL=n|O&Et*bK8rS9sl=ACH7v|qo6AFoA$ z-&^Bv>W=hBJXP}5xo~V(J&FP!cGvYHy%sM35Wna~7=WH%)@7~qKpUVQL;wqKdi$@j zcu;-@f6jl&g+ugcqnvdyidLY7>68?|E^b}x&Afu%DOV?JMfbCRINuLEdL2Oj-pAz? zc?IFGs+o!wJE>$79-e*UzCb!NXSP zE0Od)e{0v|7O18{`P2{~2*Myh%Du(BZLGpX;#pm9v|S~S7Z2BOOfOjv8~xaisCKZ0 zpNqoIG#PveeELG0X^GDv*-RHb$owUn zw-3`QYU(Mj+_p1oviiXNON(p|4YYclsz>lU$yKi{GRp_8Ao(bms6NBgmT^BCzTGqj zQfrN_TOv22)qNT4Opy9>l_d&yFmUDE4AFtpACr6XGj#f0w*>8hrX#dSyuFfnv;Ujf zlUjAnEN9=ZFHBzBd>2bt&ruu5@QV%M$)5dG4zO~TG4k-wocx}}0-!^DLZ&^)pep;a zwHIXIgc^4R2dm_|j?NTo=LB*FN&1ce03E2s+|~>7;4)4<4i61CEWJwZ?!yiT)?p5? z#5=HhLEyRG{$U}!MSxQda6klQ)Et`kiz#JZ|2;Uy4rEV<33J37g}+!pj=~Zy;JoVM zqmpKYA~2v9h>*A!5E&VTuEy)RL$k#f^=kR59)ITd?B;iCI|)b8Ig*r; z3)v7k!C_)txH4m6QZ2)L0R%u+1Mt^BP1QjpNlHoH#LlR*>Ji&l6h_sNYr^WQ7MM+yXiApSRWF#`)qi4Vt|I-ekaEg(geZGk{tYwgfL%ScxG&1~@T`bII)mlX^q ze}R*m_O}ZHwxEj(A_ED**=j!I13*i5gdhMEBG^b;Mr`)W&-!gsdFnzD)zlSf=Cyvf zy1VcW{sXyEWc|D&0Pn8m3H>u8Sug~lZxvJM9J6VqdYTV2*np})PV8DJ;V?hPF5m%vCo8#*o#vAuf+niZm^W8OutES~_x6gg2?C)jJGpovp`oUP5Myj5W#f@<|qA zN^!}f)b89fEQiWZ8;toNN`XK>xF0xe?HS3>cY89;-K!{NK?j5nLIzBk2Bc1KvI>ls zM5`F)86`Pgy>4bl!Gf<;pkO|jKjkUceEk~^Uz4HC%e%acgB|j@q|j`_3w|ouLOp>r z#fFT*b|{f}*YU0;vU%14cb_V;{4R%I$D<%-7i? z#XC|`TKwM}tB>ZQ?DllKQQzjc_OV=3+3J(VS(R{ z8R$G}$XaA3nz_yt-8TaN47i|kc5|k5Gm9%OBBeN^m}+3Bv`Ua_D#2l>$YJU&3IeZo zXLt5Lu=*^W;Hx1MdOWml9+0p4YB58c&&9P-S5OLg_bVd7f&QdjETr$0*Yw*GFdOB6 zg?*YKQv%G(U-_u!_`0RKftTKT4(Y2(B|-?)CKw9Jgmv0vGWq0;&V*gwUi1`+3491I z2U@VJB%WYo6`7pi$L>)*xw5Yl&| z$X^Z4tpxb;aEVLrRUq=M+#yTB8hSZ6dp}~f~}hJ-kNlfOuE|B-{xcfgVG{>dc-5obJa;- zDS2#v>)*coOFaM}RK2wlQtTb1mot2E=wc`PZ#$-_i#NaM<3bSp^|-gc?wa+g85$6t zRzm+$yi|s)wB9=hKP{aRW>Q$U#zKSQz+n#|ocMZjLfnOBGd#fa?{^J_miG(Vzer@M70P7n4FDwNY zSCqO^H&^jt)pcuj>UM_e6V_nda(|wdJBhKuEPk6xw%bzTZ-UU*$-0QfG7DYd#d&U( z9ZxRiZrYA(q`~G&a4wC>FQ0Oqk6z0ONmD-=9Uc9PEnz!>01!e6Q%)151PMS6ckLYi zpU@a!B-iD8Ajv>lOPRDf1TNMI&*P^*zWN}d|A_4URO1+{c1r%`m|-6NC^?^mdwZHP zL;cxnsmhbXYO0kfHBaWVt5sfqLZtpZvdMmb z&nnrIu333Bf^z0e_i48szE^%&$HBcF9ir4_5x_zGIp&T7*GRn@^nZ5d zqvvbFlbIq?a84vVWx&*aNflj%50xsbd4V`2@wvH~-?*s5DK)p*>xN_8=})PXxI?q_ z`CAeLFrns!p$$MjSb9ET9tllGi|TaY50MyHAOsumWGzEn+CFHF0SxlOiZvR>7bN5+ zDqblee+;%^Qmd^g_igsdWTlwjjCy3O4`dmT@S+8T_%F_wflv@EX^Ayg*q1s<0Y3Gr zC#$sWd|!*ZAy9SYTt{Ng*V$1z&jD`I^J29MBbRE~pMePPLxYLs3If4L<o1v+GG^q!IP0^5EzbN@^4gb6}Y_`bOrIT2*MkZx1?PH9b$dB=xBEDXT}3e^!sVcHoGDh)4QCk4@V5aW#wfWVK5yeTV~Byx1j zpiCbD81i=I$8bTN5RVAp51SScAoYjSSeGJPP?iQ#)#_ndTJpVw&!1|R)5Zv7@cP0x zv(P=R)w1s?W;DDPiTIdpAQKvKS`Q0RX`x6_zppKiG1Vw4$&p zgaT+~;`Ty-;e$VsO2J389=?|qdu~28{8aB#!A|XUUCAy9Po5hGd3L@dX>tkiPP#YE zT3iIstiBq0lmG`1h1^jDU~HD#p$oGOrxN-iQWI!)pYZ?iDT`tk5BY>{;L1^d4FQMP z-=L*sn{UkHRD=y(Yr#e{p1fl{=e^F)+HU$Oe?`0|-0KK+f8L+X+ZGlb3`nIWC}I~# z2K*sZaPo#7Qh)42fVPd38m@BT*g?t$53*Kr5*@%M{FOMFkNpg;Rr}Z@svul}RFyb7;kq zP3xPcs+;vm5;P)@I$)xuLnrTKI#;#h$3#n^HLI1=0&pXHrU^NBM!TtW{T)UNA^?!a zRI6+)v*0DMsusm4d$4;ygx42B9w7P$5mPVA>h#$XfqGHn&ASJWlJ9(f=4U!)Aye1gx0QuR@Ca6zI z1ZEF-D|PBC-_L!vU94YdC+<;u9!zC-ex(U^kLi@(@qO9Vqw&l?zOET1sCkKB^;Dzr zVttADHr=!2Ww@sAX>Pins+QdBGbDQ~x5BA9NSyh^An;*Dk~G z?E8Hp8|1C~L;_$S^9ATV{Ix;smxTGv-j67TH|T-*&pv#uiYR&>d2b$o{mLR?BWns0 zm5)ONt(HmCzTDjD!|agNZzM(pWI+IMF}VUdHi!r>Ed?r-P3!3Yjf#CIi_$b5D2DM; z<43QHlm3VY3a6Ju$YuSW^cLC+8$uC^JrBiGvWykd|NUV9t#wKH3@N-Hb>v{emgj1a?xH$B!Zgr=5WtF`%{3LX8L{V z90wMq(w+ah`@#~68&EDx1nmm8bOL~trhUwCljrjMq0eX*^0qV3WmR$zJ2Ft^mIaW( z1Yk>pzlM-vzUv}GG(Yfe`kAx8G-r4@dxXJK6~$xa*f)H(- zHm=?)W>uvfHao8lKRWv|-=pcqTV=1PpeXTeeEB}FK@u)fnCe2yf4!U2>m{JReAoMctTvo5n$@G)9g=Nu%Az4fyxJeTMf}|zG5^o^lXkc2 zp1W&3ZWej0-&+HKS4^vJI1DsuCZ(inrB7*mzdSM(m^Hq(t0{^at&aD%-NwV-{cN)h zb`h8V; z+&D{HRI06yo?mG^i$4xO47ojoS7zb{7+Ay9wUHeO9kZlD(zaIh{RsQK$ozuF=h63G zUb3>*y-c&1=s?jr+~0z?F`KK*oAz)Xi7R$i(~`t~)m9i}I9B>=yj8FaTZHZl1v!U{ z=VMG%nGT-V{F(#8smO@^LtXm}HXJtie<@valK_kP_Z!CXdqXEvYdJ#Gc>CQ)IP7XC zV-lSS;cN~a##RnkL9cEv$T@-SRVmbZhlxrO;D}BVK;h!24a*G~n0O58B^MvEi2D!4 z65yx4*W4af5dh`#7HBoW0svIx1v8gpR%+tFA@f+z67KNyU^yiz08SVm0 zY3M;GKxF95K4mvoq~tOy4g58I4a&Bzx{0aWxDDKuhAT1d`JOiR2Ab{Wjsj&zY8C~m z5sj~3`3kZFHgBE97L;!46nfDxQ{yfM6 zfex5cY6UnQu&{;U5I`s!7z}ykQ$?eah8p{g7kAqgXZ)Ki!Ss zUObwzy_@W3Ok_L;Yu}khQ8o#f2q2{Ap??2i!E&eKEM2 zCUTJ@?h1*_`7^o8!nb=WFUNj^;KXX&Y07oS znz&k5F>sWXRM({*KYu;UcU50!M{FqW-qGj2WE*nZY$-`znKNk8*Ig~DNfaIztMm^72x%s5au z-fv#w-;lfZW54Ye6ZokJ^d+OrfB^CnakBWmKd}I%2n@mu7vWQu{$5B*x)GeNO$#Hp z^KCwTu2snN(th?Yr)f{`btOjzzcf2&{CD1OhouFF0>>>gF>B>pxi&f2X=?$;b46C1dra?{2YHsE3lgq>7I~}d$a}ihsv%=sR3(~P zSHkfYqh*hNr3qRQDA^YY-2JvR-Ufa6^~t~DE`MzP7|jJnehow>3wgP#n$Uk#5`|GO&r@FKk}-wwLmmH9`X1chddV zn4keG#QN-zH(#-L20$1?27tIPb8XU>y_auk)yGK_LWw(Ig6p~&msYn)we}X z(AVkcc)R7=RKQCJfaH-o3X^Jgs=$AgFJUpGLlrPwb%XM#2h`{)V|fr7oFC|8;|g$J zAx$rhKS{9Kr!YeUZ{j-M8n{Z_DY!)Z(WOY_?eTrXV%6HjGiarUss>9^9WLP__B%GWftW zwl+MizA821F-TW$m0a5p04pvE;n$itG6@h6%zO6Agzn)AV_r-9?j}Gifucfx%jyyd z*tl`+8N^|<(`*afw!bhFZ95Fk8L3!2mAa7ul%Q>BYp^YdkF>gcg)iuzyBnUi*DN>) zYfA=)t6i-AaT;t$i`3bn_LKsS$)TfGHJT!WO}Nr5{!ca|FV@@cudpTk?l=_RT`Cg| zJf@cT#`Z_`80++G>4ttQ7BrPy9s)0RkFCC8@YuFfzV*t1HT@fU2zpWmrNL_?Fl2UWnWo1noo z&VF45DuHX#<06Rw2Io!CbDrFo;|?7V)r24DWkb_K%KNkEq81={>mohmqK4@)hK|P` z+2rW>zqkE>(I5o^0BD7liX?K2f_@8E)A%LH|?4f)&WDo;UHb$m4472*YhX y)`myvGFdVZFZu0j2=KUy%q%7C7Nje9#hiYwIaZ+}27Ukki@744C`d~X?o|LE2KQ?K literal 0 HcmV?d00001 diff --git a/man/compute_residual_metrics.Rd b/man/compute_residual_metrics.Rd new file mode 100644 index 0000000..c966227 --- /dev/null +++ b/man/compute_residual_metrics.Rd @@ -0,0 +1,87 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/compute_residual_metrics.R +\name{compute_residual_metrics} +\alias{compute_residual_metrics} +\title{Residual-based posterior predictive metrics for longitudinal antibody curves} +\usage{ +compute_residual_metrics( + model, + dataset, + ids, + antigen_iso, + scale = c("original", "log"), + summary_level = c("id_antigen", "pointwise", "antigen", "overall") +) +} +\arguments{ +\item{model}{A data frame of posterior draws in long format with columns: +\code{Subject}, \code{Iso_type}, \code{Chain}, \code{Iteration}, +\code{Parameter}, \code{value}.} + +\item{dataset}{A \code{serodynamics} case dataset produced by +\code{serodynamics::as_case_data()} (must contain \code{id} and +\code{antigen_iso} columns, and time/value attributes).} + +\item{ids}{Character vector of subject IDs to include (matched against +\code{dataset$id} and \code{model$Subject}).} + +\item{antigen_iso}{Character scalar specifying the antigen/isotype to analyze +(matched against \code{dataset$antigen_iso} and \code{model$Iso_type}).} + +\item{scale}{Scale on which to compute residuals. One of \code{"original"} +or \code{"log"}. If \code{"log"}, residuals are computed on the natural +log scale and observations/predictions \eqn{\le 0} are removed.} + +\item{summary_level}{Level at which to summarize metrics. One of: +\describe{ +\item{\code{"pointwise"}}{ +Return pointwise residuals at each observed time point. +} +\item{\code{"id_antigen"}}{ +Summarize by \code{id} and \code{antigen_iso}. +} +\item{\code{"antigen"}}{Summarize by \code{antigen_iso} only.} +\item{\code{"overall"}}{Single summary across all included observations.} +}} +} +\value{ +A tibble. For \code{summary_level = "pointwise"}, returns +per-observation residuals. Otherwise returns MAE, RMSE, SSE, and +\code{n_obs} at the requested summary level. +} +\description{ +Computes residuals between observed antibody measurements and posterior +median predictions evaluated at the observed time points. Returns pointwise +residuals or aggregated error metrics (MAE, RMSE, SSE) at multiple +summary levels. +} +\details{ +The posterior predictive summary uses the posterior median at each observed +time point, with a 95\\% credible interval computed from draw-level +predictions. Predictions are generated by +\code{\link{predict_posterior_at_times}}. +} +\examples{ +\dontrun{ +# Per-ID error metrics (original scale) +m_id <- compute_residual_metrics( + model = overall_sf2a, + dataset = dL_clean_sf2a, + ids = unique(dL_clean_sf2a$id), + antigen_iso = "IgG", + scale = "original", + summary_level = "id_antigen" +) + +# Pointwise residuals on log scale +r_pw <- compute_residual_metrics( + model = overall_sf2a, + dataset = dL_clean_sf2a, + ids = "SOSAR-22008", + antigen_iso = "IgA", + scale = "log", + summary_level = "pointwise" +) +} + +} diff --git a/man/fig2_overall_newperson.Rd b/man/fig2_overall_newperson.Rd new file mode 100644 index 0000000..c8a5701 --- /dev/null +++ b/man/fig2_overall_newperson.Rd @@ -0,0 +1,69 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/fig2_overall_newperson.R +\name{fig2_overall_newperson} +\alias{fig2_overall_newperson} +\title{Summarize population-level ("newperson") antibody trajectories +from overall models} +\usage{ +fig2_overall_newperson( + overall_models, + osps = c("IpaB", "Sf2a", "Sf3a", "Sf6", "Sonnei"), + ids = "newperson", + isotypes = c("IgG", "IgA"), + t_grid = seq(0, 210, by = 5), + cred = 0.95, + log_y = TRUE, + xlim = c(0, 210), + ylab = "Normalized MFI", + line_color = "#1f77b4", + ribbon_alpha = 0.2, + facet_scales = "fixed", + return_data = FALSE +) +} +\arguments{ +\item{overall_models}{Named list of posterior draws in long format +(one per antigen), +with columns at least: Subject, Iso_type, Chain, Iteration, +Parameter, value.} + +\item{osps}{Character vector of antigen names +(must match names in \code{overall_models}).} + +\item{ids}{Character vector of Subject IDs to include (default: "newperson").} + +\item{isotypes}{Character vector of isotypes +(default: c("IgG", "IgA")).} + +\item{t_grid}{Numeric vector of time points (days) to evaluate.} + +\item{cred}{Credible level (default 0.95).} + +\item{log_y}{Logical; if TRUE, applies log10 scale to y when returning plot.} + +\item{xlim}{Optional numeric length-2 vector for x-axis limits when +returning plot.} + +\item{ylab}{Y-axis label for plot.} + +\item{line_color}{Line color for plot.} + +\item{ribbon_alpha}{Alpha for credible ribbon in plot.} + +\item{facet_scales}{Passed to ggplot facet scales.} + +\item{return_data}{If TRUE, returns a list with \code{plot} and \code{data}.} +} +\value{ +By default, a ggplot. If \code{return_data = TRUE}, returns +list(plot = p, data = df). +} +\description{ +Computes median and credible interval bands of the antibody trajectory for a +hypothetical new individual drawn from the population distribution +("newperson"). +} +\details{ +Requires that each model draw can be pivoted to wide parameters including: +y0, y1, t1, alpha, shape. +} diff --git a/man/fmt_mci.Rd b/man/fmt_mci.Rd new file mode 100644 index 0000000..13129eb --- /dev/null +++ b/man/fmt_mci.Rd @@ -0,0 +1,25 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/fmt_mci.R +\name{fmt_mci} +\alias{fmt_mci} +\title{Format median and credible interval as a single string} +\usage{ +fmt_mci(med, lo, hi, digits = 2, sci = FALSE) +} +\arguments{ +\item{med}{Median value.} + +\item{lo}{Lower bound.} + +\item{hi}{Upper bound.} + +\item{digits}{Number of digits.} + +\item{sci}{Logical; if TRUE uses scientific notation.} +} +\value{ +A string like "1.23 (0.50ā€“2.00)". +} +\description{ +Format median and credible interval as a single string +} diff --git a/man/make_model_comparison_table.Rd b/man/make_model_comparison_table.Rd new file mode 100644 index 0000000..783e9dd --- /dev/null +++ b/man/make_model_comparison_table.Rd @@ -0,0 +1,40 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/model_comparison_table.R +\name{make_model_comparison_table} +\alias{make_model_comparison_table} +\title{Compare serotype-specific vs overall models using residual metrics} +\usage{ +make_model_comparison_table( + model_serospecific, + data_serospecific, + model_overall, + data_overall, + antigen_iso, + scale = c("original", "log"), + tie_tol = 1e-08 +) +} +\arguments{ +\item{model_serospecific}{Posterior draws (long format) +for the serotype-specific model.} + +\item{data_serospecific}{Case dataset used for the serotype-specific model.} + +\item{model_overall}{Posterior draws (long format) for the overall model.} + +\item{data_overall}{Case dataset used for the overall model.} + +\item{antigen_iso}{Character scalar antigen/isotype label.} + +\item{scale}{"original" or "log".} + +\item{tie_tol}{Numeric tolerance to declare a tie.} +} +\value{ +A tibble with per-ID MAE/RMSE for each model and +deltas/winner labels. +} +\description{ +Computes per-ID residual metrics for two models on the intersection of IDs +present in both datasets, then reports absolute and percent differences. +} diff --git a/man/mock_case_data.Rd b/man/mock_case_data.Rd new file mode 100644 index 0000000..22c3509 --- /dev/null +++ b/man/mock_case_data.Rd @@ -0,0 +1,39 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/data.R +\docType{data} +\name{mock_case_data} +\alias{mock_case_data} +\title{Mock case data for testing} +\format{ +A data frame with class \code{c("case_data", "data.frame")} +and columns: +\describe{ +\item{id}{Character. Subject ID} +\item{antigen_iso}{Character. Isotype (e.g., "IgG", "IgA")} +\item{timepoint}{Numeric. Time in days since infection} +\item{value}{Numeric. Antibody measurement} +} +} +\usage{ +mock_case_data +} +\description{ +A mock longitudinal antibody dataset compatible with +\code{serodynamics::as_case_data()}, for testing functions like +\code{\link{compute_residual_metrics}}. +} +\details{ +This is synthetic data generated for testing and examples. Real Shigella +case data will be added separately. + +The dataset has attributes: +\itemize{ +\item \code{attr(mock_case_data, "timeindays") = "timepoint"} +\item \code{attr(mock_case_data, "value_var") = "value"} +} +} +\examples{ +head(mock_case_data) +table(mock_case_data$id) +} +\keyword{datasets} diff --git a/man/mock_posterior_draws.Rd b/man/mock_posterior_draws.Rd new file mode 100644 index 0000000..136fe75 --- /dev/null +++ b/man/mock_posterior_draws.Rd @@ -0,0 +1,43 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/data.R +\docType{data} +\name{mock_posterior_draws} +\alias{mock_posterior_draws} +\title{Mock posterior draws for testing} +\format{ +A data frame with columns: +\describe{ +\item{Subject}{Character. Subject ID (e.g., "newperson", "SOSAR-22008")} +\item{Iso_type}{Character. Isotype (e.g., "IgG", "IgA")} +\item{Chain}{Integer. MCMC chain number} +\item{Iteration}{Integer. MCMC iteration number} +\item{Parameter}{Character. Parameter name (y0, y1, t1, alpha, shape)} +\item{value}{Numeric. Parameter value} +} +} +\usage{ +mock_posterior_draws +} +\description{ +A mock dataset of posterior parameter draws in long format, mimicking the +structure expected by functions like \code{\link{compute_residual_metrics}} +and \code{\link{predict_posterior_at_times}}. +} +\details{ +This is synthetic data generated for testing and examples. Real Shigella +posterior draws will be added to the package separately. + +Parameters represent: +\itemize{ +\item \code{y0}: Baseline antibody level +\item \code{y1}: Peak antibody level +\item \code{t1}: Time to peak (days) +\item \code{alpha}: Decay rate parameter +\item \code{shape}: Decay shape parameter (rho) +} +} +\examples{ +head(mock_posterior_draws) +table(mock_posterior_draws$Parameter) +} +\keyword{datasets} diff --git a/man/model_comparison_table.Rd b/man/model_comparison_table.Rd new file mode 100644 index 0000000..7fda11a --- /dev/null +++ b/man/model_comparison_table.Rd @@ -0,0 +1,35 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/model_comparison_table.R +\name{model_comparison_table} +\alias{model_comparison_table} +\title{Compare summary residual metrics between two model outputs} +\usage{ +model_comparison_table( + metrics_overall, + metrics_pointwise, + model_overall_label = "Overall Model", + model_pointwise_label = "Pointwise Model" +) +} +\arguments{ +\item{metrics_overall}{Data frame with one row containing at least +\code{MAE}, \code{RMSE}, \code{SSE}, and \code{n_obs} +for the overall model.} + +\item{metrics_pointwise}{Data frame with one row containing at least +\code{MAE}, \code{RMSE}, \code{SSE}, and \code{n_obs} +for the pointwise model.} + +\item{model_overall_label}{Label used for the overall model row.} + +\item{model_pointwise_label}{Label used for the pointwise model row.} +} +\value{ +A tibble with one row per model and comparison columns: +\code{delta_MAE}, \code{delta_RMSE}, \code{pct_improve_MAE}, +\code{pct_improve_RMSE}. +} +\description{ +Builds a compact comparison table from pre-computed summary metrics for an +overall model and a pointwise model. +} diff --git a/man/postprocess_jags_output.Rd b/man/postprocess_jags_output.Rd deleted file mode 100644 index 001ab79..0000000 --- a/man/postprocess_jags_output.Rd +++ /dev/null @@ -1,17 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/postprocess_jags_output.R -\name{postprocess_jags_output} -\alias{postprocess_jags_output} -\title{Postprocess JAGS output} -\usage{ -postprocess_jags_output(jags_output) -} -\arguments{ -\item{jags_output}{output from \code{\link[runjags:run.jags]{runjags::run.jags()}}} -} -\value{ -a \link[tibble:tbl_df-class]{tibble::tbl_df} -} -\description{ -Postprocess JAGS output -} diff --git a/man/predict_posterior_at_times.Rd b/man/predict_posterior_at_times.Rd new file mode 100644 index 0000000..c18258d --- /dev/null +++ b/man/predict_posterior_at_times.Rd @@ -0,0 +1,64 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/predict_posterior_at_times.R +\name{predict_posterior_at_times} +\alias{predict_posterior_at_times} +\title{Posterior predictions at specified times for given subjects +and antigen/isotype} +\usage{ +predict_posterior_at_times(model, ids, antigen_iso, times) +} +\arguments{ +\item{model}{A data frame of posterior draws in long format with columns: +\code{Subject}, \code{Iso_type}, \code{Chain}, \code{Iteration}, +\code{Parameter}, \code{value}.} + +\item{ids}{Character vector of subject IDs to include +(matched against \code{Subject}).} + +\item{antigen_iso}{Character scalar specifying the antigen/isotype to include +(matched against \code{Iso_type}).} + +\item{times}{Numeric vector of time points (days) at which to evaluate +predictions.} +} +\value{ +A tibble with one row per +(posterior draw \eqn{\times} time \eqn{\times} subject), including the +evaluated prediction \code{res}. Output includes at least: +\describe{ +\item{id}{Subject ID (character).} +\item{t}{Time (days) at which prediction was evaluated.} +\item{Chain}{MCMC chain index (if present in \code{model}).} +\item{Iteration}{MCMC iteration index (if present in \code{model}).} +\item{sample_id}{Row index for the draw (added if missing).} +\item{y0, y1, t1, alpha, shape}{Model parameters (wide).} +\item{res}{Predicted antibody level at time \code{t}.} +} +} +\description{ +Generates draw-level posterior predictions of the antibody trajectory at +user-specified time points, for one or more subjects and a selected +antigen/isotype. This is a low-level helper used by +residual-based posterior predictive diagnostics. +} +\details{ +This function pivots posterior draws to wide format (parameters as columns), +expands them over \code{times}, and evaluates the antibody curve via +an internal implementation of the antibody kinetics model using parameters +\code{y0}, \code{y1}, \code{t1}, \code{alpha}, and \code{shape}. +} +\examples{ +\dontrun{ +preds <- predict_posterior_at_times( + model = overall_sf2a, + ids = "newperson", + antigen_iso = "IgG", + times = c(0, 30, 90, 180) +) +} + +} +\seealso{ +\code{\link{compute_residual_metrics}} +} +\keyword{internal} diff --git a/man/prep_newperson_params.Rd b/man/prep_newperson_params.Rd new file mode 100644 index 0000000..83811be --- /dev/null +++ b/man/prep_newperson_params.Rd @@ -0,0 +1,21 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/prep_newperson_params.R +\name{prep_newperson_params} +\alias{prep_newperson_params} +\title{Extract "newperson" parameter draws and summarize for Table 2} +\usage{ +prep_newperson_params(draws_long, antigen_label) +} +\arguments{ +\item{draws_long}{Posterior draws in long format with columns: +Subject, Iso_type, Chain, Iteration, Parameter, value.} + +\item{antigen_label}{Character scalar antigen name to attach.} +} +\value{ +A tibble of newperson draws in wide format with columns: +Iteration, Chain, antigen, Iso_type, y0, y1, t1, alpha, rho. +} +\description{ +Extract "newperson" parameter draws and summarize for Table 2 +} diff --git a/man/process_shigella_data.Rd b/man/process_shigella_data.Rd new file mode 100644 index 0000000..491c184 --- /dev/null +++ b/man/process_shigella_data.Rd @@ -0,0 +1,42 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/process_shigella_data.R +\name{process_shigella_data} +\alias{process_shigella_data} +\title{Reshape Shigella longitudinal data for serodynamics workflows} +\usage{ +process_shigella_data(data, study_filter, antigen) +} +\arguments{ +\item{data}{A data frame containing longitudinal measurements.} + +\item{study_filter}{Character scalar. Value of \code{study_name} to keep +(e.g. "SOSAR").} + +\item{antigen}{Unquoted column name for the antigen measurement +(e.g. n_ipab_MFI).} +} +\value{ +A tibble with standardized columns: +\describe{ +\item{index_id}{Participant ID (copied from \code{sid}).} +\item{antigen_iso}{Isotype label (copied from \code{isotype_name}).} +\item{visit}{Visit label (copied from \code{timepoint}).} +\item{timeindays}{Time since infection (copied from \verb{Actual day}).} +\item{result}{Antibody measurement (from \code{antigen}).} +} +} +\description{ +Filters a dataset to a given study and antigen column, standardizes column +names, and returns a visit-ordered long dataset suitable for conversion to +\code{serodynamics::as_case_data()}. +} +\examples{ +\dontrun{ +dat_long <- process_shigella_data(df, "SOSAR", n_ipab_MFI) +dL <- serodynamics::as_case_data(dat_long, + id_var = "index_id", biomarker_var = "antigen_iso", + time_in_days = "timeindays", value_var = "result" +) +} + +} diff --git a/man/utils_internal.Rd b/man/utils_internal.Rd new file mode 100644 index 0000000..0baed5f --- /dev/null +++ b/man/utils_internal.Rd @@ -0,0 +1,9 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/utils_internal.R +\name{utils_internal} +\alias{utils_internal} +\title{Internal utility functions} +\description{ +Internal utility functions +} +\keyword{internal} diff --git a/tests/testthat.R b/tests/testthat.R new file mode 100644 index 0000000..8f6c457 --- /dev/null +++ b/tests/testthat.R @@ -0,0 +1,12 @@ +# This file is part of the standard setup for testthat. +# It is recommended that you do not modify it. +# +# Where should you do additional test configuration? +# Learn more about the roles of various files in: +# * https://r-pkgs.org/testing-design.html#sec-tests-files-overview +# * https://testthat.r-lib.org/articles/special-files.html + +library(testthat) +library(shigella) + +test_check("shigella") diff --git a/tests/testthat/test-compute_residual_metrics.R b/tests/testthat/test-compute_residual_metrics.R new file mode 100644 index 0000000..1d07832 --- /dev/null +++ b/tests/testthat/test-compute_residual_metrics.R @@ -0,0 +1,102 @@ +test_that("compute_residual_metrics returns expected structure", { + # Note: This test uses mock data structure + # Once mock_posterior_draws and mock_case_data are generated, + # these tests should work + + # Create minimal mock data inline for testing + mock_model <- data.frame( + Subject = rep("test_id", 15), + Iso_type = rep("IgG", 15), + Chain = rep(1, 15), + Iteration = rep(1:3, 5), + Parameter = rep(c("y0", "y1", "t1", "alpha", "shape"), each = 3), + value = c( + 200, 210, 190, # y0 + 2000, 2100, 1900, # y1 + 10, 12, 9, # t1 + 0.02, 0.025, 0.015, # alpha + 1.0, 1.1, 0.9 # shape + ) + ) + + mock_dataset <- data.frame( + id = rep("test_id", 4), + antigen_iso = rep("IgG", 4), + timepoint = c(0, 7, 30, 90), + value = c(250, 1500, 1200, 800) + ) + attr(mock_dataset, "timeindays") <- "timepoint" + attr(mock_dataset, "value_var") <- "value" + class(mock_dataset) <- c("case_data", "data.frame") + + # Test pointwise summary + result_pw <- compute_residual_metrics( + model = mock_model, + dataset = mock_dataset, + ids = "test_id", + antigen_iso = "IgG", + scale = "original", + summary_level = "pointwise" + ) + + expect_s3_class(result_pw, "data.frame") + expect_true(nrow(result_pw) > 0) + expect_true(all(c("id", "antigen_iso", "t", "obs", "pred_med", + "residual", "abs_residual", "sq_residual") %in% names(result_pw))) + expect_true(all(is.finite(result_pw$residual))) + + # Test id_antigen summary + result_id <- compute_residual_metrics( + model = mock_model, + dataset = mock_dataset, + ids = "test_id", + antigen_iso = "IgG", + scale = "original", + summary_level = "id_antigen" + ) + + expect_s3_class(result_id, "data.frame") + expect_equal(nrow(result_id), 1) + expect_true(all(c("id", "antigen_iso", "MAE", "RMSE", "SSE", "n_obs") %in% names(result_id))) + expect_true(result_id$MAE >= 0) + expect_true(result_id$RMSE >= 0) + expect_true(result_id$SSE >= 0) + expect_true(result_id$n_obs > 0) +}) + +test_that("compute_residual_metrics handles log scale", { + # Minimal mock data + mock_model <- data.frame( + Subject = rep("test_id", 15), + Iso_type = rep("IgG", 15), + Chain = rep(1, 15), + Iteration = rep(1:3, 5), + Parameter = rep(c("y0", "y1", "t1", "alpha", "shape"), each = 3), + value = c(200, 210, 190, 2000, 2100, 1900, 10, 12, 9, + 0.02, 0.025, 0.015, 1.0, 1.1, 0.9) + ) + + mock_dataset <- data.frame( + id = rep("test_id", 3), + antigen_iso = rep("IgG", 3), + timepoint = c(7, 30, 90), + value = c(1500, 1200, 800) + ) + attr(mock_dataset, "timeindays") <- "timepoint" + attr(mock_dataset, "value_var") <- "value" + class(mock_dataset) <- c("case_data", "data.frame") + + result_log <- compute_residual_metrics( + model = mock_model, + dataset = mock_dataset, + ids = "test_id", + antigen_iso = "IgG", + scale = "log", + summary_level = "overall" + ) + + expect_s3_class(result_log, "data.frame") + expect_equal(nrow(result_log), 1) + expect_true(all(is.finite(result_log$MAE))) + expect_true(all(is.finite(result_log$RMSE))) +}) diff --git a/tests/testthat/test-model_comparison_table.R b/tests/testthat/test-model_comparison_table.R new file mode 100644 index 0000000..4f8f81f --- /dev/null +++ b/tests/testthat/test-model_comparison_table.R @@ -0,0 +1,68 @@ +test_that("model_comparison_table computes differences correctly", { + # Create mock metrics for two models + metrics_overall <- data.frame( + MAE = 150, + RMSE = 200, + SSE = 40000, + n_obs = 100 + ) + + metrics_pointwise <- data.frame( + MAE = 120, + RMSE = 180, + SSE = 32400, + n_obs = 100 + ) + + result <- model_comparison_table( + metrics_overall = metrics_overall, + metrics_pointwise = metrics_pointwise, + model_overall_label = "Overall Model", + model_pointwise_label = "Pointwise Model" + ) + + expect_s3_class(result, "data.frame") + expect_equal(nrow(result), 2) + expect_true(all(c("Model", "MAE", "RMSE", "SSE", "n_obs") %in% names(result))) + expect_true("delta_MAE" %in% names(result)) + expect_true("delta_RMSE" %in% names(result)) + + # Check that differences are computed correctly + delta_mae <- result$delta_MAE[result$Model == "Pointwise Model"] + expect_equal(delta_mae, 120 - 150) + + delta_rmse <- result$delta_RMSE[result$Model == "Pointwise Model"] + expect_equal(delta_rmse, 180 - 200) +}) + +test_that("model_comparison_table computes percent improvement", { + metrics_overall <- data.frame( + MAE = 100, + RMSE = 200, + SSE = 40000, + n_obs = 100 + ) + + metrics_pointwise <- data.frame( + MAE = 80, + RMSE = 160, + SSE = 25600, + n_obs = 100 + ) + + result <- model_comparison_table( + metrics_overall = metrics_overall, + metrics_pointwise = metrics_pointwise + ) + + expect_true("pct_improve_MAE" %in% names(result)) + expect_true("pct_improve_RMSE" %in% names(result)) + + # Check percent improvements are computed correctly + # Pointwise model: (80 - 100) / 100 * 100 = -20% (20% improvement) + pct_mae <- result$pct_improve_MAE[result$Model == "Pointwise Model"] + expect_equal(pct_mae, -20) + + pct_rmse <- result$pct_improve_RMSE[result$Model == "Pointwise Model"] + expect_equal(pct_rmse, -20) +}) diff --git a/tests/testthat/test-predict_posterior_at_times.R b/tests/testthat/test-predict_posterior_at_times.R new file mode 100644 index 0000000..0cf4b08 --- /dev/null +++ b/tests/testthat/test-predict_posterior_at_times.R @@ -0,0 +1,68 @@ +test_that("predict_posterior_at_times returns expected structure", { + # Create minimal mock posterior data + mock_model <- data.frame( + Subject = rep("test_id", 15), + Iso_type = rep("IgG", 15), + Chain = rep(1, 15), + Iteration = rep(1:3, 5), + Parameter = rep(c("y0", "y1", "t1", "alpha", "shape"), each = 3), + value = c( + 200, 210, 190, # y0 + 2000, 2100, 1900, # y1 + 10, 12, 9, # t1 + 0.02, 0.025, 0.015, # alpha + 1.0, 1.1, 0.9 # shape + ) + ) + + times <- c(0, 30, 90) + + result <- predict_posterior_at_times( + model = mock_model, + ids = "test_id", + antigen_iso = "IgG", + times = times + ) + + expect_s3_class(result, "data.frame") + expect_true(nrow(result) > 0) + expect_true(all(c("id", "t", "res") %in% names(result))) + expect_true(all(result$t %in% times)) + expect_true(all(is.finite(result$res))) + expect_true(all(result$res > 0)) +}) + +test_that("predict_posterior_at_times handles multiple subjects", { + # Create mock data for multiple subjects + mock_model <- rbind( + data.frame( + Subject = rep("id1", 15), + Iso_type = rep("IgG", 15), + Chain = rep(1, 15), + Iteration = rep(1:3, 5), + Parameter = rep(c("y0", "y1", "t1", "alpha", "shape"), each = 3), + value = c(200, 210, 190, 2000, 2100, 1900, 10, 12, 9, + 0.02, 0.025, 0.015, 1.0, 1.1, 0.9) + ), + data.frame( + Subject = rep("id2", 15), + Iso_type = rep("IgG", 15), + Chain = rep(1, 15), + Iteration = rep(1:3, 5), + Parameter = rep(c("y0", "y1", "t1", "alpha", "shape"), each = 3), + value = c(250, 260, 240, 2500, 2600, 2400, 12, 14, 11, + 0.03, 0.035, 0.025, 1.2, 1.3, 1.1) + ) + ) + + result <- predict_posterior_at_times( + model = mock_model, + ids = c("id1", "id2"), + antigen_iso = "IgG", + times = c(0, 30) + ) + + expect_s3_class(result, "data.frame") + expect_true(all(c("id1", "id2") %in% result$id)) + expect_true(all(is.finite(result$res))) +}) diff --git a/tests/testthat/test-process_shigella_data.R b/tests/testthat/test-process_shigella_data.R new file mode 100644 index 0000000..cf7323e --- /dev/null +++ b/tests/testthat/test-process_shigella_data.R @@ -0,0 +1,46 @@ +test_that("process_shigella_data filters and reshapes correctly", { + # Create mock input data + mock_raw_data <- data.frame( + study_name = rep(c("SOSAR", "OTHER"), each = 6), + isotype_name = rep(c("IgG", "IgA"), 6), + sid = rep(c("ID001", "ID002"), each = 6), + timepoint = rep(c("V1", "V2", "V3"), 4), + `Actual day` = rep(c(0, 7, 30), 4), + n_ipab_MFI = runif(12, 100, 5000), + check.names = FALSE + ) + + result <- process_shigella_data( + data = mock_raw_data, + study_filter = "SOSAR", + antigen = n_ipab_MFI + ) + + expect_s3_class(result, "data.frame") + expect_true(nrow(result) > 0) + expect_true(all(c("index_id", "antigen_iso", "visit", "timeindays", "result") %in% names(result))) + expect_true(all(result$index_id %in% c("ID001", "ID002"))) + expect_true(all(result$antigen_iso %in% c("IgG", "IgA"))) + expect_true(all(!is.na(result$timeindays))) +}) + +test_that("process_shigella_data removes NA timepoints", { + mock_raw_data <- data.frame( + study_name = rep("SOSAR", 6), + isotype_name = rep("IgG", 6), + sid = rep("ID001", 6), + timepoint = c("V1", "V2", "V3", "V4", "V5", "V6"), + `Actual day` = c(0, 7, NA, 30, NA, 90), + n_ipab_MFI = runif(6, 100, 5000), + check.names = FALSE + ) + + result <- process_shigella_data( + data = mock_raw_data, + study_filter = "SOSAR", + antigen = n_ipab_MFI + ) + + expect_true(all(!is.na(result$timeindays))) + expect_equal(nrow(result), 4) # Should only have 4 rows (excluding 2 NAs) +}) diff --git a/tests/testthat/test-utils_internal.R b/tests/testthat/test-utils_internal.R new file mode 100644 index 0000000..5cd71f6 --- /dev/null +++ b/tests/testthat/test-utils_internal.R @@ -0,0 +1,45 @@ +test_that("ab() returns y0 at t = 0", { + result <- ab(t = 0, y0 = 200, y1 = 2000, t1 = 10, alpha = 0.02, shape = 1.2) + expect_equal(result, 200) +}) + +test_that("ab() returns y1 at t = t1", { + result <- ab(t = 10, y0 = 200, y1 = 2000, t1 = 10, alpha = 0.02, shape = 1.2) + expect_equal(result, 2000) +}) + +test_that("ab() shows expected rise for t < t1", { + result <- ab(t = 5, y0 = 200, y1 = 2000, t1 = 10, alpha = 0.02, shape = 1.2) + expected <- 200 * exp((log(2000 / 200) / 10) * 5) + expect_equal(result, expected) +}) + +test_that("ab() decays toward zero for large t (shape > 1)", { + result <- ab(t = 10000, y0 = 200, y1 = 2000, t1 = 10, alpha = 0.02, shape = 1.2) + expect_true(is.finite(result)) + expect_true(result >= 0) +}) + +test_that("ab() handles vector input", { + times <- c(0, 5, 10, 30, 90) + result <- ab(t = times, y0 = 200, y1 = 2000, t1 = 10, alpha = 0.02, shape = 1.2) + expect_length(result, 5) + expect_true(all(is.finite(result))) + expect_equal(result[1], 200) + expect_equal(result[3], 2000) +}) + +test_that("ab() does not error when t1 = 0", { + expect_no_error(ab(t = c(0, 10, 30), y0 = 200, y1 = 2000, t1 = 0, alpha = 0.02, shape = 1.2)) +}) + +test_that("get_timeindays_var() finds time variable from attribute", { + mock_data <- data.frame(timepoint = c(0, 7, 30), value = c(100, 200, 300)) + attr(mock_data, "timeindays") <- "timepoint" + expect_equal(get_timeindays_var(mock_data), "timepoint") +}) + +test_that("get_timeindays_var() defaults to timeindays when attribute missing", { + mock_data <- data.frame(timeindays = c(0, 7, 30), value = c(100, 200, 300)) + expect_equal(get_timeindays_var(mock_data), "timeindays") +}) diff --git a/vignettes/getting-started.Rmd b/vignettes/getting-started.Rmd new file mode 100644 index 0000000..d25dd23 --- /dev/null +++ b/vignettes/getting-started.Rmd @@ -0,0 +1,134 @@ +--- +title: "Getting Started with shigella" +output: rmarkdown::html_vignette +vignette: > + %\VignetteIndexEntry{Getting Started with shigella} + %\VignetteEngine{knitr::rmarkdown} + %\VignetteEncoding{UTF-8} +--- + +```{r, include = FALSE} +knitr::opts_chunk$set( + collapse = TRUE, + comment = "#>" +) +``` + +```{r setup} +library(shigella) +``` + +## Overview + +The `shigella` package provides tools for analyzing longitudinal antibody kinetics data from Shigella infection studies. This vignette introduces the main functionality using mock data. + +**Note:** This package currently uses mock data for testing and examples. Real Shigella datasets will be added in a future release. + +## Key Functions + +### Data Processing + +The `process_shigella_data()` function reshapes raw longitudinal data into a format compatible with the `serodynamics` package: + +```{r, eval=FALSE} +# Example with real data (not run) +dat_long <- process_shigella_data( + data = raw_data, + study_filter = "SOSAR", + antigen = n_ipab_MFI +) + +# Convert to case_data format +dL <- serodynamics::as_case_data( + dat_long, + id_var = "index_id", + biomarker_var = "antigen_iso", + time_in_days = "timeindays", + value_var = "result" +) +``` + +### Model Evaluation + +The package provides functions to evaluate and compare longitudinal antibody models: + +#### Computing Residual Metrics + +```{r, eval=FALSE} +# Compute residuals at individual level +metrics_id <- compute_residual_metrics( + model = posterior_draws, + dataset = case_data, + ids = unique(case_data$id), + antigen_iso = "IgG", + scale = "original", + summary_level = "id_antigen" +) + +# Overall summary metrics +metrics_overall <- compute_residual_metrics( + model = posterior_draws, + dataset = case_data, + ids = unique(case_data$id), + antigen_iso = "IgG", + scale = "original", + summary_level = "overall" +) +``` + +#### Model Comparison + +Compare model fit metrics between different modeling approaches: + +```{r, eval=FALSE} +comparison <- model_comparison_table( + metrics_overall = metrics_model1, + metrics_pointwise = metrics_model2, + model_overall_label = "Population Model", + model_pointwise_label = "Individual Model" +) +``` + +### Visualization + +The `fig2_overall_newperson()` function creates trajectory plots for population-level ("newperson") predictions: + +```{r, eval=FALSE} +fig <- fig2_overall_newperson( + overall_models = list( + IpaB = model_ipab, + Sf2a = model_sf2a + ), + osps = c("IpaB", "Sf2a"), + isotypes = c("IgG", "IgA"), + t_grid = seq(0, 210, by = 5), + log_y = TRUE +) +``` + +## Mock Data + +The package includes two mock datasets for testing: + +- `mock_posterior_draws`: Posterior parameter draws in long format +- `mock_case_data`: Longitudinal antibody measurements + +These datasets have the same structure as real data but contain synthetic values. + +```{r, eval=FALSE} +# Load mock data +data("mock_posterior_draws") +data("mock_case_data") + +# View structure +head(mock_posterior_draws) +head(mock_case_data) +``` + +## Next Steps + +For more detailed examples and analysis workflows, see: + +- Package documentation: `?shigella` +- Function help: `?compute_residual_metrics`, `?process_shigella_data` +- GitHub repository: https://github.com/UCD-SERG/shigella From 756544dba6196d8aee16ee518b15797ca6570c56 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Tue, 31 Mar 2026 18:08:15 +0000 Subject: [PATCH 23/30] Fix 5 CI failures: version bump, LICENSE, non-ASCII chars, VignetteBuilder, unused Imports, WORDLIST, line lengths Agent-Logs-Url: https://github.com/UCD-SERG/shigella/sessions/be474f7c-718f-491f-84a7-bf1b17ecce18 Co-authored-by: Kwan-Jenny <68584166+Kwan-Jenny@users.noreply.github.com> --- DESCRIPTION | 12 +++++++----- LICENSE | 2 ++ R/fmt_mci.R | 10 +++++++--- R/prep_newperson_params.R | 8 ++++++-- inst/WORDLIST | 16 ++++++++++++++++ 5 files changed, 38 insertions(+), 10 deletions(-) create mode 100644 LICENSE diff --git a/DESCRIPTION b/DESCRIPTION index 50e8fa8..52f29cd 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ Package: shigella Title: What the Package Does (One Line, Title Case) -Version: 0.0.0.9004 +Version: 0.0.0.9005 Authors@R: c( person("Kwan Ho", "Lee", , "ksjlee@ucdavis.edu", role = c("aut", "cre")), person("Douglas Ezra", "Morrison", , "demorrison@ucdavis.edu", role = c("aut"), @@ -13,13 +13,9 @@ RoxygenNote: 7.3.3 Imports: cli, dplyr, - future, - future.apply, ggplot2, - gridExtra, rlang, serocalculator, - serodynamics (>= 0.0.0.9011), tibble, tidyr URL: https://ucd-serg.github.io/shigella/ @@ -31,8 +27,11 @@ Suggests: coda, forcats, furrr, + future, + future.apply, ggeasy, ggmcmc, + gridExtra, gtsummary, haven, here, @@ -47,12 +46,15 @@ Suggests: purrr, quarto, readxl, + rmarkdown, runjags, scales, + serodynamics (>= 0.0.0.9011), spelling, table1, testthat (>= 3.0.0), tidyverse +VignetteBuilder: knitr Depends: R (>= 3.5) LazyData: true diff --git a/LICENSE b/LICENSE new file mode 100644 index 0000000..aec8417 --- /dev/null +++ b/LICENSE @@ -0,0 +1,2 @@ +YEAR: 2024 +COPYRIGHT HOLDER: Kwan Ho Lee, Douglas Ezra Morrison diff --git a/R/fmt_mci.R b/R/fmt_mci.R index 9982358..3cc4162 100644 --- a/R/fmt_mci.R +++ b/R/fmt_mci.R @@ -6,12 +6,16 @@ #' @param digits Number of digits. #' @param sci Logical; if TRUE uses scientific notation. #' -#' @return A string like "1.23 (0.50ā€“2.00)". +#' @return A string like "1.23 (0.50--2.00)". #' #' @export fmt_mci <- function(med, lo, hi, digits = 2, sci = FALSE) { f <- function(x) { - if (sci) formatC(x, format = "e", digits = digits) else formatC(x, format = "f", digits = digits) + if (sci) { + formatC(x, format = "e", digits = digits) + } else { + formatC(x, format = "f", digits = digits) + } } - sprintf("%s (%sā€“%s)", f(med), f(lo), f(hi)) + sprintf("%s (%s--%s)", f(med), f(lo), f(hi)) } diff --git a/R/prep_newperson_params.R b/R/prep_newperson_params.R index 7371d58..b55d0a1 100644 --- a/R/prep_newperson_params.R +++ b/R/prep_newperson_params.R @@ -14,8 +14,12 @@ prep_newperson_params <- function(draws_long, antigen_label) { dplyr::mutate(antigen = antigen_label) |> dplyr::select(.data$Iteration, .data$Chain, .data$antigen, .data$Iso_type, .data$Parameter, .data$value) |> - tidyr::pivot_wider(names_from = .data$Parameter, values_from = .data$value) |> - # The model output uses "shape" for the decay-shape parameter, which is renamed to "rho" for consistency + tidyr::pivot_wider( + names_from = .data$Parameter, + values_from = .data$value + ) |> + # The model output uses "shape" for the decay-shape parameter, + # which is renamed to "rho" for consistency dplyr::rename(rho = .data$shape) |> dplyr::select(.data$Iteration, .data$Chain, .data$antigen, .data$Iso_type, .data$y0, .data$y1, .data$t1, .data$alpha, .data$rho) diff --git a/inst/WORDLIST b/inst/WORDLIST index d314679..9cd8da1 100644 --- a/inst/WORDLIST +++ b/inst/WORDLIST @@ -1,7 +1,23 @@ CMD Codecov Github +IgA +IgG +Iso +Isotype +MFI ORCID Postprocess +SOSAR Seroresponse +df +ggplot +ipab +isotype +isotypes +newperson +pre repo +serodynamics +serotype +tibble From e4e6c4235ad1fb64c0a28c42b2707bebbe23ece1 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Tue, 31 Mar 2026 19:32:05 +0000 Subject: [PATCH 24/30] Update man/fmt_mci.Rd to match roxygen docs (devtools::document() sync) Agent-Logs-Url: https://github.com/UCD-SERG/shigella/sessions/a9ff4314-f41a-414a-955e-d0c37e3e3dd8 Co-authored-by: Kwan-Jenny <68584166+Kwan-Jenny@users.noreply.github.com> --- man/fmt_mci.Rd | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/man/fmt_mci.Rd b/man/fmt_mci.Rd index 13129eb..ad0dd2e 100644 --- a/man/fmt_mci.Rd +++ b/man/fmt_mci.Rd @@ -18,7 +18,7 @@ fmt_mci(med, lo, hi, digits = 2, sci = FALSE) \item{sci}{Logical; if TRUE uses scientific notation.} } \value{ -A string like "1.23 (0.50ā€“2.00)". +A string like "1.23 (0.50--2.00)". } \description{ Format median and credible interval as a single string From 190c42f07683d8fabd5c144eabe1291c9c48ead0 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Tue, 31 Mar 2026 20:45:25 +0000 Subject: [PATCH 25/30] Fix review feedback: rlang::sym rename, t1=0 edge case, shigella::: in tests, TODO comment Agent-Logs-Url: https://github.com/UCD-SERG/shigella/sessions/a030fe7b-7a02-4d32-a346-01e663e62bad Co-authored-by: Kwan-Jenny <68584166+Kwan-Jenny@users.noreply.github.com> --- R/compute_residual_metrics.R | 5 ++- R/model_comparison_table.R | 1 + R/utils_internal.R | 15 ++++++- .../test-predict_posterior_at_times.R | 4 +- tests/testthat/test-utils_internal.R | 41 +++++++++++++++---- 5 files changed, 52 insertions(+), 14 deletions(-) diff --git a/R/compute_residual_metrics.R b/R/compute_residual_metrics.R index 815721e..89f124a 100644 --- a/R/compute_residual_metrics.R +++ b/R/compute_residual_metrics.R @@ -83,7 +83,10 @@ compute_residual_metrics <- function(model, value_var <- serocalculator::get_values_var(dataset) observed_data <- dataset |> - dplyr::rename(t = {{ time_var }}, obs = {{ value_var }}) |> + dplyr::rename( + t = !!rlang::sym(time_var), + obs = !!rlang::sym(value_var) + ) |> dplyr::select(dplyr::all_of(c("id", "t", "obs", "antigen_iso"))) |> dplyr::mutate(id = as.character(.data$id)) |> dplyr::filter(.data$id %in% ids, .data$antigen_iso == antigen_iso) diff --git a/R/model_comparison_table.R b/R/model_comparison_table.R index 61a6cfd..4f1418a 100644 --- a/R/model_comparison_table.R +++ b/R/model_comparison_table.R @@ -85,6 +85,7 @@ make_model_comparison_table <- function(model_serospecific, data_serospecific, scale <- match.arg(scale) ids_common <- intersect(unique(data_serospecific$id), unique(data_overall$id)) + # TODO: add early error for empty ids_common m_sero <- compute_residual_metrics( model = model_serospecific, diff --git a/R/utils_internal.R b/R/utils_internal.R index 91ccd9f..cf83f6c 100644 --- a/R/utils_internal.R +++ b/R/utils_internal.R @@ -35,12 +35,23 @@ get_timeindays_var <- function(dataset) { #' @return Numeric vector of predicted antibody levels #' @keywords internal #' @noRd -bt <- function(y0, y1, t1) log(y1 / y0) / t1 +bt <- function(y0, y1, t1) { + if (t1 == 0) return(0) + log(y1 / y0) / t1 +} -# Mirrors serodynamics:::ab() ā€” see ucdavis/serodynamics/R/ab.R +# Mirrors serodynamics:::ab() -- see ucdavis/serodynamics/R/ab.R ab <- function(t, y0, y1, t1, alpha, shape) { # NOTE: this preserves serodynamics behavior; when shape == 1 the decay branch # uses 1 / (1 - shape), matching upstream implementation. + if (t1 == 0) { + # Instant peak: for t <= 0 return y1, then decay from t1 = 0 + return(ifelse( + t <= 0, + y1, + (y1^(1 - shape) - (1 - shape) * alpha * t)^(1 / (1 - shape)) + )) + } beta <- bt(y0, y1, t1) yt <- ifelse( t <= t1, diff --git a/tests/testthat/test-predict_posterior_at_times.R b/tests/testthat/test-predict_posterior_at_times.R index 0cf4b08..1cf05a7 100644 --- a/tests/testthat/test-predict_posterior_at_times.R +++ b/tests/testthat/test-predict_posterior_at_times.R @@ -17,7 +17,7 @@ test_that("predict_posterior_at_times returns expected structure", { times <- c(0, 30, 90) - result <- predict_posterior_at_times( + result <- shigella:::predict_posterior_at_times( model = mock_model, ids = "test_id", antigen_iso = "IgG", @@ -55,7 +55,7 @@ test_that("predict_posterior_at_times handles multiple subjects", { ) ) - result <- predict_posterior_at_times( + result <- shigella:::predict_posterior_at_times( model = mock_model, ids = c("id1", "id2"), antigen_iso = "IgG", diff --git a/tests/testthat/test-utils_internal.R b/tests/testthat/test-utils_internal.R index 5cd71f6..9408c86 100644 --- a/tests/testthat/test-utils_internal.R +++ b/tests/testthat/test-utils_internal.R @@ -1,28 +1,38 @@ test_that("ab() returns y0 at t = 0", { - result <- ab(t = 0, y0 = 200, y1 = 2000, t1 = 10, alpha = 0.02, shape = 1.2) + result <- shigella:::ab( + t = 0, y0 = 200, y1 = 2000, t1 = 10, alpha = 0.02, shape = 1.2 + ) expect_equal(result, 200) }) test_that("ab() returns y1 at t = t1", { - result <- ab(t = 10, y0 = 200, y1 = 2000, t1 = 10, alpha = 0.02, shape = 1.2) + result <- shigella:::ab( + t = 10, y0 = 200, y1 = 2000, t1 = 10, alpha = 0.02, shape = 1.2 + ) expect_equal(result, 2000) }) test_that("ab() shows expected rise for t < t1", { - result <- ab(t = 5, y0 = 200, y1 = 2000, t1 = 10, alpha = 0.02, shape = 1.2) + result <- shigella:::ab( + t = 5, y0 = 200, y1 = 2000, t1 = 10, alpha = 0.02, shape = 1.2 + ) expected <- 200 * exp((log(2000 / 200) / 10) * 5) expect_equal(result, expected) }) test_that("ab() decays toward zero for large t (shape > 1)", { - result <- ab(t = 10000, y0 = 200, y1 = 2000, t1 = 10, alpha = 0.02, shape = 1.2) + result <- shigella:::ab( + t = 10000, y0 = 200, y1 = 2000, t1 = 10, alpha = 0.02, shape = 1.2 + ) expect_true(is.finite(result)) expect_true(result >= 0) }) test_that("ab() handles vector input", { times <- c(0, 5, 10, 30, 90) - result <- ab(t = times, y0 = 200, y1 = 2000, t1 = 10, alpha = 0.02, shape = 1.2) + result <- shigella:::ab( + t = times, y0 = 200, y1 = 2000, t1 = 10, alpha = 0.02, shape = 1.2 + ) expect_length(result, 5) expect_true(all(is.finite(result))) expect_equal(result[1], 200) @@ -30,16 +40,29 @@ test_that("ab() handles vector input", { }) test_that("ab() does not error when t1 = 0", { - expect_no_error(ab(t = c(0, 10, 30), y0 = 200, y1 = 2000, t1 = 0, alpha = 0.02, shape = 1.2)) + expect_no_error( + shigella:::ab( + t = c(0, 10, 30), y0 = 200, y1 = 2000, + t1 = 0, alpha = 0.02, shape = 1.2 + ) + ) +}) + +test_that("ab() returns y1 at t = 0 when t1 = 0 (instant peak)", { + result <- shigella:::ab( + t = 0, y0 = 200, y1 = 2000, t1 = 0, alpha = 0.02, shape = 1.2 + ) + expect_equal(result, 2000) }) test_that("get_timeindays_var() finds time variable from attribute", { mock_data <- data.frame(timepoint = c(0, 7, 30), value = c(100, 200, 300)) attr(mock_data, "timeindays") <- "timepoint" - expect_equal(get_timeindays_var(mock_data), "timepoint") + expect_equal(shigella:::get_timeindays_var(mock_data), "timepoint") }) -test_that("get_timeindays_var() defaults to timeindays when attribute missing", { +test_that( + "get_timeindays_var() defaults to timeindays when attribute missing", { mock_data <- data.frame(timeindays = c(0, 7, 30), value = c(100, 200, 300)) - expect_equal(get_timeindays_var(mock_data), "timeindays") + expect_equal(shigella:::get_timeindays_var(mock_data), "timeindays") }) From c10194148941d527f419abbdfe5517b61575b5ed Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Wed, 1 Apr 2026 00:40:51 +0000 Subject: [PATCH 26/30] Restore ab() and bt() to exact upstream implementations; update t1=0 test Agent-Logs-Url: https://github.com/UCD-SERG/shigella/sessions/40bd4f82-f569-4387-98be-3bfec61029b2 Co-authored-by: Kwan-Jenny <68584166+Kwan-Jenny@users.noreply.github.com> --- R/utils_internal.R | 14 ++------------ tests/testthat/test-utils_internal.R | 20 ++++++++------------ 2 files changed, 10 insertions(+), 24 deletions(-) diff --git a/R/utils_internal.R b/R/utils_internal.R index cf83f6c..33e916c 100644 --- a/R/utils_internal.R +++ b/R/utils_internal.R @@ -36,22 +36,12 @@ get_timeindays_var <- function(dataset) { #' @keywords internal #' @noRd bt <- function(y0, y1, t1) { - if (t1 == 0) return(0) - log(y1 / y0) / t1 + to_return <- log(y1 / y0) / t1 + return(to_return) } # Mirrors serodynamics:::ab() -- see ucdavis/serodynamics/R/ab.R ab <- function(t, y0, y1, t1, alpha, shape) { - # NOTE: this preserves serodynamics behavior; when shape == 1 the decay branch - # uses 1 / (1 - shape), matching upstream implementation. - if (t1 == 0) { - # Instant peak: for t <= 0 return y1, then decay from t1 = 0 - return(ifelse( - t <= 0, - y1, - (y1^(1 - shape) - (1 - shape) * alpha * t)^(1 / (1 - shape)) - )) - } beta <- bt(y0, y1, t1) yt <- ifelse( t <= t1, diff --git a/tests/testthat/test-utils_internal.R b/tests/testthat/test-utils_internal.R index 9408c86..2aac02e 100644 --- a/tests/testthat/test-utils_internal.R +++ b/tests/testthat/test-utils_internal.R @@ -39,20 +39,16 @@ test_that("ab() handles vector input", { expect_equal(result[3], 2000) }) -test_that("ab() does not error when t1 = 0", { - expect_no_error( - shigella:::ab( - t = c(0, 10, 30), y0 = 200, y1 = 2000, - t1 = 0, alpha = 0.02, shape = 1.2 - ) - ) -}) - -test_that("ab() returns y1 at t = 0 when t1 = 0 (instant peak)", { +test_that("ab() handles t1 = 0 same as upstream", { + # When t1 = 0, bt() returns Inf, which is expected + # upstream behavior. We don't add special handling. result <- shigella:::ab( - t = 0, y0 = 200, y1 = 2000, t1 = 0, alpha = 0.02, shape = 1.2 + t = c(0, 10, 30), + y0 = 200, y1 = 2000, t1 = 0, + alpha = 0.02, shape = 1.2 ) - expect_equal(result, 2000) + expect_length(result, 3) + # Results may be NaN/Inf -- that's OK, matches upstream }) test_that("get_timeindays_var() finds time variable from attribute", { From 913eb3b63301cb80ad8e5941eee9c0279157d98d Mon Sep 17 00:00:00 2001 From: Kwan-Jenny <68584166+Kwan-Jenny@users.noreply.github.com> Date: Wed, 1 Apr 2026 17:09:27 -0700 Subject: [PATCH 27/30] Update data-raw/mock_data.R Co-authored-by: Copilot <175728472+Copilot@users.noreply.github.com> --- data-raw/mock_data.R | 7 +++++-- 1 file changed, 5 insertions(+), 2 deletions(-) diff --git a/data-raw/mock_data.R b/data-raw/mock_data.R index ce32604..29c7fdf 100644 --- a/data-raw/mock_data.R +++ b/data-raw/mock_data.R @@ -77,8 +77,11 @@ attr(mock_case_data, "value_var") <- "value" class(mock_case_data) <- c("case_data", "data.frame") # Save datasets ----------------------------------------------------------------- -usethis::use_data(mock_posterior_draws, overwrite = TRUE) -usethis::use_data(mock_case_data, overwrite = TRUE) +if (!dir.exists("data")) { + dir.create("data", recursive = TRUE) +} +save(mock_posterior_draws, file = file.path("data", "mock_posterior_draws.rda")) +save(mock_case_data, file = file.path("data", "mock_case_data.rda")) message("Mock datasets created successfully!") message("- mock_posterior_draws: ", nrow(mock_posterior_draws), " rows") From 1fa288b04c66c0bcff841b9b277d11bd29ed3bb0 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Thu, 2 Apr 2026 00:10:56 +0000 Subject: [PATCH 28/30] Fix roxygen docs: give bt() its own params, move trajectory docs to ab() Agent-Logs-Url: https://github.com/UCD-SERG/shigella/sessions/9c7dfb28-b082-4292-98b7-37fbefb0b36b Co-authored-by: Kwan-Jenny <68584166+Kwan-Jenny@users.noreply.github.com> --- R/utils_internal.R | 25 ++++++++++++++++++------- 1 file changed, 18 insertions(+), 7 deletions(-) diff --git a/R/utils_internal.R b/R/utils_internal.R index 33e916c..ed966a0 100644 --- a/R/utils_internal.R +++ b/R/utils_internal.R @@ -22,17 +22,15 @@ get_timeindays_var <- function(dataset) { return(var) } -#' Antibody kinetics model function +#' Compute log-linear rise rate for antibody kinetics #' -#' Evaluates the antibody trajectory model at specified time points. +#' Helper that computes the exponential rise rate beta = log(y1/y0) / t1, +#' used internally by `ab()`. #' -#' @param t Numeric vector of time points (days) #' @param y0 Baseline antibody level #' @param y1 Peak antibody level #' @param t1 Time to peak (days) -#' @param alpha Decay rate parameter -#' @param shape Decay shape parameter (rho) -#' @return Numeric vector of predicted antibody levels +#' @return Numeric scalar: the log-linear rise rate #' @keywords internal #' @noRd bt <- function(y0, y1, t1) { @@ -40,7 +38,20 @@ bt <- function(y0, y1, t1) { return(to_return) } -# Mirrors serodynamics:::ab() -- see ucdavis/serodynamics/R/ab.R +#' Antibody kinetics trajectory function +#' +#' Evaluates the antibody trajectory model at specified time points. +#' Mirrors `serodynamics:::ab()` -- see ucdavis/serodynamics/R/ab.R. +#' +#' @param t Numeric vector of time points (days) +#' @param y0 Baseline antibody level +#' @param y1 Peak antibody level +#' @param t1 Time to peak (days) +#' @param alpha Decay rate parameter +#' @param shape Decay shape parameter (rho) +#' @return Numeric vector of predicted antibody levels +#' @keywords internal +#' @noRd ab <- function(t, y0, y1, t1, alpha, shape) { beta <- bt(y0, y1, t1) yt <- ifelse( From 6b11bbf4e1cdf339647968a506675d0945870e51 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Thu, 2 Apr 2026 00:16:15 +0000 Subject: [PATCH 29/30] Fix Github -> GitHub capitalization in inst/WORDLIST Agent-Logs-Url: https://github.com/UCD-SERG/shigella/sessions/be13d809-cfdc-485c-889c-92a9d2ebd9c2 Co-authored-by: Kwan-Jenny <68584166+Kwan-Jenny@users.noreply.github.com> --- inst/WORDLIST | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/inst/WORDLIST b/inst/WORDLIST index 9cd8da1..d92948f 100644 --- a/inst/WORDLIST +++ b/inst/WORDLIST @@ -1,6 +1,6 @@ CMD Codecov -Github +GitHub IgA IgG Iso From 9059965ba5a2c52068661afe13c17fd45dec8b53 Mon Sep 17 00:00:00 2001 From: "copilot-swe-agent[bot]" <198982749+Copilot@users.noreply.github.com> Date: Thu, 2 Apr 2026 05:52:14 +0000 Subject: [PATCH 30/30] Fix Github -> GitHub capitalization in NEWS.md Agent-Logs-Url: https://github.com/UCD-SERG/shigella/sessions/9263f1b7-8152-4f90-9f85-e2bb428c4798 Co-authored-by: Kwan-Jenny <68584166+Kwan-Jenny@users.noreply.github.com> --- NEWS.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/NEWS.md b/NEWS.md index 637844b..943bd49 100644 --- a/NEWS.md +++ b/NEWS.md @@ -10,7 +10,7 @@ * added `.lintr` configuration file * Moved helper functions into `R/` subdirectory (#1) -* simplified spell-check Github Action (#3) +* simplified spell-check GitHub Action (#3) # shigella 0.0.0