Merge branch 'main' into accumulate

.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 \

.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

.github/workflows/pkgdown.yaml

@@ -36,7 +36,7 @@ jobs:
 
       - name: Deploy to GitHub pages 🚀
         if: github.event_name != 'pull_request'
-        uses: JamesIves/[email protected].8
+        uses: JamesIves/[email protected].9
         with:
           clean: false
           branch: gh-pages

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
 

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]

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

README.Rmd

@@ -163,6 +163,7 @@ you can file an issue [here](https://github.com/epiforecasts/EpiNow2/issues). We
 
 ## Contributors
 
+
 <!-- ALL-CONTRIBUTORS-LIST:START - Do not remove or modify this section -->
 <!-- prettier-ignore-start -->
 <!-- markdownlint-disable -->
@@ -230,7 +231,8 @@ All contributions to this project are gratefully acknowledged using the [`allcon
 <a href="https://github.com/epiforecasts/EpiNow2/issues?q=is%3Aissue+author%3Alorenzwalthert">lorenzwalthert</a>, 
 <a href="https://github.com/epiforecasts/EpiNow2/issues?q=is%3Aissue+author%3Anlinton">nlinton</a>, 
 <a href="https://github.com/epiforecasts/EpiNow2/issues?q=is%3Aissue+author%3Amartinamcm">martinamcm</a>, 
-<a href="https://github.com/epiforecasts/EpiNow2/issues?q=is%3Aissue+author%3Aadrian-lison">adrian-lison</a>
+<a href="https://github.com/epiforecasts/EpiNow2/issues?q=is%3Aissue+author%3Aadrian-lison">adrian-lison</a>, 
+<a href="https://github.com/epiforecasts/EpiNow2/issues?q=is%3Aissue+author%3Azsusswein">zsusswein</a>
 
 
 
@@ -248,3 +250,4 @@ All contributions to this project are gratefully acknowledged using the [`allcon
 <!-- prettier-ignore-end -->
 <!-- ALL-CONTRIBUTORS-LIST:END -->
 
+

README.md

@@ -320,7 +320,8 @@ specification. Contributions of any kind are welcome!
 <a href="https://github.com/epiforecasts/EpiNow2/issues?q=is%3Aissue+author%3Alorenzwalthert">lorenzwalthert</a>,
 <a href="https://github.com/epiforecasts/EpiNow2/issues?q=is%3Aissue+author%3Anlinton">nlinton</a>,
 <a href="https://github.com/epiforecasts/EpiNow2/issues?q=is%3Aissue+author%3Amartinamcm">martinamcm</a>,
-<a href="https://github.com/epiforecasts/EpiNow2/issues?q=is%3Aissue+author%3Aadrian-lison">adrian-lison</a>
+<a href="https://github.com/epiforecasts/EpiNow2/issues?q=is%3Aissue+author%3Aadrian-lison">adrian-lison</a>,
+<a href="https://github.com/epiforecasts/EpiNow2/issues?q=is%3Aissue+author%3Azsusswein">zsusswein</a>
 
 ### Issue Contributors
 

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)

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)

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,16 +185,16 @@ 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.
 
 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

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)$

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,21 +59,21 @@ 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
 
 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):