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

Background The global impact of COVID-19 and the country-specific responses to the pandemic provide an unparalleled opportunity to learn about different patterns of the outbreak and interventions. We model the global pattern of reported COVID-19 cases during the primary response period, with the aim of learning from the past to prepare for the future. Methods Using Bayesian methods, we analyse the response to the COVID-19 outbreak for 158 countries for the period 22 January to 9 June 2020. This encompasses the period in which many countries imposed a variety of response measures and initial relaxation strategies. Instead of modelling specific intervention types and timings for each country explicitly, we adopt a stochastic epidemiological model including a feedback mechanism on virus transmission to capture complex nonlinear dynamics arising from continuous changes in community behaviour in response to rising case numbers. We analyse the overall effect of interventions and community responses across diverse regions. This approach mitigates explicit consideration of issues such as period of infectivity and public adherence to government restrictions. Results Countries with the largest cumulative case tallies are characterised by a delayed response, whereas countries that avoid substantial community transmission during the period of study responded quickly. Countries that recovered rapidly also have a higher case identification rate and small numbers of undocumented community transmission at the early stages of the outbreak. We also demonstrate that uncertainty in numbers of undocumented infections dramatically impacts the risk of multiple waves. Our approach is also effective at pre-empting potential flare-ups. Conclusions We demonstrate the utility of modelling to interpret community behaviour in the early epidemic stages. Two lessons learnt that are important for the future are: i) countries that imposed strict containment measures early in the epidemic fared better with respect to numbers of reported cases; and ii) broader testing is required early in the epidemic to understand the magnitude of undocumented infections and recover rapidly. We conclude that clear patterns of containment are essential prior to relaxation of restrictions and show that modelling can provide insights to this end. Supplementary Information The online version contains supplementary material available at (doi:10.1186/s12889-020-09972-z).


Appendix A Stochastic model details
In the main manuscript, we present our model as a nonlinear system of ordinary differential equations (ODEs) for ease of discourse. Due the variability in daily reported cases, however, it is important to take into account stochasticity in our analysis.

A.1 Model formulation
We implement our model as a stochastic compartmental epidemiological model consisting of the compartments, S (susceptible), I (undocumented infected), A (confirmed active), R (confirmed recovered), D (confirmed death), and R u (undocumented recovered). Transition events between compartments are Let X t = [S t , I t , A t , R t , D t , R u t ] T be the state vector of sub-population counts at time t > 0, and assume a well-mixed population of size P . Conditional on the state X t , the waiting time to the next occurrence of event E j is assumed to be exponentially distributed with rate parameter h j (X t ), where h j (X t ) is the hazard function for E j . The hazard functions can be interpreted as the instantaneous rate of events conditional on the current state.

A.2 Approximate stochastic simulation
While exact realisations of this process can be generated using event-based simulation [3,4], this is prohibitive within an approximate Bayesian computational setting with large population sizes and event numbers. Therefore, we apply a first order approximation to the integral over the interval [t, t + τ ) to obtain the tau-leaping approximation [5], where Y j ∼ Poisson (h j (X t )τ ) counts the number of times event j occurs in the interval [t, t + τ ). Simulations proceed as per Algorithm 1. For our simulations we use τ = 1 (days), and initial condition where A 0 , R 0 and D 0 come from the Johns Hopkins University COVID-19 data, and κ is the relative number of unobserved cases at t = 0.

Appendix B Bayesian analysis
We apply Bayesian inference to quantify uncertainty in the model parameters associated with event rates and community response, θ = [α 0 , α, β, γ, δ, η, n, κ, w A , w R , w D ], for country i using Johns Hopkins University data The task is to sample the posterior distribution with probability density given by Bayes' Theorem, where p(θ) is the prior, p(D i | θ) is the likelihood and p(D i ) is the evidence. For the remainder of this section, we omit the country index i for notational convenience.

B.1 Approximate Bayesian computation
For our model, the likelihood function can be written in terms of the solution to the forwards Kolmogorov equation, however, this requires large matrix exponentials to evaluate. Furthermore, since the full model state vector is only partially observable, the data model is no longer Markovian and is therefore computationally intractable. To deal with this likelihood intractability, we apply approximate Bayesian computation (ABC) [8][9][10], that samples from an approximation to the posterior for each country, where D is the COVID-19 data for the country of interest, D s ∼ s(· | θ) is simulated data, ρ(D, D s ) is a discrepancy metric, is the discrepancy threshold and 1 (0, ] (ρ(D, D s )) = 1 if ρ(D, D s ) ≤ , and 1 (0, ] (ρ(D, D s )) = 0 otherwise. For our implementation, we apply the discrepancy metric,

B.2 Sequential Monte Carlo sampling
We apply a sequential Monte Carlo (SMC) scheme [1,7] to move an initial set of N p samples from the prior through a sequence of ABC approximations defined by a decreasing sequence of T discrepancy thresholds, 1 > 2 > · · · > T = . Our particular implementation (Algorithm 2), based on the work of Drovandi and Pettit [2], adaptively selects the acceptance thresholds and utilises MCMC steps using tuned Gaussian random walk proposals. For all model calibrations we apply adaptive SMC with N p = 2000 particles, tuning parameters c = 0.01 and terminate sampling when the MCMC acceptance probably p acc drops below p min = 0.001.

11:
variate Gaussian density function and dim(θ) is the number of parameters;

B.3 Point estimate computation
Given a set of particles, {θ j } Np j=1 , from the ABC posterior generated from the SMC sampler, we require a point estimate for each parameter to compare overall trends between countries. We found the posterior means often did not result is a posterior predictive distribution that was representative of the data for most countries.
Instead we take the ABC posterior particle that has the smallest average discrepancy with the data. That is, for country i, the point estimate aŝ We evaluate this through direct Monte Carlo (Algorithm 3) using M = 100 simulations per particle. This section includes box-plots for comparison of marginal posterior distributions across countries for the model parameters θ = [α 0 , α, β, γ, η, n, κ, w A ]. Results for the first analysis period (22 January-30 March) are given the Fig. 1-3, results for the second analysis period (22 January-13 April) are given the Fig. 4-6, and results for the third analysis period (22 January-9 June) are given the Fig. 7-9. For reference, Table 1 provides the ISO-3166 alpha3 codes for each country.          Vertical bars indicate daily reported cases (yellow), recoveries (green) and deaths (red). The 50% (dark shaded region) and 95% credible intervals (light shaded region) of the posterior predictive distributions are plotted against the observational data.