This study involves secondary analysis of existing aggregate research data, including SARS-CoV-1 data from Leung et al [24], MERS data from Cowling et al [7], and COVID-19 data that are publicly available [25]. All data were deidentified and only aggregate data were used. The University of Hong Kong/Hospital Authority Hong Kong West Institutional Review Board approved the secondary analysis without requiring further consent. As participants were not directly involved in the current research activities, no compensation was provided to the individuals.

Back-calculation was first designed for short-term prediction of diseases with a long incubation period [26] and has also been applied to point-source outbreaks with limited onward transmission, such as Legionnaires disease [27]. We extended the method and focused on estimating the final size of symptomatic or laboratory-confirmed cases at the early stage of a coronavirus SSE, and asymptomatic cases were not considered in the estimation if they were not confirmed and reported.

We assumed that (1) all cases in a given SSE were infected at around the same time (eg, infected by the same superspreader or environmental exposure within a short period of time); (2) all cases followed the same probability distribution function (PDF) of the incubation period, the same PDF of the onset-to-confirmation delay, and the same PDF of the generation time; and (3) there was limited onward secondary transmission when the SSEs were controlled by sufficient contact tracing and timely responses. For SARS and MERS, we assumed that the R_t of secondary transmission was 0, given their relatively lower potential of person-to-person transmission [2]. For COVID-19, which has much higher transmissibility and significant presymptomatic transmission, we assumed that the R_t of secondary transmission could be greater than 1, but the transmission was restricted within 1 disease generation. We assumed that testing and contact tracing policies were the same within 14 days, given that the generation time for SARS-CoV-1, MERS, and ancestral SAR-CoV-2 are around 7 days (Table 1).

Figure 1A-B summarizes the characteristics of coronavirus SSEs in disease surveillance. Here, the incubation period is the time interval between infection and symptom onset; in practice, there is a delay from symptom onset to case confirmation, and the time of symptom onset can be unavailable for asymptomatic infections or due to underreporting.

We formulated the model to estimate the epidemic size of SSEs as follows. Let X denote the incubation period with PDF g(x) and Y denote the onset-to-confirmation with PDF h(y). For an SSE j, let n_j be the total number of people infected at time 0 and k_j be the number of symptomatic cases that have been confirmed up to (and including) time t_j since the exposure time for this SSE.

Let x_i,j and y_i,j be the values of X and Y for the ith case in the jth SSE. Let f(X,Y|θ) be the joint probability distribution of X and Y with θ representing the parameters. We assume that X and Y are independent such that f(X,Y|θ) = g(x)h(y), and F(X,Y|θ) is the cumulative probability distribution of f(X,Y|θ). The number of symptomatic cases k_j has a binomial distribution with SSE size n_j and success probability F(X+Y>t; θ). Therefore, the likelihood function for jth SSE is:

If the time of symptom onset is known, let M_j,1 be the set of confirmed cases with known time from exposure to symptom onset and known time from symptom onset to confirmation with a size of m_j,1. If the time of symptom onset is unknown, let M_j,2 be the set of confirmed cases with unknown time from exposure to symptom onset but known time from exposure to confirmation with a size of m_j,2. Let Z denote X+Y with the probability distribution q(z|θ). The likelihood function for the jth SSE becomes:

The overall likelihood function for w SSEs is:

We estimated the SSE size {n_j} and other parameters θ using Markov chain Monte Carlo with Gibbs sampling and noninformative flat priors. We considered an estimate to be accurate if the 95% credible interval (Crl) covered the true value. An estimate was precise if the Crl had a relative error within 1.

We simulated SSE scenarios based on the epidemiological characteristics of SARS, MERS, and COVID-19 (Figure 1C). For SARS and MERS to be considered under control, we assumed that there was a single exposure and no secondary transmission. For COVID-19, we assumed that 60% of the transmission occurred during the presymptomatic phase [19,20,35,39,40]. For example, if we simulated a single-exposure COVID-19 SSE with 10 cases, 6 (60%) second-generation cases were infected by the other 4 first-generation cases in the SSE during the presymptomatic phase. However, the 6 second-generation cases would be taken as first-generation infections directly linked to the primary exposure when the presymptomatic transmission was undetected. With limited onward secondary transmission within 1 disease generation, R_t = 6/4 = 1.5 (Figure 1D). Details about the estimation of COVID-19 SSE sizes and the adjustment for presymptomatic transmission are described in Multimedia Appendix 1.

