Estimation and worldwide monitoring of the effective reproductive number of SARS-CoV-2

Abstract

The effective reproductive number R_e is a key indicator of the growth of an epidemic. Since the start of the SARS-CoV-2 pandemic, many methods and online dashboards have sprung up to monitor this number through time. However, these methods are not always thoroughly tested, correctly placed in time, or are overly confident during high incidence periods. Here, we present a method for timely estimation of R_e, applied to COVID-19 epidemic data from 170 countries. We thoroughly evaluate the method on simulated data, and present an intuitive web interface for interactive data exploration. We show that, in early 2020, in the majority of countries the estimated R_e dropped below 1 only after the introduction of major non-pharmaceutical interventions. For Europe the implementation of non-pharmaceutical interventions was broadly associated with reductions in the estimated R_e. Globally though, relaxing non-pharmaceutical interventions had more varied effects on subsequent R_e estimates. Our framework is useful to inform governments and the general public on the status of epidemics in their country, and is used as the official source of R_e estimates for SARS-CoV-2 in Switzerland. It further allows detailed comparison between countries and in relation to covariates such as implemented public health policies, mobility, behaviour, or weather data.

Keywords: SARS-CoV-2; epidemiology; global health; human; infectious disease; microbiology; public health surveillance; reproductive number; viruses.

Figure 1.. Evaluation of the pipeline on simulated data.

(A) The specified $R_e$ trajectory (black line) was used to stochastically simulate 100 trajectories of observed cases. From each trajectory, we estimated $R_e$ (yellow boxplots) and constructed a 95% confidence interval (purple boxplots of the lower/upper endpoint). (B) Fraction of simulations for which the true $R_e$ value was within the 95% confidence interval. The dashed red line indicates the nominal 95% coverage. (C) Root mean squared relative error for every time point. (D) Fraction of simulations for which $R_e$ is estimated to be significantly above or below one, depending on the true value of $R_e$.

Figure 2.. Stability of the Swiss Re estimates based on confirmed COVID-19 cases, upon adding additional days of observations.

(A) Line segments correspond to 3 weeks of estimates made with the same input data (e.g. data up to December 1st). The segments were assigned an arbitrary colour for ease of distinction. For each day, $R_e$ estimates and associated 95% confidence intervals (CIs) are overlaid, from the first possible estimate for that day up to estimates including 3 additional weeks of data. The latter, always the left end of a line segment, corresponds to the stable estimate. (B) Percentage of the first estimated CI that is contained in the stable CI based on data from 30 April 2021. This percentage was calculated as the width of the intersection of both CIs, divided by the width of the first CI. The colour indicates whether the stable $R_e$ estimate was contained in the first reported CI. In both rows, the left column shows uncertainty intervals from EpiEstim on the original data, and the right our improved 95% CIs. Both columns use the same pipeline, and differ only in the construction of the uncertainty intervals.

Figure 3.. Comparison of published Re estimates for three countries.

Point estimates are presented with a solid line and 95% confidence intervals are presented as coloured ribbons.

Figure 4.. Example panels from the online dashboard.

(A) Swiss case incidence with evidence of weekly testing patterns (top row), $R_e$ estimates with associated 95% confidence intervals from four types of observation data (middle row), and timeline of stringency index and vaccination coverage (bottom row). (B) World map of incidence per 100’000 inhabitants over the last 14 days. One can also display the worldwide $R_e$ estimates instead. (C) Comparison of $R_e$ estimates across four countries (Austria, Chile, India and Morocco), with timelines of stringency indices and vaccination coverage. All panels were extracted on May 12, 2022. Dashboard url: https://ibz-shiny.ethz.ch/covid-19-re-international.

Figure 5.. The association between the implementation or lifting of non-pharmaceutical interventions and changes in Re until May 2021.

(A) The change in the estimated $R_e$ at the same time as ($R\mathrm{}\mathrm{(}t\mathrm{)}$) or following ($R\mathrm{}\mathrm{(}t\mathrm{+}\mathrm{7}\mathrm{)}$ and $R\mathrm{}\mathrm{(}t\mathrm{+}\mathrm{14}\mathrm{)}$) the implementation (above x-axis) or lifting (below x-axis) of NPIs in a given week. (B) The change in the estimated $R_e$ related to the change in mobility in the same week. The error bars indicate the Q1 and Q3 quartiles.

