---
title: "Introduction to diseasenowcasting: Real-Time Epidemic Nowcasting"
author: "diseasenowcasting team"
output: rmarkdown::html_vignette
vignette: >
  %\VignetteIndexEntry{Introduction to diseasenowcasting: Real-Time Epidemic Nowcasting}
  %\VignetteEngine{knitr::rmarkdown}
  %\VignetteEncoding{UTF-8}
---

```{r setup, include=FALSE}
knitr::opts_chunk$set(
  message  = FALSE,
  warning  = FALSE,
  fig.width  = 7,
  fig.align='center',
  fig.height = 4,
  collapse = TRUE,
  comment  = "#>"
)
```

`diseasenowcasting` 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`](https://cran.r-project.org/web/packages/RTMB/index.html)) to infer the cases that have not yet been reported thus providing a prediction of the final number of cases.

## Native modelling and the shared workflow

Use `diseasenowcasting`'s native interface when specifying the statistical
model: `model()`, `nowcast()`, and `auto_nowcast()` expose its epidemic, delay,
likelihood, revision, and cumulative-process choices. Use `fit_check()` and
`nowcast_diagnostic()` for diagnostics that are specifically about the RTMB fit.

The returned object is already a `tbl.now::tbl_nowcast`. Plotting, tidying,
predictive scoring, ensembling, saving, and loading therefore use the common
result directly. Likewise, `backtest()` is a convenience for translating native
`model()` specifications into `tbl.now` engines; it returns a canonical
`nowcast_backtest`, not a second native backtest class. That result can be passed
directly to the scoringutils `as_forecast_quantile()`, `as_forecast_point()`, or
`as_forecast_sample()` methods. See `?diseasenowcasting_workflows` for the full
boundary and examples.

## Your data: the `tbl_now` format

`diseasenowcasting` works with data organised as a `tbl_now` object from the companion
[`tbl.now`](https://rodrigozepeda.github.io/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)


```{r load}
library(diseasenowcasting)
library(tbl.now)
library(dplyr)
library(ggplot2)
set.seed(27653)
```

As a quick example, here is how to build a `tbl_now` using the following surveillance data 
for dengue in Puerto Rico:

```{r tblnow-demo}
data(denguedat)
```

```{r, echo = FALSE}
head(denguedat, 6)
```


We 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`). 

```{r}
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
```

Once your data is a `tbl_now`, a single call to `nowcast()` does the rest.

