Real-time estimation of the effective reproduction number of SARS-CoV-2 in Aotearoa New Zealand

Introduction

An important measure in epidemiology is the reproduction number. The basic reproduction number, R₀, is the average number of people infected by a single infected individual, in a fully susceptible population and in the absence of any public health measures. For the highly transmissible Delta variant of SARS-CoV-2 (the virus that causes COVID-19), R₀ is estimated to be at least 6, compared to around 2.5–3 for the original strain of the virus. In reality, the proportion of the population who have immunity to the virus, either through infection or vaccination, changes during the course of an epidemic and transmission can be altered by introduction of public health measures, such as facemask use, contact tracing, quarantining infectious individuals, stay-at-home orders, and school or workplace closures. The effective reproduction number, R_t, is the average number of people infected by a single infected individual at a particular time t, accounting for this partial population immunity and changes in transmission levels due to public health measures.

Tracking estimates of the effective reproduction number in real time is useful for: (1) estimating how quickly an outbreak is likely to grow or shrink under current public health measures, to inform decisions on whether measures need to be escalated or scaled back, and (2) assessing the effectiveness of these measures for reducing transmission. Different models have been developed for estimating R_t from time-series data on the number of reported cases. Each model comes with its own set of assumptions and limitations, therefore it is useful for decision makers to have a range of estimates from different models to compare. One such model, implemented in the EpiNow2 package in R statistical software (Abbott et al., 2020a; Abbott et al., 2020b; R Core Team, 2021), follows current best practices (Gostic et al., 2020) to provide estimates of the time-varying instantaneous reproduction number, case growth rate, doubling time, and case numbers by date of infection. Their modelling approach is based on a method by Cori et al. (2013) and Thompson et al. (2019), extended to account for delays from infection to reporting.

From the start of the COVID-19 pandemic until the end of 2021, New Zealand pursued an elimination strategy (Baker et al., 2020). This involved a four-tier alert level system of restrictions (Unite Against COVID-19, 2022b), including nationwide and regional lockdown measures, which successfully eliminated community transmission of SARS-CoV-2 during outbreaks in March–May 2020, August–September 2020 and February–March 2021. On 17 August 2021, after several months with no known community transmission, New Zealand detected its first case of the Delta variant in the community, in a 58-year old man with no clear link to the border. At that time, only 38% of the population had received at least one vaccine dose, which would have been insufficient to prevent a large outbreak occurring without strong public health measures in place (Steyn et al., 2021; Steyn et al., 2022). The government therefore made the decision to move the entire country into Alert Level 4, the most stringent tier of restrictions (Table 1). Over the following days, contact tracing and surveillance testing efforts were ramped up in an attempt to ring-fence the outbreak as daily case numbers rapidly increased to 79 new cases on 26 August. Whole genome sequencing later identified possible links to infected travellers who arrived at a managed isolation and quarantine (MIQ) facility from New South Wales, Australia (Jelley et al., 2022).

Ten days after the introduction of Alert Level 4 restrictions, daily cases levelled off then steadily declined. Results from widespread community testing and wastewater testing suggested the outbreak was so far contained within Auckland and on 1 September regions outside of Auckland began a phased easing to Alert Level 2. Over the following weeks, daily case numbers dwindled in a long outbreak tail, with 10-20 new cases reported each day. However, elimination proved more difficult than previous outbreaks of earlier variants, particularly when transmission started to occur within unvaccinated groups in emergency and transitional housing, where public health measures including contact tracing and self-isolation were more challenging or impossible. From 22 September, restrictions in Auckland were eased to Alert Level 3 and, following detections in the neighbouring regions of Waikato and Northland, these regions moved up to Alert Level 3.

