Added functions to simulate and plot the antibody model easily

NAMESPACE

@@ -1,5 +1,6 @@
 # Generated by roxygen2: do not edit by hand
 
+export()
 export(add_measurement_shifts)
 export(add_noise)
 export(antibody_model)
@@ -68,6 +69,7 @@ export(pad_infection_model_prior_parameters)
 export(plot_2d_density)
 export(plot_antibody_data)
 export(plot_antibody_dependent_boosting)
+export(plot_antibody_model)
 export(plot_attack_rates)
 export(plot_attack_rates_pointrange)
 export(plot_cumulative_infection_histories)
@@ -94,6 +96,7 @@ export(setup_antibody_data_for_posterior_func)
 export(setup_infection_histories)
 export(setup_infection_histories_antibody_level)
 export(setup_infection_histories_prior)
+export(simulate_antibody_model)
 export(simulate_attack_rates)
 export(simulate_data)
 export(simulate_group)

R/plot_antibody_model.R

@@ -27,6 +27,56 @@ plot_antibody_dependent_boosting <- function(chain, n, titres = seq(0, 8, by = 0
   return(range)
 }
 
+#' Plot the antibody model
+#' 
+#' Plots the trajectory of the serosolver antibody model using specified parameters and optionally a specified antigenic map and infection history.
+#' @inheritParams simulate_antibody_model
+#' @return a list with two ggplot objects, one showing the simulated antibody kinetics over time, stratified by biomarker ID, the other showing simulated antibody kinetics for each biomarker ID, stratified by sample time
+#' @examples
+#' plot_antibody_model(c("boost_long"=2,"boost_short"=3,"boost_delay"=1,"wane_short"=0.2,"wane_long"=0.01, "antigenic_seniority"=0,"cr_long"=0.1,"cr_short"=0.03), times=seq(1,25,by=1),infection_history=NULL,antigenic_map=example_antigenic_map)
+#'  
+#' @export
+plot_antibody_model <- function(pars, 
+                                times=NULL, 
+                                infection_history=NULL, 
+                                antigenic_map=NULL){
+  y <- simulate_antibody_model(pars,times, infection_history,antigenic_map)
+  y$biomarker_id_label <- paste0("Biomarker ID: ", y$biomarker_ids)
+  y$sample_label <- paste0("Sample time: ", y$sample_times)
+
+  y$biomarker_id_label <- factor(y$biomarker_id_label, levels=unique(paste0("Biomarker ID: ", y$biomarker_ids)))
+  y$sample_label <- factor(y$sample_label, levels=unique(paste0("Sample time: ", y$sample_times)))
+
+  if(is.null(infection_history)){
+    infection_history <- 1
+  }
+
+  p_long <- ggplot(y) + 
+    geom_line(aes(x=sample_times,y=antibody_level)) + 
+    geom_vline(data=data.frame(x=infection_history),linetype="dashed",aes(col="Infection",xintercept=x)) +
+    facet_wrap(~biomarker_id_label) +
+    scale_color_manual(name="",values=c("Infection"="grey40")) +
+    scale_y_continuous(expand=c(0,0)) +
+    scale_x_continuous(expand=c(0,0)) +
+    xlab("Sample time") +
+    ylab("Antibody level") +
+    theme_pubr()+
+    theme(legend.position="bottom")
+
+  p_cr <- ggplot(y) + 
+    geom_ribbon(aes(x=biomarker_ids,ymax=antibody_level,ymin=0),fill='grey70',col="black") + 
+    geom_vline(data=data.frame(x=infection_history),linetype="dashed",aes(col="Infection",xintercept=x)) +
+    facet_wrap(~sample_label) +
+    scale_color_manual(name="", values=c("Infection"="grey40")) +
+    scale_y_continuous(expand=c(0,0)) +
+    scale_x_continuous(expand=c(0,0)) +
+    xlab("Biomarker ID") +
+    ylab("Antibody level") +
+    theme_pubr()+
+    theme(legend.position="bottom")
+
+  return(list(p_long,p_cr))
+}
 
 #' Plots infection histories and antibody model 
 #'

R/posteriors.R

@@ -14,7 +14,10 @@
 #' @param n_alive if not NULL, uses this as the number alive in a given year rather than calculating from the ages. This is needed if the number of alive individuals is known, but individual birth dates are not
 #' @param function_type integer specifying which version of this function to use. Specify 1 to give a posterior solving function; 2 to give the gibbs sampler for infection history proposals; otherwise just solves the titre model and returns predicted titres. NOTE that this is not the same as the attack rate prior argument, \code{version}!
 #' @param titre_before_infection TRUE/FALSE value. If TRUE, solves titre predictions, but gives the predicted titre at a given time point BEFORE any infection during that time occurs.
-#' @param data_type integer, currently accepting 1 for discrete or 2 for continuous. 
+#' @param data_type integer or vector, with an entry for each unique data type in `antibody_data`. Set to 1 for discrete data (e.g., fold dilution) or 2 for continuous (e.g., ELISA optical density). 
+#' @param biomarker_groups_weights integer or vector, giving a factor to multiply the log-likelihood contribution of this data type towards the overall likelihood.
+#' @param start_level character, to tell the model how to treat initial antibody levels. This uses the observed data to either select starting values for each unique `individual`, `biomarker_id` and `biomarker_group` combination. See \code{\link{create_start_level_data}}. One of "min", "max", "mean", "median", or "full_random". Any other entry assumes all antibody starting levels are set to 0. Can also pass a tibble or data frame of starting levels matching the output of \code{\link{create_start_level_data}}.
+#' @param start_level_randomize if TRUE, and data is discretized, then sets the starting antibody level to a random value between floor(x) and floor(x) + 1. Does nothing if using continuous data.
 #' @param verbose if TRUE, prints warning messages
 #' @param ... other arguments to pass to the posterior solving function
 #' @return a single function pointer that takes only pars and infection_histories as unnamed arguments. This function goes on to return a vector of posterior values for each individual
@@ -48,11 +51,10 @@ create_posterior_func <- function(par_tab,
                                   antibody_level_before_infection=FALSE,
                                   data_type=1,
                                   biomarker_groups_weights =1,
-                                  verbose=FALSE,
                                   start_level = "other",
                                   start_level_randomize=FALSE,
+                                  verbose=FALSE,
                                   ...) {
-  #browser()
 
     check_par_tab(par_tab, TRUE, prior_version,verbose)
     if (!("population_group" %in% colnames(antibody_data))) {
@@ -75,7 +77,7 @@ create_posterior_func <- function(par_tab,
     }
 
     if(verbose){
-      message(cat("Setting starting antibody levels based on data using command:", start_level, ";and randomizing starting antibody levels set to:", start_level_randomize))
+      message(cat("Setting starting antibody levels based on data using command:", start_level, "; and randomizing starting antibody levels set to:", start_level_randomize))
     }
 
     ## Check that antibody data is formatted correctly
@@ -194,8 +196,14 @@ create_posterior_func <- function(par_tab,
 
     ## Starting antibody levels
     births <- antibody_data_unique$birth
-    start_levels <- create_start_level_data(antibody_data,start_level,start_level_randomize) %>% 
-      arrange(individual, biomarker_group, sample_time, biomarker_id, repeat_number)
+
+    ## If specifying starting levels based on a provided data frame or tibble, set them here. Otherwise, calculate them from the antibody data.
+    if(class(start_level) %in% c("data.frame","tibble")){
+      start_levels <- start_level 
+    } else {
+      start_levels <- create_start_level_data(antibody_data,start_level,start_level_randomize) %>% 
+        arrange(individual, biomarker_group, sample_time, biomarker_id, repeat_number)
+    }
 
     start_levels <- start_levels %>% filter(repeat_number == 1)
     #browser()

R/simulate_data.R

@@ -605,6 +605,68 @@ simulate_infection_histories <- function(p_inf, possible_exposure_times, samplin
   return(list(infection_histories, ARs))
 }
 
+#' Simulate the antibody model
+#' 
+#' Simulates the trajectory of the serosolver antibody model using specified parameters and optionally a specified antigenic map and infection history.
+#' @param pars the vector of named model parameters, including `boost_long`, `boost_short`,`boost_delay`,`wane_long`,`wane_short`,`cr_long`, and `cr_short`.
+#' @param times the vector of times to solve the model over. A continuous vector of discrete timepoints. Can be left to NULL if this information is included in the `antigenic_map` argument.
+#' @param infection_history the vector of times matching entries in `times` to simulate infections in.
+#' @param antigenic_map the antigenic map to solve the model with. Can be left to NULL to ssume all biomarker IDs have the same antigenic coordinates.
+#' @return a data frame with variables `sample_times`, `biomarker_id` and `antibody_level`
+#' @examples
+#' simulate_antibody_model(c("boost_long"=2,"boost_short"=3,"boost_delay"=1,"wane_short"=0.2,"wane_long"=0.01, "antigenic_seniority"=0,"cr_long"=0.1,"cr_short"=0.03), times=seq(1,25,by=1),infection_history=NULL,antigenic_map=example_antigenic_map)
+#'  
+#' @export
+simulate_antibody_model <- function(pars, 
+                                times=NULL, 
+                                infection_history=NULL, 
+                                antigenic_map=NULL){
+  if(is.null(times) & is.null(antigenic_map)){
+    stop("Must provide one of times or antigenic_map to give the possible infection times and biomarker IDs over which to solve the model.")
+  }
+
+  ## If no vector of times provided, take from the antigenic map
+  if(is.null(times) & !is.null(antigenic_map)){
+    times <- antigenic_map$inf_times
+  }
+
+  ## If no antigenic map is provided, create a dummy antigenic map where each element has the same antigenic coordinate
+  if(is.null(antigenic_map)){
+    antigenic_map <- data.frame(x_coord=1,y_coord=1,inf_times=times)
+    biomarker_ids <- 0
+  } else {
+    biomarker_ids <- match(antigenic_map$inf_times[seq_along(times)], antigenic_map$inf_times[seq_along(times)]) -1
+  }
+
+  ## Setup antigenic map
+  use_antigenic_map <- melt_antigenic_coords(antigenic_map[seq_along(times),c("x_coord","y_coord")])
+  antigenic_map_long <- matrix(create_cross_reactivity_vector(use_antigenic_map, pars["cr_long"]),ncol=1)
+  antigenic_map_short <- matrix(create_cross_reactivity_vector(use_antigenic_map, pars["cr_short"]),ncol=1)
+
+  ## If no infection history was provided, setup a dummy infection history vector with only the first entry as an infection
+  if(is.null(infection_history)){
+    infection_history <- 1
+  }
+
+  infection_indices <- infection_history-1
+  sample_times <- seq(1, max(times),by=1)
+
+
+  start_levels <- rep(0,length(biomarker_ids))
+
+  y <- antibody_model_individual_wrapper(pars["boost_long"],pars["boost_short"],pars["boost_delay"],
+                                         pars["wane_short"],pars["wane_long"],pars["antigenic_seniority"],
+                                         0,
+                                         start_levels,
+                                         length(times),
+                                         infection_history,
+                                         infection_indices,
+                                         biomarker_ids,
+                                         sample_times,
+                                         antigenic_map_long,
+                                         antigenic_map_short)
+  data.frame(sample_times = rep(sample_times,each=length(biomarker_ids)),biomarker_ids = rep(biomarker_ids,length(sample_times)),antibody_level=y)
+}
 
 #' Generates attack rates from an SIR model with fixed beta/gamma, specified final attack rate and the number of time "buckets" to solve over ie. buckets=12 returns attack rates for 12 time periods
 generate_ar_annual <- function(AR, buckets) {

inst/extdata/scripts/adding_fixed_starting_titers.R

@@ -22,7 +22,18 @@ start_titres <- runif(length(biomarker_ids), 0,2)
 biomarker_ids <- seq(0,24,by=4)
 
 print(start_titres)
-y <- antibody_model_individual_wrapper(3,2,1,0.5,0.02,0.05,3,
+
+plot_antibody_model(c("boost_long"=2,"boost_short"=3,"boost_delay"=1,"wane_short"=0.2,"wane_long"=0.01, "antigenic_seniority"=0,"cr_long"=0.1,"cr_short"=0.03), times=seq(1,25,by=1),infection_history=NULL,
+
+                    antigenic_map=example_antigenic_map)
+
+y <- simulate_antibody_model(c("boost_long"=2,"boost_short"=3,"boost_delay"=1,"wane_short"=0.2,"wane_long"=0.01, "antigenic_seniority"=0,"cr_long"=0.1,"cr_short"=0.03), times=seq(1,25,by=1),infection_history=NULL,
+
+                        antigenic_map=NULL)
+ggplot(y) + geom_line(aes(x=sample_times,y=antibody_level)) + facet_wrap(~biomarker_ids)
+
+y <- antibody_model_individual_wrapper(3,2,1,0.5,0.02,0.05,
+                                       3,
                                        start_titres,
                                        length(exposure_times),
                                   infection_times,
@@ -46,7 +57,7 @@ example_inf_hist1 <- matrix(0, nrow=50,ncol=48)
 antibody_data <- example_antibody_data #%>% filter(repeat_number == 1)
 antibody_data <- antibody_data %>% arrange(individual, biomarker_group, sample_time, biomarker_id, repeat_number)
 f <- create_posterior_func(par_tab,antibody_data,example_antigenic_map,function_type=4,
-                           start_level="max")
+                           start_level="mean")
 
 antibody_data$y <- f(par_tab$values, example_inf_hist1)
 
@@ -72,13 +83,14 @@ library(coda)
 serosolver(par_tab, example_antibody_data, example_antigenic_map,
                        filename="readme", prior_version=2,n_chains=1,parallel=FALSE,
                        mcmc_pars=c(adaptive_iterations=10000, iterations=50000),verbose=TRUE,
-           start_level="min")
+           start_level="none")
 
-chains <- load_mcmc_chains(location=getwd(),par_tab=example_par_tab,burnin = 100000,unfixed=TRUE)
+chains <- load_mcmc_chains(location=getwd(),par_tab=par_tab,burnin = 10000,unfixed=TRUE)
 plot_model_fits(chain = chains$theta_chain,
                 infection_histories = chains$inf_chain,
-                known_infection_history = example_inf_hist,
-                antibody_data = example_antibody_data,individuals=1:4,
+                known_infection_history = NULL,
+                antibody_data = example_antibody_data,individuals=1:5,
                 antigenic_map=example_antigenic_map,
-                par_tab=example_par_tab,
+                nsamp=100,
+                par_tab=par_tab,
                 orientation="cross-sectional")