vignettes/epinow2-comparison.Rmd
epinow2-comparison.RmdEpiNow2 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.
We’ll replicate a standard EpiNow2 workflow for estimating 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 = pmax(infections + rnorm(50, 0, 5), 1) # Add noise
)
head(outbreak_data, 10)
#> date confirm
#> 1 2024-01-01 36.84301
#> 2 2024-01-02 28.87115
#> 3 2024-01-03 37.58235
#> 4 2024-01-04 22.25624
#> 5 2024-01-05 32.92307
#> 6 2024-01-06 30.61927
#> 7 2024-01-07 31.07971
#> 8 2024-01-08 56.89820
#> 9 2024-01-09 73.48838
#> 10 2024-01-10 63.33396A typical EpiNow2 analysis with a daily random walk looks like this:
library(EpiNow2)
# 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
epinow2_results <- epinow(
reported_cases = outbreak_data,
generation_time = generation_time_opts(generation_time),
delays = delay_opts(reporting_delay),
rt = rt_settings,
stan = stan_opts(cores = 4, chains = 4)
)
# View results
summary(epinow2_results)
plot(epinow2_results)The rw = 1 option creates a daily random walk for
,
which is an alternative choice to the default Gaussian Process.
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(1.0, 0.01, 0.99, 1)), # Near 1 = random walk
init_priors = list(norm(0, 1)),
std_prior = halfnorm(0.2)
)
# 2. Infection model: renewal equation with generation time
infection <- Renewal(
gen_distribution = gamma_dist(5.0, 0.4), # Mean 5, SD ~2
initialisation_prior = norm(log(30), 0.5)
)
# 3. Observation model: negative binomial with reporting delay
observation <- LatentDelay(
model = NegativeBinomialError(
cluster_factor_prior = halfnorm(0.2)
),
delay_distribution = lognorm(log(2), 0.5) # 2-day delay
)
# 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
results <- fit(
model = model,
data = outbreak_data,
method = nuts_sampler(
warmup = 500,
draws = 500,
chains = 4
)
)
# View results
print(results)
summary(results)
plot(results, type = "Rt")EpiAwareR’s compositional design makes the model structure transparent:
This explicitness enables:
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)
# Compare using LOO
library(loo)
loo_compare(loo(results_flexible$samples), loo(results_smooth$samples))
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 datasets utils methods base
#>
#> other attached packages:
#> [1] EpiAwareR_0.1.0.9000
#>
#> loaded via a namespace (and not attached):
#> [1] cli_3.6.5 knitr_1.50 rlang_1.1.6 xfun_0.54
#> [5] renv_1.1.5 textshaping_1.0.4 jsonlite_2.0.0 backports_1.5.0
#> [9] htmltools_0.5.9 ragg_1.5.0 sass_0.4.10 rmarkdown_2.30
#> [13] evaluate_1.0.5 jquerylib_0.1.4 fastmap_1.2.0 yaml_2.3.11
#> [17] lifecycle_1.0.4 compiler_4.5.2 fs_1.6.6 Rcpp_1.1.0
#> [21] JuliaCall_0.17.6 systemfonts_1.3.1 digest_0.6.39 R6_2.6.1
#> [25] bslib_0.9.0 checkmate_2.3.3 tools_4.5.2 pkgdown_2.2.0
#> [29] cachem_1.1.0 desc_1.4.3