On 6 October, despite a rise in case numbers, the government announced further easing of some restrictions for Auckland and case numbers continued to steadily grow. This decision marked a transition from New Zealand’s elimination strategy to one of suppression: using restrictions to control case numbers to low levels while continuing to increase vaccination coverage in the population through an accelerated vaccine roll-out. Modelling suggested that well over 90% of New Zealand’s eligible population, aged 12 years and over, need to be fully vaccinated (two vaccine doses) to safely control community transmission without the need for strong public health measures, such as stay-at-home orders and workplace closures (Steyn, Plank & Hendy, 2021; Steyn et al., 2022). By mid-November, with close to 80% of the eligible population vaccinated with two doses and all regions except Auckland at Alert Level 2, the Delta outbreak reached its peak at approximately 200–220 new cases per day. Cases steadily declined from this time until January 2022 and by 18 January there were a total of 11,396 reported community cases and 26 deaths related to the Delta outbreak. The four-tier Alert Level system was replaced with a 3-tier COVID-19 Protection Framework (Unite Against COVID-19, 2022a) (traffic light levels, with red being the highest restriction setting, focused on minimising spread and providing protection through high vaccination coverage) which came into effect on 2 December 2021. Following the detection of nine community cases infected with the highly transmissible Omicron variant announced on 23 January, cases numbers rapidly increased again and New Zealand experienced its largest COVID-19 outbreak to date, peaking in mid-March 2022.

