Hindsight is 2020 vision: Characterisation of the global response to the COVID-19 pandemic

Model fit

For almost all countries, the posterior predictive median lies close to the observed time series, and is almost always contained within the 95% credible intervals. Fig. 1 provides examples of the model performance for the United States (Fig. 1(A)), Italy (Fig. 1(B)), Australia (Fig. 1(C)), China (Fig. 1(D)), South Korea (Fig. 1(E)), and Taiwan (Fig. 1(F)). Fig. 1(A)-(C),(F) are representative of the accurate results we achieve across almost all 103 countries in this study. The model fits for China (Fig. 1(D)) and South Korea (Fig. 1(E)) are the only exceptions. The inaccuracy of the active trajectory for China is, at least, partly due to the changes in reporting methodology from China on 13 February [43]; however, this does not explain the overestimation of the early recovery period for both China and South Korea. We discuss potential model improvements to account for this in the Discussion section. Notwithstanding this, our model appears to be sufficiently accurate for a broad comparative analysis of global responses to COVID-19 during this period.

• Download figure
• Open in new tab

Figure 1:

Examples of model fit for (A) United States, (B) Italy, (C) Australia, (D) China, (E) South Korea, and (F) Taiwan for COVID-19 data up to 31 March. The median (solid lines), and 95% credible interval (shaded region) of the posterior predictive distributions is plotted for active cases (blue), recoveries (yellow), and deaths (red). This is compared against observations as recorded in data for active cases (circles), recoveries (squares), and deaths (crosses).

Assessing global response up to 31 March 2020

The marginal posterior distributions for the parameters α₀ (residual transmission rate), α (regulated transmission rate), β (case recovery rate), and δ (case fatality rate) (See Figs. S1 and S2 in Supplementary material) are similar across countries. This is expected as these parameters mainly describe the SARS-CoV-2 infection itself, rather than differing national responses. For the period up to 31 March, exceptions include: a higher uncertainty in death rate for Honduras, higher uncertainty in recovery rate for Nigeria and Venezuela, and increased uncertainty in the α₀ parameter of Afghanistan, Cameroon, Cuba, Honduras, Cambodia, Mauritius, Nigeria, Palestine, and Venezuela. The data also appear to be uninformative for the relative recovery rate, η since the posterior is very close to the prior (Fig. S3). The posterior marginals for the identification rate, γ, the regulatory parameter, n, and the initial infected scale parameter κ, relate to the nature of the initial outbreak and subsequent response. As a result, there is significant variation in these distributions across countries (Figs. S3 and S4).

To investigate global patterns, we compute, for each country, the marginal posterior means for the parameters γ, n, and κ. Fig. 2 shows the distribution of these point estimates across all countries, and highlights countries having the highest cumulative confirmed cases to date (31 March 2020) without a decline in active cases: the United States (USA), Spain (ESP), Italy (ITA), Germany (DEU), France (FRA), United Kingdom (GBR), Iran (IRN), and Switzerland (CHE). Fig. 2 also highlights countries that are recovering from the outbreak or prevented a severe outbreak as of 31 March: China (CHN) [23], South Korea (KOR) [30], Singapore (SGP) [41] and Taiwan (TWN) [39]. Australia (AUS) is also included as it was still undergoing an outbreak on 31 March, but significantly reduced transmission by 13 April.

• Download figure
• Open in new tab

Figure 2:

Pairwise scatter plots of point estimates of each assessed country (grey points, with lines to indicate 95% CI) for the key parameters related to the management of an COVID-19 outbreak up to 31 March 2020: regulatory parameter n; detection rate γ; and the relative initial undocumented cases κ. (A) n versus γ; (B) κ versus γ; and (C) κ versus n. The top ten countries for confirmed case counts are highlighted (red) along with countries that have managed to control the outbreak are highlighted (green). Labels identify the country by ISO-3166 alpha-3 code.

