Materials and methods

Estimating the time-varying effective reproduction number

Overview
The method used to estimate R𝑒𝑓𝑓 is described in Cori et al., 2013, as implemented in the R package, EpiNow (Abbott et al., 2020). This method is currently in development by the Centre for the Mathematical Modelling of Infectious Diseases at the London School of Hygiene and Tropical Medicine (London School of Hygiene & Tropical Medicine Mathematical Modelling of Infectious Diseases nCoV working group, 2020). Full details of their statistical analysis and code base is available via their website (https://epiforecasts.io/covid/).
The uncertainty in the R𝑒𝑓𝑓 estimates (shown in Figure 2; Figure 2—figure supplements 1, 2 and 3) represents variability in a population-level average as a result of imperfect data, rather than individual-level heterogeneity in transmission (i.e., the variation in the number of secondary cases generated by each case). This is akin to the variation represented by a confidence interval (i.e., variation in the estimate resulting from a finite sample), rather than a prediction interval (i.e., variation in individual observations).
We provide a brief overview of the method and sources of imperfect data below, focusing on how the analysis was adapted to the Australian context.

Data
We used line-lists of reported cases for each Australian state/territory extracted from the national COVID-19 database. The line-lists contain the date when the individual first exhibited symptoms, date when the case notification was received by the jurisdictional health department and where the infection was acquired (i.e., overseas or locally).

Reporting delays and under-reporting
A pre-hoc statistical analysis was conducted in order to estimate a distribution of the reporting delays from the line-lists of cases, using the code base provided by London School of Hygiene & Tropical Medicine Mathematical Modelling of Infectious Diseases nCoV working group, 2020. The estimated reporting delay is assumed to remain constant over time. These reporting delays are used to: (i) infer the time of symptom onset for those without this information, and; (ii) infer how many cases in recent days are yet to be recorded. Adjusting for reporting delays is critical for inferring when a drop in observed cases reflects a true drop in cases.
Trends identified using this approach are robust to under-reporting, assuming that it is constant. However, absolute values of R𝑒𝑓𝑓 may be biased by reporting rates. Pronounced changes in reporting rates may also impact the trends identified.
The delay from symptom onset to reporting is likely to decrease over the course of the epidemic, due to improved surveillance and reporting. We used a delay distribution estimated from observed reporting delays from the analysis period, which is therefore likely to underestimate reporting delays early in the epidemic, and overestimate them as the epidemic progressed. Underestimating the delay would result in an overestimate of R𝑒𝑓𝑓, as the inferred onset dates (for those that were unknown) and adjustment for right-truncation, would result in more concentrated inferred daily cases (i.e., the inferred cases would be more clustered in time than in reality). The converse would be true when overestimating the delay. The impact of this misspecified distribution will be greatest on the most recent estimates of R𝑒𝑓𝑓, where inference for both right-truncation and missing symptom onset dates is required.

Estimating the effective reproduction number over time
Briefly, the R𝑒𝑓𝑓 was estimated for each day from 24 February 2020 up to 5 April 2020 using line list data – date of symptom onset, date of report, and import status – for each state. The method assumes that the serial interval (i.e., time between symptom onset for an index and secondary case) is uncertain, with a mean of 4.7 days (95% CrI: 3.7, 6.0) and a standard deviation of 2.9 days (95% CrI: 1.9, 4.9), as estimated from early outbreak data in Wuhan, China (Nishiura et al., 2020). Combining the incidence over time with the uncertain distribution of serial intervals allows us to estimate R𝑒𝑓𝑓 over time.
A different choice of serial interval distribution would affect the estimated time varying R𝑒𝑓𝑓. This sensitivity is explored in detail in Flaxman et al., 2020, though we provide a brief description of the impact here. For the same daily case data, a longer average serial interval would correspond to an increased estimate of R𝑒𝑓𝑓 when R𝑒𝑓𝑓>1, and a decreased estimate when R𝑒𝑓𝑓<1. This effect can be understood intuitively by considering the epidemic dynamics in these two situations. When R𝑒𝑓𝑓>1 , daily case counts are increasing on average. The weighted average case counts (weighted by the serial interval distribution), decrease as the mean of the serial interval increases (i.e., as the support is shifted to older/lower daily case data). In order to generate the same number of observed cases in the present, R𝑒𝑓𝑓 must increase. A similar observation can be made for R𝑒𝑓𝑓<1.
In the context of our analyses (Figure 2), when the estimated R𝑒𝑓𝑓 is above 1, assuming a longer mean serial interval would further increase the R𝑒𝑓𝑓 estimates in each jurisdiction (i.e., the upper 75% of the Victorian posterior distribution for approximately the first 7–10 days, while stretching the upper tails in the other jurisdictions). When the estimated R𝑒𝑓𝑓 is below 1, a higher mean serial interval would further decrease those estimates. Qualitatively, this does not impact on the time series of R𝑒𝑓𝑓 in each Australian jurisdiction.
A prior distribution was specified for R𝑒𝑓𝑓, with mean 2.6 (informed by Imai et al., 2020) and a broad standard deviation of 2 so as to allow for a range of R𝑒𝑓𝑓 values. Finally, R𝑒𝑓𝑓 is estimated with a moving average window, selected to optimise the continuous ranked probability score, in order to smooth the curve and reduce the impact of localised events (i.e., cases clustered in time) causing large variations.
Note that up to 20% of reported cases in the Australian national COVID-19 database do not have a reported import status (see Figure 1). Conservatively, we assumed that all cases with an unknown or unconfirmed source of acquisition were locally acquired.

Accounting for imported cases
A large proportion of cases reported in Australia from January until now were imported from overseas. It is critical to account for two distinct populations in the case notification data – imported and locally acquired – in order to perform robust analyses of transmission in the early stages of this outbreak. The estimated time-varying R𝑒𝑓𝑓 value is based on cases that have been identified as a result of local transmission, whereas imported cases contribute to transmission only (Thompson et al., 2019).
Specifically, the method assumes that local and imported cases contribute equally to transmission. The results under this assumption are presented in Figure 2. However, it is likely that imported cases contributed relatively less to transmission than locally acquired cases, as a result of quarantine and other border measures which targeted these individuals (Figure 1—figure supplement 2). In the absence of data on whether the infector of local cases was themselves an imported or local case (from which we could robustly estimate the contribution of imported cases to transmission), we explored this via a sensitivity analysis. We aimed to explore the impact of a number of plausible scenarios, based on our knowledge of the timing, extent and level of enforcement of different quarantine policies enacted over time.
Prior to 15 March, returning Australian residents and citizens (and their dependents) from mainland China were advised to self-quarantine. Note that further border measures were implemented during this period, including enhanced testing and provision of advice on arrivals from selected countries based on a risk assessment tool developed in early February (Shearer et al., 2020). On 15 March, Australian authorities imposed a self-quarantine requirement on all international arrivals, and from 27 March, moved to a mandatory quarantine policy for all international arrivals.
Hence for the sensitivity analysis, we assumed two step changes in the effectiveness of quarantine of overseas arrivals (timed to coincide with the two key policy changes), resulting in three intervention phases: prior to 15 March (self-quarantine of arrivals from selected countries); 15–27 March inclusive (self-quarantine of arrivals from all countries); and 27 March onward (mandatory quarantine of overseas arrivals from all countries). We further assumed that the relative infectiousness of imported cases decreased with each intervention phase. The first two intervention phases correspond to self-quarantine policies, so we assume that they resulted in a relatively small reduction in the relative infectiousness of imported cases (the first smaller than the second, since the pre-15 March policy only applied to arrivals from selected countries). The third intervention phase corresponds to mandatory quarantine of overseas arrivals in hotels which we assume is highly effective at reducing onward transmission from imported cases, but allows for the occasional transmission event. We then varied the percentage of imported cases contributing to transmission over the three intervention phases, as detailed in Table 2.
Table 2.  Percentage of imported cases assumed to be contributing to transmission over three intervention phases for each sensitivity analysis.
We assume two step changes in the effectiveness of quarantine of overseas arrivals, resulting in three intervention phases: prior to 15 March (self-quarantine of arrivals from selected countries); 15–27 March inclusive (self-quarantine of arrivals from all countries); and 27 March onward (mandatory quarantine of overseas arrivals from all countries).
Imported cases contributing to transmission
Sensitivity analysis  Prior to 15 March  15–27 March  27 March–
1  90%  50%  1%
2  80%  50%  1%
3  50%  20%  1%
The results of these three analyses are shown in Figure 2—figure supplements 1, 2 and 3, respectively.

Forecasting short-term ward and ICU bed occupancy
We used the estimates of time-varying R𝑒𝑓𝑓 to forecast the national short-term ward/ICU occupancy due to COVID-19 patients.

Forecasting case counts
The forecasting method combines an SEEIIR (susceptible-exposed-infectious-recovered) population model of infection with daily COVID-19 case notification counts, through the use of a bootstrap particle filter (Arulampalam et al., 2002). This approach is similar to that implemented and described in Moss et al., 2019b, in the context of seasonal influenza forecasts for several major Australian cities. Briefly, the particle filter method uses post-regularisation (Doucet et al., 2001), with a deterministic resampling stage (Kitagawa, 1996). Code and documentation are available at https://epifx.readthedocs.io/en/latest/. The daily case counts by date of diagnosis were modelled using a negative binomial distribution with a fixed dispersion parameter k, and the expected number of cases was proportional to the daily incidence of symptomatic infections in the SEEIIR model; this proportion was characterised by the observation probability. Natural disease history parameters were sampled from narrow uniform priors, based on values reported in the literature for COVID-19 (Table 3), and each particle was associated with an R𝑒𝑓𝑓 trajectory that was drawn from the state/territory R𝑒𝑓𝑓 trajectories in Figure 2 up to 5 April, from which point they are assumed to be constant. The model was subsequently projected forward from April 14 to April 28, to forecast the number of reported cases, assuming a detection probability of 80%.
Table 3.  SEEIIR forecasting model parameters.
Parameter  Definition  Value/Prior distribution
σ  Inverse of the mean incubation period  U ⁢ ( 4 - 1   , 3 - 1   )
γ  Inverse of the mean infectious period  U ⁢ ( 10 - 1   , 9 - 1   )
τ  Time of first exposure (days since 2020-01-01)  U ⁢ ( 0 , 28 )
p 𝑜𝑏𝑠  Probability of observing a case  0.8
k  Dispersion parameter on Negative-Binomial  100
observation model In order to account for imported cases, we used daily counts of imported cases to construct a time-series of the expected daily importation rate and, assuming that such cases were identified one week after initial exposure, introduced exposure events into each particle trajectory by adding an extra term to the force of infection equation.
Model equations below describe the flow of individuals in the population from the susceptible class (S), through two exposed classes (E1, E2), two infectious classes (I1, I2) and finally into a removed class (R). The state variables S,E1,E2,I1,I2,R correspond to the proportion of individuals in the population (of size N) in each compartment. Given the closed population and unidirectional flow of individuals through the compartments, we evaluate the daily incidence of symptomatic individuals (at time t) as the change in cumulative incidence (the bracketed term in the expression for 𝔼⁢[yt] below). Two exposed and infectious classes are chosen such that the duration of time in the exposed or infectious period has an Erlang distribution. The corresponding parameters are given in Table 2.
Model equations:dSdt=−β(t)⋅S(I1+I2)dE1dt=β(t)⋅S(I1+I2)−2σE1dE2dt=2σE1−2σE2dI1dt=2σE2−2γI1dI2dt=2γI1−2γI2dRdt=2γI2
With initial conditions:S(0)=N−10NE1(0)=10NE2(0)=I1(0)=I2(0)=R(0)=0
Observation model:E[yt]=N⋅pobs⋅[I2(t)+R(t)−(I2(t−1)+R(t−1))]xt=[S(t),E1(t),E2(t),I1(t),I2(t),R(t),βi(t),σ,γ,τ]ℒ(yt∣xt)∼NegBin(E[yt],k)
With time-varying transmission rate corresponding to R𝑒𝑓𝑓 trajectory i:βi(t)={0,ift<τReffi(t)⋅γ,ift≥τ,for i∈{1,2,...,10}

Forecasting ward and ICU bed occupancy from observed and projected case counts
The number of new daily hospitalisations and ICU admissions were estimated from recently observed and forecasted case counts by:
Estimating the age distribution of projected case counts using data from the national COVID-19 database on the age-specific proportion of confirmed cases;
Estimating the age-specific hospitalisation and ICU admission rates using data from the national COVID-19 database. We assumed that all hospitalisations and ICU admissions were either recorded or were missing at random (31% and 58% of cases had no information recorded under hospitalisation or ICU status, respectively);
Randomly drawing the number of hospitalisations/ICU admissions in each age-group (for both the observed and projected case counts) from a binomial distribution with number of trials given by the expected number of cases in each age group (from 1), and probability given by the observed proportion of hospitalisations/ICU admissions by age group (from 2).
Finally, in order to calculate the number of occupied ward/ICU beds per day, length-of-stay in a ward bed and ICU bed were assumed to be Gamma distributed with means (SD) of 11 (3.42) days and 14 (5.22) days, respectively. We assumed ICU admissions required a ward bed prior to, and following, ICU stay for a Poisson distributed number of days with mean 2.5. Relevant Australian data were not available to parameterise a model that captures the dynamics of patient flow within the hospital system in more detail. Instead, these distributions were informed by a large study of clinical characteristics of 1099 COVID-19 patients in China (Guan et al., 2020). This model provides a useful indication of hospital bed occupancy based on limited available data and may be updated as more specific data (e.g., on COVID-19 patient length-of-stay) becomes available.