Appendix 2—figure 1.. Schematic of the effect of smoothing on the ability to estimate when Re=1.

The true $R_e$ is indicated by the black solid line, the black dashed line shows a linear approximation of the smoothed $R_e$. Instead of crossing 1 at $t_{\mathrm{0}}$, this line crosses 1 at $t_{\mathrm{0}}\mathrm{+}s$.

Appendix 3—figure 1.. Performance of our method on simulated scenarios with differing slopes.

(A) The specified $R_e$ trajectory (black line; see Methods) was used to simulate a trajectory of reported cases (with Swiss case observation noise) 100 times. From each trajectory we estimated $R_e$ (yellow boxplots), and constructed a 95% confidence interval (purple boxplots of the lower/upper endpoint). We varied the time it took to change from one $R_e$ value to the next, $t\mathrm{\in }\mathrm{\{}\mathrm{7}\mathrm{,}\mathrm{14}\mathrm{,}\mathrm{28}\mathrm{\}}$ (columns). Larger values of $t$ correspond to less abrupt changes. (B) The fraction of simulations where the true $R_e$ value was within the 95% confidence interval. The dashed red line indicates the nominal 95% coverage. (C) The root mean squared relative error for every time point. (D) The fraction of simulations where we estimate $R_e$ is significantly above or below one, depending on the true value of $R_e$. We see that the method closely tracks the true $R_e$ in all scenarios, although the error is greater for steeper slopes. In the case of steeper changes in $R_e$ the overall size of the epidemic is also smaller, which explains the larger confidence intervals.

Appendix 3—figure 2.. Performance of our method, modified to skip the smoothing step in the pipeline, on simulated scenarios with observation noise.

The specified $R_e$ trajectory (black line; see Methods) was used to simulate a trajectory of reported cases (with varying country-specific noise profiles; rows) 100 times. From each trajectory we estimated $R_e$ (yellow boxplots), and constructed a 95% confidence interval (purple boxplots of the lower/upper endpoint). Contrary to our normal pipeline, the observations were not smoothed prior to the deconvolution and $R_e$ estimation. We see that the estimates are highly variable.

Appendix 3—figure 3.. Performance of our method on simulated scenarios with observation noise.

The columns differ in the method used to construct confidence intervals: EpiEstim reports the 95% HPD of EpiEstim on the original data, estimateR refers to our method, Quantiles and Reverse Quantiles use the 5 and 95% quantiles of the estimated $R_e$ to construct the CIs. (A) The specified $R_e$ trajectory (black line; see Methods) was used to simulate a trajectory of reported cases (with varying country-specific noise profiles; rows) 100 times. From each trajectory we estimated $R_e$ (yellow ribbons represent the estimated mean ± sd across 100 simulations), and constructed a 95% confidence interval (purple ribbons represent the mean ± sd of the estimated lower/upper endpoint). (B) The fraction of simulations where the true $R_e$ value was within the 95% confidence interval. The dashed red line indicates the nominal 95% coverage.

Appendix 3—figure 4.. Performance of our method on simulated scenarios with varying population size, using confidence intervals from EpiEstim.

(A) We specified a constant $R_e\mathrm{\in }\mathrm{\{}\mathrm{0.8}\mathrm{,}\mathrm{1}\mathrm{,}\mathrm{1.5}\mathrm{\}}$ value (black line; rows) to simulate a trajectory of reported cases (with Swiss case observation noise) 100 times. From each trajectory we estimated $R_e$ (yellow boxplots), and constructed a 95% confidence interval (purple boxplots of the lower/upper endpoint). The simulated scenarios had differing initial incidence of $I_{\mathrm{0}}\mathrm{\in }\mathrm{\{}\mathrm{10}\mathrm{,}\mathrm{100}\mathrm{,}\mathrm{1000}\mathrm{,}\mathrm{5000}\mathrm{,}\mathrm{10000}\mathrm{\}}$ infections per day (columns). In the top row, ${\displaystyle R_e<1}$ so the epidemic is decreasing. In the middle row, $R_e\mathrm{=}\mathrm{1}$, the epidemic is constant, and in the bottom row, ${\displaystyle R_e>1}$, the epidemic is increasing. The bias at the start is due to the initialisation of the simulation. (B) The fraction of simulations where the true $R_e$ value was within the 95% confidence interval. The dashed red line indicates the nominal 95% coverage. We see that the EpiEstim coverage strongly declines with increased epidemic size.

