From eac7a3d53afb7bf3e420c2bb48646f972bd6b49f Mon Sep 17 00:00:00 2001 From: jamesaazam Date: Wed, 9 Oct 2024 12:43:07 +0100 Subject: [PATCH] Run styler + lintr --- R/borel.R | 34 +++- R/checks.R | 8 +- R/epichains.R | 88 ++++----- R/helpers.R | 11 +- R/likelihood.R | 70 ++++--- R/simulate.R | 16 +- R/stat_likelihoods.R | 61 ++++-- R/utils.R | 4 +- README.Rmd | 4 +- man/likelihood.Rd | 18 +- man/print.epichains.Rd | 12 +- man/print.epichains_summary.Rd | 10 +- man/summary.epichains_summary.Rd | 10 +- tests/spelling.R | 10 +- tests/testthat/helper-state.R | 34 ++-- tests/testthat/setup-options.R | 6 +- tests/testthat/test-checks.R | 300 ++++++++++++++--------------- tests/testthat/test-epichains.R | 2 +- tests/testthat/test-likelihood.R | 247 ++++++++++++------------ tests/testthat/test-simulate.R | 2 +- vignettes/interventions.Rmd | 26 ++- vignettes/projecting_incidence.Rmd | 4 +- 22 files changed, 518 insertions(+), 459 deletions(-) diff --git a/R/borel.R b/R/borel.R index f2d28959..79991922 100644 --- a/R/borel.R +++ b/R/borel.R @@ -11,10 +11,15 @@ #' dborel(1:5, 1) dborel <- function(x, mu, log = FALSE) { checkmate::assert_numeric( - x, lower = 1, upper = Inf + x, + lower = 1, + upper = Inf ) checkmate::assert_number( - mu, lower = 0, finite = TRUE, na.ok = FALSE + mu, + lower = 0, + finite = TRUE, + na.ok = FALSE ) ld <- -mu * x + (x - 1) * log(mu * x) - lgamma(x + 1) @@ -41,10 +46,16 @@ dborel <- function(x, mu, log = FALSE) { #' rborel(5, 1) rborel <- function(n, mu, censor_at = Inf) { checkmate::assert_number( - n, lower = 1, finite = TRUE, na.ok = FALSE + n, + lower = 1, + finite = TRUE, + na.ok = FALSE ) checkmate::assert_number( - mu, lower = 0, finite = TRUE, na.ok = FALSE + mu, + lower = 0, + finite = TRUE, + na.ok = FALSE ) # Run simulations out <- simulate_chain_stats( @@ -81,16 +92,22 @@ rgborel <- function(n, size, prob, mu, censor_at = Inf) { # its "correct implementation" for documentation/clarity purposes, as well as # for simulations. checkmate::assert_number( - size, finite = TRUE, lower = 0 + size, + finite = TRUE, + lower = 0 ) if (!missing(prob)) { checkmate::assert_number( - prob, lower = 0, upper = 1 + prob, + lower = 0, + upper = 1 ) } if (!missing(mu)) { checkmate::assert_number( - mu, finite = TRUE, lower = 0 + mu, + finite = TRUE, + lower = 0 ) } if (!missing(prob)) { @@ -101,6 +118,7 @@ rgborel <- function(n, size, prob, mu, censor_at = Inf) { x <- rgamma(n, shape = size, rate = size / mu) # then, sample from borel return(vapply( - x, rborel, n = 1, censor_at = censor_at, FUN.VALUE = numeric(1) + x, rborel, + n = 1, censor_at = censor_at, FUN.VALUE = numeric(1) )) } diff --git a/R/checks.R b/R/checks.R index 2fb2a0dc..e67e292b 100644 --- a/R/checks.R +++ b/R/checks.R @@ -72,9 +72,9 @@ # check that stat_threshold is an integer or Inf. checkmate::assert( checkmate::anyInfinite(stat_threshold), - checkmate::check_integerish( - stat_threshold, - lower = 1 + checkmate::check_integerish( + stat_threshold, + lower = 1 ), combine = "or" ) @@ -102,7 +102,7 @@ if (!is.null(generation_time)) { .check_generation_time_valid(generation_time) } else if (tf_specified) { - stop("If `tf` is specified, `generation_time` must be specified too.") + stop("If `tf` is specified, `generation_time` must be specified too.") } checkmate::assert_number( tf, diff --git a/R/epichains.R b/R/epichains.R index 1eed2b6a..ee2b4dc3 100644 --- a/R/epichains.R +++ b/R/epichains.R @@ -16,11 +16,11 @@ #' @author James M. Azam #' @keywords internal .new_epichains <- function(sim_df, - n_chains, - statistic, - offspring_dist, - stat_threshold, - track_pop) { + n_chains, + statistic, + offspring_dist, + stat_threshold, + track_pop) { # Assemble the elements of the object obj <- sim_df class(obj) <- c("epichains", class(obj)) @@ -55,11 +55,11 @@ #' @author James M. Azam #' @keywords internal .epichains <- function(sim_df, - n_chains, - offspring_dist, - track_pop, - statistic = c("size", "length"), - stat_threshold = Inf) { + n_chains, + offspring_dist, + track_pop, + statistic = c("size", "length"), + stat_threshold = Inf) { # Check that inputs are well specified checkmate::assert_data_frame(sim_df, min.cols = 3, min.rows = n_chains) checkmate::assert_integerish( @@ -112,10 +112,10 @@ #' @author James M. Azam #' @keywords internal .new_epichains_summary <- function(chains_summary, - n_chains, - statistic, - offspring_dist, - stat_threshold) { + n_chains, + statistic, + offspring_dist, + stat_threshold) { # Assemble the elements of the object obj <- chains_summary class(obj) <- c("epichains_summary", class(chains_summary)) @@ -141,10 +141,10 @@ #' @author James M. Azam #' @keywords internal .epichains_summary <- function(chains_summary, - n_chains, - offspring_dist, - statistic = c("size", "length"), - stat_threshold = Inf) { + n_chains, + offspring_dist, + statistic = c("size", "length"), + stat_threshold = Inf) { # chain_summary can sometimes contain infinite values, so check # that finite elements are integerish. checkmate::check_integerish( @@ -194,12 +194,12 @@ #' # population up to chain size 10. #' set.seed(32) #' chains_pois_offspring <- simulate_chains( -#' n_chains = 10, -#' statistic = "size", -#' offspring_dist = rpois, -#' stat_threshold = 10, -#' generation_time = function(n) rep(3, n), -#' lambda = 2 +#' n_chains = 10, +#' statistic = "size", +#' offspring_dist = rpois, +#' stat_threshold = 10, +#' generation_time = function(n) rep(3, n), +#' lambda = 2 #' ) #' chains_pois_offspring # Print the object print.epichains <- function(x, ...) { @@ -226,11 +226,11 @@ print.epichains <- function(x, ...) { #' # population up to chain size 10. #' set.seed(32) #' chain_summary_print_eg <- simulate_chain_stats( -#' n_chains = 10, -#' statistic = "size", -#' offspring_dist = rpois, -#' stat_threshold = 10, -#' lambda = 2 +#' n_chains = 10, +#' statistic = "size", +#' offspring_dist = rpois, +#' stat_threshold = 10, +#' lambda = 2 #' ) #' chain_summary_print_eg # Print the object print.epichains_summary <- function(x, ...) { @@ -322,9 +322,9 @@ format.epichains_summary <- function(x, ...) { "Max: %s", ifelse( is.infinite( - statistics[["max_stat"]]), - paste0(">=", attr(x, "stat_threshold") + statistics[["max_stat"]] ), + paste0(">=", attr(x, "stat_threshold")), statistics[["max_stat"]] ) ), @@ -332,9 +332,9 @@ format.epichains_summary <- function(x, ...) { "Min: %s", ifelse( is.infinite( - statistics[["min_stat"]]), - paste0(">=", attr(x, "stat_threshold") + statistics[["min_stat"]] ), + paste0(">=", attr(x, "stat_threshold")), statistics[["min_stat"]] ) ) @@ -455,11 +455,11 @@ summary.epichains <- function(object, ...) { #' # population up to chain size 10. #' set.seed(32) #' chain_stats <- simulate_chain_stats( -#' n_chains = 10, -#' statistic = "size", -#' offspring_dist = rpois, -#' stat_threshold = 10, -#' lambda = 2 +#' n_chains = 10, +#' statistic = "size", +#' offspring_dist = rpois, +#' stat_threshold = 10, +#' lambda = 2 #' ) #' summary(chain_stats) summary.epichains_summary <- function(object, ...) { @@ -526,7 +526,7 @@ summary.epichains_summary <- function(object, ...) { stopifnot( "object does not contain the correct columns" = c("chain", "infector", "infectee", "generation") %in% - colnames(x), + colnames(x), "column `chain` must be a numeric" = is.numeric(x$chain), "column `infector` must be a numeric" = @@ -643,11 +643,11 @@ tail.epichains <- function(x, ...) { #' cases_per_gen <- aggregate(chains, by = "generation") #' head(cases_per_gen) aggregate.epichains <- function(x, - by = c( - "time", - "generation" - ), - ...) { + by = c( + "time", + "generation" + ), + ...) { .validate_epichains(x) checkmate::assert_string(by) # Get grouping variable diff --git a/R/helpers.R b/R/helpers.R index bc5e49ee..62763864 100644 --- a/R/helpers.R +++ b/R/helpers.R @@ -10,8 +10,7 @@ #' @keywords internal .update_chain_stat <- function(stat_type, stat_latest, n_offspring) { return( - switch( - stat_type, + switch(stat_type, size = stat_latest + n_offspring, length = stat_latest + pmin(1, n_offspring), stop("stat_type must be 'size' or 'length'") @@ -27,8 +26,7 @@ #' @keywords internal .get_statistic_func <- function(chain_statistic) { return( - switch( - chain_statistic, + switch(chain_statistic, size = .rbinom_size, length = .rgen_length, stop("chain_statistic must be 'size' or 'length'") @@ -73,7 +71,6 @@ offspring_func_pars, n_offspring, chains) { - possible_new_offspring <- do.call( offspring_func, c( @@ -104,8 +101,8 @@ #' given the current susceptible population size. #' @keywords internal .get_susceptible_offspring <- function(new_offspring, - susc_pop, - pop) { + susc_pop, + pop) { # We first adjust for the case where susceptible can be Inf but prob can only # be maximum 1. binom_prob <- min(1, susc_pop / pop, na.rm = TRUE) diff --git a/R/likelihood.R b/R/likelihood.R index 07497803..f52b101b 100644 --- a/R/likelihood.R +++ b/R/likelihood.R @@ -61,16 +61,16 @@ #' # Example using an object #' set.seed(32) #' epichains_obj_eg <- simulate_chains( -#' n_chains = 10, -#' pop = 100, -#' percent_immune = 0, -#' statistic = "size", -#' offspring_dist = rnbinom, -#' stat_threshold = 10, -#' generation_time = function(n) rep(3, n), -#' mu = 2, -#' size = 0.2 -#') +#' n_chains = 10, +#' pop = 100, +#' percent_immune = 0, +#' statistic = "size", +#' offspring_dist = rnbinom, +#' stat_threshold = 10, +#' generation_time = function(n) rep(3, n), +#' mu = 2, +#' size = 0.2 +#' ) #' #' epichains_obj_eg_lik <- likelihood( #' chains = epichains_obj_eg, @@ -103,7 +103,7 @@ #' ) #' epichains_summary_eg_lik #' @export -#nolint start: cyclocomp_linter +# nolint start: cyclocomp_linter likelihood <- function(chains, statistic = c("size", "length"), offspring_dist, nsim_obs, obs_prob = 1, stat_threshold = Inf, log = TRUE, exclude = NULL, individual = FALSE, ...) { @@ -112,7 +112,8 @@ likelihood <- function(chains, statistic = c("size", "length"), offspring_dist, ## Input checking checkmate::assert( checkmate::check_numeric( - chains, lower = 0, upper = Inf, any.missing = FALSE + chains, + lower = 0, upper = Inf, any.missing = FALSE ), checkmate::check_class( chains, "epichains" @@ -128,16 +129,20 @@ likelihood <- function(chains, statistic = c("size", "length"), offspring_dist, ) .check_offspring_func_valid(offspring_dist) checkmate::assert_number( - obs_prob, lower = 0, upper = 1, finite = TRUE, na.ok = FALSE + obs_prob, + lower = 0, upper = 1, finite = TRUE, na.ok = FALSE ) checkmate::assert_logical( - log, any.missing = FALSE, all.missing = FALSE, len = 1 + log, + any.missing = FALSE, all.missing = FALSE, len = 1 ) checkmate::assert_logical( - individual, any.missing = FALSE, all.missing = FALSE, len = 1 + individual, + any.missing = FALSE, all.missing = FALSE, len = 1 ) checkmate::assert_numeric( - exclude, null.ok = TRUE + exclude, + null.ok = TRUE ) # likelihood is designed to work with numeric objects to objects # need to be coerced to objects (numeric vector under @@ -167,16 +172,16 @@ likelihood <- function(chains, statistic = c("size", "length"), offspring_dist, # censor the chains to be at most stat_threshold chains <- pmin(chains, stat_threshold) } else if (any(is.infinite(chains)) && !.is_epichains_summary(chains)) { - # Developer note: objects also pass the`is.numeric()` - # check, so we need to check if the object is not an - # object. - # Numeric chains can only contain Inf values based on human error/decision - # to censor it with infinite values, so we ask the user to fix that. - stop( - "`chains` must be censored with a finite value. ", - "Replace the `Inf` values with a finite value.", - call. = FALSE - ) + # Developer note: objects also pass the`is.numeric()` + # check, so we need to check if the object is not an + # object. + # Numeric chains can only contain Inf values based on human error/decision + # to censor it with infinite values, so we ask the user to fix that. + stop( + "`chains` must be censored with a finite value. ", + "Replace the `Inf` values with a finite value.", + call. = FALSE + ) } # Apply the observation process @@ -185,7 +190,8 @@ likelihood <- function(chains, statistic = c("size", "length"), offspring_dist, stop("'nsim_obs' must be specified if 'obs_prob' is < 1") } else { checkmate::assert_integerish( - nsim_obs, lower = 1 + nsim_obs, + lower = 1 ) } @@ -228,8 +234,9 @@ likelihood <- function(chains, statistic = c("size", "length"), offspring_dist, ## calculate log-likelihoods possible_internal_func <- paste0(".", ll_func) if (exists(possible_internal_func, - where = asNamespace("epichains"), - mode = "function") + where = asNamespace("epichains"), + mode = "function" + ) ) { func <- get(possible_internal_func) likelihoods[calc_sizes] <- do.call(func, c(list(x = calc_sizes), pars)) @@ -259,8 +266,7 @@ likelihood <- function(chains, statistic = c("size", "length"), offspring_dist, likelihoods[exclude] <- -Inf stat_rep_list <- lapply(stat_rep_list, function(y) { y[!(y %in% exclude)] - } - ) + }) } ## assign likelihoods @@ -285,4 +291,4 @@ likelihood <- function(chains, statistic = c("size", "length"), offspring_dist, return(chains_likelihood) } -#nolint end +# nolint end diff --git a/R/simulate.R b/R/simulate.R index 22a2c017..7d9c9fca 100644 --- a/R/simulate.R +++ b/R/simulate.R @@ -364,9 +364,9 @@ simulate_chains <- function(n_chains, #' the `likelihood()` function and for sampling from the borel #' distribution (See ?epichains::rborel). It is used extensively in the #' vignette on -#nolint start +# nolint start #' [modelling disease control](https://epiverse-trace.github.io/epichains/articles/interventions.html), -#nolint end +# nolint end #' where only data on observed chain sizes and lengths are available. #' @author James M. Azam, Sebastian Funk #' @examples @@ -395,12 +395,12 @@ simulate_chains <- function(n_chains, #' ) #' @export simulate_chain_stats <- function(n_chains, - statistic = c("size", "length"), - offspring_dist, - ..., - stat_threshold = Inf, - pop = Inf, - percent_immune = 0) { + statistic = c("size", "length"), + offspring_dist, + ..., + stat_threshold = Inf, + pop = Inf, + percent_immune = 0) { # Check offspring and population-related arguments .check_sim_args( n_chains = n_chains, diff --git a/R/stat_likelihoods.R b/R/stat_likelihoods.R index 138ebe0c..aa7f5223 100644 --- a/R/stat_likelihoods.R +++ b/R/stat_likelihoods.R @@ -7,10 +7,14 @@ #' @keywords internal .pois_size_ll <- function(x, lambda) { checkmate::assert_numeric( - x, lower = 0, any.missing = FALSE + x, + lower = 0, + any.missing = FALSE ) checkmate::assert_number( - lambda, finite = TRUE, lower = 0 + lambda, + finite = TRUE, + lower = 0 ) out <- (x - 1) * log(lambda) - lambda * x + (x - 2) * log(x) - lgamma(x) @@ -27,19 +31,27 @@ #' @keywords internal .nbinom_size_ll <- function(x, size, prob, mu) { checkmate::assert_numeric( - x, lower = 0, any.missing = FALSE + x, + lower = 0, + any.missing = FALSE ) checkmate::assert_number( - size, finite = TRUE, lower = 0 + size, + finite = TRUE, + lower = 0 ) if (!missing(prob)) { checkmate::assert_number( - prob, lower = 0, upper = 1 + prob, + lower = 0, + upper = 1 ) } if (!missing(mu)) { checkmate::assert_number( - mu, finite = TRUE, lower = 0 + mu, + finite = TRUE, + lower = 0 ) } if (!missing(prob)) { @@ -61,19 +73,27 @@ #' @keywords internal .gborel_size_ll <- function(x, size, prob, mu) { checkmate::assert_numeric( - x, lower = 0, any.missing = FALSE + x, + lower = 0, + any.missing = FALSE ) checkmate::assert_number( - size, finite = TRUE, lower = 0 + size, + finite = TRUE, + lower = 0 ) if (!missing(prob)) { checkmate::assert_number( - prob, lower = 0, upper = 1 + prob, + lower = 0, + upper = 1 ) } if (!missing(mu)) { checkmate::assert_number( - mu, finite = TRUE, lower = 0 + mu, + finite = TRUE, + lower = 0 ) } @@ -97,10 +117,14 @@ #' @keywords internal .pois_length_ll <- function(x, lambda) { checkmate::assert_numeric( - x, lower = 0, any.missing = FALSE + x, + lower = 0, + any.missing = FALSE ) checkmate::assert_number( - lambda, finite = TRUE, lower = 0 + lambda, + finite = TRUE, + lower = 0 ) ## iterated exponential function @@ -124,10 +148,13 @@ #' @keywords internal .geom_length_ll <- function(x, prob) { checkmate::assert_numeric( - x, lower = 0, any.missing = FALSE + x, + lower = 0, any.missing = FALSE ) checkmate::assert_number( - prob, lower = 0, upper = 1 + prob, + lower = 0, + upper = 1 ) lambda <- 1 / prob @@ -155,10 +182,12 @@ #' @author Sebastian Funk James M. Azam #' @keywords internal .offspring_ll <- function(x, offspring_dist, statistic, - nsim_offspring = 100, ...) { + nsim_offspring = 100, ...) { # Input checking checkmate::assert_numeric( - x, lower = 0, any.missing = FALSE + x, + lower = 0, + any.missing = FALSE ) # Remaining checks are done in simulate_chain_stats() # Simulate the chains diff --git a/R/utils.R b/R/utils.R index 7981064c..40f1bf2c 100644 --- a/R/utils.R +++ b/R/utils.R @@ -7,7 +7,9 @@ #' @keywords internal .complementary_logprob <- function(x) { checkmate::assert_numeric( - x, lower = -Inf, upper = 0 + x, + lower = -Inf, + upper = 0 ) out <- tryCatch(log1p(-sum(exp(x))), error = function(e) -Inf) diff --git a/README.Rmd b/README.Rmd index d392b413..6c764550 100644 --- a/README.Rmd +++ b/README.Rmd @@ -124,7 +124,9 @@ sim_chains <- simulate_chains( statistic = "size", offspring_dist = rpois, stat_threshold = 25, - generation_time = function(n) {rep(3, n)}, # constant generation time of 3 + generation_time = function(n) { + rep(3, n) + }, # constant generation time of 3 lambda = 1 # mean of the Poisson distribution ) # View the head of the simulation diff --git a/man/likelihood.Rd b/man/likelihood.Rd index 2e1ccb5b..244cd89a 100644 --- a/man/likelihood.Rd +++ b/man/likelihood.Rd @@ -109,15 +109,15 @@ likelihood( # Example using an object set.seed(32) epichains_obj_eg <- simulate_chains( - n_chains = 10, - pop = 100, - percent_immune = 0, - statistic = "size", - offspring_dist = rnbinom, - stat_threshold = 10, - generation_time = function(n) rep(3, n), - mu = 2, - size = 0.2 + n_chains = 10, + pop = 100, + percent_immune = 0, + statistic = "size", + offspring_dist = rnbinom, + stat_threshold = 10, + generation_time = function(n) rep(3, n), + mu = 2, + size = 0.2 ) epichains_obj_eg_lik <- likelihood( diff --git a/man/print.epichains.Rd b/man/print.epichains.Rd index 58bf795b..c1874a02 100644 --- a/man/print.epichains.Rd +++ b/man/print.epichains.Rd @@ -23,12 +23,12 @@ Print an \verb{} object # population up to chain size 10. set.seed(32) chains_pois_offspring <- simulate_chains( - n_chains = 10, - statistic = "size", - offspring_dist = rpois, - stat_threshold = 10, - generation_time = function(n) rep(3, n), - lambda = 2 + n_chains = 10, + statistic = "size", + offspring_dist = rpois, + stat_threshold = 10, + generation_time = function(n) rep(3, n), + lambda = 2 ) chains_pois_offspring # Print the object } diff --git a/man/print.epichains_summary.Rd b/man/print.epichains_summary.Rd index 7a66e409..1071b48e 100644 --- a/man/print.epichains_summary.Rd +++ b/man/print.epichains_summary.Rd @@ -28,11 +28,11 @@ limit. See \code{?epichains_summary()} for the definition of \code{stat_threshol # population up to chain size 10. set.seed(32) chain_summary_print_eg <- simulate_chain_stats( - n_chains = 10, - statistic = "size", - offspring_dist = rpois, - stat_threshold = 10, - lambda = 2 + n_chains = 10, + statistic = "size", + offspring_dist = rpois, + stat_threshold = 10, + lambda = 2 ) chain_summary_print_eg # Print the object } diff --git a/man/summary.epichains_summary.Rd b/man/summary.epichains_summary.Rd index 417dc6bb..127d1423 100644 --- a/man/summary.epichains_summary.Rd +++ b/man/summary.epichains_summary.Rd @@ -30,11 +30,11 @@ Summary method for \verb{} class # population up to chain size 10. set.seed(32) chain_stats <- simulate_chain_stats( - n_chains = 10, - statistic = "size", - offspring_dist = rpois, - stat_threshold = 10, - lambda = 2 + n_chains = 10, + statistic = "size", + offspring_dist = rpois, + stat_threshold = 10, + lambda = 2 ) summary(chain_stats) } diff --git a/tests/spelling.R b/tests/spelling.R index fc341c34..647406c2 100644 --- a/tests/spelling.R +++ b/tests/spelling.R @@ -1,7 +1,7 @@ if (requireNamespace("spelling", quietly = TRUE)) { - spelling::spell_check_test( - vignettes = TRUE, - error = TRUE, - skip_on_cran = TRUE - ) + spelling::spell_check_test( + vignettes = TRUE, + error = TRUE, + skip_on_cran = TRUE + ) } diff --git a/tests/testthat/helper-state.R b/tests/testthat/helper-state.R index 29466710..4910faa6 100644 --- a/tests/testthat/helper-state.R +++ b/tests/testthat/helper-state.R @@ -6,21 +6,21 @@ # We add a test on R >= 4.0.0 because some functions such as # `globalCallingHandlers()` did not exist before. if (getRversion() >= "4.0.0") { - testthat::set_state_inspector(function() { - list( - attached = search(), - connections = getAllConnections(), - cwd = getwd(), - envvars = Sys.getenv(), - handlers = globalCallingHandlers(), - libpaths = .libPaths(), - locale = Sys.getlocale(), - options = options(), - par = par(), - packages = .packages(all.available = TRUE), - sink = sink.number(), - timezone = Sys.timezone(), - NULL - ) - }) + testthat::set_state_inspector(function() { + list( + attached = search(), + connections = getAllConnections(), + cwd = getwd(), + envvars = Sys.getenv(), + handlers = globalCallingHandlers(), + libpaths = .libPaths(), + locale = Sys.getlocale(), + options = options(), + par = par(), + packages = .packages(all.available = TRUE), + sink = sink.number(), + timezone = Sys.timezone(), + NULL + ) + }) } diff --git a/tests/testthat/setup-options.R b/tests/testthat/setup-options.R index 0a86b247..4ca3f1c0 100644 --- a/tests/testthat/setup-options.R +++ b/tests/testthat/setup-options.R @@ -1,7 +1,7 @@ # We want to flag partial matching as part of our testing & continuous # integration process because it makes code more brittle. options( - warnPartialMatchAttr = TRUE, - warnPartialMatchDollar = TRUE, - warnPartialMatchArgs = TRUE + warnPartialMatchAttr = TRUE, + warnPartialMatchDollar = TRUE, + warnPartialMatchArgs = TRUE ) diff --git a/tests/testthat/test-checks.R b/tests/testthat/test-checks.R index 38a9682a..2469a514 100644 --- a/tests/testthat/test-checks.R +++ b/tests/testthat/test-checks.R @@ -1,169 +1,169 @@ # A function that passes all checks in .check_sim_args but returns an error # message if the supplied arguments are invalid .check_sim_args_default <- function(...) { - default_args <- list( - n_chains = 10, - statistic = "size", - offspring_dist = rpois, - stat_threshold = 10, - pop = 10, - percent_immune = 0.1 - ) - # Modify the default arguments with the user's arguments - new_args <- modifyList( - default_args, - list(...) - ) - # Run the check_sim_args function and capture the output - out <- tryCatch( - expr = { - do.call( - .check_sim_args, - new_args - ) - }, - error = function(e) { - stop(e) - } - ) - # Return the output - return(out) + default_args <- list( + n_chains = 10, + statistic = "size", + offspring_dist = rpois, + stat_threshold = 10, + pop = 10, + percent_immune = 0.1 + ) + # Modify the default arguments with the user's arguments + new_args <- modifyList( + default_args, + list(...) + ) + # Run the check_sim_args function and capture the output + out <- tryCatch( + expr = { + do.call( + .check_sim_args, + new_args + ) + }, + error = function(e) { + stop(e) + } + ) + # Return the output + return(out) } # A function that passes all checks in `.check_time_args` but returns an error # message if the supplied arguments are invalid .check_time_args_default <- function(...) { - default_args <- list( - t0 = 0, - tf_specified = FALSE, # tf is not specified but default tf = Inf below - tf = Inf, - generation_time = NULL - ) - new_args <- modifyList( - default_args, - list(...) - ) - out <- tryCatch( - expr = { - do.call( - .check_time_args, - new_args - ) - }, - error = function(e) { - stop(e) - } - ) - return(out) + default_args <- list( + t0 = 0, + tf_specified = FALSE, # tf is not specified but default tf = Inf below + tf = Inf, + generation_time = NULL + ) + new_args <- modifyList( + default_args, + list(...) + ) + out <- tryCatch( + expr = { + do.call( + .check_time_args, + new_args + ) + }, + error = function(e) { + stop(e) + } + ) + return(out) } test_that("Smaller checker functions work", { - expect_error( - .check_offspring_func_valid(rrpois), - "not found" - ) - expect_error( - .check_generation_time_valid("a"), - "Must be a function" - ) - expect_error( - .check_generation_time_valid(function(x) rep("a", 10)), - "numeric" - ) - expect_error( - .check_generation_time_valid(function(x) 3), - "Must have length" - ) - expect_no_error( - .check_statistic_args( - statistic = "size", - stat_threshold = 10 - ) - ) + expect_error( + .check_offspring_func_valid(rrpois), + "not found" + ) + expect_error( + .check_generation_time_valid("a"), + "Must be a function" + ) + expect_error( + .check_generation_time_valid(function(x) rep("a", 10)), + "numeric" + ) + expect_error( + .check_generation_time_valid(function(x) 3), + "Must have length" + ) + expect_no_error( + .check_statistic_args( + statistic = "size", + stat_threshold = 10 + ) + ) }) test_that(".check_sim_args() returns errors", { - # Checks with .check_sim_args - expect_no_error( - .check_sim_args_default() - ) - # n_chains must be >= 1 - expect_error( - .check_sim_args_default( - n_chains = 0 - ), - "Must be >= 1." - ) - # statistic can only be "size" or "length" - expect_error( - .check_sim_args_default( - statistic = "duration" - ), - "Must be element of set \\{'size','length'\\}" - ) - # offspring_dist must be a function - expect_error( - .check_sim_args_default( - offspring_dist = "rpois" - ), - "Must be a function, not 'character'" - ) - # offspring_dist must be a known function (in the environment) - expect_error( - .check_sim_args_default( - offspring_dist = r - ), - "object 'r' not found" - ) - # stat_threshold must be >= 1 - expect_error( - .check_sim_args_default( - stat_threshold = 0 - ), - "Assertion failed." - ) - # pop cannot be negative - expect_error( - .check_sim_args_default( - pop = -1 - ), - "Element 1 is not >= 1." - ) - # percent_immune must be in between 0 and 1 - expect_error( - .check_sim_args_default( - percent_immune = 1.1 - ), - "Element 1 is not <= 1." - ) + # Checks with .check_sim_args + expect_no_error( + .check_sim_args_default() + ) + # n_chains must be >= 1 + expect_error( + .check_sim_args_default( + n_chains = 0 + ), + "Must be >= 1." + ) + # statistic can only be "size" or "length" + expect_error( + .check_sim_args_default( + statistic = "duration" + ), + "Must be element of set \\{'size','length'\\}" + ) + # offspring_dist must be a function + expect_error( + .check_sim_args_default( + offspring_dist = "rpois" + ), + "Must be a function, not 'character'" + ) + # offspring_dist must be a known function (in the environment) + expect_error( + .check_sim_args_default( + offspring_dist = r + ), + "object 'r' not found" + ) + # stat_threshold must be >= 1 + expect_error( + .check_sim_args_default( + stat_threshold = 0 + ), + "Assertion failed." + ) + # pop cannot be negative + expect_error( + .check_sim_args_default( + pop = -1 + ), + "Element 1 is not >= 1." + ) + # percent_immune must be in between 0 and 1 + expect_error( + .check_sim_args_default( + percent_immune = 1.1 + ), + "Element 1 is not <= 1." + ) }) test_that(".check_time_args() returns errors", { - # Checks with .check_time_args - expect_no_error( - .check_time_args_default() - ) - # t0 cannot be negative - expect_error( - .check_time_args_default( - t0 = -1 - ), - "Element 1 is not >= 0." - ) - # tf cannot be negative - expect_error( - .check_time_args_default( - tf = -1 - ), - "Element 1 is not >= 0." - ) - # If tf is specified, generation_time must be specified too - expect_error( - .check_time_args_default( - tf_specified = TRUE, - tf = 10 - ), - "If `tf` is specified, `generation_time` must be specified too." - ) + # Checks with .check_time_args + expect_no_error( + .check_time_args_default() + ) + # t0 cannot be negative + expect_error( + .check_time_args_default( + t0 = -1 + ), + "Element 1 is not >= 0." + ) + # tf cannot be negative + expect_error( + .check_time_args_default( + tf = -1 + ), + "Element 1 is not >= 0." + ) + # If tf is specified, generation_time must be specified too + expect_error( + .check_time_args_default( + tf_specified = TRUE, + tf = 10 + ), + "If `tf` is specified, `generation_time` must be specified too." + ) }) diff --git a/tests/testthat/test-epichains.R b/tests/testthat/test-epichains.R index 42a30a7b..68d863dc 100644 --- a/tests/testthat/test-epichains.R +++ b/tests/testthat/test-epichains.R @@ -11,7 +11,7 @@ shared_args <- list( lambda = 0.9 ) -simulate_chains_default <- function(...){ +simulate_chains_default <- function(...) { default_args <- c( shared_args, generation_time = generation_time_fn diff --git a/tests/testthat/test-likelihood.R b/tests/testthat/test-likelihood.R index d3addd89..9bcd9313 100644 --- a/tests/testthat/test-likelihood.R +++ b/tests/testthat/test-likelihood.R @@ -1,138 +1,137 @@ test_that("Likelihoods can be calculated", { chains <- c(1, 1, 4, 7) set.seed(12) - expect_lt( - likelihood( - chains = chains, - statistic = "size", - offspring_dist = rpois, - lambda = 0.5 - ), - 0 - ) - expect_lt( - likelihood( - chains = chains, - statistic = "size", - offspring_dist = rpois, - lambda = 0.5, - exclude = 1 - ), - 0 - ) - expect_lt( - likelihood( - chains = chains, - statistic = "size", - offspring_dist = rpois, - lambda = 0.5, - stat_threshold = 5 - ), - 0 - ) - expect_lt( - likelihood( - chains = chains, - statistic = "size", - offspring_dist = rpois, - lambda = 0.5, - obs_prob = 0.5, - nsim_obs = 1 - ), - 0 - ) - expect_lt( - likelihood( - chains = chains, - statistic = "size", - offspring_dist = rpois, - lambda = 0.5, - stat_threshold = 5, - obs_prob = 0.5, - nsim_obs = 1 - ), - 0 - ) - expect_lt( - likelihood( - chains = chains, - statistic = "length", - offspring_dist = rbinom, - size = 1, - prob = 0.5 - ), - 0 - ) - expect_gte( - likelihood( - chains = chains, - statistic = "length", - offspring_dist = rbinom, - size = 1, - prob = 0.5, - log = FALSE - ), - 0 - ) - expect_gte( - likelihood( - chains = chains, - statistic = "length", - offspring_dist = rbinom, - size = 1, - prob = 0.5, - individual = FALSE, - log = FALSE - ), - 0 - ) - } -) + expect_lt( + likelihood( + chains = chains, + statistic = "size", + offspring_dist = rpois, + lambda = 0.5 + ), + 0 + ) + expect_lt( + likelihood( + chains = chains, + statistic = "size", + offspring_dist = rpois, + lambda = 0.5, + exclude = 1 + ), + 0 + ) + expect_lt( + likelihood( + chains = chains, + statistic = "size", + offspring_dist = rpois, + lambda = 0.5, + stat_threshold = 5 + ), + 0 + ) + expect_lt( + likelihood( + chains = chains, + statistic = "size", + offspring_dist = rpois, + lambda = 0.5, + obs_prob = 0.5, + nsim_obs = 1 + ), + 0 + ) + expect_lt( + likelihood( + chains = chains, + statistic = "size", + offspring_dist = rpois, + lambda = 0.5, + stat_threshold = 5, + obs_prob = 0.5, + nsim_obs = 1 + ), + 0 + ) + expect_lt( + likelihood( + chains = chains, + statistic = "length", + offspring_dist = rbinom, + size = 1, + prob = 0.5 + ), + 0 + ) + expect_gte( + likelihood( + chains = chains, + statistic = "length", + offspring_dist = rbinom, + size = 1, + prob = 0.5, + log = FALSE + ), + 0 + ) + expect_gte( + likelihood( + chains = chains, + statistic = "length", + offspring_dist = rbinom, + size = 1, + prob = 0.5, + individual = FALSE, + log = FALSE + ), + 0 + ) +}) test_that("likelihood() works with epichains and epichains_summary objects", { - # Simulate an object - set.seed(32) - chains_tree_eg <- simulate_chains( - n_chains = 10, - pop = 100, - percent_immune = 0, + # Simulate an object + set.seed(32) + chains_tree_eg <- simulate_chains( + n_chains = 10, + pop = 100, + percent_immune = 0, + statistic = "size", + offspring_dist = rpois, + stat_threshold = 10, + generation_time = function(n) rep(3, n), + lambda = 0.9 + ) + # Simulate an object + set.seed(32) + chains_summary_eg <- simulate_chain_stats( + n_chains = 10, + pop = 100, + percent_immune = 0, + statistic = "size", + offspring_dist = rpois, + stat_threshold = 10, + lambda = 0.9 + ) + # Use the simulated object to calculate likelihood + expect_equal( + likelihood( + chains = chains_tree_eg, statistic = "size", offspring_dist = rpois, - stat_threshold = 10, - generation_time = function(n) rep(3, n), lambda = 0.9 - ) - # Simulate an object - set.seed(32) - chains_summary_eg <- simulate_chain_stats( - n_chains = 10, - pop = 100, - percent_immune = 0, + ), + -23.538996774 + ) + # Use the simulated object to calculate likelihood + expect_equal( + likelihood( + chains = chains_summary_eg, statistic = "size", offspring_dist = rpois, - stat_threshold = 10, lambda = 0.9 - ) - # Use the simulated object to calculate likelihood - expect_equal( - likelihood( - chains = chains_tree_eg, - statistic = "size", - offspring_dist = rpois, - lambda = 0.9 - ), - -23.538996774 - ) - # Use the simulated object to calculate likelihood - expect_equal( - likelihood( - chains = chains_summary_eg, - statistic = "size", - offspring_dist = rpois, - lambda = 0.9 - ), - -23.538997 - ) + ), + -23.538997 + ) }) test_that("Likelihoods are numerically correct", { diff --git a/tests/testthat/test-simulate.R b/tests/testthat/test-simulate.R index 72940556..3eacb2cb 100644 --- a/tests/testthat/test-simulate.R +++ b/tests/testthat/test-simulate.R @@ -240,7 +240,7 @@ test_that("simulate_chain_stats throws errors", { ) }) -test_that("simulate_chain_stats is numerically correct",{ +test_that("simulate_chain_stats is numerically correct", { # Run a simulation in a small population so that # we encounter the case where we have more potential offspring than # susceptible individuals diff --git a/vignettes/interventions.Rmd b/vignettes/interventions.Rmd index b73d797d..cb8ce7cd 100644 --- a/vignettes/interventions.Rmd +++ b/vignettes/interventions.Rmd @@ -62,8 +62,10 @@ We then plot the resulting distribution of chain sizes sims[is.infinite(sims)] <- 100 # Replace infections > 99 with 100 for plotting. ggplot(data.frame(x = sims), aes(x = x)) + geom_histogram(breaks = seq(0, 100, by = 5), closed = "left") + - scale_x_continuous(breaks = c(0, 25, 50, 75, 100), - labels = c(0, 25, 50, 75, ">99")) + + scale_x_continuous( + breaks = c(0, 25, 50, 75, 100), + labels = c(0, 25, 50, 75, ">99") + ) + theme_bw() ``` @@ -84,8 +86,10 @@ sims <- simulate_chain_stats( sims[is.infinite(sims)] <- 100 # Replace infections > 99 with 100 for plotting. ggplot(data.frame(x = sims), aes(x = x)) + geom_histogram(breaks = seq(0, 100, by = 5), closed = "left") + - scale_x_continuous(breaks = c(0, 25, 50, 75, 100), - labels = c(0, 25, 50, 75, ">99")) + + scale_x_continuous( + breaks = c(0, 25, 50, 75, 100), + labels = c(0, 25, 50, 75, ">99") + ) + theme_bw() ``` @@ -122,8 +126,10 @@ sims <- simulate_chain_stats( sims[is.infinite(sims)] <- 100 # Replace infections > 99 with 100 for plotting. ggplot(data.frame(x = sims), aes(x = x)) + geom_histogram(breaks = seq(0, 100, by = 5), closed = "left") + - scale_x_continuous(breaks = c(0, 25, 50, 75, 100), - labels = c(0, 25, 50, 75, ">99")) + + scale_x_continuous( + breaks = c(0, 25, 50, 75, 100), + labels = c(0, 25, 50, 75, ">99") + ) + theme_bw() ``` @@ -147,13 +153,15 @@ This can be likened to a disease control strategy where gatherings are limited t ```{r simulate_chains_truncated} sims <- simulate_chain_stats( n_chains = 200, offspring_dist = rnbinom_truncated, stat_threshold = 99, - mu = 1.2, size = 0.5, max = 10, statistic = "size" + mu = 1.2, size = 0.5, max = 10, statistic = "size" ) sims[is.infinite(sims)] <- 100 # Replace infections > 99 with 100 for plotting. ggplot(data.frame(x = sims), aes(x = x)) + geom_histogram(breaks = seq(0, 100, by = 5), closed = "left") + - scale_x_continuous(breaks = c(0, 25, 50, 75, 100), - labels = c(0, 25, 50, 75, ">99")) + + scale_x_continuous( + breaks = c(0, 25, 50, 75, 100), + labels = c(0, 25, 50, 75, ">99") + ) + theme_bw() ``` diff --git a/vignettes/projecting_incidence.Rmd b/vignettes/projecting_incidence.Rmd index b03ad1d6..e4236343 100644 --- a/vignettes/projecting_incidence.Rmd +++ b/vignettes/projecting_incidence.Rmd @@ -23,7 +23,6 @@ knitr::opts_chunk$set( collapse = TRUE, comment = "#>" ) - ``` ## Overview @@ -137,7 +136,7 @@ first argument to determine the number of observations to sample mu <- 4.7 sgma <- 2.9 -log_mean <- log((mu^2) / (sqrt(sgma^2 + mu^2))) # log mean +log_mean <- log((mu^2) / (sqrt(sgma^2 + mu^2))) # log mean log_sd <- sqrt(log(1 + (sgma / mu)^2)) # log sd #' serial interval function @@ -308,7 +307,6 @@ head(median_daily_cases) We will now plot the individual simulation results alongside the median of the aggregated results. ```{r viz, fig.cap ="COVID-19 incidence in South Africa projected over a two week window in 2020. The light gray lines represent the individual simulations, the red line represents the median daily cases across all simulations, the black connected dots represent the observed data, and the dashed vertical line marks the beginning of the projection.", fig.width=6.0, fig.height=6} - # since all simulations may end at a different date, we will find the minimum # final date for all simulations for the purposes of visualisation. final_date <- incidence_ts_by_date %>%