From 8a46c11ed5921f7c6857becfe89f4f0b56c5745f Mon Sep 17 00:00:00 2001 From: Sebastian Funk Date: Wed, 14 Feb 2024 21:56:43 +0000 Subject: [PATCH 1/6] bring phi specification in line with scale --- R/create.R | 4 ++-- R/opts.R | 25 +++++++++++++------------ man/obs_opts.Rd | 10 ++++++---- 3 files changed, 21 insertions(+), 18 deletions(-) diff --git a/R/create.R b/R/create.R index 182781cd7..835a1f873 100644 --- a/R/create.R +++ b/R/create.R @@ -411,8 +411,8 @@ create_gp_data <- function(gp = gp_opts(), data) { create_obs_model <- function(obs = obs_opts(), dates) { data <- list( model_type = as.numeric(obs$family == "negbin"), - phi_mean = obs$phi[1], - phi_sd = obs$phi[2], + phi_mean = obs$phi$mean, + phi_sd = obs$phi$sd, week_effect = ifelse(obs$week_effect, obs$week_length, 1), obs_weight = obs$weight, obs_scale = as.integer(obs$scale$sd > 0 || obs$scale$mean != 1), diff --git a/R/opts.R b/R/opts.R index d0bf688ad..2232ed244 100644 --- a/R/opts.R +++ b/R/opts.R @@ -428,9 +428,11 @@ gp_opts <- function(basis_prop = 0.2, #' model. Custom settings can be supplied which override the defaults. #' @param family Character string defining the observation model. Options are #' Negative binomial ("negbin"), the default, and Poisson. -#' @param phi A numeric vector of length 2, defaults to 0, 1. Indicates the mean -#' and standard deviation of the normal prior used for the observation -#' process. +#' @param phi Overdispersion parameter of the reporting process, used only if +#' `familiy` is "negbin". Can be supplied either as a single numeric value +#' (fixed overdispersion) or a list with numeric elements mean (`mean`) and +#' standard deviation (`sd`) defining a normally distributed overdispersion. +#' Defaults to a list with elements `mean = 0` and `sd = 1`. #' @param weight Numeric, defaults to 1. Weight to give the observed data in the #' log density. #' @param week_effect Logical defaulting to `TRUE`. Should a day of the week @@ -470,7 +472,7 @@ gp_opts <- function(basis_prop = 0.2, #' # Scale reported data #' obs_opts(scale = list(mean = 0.2, sd = 0.02)) obs_opts <- function(family = "negbin", - phi = c(0, 1), + phi = list(mean = 0, sd = 1), weight = 1, week_effect = TRUE, week_length = 7, @@ -478,9 +480,6 @@ obs_opts <- function(family = "negbin", na = c("missing", "accumulate"), likelihood = TRUE, return_likelihood = FALSE) { - if (length(phi) != 2 || !is.numeric(phi)) { - stop("phi be numeric and of length two") - } na <- arg_match(na) if (na == "accumulate") { message( @@ -505,11 +504,13 @@ obs_opts <- function(family = "negbin", return_likelihood = return_likelihood ) - 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") + for (param in c("phi", "scale")) { + if (is.numeric(obs[[param]])) { + obs[[param]] <- list(mean = obs[[param]], sd = 0) + } + if (!(all(c("mean", "sd") %in% names(obs[[param]])))) { + stop("If specifying a ", param, " as list both a mean and sd are needed") + } } attr(obs, "class") <- c("obs_opts", class(obs)) diff --git a/man/obs_opts.Rd b/man/obs_opts.Rd index 05221d53d..c49c8bb2e 100644 --- a/man/obs_opts.Rd +++ b/man/obs_opts.Rd @@ -6,7 +6,7 @@ \usage{ obs_opts( family = "negbin", - phi = c(0, 1), + phi = list(mean = 0, sd = 1), weight = 1, week_effect = TRUE, week_length = 7, @@ -20,9 +20,11 @@ obs_opts( \item{family}{Character string defining the observation model. Options are Negative binomial ("negbin"), the default, and Poisson.} -\item{phi}{A numeric vector of length 2, defaults to 0, 1. Indicates the mean -and standard deviation of the normal prior used for the observation -process.} +\item{phi}{Overdispersion parameter of the reporting process, used only if +\code{familiy} is "negbin". Can be supplied either as a single numeric value +(fixed overdispersion) or a list with numeric elements mean (\code{mean}) and +standard deviation (\code{sd}) defining a normally distributed overdispersion. +Defaults to a list with elements \code{mean = 0} and \code{sd = 1}.} \item{weight}{Numeric, defaults to 1. Weight to give the observed data in the log density.} From 7998d04cfb6b5b684249fcccfa3adc2a13444c4c Mon Sep 17 00:00:00 2001 From: Sebastian Funk Date: Wed, 14 Feb 2024 22:36:30 +0000 Subject: [PATCH 2/6] support and deprecate legacy option --- R/opts.R | 7 +++++++ 1 file changed, 7 insertions(+) diff --git a/R/opts.R b/R/opts.R index 2232ed244..b27a67453 100644 --- a/R/opts.R +++ b/R/opts.R @@ -492,6 +492,13 @@ obs_opts <- function(family = "negbin", ) } + if (length(phi) == 2 && is.numeric(phi)) { + warning( + "Specifying `phi` as a length 2 vector is deprecated. Mean and SD ", + "should be given as list elements." + ) + phi <- list(mean = phi[1], sd = phi[2]) + } obs <- list( family = arg_match(family, values = c("poisson", "negbin")), phi = phi, From 72d4f478710d00165904f6e232d34c91404eee8d Mon Sep 17 00:00:00 2001 From: Sebastian Funk Date: Wed, 14 Feb 2024 22:36:42 +0000 Subject: [PATCH 3/6] update stan model --- inst/stan/functions/observation_model.stan | 6 ++++-- 1 file changed, 4 insertions(+), 2 deletions(-) diff --git a/inst/stan/functions/observation_model.stan b/inst/stan/functions/observation_model.stan index 6535c07f7..97ca63967 100644 --- a/inst/stan/functions/observation_model.stan +++ b/inst/stan/functions/observation_model.stan @@ -77,8 +77,10 @@ void report_lp(array[] int cases, array[] int cases_time, vector reports, obs_cases = cases; } if (model_type) { - real dispersion = 1 / pow(rep_phi[model_type], 2); - rep_phi[model_type] ~ normal(phi_mean, phi_sd) T[0,]; + real dispersion = 1 / pow(phi_sd > 0 ? rep_phi[model_type] : phi_mean, 2); + if (phi_sd > 0) { + rep_phi[model_type] ~ normal(phi_mean, phi_sd) T[0,]; + } if (weight == 1) { obs_cases ~ neg_binomial_2(obs_reports, dispersion); } else { From ad80dfe9716a9be666c255d40715ae5570408d91 Mon Sep 17 00:00:00 2001 From: Sebastian Funk Date: Wed, 14 Feb 2024 22:36:49 +0000 Subject: [PATCH 4/6] update tests --- tests/testthat/test-create_obs_model.R | 18 ++++++++++++++++-- 1 file changed, 16 insertions(+), 2 deletions(-) diff --git a/tests/testthat/test-create_obs_model.R b/tests/testthat/test-create_obs_model.R index 72cc67591..b2940a284 100644 --- a/tests/testthat/test-create_obs_model.R +++ b/tests/testthat/test-create_obs_model.R @@ -54,9 +54,23 @@ test_that("create_obs_model can be used with a custom week length", { }) test_that("create_obs_model can be used with a user set phi", { - obs <- create_obs_model(dates = dates, obs = obs_opts(phi = c(10, 0.1))) + obs <- create_obs_model( + dates = dates, obs = obs_opts(phi = list(mean = 10, sd = 0.1)) + ) expect_equal(obs$phi_mean, 10) expect_equal(obs$phi_sd, 0.1) - expect_error(obs_opts(phi = c(10))) + obs <- create_obs_model( + dates = dates, + obs = obs_opts(phi = 0.5) + ) + expect_equal(obs$phi_mean, 0.5) + expect_equal(obs$phi_sd, 0) expect_error(obs_opts(phi = c("Hi", "World"))) }) + +test_that("using a vector for phi in create_obs_model is deprecated", { + expect_warning( + create_obs_model(dates = dates, obs = obs_opts(phi = c(10, 0.1))), + "deprecated" + ) +}) From e007e5e513f3c4224458ecbbd5fa9ff351974bc2 Mon Sep 17 00:00:00 2001 From: Sebastian Funk Date: Tue, 20 Feb 2024 09:55:41 +0000 Subject: [PATCH 5/6] add news item --- NEWS.md | 1 + 1 file changed, 1 insertion(+) diff --git a/NEWS.md b/NEWS.md index bdb383159..34e16e1e2 100644 --- a/NEWS.md +++ b/NEWS.md @@ -23,6 +23,7 @@ * 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. +* Added the possibility of specifying fixed overdispersion. By @sbfnk in . ## Model changes From dc21eae8ce3983ffb205bd4ecbecf5c40e2295a1 Mon Sep 17 00:00:00 2001 From: Sebastian Funk Date: Tue, 20 Feb 2024 11:37:36 +0000 Subject: [PATCH 6/6] add PR number and reviewer --- NEWS.md | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/NEWS.md b/NEWS.md index 34e16e1e2..bdbcd6b6d 100644 --- a/NEWS.md +++ b/NEWS.md @@ -23,7 +23,7 @@ * 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. -* Added the possibility of specifying fixed overdispersion. By @sbfnk in . +* Added the possibility of specifying fixed overdispersion. By @sbfnk in #560 and reviewed by @seabbs. ## Model changes