key: cord-0152695-x04ctmpc
authors: Oliveira, Thiago de Paula; Moral, Rafael de Andrade
title: Global Short-Term Forecasting of Covid-19 Cases
date: 2020-05-29
journal: nan
DOI: nan
sha: ec56af63c9c2d05ab1fceda0eb86472b4b0bce05
doc_id: 152695
cord_uid: x04ctmpc

The continuously growing number of COVID-19 cases pressures healthcare services worldwide. Accurate short-term forecasting is thus vital to support country-level policy making. The strategies adopted by countries to combat the pandemic vary, generating different uncertainty levels about the actual number of cases. Accounting for the hierarchical structure of the data and accommodating extra-variability is therefore fundamental. We introduce a new modelling framework to describe the course of the pandemic with great accuracy, and provide short-term daily forecasts for every country in the world. We show that our model generates highly accurate forecasts up to six days ahead, and use estimated model components to cluster countries based on recent events. We introduce statistical novelty in terms of modelling the autoregressive parameter as a function of time, increasing predictive power and flexibility to adapt to each country. Our model can also be used to forecast the number of deaths, study the effects of covariates (such as lockdown policies), and generate forecasts for smaller regions within countries. Consequently, it has strong implications for global planning and decision making. We constantly update forecasts and make all results freely available to any country in the world through an online Shiny dashboard.

Outbreaks of the COVID-19 pandemic have been causing worldwide socioeconomic and health concerns since December 2019, putting high pressure on healthcare services. SARS-CoV-2, the causative agent of COVID-19, spreads efficiently and, consequently, the effectiveness of control measures depends on the relationship between epidemiological variables, human behaviour, and government intervention to control the spread of the disease 3, 4 . Different attempts to model the virus outbreaks in many countries have been made, involving mechanistic models 5, 6 , and extensions based on susceptible-infected-recovered (SIR) systems [7] [8] [9] . Countries have also put together task forces to work with COVID-19 data and their indirect impact on the population, economy, banking and insurance, and financial markets 2, 10-12 . In addition, funding agencies have put together rapid response calls worldwide for projects that can help dealing with this pandemic. However, further investment is still needed to foment priority research involving SARS-CoV-2, so as to establish a high level coordination of essential, policy-relevant, social and mental health science 13, 14 .

This pandemic is associated with high basic reproduction numbers 15, 16 , spreading with great speed since a significant number of infected individuals remain asymptomatic, while still being able to transmit the virus 17 . In the UK, for example, it is estimated that approximately 4.4% of people who tested positive for SARS-CoV-2 overall required hospitalisation, increasing to 17 − 27% in the group of persons aged 65 years or older. A pressing concern here is how to avoid bringing healthcare systems to a collapse 17, 18 . Knowing how the outbreak is progressing is crucial to predict whether or when this will happen, and therefore to plan and implement measures to reduce the number of cases so as to avoid it.

Policies for reducing the number of infected people such as social distancing and movement restrictions have been put in place in many countries, but for many others a full lockdown may be very difficult (if not impossible) to implement. This also depends heavily on the country's political leadership, socioeconomic reality, and epidemic stage 1, 19 . In this context, accurate short-term forecasting would prove itself invaluable, especially for systems on the brink of collapse and countries whose governments must consider trade-offs between lockdowns and avoiding full economical catastrophe.

The main problem is that not only is this disease new, but there are also many factors acting in concert resulting in a seemingly unpredictable outbreak progression. Forecasting with great accuracy under these circumstances is very difficult. Here we propose a new modelling framework, based on a state-space hierarchical model, that is able to generate forecasts with very good accuracy for up to seven days ahead. To aid policy making and effective implementation of restrictions or reopening measures, we provide all results as an R Shiny Dashboard, including week-long forecasts for every country in the world whose data is collected and made available by the European Centre for Disease Prevention and Control (ECDC).

