Replicating EpiNow2 Models

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 $R_t$ from case data.

Test Data

We generate simulated epidemic data with an intervention effect at day 25:

# Generate test 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

First, we run EpiNow2. Note: EpiNow2 must be loaded and run before EpiAwareR due to a runtime conflict between Stan and Julia.

library(EpiNow2)

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

# No reporting delay for this simple example (data is directly observed)
# For real data, you would add: delays = delay_opts(reporting_delay)

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

# Run estimation with timing
# Wrap in tryCatch to handle CI environments where Stan may fail
epinow2_results <- NULL
epinow2_time <- tryCatch({
  system.time({
    epinow2_results <- epinow(
      outbreak_data,
      generation_time = generation_time_opts(generation_time),
      rt = rt_settings,
      stan = stan_opts(cores = 2, chains = 2, samples = 500, warmup = 250)
    )
  })
}, error = function(e) {
  message("EpiNow2 fitting failed (may occur in some CI environments): ", e$message)
  NULL
})

if (!is.null(epinow2_results)) {
  cat("EpiNow2 runtime:", epinow2_time["elapsed"], "seconds\n")
  summary(epinow2_results)
  plot(epinow2_results)
} else {
  cat("EpiNow2 results not available - see message above\n")
}
#> EpiNow2 runtime: 26.802 seconds

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

Step 2: The EpiAwareR Equivalent

Now we load EpiAwareR and replicate the analysis. We explicitly build the model from components:

library(EpiAwareR)

# 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(log(1.5), 0.3)), # log(Rt) centered at 1.5 for growing epidemic
  std_prior = halfnorm(0.1)
)
#> Starting Julia ...

# 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(30), 0.5) # Match simulated initial infections
)

# 3. Observation model: negative binomial (no delay for this simple example)
observation <- NegativeBinomialError(
  cluster_factor_prior = halfnorm(0.5)
)

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

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
#>   Running Pathfinder initialization...
#>   Pathfinder init failed, using default initialization...
#> Processing results...

cat("EpiAwareR runtime:", epiaware_time["elapsed"], "seconds\n")
#> EpiAwareR runtime: 565.882 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_negbin 
#> 
#> Sampling:
#>   Method: NUTS
#>   Chains: 2 
#>   Draws: 500 (per chain)
#> 
#> Convergence:
#>   Max Rhat: 1.045 
#>   Min ESS (bulk): 68 
#>   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
#>    <chr>              <dbl>   <dbl>   <dbl>   <dbl>   <dbl> <dbl> <dbl>    <dbl>
#>  1 latent.ar_init.…  0.345   0.348  0.261   0.269   -0.0954 0.780 1.00      526.
#>  2 latent.damp_AR.…  0.986   0.987  0.00795 0.00794  0.972  0.998 0.999     644.
#>  3 latent.std        0.184   0.183  0.0250  0.0240   0.145  0.227 1.02      170.
#>  4 latent.ϵ_t.1.    -0.535  -0.515  0.819   0.783   -1.91   0.857 1.000     842.
#>  5 latent.ϵ_t.2.     0.0213  0.0188 0.797   0.825   -1.29   1.29  1.000    1311.
#>  6 latent.ϵ_t.3.    -0.905  -0.904  0.801   0.776   -2.25   0.386 1.00      877.
#>  7 latent.ϵ_t.4.     0.334   0.329  0.812   0.800   -1.01   1.66  1.00     1139.
#>  8 latent.ϵ_t.5.     0.0533  0.0645 0.763   0.759   -1.22   1.28  1.01      870.
#>  9 latent.ϵ_t.6.     0.519   0.526  0.761   0.735   -0.748  1.75  1.00      853.
#> 10 latent.ϵ_t.7.     1.67    1.68   0.713   0.706    0.504  2.83  1.00     1035.
#> # ℹ 56 more rows
#> # ℹ 1 more variable: ess_tail <dbl>
plot(results, type = "Rt")

Note on initial uncertainty: Renewal models often show high uncertainty in the initial period (first ~7 days) due to limited infection history for estimating $R_t$ . This “burn-in” effect is a known characteristic of these models.

Runtime Comparison

cat("Runtime comparison (2 chains each):\n")
#> Runtime comparison (2 chains each):
if (!is.null(epinow2_time)) {
  cat("  EpiNow2:   ", round(epinow2_time["elapsed"], 1), "seconds\n")
  cat("  EpiAwareR: ", round(epiaware_time["elapsed"], 1), "seconds\n")
  cat("  Speedup:   ", round(epinow2_time["elapsed"] / epiaware_time["elapsed"], 1), "x\n")
} else {
  cat("  EpiNow2:    (not available)\n")
  cat("  EpiAwareR: ", round(epiaware_time["elapsed"], 1), "seconds\n")
}
#>   EpiNow2:    26.8 seconds
#>   EpiAwareR:  565.9 seconds
#>   Speedup:    0 x

