Replicating EpiNow2 Models

Sebastian Funk

2025-12-12

Source: vignettes/epinow2-comparison.Rmd

Introduction

EpiNow2 is a popular R package for estimating time-varying reproduction numbers using renewal equation models. This vignette shows how to replicate a typical EpiNow2 analysis using EpiAwareR’s compositional approach.

Both packages use renewal equations and Bayesian inference, but EpiAwareR’s compositional design allows you to build a much wider class of models by assembling components.

Example: A Typical EpiNow2 Analysis

We’ll replicate a standard EpiNow2 workflow for estimating RtR_t from case data.

set.seed(123)
dates <- seq.Date(as.Date("2024-01-01"), by = "day", length.out = 50)

# Epidemic with intervention at day 25
infections <- numeric(50)
infections[1:7] <- 30
for (i in 8:50) {
  rt <- if (i < 25) 2.0 else 0.7 # Intervention effect
  infections[i] <- rpois(1, rt * mean(infections[max(1, i - 7):(i - 1)]))
}

outbreak_data <- data.frame(
  date = dates,
  confirm = as.integer(pmax(round(infections + rnorm(50, 0, 5)), 1)) # Integer counts
)

head(outbreak_data, 10)
#>          date confirm
#> 1  2024-01-01      37
#> 2  2024-01-02      29
#> 3  2024-01-03      38
#> 4  2024-01-04      22
#> 5  2024-01-05      33
#> 6  2024-01-06      31
#> 7  2024-01-07      31
#> 8  2024-01-08      57
#> 9  2024-01-09      73
#> 10 2024-01-10      63

Step 1: The EpiNow2 Approach

A typical EpiNow2 analysis with a daily random walk looks like this:

library(EpiNow2)
#> 
#> Attaching package: 'EpiNow2'
#> The following object is masked from 'package:stats':
#> 
#>     Gamma

# Define generation time
generation_time <- Gamma(
  mean = 5.0,
  sd = 2.0
)

# Define reporting delay
reporting_delay <- LogNormal(
  mean = 2.0,
  sd = 1.0
)

# Use daily random walk for Rt (alternative to GP)
rt_settings <- rt_opts(
  prior = Normal(mean = 1, sd = 1),
  rw = 1 # Daily random walk
)

# Run estimation with timing
epinow2_time <- system.time({
  epinow2_results <- epinow(
    outbreak_data,
    generation_time = generation_time_opts(generation_time),
    delays = delay_opts(reporting_delay),
    rt = rt_settings,
    stan = stan_opts(cores = 2, chains = 2, samples = 500, warmup = 250)
  )
})
#> Logging threshold set at INFO for the name logger
#> Writing EpiNow2 logs to the console and:
#> /tmp/Rtmpru51SC/regional-epinow/2024-02-19.log.
#> Logging threshold set at INFO for the name logger
#> Writing EpiNow2.epinow logs to the console and:
#> /tmp/Rtmpru51SC/epinow/2024-02-19.log.
#> ℹ Unconstrained distributon passed as a delay.
#> ℹ Constraining with default CDF cutoff 0.001.
#> ℹ To silence this message, specify delay distributions with `max` or
#>   `default_cdf_cutoff`.
#> ℹ Unconstrained distributon passed as a delay.
#> ℹ Constraining with default CDF cutoff 0.001.
#> ℹ To silence this message, specify delay distributions with `max` or
#>   `default_cdf_cutoff`.
#> WARN [2025-12-12 13:21:46] epinow: Bulk Effective Samples Size (ESS) is too low, indicating posterior means and medians may be unreliable.
#> Running the chains for more iterations may help. See
#> https://mc-stan.org/misc/warnings.html#bulk-ess - 
#> WARN [2025-12-12 13:21:46] epinow: Tail Effective Samples Size (ESS) is too low, indicating posterior variances and tail quantiles may be unreliable.
#> Running the chains for more iterations may help. See
#> https://mc-stan.org/misc/warnings.html#tail-ess -

cat("EpiNow2 runtime:", epinow2_time["elapsed"], "seconds\n")
#> EpiNow2 runtime: 48.18 seconds

# View results
summary(epinow2_results)
#>                         measure               estimate
#>                          <char>                 <char>
#> 1:       New infections per day          53 (33 -- 87)
#> 2:   Expected change in reports             Decreasing
#> 3:   Effective reproduction no.    0.59 (0.37 -- 0.94)
#> 4:               Rate of growth -0.091 (-0.47 -- 0.23)
#> 5: Doubling/halving time (days)       -7.7 (3 -- -1.5)
plot(epinow2_results)

The rw = 1 option creates a daily random walk for RtR_t, which is an alternative choice to the default Gaussian Process.

Step 2: The EpiAwareR Equivalent

To replicate this in EpiAwareR, we explicitly build the model from components:

# 1. Latent process: daily random walk for Rt
# AR(1) with damping ≈ 1 creates a random walk (equivalent to EpiNow2's rw=1)
latent <- AR(
  order = 1,
  damp_priors = list(truncnorm(0.99, 0.01, 0.9, 1)), # Near 1 for random walk
  init_priors = list(norm(0, 1)), # log(Rt) prior: Rt centered at 1
  std_prior = halfnorm(0.1)
)