Our model displayed an excellent predictive performance for short-term forecasting. We validated the model by fitting it to the data up to 13-May-2020, after removing the last seven time points (from 14-May-2020 until 20-May-2020), and compared the forecasted values with the observed ones ( Figure 1A ). We observe that it performs very well, especially for the first six days ahead ( Figure 1B ). Even though performance falls substantially for the seventh day ahead, with a concordance correlation between observed and forecasted values close to 0.5, there are still many countries for which the forecasted daily number of new cases is very close to the observed one. We carried out this same type of validation study using data up to 6-May-2020 and up to 29-Apr-2020, and the results were very similar, with a concordance correlation between observed and forecasted values greater than 0.75 for up to five days ahead (see Supplementary Materials). Figure 1 . A Logarithm of the observed y it versus the forecasted daily number of cases y * it for each country, for up to seven days ahead, where each day ahead constitutes one panel. The forecasts were obtained from the autoregressive state-space hierarchical negative binomial model, fitted using data up to 13-May-2020. The first day ahead corresponds to 14-May-2020, and the seventh to 20-May-2020. Each dot represents a country, and the sixteen countries shown in Figure 2 are represented by blue triangles. We add 1 to the values before taking the logarithm. B Observed accuracy, concordance correlation coefficient (CCC) and Pearson correlation (r) between observed (y it ) and forecasted (y * it ) values for each of the days ahead of 13-May-2020.

The autoregressive component in the model has a direct relationship with the pandemic behaviour over time for each country (see Supplementary Materials). It is directly proportional to the natural logarithm of the daily number of cases, given what happened in the previous day. Therefore, it is sensitive to changes and can be helpful detecting a possible second wave. See,

for example, its behaviour for Ireland, Spain, Italy and France -it shows that the outbreak is decaying, however it may still take time to subside completely ( Figure 2 ). On the other hand, in Iceland the outbreak has ended, according to the estimated autoregressive component. In China and South Korea, however, even though it appears that the outbreak has come to an end, our results show that a second outbreak could begin, and hence countries must be very cautious when relaxing restrictions. In the UK and the United States, it seems as if the outbreak is taking a long time to subside, and in Brazil and Russia it has not even reached its peak yet. The estimates for the model parameters suggest that about 10% of the number of reported cases can be viewed as contributing to extra variability, and possibly may consist of outlying observations (see the model estimates in the Supplementary Materials). These observations may be actual outliers, but it is likely that this is a feature of the data collection process. In many countries, the data that is recorded for a particular day actually reflects tests that were done over the previous week or even earlier. This generates aggregated-level data, which is prone to exhibiting overdispersion, which is accounted for in our model, but for some observations this variability is even greater, since they reflect a large number of accumulated suspected cases that were then confirmed. There is also a large variability between countries in terms of their underlying autoregressive processes (see the estimates for the variance components in the Supplementary Materials). This corroborates the fact that countries are dealing with the pandemic in different ways, and also may collect and/or report data differently.

We propose clustering the countries based on the behaviour of their estimated autoregressive parameter over the last 60 days ( Figure 3 ). This gives governments the opportunity to see which countries have had the most similar recent behaviour of the outbreak, and study similar or different measures taken by these other countries that may help determine policy. We observe, for example, that Brazil and Russia have been experiencing a very similar situation recently; the same can be said about France, Germany, Italy and Spain. Our R Shiny dashboard displays up-to-date results in terms of forecasts and also country clustering, and can be accessed at https://prof-thiagooliveira.shinyapps.io/COVIDForecast. Through the dashboard, users can choose to highlight a different number of clusters, which may provide different insights. 

Our modelling framework allows for forecasting the daily number of new COVID-19 cases for each country and territory for which data has been gathered by the ECDC. It introduces statistical novelty in terms of modelling the autoregressive parameter as a function of time. This makes it highly flexible to adapt to changes in the autoregressive structure of the data over time. In the COVID-19 pandemic, this translates directly into improved predictive power in terms of forecasting future numbers of 4/22 daily cases. Our objective here is to provide a simple, yet not overly simplistic, framework for country-level decision-making, and we understand this might be easier for smaller countries when compared to nations of continental dimensions, where state-level decision-making should be more effective 20 . Moreover, the model can be adapted to other types of data, such as number of deaths, and also be used to obtain forecasts for smaller regions within a country. A natural extension would be to include an additional hierarchical structure taking into account the nested relationship between cities, states, countries and ultimately continents, while also accommodating the spatial autocorrelation. This would allow for capturing the extra-variability introduced by aggregating counts over different cities and districts.

