Robust Rt estimation

We provide simple Julia implementations of two popular methods for reproduction number estimation: EpiEstim (Cori et al. 2013) and EpiFilter (Parag 2021). In addition to the original methods, we also provide functions for calculating model likelihoods and marginalising out uncertainty about model parameters from final estimates of \(R_t\).

These scripts were developed as part of Robust uncertainty quantification in popular estimators of the instantaneous reproduction number. All technical details can be found in this paper and the corresponding supplementary material.

While the paper above focuses on real-time estimation, we also provide a discussion and comparison of retrospective methods.

Code

using Plots, Measures

# Ensure we are working in the root directory
const rootdir = @__DIR__
cd(joinpath(rootdir, ".."))

# Load scripts
include("src/EpiEstim.jl")
include("src/support.jl")

# Set k parameter and load data
k = 7
(Ct, w) = loadData("NZCOVID")
dates = Date("2020-02-26") + Day.(0:99)

# Run the model and extract mean and credible intervals
pR = EpiEstim(k, w, Ct)
m = mean.(pR)
l = quantile.(pR, 0.025)
u = quantile.(pR, 0.975)

# Plot
plot_cases = bar(dates, Ct, xlabel="Date", ylabel="Reported cases", label=false, size=(800, 400), margins=3mm, color="#f25a2a")
plot_Rt = plot(dates, m, ribbon=(m-l,u-m), label=false, xlabel="Date", ylabel="Reproduction number", size=(800,400), margins=3mm, color="#13643f")
hline!(plot_Rt, [1], label=false, color=:black, linestyle=:dash)
plot(plot_cases, plot_Rt, layout=(2,1), size=(800,600))

Reported cases from the first 100 days of the COVID-19 pandemic in Aotearoa New Zealand (upper) and the estimated reproduction number using EpiEstim with a default smoothing value of k = 7 (lower). Data were obtained from the Ministry of Health.

Other than a Julia installation, no external software is required to run this code. All functions are provided in the /src/ directory of the GitHub repository.

References

Cori, Anne, Neil M. Ferguson, Christophe Fraser, and Simon Cauchemez. 2013. “A New Framework and Software to Estimate Time-Varying Reproduction Numbers During Epidemics.” American Journal of Epidemiology 178 (9): 1505–12. https://doi.org/10.1093/aje/kwt133.

Parag, Kris V. 2021. “Improved Estimation of Time-Varying Reproduction Numbers at Low Case Incidence and Between Epidemic Waves.” PLOS Computational Biology 17 (9): e1009347. https://doi.org/10.1371/journal.pcbi.1009347.