Skip to content
Open
Show file tree
Hide file tree
Changes from all commits
Commits
File filter

Filter by extension

Filter by extension

Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
2 changes: 1 addition & 1 deletion DESCRIPTION
Original file line number Diff line number Diff line change
Expand Up @@ -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
8 changes: 4 additions & 4 deletions R/EpiEstim-package.R
Original file line number Diff line number Diff line change
Expand Up @@ -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
#'
Expand Down
116 changes: 36 additions & 80 deletions R/gibbs_draws.R
Original file line number Diff line number Diff line change
Expand Up @@ -33,7 +33,7 @@
#' get_shape_R_flat(incid, priors)

get_shape_R_flat <- function(incid, priors, t_min = 2L, t_max = nrow(incid)) {
t <- seq(t_min, t_max, 1)

Check warning on line 36 in R/gibbs_draws.R

View workflow job for this annotation

GitHub Actions / lint-changed-files

file=R/gibbs_draws.R,line=36,col=3,[object_overwrite_linter] 't' is an exported object from package 'base'. Avoid re-using such symbols.
shape <- apply(incid[t, , , drop = FALSE], c(1, 2), sum) + priors$R$shape
as.numeric(shape)
}
Expand Down Expand Up @@ -87,7 +87,7 @@

get_shape_epsilon <- function(incid, lambda, priors,
t_min = 2L, t_max = nrow(incid)) {
t <- seq(t_min, t_max, 1)

Check warning on line 90 in R/gibbs_draws.R

View workflow job for this annotation

GitHub Actions / lint-changed-files

file=R/gibbs_draws.R,line=90,col=3,[object_overwrite_linter] 't' is an exported object from package 'base'. Avoid re-using such symbols.
vnapply(seq(2, dim(lambda)[3]), function(e)
sum(incid[t, , e])) + priors$epsilon$shape
}
Expand Down Expand Up @@ -183,15 +183,6 @@
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])) {
Expand Down Expand Up @@ -242,6 +233,11 @@
#'
#' @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
#'
Expand All @@ -263,34 +259,25 @@
#' 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")

Check warning on line 271 in R/gibbs_draws.R

View workflow job for this annotation

GitHub Actions / lint-changed-files

file=R/gibbs_draws.R,line=271,col=5,[condition_call_linter] Use stop(., call. = FALSE) not to display the call in an error message.
}
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)

Check warning on line 274 in R/gibbs_draws.R

View workflow job for this annotation

GitHub Actions / lint-changed-files

file=R/gibbs_draws.R,line=274,col=3,[object_overwrite_linter] 't' is an exported object from package 'base'. Avoid re-using such symbols.
if (is.null(shape_epsilon)) {
shape_epsilon <- get_shape_epsilon(incid, lambda, priors, t_min, t_max)
}
rate <- vnapply(seq(2, dim(lambda)[3]), function(e)
sum(R[t, ] * lambda[t, , e]) + 1 / priors$epsilon$scale)
scale <- 1 / rate

Check warning on line 280 in R/gibbs_draws.R

View workflow job for this annotation

GitHub Actions / lint-changed-files

file=R/gibbs_draws.R,line=280,col=3,[object_overwrite_linter] 'scale' is an exported object from package 'base'. Avoid re-using such symbols.
stats::rgamma(dim(lambda)[3] - 1, shape = shape_epsilon, scale = scale)
}

Expand Down Expand Up @@ -333,6 +320,8 @@
#'
#' @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
Expand All @@ -358,24 +347,17 @@
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")

Check warning on line 357 in R/gibbs_draws.R

View workflow job for this annotation

GitHub Actions / lint-changed-files

file=R/gibbs_draws.R,line=357,col=5,[condition_call_linter] Use stop(., call. = FALSE) not to display the call in an error message.
}
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)

Check warning on line 360 in R/gibbs_draws.R

View workflow job for this annotation

GitHub Actions / lint-changed-files

file=R/gibbs_draws.R,line=360,col=3,[object_overwrite_linter] 't' is an exported object from package 'base'. Avoid re-using such symbols.
if (is.null(shape_R_flat)) {
shape_R_flat <- get_shape_R_flat(incid, priors, t_min, t_max)
}
Expand All @@ -387,11 +369,11 @@
temp <- temp + epsilon[var - 1] * lambda[t, , var]
}
rate <- temp + 1 / priors$R$scale
scale <- 1 / rate

Check warning on line 372 in R/gibbs_draws.R

View workflow job for this annotation

GitHub Actions / lint-changed-files

file=R/gibbs_draws.R,line=372,col=3,[object_overwrite_linter] 'scale' is an exported object from package 'base'. Avoid re-using such symbols.
scale_flat <- as.numeric(scale)
R_flat <- stats::rgamma(length(shape_R_flat), shape = shape_R_flat, scale = scale_flat)
R_fill <- matrix(R_flat, nrow = length(t), ncol = ncol(incid))
R <- matrix(NA, nrow(incid), ncol(incid))

Check warning on line 376 in R/gibbs_draws.R

View workflow job for this annotation

GitHub Actions / lint-changed-files

file=R/gibbs_draws.R,line=376,col=3,[object_overwrite_linter] 'R' is an exported object from package 'tools'. Avoid re-using such symbols.
R[t, ] <- R_fill
R
}
Expand All @@ -412,10 +394,7 @@
#' @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(
Expand Down Expand Up @@ -445,7 +424,7 @@
function(vec) Position(function(x) x > 0, vec)
)
if (anyNA(t_min_incid)) {
warning(

Check warning on line 427 in R/gibbs_draws.R

View workflow job for this annotation

GitHub Actions / lint-changed-files

file=R/gibbs_draws.R,line=427,col=5,[condition_call_linter] Use warning(., call. = FALSE) not to display the call in an error message.
"For some variants/locations, incidence is
always zero. This will cause estimate_advantage to fail."
)
Expand Down Expand Up @@ -529,6 +508,10 @@
#' 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)
Expand Down Expand Up @@ -580,61 +563,32 @@
#' 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)
}
Comment thread
sangeetabhatia03 marked this conversation as resolved.

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)
Comment thread
sangeetabhatia03 marked this conversation as resolved.
}

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)

Expand Down Expand Up @@ -820,6 +774,8 @@
res
}



## TODO: check dimensions of objects is correct everywhere
## TODO: fix number of variants to be 2

6 changes: 6 additions & 0 deletions R/utilities.R
Original file line number Diff line number Diff line change
Expand Up @@ -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), "/")
}
154 changes: 154 additions & 0 deletions R/validate_mv_epiestim_inputs.R
Original file line number Diff line number Diff line change
@@ -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)
}
Loading
Loading