diff --git a/DESCRIPTION b/DESCRIPTION index a6368f98..e3aa12d1 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -58,6 +58,6 @@ License: GPL (>=2) LazyLoad: yes Encoding: UTF-8 Roxygen: list(markdown = TRUE) -RoxygenNote: 7.3.3 VignetteBuilder: knitr Remotes: reconverse/incidence2/pkg +Config/roxygen2/version: 8.0.0 diff --git a/R/EpiEstim-package.R b/R/EpiEstim-package.R index 695a2efd..999f1ae9 100644 --- a/R/EpiEstim-package.R +++ b/R/EpiEstim-package.R @@ -8,10 +8,10 @@ #' \code{system.file("epiestimpapers.bib", package = "EpiEstim")}. #' #' @importFrom ggplot2 last_plot ggplot aes geom_step ggtitle -#' geom_ribbon geom_line xlab ylab xlim geom_hline ylim geom_histogram -#' scale_colour_manual scale_fill_manual scale_linetype_manual lims theme -#' margin element_rect theme_light %+replace% element_blank element_line -#' element_text scale_y_continuous +#' @importFrom ggplot2 geom_ribbon geom_line xlab ylab xlim geom_hline ylim geom_histogram +#' @importFrom ggplot2 scale_colour_manual scale_fill_manual scale_linetype_manual lims theme +#' @importFrom ggplot2 margin element_rect theme_light %+replace% element_blank element_line +#' @importFrom ggplot2 element_text scale_y_continuous #' #' @importFrom rlang .data #' diff --git a/R/gibbs_draws.R b/R/gibbs_draws.R index 9bffd8fd..982929ce 100644 --- a/R/gibbs_draws.R +++ b/R/gibbs_draws.R @@ -183,15 +183,6 @@ compute_lambda <- function(incid, si_distr) { msg2 <- "Use function 'process_I_multivariant' first" stop(msg1, msg2) } - if (any(si_distr[1,] != 0)){ - stop("Values in the first row of si_distr must be 0") - } - if (any(abs(colSums(si_distr) - 1) > 0.01)) { # allow tolerance - stop("The sum of each column in si_distr should be equal to 1") - } - if (any(si_distr < 0)){ - stop("si_distr must be >=0") - } lambda <- array(NA, dim = dim(incid$local)) for(l in seq_len(dim(incid$local)[2])) { for(v in seq_len(dim(incid$local)[3])) { @@ -242,6 +233,11 @@ compute_lambda <- function(incid, si_distr) { #' #' @param seed a numeric value used to fix the random seed #' +#' @param validate_inputs a boolean value indicating whether to validate the +#' inputs; default value is FALSE as inputs are validated in +#' \code{estimate_advantage()}, where this function is called. Can be set to TRUE +#' if this function is called directly for debugging or testing. +#' #' @return A value or vector of values for epsilon for each non reference #' pathogen/strain/variant, drawn from the marginal posterior distribution #' @@ -263,26 +259,17 @@ compute_lambda <- function(incid, si_distr) { #' R <- matrix(1, nrow = T, ncol = n_loc) #' R[1, ] <- NA # no estimates of R on first time step #' draw_epsilon(R, incid$local, lambda, priors, seed = 1) - draw_epsilon <- function(R, incid, lambda, priors, shape_epsilon = NULL, t_min = 2L, t_max = nrow(incid), - seed = NULL) { - if (!is.integer(t_min) || !is.integer(t_max)){ - stop("t_min and t_max must be integers") - } - if (t_min < 2 || t_max < 2){ - stop("t_min and t_max must be >=2") - } - if(t_min > nrow(incid) || t_max > nrow(incid)){ - stop("t_min and t_max must be <= nrow(incid)") + seed = NULL, validate_inputs = FALSE) { + if (validate_inputs) { + eps_args <- as.list(environment()) + do.call(check_estimate_advantage_inputs, args = eps_args) } if(any(R[!is.na(R)] < 0)) { stop("R must be >= 0") } - if (!is.null(seed) && !is.numeric(seed)){ - stop("seed must be numeric") - } if (!is.null(seed)) set.seed(seed) t <- seq(t_min, t_max, 1) if (is.null(shape_epsilon)) { @@ -333,6 +320,8 @@ draw_epsilon <- function(R, incid, lambda, priors, #' #' @param seed a numeric value used to fix the random seed #' +#' @inheritParams draw_epsilon +#' #' @return a matrix of the instantaneous reproduction number R for the reference #' pathogen/strain/variant for each time step (row) and each location (column) #' drawn from the marginal posterior distribution @@ -358,22 +347,15 @@ draw_epsilon <- function(R, incid, lambda, priors, draw_R <- function(epsilon, incid, lambda, priors, shape_R_flat = NULL, t_min = NULL, t_max = nrow(incid), - seed = NULL) { - if (!is.integer(t_min) || !is.integer(t_max)){ - stop("t_min and t_max must be integers") - } - if (t_min < 2 || t_max < 2){ - stop("t_min and t_max must be >=2") - } - if(t_min > nrow(incid) || t_max > nrow(incid)){ - stop("t_min and t_max must be <= nrow(incid)") + seed = NULL, validate_inputs = FALSE) { + if (validate_inputs) { + R_args <- as.list(environment()) + do.call(check_estimate_advantage_inputs, args = R_args) } + if (any(epsilon < 0)){ stop("epsilon must be > 0") } - if (!is.null(seed) && !is.numeric(seed)){ - stop("seed must be numeric") - } if (!is.null(seed)) set.seed(seed) t <- seq(t_min, t_max, 1) if (is.null(shape_R_flat)) { @@ -412,10 +394,7 @@ draw_R <- function(epsilon, incid, lambda, priors, #' @export compute_si_cutoff <- function(si_distr, miss_at_most = 0.05) { - if (any(colSums(si_distr) != 1)) { - warning("Input SI distributions should sum to 1. Normalising now") - si_distr <- si_distr / colSums(si_distr) - } + cutoff <- 1 - miss_at_most cdf <- apply(si_distr, 2, cumsum) idx <- apply( @@ -529,6 +508,10 @@ compute_t_min <- function(incid, si_distr, miss_at_most) { #' is temporarily assigned to `[, , 1]` of the incidence array. We recommend the #' default value of `TRUE` as we find this to stabilise inference. #' +#' @param validate_inputs a boolean (defaulting to `TRUE`) indicating whether to +#' validate the inputs before running the estimation. We recommned that this +#' is set to `TRUE` as it will ensure any mistmatches in inputs are caught early. +#' #' @return A list with the following elements: #' - `epsilon`: a matrix containing the MCMC chain (thinned and after burnin) #' for the relative transmissibility of the "new" pathogen/strain/variant(s) @@ -580,61 +563,32 @@ compute_t_min <- function(incid, si_distr, miss_at_most) { #' abline(h = 1, col = "red") #' plot(x$R[30, 3, ], type = "l", #' xlab = "Iteration", ylab = "R time 30 location 3") - estimate_advantage <- function(incid, si_distr, priors = default_priors(), mcmc_control = default_mcmc_controls(), t_min = NULL, t_max = nrow(incid), seed = NULL, incid_imported = NULL, precompute = TRUE, - reorder_incid = TRUE) { + reorder_incid = TRUE, + validate_inputs = TRUE) { + if (any(colSums(si_distr) != 1)) { + warning( + "Input SI distributions should sum to 1. Normalising now", call. = FALSE + ) + si_distr <- normalise_si_distr(si_distr) + } + if (is.null(t_min)) { t_min <- compute_t_min(incid, si_distr) } - if (!is.integer(t_min) || !is.integer(t_max)) { - stop("t_min and t_max must be integers") - } - if (t_min < 2 || t_max < 2){ - stop("t_min and t_max must be >=2") - } - if(t_min > nrow(incid) || t_max > nrow(incid)){ - stop("t_min and t_max must be <= nrow(incid)") - } - if (any(si_distr[1,] != 0)){ - stop("Values in the first row of si_distr must be 0") - } - if (any(abs(colSums(si_distr) - 1) > 0.01)) { # allow tolerance - stop("The sum of each column in si_distr should be equal to 1") - } - if (any(si_distr < 0)){ - stop("si_distr must be >=0") - } - if (mcmc_control$n_iter < 0 || !is.integer(mcmc_control$n_iter)){ - stop("n_iter in mcmc_control must be a positive integer") - } - if (mcmc_control$burnin < 0 || !is.integer(mcmc_control$burnin)){ - stop("burnin in mcmc_control must be a positive integer") - } - if (mcmc_control$thin < 0 || !is.integer(mcmc_control$thin)){ - stop("thin in mcmc_control must be a positive integer") - } - if (mcmc_control$n_iter < mcmc_control$burnin + mcmc_control$thin){ - stop("In mcmc_control, n_iter must be greater than burnin + thin") - } - if (!is.null(seed) && !is.numeric(seed)){ - stop("seed must be numeric") - } - if (!is.null(seed)) set.seed(seed) - if (t_min > t_max) { - stop("t_min is greater than t_max. You can specify a smaller t_min or increase t_max.") + estimate_advantage_args <- as.list(environment()) + estimate_advantage_args$si_distr <- si_distr # normalized value + if (validate_inputs) { + do.call(check_estimate_advantage_inputs, args = estimate_advantage_args) } - if (!identical(priors, default_priors())) { - warning("Priors where the mean of epsilon is different from 1 are not currently supported.") - } - T <- nrow(incid) n_loc <- ncol(incid) @@ -820,6 +774,8 @@ process_I_multivariant <- function(incid, incid_imported = NULL) { res } + + ## TODO: check dimensions of objects is correct everywhere ## TODO: fix number of variants to be 2 diff --git a/R/utilities.R b/R/utilities.R index 72154d4d..84c2ecec 100644 --- a/R/utilities.R +++ b/R/utilities.R @@ -518,3 +518,9 @@ si_from_data_valid_distrs <- function(dist) { list(is_dist_valid = dist %in% valid_names, all_valid_distrs = valid_names) } +#' Normalise a serial interval distribution +#' @keywords internal +#' @noRd +normalise_si_distr <- function(si_distr) { + sweep(si_distr, 2, colSums(si_distr), "/") +} diff --git a/R/validate_mv_epiestim_inputs.R b/R/validate_mv_epiestim_inputs.R new file mode 100644 index 00000000..446a8695 --- /dev/null +++ b/R/validate_mv_epiestim_inputs.R @@ -0,0 +1,154 @@ +#' Check incidence input for MV-EpiEstim +#' +#' Check that incid is a 3-dimensional array with non-negatuve enteries. +#' +#' @inheritParams estimate_advantage +#' @returns Silently returns TRUE if the checks are passed, otherwise throws an +#' error. +#' @author Sangeeta Bhatia +#' @keywords internal +check_incidence <- function(incid) { + if (any(incid < 0)) { + stop("incid must be >=0") + } + if (!is.array(incid) || length(dim(incid)) != 3) { + stop( + "incid must be a 3-dimensional array with dimensions time, location, variant", + .call = FALSE + ) + } + invisible(TRUE) +} + + +#' Check priors for MV-EpiEstim +#' +#' Check that priors is a list of the correct format and that the mean of the +#' prior for epsilon is 1 (as currently only priors with mean of epsilon equal +#' to 1 is supported). +#' +#' @inheritParams estimate_advantage +#' @returns Silently returns TRUE if the checks are passed, otherwise throws an +#' error. +#' @author Sangeeta Bhatia +#' @keywords internal +check_priors <- function(priors) { + + if (!identical(priors, default_priors())) { + warning( + "Priors where the mean of epsilon is different from 1 are not currently supported.", + .call = FALSE + ) + } + invisible(TRUE) +} + +#' Check MCMC control parameters for MV-EpiEstim +#' +#' Check that (1) mcmc_control is a list of the correct format, (2) n_iter, burnin +#' and thin are positive integers and (3) n_iter is greater than burnin + thin. +#' +#' @inheritParams estimate_advantage +#' @returns Silently returns TRUE if the checks are passed, otherwise throws an +#' error. +#' @author Sangeeta Bhatia +#' @keywords internal +check_mcmc_control <- function(mcmc_control) { + + if (mcmc_control$n_iter < 0 || !is.integer(mcmc_control$n_iter)) { + stop("n_iter in mcmc_control must be a positive integer", .call = FALSE) + } + if (mcmc_control$burnin < 0 || !is.integer(mcmc_control$burnin)) { + stop("burnin in mcmc_control must be a positive integer", .call = FALSE) + } + if (mcmc_control$thin < 0 || !is.integer(mcmc_control$thin)) { + stop("thin in mcmc_control must be a positive integer", .call = FALSE) + } + if (mcmc_control$n_iter < mcmc_control$burnin + mcmc_control$thin) { + stop("In mcmc_control, n_iter must be greater than burnin + thin", .call = FALSE) + } + invisible(TRUE) +} + +#' Validate t_min and t_max inputs for MV-EpiEstim +#' +#' Check that (1) t_min and t_max are integers, (2) that they are >= 2, that +#' they are <= nrow(incid) and (3) t_min is not greater than t_max. +#' +#' @inheritParams estimate_advantage +#' @returns Silently returns TRUE if the checks are passed, otherwise throws an +#' error. +#' @author Sangeeta Bhatia +#' @keywords internal +check_t_min_t_max <- function(t_min, t_max, incid) { + if (!is.integer(t_min) || !is.integer(t_max)) { + stop("t_min and t_max must be integers", .call = FALSE) + } + if (t_min < 2 || t_max < 2) { + stop("t_min and t_max must be >=2", .call = FALSE) + } + if (t_min > nrow(incid) || t_max > nrow(incid)) { + stop("t_min and t_max must be <= nrow(incid)", .call = FALSE) + } + if (t_min > t_max) { + stop( + "t_min is greater than t_max. You can specify a smaller t_min or increase t_max.", + .call = FALSE + ) + } + invisible(TRUE) +} + + +#' Validate seed input for MV-EpiEstim +#' +#' Validate that the seed is not null and is numeric +#' +#' @inheritParams estimate_advantage +#' @return Silently returns TRUE if the checks are passed, otherwise throws an error. +#' @author Sangeeta Bhatia +#' @keywords internal +check_seed <- function(seed) { + if (!is.null(seed) && !is.numeric(seed)) { + stop("supplied seed is not a valid integer", .call = FALSE) + } + if (!is.null(seed)) set.seed(seed) + invisible(TRUE) +} + + +#' Validate inputs to estimate_advantage +#' +#' Runs all available input checks for \code{\link{estimate_advantage}} based on +#' the supplied arguments. Each check is called only if the corresponding +#' argument is present. +#' +#' @param ... Named arguments passed to \code{\link{estimate_advantage}}. +#' @returns Invisibly returns \code{NULL}; throws an error if any input check +#' fails. +#' @author Sangeeta Bhatia +#' @keywords internal +check_estimate_advantage_inputs <- function(...) { + estimate_advantage_args <- list(...) + arg_names <- names(estimate_advantage_args) + if ("incid" %in% arg_names) + check_incidence(estimate_advantage_args$incid) + + if ("si" %in% arg_names) + apply(estimate_advantage_args$si_distr, 2, check_si_distr) + + if ("priors" %in% arg_names) + check_priors(estimate_advantage_args$priors) + + if ("mcmc_control" %in% arg_names) + check_mcmc_control(estimate_advantage_args$mcmc_control) + + if (all(c("t_min", "t_max", "incid") %in% arg_names)) { + check_t_min_t_max( + estimate_advantage_args$t_min, estimate_advantage_args$t_max, + estimate_advantage_args$incid + ) + } + if ("seed" %in% arg_names) + check_seed(estimate_advantage_args$seed) +} diff --git a/man/EpiEstim-package.Rd b/man/EpiEstim-package.Rd index a5efc74c..b7b766bf 100644 --- a/man/EpiEstim-package.Rd +++ b/man/EpiEstim-package.Rd @@ -26,6 +26,11 @@ Useful links: \author{ \strong{Maintainer}: Anne Cori \email{a.cori@imperial.ac.uk} (\href{https://orcid.org/0000-0002-8443-9162}{ORCID}) +Authors: +\itemize{ + \item Anne Cori \email{a.cori@imperial.ac.uk} (\href{https://orcid.org/0000-0002-8443-9162}{ORCID}) +} + Other contributors: \itemize{ \item Simon Cauchemez [contributor] diff --git a/man/check_estimate_advantage_inputs.Rd b/man/check_estimate_advantage_inputs.Rd new file mode 100644 index 00000000..aaa3f3c1 --- /dev/null +++ b/man/check_estimate_advantage_inputs.Rd @@ -0,0 +1,24 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/validate_mv_epiestim_inputs.R +\name{check_estimate_advantage_inputs} +\alias{check_estimate_advantage_inputs} +\title{Validate inputs to estimate_advantage} +\usage{ +check_estimate_advantage_inputs(...) +} +\arguments{ +\item{...}{Named arguments passed to \code{\link{estimate_advantage}}.} +} +\value{ +Invisibly returns \code{NULL}; throws an error if any input check +fails. +} +\description{ +Runs all available input checks for \code{\link{estimate_advantage}} based on +the supplied arguments. Each check is called only if the corresponding +argument is present. +} +\author{ +Sangeeta Bhatia +} +\keyword{internal} diff --git a/man/check_incidence.Rd b/man/check_incidence.Rd new file mode 100644 index 00000000..a79349e9 --- /dev/null +++ b/man/check_incidence.Rd @@ -0,0 +1,24 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/validate_mv_epiestim_inputs.R +\name{check_incidence} +\alias{check_incidence} +\title{Check incidence input for MV-EpiEstim} +\usage{ +check_incidence(incid) +} +\arguments{ +\item{incid}{a multidimensional array containing values of the incidence +for each time step (1st dimension), location (2nd dimension) and +pathogen/strain/variant (3rd dimension)} +} +\value{ +Silently returns TRUE if the checks are passed, otherwise throws an +error. +} +\description{ +Check that incid is a 3-dimensional array with non-negatuve enteries. +} +\author{ +Sangeeta Bhatia +} +\keyword{internal} diff --git a/man/check_mcmc_control.Rd b/man/check_mcmc_control.Rd new file mode 100644 index 00000000..0e0972ca --- /dev/null +++ b/man/check_mcmc_control.Rd @@ -0,0 +1,24 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/validate_mv_epiestim_inputs.R +\name{check_mcmc_control} +\alias{check_mcmc_control} +\title{Check MCMC control parameters for MV-EpiEstim} +\usage{ +check_mcmc_control(mcmc_control) +} +\arguments{ +\item{mcmc_control}{a list of default MCMC control parameters, obtained by +default from \code{\link[=default_mcmc_controls]{default_mcmc_controls()}}} +} +\value{ +Silently returns TRUE if the checks are passed, otherwise throws an +error. +} +\description{ +Check that (1) mcmc_control is a list of the correct format, (2) n_iter, burnin +and thin are positive integers and (3) n_iter is greater than burnin + thin. +} +\author{ +Sangeeta Bhatia +} +\keyword{internal} diff --git a/man/check_priors.Rd b/man/check_priors.Rd new file mode 100644 index 00000000..d3109e8d --- /dev/null +++ b/man/check_priors.Rd @@ -0,0 +1,27 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/validate_mv_epiestim_inputs.R +\name{check_priors} +\alias{check_priors} +\title{Check priors for MV-EpiEstim} +\usage{ +check_priors(priors) +} +\arguments{ +\item{priors}{a list of prior parameters (shape and scale of a gamma +distribution) for epsilon and R; can be obtained from the function +\code{\link[=default_priors]{default_priors()}}. The prior for R is assumed to be the same for all +time steps and all locations} +} +\value{ +Silently returns TRUE if the checks are passed, otherwise throws an +error. +} +\description{ +Check that priors is a list of the correct format and that the mean of the +prior for epsilon is 1 (as currently only priors with mean of epsilon equal +to 1 is supported). +} +\author{ +Sangeeta Bhatia +} +\keyword{internal} diff --git a/man/check_seed.Rd b/man/check_seed.Rd new file mode 100644 index 00000000..db56ac74 --- /dev/null +++ b/man/check_seed.Rd @@ -0,0 +1,21 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/validate_mv_epiestim_inputs.R +\name{check_seed} +\alias{check_seed} +\title{Validate seed input for MV-EpiEstim} +\usage{ +check_seed(seed) +} +\arguments{ +\item{seed}{a numeric value used to fix the random seed} +} +\value{ +Silently returns TRUE if the checks are passed, otherwise throws an error. +} +\description{ +Validate that the seed is not null and is numeric +} +\author{ +Sangeeta Bhatia +} +\keyword{internal} diff --git a/man/check_t_min_t_max.Rd b/man/check_t_min_t_max.Rd new file mode 100644 index 00000000..c836b0de --- /dev/null +++ b/man/check_t_min_t_max.Rd @@ -0,0 +1,34 @@ +% Generated by roxygen2: do not edit by hand +% Please edit documentation in R/validate_mv_epiestim_inputs.R +\name{check_t_min_t_max} +\alias{check_t_min_t_max} +\title{Validate t_min and t_max inputs for MV-EpiEstim} +\usage{ +check_t_min_t_max(t_min, t_max, incid) +} +\arguments{ +\item{t_min}{an integer > 1 giving the minimum time step to consider in the +estimation. +If \code{NULL}, \code{t_min} is calculated using the function \code{\link[=compute_si_cutoff]{compute_si_cutoff()}} +which gets the maximum (across all variants) of the 95th percentile of the +SI distribution.} + +\item{t_max}{an integer > \code{t_min} and <= \code{nrow(incid)} giving the maximum time +step to consider in the estimation. Default value is \code{nrow(incid)}.} + +\item{incid}{a multidimensional array containing values of the incidence +for each time step (1st dimension), location (2nd dimension) and +pathogen/strain/variant (3rd dimension)} +} +\value{ +Silently returns TRUE if the checks are passed, otherwise throws an +error. +} +\description{ +Check that (1) t_min and t_max are integers, (2) that they are >= 2, that +they are <= nrow(incid) and (3) t_min is not greater than t_max. +} +\author{ +Sangeeta Bhatia +} +\keyword{internal} diff --git a/man/default_priors.Rd b/man/default_priors.Rd index 9ab5ae84..762c2b3e 100644 --- a/man/default_priors.Rd +++ b/man/default_priors.Rd @@ -9,8 +9,8 @@ default_priors() \value{ a list of default parameters for the priors. Values can then be manually edited as in the examples below. -Users could use functions \code{\link[epitrix:gamma_tools]{epitrix::gamma_shapescale2mucv()}} and -\code{\link[epitrix:gamma_tools]{epitrix::gamma_mucv2shapescale()}} to set the shape and scale corresponding +Users could use functions \code{\link[epitrix:gamma_shapescale2mucv]{epitrix::gamma_shapescale2mucv()}} and +\code{\link[epitrix:gamma_mucv2shapescale]{epitrix::gamma_mucv2shapescale()}} to set the shape and scale corresponding to the desired prior mean and coefficient of variation. } \description{ diff --git a/man/draw_R.Rd b/man/draw_R.Rd index 3b947147..790bb9c0 100644 --- a/man/draw_R.Rd +++ b/man/draw_R.Rd @@ -12,7 +12,8 @@ draw_R( shape_R_flat = NULL, t_min = NULL, t_max = nrow(incid), - seed = NULL + seed = NULL, + validate_inputs = FALSE ) } \arguments{ @@ -52,6 +53,11 @@ observations at time step 1 and can therefore only start at time step 2).} step to consider in the estimation. Default value is \code{nrow(incid)}.} \item{seed}{a numeric value used to fix the random seed} + +\item{validate_inputs}{a boolean value indicating whether to validate the +inputs; default value is FALSE as inputs are validated in +\code{estimate_advantage()}, where this function is called. Can be set to TRUE +if this function is called directly for debugging or testing.} } \value{ a matrix of the instantaneous reproduction number R for the reference diff --git a/man/draw_epsilon.Rd b/man/draw_epsilon.Rd index 84b0a3bd..1ee96f15 100644 --- a/man/draw_epsilon.Rd +++ b/man/draw_epsilon.Rd @@ -12,7 +12,8 @@ draw_epsilon( shape_epsilon = NULL, t_min = 2L, t_max = nrow(incid), - seed = NULL + seed = NULL, + validate_inputs = FALSE ) } \arguments{ @@ -51,6 +52,11 @@ observations at time step 1 and can therefore only start at time step 2).} step to consider in the estimation. Default value is \code{nrow(incid)}.} \item{seed}{a numeric value used to fix the random seed} + +\item{validate_inputs}{a boolean value indicating whether to validate the +inputs; default value is FALSE as inputs are validated in +\code{estimate_advantage()}, where this function is called. Can be set to TRUE +if this function is called directly for debugging or testing.} } \value{ A value or vector of values for epsilon for each non reference diff --git a/man/estimate_R.Rd b/man/estimate_R.Rd index 17c1342a..3d6b9db7 100644 --- a/man/estimate_R.Rd +++ b/man/estimate_R.Rd @@ -15,7 +15,7 @@ estimate_R( dt_out = 7L, recon_opt = "naive", iter = 10L, - tol = 1e-06, + tol = 0.000001, grid = list(precision = 0.001, min = -1, max = 1), backimputation_window = 0, date_convention = NULL diff --git a/man/estimate_R_agg.Rd b/man/estimate_R_agg.Rd index f39e7b0e..8711e794 100644 --- a/man/estimate_R_agg.Rd +++ b/man/estimate_R_agg.Rd @@ -9,7 +9,7 @@ estimate_R_agg( dt = 7L, dt_out = 7L, iter = 10L, - tol = 1e-06, + tol = 0.000001, recon_opt = "naive", config = make_config(), method = c("non_parametric_si", "parametric_si"), diff --git a/man/estimate_advantage.Rd b/man/estimate_advantage.Rd index a7566640..a08ac722 100644 --- a/man/estimate_advantage.Rd +++ b/man/estimate_advantage.Rd @@ -14,7 +14,8 @@ estimate_advantage( seed = NULL, incid_imported = NULL, precompute = TRUE, - reorder_incid = TRUE + reorder_incid = TRUE, + validate_inputs = TRUE ) } \arguments{ @@ -63,6 +64,10 @@ incidence array can be internally reordered during the estimation of the transmission advantage. If \code{TRUE}, the most transmissible pathogen/strain/variant is temporarily assigned to \verb{[, , 1]} of the incidence array. We recommend the default value of \code{TRUE} as we find this to stabilise inference.} + +\item{validate_inputs}{a boolean (defaulting to \code{TRUE}) indicating whether to +validate the inputs before running the estimation. We recommned that this +is set to \code{TRUE} as it will ensure any mistmatches in inputs are caught early.} } \value{ A list with the following elements: diff --git a/tests/testthat/test-multivariant-input.R b/tests/testthat/test-multivariant-input.R index ef1d5823..44490a69 100644 --- a/tests/testthat/test-multivariant-input.R +++ b/tests/testthat/test-multivariant-input.R @@ -19,11 +19,11 @@ sidistr_3=cbind(c(0,-0.1,0.7,0.4),c(0,0.2,0.5,0.3)) # negative value test_that("si_distr is specified correctly", { expect_error(compute_lambda(incid=incid_processed, si_distr=sidistr_1), - "Values in the first row of si_distr must be 0") - expect_error(compute_lambda(incid=incid_processed, si_distr=sidistr_2), - "The sum of each column in si_distr should be equal to 1") + "si_distr should be so that si_distr[1] = 0.", fixed = TRUE) + expect_warning(compute_lambda(incid=incid_processed, si_distr=sidistr_2), + "si_distr does not sum to 1.") expect_error(compute_lambda(incid=incid_processed, si_distr=sidistr_3), - "si_distr must be >=0") + "si_distr must be a positive vector.") }) ############################## @@ -45,15 +45,15 @@ tmin3 <- as.integer(nrow(incid)+1) # greater than nrow(incid) test_that("tmin and tmax are specified correctly", { expect_error(draw_epsilon(R=R, incid=incid, lambda=lambda, priors=priors, t_min = tmin1, t_max = nrow(incid), - seed = NULL), + seed = NULL, validate_inputs = TRUE), "t_min and t_max must be integers") expect_error(draw_epsilon(R=R, incid=incid, lambda=lambda, priors=priors, t_min = tmin2, t_max = nrow(incid), - seed = NULL), + seed = NULL, validate_inputs = TRUE), "t_min and t_max must be >=2") expect_error(draw_epsilon(R=R, incid=incid, lambda=lambda, priors=priors, t_min = tmin3, t_max = nrow(incid), - seed = NULL), + seed = NULL, validate_inputs = TRUE), "t_min and t_max must be <= nrow(incid)", fixed=TRUE) }) @@ -64,7 +64,7 @@ test_that("R is specified correctly",{ Rneg[1,1] <- -1 expect_error(draw_epsilon(R=Rneg, incid=incid, lambda=lambda, priors=priors, t_min = 2L, t_max = nrow(incid), - seed = NULL), + seed = NULL, validate_inputs = TRUE), "R must be >= 0") }) @@ -75,8 +75,8 @@ seed <- "a" test_that("seed is specified correctly",{ expect_error(draw_epsilon(R=R, incid=incid, lambda=lambda, priors=priors, t_min = 2L, t_max = nrow(incid), - seed = "a"), - "seed must be numeric") + seed = "a", validate_inputs = TRUE), + "supplied seed is not a valid integer") }) ######################## @@ -100,15 +100,15 @@ tmin3 <- as.integer(nrow(incid)+1) # greater than nrow(incid) test_that("tmin and tmax are specified correctly", { expect_error(draw_R(epsilon=epsilon, incid=incid, lambda=lambda, priors=priors, t_min = tmin1, t_max = nrow(incid), - seed = NULL), + seed = NULL, validate_inputs = TRUE), "t_min and t_max must be integers") expect_error(draw_R(epsilon=epsilon, incid=incid, lambda=lambda, priors=priors, t_min = tmin2, t_max = nrow(incid), - seed = NULL), + seed = NULL, validate_inputs = TRUE), "t_min and t_max must be >=2") expect_error(draw_R(epsilon=epsilon, incid=incid, lambda=lambda, priors=priors, t_min = tmin3, t_max = nrow(incid), - seed = NULL), + seed = NULL, validate_inputs = TRUE), "t_min and t_max must be <= nrow(incid)", fixed=TRUE) }) @@ -120,8 +120,8 @@ seed <- "a" test_that("seed is specified correctly",{ expect_error(draw_R(epsilon=epsilon, incid=incid, lambda=lambda, priors=priors, t_min = 2L, t_max = nrow(incid), - seed = seed), - "seed must be numeric") + seed = seed, validate_inputs = TRUE), + "supplied seed is not a valid integer") }) @@ -132,7 +132,7 @@ epsilon <- -1 test_that("epsilon is specified correctly",{ expect_error(draw_R(epsilon=epsilon, incid=incid, lambda=lambda, priors=priors, t_min = 2L, t_max = nrow(incid), - seed = NULL), + seed = NULL, validate_inputs = TRUE), "epsilon must be > 0") }) @@ -183,17 +183,17 @@ test_that("si_distr is specified correctly", { mcmc_control = default_mcmc_controls(), t_min = 2L, t_max = nrow(incid), seed = NULL), - "Values in the first row of si_distr must be 0") - expect_error(estimate_advantage(incid=incid, si_distr=sidistr_2, priors=priors, + "si_distr should be so that si_distr[1] = 0.", fixed=TRUE) + expect_warning(estimate_advantage(incid=incid, si_distr=sidistr_2, priors=priors, mcmc_control = default_mcmc_controls(), t_min = 2L, t_max = nrow(incid), seed = NULL), - "The sum of each column in si_distr should be equal to 1") + "Input SI distributions should sum to 1. Normalising now") expect_error(estimate_advantage(incid=incid, si_distr=sidistr_3, priors=priors, mcmc_control = default_mcmc_controls(), t_min = 2L, t_max = nrow(incid), seed = NULL), - "si_distr must be >=0") + "si_distr must be a positive vector.") }) @@ -245,7 +245,7 @@ test_that("seed is specified correctly",{ mcmc_control = default_mcmc_controls(), t_min = 2L, t_max = nrow(incid), seed = seed), - "seed must be numeric") + "supplied seed is not a valid integer") }) # convergence check @@ -327,4 +327,5 @@ test_that("the Rt of the reference variant is returned both with and without inc seed = NULL, reorder_incid = FALSE) expect_equal(stats::median(out$R, na.rm = T), stats::median(out2$R, na.rm = T), tolerance = 0.001) -}) \ No newline at end of file +}) +