diff --git a/NEWS.md b/NEWS.md index e68a8f9e5..bdb383159 100644 --- a/NEWS.md +++ b/NEWS.md @@ -22,6 +22,7 @@ * Changed all instances of arguments that refer to the maximum of a distribution to reflect the maximum. Previously this did, in some instance, refer to the length of the PMF. By @sbfnk in #468. * Fixed a bug in the bounds of delays when setting initial conditions. By @sbfnk in #474. * Added input checking to `estimate_infections()`, `estimate_secondary()`, `estimate_truncation()`, `simulate_infections()`, and `epinow()`. `check_reports_valid()` has been added to validate the reports dataset passed to these functions. Tests are added to check `check_reports_valid()`. As part of input validation, the various `*_opts()` functions now return subclasses of the same name as the functions and are tested against passed arguments to ensure the right `*_opts()` is passed to the right argument. For example, the `obs` argument in `estimate_secondary()` is expected to only receive arguments passed through `obs_opts()` and will error otherwise. By @jamesmbaazam in #476 and reviewed by @sbfnk and @seabbs. +* Added the possibility of specifying a fixed observation scaling. By @sbfnk in #550 and reviewed by @seabbs. ## Model changes diff --git a/R/create.R b/R/create.R index f58ad80d6..182781cd7 100644 --- a/R/create.R +++ b/R/create.R @@ -415,7 +415,9 @@ create_obs_model <- function(obs = obs_opts(), dates) { phi_sd = obs$phi[2], week_effect = ifelse(obs$week_effect, obs$week_length, 1), obs_weight = obs$weight, - obs_scale = as.numeric(length(obs$scale) != 0), + obs_scale = as.integer(obs$scale$sd > 0 || obs$scale$mean != 1), + obs_scale_mean = obs$scale$mean, + obs_scale_sd = obs$scale$sd, accumulate = obs$accumulate, likelihood = as.numeric(obs$likelihood), return_likelihood = as.numeric(obs$return_likelihood) @@ -423,14 +425,6 @@ create_obs_model <- function(obs = obs_opts(), dates) { data$day_of_week <- add_day_of_week(dates, data$week_effect) - data <- c(data, list( - obs_scale_mean = ifelse(data$obs_scale, - obs$scale$mean, 0 - ), - obs_scale_sd = ifelse(data$obs_scale, - obs$scale$sd, 0 - ) - )) return(data) } #' Create Stan Data Required for estimate_infections @@ -614,7 +608,7 @@ create_initial_conditions <- function(data) { out$bp_sd <- array(numeric(0)) out$bp_effects <- array(numeric(0)) } - if (data$obs_scale == 1) { + if (data$obs_scale_sd > 0) { out$frac_obs <- array(truncnorm::rtruncnorm(1, a = 0, b = 1, mean = data$obs_scale_mean, diff --git a/R/extract.R b/R/extract.R index d4b7fbc33..8d4cb1232 100644 --- a/R/extract.R +++ b/R/extract.R @@ -243,7 +243,7 @@ extract_parameter_samples <- function(stan_fit, data, reported_dates, value.V1 := NULL ] } - if (data$obs_scale == 1) { + if (data$obs_scale_sd > 0) { out$fraction_observed <- extract_static_parameter("frac_obs", samples) out$fraction_observed <- out$fraction_observed[, value := value.V1][, value.V1 := NULL diff --git a/R/opts.R b/R/opts.R index 2ad5f76a8..d0bf688ad 100644 --- a/R/opts.R +++ b/R/opts.R @@ -436,12 +436,13 @@ gp_opts <- function(basis_prop = 0.2, #' @param week_effect Logical defaulting to `TRUE`. Should a day of the week #' effect be used in the observation model. #' @param week_length Numeric assumed length of the week in days, defaulting to -#' 7 days. This can be modified if data aggregated over a period other than a -#' week or if data has a non-weekly periodicity. -#' @param scale List, defaulting to an empty list. Should an scaling factor be -#' applied to map latent infections (convolved to date of report). If none -#' empty a mean (`mean`) and standard deviation (`sd`) needs to be supplied -#' defining the normally distributed scaling factor. +#' 7 days. This can be modified if data aggregated over a period other than a +#' week or if data has a non-weekly periodicity. +#' @param scale Scaling factor to be applied to map latent infections (convolved +#' to date of report). Can be supplied either as a single numeric value (fixed +#' scale) or a list with numeric elements mean (`mean`) and standard deviation +#' (`sd`) defining a normally distributed scaling factor. Defaults to 1, i.e. +#' no scaling. #' @param na Character. Options are "missing" (the default) and "accumulate". #' This determines how NA values in the data are interpreted. If set to #' "missing", any NA values in the observation data set will be interpreted as @@ -473,7 +474,7 @@ obs_opts <- function(family = "negbin", weight = 1, week_effect = TRUE, week_length = 7, - scale = list(), + scale = 1, na = c("missing", "accumulate"), likelihood = TRUE, return_likelihood = FALSE) { @@ -504,13 +505,13 @@ obs_opts <- function(family = "negbin", return_likelihood = return_likelihood ) - if (length(obs$scale) != 0) { - scale_names <- names(obs$scale) - scale_correct <- "mean" %in% scale_names & "sd" %in% scale_names - if (!scale_correct) { - stop("If specifying a scale both a mean and sd are needed") - } + if (is.numeric(obs$scale)) { + obs$scale <- list(mean = obs$scale, sd = 0) + } + if (!(all(c("mean", "sd") %in% names(obs$scale)))) { + stop("If specifying a scale as list both a mean and sd are needed") } + attr(obs, "class") <- c("obs_opts", class(obs)) return(obs) } diff --git a/inst/stan/estimate_infections.stan b/inst/stan/estimate_infections.stan index 7eb128b2d..e1d4d0159 100644 --- a/inst/stan/estimate_infections.stan +++ b/inst/stan/estimate_infections.stan @@ -51,7 +51,7 @@ parameters{ array[delay_n_p] real delay_mean; // mean of delays array[delay_n_p] real delay_sd; // sd of delays simplex[week_effect] day_of_week_simplex;// day of week reporting effect - array[obs_scale] real frac_obs; // fraction of cases that are ultimately observed + array[obs_scale_sd > 0 ? 1 : 0] real frac_obs; // fraction of cases that are ultimately observed array[model_type] real rep_phi; // overdispersion of the reporting process } @@ -105,7 +105,7 @@ transformed parameters { } // scaling of reported cases by fraction observed if (obs_scale) { - reports = scale_obs(reports, frac_obs[1]); + reports = scale_obs(reports, obs_scale_sd > 0 ? frac_obs[1] : obs_scale_mean); } // truncate near time cases to observed reports if (trunc_id) { @@ -142,7 +142,7 @@ model { ); } // prior observation scaling - if (obs_scale) { + if (obs_scale_sd > 0) { frac_obs[1] ~ normal(obs_scale_mean, obs_scale_sd) T[0, 1]; } // observed reports from mean of reports (update likelihood) diff --git a/man/obs_opts.Rd b/man/obs_opts.Rd index 85aecaa7b..05221d53d 100644 --- a/man/obs_opts.Rd +++ b/man/obs_opts.Rd @@ -10,7 +10,7 @@ obs_opts( weight = 1, week_effect = TRUE, week_length = 7, - scale = list(), + scale = 1, na = c("missing", "accumulate"), likelihood = TRUE, return_likelihood = FALSE @@ -34,10 +34,11 @@ effect be used in the observation model.} 7 days. This can be modified if data aggregated over a period other than a week or if data has a non-weekly periodicity.} -\item{scale}{List, defaulting to an empty list. Should an scaling factor be -applied to map latent infections (convolved to date of report). If none -empty a mean (\code{mean}) and standard deviation (\code{sd}) needs to be supplied -defining the normally distributed scaling factor.} +\item{scale}{Scaling factor to be applied to map latent infections (convolved +to date of report). Can be supplied either as a single numeric value (fixed +scale) or a list with numeric elements mean (\code{mean}) and standard deviation +(\code{sd}) defining a normally distributed scaling factor. Defaults to 1, i.e. +no scaling.} \item{na}{Character. Options are "missing" (the default) and "accumulate". This determines how NA values in the data are interpreted. If set to diff --git a/tests/testthat/test-create_obs_model.R b/tests/testthat/test-create_obs_model.R index 4211c7dbe..72cc67591 100644 --- a/tests/testthat/test-create_obs_model.R +++ b/tests/testthat/test-create_obs_model.R @@ -6,9 +6,8 @@ test_that("create_obs_model works with default settings", { expect_equal(length(obs), 12) expect_equal(names(obs), c( "model_type", "phi_mean", "phi_sd", "week_effect", "obs_weight", - "obs_scale", "accumulate", "likelihood", "return_likelihood", - "day_of_week", "obs_scale_mean", - "obs_scale_sd" + "obs_scale", "obs_scale_mean", "obs_scale_sd", "accumulate", + "likelihood", "return_likelihood", "day_of_week" )) expect_equal(obs$model_type, 1) expect_equal(obs$week_effect, 7) @@ -16,7 +15,7 @@ test_that("create_obs_model works with default settings", { expect_equal(obs$likelihood, 1) expect_equal(obs$return_likelihood, 0) expect_equal(obs$day_of_week, c(7, 1, 2, 3, 4, 5, 6, 7, 1, 2, 3, 4, 5, 6, 7)) - expect_equal(obs$obs_scale_mean, 0) + expect_equal(obs$obs_scale_mean, 1) expect_equal(obs$obs_scale_sd, 0) }) @@ -34,6 +33,15 @@ test_that("create_obs_model can be used with a scaling", { expect_equal(obs$obs_scale_sd, 0.01) }) +test_that("create_obs_model can be used with fixed scaling", { + obs <- create_obs_model( + dates = dates, + obs = obs_opts(scale = 0.4) + ) + expect_equal(obs$obs_scale_mean, 0.4) + expect_equal(obs$obs_scale_sd, 0) +}) + test_that("create_obs_model can be used with no week effect", { obs <- create_obs_model(dates = dates, obs = obs_opts(week_effect = FALSE)) expect_equal(obs$week_effect, 1)