# 2. Infection model: renewal equation with generation time
# EpiNow2: Gamma(mean=5, sd=2) -> shape=6.25, scale=0.8
infection <- Renewal(
  gen_distribution = gamma_dist(6.25, 0.8),
  initialisation_prior = norm(log(100), 2)
)

# 3. Observation model: negative binomial with reporting delay
# EpiNow2: LogNormal(mean=2, sd=1) -> meanlog=0.58, sdlog=0.47
observation <- LatentDelay(
  model = NegativeBinomialError(
    cluster_factor_prior = halfnorm(0.5)
  ),
  delay_distribution = lognorm(0.58, 0.47)
)

# 4. Compose into complete model
model <- EpiProblem(
  epi_model = infection,
  latent_model = latent,
  observation_model = observation,
  tspan = c(1, 50)
)

print(model)
#> <EpiAware Epidemiological Model>
#>   Time span: 1 to 50 
#>   Components:
#>     - Infection model: epiaware_renewal 
#>     - Latent model: epiaware_ar 
#>     - Observation model: epiaware_delay

Fit the Model 

# Run estimation with timing (2 chains to match EpiNow2)
epiaware_time <- system.time({
  results <- fit(
    model = model,
    data = outbreak_data,
    method = nuts_sampler(
      warmup = 250,
      draws = 500,
      chains = 2
    )
  )
})
#> Generating Turing.jl model...
#> Running NUTS sampling...
#>   Chains: 2
#>   Warmup: 250
#>   Draws: 500
#> Processing results...

cat("EpiAwareR runtime:", epiaware_time["elapsed"], "seconds\n")
#> EpiAwareR runtime: 633.585 seconds

# View results
print(results)
#> <EpiAware Model Fit>
#> 
#> Model:
#>   Time span: 1 to 50 
#>   Infection model: epiaware_renewal 
#>   Latent model: epiaware_ar 
#>   Observation model: epiaware_delay 
#> 
#> Sampling:
#>   Method: NUTS
#>   Chains: 2 
#>   Draws: 500 (per chain)
#> 
#> Convergence:
#>   Max Rhat: 1.05 
#>   Min ESS (bulk): 2 
#>   Warning: Some parameters have ESS < 100
#> 
#> Use summary() for parameter estimates
#> Use plot() to visualize results
summary(results)
#> # A tibble: 66 × 10
#>    variable    mean  median      sd     mad     q5   q95  rhat ess_bulk ess_tail
#>    <chr>      <dbl>   <dbl>   <dbl>   <dbl>  <dbl> <dbl> <dbl>    <dbl>    <dbl>
#>  1 latent.a… 0.601  0.630   0.513   0.488   -0.318 1.40   1.00     695.     512.
#>  2 latent.d… 0.985  0.986   0.00842 0.00878  0.970 0.998  1.01    1075.     485.
#>  3 latent.s… 0.195  0.191   0.0374  0.0371   0.138 0.261  1.01     341.     655.
#>  4 latent.ϵ… 0.168  0.174   0.966   0.967   -1.46  1.79   1.00     981.     877.
#>  5 latent.ϵ… 0.111  0.0948  0.978   1.03    -1.50  1.66   1.00    1545.     748.
#>  6 latent.ϵ… 0.0194 0.00479 0.993   0.878   -1.67  1.64   1.00    2728.     493.
#>  7 latent.ϵ… 0.0626 0.0542  0.977   0.990   -1.57  1.66   1.00    1566.     702.
#>  8 latent.ϵ… 0.500  0.508   0.877   0.892   -0.885 1.98   1.00    1484.     770.
#>  9 latent.ϵ… 0.656  0.667   0.910   0.818   -0.956 2.15   1.00    1510.     573.
#> 10 latent.ϵ… 0.0121 0.0224  0.863   0.839   -1.39  1.38   1.00    1283.     845.
#> # ℹ 56 more rows
plot(results, type = "Rt")

Runtime Comparison 

cat("Runtime comparison (2 chains each):\n")
#> Runtime comparison (2 chains each):
cat("  EpiNow2:   ", round(epinow2_time["elapsed"], 1), "seconds\n")
#>   EpiNow2:    48.2 seconds
cat("  EpiAwareR: ", round(epiaware_time["elapsed"], 1), "seconds\n")
#>   EpiAwareR:  633.6 seconds
cat("  Speedup:   ", round(epinow2_time["elapsed"] / epiaware_time["elapsed"], 1), "x\n")
#>   Speedup:    0.1 x

Benefits of the explicit approach

EpiAwareR’s compositional design makes the model structure transparent:

  1. Latent model explicitly states how RtR_t evolves (random walk, AR, MA, etc.)
  2. Infection model clearly specifies generation time and renewal dynamics
  3. Observation model shows delays and overdispersion assumptions