Excluding the United Kingdom, all of the countries undergoing a severe outbreak (in terms of cumulative confirmed cases) have point estimates of n < 1 (Fig. 2(A), (D)). This indicates a high tolerance response in which the transmission rate did not decline significantly until cases increased to larger numbers (See Fig. 6 in Materials and Methods). This is consistent with reported delays in response across Europe and the United States [9]. However, it is important to restate that a high tolerance response n < 1 is not equivalent to a slow response in time. In ascending order, the five smallest values for n are: China (n = 0.471; 95% CI: [0.464,0.477]), South Korea (n = 0.559; 95% CI: [0.551,0.567]), Italy (n = 0.618; 95% CI: [0.610,0.625]), Spain (n = 0.638; 95% CI: [0.630,0.647]), then Germany (n = 0.673; 95% CI: [0.663,0.683]). Of these five, China and South Korea are characterised by a higher identification rate γ > 0.5 (China γ = 0.549; 95% CI:[0.535,0.562]; South Korea γ = 0.632; 95% [0.619,0.644]) which is indicative of their strict testing, isolation, and tracing regimes [6, 23]. It is unexpected that the value of n for China is low, since strong lockdown measures were taken in Hubei [23]; we expect this is due to the absence of clinically diagnosed cases reported before to 13 February [43]. It is also important to note that small n does not mean that NPIs were not effective, but rather that they became effective at higher tallies of confirmed cases. France, Iran, Turkey, and the United States also have comparable values for γ, but did not experience a decline in active cases in the time frame of the study. These countries have a larger value for κ, indicating a large number of latent infections at an earlier stage of the outbreak, therefore a delay in the effect is expected. With κ = 1.586 (95% CI: [1.567, 1.606]), Iran is a particularly extreme example of this, as there was significant undetected spread of SARS-CoV-2 through the volumes of pilgrims passing through the city of Qom [15]. Of the severely affected countries, the United Kingdom has the closest signature to Taiwan and Singapore, albeit with a higher κ and lower n. The larger κ likely implies early unobserved transmission prior to abandonment of “herd immunity” targets in favour of social distancing, and closing of non-essential business and schools [16,27]. Australia has a very similar characterisation to the United Kingdom, with higher κ, but also higher n that likely reflects the introduction of travel bans and social gathering restrictions introduced in the second half of March.

Overall, our modelling approach identifies features of how communities globally respond to the COVID-19 pandemic. Without explicit modelling of specific NPIs, we obtain results that are consistent with the intervention strategies of the most severely affected countries. Globally, we see a low tolerance regulatory response with many countries having high values of n. Unfortunately, most countries also have large values of κ, indicating that delays in early detection have lead to a larger unobserved infectious population. As a result, the active case curve will take longer to turn around. This downturn could be brought forward by increasing γ through increased identification and quarantine/isolation capacity.

These results are for the period 22 Jan–31 March. In the next section, we consider an additional two weeks of data (the period 1-13 Apr), provide refined parameter estimates and offer insights into the evolution of the pandemic.

Re-assessing global response up to 13 April 2020

The fast moving evolution of epidemics and pandemics such as COVID-19 demand that model outputs and the corresponding effectiveness of mitigating measures are regularly re-assessed. We re-evaluate our model estimates and inferences for 103 countries that satisfied our inclusion criteria (See Materials and Methods) over the period 22 January–13 April 2020. These results reveal that as the pandemic evolves, we can learn more about the possible response histories and latent infected populations.

