Transmission Dynamics and Effectiveness of Control Measures during COVID-19 Surge, Taiwan, April–August 2021

An unprecedented surge of COVID-19 cases in Taiwan in May 2021 led the government to implement strict nationwide control measures beginning May 15. During the surge, the government was able to bring the epidemic under control without a complete lockdown despite the cumulative case count reaching >14,400 and >780 deaths. We investigated the effectiveness of the public health and social measures instituted by the Taiwan government by quantifying the change in the effective reproduction number, which is a summary measure of the ability of the pathogen to spread through the population. The control measures that were instituted reduced the effective reproduction number from 2.0–3.3 to 0.6–0.7. This decrease was correlated with changes in mobility patterns in Taiwan, demonstrating that public compliance, active case finding, and contact tracing were effective measures in preventing further spread of the disease.

were fitted to a mixture of three unimodal (gamma, Weibull, and lognormal) distributions within the framework of a Bayesian mixture model. Each of the three component distributions was defined by the probability density function (∘; ) ( = 1,2,3) that was given a two-dimensional parameter vector = ( , ), where was the mean and was the coefficient of variation (CV).
Both parameters were common to all three component distributions.
Suppose that the time of the first event for a case ( ∈ ) was denoted by and the time of the second event was denoted by . When both times were recorded in the dataset, they were given by their calendar dates, and , respectively. Both times, oₖ and cₖ, were assigned to uniform priors of 1-day: ~ Uniform(max({ , }), + 1 day) where ≤ was held for all aggregated records, and < was imposed for simplicity for = . Because of a non-zero probability of viral RNA detection in a sample from a presymptomatic individual, the condition < can be invalid in general. However, when the alternative model with shifted distributions was verified, we found that it did not improve the data fit, and therefore chose to use the simpler model version.
When the date of the first event, , was not recorded, the case was asymptomatic at the time of testing, . According to the data collection procedures of the Taiwan CDC, the date of testing is backlogged by 1 day or occurs on the same day as the reporting day ℛ if confirmation of the case represents social significance defined based on expert opinion within the agency. Because the severity status of all confirmed cases was regularly updated by the Taiwan CDC, we noticed that some cases with an initially unknown symptom onset date became severe at later dates, denoted as . Suppose that two conditions are additionally met: ≤ and ℛ < k < . The prior (1) can then be replaced with the following prior valid for that subset of cases: where (∘; ) is the cumulative distribution function of (∘; ), and is the date of the latest release of data by Taiwan CDC, set as 2 p.m. on of August 25, 2021. The probability of selecting the best-fit distribution could be then given by categorical sampling from the three components with relative probabilities / Σ ( = 1,2,3).

