The hardware and bandwidth for this mirror is donated by dogado GmbH, the Webhosting and Full Service-Cloud Provider. Check out our Wordpress Tutorial.
If you wish to report a bug, or if you are interested in having us mirror your free-software or open-source project, please feel free to contact us at mirror[@]dogado.de.
This package uses Poisson likelihood with trend filtering penalty (a type of regularized nonparametric regression) to estimate the effective reproductive number, \(R_t\). This value roughly says “how many new infections will result from each new infection today”. Values larger than 1 indicate that an epidemic is growing while those less than 1 indicate decline.
This vignette provides a few examples to demonstrate the usage of
{rtestim}
to estimate the effective reproduction number,
\(R_t\). {rtestim}
finds a
sequence of \(R_t\), \(\{R_t\}_{t=1}^n\) of an infectious disease
by solving the following penalized Poisson regression \[\begin{equation} \label{eq:objective_fn}
\hat{\theta} = \mathop{\mathrm{argmin}}_{\theta} \frac{1}{n}
\sum_{t=1}^n \left(e^{\theta_{t}}x_t - y_t\theta_{t}\right) + \lambda
\Vert D^{(k+1)} \theta \Vert_1
\end{equation}\] where \(y_t\)
is the an epidemic signal, ideally, incident infections, but most
frequently, incident cases, on day \(t\), \(\theta_{t}
= \log(R_t)\) is the natural logarithm of \(R_t\) at time \(t\), \(D^{(k)}\) is the \(k\)-th order divided difference operator
(\(k \geq 0\)). The penalty \(\Vert D^{(k+1)} \theta \Vert_1\) imposes
smoothness on the solution and \(\lambda\) controls the level of this
smoothness, with larger \(\lambda\)
resulting in smoother estimates.
In particular \[\begin{equation} x_t = \sum_{a = 1}^m y_{t-a} w_a \end{equation}\] is the weighted sum of previous incidence at \(t\), calculated by convolving the preceding \(m\) days of new infections with the discretized serial interval distribution \(w\) of length \(m\). This delay distribution encapsulates the duration of time that a previous infection is likely to lead to future infection.
To compute \(\{R_t\}_{t=1}^n\) with
{rtestim}
, the minimal information needed is the new case
counts at days up until \(t\) and a
parametric form for the serial interval distribution (a Gamma density).
By default, \(2.5\) is used for both
the scale and shape parameters, based on the literature on contract
tracing, representing the typical delay between case onsets. This is
discretized to be supported on the integers. The order of the difference
operator, the degree of smoothness, defaults to \(3\). The sequence of smoothness penalty
\(\lambda\), if no \(\lambda\) is provided, is calculated
internally by the algorithm.
We first demonstrate the usage of the package on synthetic data, where the new daily case counts are generated from a Poisson distribution with mean parameter that roughly follows a wave. Note that the first observation must be strictly larger than 0.
set.seed(12345)
case_counts <- c(1, rpois(100, dnorm(1:100, 50, 15) * 500 + 1))
ggplot(data.frame(x = 1:101, case_counts), aes(x, case_counts)) +
geom_point(colour = "cornflowerblue") +
labs(x = "Time", y = "Case Counts")
Next, we fit the model and visualize the resulting \(\{R_t\}_{t=1}^n\):
{rtestim}
estimates a spectrum of \(\{R_t\}_{t=1}^n\)s for a range of \(\lambda\) values, where each \(\{R_t\}_{t=1}^n\) corresponds to a specific
\(\lambda\) value. If no \(\lambda\) value is supplied by the user,
{rtestim}
will automatically calculate a sequence of \(\lambda\) values. The additional parameter
nsol = 20
specifies the number of \(\lambda\)s for which the \(\{R_t\}_{t=1}^n\) is calculated
{rtestim}
also provides a cross validation procedure for
selecting the amount of smoothness to be used in the final estimate
(leave-every-k-th-out cross validation). Minimizing this metric, in
principle, balances prediction error and smoothness
(lambda.min
) though if smoother estimates are desired, one
can instead use lambda.1se
, the largest value of \(\lambda\) within one standard error of the
minimum.
The following command plots the cross validation errors for each \(\lambda\) in ascending order.
The plot above displays vertical lines that correspond to the
cross-validation scores for specific values of \(\lambda\). The blue point at the center of
each line represents the mean score for that value of \(\lambda\) across all cross-validation
folds. The top and bottom caps of each line indicate one
cross-validation standard error above and below the mean score for the
given value of \(\lambda\) across all
cross-validation folds. Two special values of \(\lambda\)’s are highlighted with dashed
lines. The one on the left represents the \(\lambda\) that gives minimum mean
cross-validated error, called lambda.min
, and the one on
the right gives the most regularized model such that the cross-validated
error is within one standard error of the minimum, called
lambda.1se
.
Users may wish to visualize the particular \(\{R_t\}_{t=1}^n\) which minimizes the cross-validation error while prioritizing smoothness.
Ideally, case counts are observed at regular intervals, such as daily
or weekly, but this is not always the case. {rtestim}
also
accommodates scenarios in which cases are reported with uneven
intervals. To demonstrate this, we generate a sequence of integers
representing the days at which we observe the case counts.
observation_incr <- rpois(101, lambda = 2)
observation_incr[observation_incr == 0] <- 1
observation_time <- cumsum(observation_incr)
We can then fit the model by passing the observation time point as
x
.
mod <- estimate_rt(observed_counts = case_counts, x = observation_time)
plot(mod) + coord_cartesian(ylim = c(0, 5))
The degree of the estimated penalized Poisson regression function \(k\) defaults to 3 for the algorithm, which corresponds to a piece-wise cubic estimate \(\{R_t\}_{t=1}^n\). To estimate \(\{R_t\}_{t=1}^n\) with piece-wise constant curves for example, use the command
Finally, we use a long history of real case counts in Canada. The data is available from opencovid.ca and the version downloaded on 4 July 2023 is included in the package. We use this data to estimate \(R_t\).
can <- estimate_rt(
observed_counts = cancovid$incident_cases,
x = cancovid$date,
korder = 2,
nsol = 20,
maxiter = 1e5
)
plot(can) + coord_cartesian(ylim = c(0.5, 2))
We also provide functionality for computing approximate confidence bands for Rt based on normal approximations and the delta method. These are intended to be fast and to provide some idea of uncertainty, but they likely don’t have guaranteed coverage.
#> An `rt_confidence_band` object.
#>
#> * type = Rt
#> * lambda = 8976.305
#> * degrees of freedom = 12
#>
#> # A tibble: 1,253 × 7
#> fit `2.5%` `10.0%` `25.0%` `75.0%` `90.0%` `97.5%`
#> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 2.05 NaN NaN NaN NaN NaN NaN
#> 2 2.03 1.95 1.97 2.00 2.05 2.08 2.10
#> 3 2.01 1.85 1.91 1.95 2.06 2.11 2.16
#> 4 1.99 1.74 1.83 1.90 2.07 2.15 2.23
#> 5 1.97 1.60 1.73 1.84 2.09 2.21 2.33
#> 6 1.95 1.44 1.62 1.77 2.12 2.28 2.45
#> 7 1.93 1.27 1.50 1.70 2.16 2.36 2.59
#> 8 1.91 1.09 1.38 1.63 2.19 2.45 2.73
#> 9 1.89 0.918 1.26 1.56 2.23 2.53 2.87
#> 10 1.88 0.756 1.14 1.49 2.26 2.61 3.00
#> # ℹ 1,243 more rows
These binaries (installable software) and packages are in development.
They may not be fully stable and should be used with caution. We make no claims about them.
Health stats visible at Monitor.