Benefits of the explicit approach

EpiAwareR’s compositional design makes the model structure transparent:

Latent model explicitly states how $R_t$ evolves (random walk, AR, MA, etc.)
Infection model clearly specifies generation time and renewal dynamics
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(0.99, 0.01, 0.9, 1)),
  init_priors = list(norm(0, 0.5)),
  std_prior = halfnorm(0.3) # Larger innovations
)

# More smooth (less variable Rt)
rw_smooth <- AR(
  order = 1,
  damp_priors = list(truncnorm(0.99, 0.01, 0.9, 1)),
  init_priors = list(norm(0, 0.5)),
  std_prior = halfnorm(0.05) # 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

EpiNow2 Documentation - The original package
Getting Started - EpiAwareR basics
Mishra Case Study - Detailed worked example

Session Info

sessionInfo()
#> R version 4.6.1 (2026-06-24)
#> Platform: x86_64-pc-linux-gnu
#> Running under: Ubuntu 24.04.4 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           LC_TIME=C.UTF-8       
#>  [4] LC_COLLATE=C.UTF-8     LC_MONETARY=C.UTF-8    LC_MESSAGES=C.UTF-8   
#>  [7] LC_PAPER=C.UTF-8       LC_NAME=C              LC_ADDRESS=C          
#> [10] LC_TELEPHONE=C         LC_MEASUREMENT=C.UTF-8 LC_IDENTIFICATION=C   
#> 
#> time zone: UTC
#> tzcode source: system (glibc)
#> 
#> attached base packages:
#> [1] stats     graphics  grDevices utils     datasets  methods   base     
#> 
#> other attached packages:
#> [1] EpiAwareR_0.2.0.9000 EpiNow2_1.9.0       
#> 
#> loaded via a namespace (and not attached):
#>  [1] gtable_0.3.6          tensorA_0.36.2.1      xfun_0.60            
#>  [4] bslib_0.12.0          ggplot2_4.0.3         QuickJSR_1.11.0      
#>  [7] inline_0.3.21         vctrs_0.7.3           tools_4.6.1          
#> [10] generics_0.1.4        curl_7.1.0            parallel_4.6.1       
#> [13] stats4_4.6.1          tibble_3.3.1          pkgconfig_2.0.3      
#> [16] R.oo_1.27.1           data.table_1.18.6.1   checkmate_2.3.4      
#> [19] RColorBrewer_1.1-3    S7_0.2.2              desc_1.4.3           
#> [22] distributional_0.8.1  RcppParallel_6.2.0    truncnorm_1.0-9      
#> [25] lifecycle_1.0.5       compiler_4.6.1        farver_2.1.2         
#> [28] textshaping_1.0.5     codetools_0.2-20      htmltools_0.5.9      
#> [31] sass_0.4.10           yaml_2.3.12           pillar_1.11.1        
#> [34] pkgdown_2.2.1         jquerylib_0.1.4       R.utils_2.13.0       
#> [37] cachem_1.1.0          StanHeaders_2.32.10   abind_1.4-8          
#> [40] posterior_1.7.0       rstan_2.32.7          tidyselect_1.2.1     
#> [43] digest_0.6.39         dplyr_1.2.1           purrr_1.2.2          
#> [46] labeling_0.4.3        fastmap_1.2.0         grid_4.6.1           
#> [49] primarycensored_1.5.1 cli_3.6.6             magrittr_2.0.5       
#> [52] loo_2.10.1            patchwork_1.3.2       utf8_1.2.6           
#> [55] pkgbuild_1.4.8        withr_3.0.3           runner_0.4.6         
#> [58] scales_1.4.0          backports_1.5.1       lubridate_1.9.5      
#> [61] timechange_0.4.0      rmarkdown_2.31        lambda.r_1.2.4       
#> [64] matrixStats_1.5.0     otel_0.2.0            gridExtra_2.3.1      
#> [67] JuliaConnectoR_1.1.5  futile.logger_1.4.9   ragg_1.5.2           
#> [70] R.methodsS3_1.8.2     juliaready_0.1.0      evaluate_1.0.5       
#> [73] knitr_1.51            V8_8.2.0              rstantools_2.7.0     
#> [76] rlang_1.3.0           futile.options_1.0.1  Rcpp_1.1.2           
#> [79] glue_1.8.1            formatR_1.14          jsonlite_2.0.0       
#> [82] R6_2.6.1              systemfonts_1.3.2     fs_2.1.0

Sebastian Funk

2026-08-25