We remark that one must be very careful when looking at the forecasted number of cases. These values must not be looked at in isolation. It is imperative that the entire context is looked at and that we understand how the data is actually generated 21 . The model will obtain forecasts based on what is being reported, and therefore will be a direct reflection of the data collection process, be it appropriate or not. When data collection is not trustworthy, model forecasts will not be either. Our estimated number of outlying observations of approximately 10% is relatively high, and when looking at each country's time series we observe that countries that display unusual behaviour in terms of more outlying observations are usually related to a poorer model predictive performance as well. This suggests that the data collection process is far from optimal. Looking at the full context of the data is therefore key to implementing policy based on model forecasts.

As self-criticism, our model may possibly be simplistic in the sense that it relies on the previous observation to predict the next one, but does not rely on mechanistic biological processes, such as SEIR-type models 5, 22, 23 . These types of models allow for a better understanding of the pandemic in terms of the disease dynamics, and are able to provide long-term estimates 6 . However, as highlighted above, our objective is short-term forecasting, which is a task that is already very difficult. Long-term forecasting is even more difficult, and we believe it could even be speculative when dealing with its implications pragmatically. A more tangible solution to this issue is to combine the short-term forecasting power of our proposed modelling framework with the long-term projections provided by mechanistic models, so as to implement policy that can solve pressing problems efficiently, while at the same time looking at how it may affect society in the long run. We acknowledge all models are wrong, including the one presented here. We have showed, however, that it provides very good forecasts for up to a week ahead, and this can help inform on country reopening policies. Even though the model performs very well, we stress the importance of constantly validating the forecasts, since not only can the underlying process change, but also data collection practices change (to either better or worse). Many countries are still not carrying out enough tests, and hence the true number of infections is sub-notified 24 . This hampers model performance significantly, especially that of biologically realistic models 5 .

There is a dire need for better quality data 21 . This is a disease with a very long incubation period (up to 24 days 25 ), with median value varying from 4.5 to 5.8 days 26 , which makes it even harder to pinpoint exactly when the infection was actually transmitted from one person to the other. Furthermore, around 97.5% of those who develop symptoms will do so in between 8.2 to 15.6 days 26 . The number of reported new cases for today reflects a mixture of infections that were transmitted at any point in time over the last few weeks. The time testing takes also influences this. One possible alternative is to model excess deaths when compared to averages over previous years 27, 28 . This is an interesting approach, since it can highlight the effects of the pandemic (see, e.g., the insighful visualisations at Our World in Data 29 , and the Finantial Times 30 ). Again, this is highly dependent on the quality of the data collection process not only for COVID-19 related deaths, but also those arising from different sources. Many different teams of data scientists and statisticians worldwide are developing different approaches to work with COVID-19 data 3, 9, 16, 31 . It is the duty of each country's government to collect and provide accurate data. This way, it can be used with the objective of improving healthcare systems and general social wellbeing.

The developed R Shiny Dashboard displays seven-day forecasts for all countries and territories whose data are collected by the ECDC, and updated clustering of countries based on the last 60 days. This can help governments currently implementing or lifting restrictions. It is possible to compare government policies between countries at similar stages of the pandemic to determine the most effective courses of action. These can then be tailored by each particular government to their own country. The efficiency of measures put in place in each country can also be studied using our modelling framework, since it easily accommodates covariates in the linear predictor. Then, the contribution of these country-specific effects to the overall number of cases can serve as an indicator of how they may be influencing the behaviour of the outbreaks over time.