We used the Latin-hypercube sampling to generate 6 simulated coronavirus scenarios. The parameters for the incubation period, onset-to-confirmation delay, and generation time were sampled from the ranges observed from major SSEs of SARS, MERS and COVID-19 (Table 1). Since most of these 3 coronavirus SSEs were reported to consist of 10 to 300 cases [8,11], we first simulated 1000 stochastic SSEs of 30, 50, 100, and 200 infections with and without onset-to-confirmation delay (ie, 8000 SSE scenarios). We assumed that prior distributions of the incubation period and onset-to-confirmation delay could be estimated by bootstrapping from 50 or 100 cases before the occurrence of the simulated SSEs.

The SSE size estimation was performed starting from the time when 5% of all cases were observed until the time when 95% of all cases were observed. To compare scenarios with different epidemic sizes, we used relative size (ie, percentage of the cases observed in an SSE) and SSE duration (ie, percentage of the duration between the time of confirmation of the first and last case) in the assessment of the accuracy of the estimation.

We applied our framework to the Amoy Gardens SSE for SARS in Hong Kong, nosocomial SSEs for MERS in South Korea, and 2 SSEs from the COVID-19 pandemic in Hong Kong. An alert would be issued when 10 or more cases were confirmed to link with the same index case, by which we assumed that we had minimal data to start the SSE size estimation.

Inferred from the model simulation, we found that the accuracy and precision of our estimation increased with the size of the SSE (Figures 2 and 3). For SARS and MERS, with and without onset-to-confirmation delay, all Crls covered the true SSE size and the relative error was less than 1 for all simulated scenarios, when the estimation was performed at the time exceeding 50% of the duration of the SSE or when more than 60% of the cases were observed (Figure 2A-B). With the presence of onset-to-confirmation delay, the range of relative error doubled when the estimation was performed before 50% of the duration of the SSE or when less than 60% of the cases were observed (Figure 2C-D). Our method tended to overestimate the size in the early stage of an SSE when less than 20% of cases were reported, especially in the scenarios with onset-to-confirmation delay.

If the presymptomatic transmission was not adjusted in the COVID-19 simulated SSEs, there would always be an underestimation even after 50% of the duration of the SSE or when more than 60% of the cases were observed (Figure 3A-D). For a typical COVID-19 SSE with 50 cases, the model would predict the final size to be around 30 when half of the cases were observed if second-generation infections cannot be distinguished from first-generation infections. If we could detect the transmission pattern of COVID-19 in time and adjust for presymptomatic transmission as per our prior information, the estimation would perform similarly well to the simulated scenarios of SARS and MERS when presymptomatic transmission was not included (Figure 3E-H).

In the case study of Amoy Gardens SARS SSE, the Crl covered the final SSE size and had a relative error within 1 in 4 days after the alert was issued, which was 25 days before the last case of the SSE was observed (Figure 4). At the time when the alert was issued, only 36.7% (114/311) of all cases were observed. The SSE size estimation became more accurate and precise thereafter.

For the series of South Korea MERS SSEs, the Crls of Cluster 1, Cluster 2, and Cluster 3 covered the final SSE sizes with a relative error within 1 in 2 days, 3 days, and 0 days after the alert was issued, respectively, which were 4 days, 15 days, and 9 days before the last case of the respective SSE was observed (Figure 5). At the time when alerts were issued, 62% (18/29), 52% (65/125), and 41% (11/27) of all cases in the respective SSEs were observed, and the estimation converged sooner for later SSEs.

For the Tao Heung COVID-19 SSE in Hong Kong, the Crl covered the final SSE size and had a relative error within 1 in 5 days after the alert was issued, which was 20 days before the last case of the SSE was observed. At this time, 76% (32/42) of all cases were observed (Figure 6). Similarly, for the Fulum COVID-19 SSE, the Crl covered the final SSE size and had a relative error within 1 in 7 days after the alert was issued and 5 days before the last case was observed. At the time of issuing the alert, 86% (38/44) of all cases had been observed (Figure 6).