> For more information about `tbl_now` check [the package's website](https://rodrigozepeda.github.io/tbl.now/index.html). 

## 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`:

```{r denguestrata}
dengue_tbl <-  dengue_tbl |>  add_strata(gender)
```

Notice that the `tbl_now` automatically prints the `strata` specification below:

```{r}
dengue_tbl
```

We can also add temporal effects for example a weekly seasonality (52 seasons) as well
as a holiday effect using the `almanac` package:

```{r, eval=requireNamespace("almanac", quietly = TRUE)}
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:

```{r dengue-fit}
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`:

```{r dengue-plot, fig.cap="*Nowcast for dengue example. The shaded bars show the median while the errorbar has the 90% credible intervals.*"}
autoplot(nc_dengue) 
```

Values can be obtained via `predict()` and `summary()`:

```{r dengue-results, eval = FALSE}
# Full posterior-predictive nowcast at every event-time
pred_dengue <- predict(nc_dengue)

#This creates a summary of mean and quantiles
summary(pred_dengue) 
```



```{r dengue-results3, echo = FALSE}
# Full posterior-predictive nowcast at every event-time
pred_dengue <- predict(nc_dengue)

#This creates a summary of mean and quantiles
summary(pred_dengue) |> tail(6)
```


Additionally 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):

```{r dengue-diagnostic, fig.height=10}
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.  

```{r mpox-prep}
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:

```{r mpoxplot, fig.height=8}
autoplot(mpox_tbl)
```


We can also set it again with `add_temporal_effects()`:

```{r}
mpox_tbl <- mpox_tbl |> 
  add_temporal_effects(temporal_effects(day_of_week = TRUE))
```

You can see that the `tbl_now` indicates its computation:

```{r, echo = FALSE}
mpox_tbl
```

One 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:

```{r mpox-fit}
#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 per `date` and `model` cell. Those cells run 
in parallel through the [future](https://future.futureverse.org/) framework.  
By default they run sequentially; to use several CPU cores, 
set a [parallel plan](https://future.futureverse.org/reference/plan.html) *before* calling `backtest()`:

```{r backtest-plan, eval=FALSE}
library(future)
plan(multisession, workers = 4)   # 4 parallel R sessions
# ... run backtest() ...
plan(sequential)                  # restore serial execution when done
```

In what follows we run a backtest in sequential mode however we
strongly recommend using as many workers in a multisession plan
as possible:

```{r backtest}
# 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,
  dates = seq(as.Date("2022-08-01"), as.Date("2022-08-15"), by = "week")
)

#This closes the plan multisession opened above
#plan(sequential)   
```

The backtest is already in the common `tbl.now` format. Convert it directly to
a scoringutils quantile forecast and compute any supported metrics:

```{r backtest-score}
backtest_scores <- backtest_mpox |>
  scoringutils::as_forecast_quantile() |>
  scoringutils::score()

relative_scores <- backtest_scores |>
  scoringutils::add_relative_skill(metric = "wis") |>
  scoringutils::summarise_scores(by = "model")
relative_scores
```

The scoringutils output includes WIS and its decomposition, median absolute
error, interval coverage, and relative WIS. A well-calibrated nowcast should
have low WIS and coverage close to the nominal levels.

> Read the table by picking the model with the lowest WIS whose coverage is
still close to nominal -- a model that wins on WIS by being overconfident is not
the one you want in production. Note that for the tutorial we only used 3
historical dates, which is far too few to reach a definite conclusion; a real
comparison would use dozens.

## When a report is not yet a case: the revision process

Everything so far assumes a report **is** a case. Many registers work
provisionally instead: a report arrives, and is later **resolved** — confirmed by
a laboratory result, or retracted when it turns out not to be a case at all.

Record that on the `tbl_now` and `diseasenowcasting` nowcasts the number that
actually settles rather than the raw report count. There is no argument to pass:
the process is **detected** from the data.

```{r resolution-demo, eval = FALSE}
# One date, plus what the result was: "confirmed", "retracted" or "pending".
dat <- tbl_now(linelist, event_date = onset, report_date = reported,
               revision_date = result, revision_type = outcome,
               data_type = "linelist")

nowcast(dat)     # works out for itself which outcomes you record
```

Whether you record only retractions, only confirmations, or both, is read from
the values in `revision_type` — and the answer targets, respectively, the cases
never retracted, the cases eventually confirmed, or the cases whose result comes
back positive.

The key point, and the reason this is not just "multiply by the confirmed
fraction": **a missing resolution date means *not resolved yet*, not "fine"**. A
report filed this morning has had no chance to be retracted; one filed two months
ago and still standing is almost certainly genuine. The model weighs every report
by how long it has had to be contradicted, so recent event times — exactly the
ones you care about — are corrected properly.

`nowcast()` prints what it is modelling and the fitted probability, and
`parameters()` gives its uncertainty (`tidy()` on a nowcast gives you the
predicted counts instead):

```{r resolution-parameters, eval = FALSE}
parameters(nc) |> filter(type == "resolution")
#>                 term estimate conf.low conf.high       type
#>  prob_not_retracted    0.8512   0.8399    0.8619 resolution
```

See the [Revision processes article](https://rodrigozepeda.github.io/diseasenowcasting/articles/Revision_processes.html) for the full treatment, including
per-stratum probabilities, censored revision dates, count-incidence data, and
the shared revision-delay assumption used when both outcomes are recorded.

## 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 the dengue data up to January 1992:

```{r auto-data}
#All dengue observed as of January 1992
dengue_92 <- denguedat |>
  filter(onset_week  <= as.Date("1992-01-01") &
         report_week <= as.Date("1992-01-01"))

dengue_tbl_92 <- tbl_now(
  dengue_92,
  event_date  = onset_week,
  report_date = report_week,
  data_type   = "linelist",
  now         = as.Date("1992-01-01")
)
```

Backtesting the whole grid is the expensive step: cost grows with the number of
candidates (`delays` x epidemic processes), the number of backtest dates
(`n_dates`), and -- most steeply -- the length of the series. The settings below
are deliberately small so this vignette builds quickly; for a real selection keep
the defaults (all three delay families, `n_dates = 6`, `n_draws = 2000`,
`K = 25`). 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):

```{r auto-fit}
# Uncomment to run candidates in parallel:
# library(future)
# plan(multisession, workers = max(parallel::detectCores() - 1, 1))
auto_ncast <- auto_nowcast(
  dengue_tbl_92,
  metric         = "wis",   # rank candidates by relative WIS (the default)
  relative_score = TRUE,
  delays         = list(lognormal_delay(), generalized_gamma_delay()),
  n_dates        = 3,       # backtest at 3 historical dates (raise for a firmer choice)
  n_draws_select = 100,     # draws while comparing (small => fast)
  n_draws        = 300,     # draws for the final fit of the winner
  K              = 5,       # delay imputations for the final fit
  K_select       = 5        # delay imputations while comparing
)
# plan(sequential)
```

We can show the scores of the models to see the best performer:

```{r auto-result}
comparison_scores(auto_ncast)  # every candidate, ranked best-first
selection_timings(auto_ncast)  # retrospective fits, refits, and total seconds
```

`best_model()` hands back the winning `model()` object, so you can reuse the same
specification on other data (or feed it to `nowcast()` / `backtest()`):

```{r auto-best-model}
winner <- best_model(auto_ncast)
winner
```

Because the result is a normal nowcast, everything else just works:

```{r auto-plot, fig.cap="_Nowcast from the model auto_nowcast() selected._"}
autoplot(auto_ncast)
```

### 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:

```{r auto-update, warning=TRUE}
# A few more weeks of dengue, observed as of 1992-04-01:
dengue_apr <- denguedat |>
  filter(onset_week  <= as.Date("1992-04-01") &
         report_week <= as.Date("1992-04-01"))

dengue_tbl_apr <- tbl_now(
  dengue_apr,
  event_date  = onset_week,
  report_date = report_week,
  data_type   = "linelist",
  now         = as.Date("1992-04-01")
)

auto_ncast_updated <- update(auto_ncast, dengue_tbl_apr)
```

Any reports with surprising delays are collected by `extreme_values()` (it
returns `NULL` when nothing looks off):

```{r auto-extreme}
extreme_values(auto_ncast_updated)
```

See the article on [Handling Outlier Delays with Censoring](https://rodrigozepeda.github.io/diseasenowcasting/articles/Handling_Outlier_Delays_with_Censoring.html) for what to do when a
delay *is* flagged.

## Example 5 -- Count-cumulative revisions and a historical origin

FluSight publishes a running level for each target week. Keep three dates
distinct: `full_data` is the retrospective extract used only to obtain scoring
truth; `historical_now` is the information cutoff; and `event_date` identifies
the weeks being nowcast. Starting near September 2023 avoids interpreting the
between-season gap as a reporting delay.

```{r flusight-data, eval = FALSE}
data(flusight, package = "tbl.now")
historical_now <- as.Date("2024-01-27")
full_data <- flusight |>
  filter(location_name == "California",
         target_end_date >= as.Date("2023-09-01"))

calendar_tbl <- tbl_now(
  full_data,
  event_date = target_end_date,
  report_date = as_of,
  case_count = observation,
  data_type = "count-cumulative",
  event_units = "weeks", report_units = "weeks",
  now = max(full_data$as_of), verbose = FALSE
)
calendar_rows <- full_data |>
  filter(target_end_date <= historical_now, as_of <= historical_now)
calendar_asof <- tbl_now(
  calendar_rows,
  event_date = target_end_date, report_date = as_of,
  case_count = observation, data_type = "count-cumulative",
  event_units = "weeks", report_units = "weeks",
  now = historical_now, verbose = FALSE
) |>
  tbl.now::complete_zeroes(max_delay = 26L, until = historical_now)

stopifnot(max(calendar_asof$target_end_date) <= historical_now,
          max(calendar_asof$as_of) <= historical_now)
```

The **calendar clock** keeps every intervening calendar week, and zero completion
makes missing cells inside the observable triangle explicit. For the
**compressed publication clock**, map the observed publication weeks to
consecutive synthetic weeks before constructing the `tbl_now`; retain the real
dates in separate columns for the as-of filter and scoring join.

```{r flusight-compressed, eval = FALSE}
publication_weeks <- sort(unique(full_data$as_of))
publication_index <- setNames(seq_along(publication_weeks) - 1L,
                              as.character(publication_weeks))
compressed_data <- full_data |>
  filter(target_end_date %in% publication_weeks) |>
  mutate(event_num = publication_index[as.character(target_end_date)],
         report_num = publication_index[as.character(as_of)],
         event_model = as.Date("2000-01-01") + 7L * event_num,
         report_model = as.Date("2000-01-01") + 7L * report_num)
compressed_tbl <- tbl_now(
  compressed_data, event_date = event_model, report_date = report_model,
  case_count = observation, data_type = "count-cumulative",
  event_units = "weeks", report_units = "weeks", verbose = FALSE
)
# The final publication of one season is immediately followed by the first of
# the next on this model clock; the original dates still enforce historical_now.
```

Count-cumulative data use their own component. `settlement = 26L` targets
$C_t(26)$, finite-horizon database retention. The cumulative-level model and the
two signed hurdle models share the same collapsed retraction kernel; the
zero-truncated Poisson magnitude has no dispersion parameter.

```{r flusight-nowcast, eval = FALSE}
level_model <- model(
  nb_likelihood(), ar1_epidemic(), lognormal_delay(),
  cumulative = cumulative_process(
    observation = "cumulative", settlement = 26L
  )
)
ztnb_model <- model(
  nb_likelihood(), ar1_epidemic(), lognormal_delay(),
  cumulative = cumulative_process(
    observation = "hurdle_ztnb", settlement = 26L
  )
)
ztp_model <- model(
  poisson_likelihood(), ar1_epidemic(), lognormal_delay(),
  cumulative = cumulative_process(
    observation = "hurdle_ztpoisson", settlement = 26L
  )
)

calendar_fits <- lapply(
  list(level = level_model, ztnb = ztnb_model, ztp = ztp_model),
  nowcast, data = calendar_tbl, now = historical_now,
  temporal_effects = "none"
)
```

An empirical multiplier is a development comparator, not part of the fitted
model. At target age $a$, calibrate it only on earlier cohorts for which both
$C_s(a)$ and $C_s(26)$ were observable by `historical_now`; discard zero
denominators and report a fallback when too few pairs remain. Retrospective
terminal values after `historical_now` may be used for scoring, never for this
calibration. See `devel/skellam_prototypes/flusight_asof_data.R` for the complete
calendar/compressed preparation and leakage-safe multiplier used by the large
backtest runner.

## 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):

```{r save-load}
saved <- tempfile(fileext = ".rds")
save_nowcast(auto_ncast, saved)

restored <- load_nowcast(saved)
# predict() / autoplot() / coef() / parameters() 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):

```{r save-refit, eval = FALSE}
nowcast(restored@data, restored@model)   # re-runs the optimisation
```

## Combining diseasenowcasting models in an ensemble

Different epidemic processes can fail in different ways. For example, an HSGP
can follow a smooth epidemic curve while an AR(1) process reacts more locally.
Because `nowcast()` returns the common `tbl.now::tbl_nowcast` result, its output
can be passed **directly** to `tbl.now::nowcast_ensemble()`; there is no need to
refit the models through `tbl.now::run_nowcast()` or convert their results.

Here both members use the same likelihood and reporting-delay model, and differ
only in their epidemic process:

```{r ensemble-models, eval = FALSE}
hsgp_fit <- nowcast(
  dengue_tbl,
  model = model(
    likelihood = nb_likelihood(),
    epidemic   = hsgp_epidemic(),
    delay      = lognormal_delay()
  ),
  n_draws = 1000,
  seed = 101
)

ar1_fit <- nowcast(
  dengue_tbl,
  model = model(
    likelihood = nb_likelihood(),
    epidemic   = ar1_epidemic(),
    delay      = lognormal_delay()
  ),
  n_draws = 1000,
  seed = 102
)
```

The default ensemble averages matching predictive quantiles from the two
models. Naming the arguments also records useful member names in the result:

```{r ensemble-quantiles, eval = FALSE}
ensemble <- tbl.now::nowcast_ensemble(
  HSGP = hsgp_fit,
  AR1  = ar1_fit
)

autoplot(ensemble)
```

Both `diseasenowcasting` fits retain posterior draws, so they can also be
combined as a mixture distribution. This preserves disagreement between the
members as part of the ensemble uncertainty and will generally give wider
intervals than averaging their quantiles:

```{r ensemble-linear-pool, eval = FALSE}
pooled_ensemble <- tbl.now::nowcast_ensemble(
  HSGP = hsgp_fit,
  AR1  = ar1_fit,
  type = "linear_pool",
  n_draws = 4000
)
```

## Next steps

This vignette covered the basics: building a `tbl_now`, fitting a nowcast with
`nowcast()`, inspecting results with `predict()` / `summary()` / `autoplot()`,
comparing models with `backtest()` / `score()`, and combining structurally
different fits with `tbl.now::nowcast_ensemble()`.

Depending on what you want to do next, check out the following
vignettes and website articles:

- **[Nowcasting at the Start of an Epidemic](Nowcasting_at_the_start_of_an_Epidemic.html)**
  —  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](Understanding_Priors.html)**
  — *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](Custom_delays_and_processes.html)**
  — 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](https://rodrigozepeda.github.io/diseasenowcasting/articles/Handling_Outlier_Delays_with_Censoring.html)**
  — *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](https://rodrigozepeda.github.io/diseasenowcasting/articles/LLM_Usage.html)**
  — *Use AI.* How to use the [`SKILL.md`](https://github.com/RodrigoZepeda/diseasenowcasting/blob/master/SKILL.md)
  to teach a Large Language Model how to you develop your nowcasts with `diseasenowcasting`. 

- **[Benchmark (diseasenowcasting vs NobBS and epinowcast)](https://rodrigozepeda.github.io/diseasenowcasting/articles/Benchmark.html)**
  — *How does it compare?* A reproducible backtest comparing `diseasenowcasting`
  against the `NobBS` and `epinowcast` packages.

- **[Mathematical Foundations of diseasenowcasting](https://rodrigozepeda.github.io/diseasenowcasting/articles/Mathematics.html)**
  — *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.