diff --git a/.devcontainer/Dockerfile b/.devcontainer/Dockerfile index b1fdb5d1f..9318a7091 100644 --- a/.devcontainer/Dockerfile +++ b/.devcontainer/Dockerfile @@ -24,11 +24,10 @@ RUN apt-get update && export DEBIAN_FRONTEND=noninteractive \ libssl-dev \ libxml2-dev \ libxt-dev \ - libmagick++-6.q16-8 \ + libmagick++-6.q16-9t64 \ pandoc \ libglpk-dev \ libv8-dev \ - pandoc-citeproc \ libpng-dev \ && apt-get autoremove -y && apt-get clean -y && rm -rf /var/lib/apt/lists/* /tmp/library-scripts \ && python3 -m pip --no-cache-dir install radian pre-commit \ diff --git a/.github/workflows/build-and-push-image.yaml b/.github/workflows/build-and-push-image.yaml index 41f85db37..576c65a82 100644 --- a/.github/workflows/build-and-push-image.yaml +++ b/.github/workflows/build-and-push-image.yaml @@ -3,6 +3,7 @@ name: Create and publish a Docker image on: push: branches: ['release', 'main'] + workflow_dispatch: env: REGISTRY: ghcr.io diff --git a/.github/workflows/pkgdown.yaml b/.github/workflows/pkgdown.yaml index 8a73de6cc..4291da077 100644 --- a/.github/workflows/pkgdown.yaml +++ b/.github/workflows/pkgdown.yaml @@ -36,7 +36,7 @@ jobs: - name: Deploy to GitHub pages 🚀 if: github.event_name != 'pull_request' - uses: JamesIves/github-pages-deploy-action@v4.6.8 + uses: JamesIves/github-pages-deploy-action@v4.6.9 with: clean: false branch: gh-pages diff --git a/NEWS.md b/NEWS.md index a55b557ae..009ba40fc 100644 --- a/NEWS.md +++ b/NEWS.md @@ -14,6 +14,13 @@ fill_missing(missing = "accumulate", initial_accumulate = 7) |> estimate_infections() ``` + + - A bug was fixed where the initial growth was never estimated (i.e. the prior mean was always zero). By @sbfnk in #853 and reviewed by @seabbs. + - A bug was fixed where an internal function for applying a default cdf cutoff failed due to a difference a vector length issue. By @jamesmbaazam in #858 and reviewed by @sbfnk. + +## Documentation + +- Brought the docs on `alpha_sd` up to date with the code change from prior PR #853. By @zsusswein in #862 and reviewed by @jamesmbaazam. # EpiNow2 1.6.1 diff --git a/R/create.R b/R/create.R index 39a189afb..2ffcfaf62 100644 --- a/R/create.R +++ b/R/create.R @@ -566,7 +566,7 @@ create_stan_data <- function(data, seeding_time, if (stan_data$seeding_time > 1 && nrow(first_week) > 1) { safe_lm <- purrr::safely(stats::lm) stan_data$prior_growth <- safe_lm(log(confirm) ~ t, - stan_data = first_week + data = first_week )[[1]] stan_data$prior_growth <- ifelse(is.null(stan_data$prior_growth), 0, stan_data$prior_growth$coefficients[2] diff --git a/R/opts.R b/R/opts.R index 8597a7ea5..d56a31a32 100644 --- a/R/opts.R +++ b/R/opts.R @@ -458,7 +458,7 @@ backcalc_opts <- function(prior = c("reports", "none", "infections"), #' deviation of the Gaussian process (logged Rt in case of the renewal model, #' logged infections in case of the nonmechanistic model). #' -#' @param alpha_sd Numeric, defaults to 0.05. The standard deviation of the +#' @param alpha_sd Numeric, defaults to 0.01. The standard deviation of the #' magnitude parameter of the Gaussian process kernel. Can be tuned to adjust #' how far alpha is allowed to deviate form its prior mean (`alpha_mean`). #' @@ -509,7 +509,7 @@ gp_opts <- function(basis_prop = 0.2, ls_min = 0, ls_max = 60, alpha_mean = 0, - alpha_sd = 0.05, + alpha_sd = 0.01, kernel = c("matern", "se", "ou", "periodic"), matern_order = 3 / 2, matern_type, @@ -1098,14 +1098,14 @@ filter_opts <- function(opts, region) { #' constrained #' @keywords internal apply_default_cdf_cutoff <- function(dist, default_cdf_cutoff, cdf_cutoff_set) { - if (!is_constrained(dist) && !is.na(sd(dist))) { + if (!is_constrained(dist) && !anyNA(sd(dist))) { #nolint start: duplicate_argument_linter cli_inform( c( "i" = "Unconstrained distributon passed as a delay. ", "i" = "Constraining with default CDF cutoff {default_cdf_cutoff}.", "i" = "To silence this message, specify delay distributions - with {.var max} or {.var tolerance}." + with {.var max} or {.var default_cdf_cutoff}." ) ) #nolint end diff --git a/README.Rmd b/README.Rmd index ef16151ca..60a489d77 100644 --- a/README.Rmd +++ b/README.Rmd @@ -163,6 +163,7 @@ you can file an issue [here](https://github.com/epiforecasts/EpiNow2/issues). We ## Contributors + @@ -230,7 +231,8 @@ All contributions to this project are gratefully acknowledged using the [`allcon lorenzwalthert, nlinton, martinamcm, -adrian-lison +adrian-lison, +zsusswein @@ -248,3 +250,4 @@ All contributions to this project are gratefully acknowledged using the [`allcon + diff --git a/README.md b/README.md index b0c1158de..170f5db45 100644 --- a/README.md +++ b/README.md @@ -320,7 +320,8 @@ specification. Contributions of any kind are welcome! lorenzwalthert, nlinton, martinamcm, -adrian-lison +adrian-lison, +zsusswein ### Issue Contributors diff --git a/man/gp_opts.Rd b/man/gp_opts.Rd index af331312d..4b21c2494 100644 --- a/man/gp_opts.Rd +++ b/man/gp_opts.Rd @@ -12,7 +12,7 @@ gp_opts( ls_min = 0, ls_max = 60, alpha_mean = 0, - alpha_sd = 0.05, + alpha_sd = 0.01, kernel = c("matern", "se", "ou", "periodic"), matern_order = 3/2, matern_type, @@ -50,7 +50,7 @@ of the Gaussian process kernel. Should be approximately the expected standard deviation of the Gaussian process (logged Rt in case of the renewal model, logged infections in case of the nonmechanistic model).} -\item{alpha_sd}{Numeric, defaults to 0.05. The standard deviation of the +\item{alpha_sd}{Numeric, defaults to 0.01. The standard deviation of the magnitude parameter of the Gaussian process kernel. Can be tuned to adjust how far alpha is allowed to deviate form its prior mean (\code{alpha_mean}).} diff --git a/tests/testthat/test-create_gp_data.R b/tests/testthat/test-create_gp_data.R index 81cd895d6..b1f7cc765 100644 --- a/tests/testthat/test-create_gp_data.R +++ b/tests/testthat/test-create_gp_data.R @@ -11,7 +11,7 @@ test_that("create_gp_data returns correct default values when GP is disabled", { expect_equal(gp_data$ls_sdlog, convert_to_logsd(21, 7)) expect_equal(gp_data$ls_min, 0) expect_equal(gp_data$ls_max, 3.54, tolerance = 0.01) - expect_equal(gp_data$alpha_sd, 0.05) + expect_equal(gp_data$alpha_sd, 0.01) expect_equal(gp_data$gp_type, 2) # Default to Matern expect_equal(gp_data$nu, 3 / 2) expect_equal(gp_data$w0, 1.0) diff --git a/tests/testthat/test-gp_opts.R b/tests/testthat/test-gp_opts.R index b34a25952..cd848b75c 100644 --- a/tests/testthat/test-gp_opts.R +++ b/tests/testthat/test-gp_opts.R @@ -6,7 +6,7 @@ test_that("gp_opts returns correct default values", { expect_equal(gp$ls_sd, 7) expect_equal(gp$ls_min, 0) expect_equal(gp$ls_max, 60) - expect_equal(gp$alpha_sd, 0.05) + expect_equal(gp$alpha_sd, 0.01) expect_equal(gp$kernel, "matern") expect_equal(gp$matern_order, 3 / 2) expect_equal(gp$w0, 1.0) diff --git a/vignettes/estimate_infections.Rmd b/vignettes/estimate_infections.Rmd index 521837013..5b5e63da5 100644 --- a/vignettes/estimate_infections.Rmd +++ b/vignettes/estimate_infections.Rmd @@ -36,8 +36,8 @@ These infections are then mapped to observations via discrete convolutions with The model is initialised before the first observed data point by assuming constant exponential growth for the mean of modelled delays from infection to case report (called `seeding_time` $t_\mathrm{seed}$ in the model): \begin{align} - I_0 &\sim \mathcal{LogNormal}(I_\mathrm{obs}, 0.2) \\ - r &\sim \mathcal{Normal}(r_\mathrm{obs}, 0.2)\\ + I_0 &\sim \mathrm{LogNormal}(I_\mathrm{obs}, 0.2) \\ + r &\sim \mathrm{Normal}(r_\mathrm{obs}, 0.2)\\ I_{0 < t \leq t_\mathrm{seed}} &= I_0 \exp \left(r t \right) \end{align} @@ -79,8 +79,8 @@ More details on the mathematical form of the GP approximation and implementation or, as a specific case of a Gaussian Process, a random walk of arbitrary length $w$. \begin{align} - \log R_{t \div w} &\sim \mathcal{Normal} (R_{t \div (w - 1)}, \sigma_R)\\ - \sigma_R &\sim \mathcal{HalfNormal}(0, 0.1) + \log R_{t \div w} &\sim \mathrm{Normal} (R_{t \div (w - 1)}, \sigma_R)\\ + \sigma_R &\sim \mathrm{HalfNormal}(0, 0.1) \end{align} where $\div$ indicates interval-valued division (i.e. the floor of the division), such that for example $w=1$ indicates a daily and $w=7$ a weekly random walk. @@ -90,7 +90,7 @@ The choice of prior for the time-varying reproduction number impact run-time, sm The initial reproduction number $R_{0}$ has a log-normal prior with a given log mean $\mu_{R}$ and log standard deviation $\sigma_{R}$, calculated from a given mean (default: 1) and standard deviation (default: 1). \begin{equation} - R_0 \sim \mathcal{LogNormal}(\mu_R, \sigma_R) + R_0 \sim \mathrm{LogNormal}(\mu_R, \sigma_R) \end{equation} The simplest possible process model option is to use no time-varying prior and rely on just the intial fixed reproduction number $R_0$. @@ -185,7 +185,7 @@ The modelled counts $D_{t}$ are related to observations $C_{t}$. By default this is assumed to follow a negative binomial distribution with overdispersion $\varphi$ (alternatively it can be modelled as a Poisson, in which case $\varphi$ is not used): \begin{align} - C_t &\sim \mathcal{NegBinom}\left(\omega_{(t \mod n_\omega)}D_t, \varphi\right) + C_t &\sim \mathrm{NegBinom}\left(\omega_{(t \mod n_\omega)}D_t, \varphi\right) \end{align} where $\omega_{t \mod n_\omega}$ is a daily reporting effect of cyclicity $n_{\omega}$. If $n_{\omega}=7$ this corresponds to a day-of-the-week reporting effect. @@ -193,8 +193,8 @@ where $\omega_{t \mod n_\omega}$ is a daily reporting effect of cyclicity $n_{\o This model uses the following priors for the observation model, \begin{align} - \frac{\omega}{n_\omega} &\sim \mathcal{Dirichlet}(1, \ldots, 1) \\ - \frac{1}{\sqrt{\varphi}} &\sim \mathcal{HalfNormal}(0, 1) + \frac{\omega}{n_\omega} &\sim \mathrm{Dirichlet}(1, \ldots, 1) \\ + \frac{1}{\sqrt{\varphi}} &\sim \mathrm{HalfNormal}(0, 1) \end{align} ## Truncation diff --git a/vignettes/estimate_truncation.Rmd b/vignettes/estimate_truncation.Rmd index 69c18213d..abd012bb1 100644 --- a/vignettes/estimate_truncation.Rmd +++ b/vignettes/estimate_truncation.Rmd @@ -24,7 +24,7 @@ Given snapshots $C^{i}_{t}$ reflecting reported counts for time $t$ where $i=1\l The model assumes that final counts $D_{t}$ are related to observed snapshots via the truncation distribution such that \begin{equation} -C^{i < S)_{t}_\sim \mathcal{NegBinom}\left(Z (T_i - t | \mu_{Z}, \sigma_{Z}) D(t) + \sigma, \varphi\right) +C^{i < S)_{t}_\sim \mathrm{NegBinom}\left(Z (T_i - t | \mu_{Z}, \sigma_{Z}) D(t) + \sigma, \varphi\right) \end{equation} where $T_i$ is the date of the final observation in snapshot $i$, $Z(\tau)$ diff --git a/vignettes/gaussian_process_implementation_details.Rmd b/vignettes/gaussian_process_implementation_details.Rmd index a92091162..e272b4edc 100644 --- a/vignettes/gaussian_process_implementation_details.Rmd +++ b/vignettes/gaussian_process_implementation_details.Rmd @@ -25,7 +25,7 @@ We make use of Gaussian Processes in several places in `EpiNow2`. For example, t # Definition -The single dimension Gaussian Processes ($\mathrm{GP}_t$) we use can be written as +The single dimension Gaussian Processes ($\mathcal{GP}_t$) we use can be written as \begin{equation} \mathcal{GP}(\mu(t), k(t, t')) @@ -44,10 +44,10 @@ with the following choices available for the kernel $k$ ## Matérn 3/2 covariance kernel (the default) \begin{equation} -k(\Delta t) = \alpha^2 \left( 1 + \frac{\sqrt{3} \Delta t}{l} \right) \exp \left( - \frac{\sqrt{3} \Delta t}{l}\right) +k(\Delta t) = \alpha^2 \left( 1 + \frac{\sqrt{3} \Delta t}{\rho} \right) \exp \left( - \frac{\sqrt{3} \Delta t}{\rho}\right) \end{equation} -with $l>0$ and $\alpha > 0$ the length scale and magnitude, respectively, of the kernel. +with $\rho>0$ and $\alpha > 0$ the length scale and magnitude, respectively, of the kernel. Note that here and later we use a slightly different definition of $\alpha$ compared to Riutort-Mayol et al. [@approxGP], where this is defined as our $\alpha^2$. ## Squared exponential kernel @@ -59,13 +59,13 @@ k(\Delta t) = \alpha^2 \exp \left( - \frac{1}{2} \frac{(\Delta t^2)}{l^2} \right ## Ornstein-Uhlenbeck (Matérn 1/2) kernel \begin{equation} -k(\Delta t) = \alpha^2 \exp{\left( - \frac{\Delta t}{2 l^2} \right)} +k(\Delta t) = \alpha^2 \exp{\left( - \frac{\Delta t}{2 \rho^2} \right)} \end{equation} ## Matérn 5/2 covariance kernel \begin{equation} -k(\Delta t) = \alpha \left( 1 + \frac{\sqrt{5} \Delta t}{l} + \frac{5}{3} \left(\frac{\Delta t}{l} \right)^2 \right) \exp \left( - \frac{\sqrt{5} \Delta t}{l}\right) +k(\Delta t) = \alpha \left( 1 + \frac{\sqrt{5} \Delta t}{\rho} + \frac{5}{3} \left(\frac{\Delta t}{l} \right)^2 \right) \exp \left( - \frac{\sqrt{5} \Delta t}{\rho}\right) \end{equation} # Hilbert space approximation @@ -73,7 +73,7 @@ k(\Delta t) = \alpha \left( 1 + \frac{\sqrt{5} \Delta t}{l} + \frac{5}{3} \left( In order to make our models computationally tractable, we approximate the Gaussian Process using a Hilbert space approximation to the Gaussian Process [@approxGP], centered around mean zero. \begin{equation} -\mathcal{GP}(0, k(\Delta t)) \approx \sum_{j=1}^m \left(S_k(\sqrt{\lambda_j}) \right)^\frac{1}{2} \phi_j(t) \beta_j +\mathrm{GP}(0, k(\Delta t)) \approx \sum_{j=1}^m \left(S_k(\sqrt{\lambda_j}) \right)^\frac{1}{2} \phi_j(t) \beta_j \end{equation} with $m$ the number of basis functions to use in the approximation, which we calculate from the number of time points $t_\mathrm{GP}$ to which the Gaussian Process is being applied (rounded up to give an integer value), as is recommended [@approxGP]. @@ -91,7 +91,7 @@ and values of $\lambda_j$ given by where $L$ is a positive number termed boundary condition, and $\beta_{j}$ are regression weights with standard normal prior \begin{equation} -\beta_j \sim \mathcal{Normal}(0, 1) +\beta_j \sim \mathrm{Normal}(0, 1) \end{equation} The function $S_k(x)$ is the spectral density relating to a particular covariance function $k$. @@ -105,7 +105,7 @@ For $\nu = 3 / 2$ (the default in `EpiNow2`) this simplifies to \begin{equation} S_k(\omega) = - \alpha^2 \frac{4 \left(\sqrt{3} / \rho \right)^3}{\left(\left(\sqrt{3} / \rho\right)^2 + \omega^2\right)^2} = + \alpha^2 \frac{4 \left(\sqrt{3} / \rho \right)^3}{\left(\left(\sqrt{3} / \rho\right)^2 + \omega^2\right)^2} = \left(\frac{2 \alpha \left(\sqrt{3} / \rho \right)^{\frac{3}{2}}}{\left(\sqrt{3} / \rho\right)^2 + \omega^2} \right)^2 \end{equation} @@ -143,8 +143,8 @@ t^* = \frac{t - \frac{1}{2}t_\mathrm{GP}}{\frac{1}{2}t_\mathrm{GP}} Relevant priors are \begin{align} -\alpha &\sim \mathcal{Normal}(\mu_\alpha, \sigma_{\alpha}) \\ -\rho &\sim \mathcal{LogNormal} (\mu_\rho, \sigma_\rho)\\ +\alpha &\sim \mathrm{Normal}(\mu_\alpha, \sigma_{\alpha}) \\ +\rho &\sim \mathrm{LogNormal} (\mu_\rho, \sigma_\rho)\\ \end{align} with $\rho$ additionally constrained to be between $\rho_\mathrm{min}$ and $\rho_\mathrm{max}$, $\mu_{\rho}$ and $\sigma_\rho$ calculated from given mean $m_{\rho}$ and standard deviation $s_\rho$, and default values (all of which can be changed by the user):