Nothing
#' Simulate secondary observations from primary observations
#'
#' Simulations are done from a given trajectory of primary observations by
#' applying any given delays and observation parameters.
#'
#' In order to simulate, all parameters that are specified such as the mean and
#' standard deviation of delays or observation scaling, must be fixed.
#' Uncertain parameters are not allowed.
#'
#' A function of the same name that was previously based on a reimplementation
#' of that model in R with potentially time-varying scalings and delays is
#' available as `convolve_and_scale()
#' @param primary a data frame of primary reports (column `primary`) by date
#' (column `date`). Column `primary` must be numeric and `date` must be in
#' date format. it will be assumed that `primary` is zero on the missing
#' days.
#' @inheritParams simulate_infections
#' @inheritParams estimate_secondary
#' @importFrom checkmate assert_data_frame assert_date assert_numeric
#' assert_subset
#' @importFrom cli cli_abort
#' @return A data.table of simulated secondary observations (column `secondary`)
#' by date.
#' @export
#' @examples
#' \donttest{
#' ## load data.table to manipulate `example_confirmed` below
#' library(data.table)
#' cases <- as.data.table(example_confirmed)[, primary := confirm]
#' sim <- simulate_secondary(
#' cases,
#' delays = delay_opts(fix_parameters(example_reporting_delay)),
#' obs = obs_opts(family = "poisson")
#' )
#' }
simulate_secondary <- function(primary,
day_of_week_effect = NULL,
secondary = secondary_opts(),
delays = delay_opts(),
truncation = trunc_opts(),
obs = obs_opts(),
CrIs = c(0.2, 0.5, 0.9),
backend = "rstan",
...) {
assert_data_frame(primary, any.missing = FALSE)
assert_subset(c("date", "primary"), colnames(primary))
assert_date(primary$date)
assert_numeric(primary$primary, lower = 0)
assert_numeric(day_of_week_effect, lower = 0, null.ok = TRUE)
assert_class(secondary, "secondary_opts")
assert_class(delays, "delay_opts")
assert_class(truncation, "trunc_opts")
assert_class(obs, "obs_opts")
## create primary values for all dates modelled
all_dates <- data.table(
date = seq.Date(min(primary$date), max(primary$date), by = "day")
)
primary <- merge.data.table(all_dates, primary, by = "date", all.x = TRUE)
primary <- primary[, primary := nafill(primary, type = "const", fill = 0)]
data <- list(
n = 1,
t = nrow(primary),
h = nrow(primary),
all_dates = 0,
obs = array(integer(0)),
primary = array(primary$primary, dim = c(1, nrow(primary))),
seeding_time = 0L
)
data <- c(data, secondary)
data <- c(data, create_stan_delays(
delay = delays,
trunc = truncation
))
if (length(data$delay_params_sd) > 0 && any(data$delay_params_sd > 0)) {
cli_abort(
c(
"!" = "Cannot simulate from uncertain parameters.",
"i" = "Use {.fn fix_parameters} to set the parameters of uncertain
distributions either using the mean or a randomly sampled value."
)
)
}
data$delay_params <- array(
data$delay_params_mean, dim = c(1, length(data$delay_params_mean))
)
data$delay_params_sd <- NULL
data <- c(data, create_obs_model(
obs, dates = primary$date
))
if (data$obs_scale_sd > 0) {
cli_abort(
c(
"!" = "Cannot simulate from uncertain observation scaling.",
"i" = "Use fixed scaling instead."
)
)
}
if (data$obs_scale) {
data$frac_obs <- array(data$obs_scale_mean, dim = c(1, 1))
} else {
data$frac_obs <- array(dim = c(1, 0))
}
data$obs_scale_mean <- NULL
data$obs_scale_sd <- NULL
if (obs$family == "negbin") {
if (data$phi_sd > 0) {
cli_abort(
c(
"!" = "Cannot simulate from uncertain overdispersion.",
"i" = "Use fixed overdispersion instead."
)
)
}
data$rep_phi <- array(data$phi_mean, dim = c(1, 1))
} else {
data$rep_phi <- array(dim = c(1, 0))
}
data$phi_mean <- NULL
data$phi_sd <- NULL
## day of week effect
if (is.null(day_of_week_effect)) {
day_of_week_effect <- rep(1, data$week_effect)
}
day_of_week_effect <- day_of_week_effect / sum(day_of_week_effect)
data$day_of_week_simplex <- array(
day_of_week_effect, dim = c(1, data$week_effect)
)
# Create stan arguments
stan <- stan_opts(backend = backend, chains = 1, samples = 1, warmup = 1)
args <- create_stan_args(
stan, data = data, fixed_param = TRUE, model = "simulate_secondary",
verbose = FALSE
)
## simulate
sim <- fit_model(args, id = "simulate_secondary")
secondary <- extract_samples(sim, "sim_secondary")$sim_secondary[1, , ]
out <- data.table(date = all_dates, secondary = secondary)
return(out[])
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.