Government policies are extremely dependent on the reality of each country. It has become clear that there are countries that are well able to withstand a complete lockdown, whereas others cannot cope with the economic downfall 19 . The issue is not only economic, but also of newly emerging health crises that are not due to COVID-19 lethality alone, but to families relying on day-to-day work to guarantee their food supplies. For these countries, there is a trade-off between avoiding a new evil versus amplifying pre-existing problems or even creating new ones. In such circumstances, it is indeed very difficult to create a one-size-fits-all plan, which makes it even more vital to strive for better data collection practices.

We hope to be able to contribute to the development of efficient short-term response to the pandemic, for countries whose healthcare systems are at capacity, and countries implementing reopening plans. By providing a means of comparing recent

behaviour of the outbreak between countries, we also hope to provide a means to opening dialogue between countries going through a similar stage, and those who have gone through similar stages before.

The data was obtained from the European Centre for Disease Prevention and Control (ECDC), and the code is implemented such that it downloads and processes the up-to-date data from https://www.ecdc.europa.eu/en/geographical-distribution-2019-ncovcases. We remove, however, the data from the current day, since it can be updated by the ECDC. We assumed non-available data to be zero prior to the first case being recorded for each country. Whenever the daily recorded data was negative (reflecting previously confirmed cases being discarded), we replaced that information with a zero. This is the case for only nine out of 29,610 observations as of the 20th of May 2020.

According to the ECDC, the data on number of cases is based on reports from health authorities worldwide (up to 500 sources), which are screened by epidemiologists for validation prior to being recorded in the ECDC dataset.

We introduce a class of state-space hierarchical models for overdispersed count time series. Let Y it be the observed number of newly infected people at country i and time t, with i = 1, . . . , N and t = 1, . . . , T . We model Y it as a Negative binomial first-order Markov process, with

The mean is modelled on the log scale as the sum of an autoregressive component (γ it ) and a component that accommodates outliers (Ω it ), i.e.

To accommodate the temporal correlation in the series, the non-stationary autoregressive process {γ it } is set up as

where η it is a Gaussian white noise process with mean 0 and variance σ 2 η . Differently from standard AR(1)-type models, here φ it is allowed to vary over time through an orthogonal polynomial regression linear predictor using the time as covariate, yielding

where P q (·) is the function that produces the orthogonal polynomial of order q, with P 0 (x) = 1 for any real number x; β q are the regression coefficients and b b b i is the vector of random effects, which are assumed to be normally distributed with mean vector 0 and variance-covariance matrix Σ b = diag σ 2 b 0 , . . . , σ 2 b Q . By assuming φ it varying by country over time, we allow for a more flexible autocorrelation function. Iterating (1) we obtain

which is equivalent to a country-specific AR(1) process. On the other hand, if φ it = φ i = β 0 , then γ it = φ t−1 γ i1 + φ t−2 η i2 + φ t−3 η i3 + . . . + φ η it−1 + η it , which is equivalent to assuming the same autocorrelation parameter for all countries.

Finally, to accommodate extra-variability we introduce the observational-level random effect

where λ it ∼ Bernoulli(π) and ω it ∼ N(0, σ 2 ω ). When λ it = 1, then observation y it is considered to be an outlier, and the extra variability is modelled by σ 2 ω . This can be seen as a mixture component that models the variance associated with outlying observations.

To forecast future observations y * i,t+1 , we use the median of the posterior distribution of Y i,t+1 |Y it . This is reasonable for short-term forecasting, since the error accumulates from one time step to the other. We produce forecasts for up to seven days ahead.

We fitted models considering different values for Q. Even though the results for Q = 3 showed that all β q parameters were different from zero when looking at the 95% credible intervals, we opted for Q = 2 for the final model fit, since it improves convergence of the model, as well as avoids overfitting by assuming a large polynomial degree, while still providing the extra flexibility introduced by the autoregressive function (2) . This can change, however, as more data becomes available for a large number of time steps.

We fitted the model without using the observations from the last seven days and obtained the forecasts y * it for each country and day. We then compared the forecasts with the true observations y it for each day ahead, by looking initially at the Pearson correlation between them. We also computed the concordance correlation coefficient 32, 33 , an agreement index that lies between −1 and 1, given by