The ability to track changes in case numbers and the reproduction number in real time is critical to inform decisions about when and where additional measures should be introduced to prevent high case numbers overwhelming healthcare and contact tracing systems, and when these measures can be safely eased. In this work, we use EpiNow2 to estimate the time-varying effective reproduction number R_t in New Zealand during the August 2021–January 2022 outbreak of the Delta variant. The EpiForecasts team who developed EpiNow2, apply the model to global, publicly available case data, aggregated using the covidregionaldata R package (Abbott et al., 2020a; Abbott et al., 2020b; Abbott et al., 2020c; R Core Team, 2021), and publish real-time estimates of the effective reproduction number on their webpage (https://epiforecasts.io/covid/) (Abbott et al., 2020b). However, there is value in applying this method using more detailed line list data held by Ministry of Health, which allows model parameters to be tailored specifically to the New Zealand setting. We compare our estimates to these EpiForecasts estimates and summarise the differences in assumptions made in our implementation of the model. Throughout New Zealand’s Delta outbreak and in the current Omicron outbreak, we continue to update our estimates regularly as new data become available and provide these to Ministry of Health to inform operational planning and decision-making in real time.

Materials & Methods

Results

Using subsampled bootstraps of the delays from symptom onset to reporting, up to a maximum of 34 days, we fitted a lognormal distribution of onset-to-report delays with a mean of 3.84 days and standard deviation of 3.94 days (parameters: μ = 0.988 (SD = 0.088) and σ = 0.847 (SD = 0.063), which was a good visual match to the data (Fig. S1). Fitting the EpiNow2 model to the 22-week time-series of daily reported case numbers, we obtained estimates for the time-varying mean reported cases, incidence, effective reproduction number and growth rate of New Zealand’s Delta outbreak up to 18 January 2022. The Rhat statistic was less than 1.05 for all parameter estimates indicating model convergence. The estimated daily number of new cases (ribbons) are shown alongside the actual reported case numbers (grey bars) in Fig. 1. Mean reported case counts (D_t ω_d,t) by date of report (Fig. 1A), were a good visual match to the data, first peaking at a median 57 (90% CrI 32, 90) cases on 27 August 2021 (cf. 79 on 26 August in actual reported cases). Plotting cases by their estimated date of infection (I_t) shows the first peak in incidence occurred on 20 August with a median 73 (48, 101) new cases infected that day (Fig. 1B).

Figure 2 shows the change in estimated effective reproduction number (R_t) over time alongside key timings of changes to levels of restrictions. On 17 August, R_t is greater than one but immediately starts to decline after introduction of Alert Level 4 restrictions from 18 August. This, along with the decline in incidence from 20 August, suggests that these restrictions had a near-immediate effect on reducing transmission, though this reduction did not filter through to reported case numbers for another 7–10 days. After twelve days under Alert Level 4 restrictions, the median R_t was reduced to 0.57 (0.42, 0.76), growth rate r_t was −0.13 (−0.18, −0.07) per day and doubling time was −5.48 (−3.89, −10.26) days (Figs. 2–3). Estimated daily cases by infection date were reduced to a minimum of 14 (11, 18) per day by 6 September (Fig. 1B) and the outbreak entered its long tail. After this time, the reproduction number and daily infections steadily increased, with R_t growing to values greater than one from 11 September. Coinciding with the further easing of some restrictions in Auckland during September and October, median R_t remained over 1 and New Zealand experienced its largest outbreak since the start of the pandemic. Transmission was at its highest on 2 October, with median R_t = 1.27 (1.15, 1.38), corresponding to a growth rate of 0.07 (0.04, 0.10) per day and doubling time of 9.66 (16.68, 6.86) days. However, by mid-November, regional Alert Level 3 restrictions remained in place and vaccination coverage had rapidly increased from 38% of the eligible population having received at least one vaccine dose on 17 August to 89% by 10 November. The median R_t subsequently decreased below 1, causing the outbreak to peak at an estimated 198 (172, 229) new cases infected on 10 November (cf. 226 cases reported on 10 November in actual data). From this point, median R_t remained relatively constant at around 0.8 to 0.9 as cases continued to decline and the COVID-19 Protection Framework was introduced. By January 2022, prior to detection of community cases with the Omicron variant, the outbreak had decreased to a median 20–30 cases infected per day.

Note that during the first 6.5 weeks of the outbreak, the vast majority of reported cases were detected only in the Auckland region. On 3 October, cases were detected in Waikato and from this time on transmission occurred within both regions but likely not between the two regions due to the Auckland travel boundary in place until 14 December. This is supported by whole genome sequencing which suggested that the Waikato cluster arose from a single introduction from Auckland (Hadfield et al., 2018; Khare et al., 2021), while clusters detected in Northland in late October—early November likely represented multiple introductions (New Zealand Ministry of Health, 2021).

Refitting the model using the mean generation times of Hart et al. (2022) resulted in similar estimates of R_t as shown in Fig. S3. With the shorter mean generation time of 3.2 days, R_t reached its minimum at a median 0.73 (0.53, 0.83) on 30 August and its maximum of 1.24 (1.15, 1.32) on 1 October. Assuming a longer mean generation time of 4.6 days resulted in a minimum R_t of 0.52 (0.40, 0.72) on 30 August and maximum on 1.32 (1.18, 1.46) on 5 October. Summary statistics for all posterior estimates are provided in Table S1.

Estimates of daily numbers of reported cases by report date and infection date, and the effective reproduction number published by EpiForecasts are shown in Fig. S2, for two 16-week periods up to 15 October 2021 and up to 28 January 2022. The EpiForecasts estimates were broadly similar to our estimates but tended to over-estimate the numbers of new infections each day, particularly in periods of low local transmission. Epiforecasts fit to data published by the World Health Organisation (World Health Organisation, 2020), which sources data from the New Zealand Ministry of Health (New Zealand Ministry of Health, 2022), but does not include information on the importation status of New Zealand’s cases. Case numbers modelled by EpiForecasts were therefore biased upwards by inclusion of internationally-imported cases detected in MIQ facilities. As the EpiForecasts model attributes all new imported cases to local transmission, this could result in over- or under-estimation of R_t. We were able to isolate the domestically-acquired cases for model fitting to provide less biased estimates of community transmission occurring in New Zealand.

Discussion

Our study demonstrates how real-time estimation of the effective reproduction number can effectively support a timely public health response during a community outbreak. New Zealand’s elimination strategy and use of stringent stay-at-home measures were successful at controlling several outbreaks of earlier variants even while vaccination coverage remained low. However, they were less effective at controlling the highly transmissible Delta variant and an accelerated vaccination rollout in combination with ongoing regional restrictions were ultimately required to reduce R_t below one. For community outbreaks where a significant proportion of cases are internationally-imported, models should distinguish between imported and domestic infections to avoid biasing estimates of transmission.

Delays from infection to reporting mean that cases reported on a particular day are really indicative of the transmissibility approximately 7–10 days prior. The EpiNow2 framework accounts for these delays, allowing case counts to be mapped from their date of reporting to date of infection to estimate an instantaneous measure of transmissibility under current conditions. While there can be larger uncertainty around R_t estimates closer to ‘now’ (due to these estimates being based on less data), the EpiNow2 approach is particularly well-suited for informing real-time decisions and facilitating a rapid public health response to changes in transmission. It also allows observed changes in R_t to be related to changes in interventions or other variables affecting transmission at certain times. More naive approaches that estimate the reproduction number directly from reported cases without accounting (or only partially accounting) for these delays provide a lagged view of an outbreak’s trajectory and may obscure the temporal relationship between changes in policy or behaviour and transmission dynamics (Gostic et al., 2020).

The instantaneous reproduction number R_t estimated by EpiNow2 is a property of the epidemic at a particular time t. It measures the average number of secondary infections per infectious individual at time t, assuming that the current conditions remain unchanged until those individuals who were infectious at time t have recovered. This differs from another widely computed measure, the case reproduction number R_c, which is the average number of people an individual infected at time t eventually infects, with no assumption of conditions being unchanged after time t. R_c is the measure estimated by Wallinga and Teunis-type methods, though Wallinga & Teunis (2004) consider individuals with symptom onset at time t as opposed to infection at t (Cori et al., 2013). In a situation where a change in conditions (such as the introduction of stay-at-home measures) brings about an abrupt change in transmissibility, the instantaneous reproduction number R_t will also change abruptly at that time. In contrast, R_c will transition smoothly because it estimates the number of people each case will actually infect, accounting for the fact that the infectious periods of cases infected prior to the change in conditions may span times before and after this change in transmission. Estimates of R_c rely on information about future incidence so can only be obtained in retrospect, and suffer from right censoring of reported case data in an epidemic that is still in progress. For estimates of the instantaneous R_t, which are calculated from current and past incidence, right censoring causes fewer issues, especially where the delays to reporting are accounted for; this measure is therefore more appropriate for tracking transmissibility in real-time (Fraser, 2007).

There are several important limitations of the EpiNow approach. First, the model assumes that the proportion of infections that are detected remains constant over the study period, i.e., consistent methods and effort for surveillance (e.g., testing, contact tracing) and following consistent case definitions. This is assumption is likely violated; in reality, the level of under-reporting varies over time as contact tracing efforts are ramped up and surveillance testing rates vary over time and over different regions. This was certainly the case over the first few weeks of the outbreak. However, such changes will only bias estimates temporarily if surveillance becomes consistent again after a change. The period of time required before estimates become unbiased depends on the generation time distribution, the delay distributions, and the correlation length scale of the Gaussian process used to model temporal variation in R_t.

Model estimates are sensitive to the choice of distribution for the delay from symptom onset to reporting. This distribution affects the upscaling of cases by date of onset that accounts for cases that have onset but have yet to be reported. We used a distribution estimated by fitting to data from New Zealand’s Delta outbreak. Reporting delays can vary over the course of an outbreak or between regions, especially early on in an outbreak or during periods of high transmission where the testing and contact tracing system may come under significant strain. The model assumes this distribution is static over time, which was a reasonable assumption for this outbreak (Table S2) but may not hold for future outbreaks in New Zealand or in other countries. If the true delay from onset to reporting at a given time is shorter than the estimated static delay distribution, then the model will overestimate onset case numbers, and vice versa for true delays longer than the distribution used. The bootstrapped subsampling approach somewhat mitigates these issues by allowing multiple delay distributions based on the observed data to be considered at the cost of increasing uncertainty in the estimates. However, if significant temporal or spatial variation in reporting delays is suspected in future outbreaks then using other estimation procedures that are robust to such changes (see e.g., Li & White, 2021; Günther et al., 2021) may be more appropriate.

Using distributions of the incubation period and onset-to-reporting delay to infer incidence from observed reported cases tends to result in over-smoothing of the incidence. The assumed distributions for generation time and incubation period sourced from Ganyani et al. (2020) and Lauer et al. (2020) were estimated by fitting to data on infections with the original strain of the virus in other countries and these distributions may differ for the Delta variant. However, the mean generation time used here sits between two estimates for Delta reported by Hart et al. (2022) and re-ftting the model using these alternative mean generation times had little impact on our estimates of R_t. Furthermore, the model incorporates uncertainty in the means and standard deviations of both distributions (as well as in the reporting delay distribution), which goes some way to accounting for variation from the expected distributions, for example due to different populations or variants. For analyses of future outbreaks, parameters for the incubation period, generation time and reporting delay distributions should continue to be calibrated using up-to-date data as it becomes available. In particular, recent studies have suggested the mean generation time for the Omicron variant may be shorter than the Delta variant (Abbott et al., 2022).

As with most models, EpiNow2 estimates are more reliable when daily case numbers are large and the time series contains at least 14 days of non-zero cases. When analysing small outbreaks, slight changes in daily case numbers can lead to large variations in model estimates. Such estimates should be treated with caution and interpreted alongside other information sources, for example the detailed case-specific information collected by contact tracers. Case numbers were too low to robustly estimate regional trends during New Zealand’s Delta outbreak using this technique. We were able to obtain reliable estimates at national scale, though the estimates based on partial data (orange bands in Figures) in the last 9 days of the time-series are relatively unstable and sensitive to the number of reported cases on the final time point. These values may change significantly as data are updated and new time points are added to the series.

We do not explicitly model interventions or vaccination so could not establish their relative contributions for reducing transmission. Other models attempting to infer the relative effectiveness of interventions and vaccination (see e.g., Li et al., 2022) generally require data from multiple countries that have implemented similar measures and with varying levels of vaccination coverage. As well as changes to levels of restrictions and vaccination, other factors can also make it challenging to attribute a change in R_t to a particular public health measure. For example, variability in the level of compliance with public health measures over time (e.g., willingness to wear facemasks, adherence to restrictions on between-household contact) and between different subgroups in the population influences R_t. Transmission levels will also vary between subgroups of the population, due to different demographics, levels of vaccination, and different environmental conditions, such as average household size and density of housing. The model masks this heterogeneity so the estimates of the time-varying R_t are an attribute of the outbreak, in whichever subgroups it is occurring at that particular time. This has important implications for small outbreaks, where the subgroups in which the virus is spreading may not be representative of the wider population. If case numbers are high enough for reliable inference, such fragmented outbreaks may be better modelled as separate local outbreaks, for example regional-level models for different district health boards (Vegvari et al., 2021). In New Zealand’s Delta outbreak, different sub-clusters have been identified (Jelley et al., 2022) but there were insufficient case numbers to fit the model for each sub-cluster separately. It is therefore not possible to disentangle the relative contributions that public health measures and compliance levels had on transmission, from other attributes of the sub-clusters dominating the outbreak at that particular point in time. For example, the increase in transmission towards the end of September could have been driven by the easing of restrictions in Auckland, changing levels of public compliance, the virus spreading into under-vaccinated groups or populations staying in emergency and transitional housing, or the combination of these factors.

Conclusions

Our analysis of New Zealand case data demonstrates the effectiveness of high vaccination coverage in combination with public health measures for controlling an outbreak of the Delta variant of SARS-CoV-2. Modelling approaches that account for delays from infection to reporting reveal important changes in underlying transmission dynamics that may not be apparent when dealing only with reported cases. Real-time estimation of the effective reproduction number remains critical for informing decision-making and operational planning in New Zealand’s ongoing pandemic response.

Supplemental Information