Estimating the time-varying reproduction number of SARS-CoV-2 using national and subnational case counts [version 2; peer review: 1 approved with reservations]

Background: Assessing temporal variations in transmission in different countries is essential for monitoring the epidemic, evaluating the effectiveness of public health interventions and estimating the impact of changes in policy. Methods: We use case and death notification data to generate daily estimates of the time-varying reproduction number globally, regionally, nationally, and subnationally over a 12-week rolling window. Our modelling framework, based on open source tooling, accounts for uncertainty in reporting delays, so that the reproduction number is estimated based on underlying latent infections. Results: Estimates of the reproduction number, trajectories of infections, and forecasts are displayed on a dedicated website as both maps and time series, and made available subnationally. This allows public health officials and policymakers to track the progress of the outbreak in near real-time using an epidemiologically valid measure. As well as providing regular updates on our website, we also provide an open source tool-set so that our approach can be used directly by researchers and policymakers on confidential data-sets. We hope that our tool will be used to support decisions in countries worldwide throughout the ongoing COVID-19 pandemic. Authors develop a method for monitoring the epidemic, evalating the effectiveness of policy intervention and estimating the impact of policy changes in public health, during the COVID-19 pandemic. They could improve the background of their paper. The authors should state clearly the aim of the paper at the end of the introduction section. Some statements at the end of the introduction could be better used for the method section. state from which countries they have drawn data and declare the