where ρ t is the Pearson correlation coefficient (a measure of precision), and C t is the bias corrector factor (a measure of accuracy) at the t−th day ahead. ρ t measures how far each observation deviated from the best-fit line while C b ∈ (0, 1] measures how far the best-fit line deviates from the identity line through the origin, defined as

is a scale shift and u = (µ 1 − µ 2 )/ √ σ 1 σ 2 is a location shift relative to the scale. When C b = 1 then there is no deviation from the identity line.

The model is estimated using a Bayesian framework, and the prior distributions used are

We used 3 MCMC chains, 2,000 adaptation iterations, 50,000 as burn-in, and 50,000 iterations per chain with a thinning of 25. We assessed convergence by looking at the trace plots, autocorrelation function plots and the Gelman-Rubin convergence diagnostic 34 .

All analyses were carried out using R software 35 and JAGS 36 . Model fitting takes approximately seven hours to run in parallel computing, in a Dell Inspiron 17 7000 with 10th Generation Intel Core i7 processor, 1.80GHz×8 processor speed, 16GB RAM plus 20GB of swap space, 64-bit integers, and the platform used is a Linux Mint 19.2 Cinnamon system version 5.2.2-050202-generic.

We used the last 60 values of the estimated autoregressive component to perform clustering so as to obtain sets of countries that presented a similar recent behaviour. First, we computed the dissimilarities between the estimated time seriesγ γ γ i between each pair of countries using the dynamic time warp (DTW) distance 37, 38 . Let M be the set of all possible sequences of m pairs preserving the order of observations in the form r = ((γ i1 ,γ i 1 ), . . . , (γ im ,γ i m )). Dynamic time warping aims to minimise the distance between the coupled observations (γ it ,γ i t ). The DTW distance may be written as

By using the DTW distance, we are able to recognise similar shapes in time series, even in the presence of shifting and/or scaling 38 .

Then, we performed hierarchical clustering using the matrix of DTW distances using Ward's method, aimed at minimising the variability within clusters 39 , and obtained ten clusters. Finally, we produced a dendrogram of the results of the hierarchical clustering analysis, with each cluster coloured differently so as to aid visualisation.

All code and datasets are made available as Supplementary Materials, and updated versions can also be downloaded from https://github.com/Prof-ThiagoOliveira/covid_forecast.

We produced a Shiny app displaying the latest forecasts for the next seven days that is updated twice a week, and can be accessed at https://prof-thiagooliveira.shinyapps.io/COVIDForecast. Table 1 displays parameter estimates for the model fitted to the data up to 19-May-2020. Figure 4 -5 depicts the two additional forecast validation procedures we have carried out. Figure 6 displays density plots for the sampled posterior distributions for each model parameter. Figures 7-15 Table 1 . Parameter estimates and associated 95% credible intervals (CI) for the fitted autoregressive hierarchical state-space negative binomial model. β 0 , β 1 and β 2 are the fixed effects associated with the time orthogonal polynomials of order 0, 1 and 2, respectively; σ b 0 , σ b 1 and σ b 2 are the standard deviations of the random effects associated with the time orthogonal polynomials of order 0, 1 and 2, respectively; σ η is the standard deviation of the random autoregressive process; π is the probability of an observation being an outlier and σ ω is the standard deviation of the mixture component ω it ; ψ is the dispersion parameter of the negative binomial distribution. Figure 4 . A Logarithm of the observed y it versus the forecasted daily number of cases y * it for each country, for up to seven days ahead, where each day ahead constitutes one panel. The forecasts were obtained from the autoregressive state-space hierarchical negative binomial model, fitted using data up to 29-April-2020. The first day ahead corresponds to 30-April-2020, and the seventh to 06-May-2020. Each dot represents a country, and the sixteen countries shown in Figure 2 are represented by blue triangles. We add 1 to the values before taking the logarithm. B Observed accuracy, concordance correlation coefficient (CCC) and Pearson correlation (r) between observed (y it ) and forecasted (y * it ) values for each of the days ahead of 29-April-2020. Figure 5 . A Logarithm of the observed y it versus the forecasted daily number of cases y * it for each country, for up to seven days ahead, where each day ahead constitutes one panel. The forecasts were obtained from the autoregressive state-space hierarchical negative binomial model, fitted using data up to 06-May-2020. The first day ahead corresponds to 07-May-2020, and the seventh to 13-May-2020. Each dot represents a country, and the sixteen countries shown in Figure 2 are represented by blue triangles. We add 1 to the values before taking the logarithm. B Observed accuracy, concordance correlation coefficient (CCC) and Pearson correlation (r) between observed (y it ) and forecasted (y * it ) values for each of the days ahead of 06-May-2020. Figure 15 . Posterior means of the autoregressive component γ it (solid lines) and associated 95% credible intervals (shaded areas) for 18 countries from the pool of 210 countries and territories in the data, from 1-Jan-2020 until 20-May-2020.

