diff --git a/DESCRIPTION b/DESCRIPTION index 58fdf6f4..3ffcc904 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -44,7 +44,6 @@ Imports: reshape2, bayesplot, loo, - cowplot, Hmisc, dplyr, gsubfn, @@ -58,7 +57,6 @@ Imports: BH, RcppEigen, RcppParallel, - pracma, purrr Suggests: knitr, diff --git a/NAMESPACE b/NAMESPACE index 1c975746..ea0f3e9d 100644 --- a/NAMESPACE +++ b/NAMESPACE @@ -1,25 +1,20 @@ # Generated by roxygen2: do not edit by hand -export(extract_seroprev_model_summary) -export(fit_seroprev_model) -export(generate_sim_data) -export(get_comparison_table) +export(extract_seromodel_summary) +export(fit_seromodel) export(get_exposure_ages) export(get_exposure_matrix) export(get_prev_expanded) -export(get_sim_counts) export(get_table_rhats) -export(group_sim_data) export(plot_foi) export(plot_info_table) export(plot_rhats) +export(plot_seromodel) export(plot_seroprev) export(plot_seroprev_fitted) -export(plot_seroprev_model) -export(plot_seroprev_models_grid) export(prepare_bin_data) -export(prepare_seroprev_data) -export(run_seroprev_model) +export(prepare_serodata) +export(run_seromodel) export(save_or_load_model) import(Rcpp) import(dplyr) diff --git a/R/model_comparison.R b/R/model_comparison.R index 0e62231b..3c62814b 100644 --- a/R/model_comparison.R +++ b/R/model_comparison.R @@ -1,105 +1,27 @@ -#' Get Table Rhats -#' -#' Function that makes the rhats table -#' @param model_object model_object +#' Method for extracting a dataframe containing the R-hat estimates for a given serological model. +#' +#' This method relies in the function \link[bayesplot]{rhat} to extract the R-hat estimates of the serological model object +#' \code{seromodel_object} and returns a table a dataframe with the estimates for each year of birth. +#' @param seromodel_object seromodel_object #' @return rhats table #' @examples #' \dontrun{ -#' seroprev_data <- prepare_seroprev_data(seroprev_data = serodata, alpha = 0.05) -#' model_object <- run_seroprev_model( -#' seroprev_data = seroprev_data, seroprev_model_name = "constant_foi_bi") -#' get_table_rhats (model_object) +#' data("serodata") +#' data_test <- prepare_serodata(serodata = serodata) +#' model_constant <- run_seromodel(serodata = data_test, +#' seromodel_name = "constant_foi_bi", +#' n_iters = 1500) +#' get_table_rhats(model_object = model_constant) #' } #' @export -get_table_rhats <- function(model_object) { - rhats <- bayesplot::rhat(model_object$fit, "foi") +get_table_rhats <- function(seromodel_object) { + rhats <- bayesplot::rhat(seromodel_object$fit, "foi") if (any(is.nan(rhats))) { rhats[which(is.nan(rhats))] <- 0 } - model_rhats <- data.frame(year = model_object$exposure_years, rhat = rhats) + model_rhats <- data.frame(year = seromodel_object$exposure_years, rhat = rhats) model_rhats$rhat[model_rhats$rhat == 0] <- NA return(model_rhats) -} - -#' Get Model Table Comparison -#' Provides a table with statistics for comparison between models and selection -#' @param model_objects_list model_objects to compare -#' @return comparison table -#' @examples -#' \dontrun{ -#' data_test <- prepare_seroprev_data(serodata) -#' model_0 <- run_seroprev_model(seroprev_data = data_test, -#' seroprev_model_name = "constant_foi_bi", -#' n_iters = 1000) -#' -#' model_1 <- run_seroprev_model(seroprev_data = data_test, -#' seroprev_model_name = "continuous_foi_normal_bi", -#' n_iters = 1000) -#' -#' model_2 <- run_seroprev_model(seroprev_data = data_test, -#' seroprev_model_name = "continuous_foi_normal_log", -#' n_iters = 1000) -#' comp_table <- get_comparison_table(model_objects_list = c(m0 = model_0, -#' m1 = model_1, -#' m2 = model_2)) -#' } -#' @export -get_comparison_table <- function(model_objects_list) { - - - dif_m0_m1 <- loo::loo_compare(model_objects_list$m0.loo_fit, - model_objects_list$m1.loo_fit) - - dif_m0_m2 <- loo::loo_compare(model_objects_list$m0.loo_fit, - model_objects_list$m2.loo_fit) - - # Aquí pendiente revisar desde la función summary_model - # No estoy segura que este parámetro venga bien desde allá ni tampoco que esté bien acá - - # model_comp$better <- NA - # model_comp$better[model_comp$difference > 0] <- 'Yes' - # model_comp$better[model_comp$difference <= 0] <-'No' - # model_comp$better[model_comp$model == 'constant_foi_bi'] <- "-" - - - model_objects_list$m0.model_summary$difference <- 0 - model_objects_list$m0.model_summary$diff_se <- 1 - - model_objects_list$m1.model_summary$difference <- dif_m0_m1[1] - model_objects_list$m1.model_summary$diff_se <- dif_m0_m1[2] - - model_objects_list$m2.model_summary$difference <- dif_m0_m2[1] - model_objects_list$m2.model_summary$diff_se <- dif_m0_m2[2] - - model_comp <- rbind(model_objects_list$m0.model_summary, - model_objects_list$m1.model_summary, - model_objects_list$m2.model_summary) - - model_comp$converged[model_comp$elpd == -1.000e+10] <- 'No' # - ds_one <- dplyr::filter(model_comp, converged == 'Yes') - print(paste0('number of converged models = ', NROW(ds_one))) - - # Ordering the best model based on elpd values - elps_order <- rev(sort(ds_one$elpd)) - best <- dplyr::filter(model_comp, elpd %in% elps_order) %>% dplyr::arrange(-.data$elpd)# This is to make sure I keep only three - best_model1 <- as.character(best$model[1]) - best_model2 <- as.character(best$model[2]) - best_model3 <- as.character(best$model[3]) - - model_comp$best_elpd <- NA - model_comp$best_elpd[model_comp$model == best_model1] <- 1 - model_comp$best_elpd[model_comp$model == best_model2] <- 2 - model_comp$best_elpd[model_comp$model == best_model3] <- 3 - model_comp <- model_comp %>% dplyr::arrange(.data$best_elpd) - - # Estimating p-values to check the difference between the models m0 and other models is actually important - model_comp <- model_comp %>% dplyr::mutate(pvalue = 1 - stats::pnorm(difference/diff_se,0,1)) - model_comp$pvalue[is.nan(model_comp$pvalue)] <- 1 - model_comp$pvalue <- model_comp$pvalue * stats::runif(NROW(model_comp), min = 1, max = 1.0001)# I make this just to ensure I get different values - model_comp$pvalue[model_comp$model == 'constant_foi_bi'] <- 0 - model_comp$pvalue <- round(model_comp$pvalue, 6) - - return(model_comp) -} +} \ No newline at end of file diff --git a/R/modelling.R b/R/modelling.R index a0c5bea4..aceb1adb 100644 --- a/R/modelling.R +++ b/R/modelling.R @@ -1,7 +1,9 @@ -# TODO For some reason, the examples cannot access the serodata variable -#' Run the specified stan model for the force-of-infection +#' Run the specified stan model for the Force-of-Infection and estimates de seroprevalence based on the result of the fit. +#' +#' This function runs the specified model for the Force-of-Infection \code{seromodel_name} using the data froma seroprevalence survey +#' \code{serodata} as the input data. See \link{fit_seromodel} for further details. #' -#' @param seroprev_data A data frame containing the data from a seroprevalence survey. +#' @param serodata A data frame containing the data from a seroprevalence survey. #' This data frame must contain the following columns: #' \tabular{ll}{ #' \code{survey} \tab survey Label of the current survey \cr \tab \cr @@ -22,8 +24,8 @@ #' \code{prev_obs_lower} \tab Lower limit of the confidence interval for the observed prevalence \cr \tab \cr #' \code{prev_obs_upper} \tab Upper limit of the confidence interval for the observed prevalence \cr \tab \cr #' } -#' The last six colums can be added to \code{seroprev_data} by means of the function \code{\link{prepare_seroprev_data}}. -#' @param seroprev_model_name Name of the selected model. Current version provides three options: +#' The last six colums can be added to \code{serodata} by means of the function \code{\link{prepare_serodata}}. +#' @param seromodel_name Name of the selected model. Current version provides three options: #' \describe{ #' \item{\code{"constant_foi_bi"}}{Runs a constant model} #' \item{\code{"continuous_foi_normal_bi"}}{Runs a normal model} @@ -36,42 +38,46 @@ #' For further details refer to the \code{control} parameter in \link[rstan]{sampling} or \href{https://mc-stan.org/rstanarm/reference/adapt_delta.html}{here}. #' @param m_treed Maximum tree depth for the binary tree used in the NUTS stan sampler. For further details refer to the \code{control} parameter in \link[rstan]{sampling}. #' @param decades Number of decades covered by the survey data. -#' @return \code{model_object}. An object containing relevant information about the implementation of the model. For further details refer to \link{fit_seroprev_model}. +#' @return \code{seromodel_object}. An object containing relevant information about the implementation of the model. For further details refer to \link{fit_seromodel}. #' @examples #' \dontrun{ -#' seroprev_data <- prepare_seroprev_data(serodata) -#' run_seroprev_model (seroprev_data, -#' seroprev_model_name = "constant_foi_bi") +#' serodata <- prepare_serodata(serodata) +#' run_seromodel (serodata, +#' seromodel_name = "constant_foi_bi") #' } #' @export -run_seroprev_model <- function(seroprev_data, - seroprev_model_name = "constant_foi_bi", - n_iters = 1000, - n_thin = 2, - delta = 0.90, - m_treed = 10, - decades = 0) { - survey <- unique(seroprev_data$survey) +run_seromodel <- function(serodata, + seromodel_name = "constant_foi_bi", + n_iters = 1000, + n_thin = 2, + delta = 0.90, + m_treed = 10, + decades = 0, + print_summary = TRUE) { + survey <- unique(serodata$survey) if (length(survey) > 1) warning("You have more than 1 surveys or survey codes") - model_object <- fit_seroprev_model(seroprev_data = seroprev_data, - seroprev_model_name = seroprev_model_name, - n_iters = n_iters, - n_thin = n_thin, - delta = delta, - m_treed = m_treed, - decades = decades); print(paste0("serofoi model ", - seroprev_model_name, - " finished running ------")) - print(t(model_object$model_summary)) - return(model_object) + seromodel_object <- fit_seromodel(serodata = serodata, + seromodel_name = seromodel_name, + n_iters = n_iters, + n_thin = n_thin, + delta = delta, + m_treed = m_treed, + decades = decades); print(paste0("serofoi model ", + seromodel_name, + " finished running ------")) + if (print_summary){ + print(t(seromodel_object$model_summary)) + } + return(seromodel_object) } -#' Save or load model +# TODO The warning 'recompiling to avoid crashing R session' still appears when the function is run for a second time. +#' Auxiliar function used to determine whether the stan model corresponding to the specified serological model has been already compiled or not. #' #' This function determines whether the corresponding .RDS file of the selected model exists or not. #' In case the .RDS file exists, it is read and returned; otherwise, the object model is created through the #' \link[rstan]{stan_model} function, saved as an .RDS file and returned as the output of the function. -#' @param seroprev_model_name Name of the selected model. Current version provides three options: +#' @param seromodel_name Name of the selected model. Current version provides three options: #' \describe{ #' \item{\code{"constant_foi_bi"}}{Runs a constant model} #' \item{\code{"continuous_foi_normal_bi"}}{Runs a normal model} @@ -79,17 +85,20 @@ run_seroprev_model <- function(seroprev_data, #' } #' @return \code{model}. The rstan model object corresponding to the selected model. #' @examples -#' save_or_load_model(seroprev_model_name="constant_foi_bi") +#' \dontrun{ +#' model <- save_or_load_model(seromodel_name="constant_foi_bi") +#' } +#' #' @export -save_or_load_model <- function(seroprev_model_name = "constant_foi_bi") { +save_or_load_model <- function(seromodel_name = "constant_foi_bi") { base_path <- config::get("stan_models_base_path", file = system.file("config.yml", package = "serofoi", mustWork = TRUE)) - rds_path <- system.file(base_path, paste(seroprev_model_name, ".rds", sep = ""), package = getPackageName()) + rds_path <- system.file(base_path, paste(seromodel_name, ".rds", sep = ""), package = getPackageName()) if (!file.exists(rds_path)) { - message(sprintf("\nNo rds file found for model %s. Compiling stan model...", seroprev_model_name)) + message(sprintf("\nNo rds file found for model %s. Compiling stan model...", seromodel_name)) } - stan_path <- system.file(base_path, paste(seroprev_model_name, ".stan", sep = ""), package = getPackageName()) + stan_path <- system.file(base_path, paste(seromodel_name, ".stan", sep = ""), package = getPackageName()) model <- rstan::stan_model(stan_path, auto_write = TRUE) @@ -97,11 +106,13 @@ save_or_load_model <- function(seroprev_model_name = "constant_foi_bi") { } -#' Fit Model -#' -#' Function that fits the selected model to the data -#' @param seroprev_data A data frame containing the data from a seroprevalence survey. For further details refer to \link{run_seroprev_model}. -#' @param seroprev_model_name Name of the selected model. Current version provides three options: +#' Function that fits the selected model to the specified seroprevalence survey data. +#' +#' This function fits the specified model \code{seromodel_name} to the serological survey data \code{serodata} +#' by means of the \link[rstan]{sampling} method. The function determines whether the corresponding stan model +#' object needs to be compiled by means of the function \link{save_or_load_model}. +#' @param serodata A data frame containing the data from a seroprevalence survey. For further details refer to \link{run_seromodel}. +#' @param seromodel_name Name of the selected model. Current version provides three options: #' \describe{ #' \item{\code{"constant_foi_bi"}}{Runs a constant model} #' \item{\code{"continuous_foi_normal_bi"}}{Runs a normal model} @@ -114,10 +125,10 @@ save_or_load_model <- function(seroprev_model_name = "constant_foi_bi") { #' For further details refer to the \code{control} parameter in \link[rstan]{sampling} or \href{https://mc-stan.org/rstanarm/reference/adapt_delta.html}{here}. #' @param m_treed Maximum tree depth for the binary tree used in the NUTS stan sampler. For further details refer to the \code{control} parameter in \link[rstan]{sampling}. #' @param decades Number of decades covered by the survey data. -#' @return \code{model_object}. An object containing relevant information about the implementation of the model. It contains the following: +#' @return \code{seromodel_object}. An object containing relevant information about the implementation of the model. It contains the following: #' \tabular{ll}{ #' \code{fit} \tab \code{stanfit} object returned by the function \link[rstan]{sampling} \cr \tab \cr -#' \code{seroprev_data} \tab A data frame containing the data from a seroprevalence survey. For further details refer to \link{run_seroprev_model}.\cr \tab \cr +#' \code{serodata} \tab A data frame containing the data from a seroprevalence survey. For further details refer to \link{run_seromodel}.\cr \tab \cr #' \code{stan_data} \tab List containing \code{Nobs}, \code{Npos}, \code{Ntotal}, \code{Age}, \code{Ymax}, \code{AgeExpoMatrix} and \code{NDecades}. #' This object is used as an input for the \link[rstan]{sampling} function \cr \tab \cr #' \code{exposure_years} \tab Integer atomic vector containing the actual exposure years (1946, ..., 2007 e.g.) \cr \tab \cr @@ -125,42 +136,43 @@ save_or_load_model <- function(seroprev_model_name = "constant_foi_bi") { #' \code{n_iters} \tab Number of interations for eah chain including the warmup. \cr \tab \cr #' \code{n_thin} \tab Positive integer specifying the period for saving samples. \cr \tab \cr #' \code{n_warmup} \tab Number of warm up iterations. Set by default as n_iters/2. \cr \tab \cr -#' \code{seroprev_model_name} \tab The name of the model\cr \tab \cr +#' \code{seromodel_name} \tab The name of the model\cr \tab \cr #' \code{delta} \tab Real number between 0 and 1 that represents the target average acceptance probability. \cr \tab \cr #' \code{m_treed} \tab Maximum tree depth for the binary tree used in the NUTS stan sampler. \cr \tab \cr #' \code{loo_fit} \tab Efficient approximate leave-one-out cross-validation. Refer to \link[loo]{loo} for further details. \cr \tab \cr #' \code{foi_cent_est} \tab A data fram e containing \code{year} (corresponding to \code{exposure_years}), \code{lower}, \code{upper}, and \code{medianv} \cr \tab \cr #' \code{foi_post_s} \tab Sample n rows from a table. Refer to \link[dplyr]{sample_n} for further details. \cr \tab \cr -#' \code{model_summary} \tab A data fram containing the summary of the model. Refer to \link{extract_seroprev_model_summary} for further details. \cr \tab \cr +#' \code{model_summary} \tab A data fram containing the summary of the model. Refer to \link{extract_seromodel_summary} for further details. \cr \tab \cr #' } #' @examples #' \dontrun{ -#' seroprev_data <- prepare_seroprev_data(serodata) -#' fit_seroprev_model (seroprev_data, -#' seroprev_model_name = "constant_foi_bi") +#' data("serodata") +#' serodata <- prepare_serodata(serodata) +#' seromodel_fit <- fit_seromodel(serodata = serodata, +#' seromodel_name = "constant_foi_bi") #' } #' #' @export -fit_seroprev_model <- function(seroprev_data, - seroprev_model_name, - n_iters = 1000, - n_thin = 2, - delta = 0.90, - m_treed = 10, - decades = 0) { - # add a warning because there are exceptions where a minimal amount of iterations need to be run - model <- save_or_load_model(seroprev_model_name) - exposure_ages <- get_exposure_ages(seroprev_data) - exposure_years <- (min(seroprev_data$birth_year):seroprev_data$tsur[1])[-1] - exposure_matrix <- get_exposure_matrix(seroprev_data) - Nobs <- nrow(seroprev_data) +fit_seromodel <- function(serodata, + seromodel_name, + n_iters = 1000, + n_thin = 2, + delta = 0.90, + m_treed = 10, + decades = 0) { + # TODO Add a warning because there are exceptions where a minimal amount of iterations is needed. + model <- save_or_load_model(seromodel_name) + exposure_ages <- get_exposure_ages(serodata) + exposure_years <- (min(serodata$birth_year):serodata$tsur[1])[-1] + exposure_matrix <- get_exposure_matrix(serodata) + Nobs <- nrow(serodata) stan_data <- list( Nobs = Nobs, - Npos = seroprev_data$counts, - Ntotal = seroprev_data$total, - Age = seroprev_data$age_mean_f, + Npos = serodata$counts, + Ntotal = serodata$total, + Age = serodata$age_mean_f, Ymax = max(exposure_ages), AgeExpoMatrix = exposure_matrix, NDecades = decades @@ -168,7 +180,7 @@ fit_seroprev_model <- function(seroprev_data, n_warmup <- floor(n_iters / 2) - if (seroprev_model_name == "continuous_foi_normal_log") { + if (seromodel_name == "continuous_foi_normal_log") { f_init <- function() { list(log_foi = rep(-3, length(exposure_ages))) } @@ -227,36 +239,36 @@ fit_seroprev_model <- function(seroprev_data, colnames(foi_post_s) <- exposure_years } - model_object <- list( + seromodel_object <- list( fit = fit, - seroprev_data = seroprev_data, + serodata = serodata, stan_data = stan_data, exposure_years = exposure_years, exposure_ages = exposure_ages, n_iters = n_iters, n_thin = n_thin, n_warmup = n_warmup, - seroprev_model_name = seroprev_model_name, + seromodel_name = seromodel_name, delta = delta, m_treed = m_treed, loo_fit = loo_fit, foi_cent_est = foi_cent_est, foi_post_s = foi_post_s ) - model_object$model_summary <- - extract_seroprev_model_summary(model_object) + seromodel_object$model_summary <- + extract_seromodel_summary(seromodel_object) } else { loo_fit <- c(-1e10, 0) - model_object <- list( + seromodel_object <- list( fit = "no model", - seroprev_data = seroprev_data, + serodata = serodata, stan_data = stan_data, exposure_years = exposure_years, exposure_ages = exposure_ages, n_iters = n_iters, n_thin = n_thin, n_warmup = n_warmup, - model = seroprev_model_name, + model = seromodel_name, delta = delta, m_treed = m_treed, loo_fit = loo_fit, @@ -264,41 +276,43 @@ fit_seroprev_model <- function(seroprev_data, ) } - return(model_object) + return(seromodel_object) } -#' Get exposure years +#' Function that generates an atomic vector containing the corresponding exposition years of a serological survey. #' -#' Function that generates an atomic vector with the exposition years in seroprev_data. The exposition years to the disease for each individual corresponds to the time from birth to the moment of the survey. -#' @param seroprev_data A data frame containing the data from a seroprevalence survey. This data frame must contain the year of birth for each individual (birth_year) and the time of the survey (tsur). birth_year can be constructed by means of the \link{prepare_seroprev_data} function. -#' @return \code{exposure_ages}. An atomic vector with the numeration of the exposition years in seroprev_data +#' This function generates an atomic vector containing the exposition years corresponding to the specified serological survey data \code{serodata}. +#' The exposition years to the disease for each individual corresponds to the time from birth to the moment of the survey. +#' @param serodata A data frame containing the data from a seroprevalence survey. This data frame must contain the year of birth for each individual (birth_year) and the time of the survey (tsur). birth_year can be constructed by means of the \link{prepare_serodata} function. +#' @return \code{exposure_ages}. An atomic vector with the numeration of the exposition years in serodata #' @examples #' \dontrun{ -#' seroprev_data <- prepare_seroprev_data(seroprev_data = serodata, alpha = 0.05) -#' exposure_ages <- get_exposure_ages(seroprev_data) +#' data("serodata") +#' serodata <- prepare_serodata(serodata = serodata, alpha = 0.05) +#' exposure_ages <- get_exposure_ages(serodata) #' } #' @export -get_exposure_ages <- function(seroprev_data) { - # TODO Verify if this change is correct - return(seq_along(min(seroprev_data$birth_year):(seroprev_data$tsur[1] - 1))) +get_exposure_ages <- function(serodata) { + return(seq_along(min(serodata$birth_year):(serodata$tsur[1] - 1))) } - -#' Get Exposure Matrix +# TODO Is necessary to explain better what we mean by the exposure matrix. +#' Function that generates the exposure matrix corresponding to a serological survey. #' -#' Function that generates the exposure matrix for a seroprevalence survey. -#' @param seroprev_data A data frame containing the data from a seroprevalence survey. This data frame must contain the year of birth for each individual (birth_year) and the time of the survey (tsur). birth_year can be constructed by means of the \link{prepare_seroprev_data} function. -#' @return \code{exposure_output}. An atomic matrix containing the expositions for each entry of \code{seroprev_data} by year. +#' Function that generates the exposure matrix corresponding to the specified serological survey data \code{serodata}. +#' @param serodata A data frame containing the data from a seroprevalence survey. This data frame must contain the year of birth for each individual (birth_year) and the time of the survey (tsur). birth_year can be constructed by means of the \link{prepare_serodata} function. +#' @return \code{exposure_output}. An atomic matrix containing the expositions for each entry of \code{serodata} by year. #' @examples #' \dontrun{ -#' seroprev_data <- prepare_seroprev_data(seroprev_data = serodata, alpha = 0.05) -#' exposure_matrix <- get_exposure_matrix(seroprev_data = seroprev_data) +#' data("serodata") +#' serodata <- prepare_serodata(serodata = serodata) +#' exposure_matrix <- get_exposure_matrix(serodata = serodata) #' } #' @export -get_exposure_matrix <- function(seroprev_data) { - age_class <- seroprev_data$age_mean_f - exposure_ages <- get_exposure_ages(seroprev_data) +get_exposure_matrix <- function(serodata) { + age_class <- serodata$age_mean_f + exposure_ages <- get_exposure_ages(serodata) ly <- length(exposure_ages) exposure <- matrix(0, nrow = length(age_class), ncol = ly) for (k in 1:length(age_class)) @@ -308,13 +322,17 @@ get_exposure_matrix <- function(seroprev_data) { } -#' Extract Model Summary +#' Method to extact a summary of the specified serological model object. #' -#' Function to generate a summary of a model. -#' @param model_object \code{model_object}. An object containing relevant information about the implementation of the model. Refer to \link{fit_seroprev_model} for further details. -#' @return \code{model_summary}. Object with a summary of \code{model_object} containing the following: +#' This method extracts a summary corresponding to a serological model object that contains information about the original serological +#' survey data used to fit the model, such as the year when the survey took place, the type of test taken and the corresponding antibody, +#' as well as information about the convergence of the model, like the expected log pointwise predictive density \code{elpd} and its +#' corresponding standar deviation. +#' @param seromodel_object \code{seromodel_object}. An object containing relevant information about the implementation of the model. +#' Refer to \link{fit_seromodel} for further details. +#' @return \code{model_summary}. Object with a summary of \code{seromodel_object} containing the following: #' \tabular{ll}{ -#' \code{seroprev_model_name} \tab Name of the selected model. For further details refer to \link{save_or_load_model}. \cr \tab \cr +#' \code{seromodel_name} \tab Name of the selected model. For further details refer to \link{save_or_load_model}. \cr \tab \cr #' \code{data_set} \tab Seroprevalence survey label.\cr \tab \cr #' \code{country} \tab Name of the country were the survey was conducted in. \cr \tab \cr #' \code{year} \tab Year in which the survey was conducted. \cr \tab \cr @@ -329,39 +347,40 @@ get_exposure_matrix <- function(seroprev_data) { #' } #' @examples #' \dontrun{ -#' seroprev_data <- prepare_seroprev_data(serodata) -#' model_object <- run_seroprev_model(seroprev_data = seroprev_data, -#' seroprev_model_name = "constant_foi_bi") -#' extract_seroprev_model_summary (model_object) +#' data("serodata") +#' serodata <- prepare_serodata(serodata) +#' seromodel_object <- run_seromodel(serodata = serodata, +#' seromodel_name = "constant_foi_bi") +#' extract_seromodel_summary(seromodel_object) #' } #' @export -extract_seroprev_model_summary <- function(model_object) { - seroprev_model_name <- model_object$seroprev_model_name - seroprev_data <- model_object$seroprev_data +extract_seromodel_summary <- function(seromodel_object) { + seromodel_name <- seromodel_object$seromodel_name + serodata <- seromodel_object$serodata #------- Loo estimates - loo_fit <- model_object$loo_fit + loo_fit <- seromodel_object$loo_fit if (sum(is.na(loo_fit)) < 1) { lll <- as.numeric((round(loo_fit$estimates[1, ], 2))) } else { lll <- c(-1e10, 0) } model_summary <- data.frame( - seroprev_model_name = seroprev_model_name, - dataset = seroprev_data$survey[1], - country = seroprev_data$country[1], - year = seroprev_data$tsur[1], - test = seroprev_data$test[1], - antibody = seroprev_data$antibody[1], - n_sample = sum(seroprev_data$total), - n_agec = length(seroprev_data$age_mean_f), - n_iter = model_object$n_iters, + seromodel_name = seromodel_name, + dataset = serodata$survey[1], + country = serodata$country[1], + year = serodata$tsur[1], + test = serodata$test[1], + antibody = serodata$antibody[1], + n_sample = sum(serodata$total), + n_agec = length(serodata$age_mean_f), + n_iter = seromodel_object$n_iters, elpd = lll[1], se = lll[2], converged = NA ) - rhats <- get_table_rhats(model_object) + rhats <- get_table_rhats(seromodel_object) if (any(rhats$rhat > 1.1) == FALSE) { model_summary$converged <- "Yes" } @@ -370,23 +389,25 @@ extract_seroprev_model_summary <- function(model_object) { } -#' Get Prevalence Expanded +#' Auxiliary function that generates an object containing the confidence interval based on a +#' Force-of-Infection fitting performed for a specified serological survey data. #' -#' Function that generates the expanded prevalence -#' @param seroprev_data A data frame containing the data from a seroprevalence survey. For further details refer to \link{run_seroprev_model}. -#' @param foi Object containing the information of the force of infection. It is obtained from \code{rstan::extract(model_object$fit, "foi", inc_warmup = FALSE)[[1]]}. +#' This function computes the corresponding binomial confidence intervals for the obtained prevalence based on a fitting +#' of the Force-of-Infection \code{foi} for plotting an analysis purposes. +#' @param serodata A data frame containing the data from a seroprevalence survey. For further details refer to \link{run_seromodel}. +#' @param foi Object containing the information of the force of infection. It is obtained from \code{rstan::extract(seromodel_object$fit, "foi", inc_warmup = FALSE)[[1]]}. #' @return \code{prev_final}. The expanded prevalence data. This is used for plotting purposes in the \code{visualization} module. #' @examples #' \dontrun{ -#' seroprev_data <- prepare_seroprev_data(serodata) -#' model_object <- run_seroprev_model(seroprev_data = seroprev_data, -#' seroprev_model_name = "constant_foi_bi") -#' foi <- rstan::extract(model_object$fit, "foi", inc_warmup = FALSE)[[1]] -#' get_prev_expanded <- function(foi, seroprev_data) +#' serodata <- prepare_serodata(serodata) +#' seromodel_object <- run_seromodel(serodata = serodata, +#' seromodel_name = "constant_foi_bi") +#' foi <- rstan::extract(seromodel_object$fit, "foi")[[1]] +#' get_prev_expanded <- function(foi, serodata) #' } #' @export get_prev_expanded <- function(foi, - seroprev_data) { + serodata) { dim_foi <- dim(foi)[2] if (dim_foi < 80) { oldest_year <- 80 - dim_foi + 1 @@ -426,7 +447,7 @@ get_prev_expanded <- function(foi, predicted_prev_upper = upper ) - observed_prev <- seroprev_data %>% + observed_prev <- serodata %>% dplyr::select(age_mean_f, prev_obs, prev_obs_lower, @@ -441,11 +462,11 @@ get_prev_expanded <- function(foi, base::merge(predicted_prev, observed_prev, by = "age", - all.x = TRUE) %>% dplyr::mutate(survey = seroprev_data$survey[1]) + all.x = TRUE) %>% dplyr::mutate(survey = serodata$survey[1]) # I added this here for those cases when binned is prefered for plotting - if (seroprev_data$age_max[1] - seroprev_data$age_min[1] < 3) { - xx <- prepare_bin_data(seroprev_data) + if (serodata$age_max[1] - serodata$age_min[1] < 3) { + xx <- prepare_bin_data(serodata) prev_final <- base::merge(prev_expanded, xx, by = "age", all.x = TRUE) } else { diff --git a/R/seroprevalence_data.R b/R/seroprevalence_data.R index e3299a4b..55327f84 100644 --- a/R/seroprevalence_data.R +++ b/R/seroprevalence_data.R @@ -1,7 +1,7 @@ -#' Prepare data +#' Function that prepares the data from a serological survey for modelling. #' -#' Function that prepares the data for modelling -#' @param seroprev_data A data frame containing the data from a seroprevalence survey. +#' This function adds the necessary additional variables to the given dataset \code{serodata} corresponding to a serological survey. +#' @param serodata A data frame containing the data from a serological survey. #' This data frame must contain the following columns: #' \tabular{ll}{ #' \code{survey} \tab survey Label of the current survey \cr \tab \cr @@ -17,7 +17,7 @@ #' \code{antibody} \tab antibody \cr \tab \cr #' } #' @param alpha probability of a type I error. For further details refer to \link{Hmisc::binconf}. -#' @return seroprev_data with additional columns necessary for the analysis. These columns are: +#' @return serodata with additional columns necessary for the analysis. These columns are: #' \tabular{ll}{ #' \code{age_mean_f} \tab Floor value of the average between age_min and age_max \cr \tab \cr #' \code{sample_size} \tab The size of the sample \cr \tab \cr @@ -28,23 +28,23 @@ #' } #' @examples #'\dontrun{ -#' data_test <- readRDS("data/data.RDS") -#' data_test <- prepare_seroprev_data(seroprev_data, alpha) +#' data("serodata") +#' data_test <- prepare_serodata(serodata) #' } #' @export -prepare_seroprev_data <- function(seroprev_data = serodata, - alpha = 0.05, - add_age_mean_f = TRUE) { +prepare_serodata <- function(serodata = serodata, + alpha = 0.05, + add_age_mean_f = TRUE) { if(add_age_mean_f){ - seroprev_data <- seroprev_data %>% + serodata <- serodata %>% dplyr::mutate(age_mean_f = floor((age_min + age_max) / 2), sample_size = sum(total)) %>% dplyr::mutate(birth_year = .data$tsur - .data$age_mean_f) } - seroprev_data <- seroprev_data %>% + serodata <- serodata %>% cbind( Hmisc::binconf( - seroprev_data$counts, - seroprev_data$total, + serodata$counts, + serodata$total, alpha = alpha, method = "exact", return.df = TRUE @@ -57,14 +57,16 @@ prepare_seroprev_data <- function(seroprev_data = serodata, ) %>% dplyr::arrange(.data$age_mean_f) - return(seroprev_data) + return(serodata) } -#' Prepare data to plot binomial confidence intervals +#' Function that prepares a pre-processed serological survey dataset to plot the binomial confidence intervals of the seroprevalence grouped by +#' age group. #' -#' Function that prepares the data for modelling -#' @param seroprev_data A data frame containing the data from a seroprevalence survey. For more information see the function \link{run_seroprev_model}. +#' This function prepapares a given pre-processed serological dataset (see \code{\link{prepare_serodata}}) to plot the binomial confidence intervals +#' of its corresponding seroprevalence grouped by age group. +#' @param serodata A data frame containing the data from a seroprevalence survey. For more information see the function \link{run_seromodel}. #' This data frame must contain the following columns: #' \tabular{ll}{ #' \code{survey} \tab survey Label of the current survey \cr \tab \cr @@ -85,19 +87,19 @@ prepare_seroprev_data <- function(seroprev_data = serodata, #' \code{prev_obs_lower} \tab Lower limit of the confidence interval for the observed prevalence \cr \tab \cr #' \code{prev_obs_upper} \tab Upper limit of the confidence interval for the observed prevalence \cr \tab \cr #' } -#' The last six colums can be added to \code{seroprev_data} by means of the function \code{\link{prepare_seroprev_data}}. +#' The last six colums can be added to \code{serodata} by means of the function \code{\link{prepare_serodata}}. #' @return data set with the binomial confidence intervals #' @examples #'\dontrun{ -#' prepare_bin_data (seroprev_data) +#' prepare_bin_data (serodata) #' } #' @export -prepare_bin_data <- function(seroprev_data) { - seroprev_data$cut_ages <- - cut(as.numeric(seroprev_data$age_mean_f), +prepare_bin_data <- function(serodata) { + serodata$cut_ages <- + cut(as.numeric(serodata$age_mean_f), seq(1, 101, by = 5), include.lowest = TRUE) - xx <- seroprev_data %>% + xx <- serodata %>% dplyr::group_by(.data$cut_ages) %>% dplyr::summarise(bin_size = sum(.data$total), bin_pos = sum(.data$counts)) @@ -118,132 +120,4 @@ prepare_bin_data <- function(seroprev_data) { p_obs_bin_u = .data$Upper ) return(xx) -} - -# TODO: Complete the documentation of get_sim_counts -#' Function that randomly generates a sample of counts for a simulated dataset -#' -#' @param sim_data A dataframe object containing the following columns: -#' \tabular{ll}{ -#' \code{birth_year} \tab List of years in which the subjects were borned \cr \tab \cr -#' \code{tsur} \tab Year of the survey\cr \tab \cr -#' \code{country} \tab Default to 'none'.\cr \tab \cr -#' \code{survey} \tab Survey label \cr \tab \cr -#' \code{age_mean_f} \tab Age \cr \tab \cr -#' } -#' @return A simulated list of counts following a binomial distribution in accordance with a given force of infection and age class sizes. -#' @examples -#'\dontrun{ -#' -#' } -#' @export -get_sim_counts <- function(sim_data, foi, size_age_class, seed = 1234){ - exposure_ages <- get_exposure_ages(sim_data) - exposure_matrix <- get_exposure_matrix(sim_data) - - set.seed(seed = seed) - sim_probabilities <- purrr::map_dbl(exposure_ages, ~1-exp(-dot(exposure_matrix[., ], foi))) - sim_counts <- purrr::map_int(sim_probabilities, ~rbinom(1, size_age_class, .)) - - return(sim_counts) -} - -# TODO: Complete the documentation of generate_sim_data -#' Function that generates simulated data from a given Force-of-Infection -#' -#' @param foi Numeric atomic vector corresponding to the desired Force-of-Infection -#' @return Dataframe object containing the simulated data generated from \code{foi} -#' @examples -#'\dontrun{ -#' size_age_class = 5 -#' foi <- rep(0.02, 50) -#' sim_data <- generate_sim_data(foi = foi, -#' size_age_class = size_age_class, -#' tsur = 2050, -#' birth_year_min = 2000, -#' survey_label = 'sim_constant_foi') -#' } -#' @export -generate_sim_data <- function(foi, - size_age_class, - tsur, - birth_year_min, - survey_label, - test = "fake", - antibody = "IgG", - seed = 1234 - ){ - sim_data <- data.frame(birth_year = c(birth_year_min:(tsur - 1))) %>% - mutate(tsur = tsur, - country = 'None', - test = test, - antibody = antibody, - survey = survey_label, - age_mean_f = tsur - birth_year) - sim_data <- sim_data %>% - mutate(counts = get_sim_counts(sim_data, foi, size_age_class, seed = seed), - total = size_age_class) %>% - prepare_seroprev_data(add_age_mean_f = FALSE) - - return(sim_data) -} - -# TODO: Complete the documentation of group_sim_data -#' Function that generates grouped simulated data from a given Force-of-Infection -#' -#' @param sim_data Dataframe with the structure of the output of \code{\linl{generate_sim_data}}. -#' @return Dataframe object containing grouped simulated data generated from \code{foi} -#' @examples -#'\dontrun{ -#' size_age_class = 5 -#' foi <- rep(0.02, 50) -#' sim_data <- generate_sim_data(foi = foi, -#' size_age_class = size_age_class, -#' tsur = 2050, -#' birth_year_min = 2000, -#' survey_label = 'sim_constant_foi') -#' sim_data_grouped <- group_sim_data(sim_data = sim_data, -#' foi = foi, -#' size_age_class = size_age_class, -#' tsur = 2050, -#' birth_year_min = 2000, -#' survey_label = 'sim_constant_foi_grouped') -#' } -#' @export -group_sim_data <- function(sim_data, - foi, - size_age_class, - tsur, - birth_year_min, - survey_label, - test = "fake", - antibody = "IgG", - seed = 1234) { - - sim_data <- sim_data %>% mutate(age_group = 'NA', age = age_mean_f) %>% arrange(age) - sim_data$age_group[sim_data$age > 0 & sim_data$age < 5] <- '01-04' - sim_data$age_group[sim_data$age > 4 & sim_data$age < 10] <- '05-09' - sim_data$age_group[sim_data$age > 9 & sim_data$age < 15] <- '10-14' - sim_data$age_group[sim_data$age > 14 & sim_data$age < 20] <- '15-19' - sim_data$age_group[sim_data$age > 19 & sim_data$age < 25] <- '20-24' - sim_data$age_group[sim_data$age > 24 & sim_data$age < 30] <- '25-29' - sim_data$age_group[sim_data$age > 29 & sim_data$age < 35] <- '30-34' - sim_data$age_group[sim_data$age > 34 & sim_data$age < 40] <- '35-39' - sim_data$age_group[sim_data$age > 39 & sim_data$age < 45] <- '40-44' - sim_data$age_group[sim_data$age > 44 & sim_data$age < 51] <- '45-50' - - - sim_data_grouped <- sim_data %>% group_by(age_group) %>% - dplyr::summarise(total = sum(total), counts = sum(counts)) %>% - mutate(tsur = sim_data$tsur[1], - country = "None", - survey = survey_label, - test = test, - antibody = antibody) %>% - mutate(age_min = as.numeric(substr(age_group, 1, 2)), - age_max = as.numeric(substr(age_group, 4, 5))) %>% - prepare_seroprev_data() - - return(sim_data_grouped) - -} +} \ No newline at end of file diff --git a/R/test_vignettes.R b/R/test_vignettes.R deleted file mode 100644 index a8deac9e..00000000 --- a/R/test_vignettes.R +++ /dev/null @@ -1,4 +0,0 @@ -#library(remotes) - -#remotes::install_github("TRACE-LAC/serofoi", ref = "dev") -#library(serofoi) diff --git a/R/visualization.R b/R/visualization.R index ffac5a0c..d94aea9b 100644 --- a/R/visualization.R +++ b/R/visualization.R @@ -1,7 +1,6 @@ -#' Generate sero-positivity plot from raw data +#' Function that generates the sero-positivity plot from a raw serological survey dataset. #' -#' Function that generates the sero positivity plot from raw data -#' @param seroprev_data A data frame containing the data from a seroprevalence survey. +#' @param serodata A data frame containing the data from a seroprevalence survey. #' This data frame must contain the following columns: #' \tabular{ll}{ #' \code{survey} \tab survey Label of the current survey \cr \tab \cr @@ -20,18 +19,18 @@ #' @return A ggplot object containing the seropositivity-vs-age graph of the raw data of a given seroprevalence survey with its corresponging binomial confidence interval. #' @examples #' \dontrun{ -#' data_test <- prepare_seroprev_data(serodata) -#' model_object <- run_seroprev_model( -#' seroprev_data = data_test, -#' seroprev_model_name = "constant_foi_bi", +#' data_test <- prepare_serodata(serodata) +#' seromodel_object <- run_seromodel( +#' serodata = data_test, +#' seromodel_name = "constant_foi_bi", #' n_iters = 1000 #') -#' plot_seroprev(model_object, size_text = 15) +#' plot_seroprev(seromodel_object, size_text = 15) #' } #' @export -plot_seroprev <- function(seroprev_data, +plot_seroprev <- function(serodata, size_text = 6) { - xx <- prepare_bin_data(seroprev_data) + xx <- prepare_bin_data(serodata) seroprev_plot <- ggplot2::ggplot(data = xx) + ggplot2::geom_errorbar(ggplot2::aes(age, ymin = p_obs_bin_l, ymax = p_obs_bin_u), width = 0.1) + @@ -45,31 +44,32 @@ plot_seroprev <- function(seroprev_data, } -#' Generate sero-positivity plot with fitted model +#' Function that generates a seropositivity plot corresponding to the specified fitted serological model. #' -#' Function that generates the seropositivity graph with fitted model. Age is located on the x axis and seropositivity on the y axis with its corresponding confidence interval. -#' @param model_object Object containing the results of fitting a model by means of \link{run_seroprev_model}. +#' This function generates a seropositivity plot of the specified serological model object. This includes the original data grouped by age +#' as well as the obtained fitting from the model implementation. Age is located on the x axis and seropositivity on the y axis with its +#' corresponding confidence interval. +#' @param seromodel_object Object containing the results of fitting a model by means of \link{run_seromodel}. #' @param size_text Text size of the graph returned by the function. #' @return A ggplot object containing the seropositivity-vs-age graph including the data, the fitted model and their corresponding confindence intervals. #' @examples #' \dontrun{ -#' data_test <- prepare_seroprev_data(serodata) -#' model_object <- run_seroprev_model( -#' seroprev_data = data_test, -#' seroprev_model_name = "constant_foi_bi", -#' n_iters = 1000 -#') -#' plot_seroprev_fitted(model_object, size_text = 15) +#' data("serodata") +#' data_test <- prepare_serodata(serodata) +#' seromodel_object <- run_seromodel(serodata = data_test, +#' seromodel_name = "constant_foi_bi", +#' n_iters = 1000) +#' plot_seroprev_fitted(seromodel_object, size_text = 15) #' } #' @export -plot_seroprev_fitted <- function(model_object, +plot_seroprev_fitted <- function(seromodel_object, size_text = 6) { - if (is.character(model_object$fit) == FALSE) { - if (class(model_object$fit@sim$samples) != "NULL" ) { + if (is.character(seromodel_object$fit) == FALSE) { + if (class(seromodel_object$fit@sim$samples) != "NULL" ) { - foi <- rstan::extract(model_object$fit, "foi", inc_warmup = FALSE)[[1]] - prev_expanded <- get_prev_expanded(foi, seroprev_data = model_object$seroprev_data) + foi <- rstan::extract(seromodel_object$fit, "foi", inc_warmup = FALSE)[[1]] + prev_expanded <- get_prev_expanded(foi, serodata = seromodel_object$serodata) prev_plot <- ggplot2::ggplot(prev_expanded) + ggplot2::geom_ribbon( @@ -121,35 +121,37 @@ plot_seroprev_fitted <- function(model_object, return(prev_plot) } -#' Generate Force-of-Infection Plot +#' Function that generates a Force-of-Infection plot corresponding to the specified fitted serological model. #' -#' Function that generates the Force-of-Infection plot. The x axis corresponds to the decades covered by the survey the y axis to the Force-of-Infection. -#' @param model_object Object containing the results of fitting a model by means of \link{run_seroprev_model}. +#' This function generates a Force-of-Infection plot from the results obtained by fitting a serological model. +#' This includes the corresponding binomial confidence interval. +#' The x axis corresponds to the decades covered by the survey the y axis to the Force-of-Infection. +#' @param seromodel_object Object containing the results of fitting a model by means of \link{run_seromodel}. #' @param size_text Text size use in the theme of the graph returned by the function. #' @return A ggplot2 object containing the Force-of-infection-vs-time including the corresponding confidence interval. #' @examples #' \dontrun{ -#' data_test <- prepare_seroprev_data(serodata) -#' model_object <- run_seroprev_model( -#' seroprev_data = data_test, -#' seroprev_model_name = "constant_foi_bi", +#' data_test <- prepare_serodata(serodata) +#' seromodel_object <- run_seromodel( +#' serodata = data_test, +#' seromodel_name = "constant_foi_bi", #' n_iters = 1000 #' ) -#' plot_foi(model_object, size_text = 15) +#' plot_foi(seromodel_object, size_text = 15) #' } #' @export -plot_foi <- function(model_object, +plot_foi <- function(seromodel_object, lambda_sim = NA, max_lambda = NA, size_text = 25) { - if (is.character(model_object$fit) == FALSE) { - if (class(model_object$fit@sim$samples) != "NULL") { - foi <- rstan::extract(model_object$fit, + if (is.character(seromodel_object$fit) == FALSE) { + if (class(seromodel_object$fit@sim$samples) != "NULL") { + foi <- rstan::extract(seromodel_object$fit, "foi", inc_warmup = FALSE)[[1]] #-------- This bit is to get the actual length of the foi data - foi_data <- model_object$foi_cent_est + foi_data <- seromodel_object$foi_cent_est if (!is.na(lambda_sim)) { lambda_mod_length <- NROW(foi_data) @@ -213,29 +215,32 @@ plot_foi <- function(model_object, return(foi_plot) } -#' Generate Rhats-Convergence Plot +#' Function that generates a plot of the R-hat estimates of the specified fitted serological model. #' -#' Function that generates the convergence graph of a model. The x axis corresponds to the decades covered by the survey and the y axis to the value of the rhats. -#' All rhats must be smaller than 1 to ensure convergence. -#' @param model_object Object containing the results of fitting a model by means of \link{run_seroprev_model}. +#' This function generates a plot of the R-hat estimates obtained for a specified fitted serological model \code{seromodel_object}. +#' The x axis corresponds to the decades covered by the survey and the y axis to the value of the rhats. +#' All rhats must be smaller than 1 to ensure convergence (for further details check \link[bayesplot]{rhat}). +#' @param seromodel_object Object containing the results of fitting a model by means of \link{run_seromodel}. #' @param size_text Text size use in the theme of the graph returned by the function. #' @return The rhats-convergence plot of the selected model. #' @examples #' \dontrun{ -#' data_test <- prepare_seroprev_data(serodata) -#' model_object <- run_seroprev_model( -#' seroprev_data = data_test, -#' seroprev_model_name = "constant_foi_bi", +#' data("serodata") +#' data_test <- prepare_serodata(serodata) +#' seromodel_object <- run_seromodel( +#' serodata = data_test, +#' seromodel_name = "constant_foi_bi", #' n_iters = 1000 #') -#' plot_rhats(model_object, size_text = 15) +#' plot_rhats(seromodel_object, +#' size_text = 15) #' } #' @export -plot_rhats <- function(model_object, +plot_rhats <- function(seromodel_object, size_text = 25) { - if (is.character(model_object$fit) == FALSE) { - if (class(model_object$fit@sim$samples) != "NULL") { - rhats <- get_table_rhats(model_object) + if (is.character(seromodel_object$fit) == FALSE) { + if (class(seromodel_object$fit@sim$samples) != "NULL") { + rhats <- get_table_rhats(seromodel_object) rhats_plot <- ggplot2::ggplot(rhats, ggplot2::aes(year, rhat)) + @@ -275,45 +280,45 @@ plot_rhats <- function(model_object, } -#' Generate a vertical arrange of plots summarizing the results of the model implementation +#' Function that generates a vertical arrange of plots showing a summary of a given implemented model, as well as its corresponding +#' seroprevalence, Force-of-Infection and R-hat estimates plots. #' -#' Function that generates the combined plots summarizing the results of the model implementation -#' @param model_object Object containing the results of fitting a model by means of \link{run_seroprev_model}. +#' @param seromodel_object Object containing the results of fitting a model by means of \link{run_seromodel}. #' @param size_text Text size use in the theme of the graph returned by the function. #' @return A ggplot object with a vertical arrange containing the seropositivity, force of infection, and convergence plots. #' @examples #' \dontrun{ -#' data_test <- prepare_seroprev_data(serodata) -#' model_object <- run_seroprev_model( -#' seroprev_data = data_test, -#' seroprev_model_name = "constant_foi_bi", +#' data_test <- prepare_serodata(serodata) +#' seromodel_object <- run_seromodel( +#' serodata = data_test, +#' seromodel_name = "constant_foi_bi", #' n_iters = 1000 #') -#' plot_seroprev_model(model_object, size_text = 15) +#' plot_seromodel(seromodel_object, size_text = 15) #' } #' @export -plot_seroprev_model <- function(model_object, +plot_seromodel <- function(seromodel_object, lambda_sim = NA, max_lambda = NA, size_text = 25) { - if (is.character(model_object$fit) == FALSE) { - if (class(model_object$fit@sim$samples) != "NULL") { - prev_plot <- plot_seroprev_fitted(model_object = model_object, + if (is.character(seromodel_object$fit) == FALSE) { + if (class(seromodel_object$fit@sim$samples) != "NULL") { + prev_plot <- plot_seroprev_fitted(seromodel_object = seromodel_object, size_text = size_text) foi_plot <- plot_foi( - model_object = model_object, + seromodel_object = seromodel_object, lambda_sim = lambda_sim, max_lambda = max_lambda, size_text = size_text ) - rhats_plot <- plot_rhats(model_object = model_object, + rhats_plot <- plot_rhats(seromodel_object = seromodel_object, size_text = size_text) summary_table <- t( - dplyr::select(model_object$model_summary, - c('seroprev_model_name', 'dataset', 'elpd', 'se', 'converged'))) + dplyr::select(seromodel_object$model_summary, + c('seromodel_name', 'dataset', 'elpd', 'se', 'converged'))) summary_plot <- plot_info_table(summary_table, size_text = size_text) @@ -348,7 +353,7 @@ plot_seroprev_model <- function(model_object, ggplot2::ylab(" ") + ggplot2::xlab(" ") g1 <- g0 - g0 <- g0 + ggplot2::labs(subtitle = model_object$model) + + g0 <- g0 + ggplot2::labs(subtitle = seromodel_object$model) + ggplot2::theme(plot.title = ggplot2::element_text(size = 10)) plot_arrange <- @@ -359,21 +364,20 @@ plot_seroprev_model <- function(model_object, } -#' Plot Info Table +#' Auxiliary function that generates a plot of a given table #' -#' Function that generates the information table #' @param info the information that will be contained in the table #' @param size_text Text size of the graph returned by the function #' @return p, a variable that will be used in the \link{visualisation} module #' @examples #' \dontrun{ -#' data_test <- prepare_seroprev_data(serodata) -#' model_object <- run_seroprev_model( -#' seroprev_data = data_test, -#' seroprev_model_name = "constant_foi_bi", +#' data_test <- prepare_serodata(serodata) +#' seromodel_object <- run_seromodel( +#' serodata = data_test, +#' seromodel_name = "constant_foi_bi", #' n_iters = 1000 #') -#' info = t(model_object$model_summary) +#' info = t(seromodel_object$model_summary) #' plot_info_table (info, size_text = 15) #' } #' @export @@ -388,30 +392,4 @@ plot_info_table <- function(info, size_text) { fontface = "bold") return(p) -} - - -#' Function that generates a plot showing the results for several models. -#' -#' @param ... The first n arguments of the function are taken as plots ti be arranged in a single plot. -#' This objects can be obtained by means of \link{run_seroprev_model}. -#' @param n_row Number of rows of the plot array. -#' @param n_col Number of columns of the plot array. -#' @return A ggplot object containing an array with the plots in \code{models_list}. -#' @examples -#' \dontrun{ -#' data_test <- prepare_seroprev_data(serodata) -#' model_object_constant <- run_seroprev_model(seroprev_data = data_test, -#' seroprev_model_name = "constant_foi_bi") -#' model_object_normalbi <- run_seroprev_model(seroprev_data = data_test, -#' seroprev_model_name = "continuous_foi_normal_bi") -#' model_object_normallog <- run_seroprev_model(seroprev_data = data_test, -#' seroprev_model_name = "continuous_foi_normal_log") -#' plot_seroprev_models_list(model_object_constant, model_object_normalbi, model_object_normallog, n_row = 1, n_col = 3) -#' } -#' @export -plot_seroprev_models_grid <- function(..., n_row = NULL, n_col = NULL) { - # TODO Distribute rows and columns according to the number of rows. - plot_seroprev_models <- cowplot::plot_grid(..., nrow = n_row, ncol = n_col) - return(plot_seroprev_models) } \ No newline at end of file diff --git a/README.Rmd b/README.Rmd index 745c0010..cf53b1af 100644 --- a/README.Rmd +++ b/README.Rmd @@ -72,14 +72,14 @@ head(serodata, 5) ``` -The function `prepare_seroprev_data` will prepare the entry data for entering the +The function `prepare_serodata` will prepare the entry data for entering the modelling functions. The seroprevalence *prepared data* can be visualised with the `plot_seroprev` function. This function also plots the binomial confidence interval of the observed data. ```{r data_test, include = TRUE, fig.height=6, fig.width=6, message=FALSE} -data_test <- prepare_seroprev_data(serodata) +data_test <- prepare_serodata(serodata) plot_seroprev(data_test, size_text = 15) @@ -87,7 +87,7 @@ plot_seroprev(data_test, size_text = 15) #### Current version of the package runs ***three*** different FoI models -The `run_seroprev_model` function allows specifying the Bayesian model from *R*, +The `run_seromodel` function allows specifying the Bayesian model from *R*, while running in the back from `rstan`. The number of iterations, thinning, and other parameters can be customised. @@ -104,8 +104,8 @@ achieving convergence, as it only fits one parameter (the constant FoI) from a binomial distribution. ```{r model_1, include = TRUE, errors = FALSE, warning = FALSE, message = FALSE } -model_1 <- run_seroprev_model(seroprev_data = data_test, - seroprev_model_name = "constant_foi_bi", +model_1 <- run_seromodel(serodata = data_test, + seromodel_name = "constant_foi_bi", n_iters = 500, n_thin = 2) @@ -120,8 +120,8 @@ of years, reflected by the difference between year of the serosurvey and the maximum age-class sampled. ```{r model_2, include = TRUE, errors = FALSE, warning = FALSE, message = FALSE } -model_2 <- run_seroprev_model(seroprev_data = data_test, - seroprev_model_name = "continuous_foi_normal_bi", +model_2 <- run_seromodel(serodata = data_test, + seromodel_name = "continuous_foi_normal_bi", n_iters = 1500, n_thin = 2) ``` @@ -132,34 +132,34 @@ For the *fast* *epidemic model,* a larger number of iterations is required for achieving convergence, compared to the previous models. ```{r model_3, include = TRUE, errors = FALSE, warning = FALSE, message = FALSE } -model_3 <- run_seroprev_model(seroprev_data = data_test, - seroprev_model_name = "continuous_foi_normal_log", +model_3 <- run_seromodel(serodata = data_test, + seromodel_name = "continuous_foi_normal_log", n_iters = 1500, n_thin = 2) ``` -For each model, the plot_seroprev_model function generate a vertical arrange of +For each model, the plot_seromodel function generate a vertical arrange of plots summarising the results of the model implementation plotting functions. Crucially, it shows the (expected) log-predictive density `elpd`, standard error `se`, and allows to check convergence based on `R-hat` convergence diagnostics. -```{r plot_seroprev_model, include = FALSE, echo = TRUE, errors = FALSE, warning = FALSE, message = FALSE, fig.height=12, fig.width=6} +```{r plot_seromodel, include = FALSE, echo = TRUE, errors = FALSE, warning = FALSE, message = FALSE, fig.height=12, fig.width=6} -plot_seroprev_model(model_2, size_text = 15) +plot_seromodel(model_2, size_text = 15) ``` -Also, the `plot_seroprev_models_grid` allows a visual a comparison of +Also, the `plot_seromodels_grid` allows a visual a comparison of the models based on the (expected) log-predictive density `elpd`, standard error `se`, and allows to check convergence based on `R-hat` convergence diagnostics. ```{r, model_comp, include=FALSE, echo=TRUE, warning=FALSE, message=FALSE, fig.height=12, message=FALSE} -p1 <- plot_seroprev_model(model_1, size_text = 15) -p2 <- plot_seroprev_model(model_2, size_text = 15) -p3 <- plot_seroprev_model(model_3, size_text = 15) +p1 <- plot_seromodel(model_1, size_text = 15) +p2 <- plot_seromodel(model_2, size_text = 15) +p3 <- plot_seromodel(model_3, size_text = 15) -plot_seroprev_models_grid(p1, p2, p3, n_row = 1, n_col = 3) +plot_seromodels_grid(p1, p2, p3, n_row = 1, n_col = 3) ``` For more detailed information and examples, please check the [online diff --git a/README.md b/README.md index 14212826..9123bd31 100644 --- a/README.md +++ b/README.md @@ -60,14 +60,14 @@ head(serodata, 5) #> 5 IgG anti-T.cruzi ``` -The function `prepare_seroprev_data` will prepare the entry data for +The function `prepare_serodata` will prepare the entry data for entering the modelling functions. The seroprevalence *prepared data* can be visualised with the `plot_seroprev` function. This function also plots the binomial confidence interval of the observed data. ``` r -data_test <- prepare_seroprev_data(serodata) +data_test <- prepare_serodata(serodata) plot_seroprev(data_test, size_text = 15) ``` @@ -76,7 +76,7 @@ plot_seroprev(data_test, size_text = 15) #### Current version of the package runs ***three*** different FoI models -The `run_seroprev_model` function allows specifying the Bayesian model +The `run_seromodel` function allows specifying the Bayesian model from *R*, while running in the back from `rstan`. The number of iterations, thinning, and other parameters can be customised. @@ -95,13 +95,13 @@ achieving convergence, as it only fits one parameter (the constant FoI) from a binomial distribution. ``` r -model_1 <- run_seroprev_model(seroprev_data = data_test, - seroprev_model_name = "constant_foi_bi", +model_1 <- run_seromodel(serodata = data_test, + seromodel_name = "constant_foi_bi", n_iters = 500, n_thin = 2) #> [1] "serofoi model constant_foi_bi finished running ------" #> [,1] -#> seroprev_model_name "constant_foi_bi" +#> seromodel_name "constant_foi_bi" #> dataset "COL-035-93" #> country "COL" #> year "2012" @@ -124,13 +124,13 @@ of years, reflected by the difference between year of the serosurvey and the maximum age-class sampled. ``` r -model_2 <- run_seroprev_model(seroprev_data = data_test, - seroprev_model_name = "continuous_foi_normal_bi", +model_2 <- run_seromodel(serodata = data_test, + seromodel_name = "continuous_foi_normal_bi", n_iters = 1500, n_thin = 2) #> [1] "serofoi model continuous_foi_normal_bi finished running ------" #> [,1] -#> seroprev_model_name "continuous_foi_normal_bi" +#> seromodel_name "continuous_foi_normal_bi" #> dataset "COL-035-93" #> country "COL" #> year "2012" @@ -150,13 +150,13 @@ For the *fast* *epidemic model,* a larger number of iterations is required for achieving convergence, compared to the previous models. ``` r -model_3 <- run_seroprev_model(seroprev_data = data_test, - seroprev_model_name = "continuous_foi_normal_log", +model_3 <- run_seromodel(serodata = data_test, + seromodel_name = "continuous_foi_normal_log", n_iters = 1500, n_thin = 2) #> [1] "serofoi model continuous_foi_normal_log finished running ------" #> [,1] -#> seroprev_model_name "continuous_foi_normal_log" +#> seromodel_name "continuous_foi_normal_log" #> dataset "COL-035-93" #> country "COL" #> year "2012" @@ -176,7 +176,7 @@ plotting functions. Crucially, it shows the (expected) log-predictive density `elpd`, standard error `se`, and allows to check convergence based on `R-hat` convergence diagnostics. -Also, the `plot_seroprev_models_grid` allows a visual a comparison of +Also, the `plot_seromodels_grid` allows a visual a comparison of the models based on the (expected) log-predictive density `elpd`, standard error `se`, and allows to check convergence based on `R-hat` convergence diagnostics. diff --git a/data/data.RData b/data/data.RData deleted file mode 100644 index d6cbc1f0..00000000 Binary files a/data/data.RData and /dev/null differ diff --git a/data/data_test.R b/data/data_test.R deleted file mode 100644 index 0104ef0a..00000000 --- a/data/data_test.R +++ /dev/null @@ -1,15 +0,0 @@ -#' Seroprevalence data on serofoi -#' -#' Data from a serological surveys -#' -#' @docType data -#' -#' @usage data_test -#' -#' @format An object of class \code{"cross"}; see \code{\link[qtl]{read.cross}}. -#' -#' @keywords datasets -#' -#' @examples -#' data_test -"serofoi" diff --git a/data/fake_data_1569526695.RDS b/data/fake_data_1569526695.RDS deleted file mode 100644 index 9d0ad0b3..00000000 Binary files a/data/fake_data_1569526695.RDS and /dev/null differ diff --git a/data/mydata.RData b/data/mydata.RData deleted file mode 100644 index 2e8a51fe..00000000 Binary files a/data/mydata.RData and /dev/null differ diff --git a/man/extract_seromodel_summary.Rd b/man/extract_seromodel_summary.Rd new file mode 100644 index 00000000..b9f4424b --- /dev/null +++ b/man/extract_seromodel_summary.Rd @@ -0,0 +1,44 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/modelling.R +\name{extract_seromodel_summary} +\alias{extract_seromodel_summary} +\title{Method to extact a summary of the specified serological model object.} +\usage{ +extract_seromodel_summary(seromodel_object) +} +\arguments{ +\item{seromodel_object}{\code{seromodel_object}. An object containing relevant information about the implementation of the model. +Refer to \link{fit_seromodel} for further details.} +} +\value{ +\code{model_summary}. Object with a summary of \code{seromodel_object} containing the following: +\tabular{ll}{ +\code{seromodel_name} \tab Name of the selected model. For further details refer to \link{save_or_load_model}. \cr \tab \cr +\code{data_set} \tab Seroprevalence survey label.\cr \tab \cr +\code{country} \tab Name of the country were the survey was conducted in. \cr \tab \cr +\code{year} \tab Year in which the survey was conducted. \cr \tab \cr +\code{test} \tab Type of test of the survey. \cr \tab \cr +\code{antibody} \tab Antibody \cr \tab \cr +\code{n_sample} \tab Total number of samples in the survey. \cr \tab \cr +\code{n_agec} \tab Number of age groups considered. \cr \tab \cr +\code{n_iter} \tab Number of interations for eah chain including the warmup. \cr \tab \cr +\code{elpd} \tab elpd \cr \tab \cr +\code{se} \tab se \cr \tab \cr +\code{converged} \tab convergence \cr \tab \cr +} +} +\description{ +This method extracts a summary corresponding to a serological model object that contains information about the original serological +survey data used to fit the model, such as the year when the survey took place, the type of test taken and the corresponding antibody, +as well as information about the convergence of the model, like the expected log pointwise predictive density \code{elpd} and its +corresponding standar deviation. +} +\examples{ +\dontrun{ +data("serodata") +serodata <- prepare_serodata(serodata) +seromodel_object <- run_seromodel(serodata = serodata, + seromodel_name = "constant_foi_bi") +extract_seromodel_summary(seromodel_object) +} +} diff --git a/man/extract_seroprev_model_summary.Rd b/man/extract_seroprev_model_summary.Rd deleted file mode 100644 index 498b3b1f..00000000 --- a/man/extract_seroprev_model_summary.Rd +++ /dev/null @@ -1,39 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/modelling.R -\name{extract_seroprev_model_summary} -\alias{extract_seroprev_model_summary} -\title{Extract Model Summary} -\usage{ -extract_seroprev_model_summary(model_object) -} -\arguments{ -\item{model_object}{\code{model_object}. An object containing relevant information about the implementation of the model. Refer to \link{fit_seroprev_model} for further details.} -} -\value{ -\code{model_summary}. Object with a summary of \code{model_object} containing the following: -\tabular{ll}{ -\code{seroprev_model_name} \tab Name of the selected model. For further details refer to \link{save_or_load_model}. \cr \tab \cr -\code{data_set} \tab Seroprevalence survey label.\cr \tab \cr -\code{country} \tab Name of the country were the survey was conducted in. \cr \tab \cr -\code{year} \tab Year in which the survey was conducted. \cr \tab \cr -\code{test} \tab Type of test of the survey. \cr \tab \cr -\code{antibody} \tab Antibody \cr \tab \cr -\code{n_sample} \tab Total number of samples in the survey. \cr \tab \cr -\code{n_agec} \tab Number of age groups considered. \cr \tab \cr -\code{n_iter} \tab Number of interations for eah chain including the warmup. \cr \tab \cr -\code{elpd} \tab elpd \cr \tab \cr -\code{se} \tab se \cr \tab \cr -\code{converged} \tab convergence \cr \tab \cr -} -} -\description{ -Function to generate a summary of a model. -} -\examples{ -\dontrun{ -seroprev_data <- prepare_seroprev_data(serodata) -model_object <- run_seroprev_model(seroprev_data = seroprev_data, - seroprev_model_name = "constant_foi_bi") -extract_seroprev_model_summary (model_object) -} -} diff --git a/man/fit_seroprev_model.Rd b/man/fit_seromodel.Rd similarity index 68% rename from man/fit_seroprev_model.Rd rename to man/fit_seromodel.Rd index 34f23b2d..da2d2d0c 100644 --- a/man/fit_seroprev_model.Rd +++ b/man/fit_seromodel.Rd @@ -1,12 +1,12 @@ % Generated by roxygen2: do not edit by hand % Please edit documentation in R/modelling.R -\name{fit_seroprev_model} -\alias{fit_seroprev_model} -\title{Fit Model} +\name{fit_seromodel} +\alias{fit_seromodel} +\title{Function that fits the selected model to the specified seroprevalence survey data.} \usage{ -fit_seroprev_model( - seroprev_data, - seroprev_model_name, +fit_seromodel( + serodata, + seromodel_name, n_iters = 1000, n_thin = 2, delta = 0.9, @@ -15,9 +15,9 @@ fit_seroprev_model( ) } \arguments{ -\item{seroprev_data}{A data frame containing the data from a seroprevalence survey. For further details refer to \link{run_seroprev_model}.} +\item{serodata}{A data frame containing the data from a seroprevalence survey. For further details refer to \link{run_seromodel}.} -\item{seroprev_model_name}{Name of the selected model. Current version provides three options: +\item{seromodel_name}{Name of the selected model. Current version provides three options: \describe{ \item{\code{"constant_foi_bi"}}{Runs a constant model} \item{\code{"continuous_foi_normal_bi"}}{Runs a normal model} @@ -37,10 +37,10 @@ For further details refer to the \code{control} parameter in \link[rstan]{sampli \item{decades}{Number of decades covered by the survey data.} } \value{ -\code{model_object}. An object containing relevant information about the implementation of the model. It contains the following: +\code{seromodel_object}. An object containing relevant information about the implementation of the model. It contains the following: \tabular{ll}{ \code{fit} \tab \code{stanfit} object returned by the function \link[rstan]{sampling} \cr \tab \cr -\code{seroprev_data} \tab A data frame containing the data from a seroprevalence survey. For further details refer to \link{run_seroprev_model}.\cr \tab \cr +\code{serodata} \tab A data frame containing the data from a seroprevalence survey. For further details refer to \link{run_seromodel}.\cr \tab \cr \code{stan_data} \tab List containing \code{Nobs}, \code{Npos}, \code{Ntotal}, \code{Age}, \code{Ymax}, \code{AgeExpoMatrix} and \code{NDecades}. This object is used as an input for the \link[rstan]{sampling} function \cr \tab \cr \code{exposure_years} \tab Integer atomic vector containing the actual exposure years (1946, ..., 2007 e.g.) \cr \tab \cr @@ -48,23 +48,26 @@ This object is used as an input for the \link[rstan]{sampling} function \cr \tab \code{n_iters} \tab Number of interations for eah chain including the warmup. \cr \tab \cr \code{n_thin} \tab Positive integer specifying the period for saving samples. \cr \tab \cr \code{n_warmup} \tab Number of warm up iterations. Set by default as n_iters/2. \cr \tab \cr -\code{seroprev_model_name} \tab The name of the model\cr \tab \cr +\code{seromodel_name} \tab The name of the model\cr \tab \cr \code{delta} \tab Real number between 0 and 1 that represents the target average acceptance probability. \cr \tab \cr \code{m_treed} \tab Maximum tree depth for the binary tree used in the NUTS stan sampler. \cr \tab \cr \code{loo_fit} \tab Efficient approximate leave-one-out cross-validation. Refer to \link[loo]{loo} for further details. \cr \tab \cr \code{foi_cent_est} \tab A data fram e containing \code{year} (corresponding to \code{exposure_years}), \code{lower}, \code{upper}, and \code{medianv} \cr \tab \cr \code{foi_post_s} \tab Sample n rows from a table. Refer to \link[dplyr]{sample_n} for further details. \cr \tab \cr -\code{model_summary} \tab A data fram containing the summary of the model. Refer to \link{extract_seroprev_model_summary} for further details. \cr \tab \cr +\code{model_summary} \tab A data fram containing the summary of the model. Refer to \link{extract_seromodel_summary} for further details. \cr \tab \cr } } \description{ -Function that fits the selected model to the data +This function fits the specified model \code{seromodel_name} to the serological survey data \code{serodata} +by means of the \link[rstan]{sampling} method. The function determines whether the corresponding stan model +object needs to be compiled by means of the function \link{save_or_load_model}. } \examples{ \dontrun{ -seroprev_data <- prepare_seroprev_data(serodata) -fit_seroprev_model (seroprev_data, - seroprev_model_name = "constant_foi_bi") +data("serodata") +serodata <- prepare_serodata(serodata) +seromodel_fit <- fit_seromodel(serodata = serodata, + seromodel_name = "constant_foi_bi") } } diff --git a/man/generate_sim_data.Rd b/man/generate_sim_data.Rd deleted file mode 100644 index 449564c7..00000000 --- a/man/generate_sim_data.Rd +++ /dev/null @@ -1,37 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/seroprevalence_data.R -\name{generate_sim_data} -\alias{generate_sim_data} -\title{Function that generates simulated data from a given Force-of-Infection} -\usage{ -generate_sim_data( - foi, - size_age_class, - tsur, - birth_year_min, - survey_label, - test = "fake", - antibody = "IgG", - seed = 1234 -) -} -\arguments{ -\item{foi}{Numeric atomic vector corresponding to the desired Force-of-Infection} -} -\value{ -Dataframe object containing the simulated data generated from \code{foi} -} -\description{ -Function that generates simulated data from a given Force-of-Infection -} -\examples{ -\dontrun{ -size_age_class = 5 -foi <- rep(0.02, 50) -sim_data <- generate_sim_data(foi = foi, - size_age_class = size_age_class, - tsur = 2050, - birth_year_min = 2000, - survey_label = 'sim_constant_foi') -} -} diff --git a/man/get_comparison_table.Rd b/man/get_comparison_table.Rd deleted file mode 100644 index e75c8333..00000000 --- a/man/get_comparison_table.Rd +++ /dev/null @@ -1,38 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/model_comparison.R -\name{get_comparison_table} -\alias{get_comparison_table} -\title{Get Model Table Comparison -Provides a table with statistics for comparison between models and selection} -\usage{ -get_comparison_table(model_objects_list) -} -\arguments{ -\item{model_objects_list}{model_objects to compare} -} -\value{ -comparison table -} -\description{ -Get Model Table Comparison -Provides a table with statistics for comparison between models and selection -} -\examples{ -\dontrun{ -data_test <- prepare_seroprev_data(serodata) -model_0 <- run_seroprev_model(seroprev_data = data_test, - seroprev_model_name = "constant_foi_bi", - n_iters = 1000) - -model_1 <- run_seroprev_model(seroprev_data = data_test, - seroprev_model_name = "continuous_foi_normal_bi", - n_iters = 1000) - -model_2 <- run_seroprev_model(seroprev_data = data_test, - seroprev_model_name = "continuous_foi_normal_log", - n_iters = 1000) -comp_table <- get_comparison_table(model_objects_list = c(m0 = model_0, - m1 = model_1, - m2 = model_2)) -} -} diff --git a/man/get_exposure_ages.Rd b/man/get_exposure_ages.Rd index 71eeb3ed..de2cec77 100644 --- a/man/get_exposure_ages.Rd +++ b/man/get_exposure_ages.Rd @@ -2,22 +2,24 @@ % Please edit documentation in R/modelling.R \name{get_exposure_ages} \alias{get_exposure_ages} -\title{Get exposure years} +\title{Function that generates an atomic vector containing the corresponding exposition years of a serological survey.} \usage{ -get_exposure_ages(seroprev_data) +get_exposure_ages(serodata) } \arguments{ -\item{seroprev_data}{A data frame containing the data from a seroprevalence survey. This data frame must contain the year of birth for each individual (birth_year) and the time of the survey (tsur). birth_year can be constructed by means of the \link{prepare_seroprev_data} function.} +\item{serodata}{A data frame containing the data from a seroprevalence survey. This data frame must contain the year of birth for each individual (birth_year) and the time of the survey (tsur). birth_year can be constructed by means of the \link{prepare_serodata} function.} } \value{ -\code{exposure_ages}. An atomic vector with the numeration of the exposition years in seroprev_data +\code{exposure_ages}. An atomic vector with the numeration of the exposition years in serodata } \description{ -Function that generates an atomic vector with the exposition years in seroprev_data. The exposition years to the disease for each individual corresponds to the time from birth to the moment of the survey. +This function generates an atomic vector containing the exposition years corresponding to the specified serological survey data \code{serodata}. +The exposition years to the disease for each individual corresponds to the time from birth to the moment of the survey. } \examples{ \dontrun{ -seroprev_data <- prepare_seroprev_data(seroprev_data = serodata, alpha = 0.05) -exposure_ages <- get_exposure_ages(seroprev_data) +data("serodata") +serodata <- prepare_serodata(serodata = serodata, alpha = 0.05) +exposure_ages <- get_exposure_ages(serodata) } } diff --git a/man/get_exposure_matrix.Rd b/man/get_exposure_matrix.Rd index e988a53a..8be080d2 100644 --- a/man/get_exposure_matrix.Rd +++ b/man/get_exposure_matrix.Rd @@ -2,22 +2,23 @@ % Please edit documentation in R/modelling.R \name{get_exposure_matrix} \alias{get_exposure_matrix} -\title{Get Exposure Matrix} +\title{Function that generates the exposure matrix corresponding to a serological survey.} \usage{ -get_exposure_matrix(seroprev_data) +get_exposure_matrix(serodata) } \arguments{ -\item{seroprev_data}{A data frame containing the data from a seroprevalence survey. This data frame must contain the year of birth for each individual (birth_year) and the time of the survey (tsur). birth_year can be constructed by means of the \link{prepare_seroprev_data} function.} +\item{serodata}{A data frame containing the data from a seroprevalence survey. This data frame must contain the year of birth for each individual (birth_year) and the time of the survey (tsur). birth_year can be constructed by means of the \link{prepare_serodata} function.} } \value{ -\code{exposure_output}. An atomic matrix containing the expositions for each entry of \code{seroprev_data} by year. +\code{exposure_output}. An atomic matrix containing the expositions for each entry of \code{serodata} by year. } \description{ -Function that generates the exposure matrix for a seroprevalence survey. +Function that generates the exposure matrix corresponding to the specified serological survey data \code{serodata}. } \examples{ \dontrun{ -seroprev_data <- prepare_seroprev_data(seroprev_data = serodata, alpha = 0.05) -exposure_matrix <- get_exposure_matrix(seroprev_data = seroprev_data) +data("serodata") +serodata <- prepare_serodata(serodata = serodata) +exposure_matrix <- get_exposure_matrix(serodata = serodata) } } diff --git a/man/get_prev_expanded.Rd b/man/get_prev_expanded.Rd index 94549724..12b5e071 100644 --- a/man/get_prev_expanded.Rd +++ b/man/get_prev_expanded.Rd @@ -2,27 +2,29 @@ % Please edit documentation in R/modelling.R \name{get_prev_expanded} \alias{get_prev_expanded} -\title{Get Prevalence Expanded} +\title{Auxiliary function that generates an object containing the confidence interval based on a +Force-of-Infection fitting performed for a specified serological survey data.} \usage{ -get_prev_expanded(foi, seroprev_data) +get_prev_expanded(foi, serodata) } \arguments{ -\item{foi}{Object containing the information of the force of infection. It is obtained from \code{rstan::extract(model_object$fit, "foi", inc_warmup = FALSE)[[1]]}.} +\item{foi}{Object containing the information of the force of infection. It is obtained from \code{rstan::extract(seromodel_object$fit, "foi", inc_warmup = FALSE)[[1]]}.} -\item{seroprev_data}{A data frame containing the data from a seroprevalence survey. For further details refer to \link{run_seroprev_model}.} +\item{serodata}{A data frame containing the data from a seroprevalence survey. For further details refer to \link{run_seromodel}.} } \value{ \code{prev_final}. The expanded prevalence data. This is used for plotting purposes in the \code{visualization} module. } \description{ -Function that generates the expanded prevalence +This function computes the corresponding binomial confidence intervals for the obtained prevalence based on a fitting +of the Force-of-Infection \code{foi} for plotting an analysis purposes. } \examples{ \dontrun{ -seroprev_data <- prepare_seroprev_data(serodata) -model_object <- run_seroprev_model(seroprev_data = seroprev_data, - seroprev_model_name = "constant_foi_bi") -foi <- rstan::extract(model_object$fit, "foi", inc_warmup = FALSE)[[1]] -get_prev_expanded <- function(foi, seroprev_data) +serodata <- prepare_serodata(serodata) +seromodel_object <- run_seromodel(serodata = serodata, + seromodel_name = "constant_foi_bi") +foi <- rstan::extract(seromodel_object$fit, "foi")[[1]] +get_prev_expanded <- function(foi, serodata) } } diff --git a/man/get_sim_counts.Rd b/man/get_sim_counts.Rd deleted file mode 100644 index b75538d1..00000000 --- a/man/get_sim_counts.Rd +++ /dev/null @@ -1,29 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/seroprevalence_data.R -\name{get_sim_counts} -\alias{get_sim_counts} -\title{Function that randomly generates a sample of counts for a simulated dataset} -\usage{ -get_sim_counts(sim_data, foi, size_age_class, seed = 1234) -} -\arguments{ -\item{sim_data}{A dataframe object containing the following columns: -\tabular{ll}{ -\code{birth_year} \tab List of years in which the subjects were borned \cr \tab \cr -\code{tsur} \tab Year of the survey\cr \tab \cr -\code{country} \tab Default to 'none'.\cr \tab \cr -\code{survey} \tab Survey label \cr \tab \cr -\code{age_mean_f} \tab Age \cr \tab \cr -}} -} -\value{ -A simulated list of counts following a binomial distribution in accordance with a given force of infection and age class sizes. -} -\description{ -Function that randomly generates a sample of counts for a simulated dataset -} -\examples{ -\dontrun{ - -} -} diff --git a/man/get_table_rhats.Rd b/man/get_table_rhats.Rd index 5a5ea6d7..92bd292a 100644 --- a/man/get_table_rhats.Rd +++ b/man/get_table_rhats.Rd @@ -2,24 +2,27 @@ % Please edit documentation in R/model_comparison.R \name{get_table_rhats} \alias{get_table_rhats} -\title{Get Table Rhats} +\title{Method for extracting a dataframe containing the R-hat estimates for a given serological model.} \usage{ -get_table_rhats(model_object) +get_table_rhats(seromodel_object) } \arguments{ -\item{model_object}{model_object} +\item{seromodel_object}{seromodel_object} } \value{ rhats table } \description{ -Function that makes the rhats table +This method relies in the function \link[bayesplot]{rhat} to extract the R-hat estimates of the serological model object +\code{seromodel_object} and returns a table a dataframe with the estimates for each year of birth. } \examples{ \dontrun{ -seroprev_data <- prepare_seroprev_data(seroprev_data = serodata, alpha = 0.05) -model_object <- run_seroprev_model( -seroprev_data = seroprev_data, seroprev_model_name = "constant_foi_bi") -get_table_rhats (model_object) +data("serodata") +data_test <- prepare_serodata(serodata = serodata) +model_constant <- run_seromodel(serodata = data_test, + seromodel_name = "constant_foi_bi", + n_iters = 1500) +get_table_rhats(model_object = model_constant) } } diff --git a/man/group_sim_data.Rd b/man/group_sim_data.Rd deleted file mode 100644 index fcde7ad0..00000000 --- a/man/group_sim_data.Rd +++ /dev/null @@ -1,44 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/seroprevalence_data.R -\name{group_sim_data} -\alias{group_sim_data} -\title{Function that generates grouped simulated data from a given Force-of-Infection} -\usage{ -group_sim_data( - sim_data, - foi, - size_age_class, - tsur, - birth_year_min, - survey_label, - test = "fake", - antibody = "IgG", - seed = 1234 -) -} -\arguments{ -\item{sim_data}{Dataframe with the structure of the output of \code{\linl{generate_sim_data}}.} -} -\value{ -Dataframe object containing grouped simulated data generated from \code{foi} -} -\description{ -Function that generates grouped simulated data from a given Force-of-Infection -} -\examples{ -\dontrun{ -size_age_class = 5 -foi <- rep(0.02, 50) -sim_data <- generate_sim_data(foi = foi, - size_age_class = size_age_class, - tsur = 2050, - birth_year_min = 2000, - survey_label = 'sim_constant_foi') -sim_data_grouped <- group_sim_data(sim_data = sim_data, - foi = foi, - size_age_class = size_age_class, - tsur = 2050, - birth_year_min = 2000, - survey_label = 'sim_constant_foi_grouped') -} -} diff --git a/man/plot_foi.Rd b/man/plot_foi.Rd index 45d0e987..da19fde2 100644 --- a/man/plot_foi.Rd +++ b/man/plot_foi.Rd @@ -2,12 +2,12 @@ % Please edit documentation in R/visualization.R \name{plot_foi} \alias{plot_foi} -\title{Generate Force-of-Infection Plot} +\title{Function that generates a Force-of-Infection plot corresponding to the specified fitted serological model.} \usage{ -plot_foi(model_object, lambda_sim = NA, max_lambda = NA, size_text = 25) +plot_foi(seromodel_object, lambda_sim = NA, max_lambda = NA, size_text = 25) } \arguments{ -\item{model_object}{Object containing the results of fitting a model by means of \link{run_seroprev_model}.} +\item{seromodel_object}{Object containing the results of fitting a model by means of \link{run_seromodel}.} \item{size_text}{Text size use in the theme of the graph returned by the function.} } @@ -15,16 +15,18 @@ plot_foi(model_object, lambda_sim = NA, max_lambda = NA, size_text = 25) A ggplot2 object containing the Force-of-infection-vs-time including the corresponding confidence interval. } \description{ -Function that generates the Force-of-Infection plot. The x axis corresponds to the decades covered by the survey the y axis to the Force-of-Infection. +This function generates a Force-of-Infection plot from the results obtained by fitting a serological model. +This includes the corresponding binomial confidence interval. +The x axis corresponds to the decades covered by the survey the y axis to the Force-of-Infection. } \examples{ \dontrun{ - data_test <- prepare_seroprev_data(serodata) - model_object <- run_seroprev_model( - seroprev_data = data_test, - seroprev_model_name = "constant_foi_bi", + data_test <- prepare_serodata(serodata) + seromodel_object <- run_seromodel( + serodata = data_test, + seromodel_name = "constant_foi_bi", n_iters = 1000 ) -plot_foi(model_object, size_text = 15) +plot_foi(seromodel_object, size_text = 15) } } diff --git a/man/plot_info_table.Rd b/man/plot_info_table.Rd index f30f7e5e..56791995 100644 --- a/man/plot_info_table.Rd +++ b/man/plot_info_table.Rd @@ -2,7 +2,7 @@ % Please edit documentation in R/visualization.R \name{plot_info_table} \alias{plot_info_table} -\title{Plot Info Table} +\title{Auxiliary function that generates a plot of a given table} \usage{ plot_info_table(info, size_text) } @@ -15,17 +15,17 @@ plot_info_table(info, size_text) p, a variable that will be used in the \link{visualisation} module } \description{ -Function that generates the information table +Auxiliary function that generates a plot of a given table } \examples{ \dontrun{ - data_test <- prepare_seroprev_data(serodata) - model_object <- run_seroprev_model( - seroprev_data = data_test, - seroprev_model_name = "constant_foi_bi", + data_test <- prepare_serodata(serodata) + seromodel_object <- run_seromodel( + serodata = data_test, + seromodel_name = "constant_foi_bi", n_iters = 1000 ) -info = t(model_object$model_summary) +info = t(seromodel_object$model_summary) plot_info_table (info, size_text = 15) } } diff --git a/man/plot_rhats.Rd b/man/plot_rhats.Rd index e79241c8..0f0563db 100644 --- a/man/plot_rhats.Rd +++ b/man/plot_rhats.Rd @@ -2,12 +2,12 @@ % Please edit documentation in R/visualization.R \name{plot_rhats} \alias{plot_rhats} -\title{Generate Rhats-Convergence Plot} +\title{Function that generates a plot of the R-hat estimates of the specified fitted serological model.} \usage{ -plot_rhats(model_object, size_text = 25) +plot_rhats(seromodel_object, size_text = 25) } \arguments{ -\item{model_object}{Object containing the results of fitting a model by means of \link{run_seroprev_model}.} +\item{seromodel_object}{Object containing the results of fitting a model by means of \link{run_seromodel}.} \item{size_text}{Text size use in the theme of the graph returned by the function.} } @@ -15,17 +15,20 @@ plot_rhats(model_object, size_text = 25) The rhats-convergence plot of the selected model. } \description{ -Function that generates the convergence graph of a model. The x axis corresponds to the decades covered by the survey and the y axis to the value of the rhats. -All rhats must be smaller than 1 to ensure convergence. +This function generates a plot of the R-hat estimates obtained for a specified fitted serological model \code{seromodel_object}. +The x axis corresponds to the decades covered by the survey and the y axis to the value of the rhats. +All rhats must be smaller than 1 to ensure convergence (for further details check \link[bayesplot]{rhat}). } \examples{ \dontrun{ -data_test <- prepare_seroprev_data(serodata) -model_object <- run_seroprev_model( - seroprev_data = data_test, - seroprev_model_name = "constant_foi_bi", +data("serodata") +data_test <- prepare_serodata(serodata) +seromodel_object <- run_seromodel( + serodata = data_test, + seromodel_name = "constant_foi_bi", n_iters = 1000 ) -plot_rhats(model_object, size_text = 15) +plot_rhats(seromodel_object, + size_text = 15) } } diff --git a/man/plot_seromodel.Rd b/man/plot_seromodel.Rd new file mode 100644 index 00000000..4275d59c --- /dev/null +++ b/man/plot_seromodel.Rd @@ -0,0 +1,37 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/visualization.R +\name{plot_seromodel} +\alias{plot_seromodel} +\title{Function that generates a vertical arrange of plots showing a summary of a given implemented model, as well as its corresponding +seroprevalence, Force-of-Infection and R-hat estimates plots.} +\usage{ +plot_seromodel( + seromodel_object, + lambda_sim = NA, + max_lambda = NA, + size_text = 25 +) +} +\arguments{ +\item{seromodel_object}{Object containing the results of fitting a model by means of \link{run_seromodel}.} + +\item{size_text}{Text size use in the theme of the graph returned by the function.} +} +\value{ +A ggplot object with a vertical arrange containing the seropositivity, force of infection, and convergence plots. +} +\description{ +Function that generates a vertical arrange of plots showing a summary of a given implemented model, as well as its corresponding +seroprevalence, Force-of-Infection and R-hat estimates plots. +} +\examples{ +\dontrun{ +data_test <- prepare_serodata(serodata) +seromodel_object <- run_seromodel( + serodata = data_test, + seromodel_name = "constant_foi_bi", + n_iters = 1000 +) +plot_seromodel(seromodel_object, size_text = 15) +} +} diff --git a/man/plot_seroprev.Rd b/man/plot_seroprev.Rd index 90fff635..e77fc305 100644 --- a/man/plot_seroprev.Rd +++ b/man/plot_seroprev.Rd @@ -2,12 +2,12 @@ % Please edit documentation in R/visualization.R \name{plot_seroprev} \alias{plot_seroprev} -\title{Generate sero-positivity plot from raw data} +\title{Function that generates the sero-positivity plot from a raw serological survey dataset.} \usage{ -plot_seroprev(seroprev_data, size_text = 6) +plot_seroprev(serodata, size_text = 6) } \arguments{ -\item{seroprev_data}{A data frame containing the data from a seroprevalence survey. +\item{serodata}{A data frame containing the data from a seroprevalence survey. This data frame must contain the following columns: \tabular{ll}{ \code{survey} \tab survey Label of the current survey \cr \tab \cr @@ -29,16 +29,16 @@ This data frame must contain the following columns: A ggplot object containing the seropositivity-vs-age graph of the raw data of a given seroprevalence survey with its corresponging binomial confidence interval. } \description{ -Function that generates the sero positivity plot from raw data +Function that generates the sero-positivity plot from a raw serological survey dataset. } \examples{ \dontrun{ - data_test <- prepare_seroprev_data(serodata) - model_object <- run_seroprev_model( - seroprev_data = data_test, - seroprev_model_name = "constant_foi_bi", + data_test <- prepare_serodata(serodata) + seromodel_object <- run_seromodel( + serodata = data_test, + seromodel_name = "constant_foi_bi", n_iters = 1000 ) -plot_seroprev(model_object, size_text = 15) +plot_seroprev(seromodel_object, size_text = 15) } } diff --git a/man/plot_seroprev_fitted.Rd b/man/plot_seroprev_fitted.Rd index 72d84c52..c793190e 100644 --- a/man/plot_seroprev_fitted.Rd +++ b/man/plot_seroprev_fitted.Rd @@ -2,12 +2,12 @@ % Please edit documentation in R/visualization.R \name{plot_seroprev_fitted} \alias{plot_seroprev_fitted} -\title{Generate sero-positivity plot with fitted model} +\title{Function that generates a seropositivity plot corresponding to the specified fitted serological model.} \usage{ -plot_seroprev_fitted(model_object, size_text = 6) +plot_seroprev_fitted(seromodel_object, size_text = 6) } \arguments{ -\item{model_object}{Object containing the results of fitting a model by means of \link{run_seroprev_model}.} +\item{seromodel_object}{Object containing the results of fitting a model by means of \link{run_seromodel}.} \item{size_text}{Text size of the graph returned by the function.} } @@ -15,16 +15,17 @@ plot_seroprev_fitted(model_object, size_text = 6) A ggplot object containing the seropositivity-vs-age graph including the data, the fitted model and their corresponding confindence intervals. } \description{ -Function that generates the seropositivity graph with fitted model. Age is located on the x axis and seropositivity on the y axis with its corresponding confidence interval. +This function generates a seropositivity plot of the specified serological model object. This includes the original data grouped by age +as well as the obtained fitting from the model implementation. Age is located on the x axis and seropositivity on the y axis with its +corresponding confidence interval. } \examples{ \dontrun{ - data_test <- prepare_seroprev_data(serodata) - model_object <- run_seroprev_model( - seroprev_data = data_test, - seroprev_model_name = "constant_foi_bi", - n_iters = 1000 -) -plot_seroprev_fitted(model_object, size_text = 15) +data("serodata") +data_test <- prepare_serodata(serodata) +seromodel_object <- run_seromodel(serodata = data_test, + seromodel_name = "constant_foi_bi", + n_iters = 1000) +plot_seroprev_fitted(seromodel_object, size_text = 15) } } diff --git a/man/plot_seroprev_model.Rd b/man/plot_seroprev_model.Rd deleted file mode 100644 index ae613bb5..00000000 --- a/man/plot_seroprev_model.Rd +++ /dev/null @@ -1,35 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/visualization.R -\name{plot_seroprev_model} -\alias{plot_seroprev_model} -\title{Generate a vertical arrange of plots summarizing the results of the model implementation} -\usage{ -plot_seroprev_model( - model_object, - lambda_sim = NA, - max_lambda = NA, - size_text = 25 -) -} -\arguments{ -\item{model_object}{Object containing the results of fitting a model by means of \link{run_seroprev_model}.} - -\item{size_text}{Text size use in the theme of the graph returned by the function.} -} -\value{ -A ggplot object with a vertical arrange containing the seropositivity, force of infection, and convergence plots. -} -\description{ -Function that generates the combined plots summarizing the results of the model implementation -} -\examples{ -\dontrun{ -data_test <- prepare_seroprev_data(serodata) -model_object <- run_seroprev_model( - seroprev_data = data_test, - seroprev_model_name = "constant_foi_bi", - n_iters = 1000 -) -plot_seroprev_model(model_object, size_text = 15) -} -} diff --git a/man/plot_seroprev_models_grid.Rd b/man/plot_seroprev_models_grid.Rd deleted file mode 100644 index 6205ccd4..00000000 --- a/man/plot_seroprev_models_grid.Rd +++ /dev/null @@ -1,34 +0,0 @@ -% Generated by roxygen2: do not edit by hand -% Please edit documentation in R/visualization.R -\name{plot_seroprev_models_grid} -\alias{plot_seroprev_models_grid} -\title{Function that generates a plot showing the results for several models.} -\usage{ -plot_seroprev_models_grid(..., n_row = NULL, n_col = NULL) -} -\arguments{ -\item{...}{The first n arguments of the function are taken as plots ti be arranged in a single plot. -This objects can be obtained by means of \link{run_seroprev_model}.} - -\item{n_row}{Number of rows of the plot array.} - -\item{n_col}{Number of columns of the plot array.} -} -\value{ -A ggplot object containing an array with the plots in \code{models_list}. -} -\description{ -Function that generates a plot showing the results for several models. -} -\examples{ -\dontrun{ - data_test <- prepare_seroprev_data(serodata) - model_object_constant <- run_seroprev_model(seroprev_data = data_test, - seroprev_model_name = "constant_foi_bi") - model_object_normalbi <- run_seroprev_model(seroprev_data = data_test, - seroprev_model_name = "continuous_foi_normal_bi") - model_object_normallog <- run_seroprev_model(seroprev_data = data_test, - seroprev_model_name = "continuous_foi_normal_log") - plot_seroprev_models_list(model_object_constant, model_object_normalbi, model_object_normallog, n_row = 1, n_col = 3) -} -} diff --git a/man/prepare_bin_data.Rd b/man/prepare_bin_data.Rd index 4bd20caa..91934657 100644 --- a/man/prepare_bin_data.Rd +++ b/man/prepare_bin_data.Rd @@ -2,12 +2,13 @@ % Please edit documentation in R/seroprevalence_data.R \name{prepare_bin_data} \alias{prepare_bin_data} -\title{Prepare data to plot binomial confidence intervals} +\title{Function that prepares a pre-processed serological survey dataset to plot the binomial confidence intervals of the seroprevalence grouped by +age group.} \usage{ -prepare_bin_data(seroprev_data) +prepare_bin_data(serodata) } \arguments{ -\item{seroprev_data}{A data frame containing the data from a seroprevalence survey. For more information see the function \link{run_seroprev_model}. +\item{serodata}{A data frame containing the data from a seroprevalence survey. For more information see the function \link{run_seromodel}. This data frame must contain the following columns: \tabular{ll}{ \code{survey} \tab survey Label of the current survey \cr \tab \cr @@ -28,16 +29,17 @@ This data frame must contain the following columns: \code{prev_obs_lower} \tab Lower limit of the confidence interval for the observed prevalence \cr \tab \cr \code{prev_obs_upper} \tab Upper limit of the confidence interval for the observed prevalence \cr \tab \cr } -The last six colums can be added to \code{seroprev_data} by means of the function \code{\link{prepare_seroprev_data}}.} +The last six colums can be added to \code{serodata} by means of the function \code{\link{prepare_serodata}}.} } \value{ data set with the binomial confidence intervals } \description{ -Function that prepares the data for modelling +This function prepapares a given pre-processed serological dataset (see \code{\link{prepare_serodata}}) to plot the binomial confidence intervals +of its corresponding seroprevalence grouped by age group. } \examples{ \dontrun{ -prepare_bin_data (seroprev_data) +prepare_bin_data (serodata) } } diff --git a/man/prepare_seroprev_data.Rd b/man/prepare_serodata.Rd similarity index 73% rename from man/prepare_seroprev_data.Rd rename to man/prepare_serodata.Rd index 208264ca..81da358a 100644 --- a/man/prepare_seroprev_data.Rd +++ b/man/prepare_serodata.Rd @@ -1,17 +1,13 @@ % Generated by roxygen2: do not edit by hand % Please edit documentation in R/seroprevalence_data.R -\name{prepare_seroprev_data} -\alias{prepare_seroprev_data} -\title{Prepare data} +\name{prepare_serodata} +\alias{prepare_serodata} +\title{Function that prepares the data from a serological survey for modelling.} \usage{ -prepare_seroprev_data( - seroprev_data = serodata, - alpha = 0.05, - add_age_mean_f = TRUE -) +prepare_serodata(serodata = serodata, alpha = 0.05, add_age_mean_f = TRUE) } \arguments{ -\item{seroprev_data}{A data frame containing the data from a seroprevalence survey. +\item{serodata}{A data frame containing the data from a serological survey. This data frame must contain the following columns: \tabular{ll}{ \code{survey} \tab survey Label of the current survey \cr \tab \cr @@ -30,7 +26,7 @@ This data frame must contain the following columns: \item{alpha}{probability of a type I error. For further details refer to \link{Hmisc::binconf}.} } \value{ -seroprev_data with additional columns necessary for the analysis. These columns are: +serodata with additional columns necessary for the analysis. These columns are: \tabular{ll}{ \code{age_mean_f} \tab Floor value of the average between age_min and age_max \cr \tab \cr \code{sample_size} \tab The size of the sample \cr \tab \cr @@ -41,11 +37,11 @@ seroprev_data with additional columns necessary for the analysis. These columns } } \description{ -Function that prepares the data for modelling +This function adds the necessary additional variables to the given dataset \code{serodata} corresponding to a serological survey. } \examples{ \dontrun{ -data_test <- readRDS("data/data.RDS") -data_test <- prepare_seroprev_data(seroprev_data, alpha) +data("serodata") +data_test <- prepare_serodata(serodata) } } diff --git a/man/run_seroprev_model.Rd b/man/run_seromodel.Rd similarity index 70% rename from man/run_seroprev_model.Rd rename to man/run_seromodel.Rd index 89923bc5..ff3d1bca 100644 --- a/man/run_seroprev_model.Rd +++ b/man/run_seromodel.Rd @@ -1,21 +1,22 @@ % Generated by roxygen2: do not edit by hand % Please edit documentation in R/modelling.R -\name{run_seroprev_model} -\alias{run_seroprev_model} -\title{Run the specified stan model for the force-of-infection} +\name{run_seromodel} +\alias{run_seromodel} +\title{Run the specified stan model for the Force-of-Infection and estimates de seroprevalence based on the result of the fit.} \usage{ -run_seroprev_model( - seroprev_data, - seroprev_model_name = "constant_foi_bi", +run_seromodel( + serodata, + seromodel_name = "constant_foi_bi", n_iters = 1000, n_thin = 2, delta = 0.9, m_treed = 10, - decades = 0 + decades = 0, + print_summary = TRUE ) } \arguments{ -\item{seroprev_data}{A data frame containing the data from a seroprevalence survey. +\item{serodata}{A data frame containing the data from a seroprevalence survey. This data frame must contain the following columns: \tabular{ll}{ \code{survey} \tab survey Label of the current survey \cr \tab \cr @@ -36,9 +37,9 @@ This data frame must contain the following columns: \code{prev_obs_lower} \tab Lower limit of the confidence interval for the observed prevalence \cr \tab \cr \code{prev_obs_upper} \tab Upper limit of the confidence interval for the observed prevalence \cr \tab \cr } -The last six colums can be added to \code{seroprev_data} by means of the function \code{\link{prepare_seroprev_data}}.} +The last six colums can be added to \code{serodata} by means of the function \code{\link{prepare_serodata}}.} -\item{seroprev_model_name}{Name of the selected model. Current version provides three options: +\item{seromodel_name}{Name of the selected model. Current version provides three options: \describe{ \item{\code{"constant_foi_bi"}}{Runs a constant model} \item{\code{"continuous_foi_normal_bi"}}{Runs a normal model} @@ -58,15 +59,16 @@ For further details refer to the \code{control} parameter in \link[rstan]{sampli \item{decades}{Number of decades covered by the survey data.} } \value{ -\code{model_object}. An object containing relevant information about the implementation of the model. For further details refer to \link{fit_seroprev_model}. +\code{seromodel_object}. An object containing relevant information about the implementation of the model. For further details refer to \link{fit_seromodel}. } \description{ -Run the specified stan model for the force-of-infection +This function runs the specified model for the Force-of-Infection \code{seromodel_name} using the data froma seroprevalence survey +\code{serodata} as the input data. See \link{fit_seromodel} for further details. } \examples{ \dontrun{ -seroprev_data <- prepare_seroprev_data(serodata) -run_seroprev_model (seroprev_data, - seroprev_model_name = "constant_foi_bi") +serodata <- prepare_serodata(serodata) +run_seromodel (serodata, + seromodel_name = "constant_foi_bi") } } diff --git a/man/save_or_load_model.Rd b/man/save_or_load_model.Rd index 49832931..2866e6fb 100644 --- a/man/save_or_load_model.Rd +++ b/man/save_or_load_model.Rd @@ -2,12 +2,12 @@ % Please edit documentation in R/modelling.R \name{save_or_load_model} \alias{save_or_load_model} -\title{Save or load model} +\title{Auxiliar function used to determine whether the stan model corresponding to the specified serological model has been already compiled or not.} \usage{ -save_or_load_model(seroprev_model_name = "constant_foi_bi") +save_or_load_model(seromodel_name = "constant_foi_bi") } \arguments{ -\item{seroprev_model_name}{Name of the selected model. Current version provides three options: +\item{seromodel_name}{Name of the selected model. Current version provides three options: \describe{ \item{\code{"constant_foi_bi"}}{Runs a constant model} \item{\code{"continuous_foi_normal_bi"}}{Runs a normal model} @@ -23,5 +23,8 @@ In case the .RDS file exists, it is read and returned; otherwise, the object mod \link[rstan]{stan_model} function, saved as an .RDS file and returned as the output of the function. } \examples{ -save_or_load_model(seroprev_model_name="constant_foi_bi") +\dontrun{ +model <- save_or_load_model(seromodel_name="constant_foi_bi") +} + } diff --git a/tests/testthat/_snaps/comparison/comp_table.csv b/tests/testthat/_snaps/comparison/comp_table.csv deleted file mode 100644 index d70bd29d..00000000 --- a/tests/testthat/_snaps/comparison/comp_table.csv +++ /dev/null @@ -1,4 +0,0 @@ -"","seroprev_model_name","dataset","country","year","test","antibody","n_sample","n_agec","n_iter","elpd","se","converged","difference","diff_se","best_elpd","pvalue" -"1","constant_foi_bi","COL-035-93","COL",2012,"ELISA","IgG anti-T.cruzi",747,72,1000,-92.75,6.4,"Yes",0,1,NA,0.500041 -"2","continuous_foi_normal_bi","COL-035-93","COL",2012,"ELISA","IgG anti-T.cruzi",747,72,1000,-74.46,5.88,"Yes",0,-18.2967809927967,NA,0.500009 -"3","continuous_foi_normal_log","COL-035-93","COL",2012,"ELISA","IgG anti-T.cruzi",747,72,1000,-71.93,7.11,"Yes",0,-20.8238925781359,NA,0.500029 diff --git a/tests/testthat/_snaps/individual_models/seroprev-fitted-individual-plot.svg b/tests/testthat/_snaps/individual_models/seroprev-fitted-individual-plot.svg deleted file mode 100644 index 7f5d9fab..00000000 --- a/tests/testthat/_snaps/individual_models/seroprev-fitted-individual-plot.svg +++ /dev/null @@ -1,139 +0,0 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -0.00 -0.25 -0.50 -0.75 -1.00 - - - - - - - - - -0 -20 -40 -60 -Age -Sero-positivity -seroprev_fitted_individual_plot - - diff --git a/tests/testthat/_snaps/models/constant-foi-plot.svg b/tests/testthat/_snaps/models/constant-foi-plot.svg new file mode 100644 index 00000000..58d67bf1 --- /dev/null +++ b/tests/testthat/_snaps/models/constant-foi-plot.svg @@ -0,0 +1,80 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +0.000 +0.001 +0.002 +0.003 +0.004 +0.005 + + + + + + + + + + +1940 +1960 +1980 +2000 +Year +Force-of-Infection +constant_foi_plot + + diff --git a/tests/testthat/_snaps/individual_models/model-0-plot.svg b/tests/testthat/_snaps/models/constant-model-plot.svg similarity index 99% rename from tests/testthat/_snaps/individual_models/model-0-plot.svg rename to tests/testthat/_snaps/models/constant-model-plot.svg index 2f930530..7bcc61d5 100644 --- a/tests/testthat/_snaps/individual_models/model-0-plot.svg +++ b/tests/testthat/_snaps/models/constant-model-plot.svg @@ -25,7 +25,7 @@ -seroprev_model_name: constant_foi_bi +seromodel_name: constant_foi_bi dataset: COL-035-93 elpd: -92.75 se: 6.4 diff --git a/tests/testthat/_snaps/models/constant-rhats-plot.svg b/tests/testthat/_snaps/models/constant-rhats-plot.svg new file mode 100644 index 00000000..975f6bbf --- /dev/null +++ b/tests/testthat/_snaps/models/constant-rhats-plot.svg @@ -0,0 +1,144 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +1.0 +1.5 +2.0 + + + + + + + +1940 +1960 +1980 +2000 +year +Convergence (R^) +constant_rhats_plot + + diff --git a/tests/testthat/_snaps/plot_functions/plot-seroprev-fitted-constant.svg b/tests/testthat/_snaps/models/constant-sp-fitted-plot.svg similarity index 99% rename from tests/testthat/_snaps/plot_functions/plot-seroprev-fitted-constant.svg rename to tests/testthat/_snaps/models/constant-sp-fitted-plot.svg index be7f9fac..3b34e920 100644 --- a/tests/testthat/_snaps/plot_functions/plot-seroprev-fitted-constant.svg +++ b/tests/testthat/_snaps/models/constant-sp-fitted-plot.svg @@ -134,6 +134,6 @@ 60 Age Sero-positivity -plot_seroprev_fitted_constant +constant_sp_fitted_plot diff --git a/tests/testthat/_snaps/models/normal-foi-plot.svg b/tests/testthat/_snaps/models/normal-foi-plot.svg new file mode 100644 index 00000000..4df3bb59 --- /dev/null +++ b/tests/testthat/_snaps/models/normal-foi-plot.svg @@ -0,0 +1,77 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +0.000 +0.005 +0.010 +0.015 +0.020 + + + + + + + + + +1940 +1960 +1980 +2000 +Year +Force-of-Infection +normal_foi_plot + + diff --git a/tests/testthat/_snaps/individual_models/foi-individual-plot.svg b/tests/testthat/_snaps/models/normal-log-foi-plot.svg similarity index 99% rename from tests/testthat/_snaps/individual_models/foi-individual-plot.svg rename to tests/testthat/_snaps/models/normal-log-foi-plot.svg index 7b42cfb5..26e211a0 100644 --- a/tests/testthat/_snaps/individual_models/foi-individual-plot.svg +++ b/tests/testthat/_snaps/models/normal-log-foi-plot.svg @@ -68,6 +68,6 @@ 2000 Year Force-of-Infection -foi_individual_plot +normal_log_foi_plot diff --git a/tests/testthat/_snaps/individual_models/model-2-plot.svg b/tests/testthat/_snaps/models/normal-log-model-plot.svg similarity index 99% rename from tests/testthat/_snaps/individual_models/model-2-plot.svg rename to tests/testthat/_snaps/models/normal-log-model-plot.svg index 718153a9..ba2f317f 100644 --- a/tests/testthat/_snaps/individual_models/model-2-plot.svg +++ b/tests/testthat/_snaps/models/normal-log-model-plot.svg @@ -25,7 +25,7 @@ -seroprev_model_name: continuous_foi_normal_log +seromodel_name: continuous_foi_normal_log dataset: COL-035-93 elpd: -71.93 se: 7.11 diff --git a/tests/testthat/_snaps/individual_models/rhats-individual-plot.svg b/tests/testthat/_snaps/models/normal-log-rhats-plot.svg similarity index 99% rename from tests/testthat/_snaps/individual_models/rhats-individual-plot.svg rename to tests/testthat/_snaps/models/normal-log-rhats-plot.svg index 920171f5..e3e2acde 100644 --- a/tests/testthat/_snaps/individual_models/rhats-individual-plot.svg +++ b/tests/testthat/_snaps/models/normal-log-rhats-plot.svg @@ -139,6 +139,6 @@ 2000 year Convergence (R^) -rhats_individual_plot +normal_log_rhats_plot diff --git a/tests/testthat/_snaps/plot_functions/plot-seroprev-fitted-normallog.svg b/tests/testthat/_snaps/models/normal-log-sp-fitted-plot.svg similarity index 99% rename from tests/testthat/_snaps/plot_functions/plot-seroprev-fitted-normallog.svg rename to tests/testthat/_snaps/models/normal-log-sp-fitted-plot.svg index b2bed6c4..2045488a 100644 --- a/tests/testthat/_snaps/plot_functions/plot-seroprev-fitted-normallog.svg +++ b/tests/testthat/_snaps/models/normal-log-sp-fitted-plot.svg @@ -134,6 +134,6 @@ 60 Age Sero-positivity -plot_seroprev_fitted_normallog +normal_log_sp_fitted_plot diff --git a/tests/testthat/_snaps/individual_models/model-1-plot.svg b/tests/testthat/_snaps/models/normal-model-plot.svg similarity index 99% rename from tests/testthat/_snaps/individual_models/model-1-plot.svg rename to tests/testthat/_snaps/models/normal-model-plot.svg index efbdbf9c..3035fb0a 100644 --- a/tests/testthat/_snaps/individual_models/model-1-plot.svg +++ b/tests/testthat/_snaps/models/normal-model-plot.svg @@ -25,7 +25,7 @@ -seroprev_model_name: continuous_foi_normal_bi +seromodel_name: continuous_foi_normal_bi dataset: COL-035-93 elpd: -74.46 se: 5.88 diff --git a/tests/testthat/_snaps/models/normal-rhats-plot.svg b/tests/testthat/_snaps/models/normal-rhats-plot.svg new file mode 100644 index 00000000..9ac7fe59 --- /dev/null +++ b/tests/testthat/_snaps/models/normal-rhats-plot.svg @@ -0,0 +1,144 @@ + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + + +1.0 +1.5 +2.0 + + + + + + + +1940 +1960 +1980 +2000 +year +Convergence (R^) +normal_rhats_plot + + diff --git a/tests/testthat/_snaps/plot_functions/plot-seroprev-fitted-normalbi.svg b/tests/testthat/_snaps/models/normal-sp-fitted-plot.svg similarity index 99% rename from tests/testthat/_snaps/plot_functions/plot-seroprev-fitted-normalbi.svg rename to tests/testthat/_snaps/models/normal-sp-fitted-plot.svg index cd181fcd..24913f03 100644 --- a/tests/testthat/_snaps/plot_functions/plot-seroprev-fitted-normalbi.svg +++ b/tests/testthat/_snaps/models/normal-sp-fitted-plot.svg @@ -134,6 +134,6 @@ 60 Age Sero-positivity -plot_seroprev_fitted_normalbi +normal_sp_fitted_plot diff --git a/tests/testthat/_snaps/plot_functions/plot-seroprev.svg b/tests/testthat/_snaps/models/serodata-plot.svg similarity index 99% rename from tests/testthat/_snaps/plot_functions/plot-seroprev.svg rename to tests/testthat/_snaps/models/serodata-plot.svg index 603ce49d..c223291d 100644 --- a/tests/testthat/_snaps/plot_functions/plot-seroprev.svg +++ b/tests/testthat/_snaps/models/serodata-plot.svg @@ -130,6 +130,6 @@ 60 Age Sero-positivity -plot_seroprev +serodata_plot diff --git a/tests/testthat/_snaps/plot_functions/plot-arrange-models.svg b/tests/testthat/_snaps/plot_functions/plot-arrange-models.svg deleted file mode 100644 index 321b8ad8..00000000 --- a/tests/testthat/_snaps/plot_functions/plot-arrange-models.svg +++ /dev/null @@ -1,1044 +0,0 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - -seroprev_model_name: constant_foi_bi -dataset: COL-035-93 -elpd: -92.75 -se: 6.4 -converged: Yes - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -0.00 -0.25 -0.50 -0.75 -1.00 - - - - - - - - - -0 -20 -40 -60 -Age -Sero-positivity - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -0.000 -0.001 -0.002 -0.003 -0.004 -0.005 - - - - - - - - - - -1940 -1960 -1980 -2000 -Year -Force-of-Infection - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -1.0 -1.5 -2.0 - - - - - - - -1940 -1960 -1980 -2000 -year -Convergence (R^) - - - - - - - - - - - - - - - - -seroprev_model_name: continuous_foi_normal_bi -dataset: COL-035-93 -elpd: -74.46 -se: 5.88 -converged: Yes - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -0.00 -0.25 -0.50 -0.75 -1.00 - - - - - - - - - -0 -20 -40 -60 -Age -Sero-positivity - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -0.000 -0.005 -0.010 -0.015 -0.020 - - - - - - - - - -1940 -1960 -1980 -2000 -Year -Force-of-Infection - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -1.0 -1.5 -2.0 - - - - - - - -1940 -1960 -1980 -2000 -year -Convergence (R^) - - - - - - - - - - - - - - - - -seroprev_model_name: continuous_foi_normal_log -dataset: COL-035-93 -elpd: -71.93 -se: 7.11 -converged: Yes - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -0.00 -0.25 -0.50 -0.75 -1.00 - - - - - - - - - -0 -20 -40 -60 -Age -Sero-positivity - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -0.00 -0.05 -0.10 -0.15 - - - - - - - - -1940 -1960 -1980 -2000 -Year -Force-of-Infection - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -1.0 -1.5 -2.0 - - - - - - - -1940 -1960 -1980 -2000 -year -Convergence (R^) - - diff --git a/tests/testthat/_snaps/plot_functions/plot-seroprev-model-constant.svg b/tests/testthat/_snaps/plot_functions/plot-seroprev-model-constant.svg deleted file mode 100644 index 2f930530..00000000 --- a/tests/testthat/_snaps/plot_functions/plot-seroprev-model-constant.svg +++ /dev/null @@ -1,357 +0,0 @@ - - - - - - - - - - - - - - - - - - - -seroprev_model_name: constant_foi_bi -dataset: COL-035-93 -elpd: -92.75 -se: 6.4 -converged: Yes - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -0.00 -0.25 -0.50 -0.75 -1.00 - - - - - - - - - -0 -20 -40 -60 -Age -Sero-positivity - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -0.000 -0.001 -0.002 -0.003 -0.004 -0.005 - - - - - - - - - - -1940 -1960 -1980 -2000 -Year -Force-of-Infection - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -1.0 -1.5 -2.0 - - - - - - - -1940 -1960 -1980 -2000 -year -Convergence (R^) - - diff --git a/tests/testthat/_snaps/plot_functions/plot-seroprev-model-normalbi.svg b/tests/testthat/_snaps/plot_functions/plot-seroprev-model-normalbi.svg deleted file mode 100644 index efbdbf9c..00000000 --- a/tests/testthat/_snaps/plot_functions/plot-seroprev-model-normalbi.svg +++ /dev/null @@ -1,354 +0,0 @@ - - - - - - - - - - - - - - - - - - - -seroprev_model_name: continuous_foi_normal_bi -dataset: COL-035-93 -elpd: -74.46 -se: 5.88 -converged: Yes - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -0.00 -0.25 -0.50 -0.75 -1.00 - - - - - - - - - -0 -20 -40 -60 -Age -Sero-positivity - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -0.000 -0.005 -0.010 -0.015 -0.020 - - - - - - - - - -1940 -1960 -1980 -2000 -Year -Force-of-Infection - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -1.0 -1.5 -2.0 - - - - - - - -1940 -1960 -1980 -2000 -year -Convergence (R^) - - diff --git a/tests/testthat/_snaps/plot_functions/plot-seroprev-model-normallog.svg b/tests/testthat/_snaps/plot_functions/plot-seroprev-model-normallog.svg deleted file mode 100644 index 718153a9..00000000 --- a/tests/testthat/_snaps/plot_functions/plot-seroprev-model-normallog.svg +++ /dev/null @@ -1,350 +0,0 @@ - - - - - - - - - - - - - - - - - - - -seroprev_model_name: continuous_foi_normal_log -dataset: COL-035-93 -elpd: -71.93 -se: 7.11 -converged: Yes - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -0.00 -0.25 -0.50 -0.75 -1.00 - - - - - - - - - -0 -20 -40 -60 -Age -Sero-positivity - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -0.00 -0.05 -0.10 -0.15 - - - - - - - - -1940 -1960 -1980 -2000 -Year -Force-of-Infection - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -1.0 -1.5 -2.0 - - - - - - - -1940 -1960 -1980 -2000 -year -Convergence (R^) - - diff --git a/tests/testthat/_snaps/simdata_cases/case-a-group-foi-1234.svg b/tests/testthat/_snaps/simdata_cases/case-a-group-foi-1234.svg deleted file mode 100644 index 47facbd7..00000000 --- a/tests/testthat/_snaps/simdata_cases/case-a-group-foi-1234.svg +++ /dev/null @@ -1,371 +0,0 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -0.00 -0.01 -0.02 -0.03 - - - - - - - - - - -2000 -2010 -2020 -2030 -2040 -2050 -Year -Force-of-Infection -case_A_group_1234 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -0.00 -0.01 -0.02 -0.03 - - - - - - - - - - -2000 -2010 -2020 -2030 -2040 -2050 -Year -Force-of-Infection - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -0.00 -0.01 -0.02 -0.03 - - - - - - - - - - -2000 -2010 -2020 -2030 -2040 -2050 -Year -Force-of-Infection - - diff --git a/tests/testthat/_snaps/simdata_cases/case-a-group-models-1234.svg b/tests/testthat/_snaps/simdata_cases/case-a-group-models-1234.svg deleted file mode 100644 index 0b9419db..00000000 --- a/tests/testthat/_snaps/simdata_cases/case-a-group-models-1234.svg +++ /dev/null @@ -1,895 +0,0 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - -seroprev_model_name: constant_foi_bi -dataset: sim_foi -elpd: -20.79 -se: 1.66 -converged: Yes - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -0.00 -0.25 -0.50 -0.75 -1.00 - - - - - - - - - -0 -20 -40 -60 -Age -Sero-positivity - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -0.00 -0.01 -0.02 -0.03 - - - - - - - - - -2010 -2020 -2030 -2040 -2050 -Year -Force-of-Infection - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -1.0 -1.5 -2.0 - - - - - - - - -2010 -2020 -2030 -2040 -2050 -year -Convergence (R^) - - - - - - - - - - - - - - - - -seroprev_model_name: continuous_foi_normal_bi -dataset: sim_foi -elpd: -20.82 -se: 1.53 -converged: Yes - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -0.00 -0.25 -0.50 -0.75 -1.00 - - - - - - - - - -0 -20 -40 -60 -Age -Sero-positivity - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -0.00 -0.01 -0.02 -0.03 - - - - - - - - - -2010 -2020 -2030 -2040 -2050 -Year -Force-of-Infection - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -1.0 -1.5 -2.0 - - - - - - - - -2010 -2020 -2030 -2040 -2050 -year -Convergence (R^) - - - - - - - - - - - - - - - - -seroprev_model_name: continuous_foi_normal_log -dataset: sim_foi -elpd: -21.27 -se: 1.53 -converged: Yes - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -0.00 -0.25 -0.50 -0.75 -1.00 - - - - - - - - - -0 -20 -40 -60 -Age -Sero-positivity - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -0.00 -0.01 -0.02 -0.03 - - - - - - - - - -2010 -2020 -2030 -2040 -2050 -Year -Force-of-Infection - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -1.0 -1.5 -2.0 - - - - - - - - -2010 -2020 -2030 -2040 -2050 -year -Convergence (R^) - - diff --git a/tests/testthat/_snaps/simdata_cases/case-a-no-group-foi-1234.svg b/tests/testthat/_snaps/simdata_cases/case-a-no-group-foi-1234.svg deleted file mode 100644 index 86bf23a2..00000000 --- a/tests/testthat/_snaps/simdata_cases/case-a-no-group-foi-1234.svg +++ /dev/null @@ -1,371 +0,0 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -0.00 -0.01 -0.02 -0.03 - - - - - - - - - - -2000 -2010 -2020 -2030 -2040 -2050 -Year -Force-of-Infection -case_A_no_group_1234 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -0.00 -0.01 -0.02 -0.03 - - - - - - - - - - -2000 -2010 -2020 -2030 -2040 -2050 -Year -Force-of-Infection - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -0.00 -0.01 -0.02 -0.03 - - - - - - - - - - -2000 -2010 -2020 -2030 -2040 -2050 -Year -Force-of-Infection - - diff --git a/tests/testthat/_snaps/simdata_cases/case-a-no-group-models-1234.svg b/tests/testthat/_snaps/simdata_cases/case-a-no-group-models-1234.svg deleted file mode 100644 index aa1d25e1..00000000 --- a/tests/testthat/_snaps/simdata_cases/case-a-no-group-models-1234.svg +++ /dev/null @@ -1,922 +0,0 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - -seroprev_model_name: constant_foi_bi -dataset: sim_foi -elpd: -66.18 -se: 5.5 -converged: Yes - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -0.00 -0.25 -0.50 -0.75 -1.00 - - - - - - - - - -0 -20 -40 -60 -Age -Sero-positivity - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -0.00 -0.01 -0.02 -0.03 - - - - - - - - - - -2000 -2010 -2020 -2030 -2040 -2050 -Year -Force-of-Infection - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -1.0 -1.5 -2.0 - - - - - - - - - -2000 -2010 -2020 -2030 -2040 -2050 -year -Convergence (R^) - - - - - - - - - - - - - - - - -seroprev_model_name: continuous_foi_normal_bi -dataset: sim_foi -elpd: -66.45 -se: 5.63 -converged: Yes - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -0.00 -0.25 -0.50 -0.75 -1.00 - - - - - - - - - -0 -20 -40 -60 -Age -Sero-positivity - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -0.00 -0.01 -0.02 -0.03 - - - - - - - - - - -2000 -2010 -2020 -2030 -2040 -2050 -Year -Force-of-Infection - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -1.0 -1.5 -2.0 - - - - - - - - - -2000 -2010 -2020 -2030 -2040 -2050 -year -Convergence (R^) - - - - - - - - - - - - - - - - -seroprev_model_name: continuous_foi_normal_log -dataset: sim_foi -elpd: -66.8 -se: 5.59 -converged: Yes - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -0.00 -0.25 -0.50 -0.75 -1.00 - - - - - - - - - -0 -20 -40 -60 -Age -Sero-positivity - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -0.00 -0.01 -0.02 -0.03 - - - - - - - - - - -2000 -2010 -2020 -2030 -2040 -2050 -Year -Force-of-Infection - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -1.0 -1.5 -2.0 - - - - - - - - - -2000 -2010 -2020 -2030 -2040 -2050 -year -Convergence (R^) - - diff --git a/tests/testthat/_snaps/simdata_cases/case-b-group-foi-1234.svg b/tests/testthat/_snaps/simdata_cases/case-b-group-foi-1234.svg deleted file mode 100644 index 276dddd2..00000000 --- a/tests/testthat/_snaps/simdata_cases/case-b-group-foi-1234.svg +++ /dev/null @@ -1,368 +0,0 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -0.0 -0.1 -0.2 -0.3 - - - - - - - - - - -2000 -2010 -2020 -2030 -2040 -2050 -Year -Force-of-Infection -case_B_group_1234 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -0.0 -0.1 -0.2 -0.3 - - - - - - - - - - -2000 -2010 -2020 -2030 -2040 -2050 -Year -Force-of-Infection - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -0.0 -0.1 -0.2 -0.3 - - - - - - - - - - -2000 -2010 -2020 -2030 -2040 -2050 -Year -Force-of-Infection - - diff --git a/tests/testthat/_snaps/simdata_cases/case-b-group-models-1234.svg b/tests/testthat/_snaps/simdata_cases/case-b-group-models-1234.svg deleted file mode 100644 index 0bcc837d..00000000 --- a/tests/testthat/_snaps/simdata_cases/case-b-group-models-1234.svg +++ /dev/null @@ -1,892 +0,0 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - -seroprev_model_name: constant_foi_bi -dataset: sim_foi -elpd: -66.85 -se: 11.13 -converged: Yes - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -0.00 -0.25 -0.50 -0.75 -1.00 - - - - - - - - - -0 -20 -40 -60 -Age -Sero-positivity - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -0.0 -0.1 -0.2 -0.3 - - - - - - - - - -2010 -2020 -2030 -2040 -2050 -Year -Force-of-Infection - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -1.0 -1.5 -2.0 - - - - - - - - -2010 -2020 -2030 -2040 -2050 -year -Convergence (R^) - - - - - - - - - - - - - - - - -seroprev_model_name: continuous_foi_normal_bi -dataset: sim_foi -elpd: -23.35 -se: 4.93 -converged: Yes - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -0.00 -0.25 -0.50 -0.75 -1.00 - - - - - - - - - -0 -20 -40 -60 -Age -Sero-positivity - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -0.0 -0.1 -0.2 -0.3 - - - - - - - - - -2010 -2020 -2030 -2040 -2050 -Year -Force-of-Infection - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -1.0 -1.5 -2.0 - - - - - - - - -2010 -2020 -2030 -2040 -2050 -year -Convergence (R^) - - - - - - - - - - - - - - - - -seroprev_model_name: continuous_foi_normal_log -dataset: sim_foi -elpd: -14.73 -se: 3.82 -converged: Yes - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -0.00 -0.25 -0.50 -0.75 -1.00 - - - - - - - - - -0 -20 -40 -60 -Age -Sero-positivity - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -0.0 -0.1 -0.2 -0.3 - - - - - - - - - -2010 -2020 -2030 -2040 -2050 -Year -Force-of-Infection - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -1.0 -1.5 -2.0 - - - - - - - - -2010 -2020 -2030 -2040 -2050 -year -Convergence (R^) - - diff --git a/tests/testthat/_snaps/simdata_cases/case-b-no-group-foi-1234.svg b/tests/testthat/_snaps/simdata_cases/case-b-no-group-foi-1234.svg deleted file mode 100644 index fe0652bb..00000000 --- a/tests/testthat/_snaps/simdata_cases/case-b-no-group-foi-1234.svg +++ /dev/null @@ -1,368 +0,0 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -0.0 -0.1 -0.2 -0.3 - - - - - - - - - - -2000 -2010 -2020 -2030 -2040 -2050 -Year -Force-of-Infection -case_B_no_group_1234 - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -0.0 -0.1 -0.2 -0.3 - - - - - - - - - - -2000 -2010 -2020 -2030 -2040 -2050 -Year -Force-of-Infection - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -0.0 -0.1 -0.2 -0.3 - - - - - - - - - - -2000 -2010 -2020 -2030 -2040 -2050 -Year -Force-of-Infection - - diff --git a/tests/testthat/_snaps/simdata_cases/case-b-no-group-models-1234.svg b/tests/testthat/_snaps/simdata_cases/case-b-no-group-models-1234.svg deleted file mode 100644 index 69d37da0..00000000 --- a/tests/testthat/_snaps/simdata_cases/case-b-no-group-models-1234.svg +++ /dev/null @@ -1,919 +0,0 @@ - - - - - - - - - - - - - - - - - - - - - - - - - - - - -seroprev_model_name: constant_foi_bi -dataset: sim_foi -elpd: -78.52 -se: 5.81 -converged: Yes - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -0.00 -0.25 -0.50 -0.75 -1.00 - - - - - - - - - -0 -20 -40 -60 -Age -Sero-positivity - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -0.0 -0.1 -0.2 -0.3 - - - - - - - - - - -2000 -2010 -2020 -2030 -2040 -2050 -Year -Force-of-Infection - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -1.0 -1.5 -2.0 - - - - - - - - - -2000 -2010 -2020 -2030 -2040 -2050 -year -Convergence (R^) - - - - - - - - - - - - - - - - -seroprev_model_name: continuous_foi_normal_bi -dataset: sim_foi -elpd: -35.07 -se: 4.6 -converged: Yes - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -0.00 -0.25 -0.50 -0.75 -1.00 - - - - - - - - - -0 -20 -40 -60 -Age -Sero-positivity - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -0.0 -0.1 -0.2 -0.3 - - - - - - - - - - -2000 -2010 -2020 -2030 -2040 -2050 -Year -Force-of-Infection - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -1.0 -1.5 -2.0 - - - - - - - - - -2000 -2010 -2020 -2030 -2040 -2050 -year -Convergence (R^) - - - - - - - - - - - - - - - - -seroprev_model_name: continuous_foi_normal_log -dataset: sim_foi -elpd: -24.23 -se: 5.18 -converged: Yes - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -0.00 -0.25 -0.50 -0.75 -1.00 - - - - - - - - - -0 -20 -40 -60 -Age -Sero-positivity - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -0.0 -0.1 -0.2 -0.3 - - - - - - - - - - -2000 -2010 -2020 -2030 -2040 -2050 -Year -Force-of-Infection - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - - -1.0 -1.5 -2.0 - - - - - - - - - -2000 -2010 -2020 -2030 -2040 -2050 -year -Convergence (R^) - - diff --git a/tests/testthat/test_comparison.R b/tests/testthat/test_comparison.R deleted file mode 100644 index 43a779aa..00000000 --- a/tests/testthat/test_comparison.R +++ /dev/null @@ -1,63 +0,0 @@ -test_that("comparison", { - # So far we are skipping tests on these platforms until - # we find an efficient way to update rstan testthat snapshots on all of them - skip_on_os(c("windows", "mac")) - source("testing_utils.R") - set.seed(1234) # For reproducibility - - package <- "serofoi" - - # TODO For some reason it is not recognizing the global `serodata` variable, so we need to explicitly load it - serodata <- readRDS(testthat::test_path("extdata", "data.RDS")) - - data_test <- prepare_seroprev_data(serodata) - - model_0 <- run_seroprev_model( - seroprev_data = data_test, - seroprev_model_name = "constant_foi_bi", - n_iters = 1000 - ) - - model_1 <- run_seroprev_model( - seroprev_data = data_test, - seroprev_model_name = "continuous_foi_normal_bi", - n_iters = 1000 - ) - - model_2 <- run_seroprev_model( - seroprev_data = data_test, - seroprev_model_name = "continuous_foi_normal_log", - n_iters = 1000 - ) - - comp_table <- get_comparison_table( - model_objects_list = c( - m0 = model_0, - m1 = model_1, - m2 = model_2 - ) - ) - - column_comparation_functions <- list( - model = equal_exact(), - dataset = equal_exact(), - country = equal_exact(), - year = equal_exact(), - test = equal_exact(), - antibody = equal_exact(), - n_sample = equal_exact(), - n_agec = equal_exact(), - n_iter = equal_exact(), - elpd = equal_with_tolerance(), - se = equal_with_tolerance(), - converged = equal_exact(), - difference = equal_exact(), - diff_se = equal_with_tolerance(), - best_elpd = equal_exact(), - pvalue = equal_with_tolerance() - ) - - expect_similar_dataframes( - "comp_table", comp_table, column_comparation_functions - ) -}) diff --git a/tests/testthat/test_individual_models.R b/tests/testthat/test_individual_models.R deleted file mode 100644 index adcd1710..00000000 --- a/tests/testthat/test_individual_models.R +++ /dev/null @@ -1,58 +0,0 @@ -test_that("individual models", { - # So far we are skipping tests on these platforms until - # we find an efficient way to update rstan testthat snapshots on all of them - skip_on_os(c("windows", "mac")) - skip_on_ci() - - library(devtools) - library(dplyr) - library(vdiffr) - - #----- Read and prepare data - data_test_path <- testthat::test_path( - "extdata", "data.RDS" - ) - data_test <- readRDS(data_test_path) %>% prepare_seroprev_data(alpha = 0.05) - - #----- Test each model - model_0_object <- run_seroprev_model( - seroprev_data = data_test, - seroprev_model_name = "constant_foi_bi", - n_iters = 1000 - ) - - model_1_object <- run_seroprev_model( - seroprev_data = data_test, - seroprev_model_name = "continuous_foi_normal_bi", - n_iters = 1000 - ) - - model_2_object <- run_seroprev_model( - seroprev_data = data_test, - seroprev_model_name = "continuous_foi_normal_log", - n_iters = 1000 - ) - - #----- Generate all plots for each model - - model_0_plot <- plot_seroprev_model(model_0_object, size_text = 6) - vdiffr::expect_doppelganger("model_0_plot", model_0_plot) - model_1_plot <- plot_seroprev_model(model_1_object, size_text = 6) - vdiffr::expect_doppelganger("model_1_plot", model_1_plot) - model_2_plot <- plot_seroprev_model(model_2_object, size_text = 6) - vdiffr::expect_doppelganger("model_2_plot", model_2_plot) - - #----- Generate each individual plot - - sp_fitted_individual_plot <- plot_seroprev_fitted(model_2_object, size_text = 15) - vdiffr::expect_doppelganger("seroprev_fitted_individual_plot", sp_fitted_individual_plot) - - foi_individual_plot <- plot_foi(model_2_object, size_text = 15) - vdiffr::expect_doppelganger("foi_individual_plot", foi_individual_plot) - - rhats_individual_plot <- plot_rhats(model_2_object, size_text = 15) - vdiffr::expect_doppelganger("rhats_individual_plot", rhats_individual_plot) - - - # bayesplot::mcmc_trace(model_1_object$fit, pars="lambda0") -}) diff --git a/tests/testthat/test_models.R b/tests/testthat/test_models.R new file mode 100644 index 00000000..9b8d3a20 --- /dev/null +++ b/tests/testthat/test_models.R @@ -0,0 +1,48 @@ +test_that("individual models", { + # So far we are skipping tests on these platforms until + # we find an efficient way to update rstan testthat snapshots on all of them + skip_on_os(c("windows", "mac")) + skip_on_ci() + source("testing_utils.R") + set.seed(1234) # For reproducibility + + library(devtools) + library(dplyr) + library(vdiffr) + + #----- Read and prepare data + data_test_path <- testthat::test_path( + "extdata", "data.RDS" + ) + data("serodata") + data_test <- serodata %>% prepare_serodata(alpha = 0.05) + + #----- Plot raw data + data_test_plot <- plot_seroprev(data_test, size_text = 15) + vdiffr::expect_doppelganger("serodata_plot", data_test_plot) + + #----- Test each model + models_to_run <- c("constant_foi_bi", + "continuous_foi_normal_bi", + "continuous_foi_normal_log") + models_short_names <- c("constant", "normal", "normal_log") + + models_list <- lapply(models_to_run, run_seromodel, serodata = data_test, n_iters = 1000) + + #----- Generate plots for each model + i = 1 + for (model in models_list) { + model_plot <- plot_seromodel(model, size_text = 6) + vdiffr::expect_doppelganger(paste0(models_short_names[i], "_model_plot"), model_plot) + + model_seroprev_plot <- plot_seroprev_fitted(model, size_text = 15) + vdiffr::expect_doppelganger(paste0(models_short_names[i], "_sp_fitted_plot"), model_seroprev_plot) + + model_foi_plot <- plot_foi(model, size_text = 15) + vdiffr::expect_doppelganger(paste0(models_short_names[i], "_foi_plot"), model_foi_plot) + + model_rhats_plot <- plot_rhats(model, size_text = 15) + vdiffr::expect_doppelganger(paste0(models_short_names[i], "_rhats_plot"), model_rhats_plot) + i = i + 1 + } +}) diff --git a/tests/testthat/test_plot_functions.R b/tests/testthat/test_plot_functions.R deleted file mode 100644 index 8c3f66eb..00000000 --- a/tests/testthat/test_plot_functions.R +++ /dev/null @@ -1,68 +0,0 @@ -test_that("plot_seroprev_fitted", { - # So far we are skipping tests on these platforms until - # we find an efficient way to update rstan testthat snapshots on all of them - skip_on_os(c("windows", "mac")) - skip_on_ci() - print("*** Test info ****") - print(R.Version()) - cat("Interactive: ", interactive()) - library(devtools) - library(dplyr) - library(vdiffr) - set.seed(1234) # For reproducibility - data_test <- readRDS(testthat::test_path("extdata", "data.RDS")) %>% prepare_seroprev_data() - - actual_plot_seroprev <- plot_seroprev(data_test, size_text = 15) - - vdiffr::expect_doppelganger("plot_seroprev", actual_plot_seroprev) - - # Constant Model - model_object_constant <- run_seroprev_model( - seroprev_data = data_test, - seroprev_model_name = "constant_foi_bi", - n_iters = 1000 - ) - - plot_seroprev_model_constant <- plot_seroprev_model(model_object_constant, size_text = 6) - - vdiffr::expect_doppelganger("plot_seroprev_model_constant", plot_seroprev_model_constant) - - plot_seroprev_fitted_constant <- plot_seroprev_fitted(model_object_constant, size_text = 15) - - vdiffr::expect_doppelganger("plot_seroprev_fitted_constant", plot_seroprev_fitted_constant) - - # Normal Bi Model - model_object_normalbi <- run_seroprev_model( - seroprev_data = data_test, - seroprev_model_name = "continuous_foi_normal_bi", - n_iters = 1000 - ) - - plot_seroprev_model_normalbi <- plot_seroprev_model(model_object_normalbi, size_text = 6) - - vdiffr::expect_doppelganger("plot_seroprev_model_normalbi", plot_seroprev_model_normalbi) - - plot_seroprev_fitted_normalbi <- plot_seroprev_fitted(model_object_normalbi, size_text = 15) - - vdiffr::expect_doppelganger("plot_seroprev_fitted_normalbi", plot_seroprev_fitted_normalbi) - - # Normal Log Model - model_object_normallog <- run_seroprev_model( - seroprev_data = data_test, - seroprev_model_name = "continuous_foi_normal_log", - n_iters = 1000 - ) - - plot_seroprev_model_normallog <- plot_seroprev_model(model_object_normallog, size_text = 6) - - vdiffr::expect_doppelganger("plot_seroprev_model_normallog", plot_seroprev_model_normallog) - - plot_seroprev_fitted_normallog <- plot_seroprev_fitted(model_object_normallog, size_text = 15) - - vdiffr::expect_doppelganger("plot_seroprev_fitted_normallog", plot_seroprev_fitted_normallog) - - # Models Comparison Plot - plot_arrange_models <- plot_seroprev_models_grid(plot_seroprev_model_constant, plot_seroprev_model_normalbi, plot_seroprev_model_normallog, n_col = 3) - - vdiffr::expect_doppelganger("plot_arrange_models", plot_arrange_models) -}) diff --git a/tests/testthat/test_simdata_cases.R b/tests/testthat/test_simdata_cases.R deleted file mode 100644 index ad707716..00000000 --- a/tests/testthat/test_simdata_cases.R +++ /dev/null @@ -1,130 +0,0 @@ -test_that("simulated data", { - # So far we are skipping tests on these platforms until - # we find an efficient way to update rstan testthat snapshots on all of them - skip_on_os(c("windows", "mac")) - skip_on_ci() - - library(dplyr) - library(serofoi) - library(pracma) - - seed = 1234 - size_age_class = 5 - #----- Test foi case A - sim_foi <- rep(0.02, 50) - case_label <- "case_A_" - max_lambda <- 0.035 - - #----- Test foi case B - # no_transm <- 0.0000000001 - # sim_foi <- c(rep(0.2, 25), rep(0.1, 10), rep(no_transm, 15)) - # case_label <- "case_B_" - # max_lambda <- 0.3 - - #----- Results paths - data_path_no_grouped <- testthat::test_path( - "extdata", paste0(case_label, "sim_data_no_grouped.csv") - ) - - data_path_grouped <- testthat::test_path( - "extdata", paste0(case_label, "sim_data_grouped.csv") - ) - - #----- Data simulation - sim_data <- generate_sim_data(foi = sim_foi, - size_age_class = size_age_class, - tsur = 2050, - birth_year_min = 2000, - survey_label = 'sim_foi', - seed = seed) - - sim_data <- sim_data %>% mutate(age_min = age_mean_f, age_max = age_mean_f) - write.csv(sim_data, data_path_no_grouped, row.names = FALSE) - - model_constant <- run_seroprev_model(seroprev_data = sim_data, - seroprev_model_name = "constant_foi_bi") - - model_normal <- run_seroprev_model(seroprev_data = sim_data, - seroprev_model_name = "continuous_foi_normal_bi") - - model_log <- run_seroprev_model(seroprev_data = sim_data, - seroprev_model_name = "continuous_foi_normal_log") - - model_plot_constant <- plot_seroprev_model(model_constant, size_text = 6, , max_lambda = max_lambda) - model_plot_normal <- plot_seroprev_model(model_normal, size_text = 6, , max_lambda = max_lambda) - model_plot_log <- plot_seroprev_model(model_log, size_text = 6, max_lambda = max_lambda) - model_plot_arrange <- plot_seroprev_models_grid(model_plot_constant, - model_plot_normal, - model_plot_log, - size_text = 6, - n_col = 3, - n_row = 1) - vdiffr::expect_doppelganger(paste0(case_label, "no_group_models_", seed), model_plot_arrange) - - sim_foi_plot_constant <- plot_foi(model_constant, size_text = 10, max_lambda = max_lambda) + - ggplot2::geom_point(data = data.frame(year = model_constant$exposure_years, - foi = sim_foi), - ggplot2::aes(year, foi)) + - ggplot2::ggtitle(paste0(case_label, "no_group_", seed)) - - sim_foi_plot_normal <- plot_foi(model_normal, size_text = 10, max_lambda = max_lambda) + - ggplot2::geom_point(data = data.frame(year = model_normal$exposure_years, - foi = sim_foi), - ggplot2::aes(year, foi)) - sim_foi_plot_log <- plot_foi(model_log, size_text = 10, max_lambda = max_lambda) + - ggplot2::geom_point(data = data.frame(year = model_log$exposure_years, - foi = sim_foi), - ggplot2::aes(year, foi)) - sim_foi_plot_arrange <- plot_seroprev_models_grid(sim_foi_plot_constant, - sim_foi_plot_normal, - sim_foi_plot_log, - n_col = 1, n_row = 3) - - vdiffr::expect_doppelganger(paste0(case_label, "no_group_foi_", seed), sim_foi_plot_arrange) - - - sim_data_grouped <- group_sim_data(sim_data = sim_data, - foi = sim_foi, - size_age_class = size_age_class, - tsur = 2050, - birth_year_min = 2000, - survey_label = 'sim_foi') - write.csv(sim_data_grouped, data_path_grouped, row.names = FALSE) - - model_grouped_constant <- run_seroprev_model(seroprev_data = sim_data_grouped, - seroprev_model_name = "constant_foi_bi") - model_grouped_normal <- run_seroprev_model(seroprev_data = sim_data_grouped, - seroprev_model_name = "continuous_foi_normal_bi") - model_grouped_log <- run_seroprev_model(seroprev_data = sim_data_grouped, - seroprev_model_name = "continuous_foi_normal_log") - - model_grouped_constant_plot <- plot_seroprev_model(model_grouped_constant, size_text = 6, max_lambda = max_lambda) - model_grouped_normal_plot <- plot_seroprev_model(model_grouped_normal, size_text = 6, max_lambda = max_lambda) - model_grouped_log_plot <- plot_seroprev_model(model_grouped_log, size_text = 6, max_lambda = max_lambda) - model_grouped_plot_arrange <- plot_seroprev_models_grid(model_grouped_constant_plot, - model_grouped_normal_plot, - model_grouped_log_plot, - n_col = 3, n_row = 1) - vdiffr::expect_doppelganger(paste0(case_label, "group_models_", seed), model_grouped_plot_arrange) - rm(list = ls(pat = "*_plot")) - - sim_foi_grouped_constant_plot <- plot_foi(model_grouped_constant, size_text = 10, max_lambda = max_lambda) + - ggplot2::geom_point(data = data.frame(year = model_constant$exposure_years, - foi = sim_foi), - ggplot2::aes(year, foi)) + - ggplot2::ggtitle(paste0(case_label, "group_", seed)) - sim_foi_grouped_normal_plot <- plot_foi(model_grouped_normal, size_text = 10, max_lambda = max_lambda) + - ggplot2::geom_point(data = data.frame(year = model_normal$exposure_years, - foi = sim_foi), - ggplot2::aes(year, foi)) - sim_foi_grouped_log_plot <- plot_foi(model_grouped_log, size_text = 10, max_lambda = max_lambda) + - ggplot2::geom_point(data = data.frame(year = model_log$exposure_years, - foi = sim_foi), - ggplot2::aes(year, foi)) - sim_foi_grouped_plot_arrange <- plot_seroprev_models_grid(sim_foi_grouped_constant_plot, - sim_foi_grouped_normal_plot, - sim_foi_grouped_log_plot, - n_col = 1, n_row = 3) - vdiffr::expect_doppelganger(paste0(case_label, "group_foi_", seed), sim_foi_grouped_plot_arrange) -}) -# TODO: solve error 'address 0x18, cause 'memory not mapped' when testing. diff --git a/vignettes/serofoi.Rmd b/vignettes/serofoi.Rmd index 7d8b042a..3c7ddc89 100644 --- a/vignettes/serofoi.Rmd +++ b/vignettes/serofoi.Rmd @@ -20,5 +20,5 @@ library(serofoi) # Prepare the data for using serofoi ```{r} -data_test <- prepare_seroprev_data(serodata) +data_test <- prepare_serodata(serodata) ``` diff --git a/vignettes/simulated_data.Rmd b/vignettes/simulated_data.Rmd index 620d6da1..8a2c3238 100644 --- a/vignettes/simulated_data.Rmd +++ b/vignettes/simulated_data.Rmd @@ -4,7 +4,7 @@ output: rmarkdown::html_vignette bibliography: references.bib link-citations: true vignette: > - %\VignetteIndexEntry{simdata} + %\VignetteIndexEntry{Simulated data} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- @@ -20,7 +20,7 @@ knitr::opts_chunk$set( library(serofoi) library(tidyverse) rownames(serodata) <- NULL -data_test <- prepare_seroprev_data(serodata) +data_test <- prepare_serodata(serodata) ``` @@ -74,8 +74,8 @@ $$ P(a,t) = 1−exp(−\lambda{a})$$ The corresponding code for the constant-FoI model with binomial distribution is: ```{r model_1, include = FALSE, echo=TRUE, errors = FALSE, warning = FALSE, message = FALSE } -model_1 <- run_seroprev_model(seroprev_data = data_test, - seroprev_model_name = "constant_foi_bi", +model_1 <- run_seromodel(serodata = data_test, + seromodel_name = "constant_foi_bi", n_iters = 500, n_thin = 2) @@ -94,15 +94,15 @@ Therefore, a serosurvey completed at time $t$ and including ages $a$ from {$age_ Currently, two models within *serofoi* allow time-varying FoI, with corresponding code: ```{r model_2, include = FALSE, echo = TRUE, errors = FALSE, warning = FALSE, message = FALSE } -model_2 <- run_seroprev_model(seroprev_data = data_test, - seroprev_model_name = "continuous_foi_normal_bi", +model_2 <- run_seromodel(serodata = data_test, + seromodel_name = "continuous_foi_normal_bi", n_iters = 1500, n_thin = 2) ``` ```{r model_3, include = FALSE, echo = TRUE, errors = FALSE, warning = FALSE, message = FALSE } -model_3 <- run_seroprev_model(seroprev_data = data_test, - seroprev_model_name = "continuous_foi_normal_log", +model_3 <- run_seromodel(serodata = data_test, + seromodel_name = "continuous_foi_normal_log", n_iters = 1500, n_thin = 2) ``` @@ -125,15 +125,15 @@ After preparing the data we can run two models (constant and time varying) avail ```{r three models, include = TRUE, errors = FALSE, warning = FALSE, message=FALSE } -chagas2012p <- prepare_seroprev_data(chagas2012) +chagas2012p <- prepare_serodata(chagas2012) -m1_cha <- run_seroprev_model(seroprev_data = chagas2012p, - seroprev_model_name = "constant_foi_bi", +m1_cha <- run_seromodel(serodata = chagas2012p, + seromodel_name = "constant_foi_bi", n_iters = 500, n_thin = 2) -m2_cha <- run_seroprev_model(seroprev_data = chagas2012p, - seroprev_model_name = "continuous_foi_normal_bi", +m2_cha <- run_seromodel(serodata = chagas2012p, + seromodel_name = "continuous_foi_normal_bi", n_iters = 1500, n_thin = 2) @@ -142,10 +142,10 @@ m2_cha <- run_seroprev_model(seroprev_data = chagas2012p, Now, we can plot the results of the two models to compare (Figure 1). As shown in @Cucunubá2017, interventions for Chagas control have been ongoing from the 1980s in Colombia having heterogeneous impact depending on the type of population. For this serosurvey, which is from tradittional indigenous rural area, the serofoi models are able to detect a modest still relevant slow decreasing trend consistent (model 2) as slightly better supported than the constant model. Notice that model 3 does not converge despite a high number of iterations. ```{r, model_comp_cha, include=TRUE, echo=TRUE, fig.width = 8, fig.height=10, warning=FALSE, message=FALSE, message=FALSE} -p1_cha <- plot_seroprev_model(m1_cha, size_text = 12, max_lambda = 0.02) -p2_cha <- plot_seroprev_model(m2_cha, size_text = 12, max_lambda = 0.02) +p1_cha <- plot_seromodel(m1_cha, size_text = 12, max_lambda = 0.02) +p2_cha <- plot_seromodel(m2_cha, size_text = 12, max_lambda = 0.02) -plot_seroprev_models_grid(p1_cha, p2_cha, n_row = 1, n_col = 2) +plot_seromodels_grid(p1_cha, p2_cha, n_row = 1, n_col = 2) ``` ## Case study 2. Hidden Alphaviruses epidemics in Panama @@ -153,28 +153,27 @@ plot_seroprev_models_grid(p1_cha, p2_cha, n_row = 1, n_col = 2) As shown in [@Carrera2020], hidden epidemic and endemic transmission of alphaviruses in Eastern Panama have been around for decades. From this paper we use a dataset measuring IgG antibodies againts Venezuelan Equine Encephalitis Virus (VEEV) in a rural village in Panamá in 2017. This dataset, `veev2017` is included in serofoi. ```{r data_veev, include = TRUE, errors = FALSE, warning = FALSE, message = FALSE } - data("veev2012") head(veev2012, 5) ``` -```{r models_veee, include = TRUE, errors = FALSE, warning = FALSE, message=FALSE } +```{r models_veev, include = TRUE, errors = FALSE, warning = FALSE, message=FALSE } -veev2012p <- prepare_seroprev_data(veev2012) +veev2012p <- prepare_serodata(veev2012) -m1_veev <- run_seroprev_model(seroprev_data = veev2012p, - seroprev_model_name = "constant_foi_bi", +m1_veev <- run_seromodel(serodata = veev2012p, + seromodel_name = "constant_foi_bi", n_iters = 500, n_thin = 2) -m2_veev <- run_seroprev_model(seroprev_data = veev2012p, - seroprev_model_name = "continuous_foi_normal_bi", +m2_veev <- run_seromodel(serodata = veev2012p, + seromodel_name = "continuous_foi_normal_bi", n_iters = 500, n_thin = 2) -m3_veev <- run_seroprev_model(seroprev_data = veev2012p, - seroprev_model_name = "continuous_foi_normal_log", +m3_veev <- run_seromodel(serodata = veev2012p, + seromodel_name = "continuous_foi_normal_log", n_iters = 500, n_thin = 2) ``` @@ -182,11 +181,11 @@ m3_veev <- run_seroprev_model(seroprev_data = veev2012p, Now, we can plot the results of the three models to compare. On Figure 3, we can observe a large increase in the estimated FOI. As suggested in @Carrera2020, an important increase in the transmission of VEEV in this region is inferred. Although the three fitted models converge well, we see much larger support for a more sudden increase in recent years with highest FoI values at 0.25 (consistent with an epidemic). ```{r, model_comp_veev, include=TRUE, echo=TRUE, fig.width = 10, fig.height=10, warning=FALSE, message=FALSE, message=FALSE} -p1_veev <- plot_seroprev_model(m1_veev, size_text = 10, max_lambda = 0.6) -p2_veev <- plot_seroprev_model(m2_veev, size_text = 10, max_lambda = 0.6) -p3_veev <- plot_seroprev_model(m3_veev, size_text = 10, max_lambda = 0.6) +p1_veev <- plot_seromodel(m1_veev, size_text = 10, max_lambda = 0.6) +p2_veev <- plot_seromodel(m2_veev, size_text = 10, max_lambda = 0.6) +p3_veev <- plot_seromodel(m3_veev, size_text = 10, max_lambda = 0.6) -plot_seroprev_models_grid(p1_veev, p2_veev, p3_veev, n_row = 1, n_col = 3) +plot_seromodels_grid(p1_veev, p2_veev, p3_veev, n_row = 1, n_col = 3) ``` @@ -197,29 +196,27 @@ Figure 3. Results of fitted models, Force-of-Infection estimates and convergence Chikungunya outbreaks ocurred rapidly after the introduction of the virus to Brazil in 2013-2014. We use the dataset from [@dias2018] that conducts a population-based study through household interviews and serologic surveys (measuring IgG antibodies against Chikungunya virus) in Bahia, Brazil during October-December 2015, right after a large CHIKV epidemic that occurred in that area. ```{r data_chik, include = TRUE, errors = FALSE, warning = FALSE, message = FALSE } - data("chik2015") - head(chik2015, 5) ``` ```{r models_chik, include = TRUE, errors = FALSE, warning = FALSE, message=FALSE } -chik2015p <- prepare_seroprev_data(chik2015) +chik2015p <- prepare_serodata(chik2015) -mod1_chik <- run_seroprev_model(seroprev_data = chik2015p, - seroprev_model_name = "constant_foi_bi", +mod1_chik <- run_seromodel(serodata = chik2015p, + seromodel_name = "constant_foi_bi", n_iters = 1000, n_thin = 2) -mod2_chik <- run_seroprev_model(seroprev_data = chik2015p, - seroprev_model_name = "continuous_foi_normal_bi", +mod2_chik <- run_seromodel(serodata = chik2015p, + seromodel_name = "continuous_foi_normal_bi", n_iters = 1500, n_thin = 2) -mod3_chik <- run_seroprev_model(seroprev_data = chik2015p, - seroprev_model_name = "continuous_foi_normal_log", +mod3_chik <- run_seromodel(serodata = chik2015p, + seromodel_name = "continuous_foi_normal_log", n_iters = 1500, n_thin = 2) ``` @@ -228,11 +225,11 @@ In Figure 4, we can observe the comparison between the three serofoi models. Her ```{r plot_chik, include=TRUE, echo=TRUE, fig.width = 10, fig.height = 10, warning = FALSE, message = FALSE} -p1_chik <- plot_seroprev_model(mod1_chik, size_text = 10, max_lambda = 0.08) -p2_chik <- plot_seroprev_model(mod2_chik, size_text = 10, max_lambda = 0.08) -p3_chik <- plot_seroprev_model(mod3_chik, size_text = 10, max_lambda = 0.08) +p1_chik <- plot_seromodel(mod1_chik, size_text = 10, max_lambda = 0.08) +p2_chik <- plot_seromodel(mod2_chik, size_text = 10, max_lambda = 0.08) +p3_chik <- plot_seromodel(mod3_chik, size_text = 10, max_lambda = 0.08) -plot_seroprev_models_grid(p1_chik, p2_chik, p3_chik, n_row = 1, n_col = 3) +plot_seromodels_grid(p1_chik, p2_chik, p3_chik, n_row = 1, n_col = 3) ``` diff --git a/vignettes/use_cases.Rmd b/vignettes/use_cases.Rmd index 2455f81f..3cba71eb 100644 --- a/vignettes/use_cases.Rmd +++ b/vignettes/use_cases.Rmd @@ -4,7 +4,7 @@ output: rmarkdown::html_vignette bibliography: references.bib link-citations: true vignette: > - %\VignetteIndexEntry{use_cases} + %\VignetteIndexEntry{Use cases} %\VignetteEngine{knitr::rmarkdown} %\VignetteEncoding{UTF-8} --- @@ -20,7 +20,7 @@ knitr::opts_chunk$set( library(serofoi) library(tidyverse) rownames(serodata) <- NULL -data_test <- prepare_seroprev_data(serodata) +data_test <- prepare_serodata(serodata) ``` @@ -74,10 +74,10 @@ $$ P(a,t) = 1−exp(−\lambda{a})$$ The corresponding code for the constant-FoI model with binomial distribution is: ```{r model_1, include = FALSE, echo=TRUE, errors = FALSE, warning = FALSE, message = FALSE } -model_1 <- run_seroprev_model(seroprev_data = data_test, - seroprev_model_name = "constant_foi_bi", - n_iters = 500, - n_thin = 2) +model_1 <- run_seromodel(serodata = data_test, + seromodel_name = "constant_foi_bi", + n_iters = 500, + n_thin = 2) ``` @@ -94,15 +94,17 @@ Therefore, a serosurvey completed at time $t$ and including ages $a$ from {$age_ Currently, two models within *serofoi* allow time-varying FoI, with corresponding code: ```{r model_2, include = FALSE, echo = TRUE, errors = FALSE, warning = FALSE, message = FALSE } -model_2 <- run_seroprev_model(seroprev_data = data_test, - seroprev_model_name = "continuous_foi_normal_bi", - n_iters = 1500, - n_thin = 2) +model_2 <- run_seromodel(serodata = data_test, + seromodel_name = "continuous_foi_normal_bi", + n_iters = 1500, + n_thin = 2) +``` -model_3 <- run_seroprev_model(seroprev_data = data_test, - seroprev_model_name = "continuous_foi_normal_log", - n_iters = 1500, - n_thin = 2) +```{r model_3, include = FALSE, echo = TRUE, errors = FALSE, warning = FALSE, message = FALSE } +model_3 <- run_seromodel(serodata = data_test, + seromodel_name = "continuous_foi_normal_log", + n_iters = 1500, + n_thin = 2) ``` ## Fitting process @@ -114,7 +116,7 @@ Models are implemented in **RStan** using the No-U-Turn sampler, a type of Hamil Based on the data and analysis shown in [@Cucunubá2017], we use one of the datasets for measuring the sero-prevalence of IgG antibodies against *Trypanosoma cruzi* infection in rural area of Colombia in 2012. The dataset is part of the serofoi package as `chagas2012`. ```{r data_chagas, include = TRUE, errors = FALSE, warning = FALSE, message = FALSE } -chagas2012 <- readRDS("../data/chagas2012.RDS") # This can be removed once data is incorporated +data("chagas2012") head(chagas2012, 5) ``` @@ -123,15 +125,15 @@ After preparing the data we can run two models (constant and time varying) avail ```{r three models, include = TRUE, errors = FALSE, warning = FALSE, message=FALSE } -chagas2012p <- prepare_seroprev_data(chagas2012) +chagas2012p <- prepare_serodata(chagas2012) -m1_cha <- run_seroprev_model(seroprev_data = chagas2012p, - seroprev_model_name = "constant_foi_bi", +m1_cha <- run_seromodel(serodata = chagas2012p, + seromodel_name = "constant_foi_bi", n_iters = 500, n_thin = 2) -m2_cha <- run_seroprev_model(seroprev_data = chagas2012p, - seroprev_model_name = "continuous_foi_normal_bi", +m2_cha <- run_seromodel(serodata = chagas2012p, + seromodel_name = "continuous_foi_normal_bi", n_iters = 1500, n_thin = 2) @@ -140,10 +142,10 @@ m2_cha <- run_seroprev_model(seroprev_data = chagas2012p, Now, we can plot the results of the two models to compare (Figure 1). As shown in @Cucunubá2017, interventions for Chagas control have been ongoing from the 1980s in Colombia having heterogeneous impact depending on the type of population. For this serosurvey, which is from tradittional indigenous rural area, the serofoi models are able to detect a modest still relevant slow decreasing trend consistent (model 2) as slightly better supported than the constant model. Notice that model 3 does not converge despite a high number of iterations. ```{r, model_comp_cha, include=TRUE, echo=TRUE, fig.width = 8, fig.height=10, warning=FALSE, message=FALSE, message=FALSE} -p1_cha <- plot_seroprev_model(m1_cha, size_text = 12, max_lambda = 0.02) -p2_cha <- plot_seroprev_model(m2_cha, size_text = 12, max_lambda = 0.02) +p1_cha <- plot_seromodel(m1_cha, size_text = 12, max_lambda = 0.02) +p2_cha <- plot_seromodel(m2_cha, size_text = 12, max_lambda = 0.02) -plot_seroprev_models_grid(p1_cha, p2_cha, n_row = 1, n_col = 2) +plot_seromodels_grid(p1_cha, p2_cha, n_row = 1, n_col = 2) ``` ## Case study 2. Hidden Alphaviruses epidemics in Panama @@ -151,28 +153,27 @@ plot_seroprev_models_grid(p1_cha, p2_cha, n_row = 1, n_col = 2) As shown in [@Carrera2020], hidden epidemic and endemic transmission of alphaviruses in Eastern Panama have been around for decades. From this paper we use a dataset measuring IgG antibodies againts Venezuelan Equine Encephalitis Virus (VEEV) in a rural village in Panamá in 2017. This dataset, `veev2017` is included in serofoi. ```{r data_veev, include = TRUE, errors = FALSE, warning = FALSE, message = FALSE } - -veev2012 <- readRDS("../data/veev2012.RDS") +data("veev2012") head(veev2012, 5) ``` -```{r models_veee, include = TRUE, errors = FALSE, warning = FALSE, message=FALSE } +```{r models_veev, include = TRUE, errors = FALSE, warning = FALSE, message=FALSE } -veev2012p <- prepare_seroprev_data(veev2012) +veev2012p <- prepare_serodata(veev2012) -m1_veev <- run_seroprev_model(seroprev_data = veev2012p, - seroprev_model_name = "constant_foi_bi", +m1_veev <- run_seromodel(serodata = veev2012p, + seromodel_name = "constant_foi_bi", n_iters = 500, n_thin = 2) -m2_veev <- run_seroprev_model(seroprev_data = veev2012p, - seroprev_model_name = "continuous_foi_normal_bi", +m2_veev <- run_seromodel(serodata = veev2012p, + seromodel_name = "continuous_foi_normal_bi", n_iters = 500, n_thin = 2) -m3_veev <- run_seroprev_model(seroprev_data = veev2012p, - seroprev_model_name = "continuous_foi_normal_log", +m3_veev <- run_seromodel(serodata = veev2012p, + seromodel_name = "continuous_foi_normal_log", n_iters = 500, n_thin = 2) ``` @@ -180,11 +181,11 @@ m3_veev <- run_seroprev_model(seroprev_data = veev2012p, Now, we can plot the results of the three models to compare. On Figure 3, we can observe a large increase in the estimated FOI. As suggested in @Carrera2020, an important increase in the transmission of VEEV in this region is inferred. Although the three fitted models converge well, we see much larger support for a more sudden increase in recent years with highest FoI values at 0.25 (consistent with an epidemic). ```{r, model_comp_veev, include=TRUE, echo=TRUE, fig.width = 10, fig.height=10, warning=FALSE, message=FALSE, message=FALSE} -p1_veev <- plot_seroprev_model(m1_veev, size_text = 10, max_lambda = 0.6) -p2_veev <- plot_seroprev_model(m2_veev, size_text = 10, max_lambda = 0.6) -p3_veev <- plot_seroprev_model(m3_veev, size_text = 10, max_lambda = 0.6) +p1_veev <- plot_seromodel(m1_veev, size_text = 10, max_lambda = 0.6) +p2_veev <- plot_seromodel(m2_veev, size_text = 10, max_lambda = 0.6) +p3_veev <- plot_seromodel(m3_veev, size_text = 10, max_lambda = 0.6) -plot_seroprev_models_grid(p1_veev, p2_veev, p3_veev, n_row = 1, n_col = 3) +plot_seromodels_grid(p1_veev, p2_veev, p3_veev, n_row = 1, n_col = 3) ``` @@ -195,28 +196,27 @@ Figure 3. Results of fitted models, Force-of-Infection estimates and convergence Chikungunya outbreaks ocurred rapidly after the introduction of the virus to Brazil in 2013-2014. We use the dataset from [@dias2018] that conducts a population-based study through household interviews and serologic surveys (measuring IgG antibodies against Chikungunya virus) in Bahia, Brazil during October-December 2015, right after a large CHIKV epidemic that occurred in that area. ```{r data_chik, include = TRUE, errors = FALSE, warning = FALSE, message = FALSE } - -chik2015 <- readRDS("../data/chik2015.RDS") +data("chik2015") head(chik2015, 5) ``` ```{r models_chik, include = TRUE, errors = FALSE, warning = FALSE, message=FALSE } -chik2015p <- prepare_seroprev_data(chik2015) +chik2015p <- prepare_serodata(chik2015) -mod1_chik <- run_seroprev_model(seroprev_data = chik2015p, - seroprev_model_name = "constant_foi_bi", +mod1_chik <- run_seromodel(serodata = chik2015p, + seromodel_name = "constant_foi_bi", n_iters = 1000, n_thin = 2) -mod2_chik <- run_seroprev_model(seroprev_data = chik2015p, - seroprev_model_name = "continuous_foi_normal_bi", +mod2_chik <- run_seromodel(serodata = chik2015p, + seromodel_name = "continuous_foi_normal_bi", n_iters = 1500, n_thin = 2) -mod3_chik <- run_seroprev_model(seroprev_data = chik2015p, - seroprev_model_name = "continuous_foi_normal_log", +mod3_chik <- run_seromodel(serodata = chik2015p, + seromodel_name = "continuous_foi_normal_log", n_iters = 1500, n_thin = 2) ``` @@ -224,11 +224,11 @@ mod3_chik <- run_seroprev_model(seroprev_data = chik2015p, In Figure 4, we can observe the comparison between the three serofoi models. Here serofoi shows strong statistical support for a sudden increase in the transmission of CHIKV close to the year of the serosurvey (2015). The exact year is not possible to estimate, mainly given the data used is largely aggregated by 20-years age groups. Despite that, the results are consistent with the empirical evidence shown by [@dias2018] with both interviews, and IgM testing. ```{r, plot_chik, include=TRUE, echo=TRUE, fig.width = 10, fig.height=10, warning=FALSE, message=FALSE, message=FALSE} -p1_chik <- plot_seroprev_model(mod1_chik, size_text = 10, max_lambda = 0.08) -p2_chik <- plot_seroprev_model(mod2_chik, size_text = 10, max_lambda = 0.08) -p3_chik <- plot_seroprev_model(mod3_chik, size_text = 10, max_lambda = 0.08) +p1_chik <- plot_seromodel(mod1_chik, size_text = 10, max_lambda = 0.08) +p2_chik <- plot_seromodel(mod2_chik, size_text = 10, max_lambda = 0.08) +p3_chik <- plot_seromodel(mod3_chik, size_text = 10, max_lambda = 0.08) -plot_seroprev_models_grid(p1_chik, p2_chik, p3_chik, n_row = 1, n_col = 3) +plot_seromodels_grid(p1_chik, p2_chik, p3_chik, n_row = 1, n_col = 3) ```