This explicitness enables:

  • Testing alternatives: Easily compare different smoothness levels, or swap e.g. for AR(2), moving average or mechanistic models (e.g. with susceptible depletion)
  • Understanding priors: See exactly what distributions drive each component
  • Model comparison: Systematically compare competing assumptions

For example, comparing different levels of smoothness in the random walk:

# Less smooth (more variable Rt)
rw_flexible <- AR(
  order = 1,
  damp_priors = list(truncnorm(1.0, 0.01, 0.99, 1)),
  init_priors = list(norm(0, 1)),
  std_prior = halfnorm(0.3) # Larger innovations
)

# More smooth (less variable Rt)
rw_smooth <- AR(
  order = 1,
  damp_priors = list(truncnorm(1.0, 0.01, 0.99, 1)),
  init_priors = list(norm(0, 1)),
  std_prior = halfnorm(0.1) # Smaller innovations
)

# Fit both and compare
model_flexible <- EpiProblem(..., latent_model = rw_flexible, ...)
model_smooth <- EpiProblem(..., latent_model = rw_smooth, ...)

results_flexible <- fit(model_flexible, outbreak_data)
results_smooth <- fit(model_smooth, outbreak_data)

# Visual inspection of results
# Formal model comparison would require LOO cross-validation (for fit)
# or leave-future-out validation (for forecasting)
print(results_flexible)
print(results_smooth)

Learn More

Session Info 

sessionInfo()
#> R version 4.5.2 (2025-10-31)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 24.04.3 LTS
#> 
#> Matrix products: default
#> BLAS:   /usr/lib/x86_64-linux-gnu/openblas-pthread/libblas.so.3 
#> LAPACK: /usr/lib/x86_64-linux-gnu/openblas-pthread/libopenblasp-r0.3.26.so;  LAPACK version 3.12.0
#> 
#> locale:
#>  [1] LC_CTYPE=C.UTF-8          LC_NUMERIC=C             
#>  [3] LC_TIME=C.UTF-8           LC_COLLATE=C.UTF-8       
#>  [5] LC_MONETARY=C.UTF-8       LC_MESSAGES=C.UTF-8      
#>  [7] LC_PAPER=C.UTF-8          LC_NAME=C.UTF-8          
#>  [9] LC_ADDRESS=C.UTF-8        LC_TELEPHONE=C.UTF-8     
#> [11] LC_MEASUREMENT=C.UTF-8    LC_IDENTIFICATION=C.UTF-8
#> 
#> time zone: UTC
#> tzcode source: system (glibc)
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] EpiNow2_1.7.1        EpiAwareR_0.1.0.9000
#> 
#> loaded via a namespace (and not attached):
#>  [1] gtable_0.3.6          tensorA_0.36.2.1      xfun_0.54            
#>  [4] bslib_0.9.0           ggplot2_4.0.1         QuickJSR_1.8.1       
#>  [7] inline_0.3.21         vctrs_0.6.5           tools_4.5.2          
#> [10] generics_0.1.4        parallel_4.5.2        stats4_4.5.2         
#> [13] tibble_3.3.0          pkgconfig_2.0.3       R.oo_1.27.1          
#> [16] data.table_1.17.8     checkmate_2.3.3       RColorBrewer_1.1-3   
#> [19] S7_0.2.1              desc_1.4.3            distributional_0.5.0 
#> [22] RcppParallel_5.1.11-1 truncnorm_1.0-9       lifecycle_1.0.4      
#> [25] compiler_4.5.2        farver_2.1.2          textshaping_1.0.4    
#> [28] codetools_0.2-20      htmltools_0.5.9       sass_0.4.10          
#> [31] yaml_2.3.12           pillar_1.11.1         pkgdown_2.2.0        
#> [34] jquerylib_0.1.4       R.utils_2.13.0        cachem_1.1.0         
#> [37] StanHeaders_2.32.10   JuliaCall_0.17.6      abind_1.4-8          
#> [40] posterior_1.6.1       rstan_2.32.7          tidyselect_1.2.1     
#> [43] digest_0.6.39         dplyr_1.1.4           purrr_1.2.0          
#> [46] labeling_0.4.3        fastmap_1.2.0         grid_4.5.2           
#> [49] cli_3.6.5             magrittr_2.0.4        loo_2.8.0            
#> [52] patchwork_1.3.2       utf8_1.2.6            pkgbuild_1.4.8       
#> [55] withr_3.0.2           runner_0.4.5          scales_1.4.0         
#> [58] backports_1.5.0       lubridate_1.9.4       timechange_0.3.0     
#> [61] rmarkdown_2.30        lambda.r_1.2.4        matrixStats_1.5.0    
#> [64] gridExtra_2.3         futile.logger_1.4.3   ragg_1.5.0           
#> [67] R.methodsS3_1.8.2     evaluate_1.0.5        knitr_1.50           
#> [70] rstantools_2.5.0      rlang_1.1.6           futile.options_1.0.1 
#> [73] Rcpp_1.1.0            glue_1.8.0            formatR_1.14         
#> [76] jsonlite_2.0.0        R6_2.6.1              systemfonts_1.3.1    
#> [79] fs_1.6.6