Appendix 3—figure 5.. Performance of our method on simulated scenarios with varying population size, using the Union of EpiEstim and Block bootstrap 95% confidence intervals.

(A) We specified a constant $R_e\mathrm{\in }\mathrm{\{}\mathrm{0.8}\mathrm{,}\mathrm{1}\mathrm{,}\mathrm{1.5}\mathrm{\}}$ value (black line; rows) to simulate a trajectory of reported cases (with Swiss case observation noise) 100 times. From each trajectory we estimated $R_e$ (yellow boxplots), and constructed a 95% confidence interval (purple boxplots of the lower/upper endpoint). The simulated scenarios had differing initial incidence of $I_{\mathrm{0}}\mathrm{\in }\mathrm{\{}\mathrm{10}\mathrm{,}\mathrm{100}\mathrm{,}\mathrm{1000}\mathrm{,}\mathrm{5000}\mathrm{,}\mathrm{10000}\mathrm{\}}$ infections per day (columns). In the top row, ${\displaystyle R_e<1}$ so the epidemic is decreasing. In the middle row, $R_e\mathrm{=}\mathrm{1}$, the epidemic is constant, and in the bottom row, ${\displaystyle R_e>1}$, the epidemic is increasing. The bias at the start is due to the initialisation of the simulation. (B) The fraction of simulations where the true $R_e$ value was within the 95% confidence interval. The dashed red line indicates the nominal 95% coverage. We see that for a wide range of infection incidences, our 95% confidence interval is informative and covers the true value of $R_e$.

Appendix 3—figure 6.. Performance of our method on simulated scenarios using a fixed shift versus the deconvolution to infer infection incidence.

The fixed shift method shifts the observations back by the mean of the delay distribution (here assumed to correspond to confirmed cases). (A) The specified $R_e$ trajectory (black line; see Methods) was used to simulate a trajectory of reported cases (with Swiss case observation noise) 100 times. From each trajectory we estimated $R_e$ (yellow boxplots), and constructed a 95% confidence interval (purple boxplots of the lower/upper endpoint). (B) The fraction of simulations where the true $R_e$ value was within the 95% confidence interval. The dashed red line indicates the nominal 95% coverage. The average coverage in this scenario was 0.90 with deconvolution and 0.78 with the fixed shift. (C) The root mean squared relative error for every time point. The average (cumulative) RMSE in this scenario was 0.0706 (25.4) with deconvolution and 0.0726 (26.6) with the fixed shift.

Appendix 3—figure 7.. Performance of our method on simulated scenarios with misspecified delay distributions.

When estimating $R_e$, we misspecified the mean of the delay distribution (5.5 for symptom-onset to case confirmation) by the numbers above the columns. (A) The specified $R_e$ trajectory (black line; see Methods) was used to simulate a trajectory of reported cases (with Swiss case observation noise) 100 times. From each trajectory we estimated $R_e$ (yellow boxplots), and a 95% confidence interval (purple boxplots of the lower/upper endpoint). (B) The fraction of simulations where the true $R_e$ value was within the 95% confidence interval. The dashed red line indicates the nominal 95% coverage. The average coverage in this scenario was 0.72, 0.85, 0.88, 0.90, 0.90, 0.89, 0.82, 0.69 from –5 to 10. (C) The root mean squared relative error for every time point. The average (cumulative) RMSE in this scenario was 0.0989 (36.0), 0.0807 (29.2), 0.0746 (26.9), 0.0706 (25.4), 0.0701 (25.2), 0.0726 (26.0), 0.0909 (32.3), 0.134 (47.0) from –5 to 10.

Appendix 3—figure 8.. Performance of our method on simulated scenarios with time-varying delay distributions.

The observations were simulated with a time-varying delay distribution for (A,B,C) confirmed cases, or (D,E,F) deaths (see Methods), and then estimated with (right column) or without (left column) taking the time-varying distributions into account. (A, D) The specified $R_e$ trajectory (black line; see Methods) was used to simulate a trajectory of reported cases or deaths (with Swiss case observation noise) 100 times. From each trajectory we estimated $R_e$ (yellow boxplots), and a 95% confidence interval (purple boxplots of the lower/upper endpoint). (B, E) The fraction of simulations where the true $R_e$ value was within the 95% confidence interval. The dashed red line indicates the nominal 95% coverage. For the cumulative cases, the average coverage in this scenario was 0.89 without and 0.89 with time variation. For the deaths, the average coverage was 0.83 without and 0.92 with time variation. (C, F) The root mean squared relative error for every time point. For the cumulative cases, the average (cumulative) RMSE was 0.0728 (26.2) without and 0.0783 (28.5) with time variation. For the deaths, the average (cumulative) RMSE was 0.113 (39.5) without and 0.0861 (30.8) with time variation.