Updated point estimates, based on marginal posterior means, are computed for the parameters γ, n and k. The distribution of these point estimates is shown in Fig. 3(A)-(C) with the same countries highlighted as in Fig. 2. For the top eight countries still experiencing severe outbreaks on 31 March, there is a significant shift in these parameter estimates using data for the period 22 January-13 April. There is also a reduction in the uncertainty for some parameters overall (See Supplementary Material Fig. S5-S8). These shifts, demonstrated in Fig. 3(D)-(F), have some broad trends. In Fig. 3(D) and (F), all eight countries show a decrease in the parameter n, with France, Iran, Turkey, the United Kingdom, and the United States having the most substantial declines. This indicates that, in light of more recent data, the community response to the outbreak had a higher tolerance to growing case numbers than earlier data indicated. Fig. 3(D) and (E) indicate that most of these countries have a lower identification rate. However, it is promising to see a small increase in γ for Germany and Switzerland, and a large increase in γ for Italy and Spain; this is possibly a reflection of increased testing capabilities within these countries between 31 March to 13 April (https://ourworldindata.org/coronavirus; Accessed 14 April). Overall, there is also a decrease in κ, indicating that the number of early undocumented cases could be lower than previously thought. One concerning exception is Italy, with an increase in the point estimate with κ = 0.956 (95% CI: [0.924, 0.988]). Australia also shows a decrease in the value of n = 0.703 (95% CI: [0.693, 0.713]), but also a decrease in κ = 0.259 (95% CI: [0.240, 0.278]). This suggests the Australian community was overall more tolerant of higher case numbers, but the undocumented case numbers were lower than previously expected. Australia shows an increase in γ = 0.456 (95% CI:[0.440,0.472]) which is consistent with a high identification rate of symptomatic cases [32]. Comparing Fig 2(A)-(C) with Fig 3(A)-(C) highlights the dynamical nature on the COVID-19 pandemic. In particular, changes in the parameters n and κ indicate that our understanding of the past trajectory of the outbreak is also evolving. This is a particular strength of our modelling approach in that we can characterise how sensitive a community response is to growing case numbers. Our framework is also flexible enough to include community responses to other factors that may be prescribed through an arbitrary utility metric (see Materials and Methods, and Supplementary Material).

• Download figure
• Open in new tab

Figure 3:

Pairwise scatter plots of point estimates of each assessed country (blue points, with lines to indicate 95% CI) for the key parameters related to the management of an COVID-19 outbreak up to 13 April 2020: regulatory parameter n; detection rate γ; and the relative initial undocumented cases κ. (A) n versus 7; (B) κ versus 7; and (C) κ versus n. Highlighted countries are the same as in Fig. 2. Labels identify the country by ISO-3166 alpha-3 code. (D)-(F) Indication of parameter estimate shifts in light of new data up to 13 April 2020, for the severe outbreaks identified on 31 March 2020.

Assessment of outbreak management

As stated above, the regulatory parameter n does not provide insight into the full management of the outbreak. It simply provides a measure of how the transmission rate declines as a function of confirmed case counts (see Materials and Methods). We compute several derived quantities that provide a more complete view of how countries are managing the pandemic. Specifically, we consider the following quantities: the case death rate relative to the case recovery rate, δ/β; the case removal rate relative to identification rate, (β + δ)/γ; the relative strength of transmission rate regulation, , where f (·) is the transmission rate regulatory function (see Materials and Methods), and respectively C₀, and C_T, are the number of confirmed cases at the start of the outbreak (first day exceeding 100 confirmed cases) and at the end of the data period (13 April); and the relative reduction in the basic reproduction number [10], , where unfigure is the basic reproductive number at the start of the outbreak and unfigure is the basic reproduction number at the end of the data period (See Supplementary Material). Point estimates of these ratios are computed through expectations with respect to the full joint parameter posterior distributions (See Materials and Methods, and Supplementary Material) for each country.

Fig. 4(A) shows the distribution of point estimates related to management and treatment of confirmed cases. The effective management of COVID-19 cases within Germany can be inferred by a low value for the relative case fatality rate δ/β = 0.083 (95% CI [0.070,0.086]) and high value for the case removal rate relative to identification rate (β + δ)/γ = 1.12 (95% CI [0.76,1.48]). The United Kingdom is an outlier with a very high δ/β = 41.32 (95% CI:[2.64, 80.01]); we expect this is due to a discrepancy in the data for recoveries in the United Kingdom (with only 626 recoveries recorded). China is another outlier with (β + δ)/γ = 6.66 (95% CI: [0, 13.73]) due to relatively poor model fit for China as discussed above. Fig. 4(B) shows that globally there is a high level of regulation with (f (C₀) — f (C_T))/f (C₀) > 0.7 (that is, regulatory strength greater than 70% of maximum) for nearly all countries. Despite this, there are many countries with only moderate reductions in the basic reproductive number, indicating that more testing is required or a higher level of residual transmission a₀. In particular, the second outbreak occurring within Singapore in the later period of study is identified by a high level regulation but very little decrease in the basic reproductive number. Fortunately, the basic reproduction number reduced by at least 50% for all of the most severely affected countries.

• Download figure
• Open in new tab

Figure 4:

Scatter plots of point estimates of each assessed country (blue points, with lines to indicate 95% CI) for derived quantities from epidemiological rate parameters. (A) Comparison of the ratio δ/β with the ratio (β + δ)/γ, outliers not shown (China with (β + δ)/γ = 6.66 (95% CI: [0, 13.73]) and United Kingdom with δ/β = 41.32 (95% CI:[2.64, 80.01])). (B) reduction in effective reproduction number R_e on 13 April 2020 relative to the basic reproduction number R₀ compared with the relative strength of transmission rate regulation f (C_T) on 13 April 2020 relative to the transmission rate at the start of the outbreak. Highlighted countries are the same as in Fig. 2. Labels identify the country by ISO-3166 alpha-3 code.

Uncertainty quantification for predictions

The primary purpose of our model is to analyse the historical trajectory and community response to the COVID-19 pandemic. However, our approach also enables predictions that account for uncertainty in the undocumented infected population and in the community tolerance of rising cases. There are many unknowns within the current pandemic that affect the accuracy of model prediction [20,28]. In particular, we have no estimates for the number of undocumented cases world wide. Furthermore, inadequate uncertainty quantification can lead to unrealistic credible prediction intervals [28]. We use our Bayesian framework to handle both uncertainty in model parameters and uncertainty in undocumented infected population. Regarding prediction uncertainty, we take a conservative approach that provides predictions accounting for uncertainty in the evolution up to the current time. We perform this by sampling from the posterior predictive distribution and continue model evolution to the prediction horizon. In this way, we consider uncertainty in all possible histories (parameter values and latent population evolution) to quantify uncertainty in predictions conditional on the model formulation.

Using posterior distributions of parameters computed using data up to 13 April, we quantify prediction uncertainty over a 7 day prediction horizon. Example predictions for active cases, recoveries and deaths (median and 95% credible interval) are shown in Fig. 5 for the United States (Fig. 5(A)), Italy (Fig. 5(B)), Spain (Fig. 5(C)), Germany (Fig. 5(D)), Taiwan (Fig. 5(E)) and Australia (Fig. 5(F)). Our approach produced wide credible intervals as a result of parameter uncertainty and uncertainty in undocumented infections. There is a wide range of possible scenarios for undocumented infection populations sizes that are completely consistent with the data [26]. By continuing simulations from the posterior predictive distribution we take an ensemble of these scenarios into account.

• Download figure
• Open in new tab

Figure 5:

Examples of uncertainty quantification over γ day prediction horizon for (A) the United States, (B) Italy, (C) Spain, (D) Germany, (E) Taiwan and (F) Australia using COVID-19 data up to 13 April 2020. The median (solid lines), and 95% credible interval (shaded region) of the posterior predictive distributions is plotted for active cases (blue), recoveries (yellow), and deaths (red). Data used for predictions are also shown: active cases (circles), recoveries (squares), and deaths (crosses).

Our predictions for the United States show a continuing increase in the number of active cases (Fig 5(A)). An increase in active cases is predicted for Italy (Fig 5(B)), although the decline in the lower bound of the credible interval indicates the possibility of a flattening of active case counts in future weeks. Predictions for Spain and Germany (Fig 5(C),(D)) also indicate a flattening out in the number of active cases over the prediction horizon. For Taiwan and Australia (Fig 5(E), (F)) the median predictions show a decline in active cases; however, the uncertainty is still substantial.

The key message from this analysis is to highlight the importance of conservative uncertainties in the future evolution of the pandemic. It is encouraging to see predictions that indicate potential future declines in active cases counts. However, uncertainties in the number of undocumented cases and model parameters indicate caution in predicting the timing of consistent declines in active case numbers. Therefore, it is essential that communities remain vigilant in fast-evolving situations such as this.

Data summary

Daily counts of reported confirmed COVID-19 cases, recoveries and deaths for each country are obtained from the Johns Hopkins University coronavirus resource center (coronavirus.jhu.edu/; Accessed 14 April 2020) for the period of 22 January 2020 to 13 April 2020. Population data for each country for 2020 were obtained from United Nations Population Division estimates (www.worldometers.info/world-population/population-by-country/; Accessed 14 April 2020).

Our analysis is performed over two different periods: i) 22 January-31 March 2020; and ii) 22 January-13 April. Our inclusion criteria for period i) are countries for which there was at least one day exceeding 100 confirmed cases of COVID-19 between before 30 March 2020. Of all countries with confirmed cases, N = 103 countries satisfied our criteria for period i). Similarly, for period ii) we include countries with at least one day exceeding 100 confirmed cases of COVID-19 before 12 April 2020.

