Introduction to diseasenowcasting: Real-Time Epidemic Nowcasting
diseasenowcasting team
Source:vignettes/introduction.Rmd
introduction.Rmddiseasenowcasting is an R package for nowcasting time
series of epidemiological cases. Epidemiologic surveillance tools
usually have an intrinsic delay between the true date of an
event (event_date) and the report date for
that event (report_date). Some examples include
the true date being symptom onset or testing time and the report date
corresponds to when the case was registered in the system.
diseasenowcasting uses censored Bayesian models (via R’s
Template Model Builder RTMB)
to infer the cases that have not yet been reported thus providing a
prediction of the final number of cases.
Your data: the tbl_now format
diseasenowcasting works with data organised as a
tbl_now object from the companion tbl.now
package. A tbl_now is simply a data frame that has been
annotated with the roles of its columns:
event date when the event happened (e.g. symptom onset)
report date when the event was reported (e.g. date entered into the database)
strata (optional) columns that defined all the strata (e.g. sex and region)
now (optional) the date until which to nowcast (assumes all events and reports before the now have been observed and missing observations correspond to no observations - i.e. if one day there were not cases the missingness can be translated into zero cases)
As a quick example, here is how to build a tbl_now using
the following surveillance data for dengue in Puerto Rico:
data(denguedat)#> onset_week report_week gender
#> 1 1990-01-01 1990-01-01 Male
#> 2 1990-01-01 1990-01-01 Female
#> 3 1990-01-01 1990-01-01 Female
#> 4 1990-01-01 1990-01-08 Female
#> 5 1990-01-01 1990-01-08 Male
#> 6 1990-01-01 1990-01-15 FemaleWe can transform the data.frame to a
tbl_now by specifying the event and report dates
(onset and report weeks respectively) as well
as the data_type and the strata (in this case,
gender).
dengue_tbl <- tbl_now(
denguedat,
event_date = onset_week, # symptom onset date
report_date = report_week, # when the record was reported
data_type = "linelist", # another option is "count-incidence" if data is aggregated
now = as.Date("1991-01-01") #When is the now of the nowcast
)
dengue_tbl
#> # A tibble: 52,987 × 6
#> # Data type: "linelist"
#> # Frequency: Event: `weeks` | Report: `weeks`
#> onset_week report_week gender .event_num .report_num .delay
#> <date> <date> <chr> <dbl> <dbl> <dbl>
#> [event_date] [report_date] [...] [...] [...] [...]
#> 1 1990-01-01 1990-01-01 Male 0 0 0
#> 2 1990-01-01 1990-01-01 Female 0 0 0
#> 3 1990-01-01 1990-01-01 Female 0 0 0
#> 4 1990-01-01 1990-01-08 Female 0 1 1
#> 5 1990-01-01 1990-01-08 Male 0 1 1
#> 6 1990-01-01 1990-01-15 Female 0 2 2
#> 7 1990-01-01 1990-01-15 Female 0 2 2
#> 8 1990-01-01 1990-01-15 Female 0 2 2
#> 9 1990-01-01 1990-01-22 Female 0 3 3
#> 10 1990-01-01 1990-01-08 Female 0 1 1
#> # ────────────────────────────────────────────────────────────────────────────────
#> # Now: 1991-01-01 | Event date: "onset_week" | Report date: "report_week"
#> # ────────────────────────────────────────────────────────────────────────────────
#> # ℹ 52,977 more rowsOnce your data is a tbl_now, a single call to
nowcast() does the rest.
For more information about
tbl_nowcheck the package’s website.
Example 1 – Dengue fever (setting up a stratified nowcast)
We fit a nowcast stratified by gender to illustrate the basic
workflow. First we add the column gender as strata to the
tbl_now:
dengue_tbl <- dengue_tbl |> add_strata(gender)Notice that the tbl_now automatically prints the
strata specification below:
dengue_tbl
#> # A tibble: 52,987 × 6
#> # Data type: "linelist"
#> # Frequency: Event: `weeks` | Report: `weeks`
#> onset_week report_week gender .event_num .report_num .delay
#> <date> <date> <chr> <dbl> <dbl> <dbl>
#> [event_date] [report_date] [strata] [...] [...] [...]
#> 1 1990-01-01 1990-01-01 Male 0 0 0
#> 2 1990-01-01 1990-01-01 Female 0 0 0
#> 3 1990-01-01 1990-01-01 Female 0 0 0
#> 4 1990-01-01 1990-01-08 Female 0 1 1
#> 5 1990-01-01 1990-01-08 Male 0 1 1
#> 6 1990-01-01 1990-01-15 Female 0 2 2
#> 7 1990-01-01 1990-01-15 Female 0 2 2
#> 8 1990-01-01 1990-01-15 Female 0 2 2
#> 9 1990-01-01 1990-01-22 Female 0 3 3
#> 10 1990-01-01 1990-01-08 Female 0 1 1
#> # ────────────────────────────────────────────────────────────────────────────────
#> # Now: 1991-01-01 | Event date: "onset_week" | Report date: "report_week"
#> # Strata: "gender"
#> # ────────────────────────────────────────────────────────────────────────────────
#> # ℹ 52,977 more rowsWe can also add temporal effects for example a weekly seasonality (52
seasons) as well as a holiday effect using the almanac
package:
library(almanac)
#Specify 52 seasons (weekly) and holidays from the US
t_effects <- temporal_effects(seasons = 52, holidays = cal_us_federal())
#Add the temporal effects
dengue_tbl <- dengue_tbl |>
add_temporal_effects(t_effects)Finally we fit the model:
nc_dengue <- nowcast(dengue_tbl)The fitted model can be visualized with autoplot(). Note
that the nowcast was already stratified by the strata specified in the
tbl_now:
autoplot(nc_dengue) 
Nowcast for dengue example. The shaded bars show the median while the errorbar has the 90% credible intervals.
Values can be obtained via predict() and
summary():
# Full posterior-predictive nowcast at every event-time
pred_dengue <- predict(nc_dengue)
#This creates a summary of mean and quantiles
summary(pred_dengue) #> mean median sd mad q2.5 q5 q10 q25 q50 q75 q90 q95
#> 154 108.6480 108 1.538597 1.4826 107 107 107 107.75 108 109 111 111.00
#> 155 89.1330 89 2.327664 1.4826 86 86 87 87.00 89 90 92 93.00
#> 156 68.1615 67 3.834680 2.9652 63 63 64 65.00 67 70 73 75.00
#> 157 45.2905 44 7.167901 5.9304 36 37 38 40.00 44 49 54 59.00
#> 158 40.6945 38 13.773635 10.3782 24 25 27 32.00 38 46 56 65.05
#> 159 36.4080 33 18.350747 16.3086 12 14 17 23.00 33 46 61 71.00
#> q97.5 .event_num stratum event_date
#> 154 112.000 47 Total 1990-11-26
#> 155 95.000 48 Total 1990-12-03
#> 156 78.000 49 Total 1990-12-10
#> 157 62.025 50 Total 1990-12-17
#> 158 73.025 51 Total 1990-12-24
#> 159 81.000 52 Total 1990-12-31Additionally the nowcast_diagnostic() shows the fitted
distribution for the delay, the smoothed epidemic process as well as the
aggregated nowcast (for the sum of all strata):
nowcast_diagnostic(nc_dengue) 
Example 2 – Mpox (modifying the nowcast model)
The mpoxdat dataset (also in tbl.now)
covers the 2022 mpox outbreak in New York City with daily case counts
stratified by race.
data(mpoxdat)
mpox_tbl <- tbl_now(
mpoxdat,
event_date = dx_date,
report_date = dx_report_date,
case_count = n,
data_type = "count-incidence",
now = as.Date("2022-08-15")
) A simple plot of the data shows that we should be taking into account day-of-the-week effects:
autoplot(mpox_tbl)
We can also set it again with
add_temporal_effects():
mpox_tbl <- mpox_tbl |>
add_temporal_effects(temporal_effects(day_of_week = TRUE))You can see that the tbl_now indicates its
computation:
#> # A tibble: 1,417 × 7
#> # Data type: "count-incidence"
#> # Frequency: Event: `days` | Report: `days`
#> dx_date dx_report_date race n .event_num .report_num .delay
#> <date> <date> <chr> <int> <dbl> <dbl> <dbl>
#> [event_date] [report_date] [...] [cas… [...] [...] [...]
#> 1 2022-07-08 2022-07-12 Asian 4 0 4 4
#> 2 2022-07-08 2022-07-12 Black 6 0 4 4
#> 3 2022-07-08 2022-07-12 Hispanic 6 0 4 4
#> 4 2022-07-08 2022-07-12 Non-Hispanic… 6 0 4 4
#> 5 2022-07-08 2022-07-13 Asian 2 0 5 5
#> 6 2022-07-08 2022-07-13 Black 3 0 5 5
#> 7 2022-07-08 2022-07-13 Hispanic 8 0 5 5
#> 8 2022-07-08 2022-07-13 Non-Hispanic… 5 0 5 5
#> 9 2022-07-08 2022-07-14 Black 1 0 6 6
#> 10 2022-07-08 2022-07-14 Hispanic 3 0 6 6
#> # ────────────────────────────────────────────────────────────────────────────────
#> # Now: 2022-08-15 | Event date: "dx_date" | Report date: "dx_report_date"
#> # T. effects (lazy): [event_date] day_of_week
#> # ────────────────────────────────────────────────────────────────────────────────
#> # ℹ 1,407 more rowsOne can also choose between several likelihoods, epidemic processess
and delay distributions and feed it into the model(). Here
we use a Susceptible-Infected-Recovered model (SIR) with a delay that
follows a lognormal distribution:
#Models can be modified via the model()
mpox_model <- model(likelihood = nb_likelihood(), #Negative binomial
epidemic = sir_epidemic(), #SIR model
delay = lognormal_delay()) #Delay distribution
#We can then fit the model
nc_mpox <- nowcast(mpox_tbl, model = mpox_model)
#And show the nowcast
autoplot(nc_mpox) 
Example 3 – Comparing models with a backtest
A backtest reruns nowcasts at multiple historical dates and scores them against the eventually-observed totals. This lets you compare between models before committing to one for real-time monitoring.
The
backtest()function fits one nowcast perdateandmodelcell. Those cells run in parallel through the future framework.
By default they run sequentially; to use several CPU cores, set a parallel plan before callingbacktest():
library(future)
plan(multisession, workers = 4) # 4 parallel R sessions
# ... run backtest() ...
plan(sequential) # restore serial execution when doneIn what follows we run a backtest in sequential mode however we strongly recommend using as many workers in a multisession plan as possible:
# Compare HSGP (flexible GP trend) vs AR1 (autoregressive trend)
# and SIR (susceptible, infected, recovered) on mpox
models_to_compare <- list(
model(nb_likelihood(), hsgp_epidemic(), lognormal_delay()),
model(nb_likelihood(), ar1_epidemic(), lognormal_delay()),
model(nb_likelihood(), sir_epidemic(), lognormal_delay())
)
#Uncomment this line to use several of your cores
#plan(multisession, workers = 4)
backtest_mpox <- backtest(
mpox_tbl,
models = models_to_compare,
n_dates = 5 #Test 3 dates at random from the mpox data
)
#This closes the plan multisession opened above
#plan(sequential) We can then plot the Weighted Interval Score (WIS) and interval coverage:
# Plot scoring metrics
autoplot(backtest_mpox)
Or obtain it via score():
# Rank by Weighted Interval Score (WIS) -- lower is better
score(backtest_mpox)
#> model wis overprediction underprediction dispersion
#> 1 SIR/nb/LogNormal 12.29247 0.000000 10.105556 2.186910
#> 2 HSGP/nb/LogNormal 12.89510 2.013889 4.527778 6.353437
#> 3 AR1/nb/LogNormal 13.28417 3.138889 2.319444 7.825833
#> coverage_50 coverage_90 ape mse n
#> 1 0.50 0.75 0.4329608 1766.000 4
#> 2 0.50 1.00 1.4150650 1679.062 4
#> 3 0.25 1.00 6.4057984 1215.812 4The score() output shows WIS, absolute percentage error
(APE), and empirical coverage at 50 % and 90 %. A well-calibrated
nowcast should have the lowest WIS and coverage close to these
levels.
In this specific test we would choose the SIR for having the lowest WIS at essentially the same coverage as HGSP. Note however that for the tutorial we only used 5 historical dates which is too low to reach a definite conslusion.
Example 4 – Letting the package choose the model
(auto_nowcast())
Doing the backtest-and-compare loop by hand (Example 3) is exactly
what auto_nowcast() automates. Given a
tbl_now, it builds a grid of candidate models sized to
how much data you have (which epidemic processes are even feasible,
crossed with the delay families), backtests them, scores them, and
refits the single best one on the full data. The result
is an ordinary nowcast, with the ranked comparison stored alongside
it.
Here we use all the dengue data up to 1994:
#All dengue observed as of January 1994
dengue_94 <- denguedat |>
filter(onset_week <= as.Date("1994-01-01") &
report_week <= as.Date("1994-01-01"))
dengue_tbl_94 <- tbl_now(
dengue_94,
event_date = onset_week,
report_date = report_week,
data_type = "linelist",
now = as.Date("1994-01-01")
)Backtesting the whole grid is the expensive step. To run the
candidates in parallel, set a future::plan() before the
call and restore it afterwards (left commented here so the vignette
stays single-process):
# Uncomment to run candidates in parallel:
# library(future)
# plan(multisession, workers = max(parallel::detectCores() - 1, 1))
auto_ncast <- auto_nowcast(
dengue_tbl_94,
metric = "wis", # rank candidates by Weighted Interval Score
n_dates = 10, # backtest at 10 historical dates (raise for a firmer choice)
n_draws_select = 150, # draws while comparing (small => fast)
n_draws = 500, # draws for the final fit of the winner
K = 8 # delay imputations (small => fast)
)
# plan(sequential)We can show the scores of the models to see the best performer:
comparison_scores(auto_ncast) # every candidate, ranked best-first
#> model wis overprediction underprediction dispersion
#> 1 HSGP/nb/GeneralizedGamma 8.502125 0.01666667 4.264444 4.221014
#> 2 HSGP/nb/LogNormal 8.540486 0.04444444 3.810000 4.686042
#> 3 HSGP/nb/Dirichlet 8.900236 0.20555556 3.876667 4.818014
#> 4 AR1/nb/Dirichlet 9.515236 0.16666667 4.177778 5.170792
#> 5 AR1/nb/GeneralizedGamma 9.730778 0.07777778 4.780000 4.873000
#> 6 AR1/nb/LogNormal 10.358556 0.08888889 5.391111 4.878556
#> 7 SIR/nb/GeneralizedGamma 21.316139 0.00000000 18.050000 3.266139
#> 8 SIR/nb/Dirichlet 23.516347 0.00000000 20.546111 2.970236
#> 9 SIR/nb/LogNormal 23.636236 0.00000000 20.493889 3.142347
#> coverage_50 coverage_90 ape mse n
#> 1 0.5 1.0 0.3484614 480.275 10
#> 2 0.3 1.0 0.3989732 582.225 10
#> 3 0.3 1.0 0.3821182 568.050 10
#> 4 0.4 0.9 0.4212942 637.150 10
#> 5 0.4 0.9 0.4103526 704.400 10
#> 6 0.4 0.9 0.4405807 785.025 10
#> 7 0.2 0.6 0.7496154 2819.350 10
#> 8 0.1 0.6 0.7324501 2765.500 10
#> 9 0.2 0.7 0.7379645 2790.175 10best_model() hands back the winning model()
object, so you can reuse the same specification on other data (or feed
it to nowcast() / backtest()):
winner <- best_model(auto_ncast)
winnerBecause the result is a normal nowcast, everything else just works:
autoplot(auto_ncast)
Nowcast from the model auto_nowcast() selected.
Updating the chosen model as new data arrive
The selected nowcast is an ordinary nowcast_class, so
once auto_nowcast() has picked a model you keep it and feed
it the next batch of reports with update() – no need to
re-run the selection. update() warm-refits the
same winning model and scores the incoming reports against the
fitted delay, warning if any arrive far later than expected:
# A few more weeks of dengue, observed as of 1994-04-01:
dengue_apr <- denguedat |>
filter(onset_week <= as.Date("1994-04-01") &
report_week <= as.Date("1994-04-01"))
dengue_tbl_apr <- tbl_now(
dengue_apr,
event_date = onset_week,
report_date = report_week,
data_type = "linelist",
now = as.Date("1994-04-01")
)
auto_ncast_updated <- update(auto_ncast, dengue_tbl_apr)
#> Warning: ! Surprising reporting delay of 11 weeks (1 report): longer than the model
#> expects (P(D >= d) = 3e-05).
#> ! Surprising reporting delay of 10 weeks (1 report): longer than the model
#> expects (P(D >= d) = 9.1e-05).
#> ! Surprising reporting delay of 9 weeks (2 reports): longer than the model
#> expects (P(D >= d) = 0.00028).
#> ! Surprising reporting delay of 8 weeks (5 reports): longer than the model
#> expects (P(D >= d) = 0.00085).
#> ! Surprising reporting delay of 7 weeks (3 reports): longer than the model
#> expects (P(D >= d) = 0.0026).
#> ! Surprising reporting delay of 6 weeks (2 reports): longer than the model
#> expects (P(D >= d) = 0.0079).
#> ℹ If these are outliers, treat them as censored with `censor_delays_above()`
#> and re-fit.
#> ℹ See all flagged delays with `extreme_values(nc)`.Any reports with surprising delays are collected by
extreme_values() (it returns NULL when nothing
looks off):
extreme_values(auto_ncast_updated)
#> delay weight mean_tail_prob cdf_prob lpd relative_surprise direction
#> 1 6 2 0.007928 0.992072 -4.6875 0.0260 long
#> 2 7 3 0.002606 0.997394 -5.7892 0.0086 long
#> 3 8 5 0.000852 0.999148 -6.9041 0.0028 long
#> 4 9 2 0.000278 0.999722 -8.0240 0.0009 long
#> 5 10 1 0.000091 0.999909 -9.1439 0.0003 long
#> 6 11 1 0.000030 0.999970 -10.2609 0.0001 long
#> surprise level
#> 1 delay 0.99
#> 2 delay 0.99
#> 3 delay 0.99
#> 4 delay 0.99
#> 5 delay 0.99
#> 6 delay 0.99See the vignette on Handling Outlier Delays with Censoring for what to do when a delay is flagged.
Example 5 – Count-cumulative data that revises up and down (FluSight)
So far every example has used incident data: each
case is counted once, and counts only ever grow as late reports arrive.
Some surveillance systems instead publish a running cumulative
total for each event-time that is re-reported week
after week – and those totals can be revised downward
as well as upward (for example when a suspected case is later
re-classified as negative). The FluSight
influenza hospitalisation data shipped with tbl.now is
exactly this kind of stream: for each target_end_date (the
epiweek being counted), the reported cumulative observation
changes across as_of report dates.
data(flusight, package = "tbl.now")
head(flusight)
#> # A tibble: 6 × 4
#> as_of target_end_date location_name observation
#> <date> <date> <chr> <dbl>
#> 1 2023-09-23 2022-02-12 Alabama 10
#> 2 2023-09-23 2022-02-12 Alaska 0
#> 3 2023-09-23 2022-02-12 Arizona 64
#> 4 2023-09-23 2022-02-12 Arkansas 29
#> 5 2023-09-23 2022-02-12 California 36
#> 6 2023-09-23 2022-02-12 Colorado 29We tell tbl_now() that these are cumulative counts with
data_type = "count-cumulative". Here we nowcast a single
location (California), keeping a recent window of the season:
california <- flusight |>
filter(location_name == "California", target_end_date >= as.Date("2023-10-01"))
flu_tbl <- tbl_now(
california,
event_date = target_end_date, # the epiweek being counted
report_date = as_of, # when that cumulative count was known
case_count = observation, # cumulative admissions (can go up OR down)
data_type = "count-cumulative",
now = as.Date("2024-01-27")
)To model the down-revisions we attach a confirmation
process to the model. confirmation_process()
describes the retraction side of the stream through a retraction delay
and a confirmation probability p – the probability
that a report is genuine and never retracted. Left unset, p
gets a strong data-informed prior (just like the lognormal delay’s
mean); pass confirmation_process(p = 0.98) to hold it
fixed, or a beta_prior() to set your own. When the data are
count-cumulative, nowcast() automatically switches from the
censored count likelihood to the signed-increment Skellam /
SkNB likelihood that the confirmation process needs:
flu_model <- model(
likelihood = nb_likelihood(),
epidemic = ar1_epidemic(),
delay = lognormal_delay(),
confirmation = confirmation_process() # models the up- and down-revisions
)
flu_ncast <- nowcast(flu_tbl, flu_model, n_draws = 1000)
autoplot(flu_ncast)
Confirmation nowcast for cumulative influenza hospitalisations in California.
Everything else works exactly as before – predict(),
summary(), coef() and the different epidemic
processes, delay families, covariates and temporal effects are all
available for count-cumulative data too:
summary(predict(flu_ncast)) |> dplyr::as_tibble() |> tail(4)
#> # A tibble: 4 × 15
#> mean median sd mad q2.5 q5 q10 q25 q50 q75 q90 q95 q97.5
#> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
#> 1 1425. 1421 9.69 2.97 1417 1418 1419 1420 1421 1426 1433 1441. 1451.
#> 2 1029. 1026 8.74 2.97 1020 1022 1023 1025 1026 1030 1040 1046. 1056
#> 3 1046. 1044 8.70 2.97 1031. 1035 1039 1043 1044 1047 1054 1061 1066.
#> 4 703. 694 25.3 4.45 685. 688 691 692 694 705 724 744 764
#> # ℹ 2 more variables: .event_num <int>, event_date <date>To nowcast several locations at once, declare the
location column as strata in tbl_now(). A
single stratified nowcast() then shares the delay and
confirmation structure across locations, which is both faster than a
separate fit per location and – on FluSight – sharper (lower Weighted
Interval Score).
Saving and loading a fitted nowcast
Fitting can take a while, so you will often want to
save a fitted nowcast and reload it later – in a
report, a dashboard, or a scheduled job – instead of re-fitting.
save_nowcast() writes it to a single .rds file
and load_nowcast() brings it back.
The autodiff engine (RTMB) cannot itself be written to
disk, so what is stored is everything needed to reuse the fit:
the model() specification, the input tbl_now,
and each fit’s parameters together with its Laplace mode and precision.
That is all predict() needs, so a reloaded nowcast behaves
exactly like the original (any number of draws, no re-fitting):
saved <- tempfile(fileext = ".rds")
save_nowcast(auto_ncast, saved)
restored <- load_nowcast(saved)
# predict() / autoplot() / coef() / tidy() all work just as before:
autoplot(restored)
Because the input data travels in the bundle, you can also re-fit the saved model later (on the original data, or on a newer extract):
nowcast(restored@data, restored@model) # re-runs the optimisationNext steps
This vignette covered the basics: building a tbl_now,
fitting a nowcast with nowcast(), inspecting results with
predict() / summary() /
autoplot(), and comparing models with
backtest() / score().
Depending on what you want to do next, check out the following vignettes:
Nowcasting at the Start of an Epidemic — A worked, end-to-end case study of monitoring an outbreak in real time: choosing a model when data are scarce, reading the nowcast as the epidemic grows, and experimenting with the prior to encode what you already believe about the epidemic before much data arrives.
Understanding Priors in diseasenowcasting — Make the model say what you mean. Shows the package’s default priors, how to tighten or loosen them, and how the prior trades off against the data. Uses the prior-predictive tools (
nowcast(..., prior_only = TRUE)) to see what a prior implies before fitting.Custom delays and epidemic processes — Two examples on how to set your own delays and epidemic processes. Includes how to use ordinary differential equation models.
Handling Outlier Delays with Censoring — Robustness to reporting glitches. How the censored likelihood copes with unusually long reporting delays, and how to flag extreme delays in your surveillance stream.
Using alongside an LLM — Use AI. How to use the
SKILL.mdto teach a Large Language Model how to you develop your nowcasts withdiseasenowcasting.Benchmark (diseasenowcasting vs NobBS and epinowcast) — How does it compare? A reproducible backtest comparing
diseasenowcastingagainst theNobBSandepinowcastpackages.Mathematical Foundations of diseasenowcasting — Under the hood. The censored likelihood, the epidemic processes (HSGP, AR(1), SIR), the delay families, and the Laplace-approximation inference that powers
RTMB.
See ?nowcast, ?backtest,
?score, and ?model for full documentation.