Time-varied reporting delay
The effective reproduction number, Rₜ, and case-fatality ratio are often estimated retrospectively when an epidemic has already been declared over, but they can also be estimated in real-time. However, real-time estimation presents unavoidable challenges to researchers (1,2 Taiwan's COVID-19 epidemic wave in mid-2021. In our framework, the coefficient of variation (CV) of the reporting delay was set constant, but the mean reporting delay varied over multiple lengthscales similarly to previous work (3). Two lengthscales were used, each characterized by two parameters: , which defined a long-range lengthscale and , which defined a short-range lengthscale. Their implementation in our statistical framework is described in more detail below. The baseline values for and were set to 10 days and 7 days, respectively (Appendix Figure 1). The sensitivity analysis included other values, = {5 days, 15 days} and = {3 days, 11 days}, demonstrating relative robustness (Appendix Figure 2). A comparison with a piecewise constant reporting delay was also performed (Appendix Figure 3).
For the long-range variation in the reporting delay controlled by the parameter , the temporal dynamics of its mean were modeled using a cubic B-spline, ( ) = , ( ), such that: The likelihood for a short-range variation in the reporting delay controlled by the parameter was defined as follows: where the set ( ) consisted of all cases whose symptom onset date (or confirmation date ) was within the time window : -≤ ≤ + (or -≤ ≤ + ), and are the inferred symptom onset and confirmation time: Because the reporting delay could be specified either by the date of symptom onset (a "forward-looking" reporting delay) or by the date of confirmation (a "backward-looking" reporting delay according to notations (5)), the two possible reporting delay distributions could lead to different results. The forward reporting delay, indicated by "(o)" in our work, e.g., ( ) , describes the situation in which the mean reporting delay was estimated using the left-end symptom onset date. By comparison, the backward reporting delay, indicated by "(c)", e.g., ( ) , describes the situation in which the mean reported delay was estimated using the right-end confirmation date. Depending on the concrete application, one form can be preferred over another and would allow different interpretations. For instance, when cases with unknown symptom onset dates were backprojected to their symptom onset dates, the form ( ) could be appropriately used. When the cases were nowcasted to a given confirmation date to predict the number of cases that had not yet been reported, the form ( ) was appropriate. Park and colleagues (5) argued that the backward-looking delay distribution is prone to a bias toward lower values when an exponential growth in cases is observed during an epidemic, whereas the forward-looking delay distribution is likely to be representative of the unbiased (intrinsic) delay distribution.
Thus, the reporting delay distribution was modeled based on a mixture of three unimodal distributions (gamma, Weibull, and lognormal). The plausibility of each distribution was determined by its relative weight in the likelihood Σ , when the probability of selecting the distribution was equal to / Σ (see equations (7) and (8)). The temporary change in the reporting delay was modeled by a time-varied mean reporting delay, whereas the CV of the reporting delay was kept unchanged. The alternative formulation, where the standard deviation (SD) was used instead of the CV and kept constant, did not fit the data, so it was omitted from our analysis. The time variation of the mean reporting delay was modeled at two timescales: for the long timescale, the change in ( ) was defined by a cubic B-spline function, while for the short timescale the reporting delay distribution was adjusted locally using the likelihood (7). The

Effective reproduction number by date of symptom onset
To estimate the effective reproduction number by the date of symptom onset, one needs to know whether the case is asymptomatic or symptomatic. If the case is symptomatic, the symptom onset date should be recorded. However, the extracted dataset was composed of two types of case records-some records contained the symptom onset date and the confirmation date, while others had a blank symptom onset date. For the latter, the asymptomatic status was not definitive, and the record could not differentiate truly asymptomatic and pre-symptomatic cases. The unknown symptom onset dates were backprojected from the case confirmation dates using the reporting delay distribution ( ) .
Let be the daily counts of cases with known symptom onset dates at day of symptom onset , and be the daily counts of cases with unknown symptom onset dates at day of confirmation . The expected daily count * that combines both was defined as follows: To account for right truncation, the expected counts * were nowcasted according to (6): where the function ( ) (∘, ) is a cumulative distribution function of the forward reporting delay distribution ( ) (∘, ).
The effective reproduction number by date of symptom onset, ( ) , was estimated using the renewal process written within the negative binomial likelihood: where is the serial interval distribution estimated for Taiwan consisting of the best-fit lognormal distribution with a mean (±SD) of 4.6 ± 3.4 days, which is in line with previous reports (7,8). The overdispersion parameter followed one of the four functional forms: (i) 1 ≡ const, which is a commonly used form of the negative binomial distribution, also known as NB2, (ii) 2 ≡ const ⋅ ( ), which resembles a quasi-Poisson distribution, also known as NB1, (iii) 3 ≡ const ⋅ ( ( )) 1/2 , which presents an intermediate situation between (i) and (ii), and (iv) 4 ≫ 1, which is a limiting case leading to the Poisson likelihood. Green and colleagues (9) previously indicated a better fit of the data when negative binomial likelihood (i) is used instead of Poisson likelihood (iv). However, other studies (10,11) have shown that the overdispersion parameter given by (iii) is a better choice than (i) or (ii). In our framework, all four likelihoods (i)-(iv) were implemented within the Bayesian mixture model with a similar likelihood as in (8). The best-fit configuration was directly selected from the results of statistical inference.
Finally, the effective reproduction number was smoothed by calculating the centered rolling average over a 5-day time window (13), which implied that the value ( ) on day t was estimated using five sequential values of the nowcasted counts on days { − 2, − 1, , + 1, + 2}.

Effective reproduction number by date of infection
The number of infections on day was generated by nowcasted cases with symptom onset dates on either side of that day and was proportional to the effective reproduction number by the date of infection ( ) . Framed within the renewal process framework (14), the equation holds the form: where denotes the profile of infectiousness, which is the probability of a secondary transmission from a primary case at time measured with respect to the symptom onset date.
Notably, covers both positive and negative integers. Negative integers imply pre-symptomatic transmission, which is characteristic of SARS-CoV-2 infection. The profile of infectiousness was previously estimated using a gamma distribution shifted 12.3 days to the left that peaked at time zero and indicated 44% of pre-symptomatic infections (15).
Because newly infected cases experience symptom onsets that are time-lagged by the incubation period, we can write: where ℎ defines the incubation period fitted by a lognormal distribution with a mean of 5.2 days (16). Combining equations (13)- (14) results in the renewal equation previously derived by Nakajo and Nishiura (14): where = 13 days is the least integer for the shift of 12.3 days. The observation model was implemented using a negative binomial likelihood function analogous to equation (13), however, ( ) was replaced by a double sum (15).
The effective reproduction number was first modeled by a piecewise constant function of time t. The change points were set to every 7 days, which was a rough approximation of both the generation time (17)(18)(19) and a calendar week. Two alternative approaches were also explored. In the first approach, the change in the effective reproduction number was attributed to PHSMs (Appendix Table 2). In this case, the effective reproduction number remained constant during the intervening time periods but changed abruptly when the PHSMs were implemented. In the second approach, the effective reproduction number was correlated with a change in mobility metrics and involved a sigmoidal, monotonically decreasing change in the baseline (basic) reproduction number, for example, owing to behavioral changes or more efficient case finding and contact tracing toward the end of the epidemic in August-September 2021. The effect of overall mobility reduction was modeled by a combination of mobility metrics in different settings. In this case, the effective reproduction number changed on the log scale according to the following equation: where the intercept log 0, is the logarithm of the baseline reproduction number given by a sigmoidal function over time (20): where each of five parameters, ₀, , , , and , were assigned to a weakly informative prior.
The parameter 0 defined the baseline reproduction number at the beginning of the epidemic.
The parameter measured the reduction in transmission during the later stage of the epidemic, The restricted number of mobility metrics min ≤ was identified by comparing the "leave-one-out information criterion" LOOIC values for the respective models with the total min metrics and min combinations among all metrics. The LOOIC is commonly used in Bayesian frameworks for model selection (21). Specifically, the fit of the model with only min mobility metrics was compared to the fit of other models by selecting the model with the smallest LOOIC value. The negative binomial likelihood used for the data fit was analogous to (11) and was sequentially replaced by the gamma likelihood (12).

Counterfactual scenarios
The statistical model for the effective reproduction number based on mobility metrics was applied to predict the incidence of COVID-19 under different counterfactual scenarios. The model with smaller values of DSS or RMSE was preferred over the others.
To predict the case counts ̃ and their posterior distribution for a given counterfactual or baseline scenario, the following two-step procedure was implemented. In the first step, the number of infections ̃ was obtained using the renewal process and the Poisson count model: for any > 0 . Otherwise, ̃= 0 = 0 and ̃< 0 = 0. The effective reproduction number

Summary
Our framework has several important implications to advance the methods used for the real-time estimation of epidemiological parameters. First, we estimated a time-varied reporting delay distribution that measured the time between symptom onset and case confirmation. Similar to previous works (3,24,25), we performed statistical inference within a Bayesian framework and conducted Bayesian nowcasting of cases not yet reported. Our approach was based on using cubic B-splines, which is less computationally demanding than Gaussian processes (3). Our method shows a similar performance to approximating the reporting delay by a piecewise constant function, but it forces the coefficient of variation to remain constant (Appendix Figure   3). It appeared that the method was robust to varying the parameters (Appendix Figure 2).
Second, we statistically inferred the effective reproduction number according to the date of infection, as proposed by Nakajo and Nishiura (14). • Suspend all outdoor gatherings with 10+ participants and indoor gatherings with 5+ participants. • Apart from essential services, which include maintenance venues, medical and public services, all other business and public venues must be closed. • In the community with an ongoing active transmission where the rapid containment is required, residents must comply with COVID-19 testing protocols, do not leave the pre-defined control places and also suspend all gatherings and close schools. Strengthened level 3 Five additional measures on top of that in Level 3: • Violation of mask wearing regulation at all times will be fined.
• Strict inspection of entertainment related venues which have already been announced to be closed. Illegal operation, which involves operators, on-site practitioners, consumers and participants will be penalized accordingly with the law. • Only take-out for all food service is possible. The enhanced crowd control is required for supermarkets and all other markets. • Cancellation of all banquets for wedding ceremonies and public memorial ceremonies for funerals. • Suspension of all religious gatherings. Religious venues must be closed for the public. Level 4 Rapid increase of local cases (average number of cases more than 100 per day within 14 days) and half of them with unknown source of infection • Residents can leave their home only for essential activities, for example, to buy food, receive medical treatment, or to do essential work. The social distance must be maintained and wearing a mask must be done at all times outdoors. • Residents must wear a mask and maintain social distance when at home. • Suspension of all gatherings. • Apart from essential services, which include maintenance venues, necessary medical and public services, all other venues and schools must be closed. • Implementation of a lockdown policy in cities/counties or districts that reported severe outbreaks. The lockdown areas must be defined precisely and clearly, which will include the definition of restriction zones for people's entrance and exit. Residents must stay home.