| Title: | Nowcast Time Series of Epidemiological Cases |
|---|---|
| Description: | Implements Bayesian nowcasting models for incident disease cases using a censored point-process approach. The reporting delay is modelled directly as a stochastic process, jointly with the epidemic dynamics, via a censored likelihood. Inference uses RTMB (CppAD autodiff + Laplace approximation). Supports arbitrary epidemic models (including HSGP, AR1, SIR), arbitrary delay families (including LogNormal, GenGamma, Dirichlet), stratified data, extreme-delay detection, and backtesting. |
| Authors: | Rodrigo Zepeda-Tello [aut, cre] (ORCID: <https://orcid.org/0000-0003-4471-5270>), Rami Yaari [aut] (ORCID: <https://orcid.org/0000-0002-8808-8937>), Matteo Perini [aut] (ORCID: <https://orcid.org/0000-0002-9465-6216>), Tal Robin [ctb] (ORCID: <https://orcid.org/0000-0002-3128-2499>), Teresa Yamana [ctb] (ORCID: <https://orcid.org/0000-0001-8349-3151>), Jeffrey Shaman [ctb], Columbia University in the City of New York [cph, fnd] |
| Maintainer: | Rodrigo Zepeda-Tello <[email protected]> |
| License: | MIT + file LICENSE |
| Version: | 2.4.0 |
| Built: | 2026-10-16 08:53:26 UTC |
| Source: | https://github.com/rodrigozepeda/diseasenowcasting |
Takes a tbl_now and chooses a model for you: it builds a grid of
candidate models (epidemic process x reporting-delay family) sized to how much
data you have, backtest()s them over several historical dates, converts the
canonical backtest to a scoringutils forecast, keeps the best one, and refits
it on the full data. The returned object
is an ordinary nowcast() result (so autoplot(), predict(), etc. work),
with the ranked scoreboard attached in its comparison slot.
auto_nowcast( data, metric = "wis", relative_score = TRUE, tie_break = c("epidemic_priority", "fastest"), type = c("auto", "two_stage", "one_stage"), sir = NULL, ar = NULL, hsgp = NULL, delays = NULL, likelihood = nb_likelihood(), models = NULL, n_dates = 6L, n_draws_select = 500L, n_draws = 2000L, K = 25L, K_select = 10L, min_ar = 15L, min_hsgp = 30L, now = NULL, seed = sample.int(.Machine$integer.max, 1), verbose = TRUE, ... )auto_nowcast( data, metric = "wis", relative_score = TRUE, tie_break = c("epidemic_priority", "fastest"), type = c("auto", "two_stage", "one_stage"), sir = NULL, ar = NULL, hsgp = NULL, delays = NULL, likelihood = nb_likelihood(), models = NULL, n_dates = 6L, n_draws_select = 500L, n_draws = 2000L, K = 25L, K_select = 10L, min_ar = 15L, min_hsgp = 30L, now = NULL, seed = sample.int(.Machine$integer.max, 1), verbose = TRUE, ... )
data |
A |
metric |
A single score column produced by |
relative_score |
Logical. When |
tie_break |
How to break effectively equal selection scores.
|
type |
Stage strategy used for both the backtest and the final fit:
|
sir, ar, hsgp
|
Optional epidemic-process components (e.g.
|
delays |
A list of delay components to compare. Default:
|
likelihood |
Either a single likelihood used for every candidate
(default |
models |
Optional |
n_dates |
Number of historical dates to backtest over (default 6). |
n_draws_select |
Posterior draws during the selection backtest (default 500 – kept small for speed). |
n_draws |
Posterior draws for the final fit of the winning model (default 2000). |
K |
Delay imputations for the final two-stage fit of the winning model (default 25). |
K_select |
Delay imputations during the selection backtest (default
10 – kept small for speed, like |
min_ar, min_hsgp
|
Series-length thresholds (in event-times) at which AR(1) and HSGP become candidates (defaults 15 and 30). |
now |
As-of date for the final fit (default: the |
seed |
RNG seed. |
verbose |
Print progress and the chosen model (default |
... |
Passed through to |
Candidate epidemic processes are chosen by series length (max_time, the
number of event-times): a process becomes a candidate as soon as the series is
long enough to support it (SIR needs the least data, the HSGP the most) and is
never dropped for being too long, so the comparison always spans every
process the data can support. With the default thresholds:
max_time < min_ar -> compares {SIR};
min_ar <= max_time < min_hsgp -> compares {SIR, AR(1)};
max_time >= min_hsgp -> compares {SIR, AR(1), HSGP}.
Any process you pass explicitly via sir / ar / hsgp is always included
(regardless of length), which is how you make a prior compete: e.g. pass
sir = sir_epidemic(R0 = lognormal_prior(log(3), 0.2)) and the SIR candidate
will use that R0 prior throughout the comparison.
Robustness. A candidate that fails to converge on a backtest date simply
drops out of the comparison there (it never aborts the search), and candidates
are scored on the common set of dates where they all produced a forecast so a
model cannot "win" on a lucky subset. The winner is then refit on the full
data; if that refit fails, auto_nowcast() falls through to the next-best
candidate (and so on), so it converges whenever any candidate would.
Candidate delays default to LogNormal, Generalized-Gamma and Dirichlet;
override with delays.
Speed. The grid is backtested with a fast configuration
(n_draws_select posterior draws over n_dates dates spread across the
history); only the winning model is refit with the full n_draws. Backtesting
is the expensive step – set a future::plan() (e.g.
future::plan(multisession)) for parallel speed-up.
A diseasenowcasting subclass of tbl.now::tbl_nowcast (as from
nowcast()) for the selected model, with the model-selection scoreboard
retained on the diseasenowcasting subclass and its native fit:
list(scores, chosen, metric, relative_score, tie_break, timings, max_time).
diseasenowcasting_workflows, nowcast(), backtest(), scoringutils::score(),
scoringutils::add_relative_skill()
library(tbl.now) data(denguedat) # A short window keeps this example quick (auto_nowcast fits a whole grid): dn <- subset(denguedat, onset_week >= as.Date("1990-06-01") & onset_week <= as.Date("1990-12-01")) tn <- tbl_now(dn, event_date = onset_week, report_date = report_week, data_type = "linelist", verbose = FALSE) # Backtesting the grid is the expensive step -- uncomment to run candidates # in parallel (then restore sequential afterwards): # future::plan(future::multisession, workers = 4) # Compare a couple of delays; make the SIR candidate use a custom R0 prior: nc <- auto_nowcast(tn, sir = sir_epidemic(R0 = lognormal_prior(log(2), 0.3)), delays = list(lognormal_delay(), dirichlet_delay()), n_dates = 2, n_draws_select = 150, n_draws = 300, temporal_effects = "none") # future::plan(future::sequential) best_model_name(nc) # the winning model's label comparison_scores(nc) # the ranked scoreboard best_score(nc) # just the winner's row selection_metric(nc) # which metric chose it winner <- best_model(nc) # the model() object, to reuse elsewherelibrary(tbl.now) data(denguedat) # A short window keeps this example quick (auto_nowcast fits a whole grid): dn <- subset(denguedat, onset_week >= as.Date("1990-06-01") & onset_week <= as.Date("1990-12-01")) tn <- tbl_now(dn, event_date = onset_week, report_date = report_week, data_type = "linelist", verbose = FALSE) # Backtesting the grid is the expensive step -- uncomment to run candidates # in parallel (then restore sequential afterwards): # future::plan(future::multisession, workers = 4) # Compare a couple of delays; make the SIR candidate use a custom R0 prior: nc <- auto_nowcast(tn, sir = sir_epidemic(R0 = lognormal_prior(log(2), 0.3)), delays = list(lognormal_delay(), dirichlet_delay()), n_dates = 2, n_draws_select = 150, n_draws = 300, temporal_effects = "none") # future::plan(future::sequential) best_model_name(nc) # the winning model's label comparison_scores(nc) # the ranked scoreboard best_score(nc) # just the winner's row selection_metric(nc) # which metric chose it winner <- best_model(nc) # the model() object, to reuse elsewhere
backtest() translates native model() specifications into labelled
tbl.now::engine_diseasenowcasting() specifications and delegates the full
retrospective workflow to tbl.now::nowcast_backtest(). The returned object
is therefore the common nowcast_backtest result used by tbl.now for
tidying, forecast conversion, scoring, weighting, and ensembling.
backtest( data, models = diseasenowcasting::model(), dates = NULL, type = c("two_stage", "one_stage", "auto"), horizon = NULL, n_dates = 4L, n_draws = 1000L, K = 25L, np_spread = 1, seed = NULL, keep_draws = FALSE, on_error = c("warn", "abort"), verbose = TRUE, truth_axis = NULL, truth_type = NULL, quantile_levels = tbl.now::nowcast_quantile_levels(), ... )backtest( data, models = diseasenowcasting::model(), dates = NULL, type = c("two_stage", "one_stage", "auto"), horizon = NULL, n_dates = 4L, n_draws = 1000L, K = 25L, np_spread = 1, seed = NULL, keep_draws = FALSE, on_error = c("warn", "abort"), verbose = TRUE, truth_axis = NULL, truth_type = NULL, quantile_levels = tbl.now::nowcast_quantile_levels(), ... )
data |
A tbl.now::tbl_now holding the full data, including observations that arrived after the retrospective nowcast dates. |
models |
A |
dates |
Retrospective nowcast origins, passed as |
type |
|
horizon |
Number of time units of hindsight used by |
n_dates |
Number of automatic retrospective origins. Ignored when
|
n_draws |
Posterior draws per fit. |
K, np_spread
|
Native two-stage fitting controls passed to every engine. |
seed |
Optional base seed. |
keep_draws |
Whether the canonical backtest retains posterior draws. |
on_error |
Either |
verbose |
Whether to report progress. |
truth_axis, truth_type
|
Canonical scoring truth controls passed to
|
quantile_levels |
Quantile probabilities requested from every model. |
... |
Additional arguments passed to every |
A tbl.now::nowcast_backtest object.
diseasenowcasting_workflows for the native/common ownership
boundary; tbl.now::score_nowcast(), tbl.now::nowcast_weights(),
tbl.now::nowcast_ensemble(), fit_check()
model() object from a nowcastReturns the fitted nowcast's model() specification. For an auto_nowcast()
result this is the selected model, so you can reuse it elsewhere, e.g.
nowcast(other_data, model = best_model(nc)) or pass it to backtest().
best_model(nc)best_model(nc)
nc |
A |
A model_class object.
auto_nowcast(), best_model_name()
auto_nowcast()
Name of the model chosen by auto_nowcast()
best_model_name(nc)best_model_name(nc)
nc |
A |
The winning model's label, a string of the form
"epidemic/likelihood/delay" (e.g. "HSGP/nb/Dirichlet").
auto_nowcast(), best_model(), comparison_scores()
auto_nowcast() choseThe single comparison_scores() row belonging to the winning model,
including its scoringutils metrics, selection score, and retrospective fit
timing, rather than the whole table.
best_score(nc)best_score(nc)
nc |
A |
A one-row data.frame.
auto_nowcast(), comparison_scores(), best_model_name()
auto_nowcast()
The ranked table of candidate models that auto_nowcast() backtested, one row
per model, best-first by the selected raw or relative scoringutils metric.
In addition to scoringutils metrics, it includes selection_score, median
and total retrospective fit seconds, successful-fit count, epidemic priority,
and original grid order.
comparison_scores(nc)comparison_scores(nc)
nc |
A |
A data.frame, one row per candidate model.
auto_nowcast(), best_score(), selection_timings(),
scoringutils::score()
Daily case counts of COVID-19 from Colombia's national epidemiological
surveillance system (INS), aggregated by notification date, diagnosis date,
and sex. Each row is a unique combination of these three variables together
with the number of cases n.
covid_colombiacovid_colombia
A data frame with 35,501 rows and 4 variables:
Date. The event date – when the case was
notified / symptom onset was recorded.
Date. The report date – when the laboratory
diagnosis was registered in the national system.
character. Biological sex of the case: "Female" or
"Male".
integer. Number of cases with this
(notification_date, diagnosis_date, sex) combination.
In the nowcasting context the event date is notification_date (when the
case symptom onset was recorded in the system) and the report date is
diagnosis_date (when the laboratory result was entered). The delay between
the two reflects the time from symptom onset to laboratory confirmation and
data entry – the quantity the nowcasting model estimates and corrects for.
The dataset covers the full first three years of the Colombian epidemic (2020-03-02 to 2023-03-03) and was used in the NobBS comparison study to benchmark the diseasenowcasting / diseasenowcast2 engine.
Instituto Nacional de Salud (INS), Colombia. Data accessed via the SIVIGILA open-data platform and pre-processed for the diseasenowcasting benchmarking study.
# Build a stratified tbl_now (event = notification, report = diagnosis, # strata = sex). Pipe through temporal effects for day-of-week covariates. if (requireNamespace("tbl.now", quietly = TRUE)) { tn <- tbl.now::tbl_now( covid_colombia, event_date = notification_date, report_date = diagnosis_date, strata = sex, case_count = n, data_type = "count-incidence", verbose = FALSE ) print(tn) }# Build a stratified tbl_now (event = notification, report = diagnosis, # strata = sex). Pipe through temporal effects for day-of-week covariates. if (requireNamespace("tbl.now", quietly = TRUE)) { tn <- tbl.now::tbl_now( covid_colombia, event_date = notification_date, report_date = diagnosis_date, strata = sex, case_count = n, data_type = "count-incidence", verbose = FALSE ) print(tn) }
Configures models for revision streams that publish cumulative levels. The
target is finite-horizon database retention C_t(H), not biological truth.
The retraction mechanism is the collapsed kernel
h_R(l) = retraction_mass * g_R(l); it does not separately identify a truth
probability and a conditional revision-delay law.
cumulative_process( observation = c("hurdle_ztnb", "hurdle_ztpoisson", "cumulative"), retraction_delay = lognormal_delay(), settlement = 26L, retraction_mass = beta_prior(1.5, 20), movement_intercept = normal_prior(-1, 2), movement_age = normal_prior(0, 1), movement_previous = normal_prior(0, 1), magnitude_size = NULL )cumulative_process( observation = c("hurdle_ztnb", "hurdle_ztpoisson", "cumulative"), retraction_delay = lognormal_delay(), settlement = 26L, retraction_mass = beta_prior(1.5, 20), movement_intercept = normal_prior(-1, 2), movement_age = normal_prior(0, 1), movement_previous = normal_prior(0, 1), magnitude_size = NULL )
observation |
Observation composite likelihood. |
retraction_delay |
Parametric delay family for retraction ages |
settlement |
Positive integer settlement horizon |
retraction_mass |
Prior or fixed value in |
movement_intercept, movement_age, movement_previous
|
Priors or fixed
values for the bounded movement-probability regression. Set
|
magnitude_size |
Positive prior or fixed value for the ZTNB magnitude
size. It is used only by |
A cumulative_process_class object for
model(cumulative = ).
cumulative_process() cumulative_process(observation = "cumulative", settlement = 52L) cumulative_process(observation = "hurdle_ztpoisson", settlement = 6L)cumulative_process() cumulative_process(observation = "cumulative", settlement = 52L) cumulative_process(observation = "hurdle_ztpoisson", settlement = 6L)
custom_delay( cdf, log_cdf = function(theta) function(d) log(cdf(theta)(d)), log_survival = function(theta) function(d) log(1 - cdf(theta)(d)), priors = list(), name = "Custom", param_names = NULL, inits = NULL, num_delay_seasons = 1L, season_distribution = std_normal_prior() )custom_delay( cdf, log_cdf = function(theta) function(d) log(cdf(theta)(d)), log_survival = function(theta) function(d) log(1 - cdf(theta)(d)), priors = list(), name = "Custom", param_names = NULL, inits = NULL, num_delay_seasons = 1L, season_distribution = std_normal_prior() )
cdf |
A function |
log_cdf |
Optional function |
log_survival |
Optional function |
priors |
A list with one element per parameter. Each element is either a
|
name |
A display name for the distribution (used in print output). |
param_names |
Optional character vector naming each parameter (used in
diagnostics); its length sets the number of parameters if |
inits |
Numeric vector of initial values for the parameters (on the
unconstrained scale passed to the functions); its length sets the number of
parameters if |
num_delay_seasons |
Number of periodic delay seasons. Default 1. |
season_distribution |
Prior for the delay-season effects. |
Allows any reporting-delay distribution that can be expressed with
RTMB-traceable functions. You supply the cumulative distribution function
cdf; log_cdf and log_survival are optional and default to the obvious
transforms of cdf. Each is given as a function of the parameter vector
theta that returns a function of the delay d:
cdf(theta) -> function(d) returning F(d) log_cdf(theta) -> function(d) returning log F(d) log_survival(theta) -> function(d) returning log(1 - F(d))
Inside the RTMB tape theta is an advector, so every operation must be
AD-traceable (see the RTMB traceability section).
A custom_delay_class object (a delay_process_class).
All operations inside cdf / log_cdf / log_survival must be differentiable
under RTMB's CppAD tape:
Allowed: +, -, *, /, exp, log, sqrt, pnorm,
pgamma, lgamma, abs, sum, scalar multiplication, fixed-length
for loops (length must not depend on parameter values).
Not allowed: if/ifelse on parameter values, pmax/pmin on AD
types (use (x + abs(x)) / 2 for pmax(x, 0)), external solvers,
anything non-differentiable.
Use validate_custom_delay() to check your functions before fitting.
validate_custom_delay(), delay_process
# Weibull delay: theta = c(log_shape, log_scale). # Only `cdf` is required; we also supply `log_survival` (exact in the tail). weibull_cdf <- function(theta) { shape <- exp(theta[1]); scale <- exp(theta[2]) function(d) 1 - exp(-(d / scale)^shape) } weibull_log_survival <- function(theta) { shape <- exp(theta[1]); scale <- exp(theta[2]) function(d) -(d / scale)^shape } custom_delay( cdf = weibull_cdf, log_survival = weibull_log_survival, priors = list(normal_prior(0, 1), normal_prior(log(7), 1)), name = "Weibull", param_names = c("log_shape", "log_scale"), inits = c(0, log(7)) )# Weibull delay: theta = c(log_shape, log_scale). # Only `cdf` is required; we also supply `log_survival` (exact in the tail). weibull_cdf <- function(theta) { shape <- exp(theta[1]); scale <- exp(theta[2]) function(d) 1 - exp(-(d / scale)^shape) } weibull_log_survival <- function(theta) { shape <- exp(theta[1]); scale <- exp(theta[2]) function(d) -(d / scale)^shape } custom_delay( cdf = weibull_cdf, log_survival = weibull_log_survival, priors = list(normal_prior(0, 1), normal_prior(log(7), 1)), name = "Weibull", param_names = c("log_shape", "log_scale"), inits = c(0, log(7)) )
Lets you supply any RTMB-traceable function intensity_fn(theta) that
returns the log expected-incidence trajectory log_mean[n_time x n_strata]
as the latent epidemic process for the nowcast. This makes the framework
epidemic-process agnostic: random walks, ODE models, regression surfaces,
and anything else that can be written in AD-safe arithmetic are all valid.
custom_epidemic( intensity_fn, priors = list(), name = "Custom", param_names = character(0), inits = numeric(0) )custom_epidemic( intensity_fn, priors = list(), name = "Custom", param_names = character(0), inits = numeric(0) )
intensity_fn |
A function |
priors |
A list with one element per parameter in |
name |
Character label shown in print and diagnostic output. |
param_names |
Character vector naming the parameters; its length sets the
number of parameters if |
inits |
Numeric vector of starting values; its length sets the number of
parameters if |
A custom_epidemic_class object (subclass of
epidemic_process_class).
validate_custom_epidemic(), epidemic_process
# Pure random walk on log-incidence (max_time = 20) # Uses cumsum() — no [<- assignment needed, so fully AD-safe. max_time <- 20L rw_fn <- function(theta) { log_mu0 <- theta[1] sigma_rw <- exp(theta[2]) eps <- theta[3:(2L + max_time)] lm <- log_mu0 + cumsum(sigma_rw * eps) matrix(lm, max_time, 1L) } # n_params is inferred (here 2 + max_time) from the priors list: custom_epi <- custom_epidemic( rw_fn, priors = c(list(normal_prior(2, 1), normal_prior(-2, 0.5)), rep(list(std_normal_prior()), max_time)), name = "RandomWalk", param_names = c("log_mu0", "log_sigma", paste0("eps_", seq_len(max_time))), inits = c(2, -2, rep(0, max_time)) )# Pure random walk on log-incidence (max_time = 20) # Uses cumsum() — no [<- assignment needed, so fully AD-safe. max_time <- 20L rw_fn <- function(theta) { log_mu0 <- theta[1] sigma_rw <- exp(theta[2]) eps <- theta[3:(2L + max_time)] lm <- log_mu0 + cumsum(sigma_rw * eps) matrix(lm, max_time, 1L) } # n_params is inferred (here 2 + max_time) from the priors list: custom_epi <- custom_epidemic( rw_fn, priors = c(list(normal_prior(2, 1), normal_prior(-2, 0.5)), rep(list(std_normal_prior()), max_time)), name = "RandomWalk", param_names = c("log_mu0", "log_sigma", paste0("eps_", seq_len(max_time))), inits = c(2, -2, rep(0, max_time)) )
Build the default prior bundle for an RTMB nowcast model
default_priors(mod, data = NULL, ...)default_priors(mod, data = NULL, ...)
mod |
A |
data |
Optional prepared-data list from |
... |
Per-key overrides (e.g. |
The default prior for each parameter is documented on the component
constructor: see the Default priors section of epidemic_process
(HSGP / AR(1) / SIR), delay_process (LogNormal / Gamma / GenGamma /
Dirichlet) and likelihood (the NB overdispersion phi). To see what a
set of priors implies for the epidemic curve before fitting, use
nowcast(prior_only = TRUE).
A named list of prior specs.
epidemic_process, delay_process, likelihood
Specify the reporting-delay distribution. Parameter slots accept a fixed
numeric, a prior_class (e.g. normal_prior()), or numeric(0) for the
default prior.
lognormal_delay( mu = numeric(0), sigma = numeric(0), num_delay_seasons = 1L, season_distribution = std_normal_prior() ) gamma_delay( shape = numeric(0), rate = numeric(0), num_delay_seasons = 1L, season_distribution = std_normal_prior() ) generalized_gamma_delay( mu = numeric(0), sigma = numeric(0), Q = numeric(0), num_delay_seasons = 1L, season_distribution = std_normal_prior() ) dirichlet_delay(alpha = numeric(0), bins = numeric(0))lognormal_delay( mu = numeric(0), sigma = numeric(0), num_delay_seasons = 1L, season_distribution = std_normal_prior() ) gamma_delay( shape = numeric(0), rate = numeric(0), num_delay_seasons = 1L, season_distribution = std_normal_prior() ) generalized_gamma_delay( mu = numeric(0), sigma = numeric(0), Q = numeric(0), num_delay_seasons = 1L, season_distribution = std_normal_prior() ) dirichlet_delay(alpha = numeric(0), bins = numeric(0))
mu |
Log-mean intercept ( |
sigma |
Log-scale / SD parameter > 0. |
num_delay_seasons |
Number of periodic delay seasons. Default 1. |
season_distribution |
Prior for the delay-season effects. |
shape, rate
|
Gamma delay parameters (the |
Q |
GenGamma shape ( |
alpha |
Dirichlet concentration (scalar broadcast to all bins). |
bins |
Dirichlet: number of explicit delay bins (geometric tail beyond). |
A delay_process_class object.
When a prior argument is left empty, default_priors() supplies these
defaults. The delay mean prior is data-informed: a normal_prior() on
the log scale, centred at the log of the median observed delay.
LogNormal: mu ~ data-informed normal_prior(log median delay);
sigma ~ gamma_prior(2, 2).
Gamma: shape ~ data-informed normal_prior(...);
rate ~ gamma_prior(2, 2 / sd) (data-informed SD).
Generalised Gamma: mu ~ data-informed normal_prior(...);
sigma ~ gamma_prior(2, 0.1); Q (shape) ~ normal_prior(0, 0.5).
Dirichlet: alpha ~ a per-bin concentration vector, data-informed
from the empirical delay pmf (0.05 + (bins+1) * pmf), else
rep(1, bins + 1).
lognormal_delay() lognormal_delay(mu = normal_prior(log(7), 0.5)) gamma_delay() generalized_gamma_delay(Q = 1) dirichlet_delay(alpha = 1, bins = 21)lognormal_delay() lognormal_delay(mu = normal_prior(log(7), 0.5)) gamma_delay() generalized_gamma_delay(Q = 1) dirichlet_delay(alpha = 1, bins = 21)
diseasenowcasting owns the statistical model and RTMB fit;
tbl.now::tbl_nowcast owns the common result grammar. The public result of
nowcast() and auto_nowcast() therefore supports both layers without an
explicit conversion.
Use model(), nowcast(), and auto_nowcast() when you need this package's
epidemic, delay, likelihood, revision, or cumulative-process components.
Use fit_check() and nowcast_diagnostic() for RTMB-specific convergence and
fit diagnostics. These functions automatically unwrap the retained native
fit from the common result.
backtest() is a native model-specification convenience, not a second
evaluation system: it translates each model() into a
tbl.now::engine_diseasenowcasting() and returns
tbl.now::nowcast_backtest() directly.
Use tbl.now::run_nowcast() with
tbl.now::engine_diseasenowcasting() when diseasenowcasting is one of several
engines. Once a fit exists, use the common methods for tidy(), autoplot(),
predictive scoring, forecast conversion, ensembling, and persistence.
In particular, a backtest can be passed directly to
scoringutils::as_forecast_quantile(),
scoringutils::as_forecast_point(), or, when draws were retained,
scoringutils::as_forecast_sample(). Relative WIS and additional scores then
come from scoringutils::score() and
scoringutils::add_relative_skill(), rather than a package-local scorer.
## Not run: fit <- nowcast(data, model = model()) autoplot(fit) # common result plot fit_check(fit) # native RTMB diagnostics bt <- backtest(data, models = list(default = model())) relative <- bt |> scoringutils::as_forecast_quantile() |> scoringutils::score() |> scoringutils::add_relative_skill(metric = "wis") ## End(Not run)## Not run: fit <- nowcast(data, model = model()) autoplot(fit) # common result plot fit_check(fit) # native RTMB diagnostics bt <- backtest(data, models = list(default = model())) relative <- bt |> scoringutils::as_forecast_quantile() |> scoringutils::score() |> scoringutils::add_relative_skill(metric = "wis") ## End(Not run)
Returns hex colour codes matching the package's visual identity (derived from the hex-sticker logo: dark charcoal background, sage green, burnt orange accent).
dn_palette(n = NULL)dn_palette(n = NULL)
n |
Number of colours to return (1-8). If |
A named character vector of hex colour codes.
Specify the latent epidemic process. Parameter slots accept a fixed
numeric, a prior_class, or numeric(0) for the default prior. The
log-incidence mean intercept is inherited from the likelihood (mu).
hsgp_epidemic( alpha = numeric(0), ell = numeric(0), gp_kernel = "matern32", gp_basis = "dirichlet", num_basis = 0, tmax_model = 0 ) ar1_epidemic(phi = numeric(0), sigma = numeric(0), error = numeric(0)) sir_epidemic( R0 = numeric(0), gamma = numeric(0), N_eff = numeric(0), N_pop = 10000, use_beta_rw_trend = TRUE )hsgp_epidemic( alpha = numeric(0), ell = numeric(0), gp_kernel = "matern32", gp_basis = "dirichlet", num_basis = 0, tmax_model = 0 ) ar1_epidemic(phi = numeric(0), sigma = numeric(0), error = numeric(0)) sir_epidemic( R0 = numeric(0), gamma = numeric(0), N_eff = numeric(0), N_pop = 10000, use_beta_rw_trend = TRUE )
alpha |
HSGP GP amplitude prior (> 0). |
ell |
HSGP GP length-scale prior (> 0). |
gp_kernel |
HSGP kernel: |
gp_basis |
HSGP eigenbasis: |
num_basis |
HSGP basis count; |
tmax_model |
HSGP time normalisation; |
phi |
AR(1) autocorrelation prior in (-1, 1). |
sigma |
AR(1) innovation SD prior (> 0). |
error |
AR(1) standardised innovation prior. |
R0 |
SIR basic reproduction number prior (> 0). |
gamma |
SIR recovery rate prior in (0, 1). |
N_eff |
SIR effective susceptible fraction prior in (0, 1). |
N_pop |
SIR total population (default 10000). |
use_beta_rw_trend |
SIR: beta follows an AR(1) walk if TRUE (default). |
An epidemic_process_class object.
When a prior argument is left empty, default_priors() supplies these
defaults (see also nowcast(prior_only = TRUE) to visualise them):
HSGP (hsgp_epidemic):
alpha (GP amplitude): half_normal_prior(0, 1)
ell (GP length-scale): inv_gamma_prior(3, 1)
AR(1) (ar1_epidemic):
phi (autocorrelation): std_normal_prior()
sigma (innovation SD): exponential_prior(100)
innovations: std_normal_prior()
SIR (sir_epidemic):
R0: lognormal_prior(log(2), 0.5)
gamma (recovery rate): lognormal_prior(log(1/5), 0.5)
N_eff (susceptible fraction): beta_prior(2, 5)
The log-incidence intercept comes from the likelihood (mu), defaulting to a
data-informed normal_prior() centred at the log median daily count.
hsgp_epidemic() hsgp_epidemic(gp_kernel = "sq_exp", num_basis = 20) ar1_epidemic(phi = 0.9) sir_epidemic(R0 = 2.5, use_beta_rw_trend = FALSE)hsgp_epidemic() hsgp_epidemic(gp_kernel = "sq_exp", num_basis = 20) ar1_epidemic(phi = 0.9) sir_epidemic(R0 = 2.5, use_beta_rw_trend = FALSE)
update()
Returns a tidy data.frame with one row per flagged surprise – currently the
reporting delays that arrived later than the model expects (surprise = "delay", direction = "long") – together with their tail probability and
the level used. Returns NULL when nothing was flagged (or surprise was
not computed).
extreme_values(nc)extreme_values(nc)
nc |
A |
A data.frame of flagged surprises, or NULL.
Optimises the negative log-posterior built from model + data. For
delay_only data this fits the reporting-delay process alone (no epidemic);
the joint epidemic fit is added in later phases.
fit( model, data, priors = NULL, init = NULL, control = list(iter.max = 500, eval.max = 1000, rel.tol = 1e-09), warn = TRUE )fit( model, data, priors = NULL, init = NULL, control = list(iter.max = 500, eval.max = 1000, rel.tol = 1e-09), warn = TRUE )
model |
A |
data |
Prepared-data list from |
priors |
Optional prior bundle; defaults to |
init |
Optional named init list. |
control |
|
warn |
If |
A list with par (named estimates), obj, opt, data, priors,
model, convergence, and (delay-only) delay_mu / delay_sigma.
Predictive accuracy belongs to tbl.now::score_nowcast() and
tbl.now::nowcast_backtest(). fit_check() deliberately reports only
diagnostics specific to the RTMB optimization performed by
diseasenowcasting.
fit_check(object, warn = TRUE)fit_check(object, warn = TRUE)
object |
A result from |
warn |
If |
A data frame with one row per retained RTMB fit and columns fit,
rung, convergence, objective, raw and projected gradients, quadratic
objective gap, Hessian status, any Laplace-precision regularization used
for prediction, overall status, and diagnostic reasons.
diseasenowcasting_workflows for the distinction between native fit
diagnostics and predictive scoring; nowcast_diagnostic(),
tbl.now::score_nowcast(), tbl.now::nowcast_backtest()
Hard-fix a parameter in a prior bundle (treat as data, drop from estimation)
fix_param(priors, key, value)fix_param(priors, key, value)
priors |
A prior bundle from |
key |
Parameter key to fix. |
value |
Fixed numeric value. |
The modified prior bundle.
max_time for custom epidemic processesInfers max_time from a model and a tbl_now. max_time is the number of
event-time points the model spans: the epidemic process runs from t = 0 to
t = max_time - 1.
infer_max_time(data, model = diseasenowcasting::model())infer_max_time(data, model = diseasenowcasting::model())
data |
A |
model |
A |
This function is useful for building custom epidemic processes whose intensity function loops over time and therefore needs to know the number of event-times before it is written (e.g. an SIR recursion or a random walk with one innovation per time point).
The number of event-time points in the model (an integer).
library(tbl.now) library(dplyr) mpox_pooled <- mpoxdat |> ungroup() |> summarise(n = sum(n), .by = c(dx_date, dx_report_date)) mpox_tn <- tbl_now(mpox_pooled, event_date = dx_date, report_date = dx_report_date, case_count = n, data_type = "count-incidence", verbose = FALSE) # The result is 316: the model spans from time t = 0 (min event_num) # to time t = 315 (max report_num). infer_max_time(mpox_tn)library(tbl.now) library(dplyr) mpox_pooled <- mpoxdat |> ungroup() |> summarise(n = sum(n), .by = c(dx_date, dx_report_date)) mpox_tn <- tbl_now(mpox_pooled, event_date = dx_date, report_date = dx_report_date, case_count = n, data_type = "count-incidence", verbose = FALSE) # The result is 316: the model spans from time t = 0 (min event_num) # to time t = 315 (max report_num). infer_max_time(mpox_tn)
Count observation model for the (truncation-corrected) case counts.
poisson_likelihood(mu = numeric(0)) nb_likelihood(mu = numeric(0), phi = lognormal_prior(log(20), 0.5))poisson_likelihood(mu = numeric(0)) nb_likelihood(mu = numeric(0), phi = lognormal_prior(log(20), 0.5))
mu |
Log-scale mean intercept prior (or fixed numeric). |
phi |
Negative-binomial overdispersion prior (or fixed numeric);
NB only. Defaults to |
A likelihood_class object.
poisson_likelihood() nb_likelihood() # Wider overdispersion (heavier-tailed counts) -- set via the likelihood: nb_likelihood(phi = lognormal_prior(log(5), 0.5))poisson_likelihood() nb_likelihood() # Wider overdispersion (heavier-tailed counts) -- set via the likelihood: nb_likelihood(phi = lognormal_prior(log(5), 0.5))
save_nowcast()
Restores the common tbl_nowcast result from a bundle written by
save_nowcast(). The native fit is restored into @fit. The
result works with predict(), autoplot(), coef(), parameters(),
mean()/median()/quantile() straight away (sampling from the stored
Laplace mode + precision). To re-fit it – on the same or new data – pass the
loaded object's model to nowcast() (the saved tbl_now is in the data
slot).
load_nowcast(file, rebuild = FALSE)load_nowcast(file, rebuild = FALSE)
file |
Path to a |
rebuild |
If |
A diseasenowcasting subclass of tbl.now::tbl_nowcast.
Combines a likelihood, an epidemic process, a reporting-delay distribution, and an optional revision process into a model object. Arguments are positional: likelihood, epidemic process, reporting delay, then revision process. Any argument can be omitted by naming the later components.
model( likelihood = nb_likelihood(), epidemic = hsgp_epidemic(), delay = lognormal_delay(), revision = no_revision(), covariate_prior = std_normal_prior(), strata_pooling = "independent", cumulative = no_cumulative() )model( likelihood = nb_likelihood(), epidemic = hsgp_epidemic(), delay = lognormal_delay(), revision = no_revision(), covariate_prior = std_normal_prior(), strata_pooling = "independent", cumulative = no_cumulative() )
likelihood |
A |
epidemic |
An |
delay |
A |
revision |
A |
covariate_prior |
A |
strata_pooling |
|
cumulative |
Dedicated count-cumulative observation configuration
from |
A model_class object.
model() model(poisson_likelihood(), hsgp_epidemic(gp_kernel = "matern52")) model(nb_likelihood(), ar1_epidemic(), lognormal_delay()) model(nb_likelihood(), ar1_epidemic(), lognormal_delay(), revision_process()) model(nb_likelihood(), hsgp_epidemic(), lognormal_delay(), strata_pooling = "hierarchical")model() model(poisson_likelihood(), hsgp_epidemic(gp_kernel = "matern52")) model(nb_likelihood(), ar1_epidemic(), lognormal_delay()) model(nb_likelihood(), ar1_epidemic(), lognormal_delay(), revision_process()) model(nb_likelihood(), hsgp_epidemic(), lognormal_delay(), strata_pooling = "hierarchical")
The main entry point. Takes a tbl_now (from the tbl.now package) and a
model(), fits the latent-epidemic + reporting-delay model as of a given
date, and returns the common tbl.now::tbl_nowcast result. Predictive draws
and quantiles are materialised in that result; the untouched native fit is
retained in @fit for predict(), mean()/median()/quantile(),
coef(), and model-specific diagnostics.
nowcast( data, model = diseasenowcasting::model(), type = c("two_stage", "one_stage", "auto"), now = NULL, K = 25L, n_draws = 2000L, delay_window = 120L, np_spread = 1, floor_mu = 0.08, floor_sig_frac = 0.08, temporal_effects = "auto", prior_only = FALSE, quantile_levels = tbl.now::nowcast_quantile_levels(), seed = sample.int(.Machine$integer.max, 1), ... )nowcast( data, model = diseasenowcasting::model(), type = c("two_stage", "one_stage", "auto"), now = NULL, K = 25L, n_draws = 2000L, delay_window = 120L, np_spread = 1, floor_mu = 0.08, floor_sig_frac = 0.08, temporal_effects = "auto", prior_only = FALSE, quantile_levels = tbl.now::nowcast_quantile_levels(), seed = sample.int(.Machine$integer.max, 1), ... )
data |
A |
model |
A |
type |
|
now |
As-of date; only events/reports up to |
K |
Number of delay imputations for the two-stage path. |
n_draws |
Default number of posterior draws used by |
delay_window |
Recent window length for the parametric Stage-1 delay fit. |
np_spread |
Dirichlet simplex imputation covariance inflation (default 1). |
floor_mu, floor_sig_frac
|
Imputation-spread floors (parametric families). |
temporal_effects |
Controls automatic seasonal / day-of-week covariates.
|
prior_only |
If |
quantile_levels |
Probabilities at which to summarise the predictive draws in the returned tbl.now::tbl_nowcast. |
seed |
Optional RNG seed (imputation draws). |
... |
Passed to |
A diseasenowcasting subclass of tbl.now::tbl_nowcast. The native
fitted model is retained in @fit; diseasenowcasting operations unwrap it
automatically. The @fit_diagnostics property (also available at
@metadata$diseasenowcasting$fit_diagnostics) records the resolved fitting
stage, imputation retention, optimizer adequacy, curvature, and any
Laplace-precision regularization used for prediction.
A revision process – reports that are later confirmed or retracted –
is a row-level mechanism for linelist and count-incidence data. It is detected
from the data, not requested by an argument. nowcast() attaches one when the
tbl_now carries revision_date / revision_type (see
tbl.now::add_revision_date()). The mode
(confirmation_only / retraction_only / both) is read from
unique(revision_type) over the full data, so it is stable across as-of
dates. Configure p and the lag with
model(revision = revision_process(...)), which always wins over the
detected default; assert the mode with revision_process(mode = ).
Revision-date censoring is likewise data metadata: set
is_censored_revision when constructing the tbl_now. There is no
nowcast() column-name argument for it.
Count-cumulative revisions instead use cumulative_process() and do
not estimate a separate revision probability p.
With type = "one_stage", the epidemic process, event-to-report delay,
report-to-revision delay, and revision probability p are estimated in
one joint objective. With type = "two_stage", Stage 1 estimates the
event-to-report delay and draws K imputations from its Laplace
approximation. Each Stage-2 fit conditions on one imputed reporting-delay
distribution and jointly estimates the epidemic process, revision delay,
and p. Posterior draws within a Stage-2 fit propagate revision uncertainty;
pooling across the fits additionally propagates reporting-delay uncertainty.
This is the same stepwise boundary used by the original event-to-report model:
two_stage separates the reporting process from the downstream model, while
the revision block remains downstream of a report.
phi)The negative-binomial overdispersion prior is not an argument of
nowcast(). Set it on the likelihood instead, e.g.
model(nb_likelihood(phi = lognormal_prior(log(5), 0.5)), ...). The default
nb_likelihood() already uses lognormal_prior(log(20), 0.5).
diseasenowcasting_workflows for when to use native modelling
operations versus the shared tbl.now result workflow; backtest() and
auto_nowcast() for retrospective comparison and automatic selection.
if (requireNamespace("tbl.now", quietly = TRUE)) { # data <- tbl.now::tbl_now(my_linelist, event_date = onset, report_date = reported) # nc <- nowcast(data, model(nb_likelihood(), hsgp_epidemic(), lognormal_delay())) # predict(nc); median(nc); coef(nc) }if (requireNamespace("tbl.now", quietly = TRUE)) { # data <- tbl.now::tbl_now(my_linelist, event_date = onset, report_date = reported) # nc <- nowcast(data, model(nb_likelihood(), hsgp_epidemic(), lognormal_delay())) # predict(nc); median(nc); coef(nc) }
Produces three stacked panels: (1) reporting-delay histogram with the
fitted density, (2) smoothed latent incidence lambda with 50/90% CI, and (3)
posterior-predictive nowcast with observed counts. Requires the
patchwork package for a combined output.
nowcast_diagnostic( object, n_draws = NULL, seed = sample.int(.Machine$integer.max, 1), previous_times = 30 )nowcast_diagnostic( object, n_draws = NULL, seed = sample.int(.Machine$integer.max, 1), previous_times = 30 )
object |
A |
n_draws |
Number of posterior draws (default: min( |
seed |
Optional RNG seed. |
previous_times |
Number of most recent event-times to display in panels 2
(smoothed incidence) and 3 (nowcast): only the last |
A patchwork combined plot (or a named list of three ggplots).
Two-stage multiple-imputation nowcast
nowcast_twostage( model, m, X = NULL, d_star = NULL, max_time = NULL, target = NULL, delay_window = 120L, K = 25L, floor_mu = 0.08, floor_sig_frac = 0.08, np_spread = 1, n_draws_per = 200L, phi = lognormal_prior(log(20), 0.5), probs = c(0.025, 0.05, 0.1, 0.25, 0.5, 0.75, 0.9, 0.95, 0.975), seed = sample.int(.Machine$integer.max, 1) )nowcast_twostage( model, m, X = NULL, d_star = NULL, max_time = NULL, target = NULL, delay_window = 120L, K = 25L, floor_mu = 0.08, floor_sig_frac = 0.08, np_spread = 1, n_draws_per = 200L, phi = lognormal_prior(log(20), 0.5), probs = c(0.025, 0.05, 0.1, 0.25, 0.5, 0.75, 0.9, 0.95, 0.975), seed = sample.int(.Machine$integer.max, 1) )
model |
A |
m |
Observation matrix |
X |
Optional covariate matrix ( |
d_star |
Optional max-observable-delay vector. |
max_time |
Time-window length; defaults to |
target |
Event-time to nowcast (default newest). |
delay_window |
Recent window length for the Stage-1 delay fit. |
K |
Number of delay imputations. |
floor_mu |
Floor on the log-mean imputation SD (parametric families). |
floor_sig_frac |
Floor on the delay-SD imputation SD (fraction of sigma). |
np_spread |
Dirichlet only: covariance-inflation factor for the simplex
imputation (samples |
n_draws_per |
Posterior nowcast draws per imputation. |
phi |
NB overdispersion prior (default |
probs |
Quantile probabilities to report. |
seed |
Optional RNG seed. |
A list with quantiles, median, pooled draws, the rung used
("multi", "anchored", or "onestage"), n_samp (imputations pooled),
and fit_diagnostics. The latter records the requested and retained
imputation counts plus warm, Stage-1, fallback, exclusion, and retained-fit
diagnostics.
Returns all estimated parameters as a long data frame with credible intervals derived from the Laplace approximation posterior precision matrix.
parameters(x, conf.level = 0.95, ...)parameters(x, conf.level = 0.95, ...)
x |
A |
conf.level |
Credible level for the interval (default 0.95). |
... |
Unused. |
tidy() on a nowcast gives you the nowcast – the predicted counts, via
tbl.now's method for the broom generic. This function gives you the
parameters. They are different questions, and this package used to answer
the second one under the first one's name.
A data.frame with columns term, estimate, std.error,
conf.low, conf.high, type.
if (requireNamespace("tbl.now", quietly = TRUE)) { # nc <- nowcast(data, model()) # parameters(nc) }if (requireNamespace("tbl.now", quietly = TRUE)) { # nc <- nowcast(data, model()) # parameters(nc) }
Prepare data for the RTMB nowcast engine
prepare_data( model, m, m_censored = NULL, X = NULL, d_star = NULL, delay_only = FALSE, max_time = NULL, num_strata = NULL, gp_L = 1.5, gp_boundary_frac = 0.62, ar_sigma_max = 1, is_confirmation = FALSE, cumulative_levels = NULL, cumulative_previous_nonzero = NULL, cumulative_settlement = NULL, retraction = NULL, resolution_mode = 0L, ... )prepare_data( model, m, m_censored = NULL, X = NULL, d_star = NULL, delay_only = FALSE, max_time = NULL, num_strata = NULL, gp_L = 1.5, gp_boundary_frac = 0.62, ar_sigma_max = 1, is_confirmation = FALSE, cumulative_levels = NULL, cumulative_previous_nonzero = NULL, cumulative_settlement = NULL, retraction = NULL, resolution_mode = 0L, ... )
model |
A |
m |
Observation matrix: columns |
m_censored |
Optional censored-observation matrix (same layout). |
X |
Optional covariate matrix ( |
d_star |
Optional max-observable-delay vector; if NULL, computed as
|
delay_only |
If TRUE, only the delay process is prepared/fit. |
max_time |
Time-window length; defaults to |
num_strata |
Number of stratum cells (the K-way product of the strata
levels). If |
gp_L |
HSGP boundary factor (> 1). Default 1.5. |
gp_boundary_frac |
Fraction of the HSGP domain placed left of the data. Default 0.62. |
ar_sigma_max |
Upper bound on the AR/beta RW innovation SD. Default 1. |
is_confirmation |
Deprecated legacy switch. |
cumulative_levels |
Optional cumulative levels aligned row-for-row with
|
cumulative_previous_nonzero |
Optional indicators, aligned with |
cumulative_settlement |
Optional positive integer settlement horizon
|
retraction |
Optional list of linelist-retraction sufficient statistics
from |
resolution_mode |
|
... |
Reserved. |
A named list of engine inputs.
Specify a prior distribution for any estimated parameter. The result is a
prior_class object that can be passed to any parameter slot in
poisson_likelihood(), nb_likelihood(), the delay process constructors,
or the epidemic process constructors.
std_normal_prior() normal_prior(mean, sd) cauchy_prior(location, scale) student_t_prior(df, location, scale) double_exponential_prior(mean, sd) flat_prior() positive_flat_prior() half_std_normal_prior() half_normal_prior(location, scale) half_cauchy_prior(location, scale) half_student_t_prior(df, scale) half_double_exponential_prior(mean, sd) gamma_prior(shape, rate) weibull_prior(shape, scale) inv_gamma_prior(shape, scale) lognormal_prior(meanlog, sdlog) chi_square_prior(df) exponential_prior(rate) logistic_prior(location, scale) beta_prior(alpha, beta)std_normal_prior() normal_prior(mean, sd) cauchy_prior(location, scale) student_t_prior(df, location, scale) double_exponential_prior(mean, sd) flat_prior() positive_flat_prior() half_std_normal_prior() half_normal_prior(location, scale) half_cauchy_prior(location, scale) half_student_t_prior(df, scale) half_double_exponential_prior(mean, sd) gamma_prior(shape, rate) weibull_prior(shape, scale) inv_gamma_prior(shape, scale) lognormal_prior(meanlog, sdlog) chi_square_prior(df) exponential_prior(rate) logistic_prior(location, scale) beta_prior(alpha, beta)
mean, location
|
Location parameter. |
sd, scale
|
Scale parameter (must be > 0). |
df |
Degrees of freedom (must be > 0). |
shape |
Shape parameter (must be > 0). |
rate |
Rate parameter (must be > 0). |
meanlog |
Mean on the log scale (lognormal). |
sdlog |
SD on the log scale (lognormal, must be > 0). |
alpha, beta
|
Beta distribution shape parameters (must be > 0). |
A prior_class object.
normal_prior(0, 1) gamma_prior(2, 0.1) exponential_prior(1) flat_prior()normal_prior(0, 1) gamma_prior(2, 0.1) exponential_prior(1) flat_prior()
The distribution g_C of the revision lag: how long after a case is
reported its result comes back. Pass one to revision_process() as
revision_delay.
lognormal_revision(mu = numeric(0), sigma = numeric(0)) gamma_revision(shape = numeric(0), rate = numeric(0)) generalized_gamma_revision(mu = numeric(0), sigma = numeric(0), Q = numeric(0)) dirichlet_revision(alpha = numeric(0), bins = numeric(0))lognormal_revision(mu = numeric(0), sigma = numeric(0)) gamma_revision(shape = numeric(0), rate = numeric(0)) generalized_gamma_revision(mu = numeric(0), sigma = numeric(0), Q = numeric(0)) dirichlet_revision(alpha = numeric(0), bins = numeric(0))
mu |
Log-mean intercept ( |
sigma |
Log-scale / SD parameter > 0. |
shape, rate
|
Gamma delay parameters (the |
Q |
GenGamma shape ( |
alpha |
Dirichlet concentration (scalar broadcast to all bins). |
bins |
Dirichlet: number of explicit delay bins (geometric tail beyond). |
These are aliases of the corresponding delay_process constructors – a
revision lag is an ordinary non-negative delay, only measured from the
report rather than from the event – so the parameters, priors and
behaviour are identical. The one
difference is the support, and it depends on the mode: under
retraction_only the lag lives on {1, 2, ...} (a retraction lands strictly
after the report it withdraws, so a case retracted in the same period is
dropped by nowcast() – it was never visible in any data vintage), while
under confirmation_only and both it lives on {0, 1, ...}, since a test
can come back the day it was ordered.
A delay_process_class object, for the revision_delay slot of
revision_process().
dirichlet_revision() is the safest default when the counts are large. The
correction applied to a pending report of age j is
rho(j) = p / (p + (1 - p) * (1 - G_C(j))), so at high counts a shape error
in g_C biases the nowcast by more than its Monte-Carlo noise: on a COVID
series of ~8000 cases/day a lognormal g_C fitted to a 1 + Poisson(2) lag
left a ~0.9% bias and lost nominal coverage, while the Dirichlet recovered
rho to four decimals. At low counts the parametric families are fine and
estimate fewer parameters.
revision_process(), delay_process, nowcast()
revision_process(lognormal_revision()) revision_process(dirichlet_revision(bins = 10)) # Held-fixed revision lag, e.g. from an external study revision_process(gamma_revision(shape = log(3), rate = 2))revision_process(lognormal_revision()) revision_process(dirichlet_revision(bins = 10)) # Held-fixed revision lag, e.g. from an external study revision_process(gamma_revision(shape = log(3), rate = 2))
A report is rarely a case outright. It is provisional, and resolved exactly once: a test comes back, and the report is either confirmed (a real case) or retracted (removed from the register). A revision process models that second step, so the nowcast targets the settled count rather than the raw report count.
revision_process( revision_delay = lognormal_delay(), p = numeric(0), stratified_p = FALSE, mode = c("auto", "confirmation_only", "retraction_only", "both") )revision_process( revision_delay = lognormal_delay(), p = numeric(0), stratified_p = FALSE, mode = c("auto", "confirmation_only", "retraction_only", "both") )
revision_delay |
A |
p |
Probability that a report resolves positive – that it is a real
case and is never retracted. Either a
|
stratified_p |
If |
mode |
Which outcomes the data record. |
Attach it with model(revision = revision_process(...)). nowcast()
switches it on automatically when the data carry it – see Detection below.
A revision_process_class object, for model(revision = ).
What differs between surveillance systems is only which resolutions get recorded, and all three possibilities share one likelihood:
retraction_only |
confirmation_only |
both |
|
| outcomes recorded | the negatives | the positives | both signs |
| a missing outcome means | not retracted yet | not confirmed yet | not resolved yet |
| target | reports never retracted | reports eventually confirmed | reports resolving positive |
| lag support | {1, 2, ...} |
{0, 1, ...} |
{0, 1, ...}
|
The lag support differs only because a retraction in the same period as its report describes a case never visible in any data vintage, whereas a test coming back the day it was ordered is ordinary.
In mode = "both", revision_delay is one shared law for positive and
negative resolutions. This first prototype deliberately does not fit separate
competing-risk lag laws. In a one-outcome mode, the same argument denotes the
lag for the outcome that is recorded: confirmation or retraction respectively.
"auto" reads unique(revision_type) over the full data, not the as-of
view, so the mode is a stable property of the data source and does not flip
between backtest dates.
nowcast() attaches a revision process when the tbl_now carries
revision_date / revision_type (see
tbl.now::add_revision_date()).
Count-cumulative data instead use cumulative_process(), whose primitive
retraction object is the finite-age kernel h_R; it does not estimate p.
Do not use this component for an aggregate cumulative stream. Configure its
down-revisions with cumulative_process(). Without individual report
outcomes, p and a conditional revision-delay law are not separately
identified.
revision_delay inherits the default priors of its delay family (see
delay_process). p is described under its argument above. At p = 1 the
revision layer is inert and the model is the ordinary count model.
model(), revision_delay, delay_process, nowcast()
# Results come back on one shared timescale, whatever the answer: revision_process(lognormal_revision()) # Attach to a model: model(nb_likelihood(), hsgp_epidemic(), lognormal_delay(), revision = revision_process( revision_delay = dirichlet_revision(bins = 10), p = beta_prior(20, 3), stratified_p = TRUE))# Results come back on one shared timescale, whatever the answer: revision_process(lognormal_revision()) # Attach to a model: model(nb_likelihood(), hsgp_epidemic(), lognormal_delay(), revision = revision_process( revision_delay = dirichlet_revision(bins = 10), p = beta_prior(20, 3), stratified_p = TRUE))
base::sample())A generic that draws from a prior_class with the appropriate random
number generator. For any other object (numeric, character, factor,
Date, list, ...) it dispatches to base::sample(), so sample() keeps its
usual base behaviour outside the package.
sample(object, size, ...) ## S7 method for class <any> sample(object, size, ...) ## S7 method for class <diseasenowcasting::prior_class> sample(object, size, ...)sample(object, size, ...) ## S7 method for class <any> sample(object, size, ...) ## S7 method for class <diseasenowcasting::prior_class> sample(object, size, ...)
object |
A |
size |
Number of draws. |
... |
Passed through to the underlying sampler. |
For a prior, a numeric vector of length size; otherwise whatever
base::sample() returns.
sample(normal_prior(log(7), 0.5), 10) sample(gamma_prior(2, 0.1), 10) sample(as.Date("2020-01-01") + 0:9, 3) # falls back to base::sample()sample(normal_prior(log(7), 0.5), 10) sample(gamma_prior(2, 0.1), 10) sample(as.Date("2020-01-01") + 0:9, 3) # falls back to base::sample()
Serializes a nowcast() (or auto_nowcast()) result to a single .rds
bundle that can be restored later with load_nowcast(). The RTMB autodiff
tape is not saved (its external pointers cannot be serialized); instead
the fitted parameters and the Laplace mode + precision are stored, which is
all predict() needs. The input tbl_now and the model() specification
are saved too, so the loaded object can also be re-fit.
save_nowcast(object, file)save_nowcast(object, file)
object |
A result from |
file |
Path to write (a single |
file, invisibly.
if (requireNamespace("tbl.now", quietly = TRUE)) { library(tbl.now) data(denguedat) dn <- subset(denguedat, onset_week <= as.Date("1990-12-01") & report_week <= as.Date("1990-12-01")) tn <- tbl_now(dn, event_date = onset_week, report_date = report_week, data_type = "linelist", verbose = FALSE) nc <- nowcast(tn, type = "one_stage", n_draws = 200) f <- tempfile(fileext = ".rds") save_nowcast(nc, f) nc2 <- load_nowcast(f) predict(nc2, summary = TRUE) # works without the original RTMB tape }if (requireNamespace("tbl.now", quietly = TRUE)) { library(tbl.now) data(denguedat) dn <- subset(denguedat, onset_week <= as.Date("1990-12-01") & report_week <= as.Date("1990-12-01")) tn <- tbl_now(dn, event_date = onset_week, report_date = report_week, data_type = "linelist", verbose = FALSE) nc <- nowcast(tn, type = "one_stage", n_draws = 200) f <- tempfile(fileext = ".rds") save_nowcast(nc, f) nc2 <- load_nowcast(f) predict(nc2, summary = TRUE) # works without the original RTMB tape }
auto_nowcast() used to pick the winnerThe metric auto_nowcast() used to pick the winner
selection_metric(nc)selection_metric(nc)
nc |
A |
The scoringutils score-column name supplied to auto_nowcast().
Consult the result's
relative_score comparison field to determine whether its relative skill
was used.
auto_nowcast(), comparison_scores()
auto_nowcast()
Fitting times recorded by auto_nowcast()
selection_timings(nc)selection_timings(nc)
nc |
A result returned by |
A list with backtest (one row per attempted retrospective fit),
refit (the full-data refit attempts), and total_seconds for the complete
automatic-selection call.
auto_nowcast(), comparison_scores()
Matches diseasenowcast2's .summarise_nowcast_matrix() exactly: one row per
event (.event_num 0-indexed, ascending) with columns
mean, median, sd, mad, q2.5, q5, q10, q25, q50, q75, q90, q95, q97.5.
summarise_nowcast_matrix(draws_matrix)summarise_nowcast_matrix(draws_matrix)
draws_matrix |
Pooled draws matrix |
A data.frame in the scoring-pipeline format.
Given a fitted nowcast and new data that arrived after the fit, compute how surprising the new data is under the fitted posterior.
Two types of surprise are supported (select via the type argument):
"count": How surprising is a new case count at a given event time
"delay": How surprising is a report with a very long delay
"both": Compute both (requires appropriate new_data format).
surprise( object, new_data, type = c("both", "count", "delay"), level = 0.99, n_draws = 500L, seed = sample.int(.Machine$integer.max, 1) )surprise( object, new_data, type = c("both", "count", "delay"), level = 0.99, n_draws = 500L, seed = sample.int(.Machine$integer.max, 1) )
object |
A |
new_data |
For |
type |
Which surprise type(s) to compute. One of |
level |
Credible level used to flag surprises (default |
n_draws |
Number of posterior draws (default 500). |
seed |
Optional RNG seed. |
A list with:
$count_surprise (data.frame): one row per event time with
event_index, observed, posterior_mean, lpd, ppp_right,
ppp_left, relative_surprise.
$delay_surprise (data.frame): one row per delay value with
delay, tail_prob (= 1 - G_D(delay)), lpd, relative_surprise.
Surprise score on a raw fit() result
## S3 method for class 'list' surprise( object, new_data, type = c("both", "count", "delay"), level = 0.99, n_draws = 500L, seed = sample.int(.Machine$integer.max, 1) )## S3 method for class 'list' surprise( object, new_data, type = c("both", "count", "delay"), level = 0.99, n_draws = 500L, seed = sample.int(.Machine$integer.max, 1) )
object |
A list returned by |
new_data |
See |
type |
See |
level |
See |
n_draws |
See |
seed |
See |
ggplot2 theme matching the diseasenowcasting visual identity
theme_diseasenowcasting(base_size = 11)theme_diseasenowcasting(base_size = 11)
base_size |
Base font size. |
validate_custom_delay(delay, test_theta = NULL, test_delays = 1:14)validate_custom_delay(delay, test_theta = NULL, test_delays = 1:14)
delay |
A |
test_theta |
Optional numeric vector of length |
test_delays |
Integer vector of delay values to evaluate (default 1:14). |
Tests the user-supplied cdf / log_cdf / log_survival on a dummy
numeric vector, checking that obj$fn() and obj$gr() are finite. Emits a
clear error (including common causes) on failure.
Invisibly returns TRUE on success.
# Custom components tape USER functions, so RTMB must be attached # (it is kept in Imports, not Depends, so attach it yourself): library(RTMB) weibull_cdf <- function(theta) { shape <- exp(theta[1]); scale <- exp(theta[2]) function(d) 1 - exp(-(d / scale)^shape) } dly <- custom_delay(weibull_cdf, priors = list(normal_prior(0, 1), normal_prior(log(7), 1)), name = "Weibull", inits = c(0, log(7))) validate_custom_delay(dly)# Custom components tape USER functions, so RTMB must be attached # (it is kept in Imports, not Depends, so attach it yourself): library(RTMB) weibull_cdf <- function(theta) { shape <- exp(theta[1]); scale <- exp(theta[2]) function(d) 1 - exp(-(d / scale)^shape) } dly <- custom_delay(weibull_cdf, priors = list(normal_prior(0, 1), normal_prior(log(7), 1)), name = "Weibull", inits = c(0, log(7))) validate_custom_delay(dly)
validate_custom_epidemic(epidemic, test_theta = NULL)validate_custom_epidemic(epidemic, test_theta = NULL)
epidemic |
A |
test_theta |
Optional numeric vector of length |
Passes intensity_fn through RTMB::MakeADFun at the supplied (or default)
initial values and confirms that the objective value and gradient are both
finite. Emits an informative error if the function is not AD-safe.
epidemic, invisibly. Emits a success message if the check passes.
# Custom components tape USER functions, so RTMB must be attached # (it is kept in Imports, not Depends, so attach it yourself): library(RTMB) max_time <- 15L rw_fn <- function(theta) matrix(theta[1] + cumsum(exp(theta[2]) * theta[3:(2L + max_time)]), max_time, 1L) custom_epi <- custom_epidemic( rw_fn, priors = c(list(normal_prior(2, 1), normal_prior(-2, 0.5)), rep(list(std_normal_prior()), max_time)), inits = c(2, -2, rep(0, max_time)) ) validate_custom_epidemic(custom_epi)# Custom components tape USER functions, so RTMB must be attached # (it is kept in Imports, not Depends, so attach it yourself): library(RTMB) max_time <- 15L rw_fn <- function(theta) matrix(theta[1] + cumsum(exp(theta[2]) * theta[3:(2L + max_time)]), max_time, 1L) custom_epi <- custom_epidemic( rw_fn, priors = c(list(normal_prior(2, 1), normal_prior(-2, 0.5)), rep(list(std_normal_prior()), max_time)), inits = c(2, -2, rep(0, max_time)) ) validate_custom_epidemic(custom_epi)