COVID-19: Mitigation or suppression? Arab

COVID-19 and finance: Agendas for future research

Temporal dynamics in viral shedding and transmissibility of COVID-19

Applying principles of behaviour change to reduce SARS-CoV-2 transmission

Mechanistic models versus machine learning, a fight worth fighting for the biological community?

Modelling the COVID-19 epidemic and implementation of population-wide interventions in Italy

Modeling infectious epidemics

Susceptible-Infected-Recovered (SIR) Dynamics of COVID-19 and Economic Impact

A Modified SIR Model for the COVID-19 Contagion in Italy

Dealing with sleep problems during home confinement due to the COVID-19 outbreak: practical recommendations from a task force of the European CBT-I Academy

Countries test tactics in 'war' against COVID-19

International Liaison Committee on Resuscitation: COVID-19 Consensus on Science, Treatment Recommendations and Task Force Insights

Multidisciplinary research priorities for the COVID-19 pandemic: a call for action for mental health science

Coronavirus disease (COVID-19) in neurology and neurosurgery: A scoping review of the early literature

Severe acute respiratory syndrome coronavirus 2 (SARS-CoV-2) and coronavirus disease-2019 (COVID-19): The epidemic and the challenges

Phase-adjusted estimation of the number of Coronavirus Disease 2019 cases in Wuhan

COVID-19 Pandemic, Corona Viruses, and Diabetes Mellitus

Impact of non-pharmaceutical interventions (NPIs) to reduce COVID-19 mortality and healthcare demand

Limiting the spread of COVID-19 in Africa: one size mitigation strategies do not fit all countries

State-level variation of initial COVID-19 dynamics in the United States: The role of local government interventions. medRxiv 1-14

Modelling COVID-19

Extending the SIR epidemic model

On the stability of an SEIR epidemic model with distributed time-delay and a general class of feedback vaccination rules

Substantial undocumented infection facilitates the rapid dissemination of novel coronavirus (SARS-CoV2)

Clinical characteristics of coronavirus disease 2019 in China

The Incubation Period of Coronavirus Disease 2019 (COVID-19) From Publicly Reported Confirmed Cases: Estimation and Application

Excess cardiovascular mortality associated with cold spells in the Czech Republic

Excess mortality associated with influenza epidemics in Portugal

FT Visual & Data Journalism team. Coronavirus tracked: the latest figures as countries fight to contain the pandemic

Forecasting the novel coronavirus COVID-19

A Concordance Correlation Coefficient to Evaluate Reproducibility

Longitudinal Concordance Correlation Function Based on Variance Components: An Application in Fruit Color Analysis

Inference from iterative simulation using multiple sequences

An R package providing interface utilities, model templates, parallel computing methods and additional distributions for MCMC models in JAGS

Dynamic Time Warping

TSclust: An R package for time series clustering

Ward's Hierarchical Agglomerative Clustering Method: Which Algorithms Implement Ward's Criterion?

We would like to thank Prof. John Hinde for helpful comments when preparing this manuscript.

T.P.O. and R.A.M. conceived and implemented the modelling framework, and wrote the manuscript. T.P.O. produced the Shiny visualisation app.

The authors declare no competing interests.