For each country, i = 1, 2, …, N, we process the raw data to produce a time-series, , where A_t,_i, R_t,i, and D_t,i are the active cases, recoveries and deaths on day t in country i. Here, day t = 0 is the first day to exceed 100 active cases. While it is known that there are variations in reporting protocols across countries and time, we make the necessary simplifying assumption in our analysis that the quantities A_t,i, R_t,i, and D_t,i have an equivalent interpretation across all countries over time.

Mathematical model

We apply a stochastic epidemiological model to describe the spread of COVID-19 within a single country that is a assumed to have a well-mixed population. Here, S is the population of those susceptible to SARS-CoV-2 infection, I is the latent infected population, A are the active confirmed cases of COVID-19, R are confirmed recoveries, D are confirmed deaths caused by COVID-19, and R^u are those that recover from COVID-19 but where never confirmed. The cumulative confirmed cases are given by C = A + R + D. The individual counts within each population are updated according to the following transition events, where is the effective transmission rate function, γ > 0 is the identification rate, β > 0 is the case recovery rate, δ > 0 is the case fatality rate, and n is the latent recovery rate relative to the case recovery rate. At time t > 0 (days), the system state is given by the population counts S_t, I_t, A_t, R_t, D_t, and and evolves according to a discrete-state continuous-time Markov process (Supplementary Material). We apply the tau-leaping numerical scheme [18] (Supplementary Material), as exact stochastic simulation was computationally infeasible. At time t = 0 (first day exceeding 100 active cases), A₀, R₀, and D₀ are obtained from time-series data. For unobserved states we initially set, I₀ = κA₀, , and S₀ = P — C₀ — I₀ where, P is total population record for the country and κ ≥ 0 is the number of unobserved cases relative to active cases at t = 0.