Appendix 4—figure 1.. Stability of Re estimates at present.

We estimated $R_e$ through time repeatedly on 4 scenarios. With each new iteration, we added one new data point at present. Each $R_e$ trajectory is presented with its own colour. Purple trajectories are iterations for which the last known case data point was furthest in the past, yellow trajectories are trajectories for which the last known case data point was closest to the present. In each panel and for each $R_e$ trajectory, 100 simulation replicates were aggregated. The median of mean estimates are presented with lines and medians of upper and lower bounds of 95% confidence intervals are shown with translucent ribbons. For each scenario, two panels are presented. Each time the right panel correspond to raw estimates and the left panel corresponds to the same estimates with the last 4 $R_e$ values removed from each $R_e$ trajectory. (A) Stable $R_e$. (B) Gradual increase in $R_e$. (C) Gradual decrease in $R_e$. (D) Abrupt decrease in $R_e$.

Appendix 4—figure 2.. Mean delay in Switzerland between onset of symptoms and reporting.

For each date, the mean is taken over the last 300 reports with known symptom onset date, based on line list data from the FOPH. For early dates, before 300 reports are available, the mean is taken over the first 300 reports.

Appendix 4—figure 3.. Comparison of the Re estimates with or without accounting for known symptom onset dates and for time-variability on reporting delays.

The comparison is based on time series of confirmed cases in Switzerland, from line list data provided by the FOPH. Both the inclusion of known symptom onset dates and of the time-variability of reporting delay distributions have an effect on the Re estimates, in particular for early estimates in this case. The fraction of cases with known symptom onset date has drastically reduced since November 2020, hence the overlap in curves with and without symptom onset data for later dates. For each trajectory the point estimate is shown with a line, and the translucent ribbon indicates the 95% confidence interval.

Appendix 4—figure 4.. Effective reproductive number estimates with or without accounting for known imports.

The comparison is based on time series of confirmed cases in Switzerland, from line list data provided by the FOPH. The analysis ignoring imports is unbiased if the number of imports equals the number of exports. Since the analysis accounting for imports is not accounting for exports, the results are a lower limit for the effective reproductive number. Very few imported cases were reported since November 2020, hence the complete overlap in the curves after that date. For each trajectory the point estimate is shown with a line, and the translucent ribbon indicates the 95% confidence interval.

Appendix 5—figure 1.. Square root and log transformations to stabilise the variance of residuals.

Each row corresponds to the results of observations from Switzerland (CHE), China (CHN), France (FRA), New Zealand (NZL) and the United States (USA), respectively. The first and last two plots correspond to the result of square root transformation and log transformation, respectively.

Appendix 5—figure 2.. Autocorrelation function (ACF) and partial autocorrelation function (PACF) plots of the observations from five different countries.

In each row, the two plots are the ACF and PACF plots of the observations from Switzerland (CHE), China (CHN), France (FRA), New Zealand (NZL) and the United States (USA), respectively.

Appendix 5—figure 3.. Autocorrelation function (ACF) and partial autocorrelation function (PACF) plots of observations from Switzerland and the fitted ARIMA model.

The two plots on the upper row are the ACF and PACF plots of the observations from Switzerland. The two plots on the lower row are the ACF and PACF plots of the residuals of the fitted ARIMA model.

Appendix 5—figure 4.. Simulated observations with and without noise.

The upper row shows the real observations from Switzerland (twice the same). The other four rows show simulated observations, the left column shows simulations without the noise term (${\stackrel{\mathrm{~}}{D}}_t$ in Section 4.7), and the right the simulated observations with the noise term ($D_t$ in Section 4.7).

Appendix 5—figure 5.. Estimated Re with different smoothing windows.

For each trajectory the point estimate is shown with a yellow line, and the purple ribbon indicates the 95% confidence interval.

Appendix 5—figure 6.. Estimated Re with different block sizes.

The ribbons indicate the 95% confidence interval.