Introduction
The coronavirus disease 2019  pandemic that emerged in December 2019 has since spread to over 100 countries in every continent except Antarctica. While some information on the progress of an outbreak in a given country can be gained from the reported numbers of confirmed cases and deaths, these numbers can obscure changes in the underlying dynamics of the outbreak due to delays between infection and the eventual reporting of a case or death. Accounting for the uncertain delays from infection to symptom onset, and the uncertain delays from symptom onset to hospital admission, diagnostic testing or potential death, followed by further delays until data are recorded in official statistics, requires the use of specific statistical methods for handling right-truncated data [1][2][3] , uncertainty, and the creation of a "nowcast" 4,5 (an estimate of the current number of newly infected or symptomatic cases).
A method for tracking the progress of an outbreak is to measure changes in the time-varying reproduction number (effective reproduction number), which represents the average number of secondary infections generated by each new infectious case [6][7][8] . This approach can be advantageous compared to monitoring numbers of newly reported or symptomatic cases since, in principle, reproduction number estimates reflect variations in transmission intensity. Due to the delays in disease progression, recorded numbers of newly notified or symptomatic cases will increase or decrease for a period after transmissibility has reduced or increased, respectively. Monitoring changes in the time-varying reproduction can account for this delay and reveals variations in transmissibility that are not clear when using only reported cases. This paper outlines the methods used to produce the website (https://epiforecasts.io/covid/), and data resource 9 , we have developed that presents real-time estimates and forecasts of reported cases by date of infection and the respective timevarying reproduction numbers globally, regionally, nationally and subnationally for Covid-19. This website relies on methods implemented in the EpiNow2 R package and data aggregated in the covidregionaldata R package, both developed by the authors 10,11 . Our estimates overcome some of the limitations of naive implementations that derive estimates for the reproduction number directly from numbers of reported cases without adjusting (or with only partial adjustments) for the delay from infection to symptom onset or from onset to notification. Our approach also incorporates multiple sources of uncertainty that if excluded can bias estimates. The code that creates and updates the website is open source, and documented for use by others, allowing policymakers and researchers to run analyses using confidential data. The methods outlined in this paper and corresponding code base are under development, and new versions of this live article will be released alongside changes to the methods to create a record of the methodology used throughout the pandemic.

Data
We use daily counts of confirmed cases and deaths reported by the European Centre for Disease Control from the last 12 weeks for all analyses conducted at the national level 11,12 . To estimate the delay from symptom onset to reporting (once confirmed with a positive laboratory test), we use all cases from a publicly available linelist for which onset and notification dates are available 11,13 . This linelist combines all known linelist data from over 100 countries at the time of writing. Countries are only included in the reported estimates if within the last 12 weeks they have fewer than 14 days with non-zero case counts. This restriction reduces the likelihood of spurious estimates for countries with limited transmission or case ascertainment.
For sub-national analyses, the data is aggregated using the covidregionaldata R package developed by the authors. Individual data sources are reported on the respective pages of our website. The data are fetched from government departments or from individuals who maintain a data source if no official data are available. Similarly to national estimates, subnational areas are only included if they report at least 14 days with non-zero cases in the last 12 weeks.
All analyses described below are run daily for each national or subnational entity under consideration. An automated timestamp is used to evaluate if data has been updated since the last time estimates were made in order to avoid repeatedly estimating based on the same data.

Delays between case onset and report
To estimate the reporting delay (i.e the delay between onset and case report or death) with appropriate uncertainty, we fit a log-normal distribution, using use the statistical modelling program stan 10,14 , to 100 subsampled bootstraps (each with 250 samples drawn with replacement) of the available delay data. Accounting for left and right censoring occurring in the data as each date is rounded to the nearest day and truncated to the maximum observed delay. There was insufficient data available on the various reporting delays to estimate spatially-or temporally-varying delays whilst also accounting for the biases induced by the growth rate of reported cases,

Amendments from Version 1
In this update we include details of our new open-source time-varying reproduction method that is based on inferring latent infections rather than attempting to reconstruct them via backsampling as discussed in the previous version of this article. This approach reduces bias in estimates and increases the potential for rapid changes over time. We also discuss, and provide links to, our flexible scheduling framework, which is developed in partnership with the Met office, and our dataverse, where version controlled reproduction number estimates can be found.
Any further responses from the reviewers can be found at the end of the article UPDATE so they were considered to be static over the 12 weeks of data considered each day.
This results in an onset to case report delay distribution with a mean of 6.5 days and a standard deviation of 17 days and an onset to death report delay distribution with a mean of 13.1 days and a standard deviation of 11.7 days. For computational reasons the maximum allowed delay is set to be 30 days. Dataset specific estimates are detailed on the respective country pages. Estimated delays are routinely updated as new data becomes available.
As data may also be right truncated due to unrecorded delays (i.e the delay between a case report and its appearance in an aggregated data set) we truncate all time-series to exclude the last 3 days of data, based on qualitative inspection of the stability of case counts in the datasets used.

Estimating the time-varying reproduction number and nowcasting reported infections
We estimated the instantaneous reproduction number (R t ) using the EpiNow2 R package (version 1.2.1) 10 on the last 12 weeks of available data, discarding estimates from the first 14 days globally, for United Nation regions, nationally, and subnationally for 10 countries. The instantaneous reproduction number represents the number of secondary cases arising from an individual showing symptoms at a particular time, assuming that conditions remain identical after that time, and is therefore a measure of the instantaneous transmissibility (in contrast to the case reproduction number -see Fraser (2007) 8 for a full discussion). EpiNow2 implements a Bayesian latent variable approach using the probabilistic programming language Stan 14 , which works as follows. The initial number of infections were estimated as a free parameter with a prior based on the initial number of cases, or deaths, respectively. For each subsequent time step, previous imputed infections (I t-1 ) were summed, weighted by an uncertain generation time probability mass function (w), and combined with an estimate of R t to give the incidence at time t (I t ) 6,7,10 . We used a log normal prior for the reproduction number (R 0 ) with mean 1 and standard deviation 1 reflecting our current belief that R t is likely to be centered around 1 in most of the world, with public health interventions and individual behaviour combining to prevent it from growing significantly larger for sustained periods. This contrasts with our earlier approach which was to use a gamma prior with a of mean 2.6 and standard deviation 2. This was based on early estimates for the basic reproduction number from the initial stages of the outbreak in Wuhan 15,16 with long tails to allow for differences in the reproduction number between countries. The infection trajectories were then mapped to mean reported case counts (D t ) by convolving over an uncertain incubation period and report delay distribution (convolved into ξ). Observed reported case counts (C t ) were then assumed to be generated from a negative binomial observation model with overdispersion ϕ (using an exponential prior with mean 1) and mean D t , multiplied by a day of the week effect with an independent parameter for each day of the week (ω (tmod7) ). Temporal variation was controlled using an approximate Gaussian process 17 with a squared exponential kernel (GP). In mathematical notation, The parameters of the Gaussian process kernnl were estimated during model fitting with the following priors. The length scale was given an inverse gamma prior with shape and scale values optimised to give a distribution with 98% of the density between 2 days and 21 days. The prior on the magnitude was standard normal. Each timeseries was fit independently using Markov-chain Monte Carlo (MCMC). A minimum of 4 chains were used with a warmup of 500 each and 4000 samples post warmup. Convergence was assessed using the R hat diagnostic 14 .
We used an estimate of the generation time sourced from 18 but refit using a log-normal incubation period with a mean of 5.2 days (SD 1.1) and SD of 1.52 days (SD 1.1) 19 rather than the incubation period used in the original study (code available here: https://github.com/seabbs/COVID19). This resulted in a distributed generation time with mean 3.6 days (standard deviation (SD) 0.7), and SD of 3.1 days (SD 0.8) for all estimates. The incubation period estimate was also used to convolve from unobserved infections to unobserved onsets in the model. See 10 for further details on the approach.
Estimating the daily growth rate and doubling time We estimated the rate of spread (r) by converting our R t estimates using an approximation derived in 20. The doubling time was then estimated by calculating ln(2) 1 r for each estimate of the rate of spread.

Estimated change in daily cases
We defined the estimated change in daily cases to correspond to the proportion of reproduction number estimates for the current day that are below 1 (the value at which an outbreak is in decline). It was assumed that if less than 5% of samples were subcritical then an increase in cases was definite, if less than 20% of samples were subcritical then an increase in cases was likely, if more than 80% of samples were subcritical then a decrease in cases was likely and if more than 95% of samples were subcritical then a decrease in cases was definite. For countries/regions with between 20% and 80% of samples being subcritical we could not make a statement about the likely change in cases (defined as unsure).

The effect of changes in testing procedure
The results presented here are sensitive to changes in COVID-19 testing practices and the level of effort put into detecting COVID-19 cases, e.g. through contact tracing. For example, if numbers of incident infections remain constant but a country begins to find and report a higher proportion of cases, then an increasing value of the reproduction number will be inferred. This is because all changes in the number of cases are attributed to changes in the number of infections resulting from previously reported cases and are not assumed to be a result of improved testing and surveillance. On the other hand, if a country reports a lower proportion of cases because a lower number of tests are performed (which can happen if reagents required for testing are no longer available, for example) or the surveillance system captures a lower proportion of infections, then the model will attribute this to a drop in the reproduction number that may not be a true reduction. In order for our estimates to be unbiased not all cases have to be reported, but the level of testing effort (and therefore the proportion of detected cases) must be constant 21 . This means that, whilst a change in testing effort will initially introduce bias, this will be reduced over time as long as the testing effort remains consistent from this point onwards.
Countries may also change the focus of their surveillance over the course of the outbreak. They may initially focus on identifying travellers returning from areas of known COVID-19 transmission and performing contact tracing on the contacts of known cases. As the outbreak evolves this may change to passive surveillance at hospitals. Here, the case definition may also change from tests based on polymerase chain reaction (PCR) to diagnoses based on symptoms and computed tomography (CT) scans. In the future, different kinds of COVID-19 tests may be deployed that could influence results, such as tests that detect both active and past infections.

Forecasting the reproduction number and case counts by date of infection
We forecast the time-varying effective reproduction number over a 14-day time horizon by assuming it remains the same as the last estimated R t . The reproduction number forecast is then transformed into a case forecast using the EpiNow2 model outlined in the previous section 10 . These forecasts are indicative only and should not be considered with a weight equal to the real-time estimates. Changes in contact rates, mobility, and public health interventions are not accounted for which may lead to significant inaccuracy.

Reporting
We report the median and 90% credible intervals for all measures with 20%, 50% and 90% credible intervals shown in figures. The analysis was conducted independently for all regions and is updated daily as new data becomes available. To highlight the proportion of cases that have yet to be reported (due to correcting for right truncation), we show a cut-off in figures based on the mean of all delays. Values prior to this point are defined as estimates with values past this point being defined as estimates based on partial data. In reality, this is a continuum with estimates closer to now progressively being based on less data and therefore becoming increasing uncertain. All estimates are available as downloadable files in csv format under an open-source license for use elsewhere 9 . The scheduling framework used to update our estimates is also available under an open-source license (https://github.com/ epiforecasts/covid-rt-estimate).

Website, summarised estimates, and interactivity
We use Rmarkdown templates and the distill framework to generate webpages summarising these estimates 22,23 . The RtD3 package is used to provide interactive visualisations of all estimates 24 . Estimates by country are provided on a dedicated static page along with global, and regional, summaries. More detailed subnational estimates are available for over 10 countries in an flexible framework into which additional subnational estimates will be added as more data becomes available.

Discussion
We provide a centralised resource which generates comparable daily estimates of the time-varying reproduction number and a daily nowcast of the number of cases newly infected derived using a standardised method. We account for the delay between infection and case notification and include all sources of quantifiable uncertainty. This resource may be useful for policymakers to track the progression of the COVID-19 outbreak and evaluate the effectiveness of intervention measures. As new data become available, we will include subnational estimates for additional countries, and provide additional support for public health agencies or researchers interested in applying our methods to their data.
There are several advantages associated with our approach. Firstly, reported counts are the only data required, which allows our approach to be used in a wide variety of contexts. It can be applied separately to counts of cases, hospital admissions, deaths or other metrics as long as appropriate delay distributions are used 25 . As our methodology is applied across a range of geographies our estimates can be compared without having to consider differences in the underlying approach (even if differences in testing should still be accounted for as discussed below). Finally, we have constructed our approach using open source tools and all of our code, raw data, and results are available online and developed with other users in mind. This means our methods can be readily applied by others to nonpublic data and be fully evaluated by end users.
Our approach is also subject to several limitations. Firstly, the model requires that the proportion of infections that are notified is constant over the 12 weeks considered. In other words, it requires consistency in the focus of the surveillance method, level of effort spent on testing, and case definition. Yet it is often the case that the level of under-reporting in a country changes over the course of an outbreak 21 . However, it should be noted that any changes in surveillance testing procedures will only bias the estimates temporarily if they begin to remain consistent again after they have changed. How long the bias remains in the reproduction number estimates will depend on the serial generation time and delay distributions, as well as the length scale of the Gaussian process used in the reproduction number estimation process. The impact of testing and other reporting biases vary between measures of transmission (test positive cases, hospital admissions, test positive deaths) 25 . For this reason we include estimates based on reported deaths and provide tooling to allow estimates to be produced for alternative datasets. In theory, estimates from disparate sources should be comparable using our approach, however if they in fact represent different sub-populations then there may be variation between them that can potentially be usefully interpreted.
In addition, the model is limited by how representative the delay that we use from infection to notification distribution is for a given location. As there is limited data to assess this, we estimate a bootstrapped global delay distribution using the combined data from every country. In particular, the delay from onset to notification can especially impact the upscaling of cases by date of onset that accounts for cases that have onset but not yet been reported. If the true delay from onset to notification for a given country is shorter than our global delay, then we will overestimate onset case numbers, and vice versa for true delays longer than the distribution we used. Additionally, estimates of the reporting delay distribution are known to be biased early in an epidemic and may vary over time 26 . However, our use of a bootstrapped subsampling approach mitigates these issues by allowing multiple delay distributions based on the observed data to be considered at the cost of increasing uncertainty in our estimates.
Our model is also limited by the data available to us. For example, the publicly available linelists contain little data on the importation status of cases. This means that cases counts may be biased upwards by attributing imported cases to local transmission. This bias is particularly problematic when case counts are low. Unfortunately, in the absence of data, this issue can only be explored via scenario analysis.
As more data becomes available, future work should look to refine the distributions used for generation time, incubation period, and the report delay. There is also the potential to extend the present model to account for changes in the delay from onset to notification over the course of an outbreak though additional data would need to be available for this to be possible. Finally, there is scope to explore how outbreak dynamics that differ among particular sub-populations, such as high-risk COVID-19 patients, can bias overall reproduction number estimates. This may be achieved by comparing reproduction number estimates from disparate data sources such as test positive cases, hospital admissions, and test positive deaths.
Our approach, providing real-time estimates of the reproduction number, serves as a valuable tool for decision makers looking to track the course of COVID-19 outbreaks. The nowcasts explicitly account for delays, using the same methodology across all countries and sub-national regions. These reproduction number estimates may also be used to ascertain the likely outbreak trajectory if no policy interventions are made. They can also provide real-time feedback on whether transmission is decreasing following a particular intervention, or whether it is increasing following the relaxing or lifting of current intervention measures. We hope that our website and the related toolkit will provide a valuable resource for devising strategies to contain COVID-19 outbreaks worldwide.

Francesco Chirico
Post-graduate School of Occupational Health, Università Cattolica del Sacro Cuore, Roma, Italy Authors develop a method for monitoring the epidemic, evalating the effectiveness of policy intervention and estimating the impact of policy changes in public health, during the COVID-19 pandemic. They could improve the background of their paper. The authors should state clearly the aim of the paper at the end of the introduction section. Some statements at the end of the introduction could be better used for the method section. With regard to pitfalls and alternatives in the use of case fatality ratio, see the paper entitled: In Methods section, authors should state from which countries they have drawn data and declare the study design.