The transmission rate is described as a function of observables of the form . Here, f (A_t,R_t,D_t) is a decreasing function in a measure of utility U(A_t,R_t,D_t), a is a scale factor such that αf(A_t,R_t,D_t) is a transmission rate that decreases as the utility function increases and α₀ is the residual transmission rate as f (A_t, R_t, D_t) → 0. We assume the form f (A_t, R_t, D_t) = 1/(1 + U(A_t, R_t, D_t)ⁿ).

We find a utility function based on cumulative confirmed case numbers, U(A_t, R_t, D_t) = A_t + R_t + D_t = C_t, provides a good fit across all countries that we consider (excepting changes in reporting methodology for China on 13 Feb [43]). However, we also investigate a utility function that only considers deaths, U(A_t,R_t,D_t) = D_t, but model fitting was inferior (See Supplementary Material). There is a wide variety of utility functions that could be considered to model specific country NPIs and behaviours.

Therefore, we simplify our notation and consider where the parameter n ≥ 0 controls the rate of decrease as shown in Fig. 6. Given an initial number of confirmed cases, C₀, we have . Note the shape of f (C_t) starts to decline rapidly before starting to level out. Increasing n shifts the initial decline to the left, effectively altering the critical values of C_t in which a population begins to maximise efforts to “flatten the curve”. For example, small n describes a population that does not significantly reduce the transmission rate until the confirmed case tally is large. Conversely, larger n describes a population that acts decisively as a response to very early stages of the outbreak. Importantly, our approach does not distinguish between different NPIs and voluntary population behaviour, but rather models the net effect of rising confirmed cases on transmission rates.

• Download figure
• Open in new tab

Figure 6:

The effect of the regulatory parameter n on the regulation function. Higher values of n lead to a rapid decrease in transmission as C_t increases, whereas smaller n requires many more cases to reduce the transmission rate. (A) linear scale; (B) logarithmic scale.

The parameter κ is used to initialise the unobserved infected population, I₀ = κA₀, and indicates the outbreak extent within the community at time t = 0. Small κ indicates that efforts of testing and contact tracing have been very effective at early time, whereas large κ implies the virus has spread extensively without early intervention.

Bayesian Analysis

Using time-series data unfigure for each country i ϵ [1, 2, …, N], we infer model parameters θ_i within a Bayesian framework. For notational convenience, we will omit the subscript i for the remainder of this section. We sample the joint posterior distribution for θ = [α₀, α, β, γ, δ, η, n, κ], where is the likelihood function and p(θ) is the prior distribution. The likelihood function is intractable, due to the unobserved populations S_t, I_t, and . Therefore, we rely on adaptive sequential Monte Carlo for approximate Bayesian computation (SMC-ABC) [12,35,40] to obtain accurate approximate posterior samples (Supplementary Material). We use independent uniform priors, , , , , , , and .

Model fit is assessed though sampling the posterior predictive distribution where is simulated data as generated by the model. We compute the mean and the 95% credible region of unfigure through generating a single model simulation for each posterior sample generated by the SMC-ABC sampler and computing quantiles of the simulation state distribution at each observation time. For prediction, each stochastic simulation is continued in time up to the prediction horizon.

Parameter point estimates are obtained using the marginal posterior means and uncertainty is quantified using the marginal posterior variance. Point estimates of derived quantities, δ/β, (β + β)/γ, (f(C₀) — /(C_T))/f(C₀), and unfigure are computed through evaluating expectations with respect to the full joint posterior distribution.