key: cord-0847835-kmnb4qcp
authors: Jalilian, Abdollah; Mateu, Jorge
title: A hierarchical spatio-temporal model to analyze relative risk variations of COVID-19: a focus on Spain, Italy and Germany
date: 2021-03-23
journal: Stoch Environ Res Risk Assess
DOI: 10.1007/s00477-021-02003-2
sha: e902944daabfe50dba4869655cbf581634e97f03
doc_id: 847835
cord_uid: kmnb4qcp

The novel coronavirus disease (COVID-19) has spread rapidly across the world in a short period of time and with a heterogeneous pattern. Understanding the underlying temporal and spatial dynamics in the spread of COVID-19 can result in informed and timely public health policies. In this paper, we use a spatio-temporal stochastic model to explain the temporal and spatial variations in the daily number of new confirmed cases in Spain, Italy and Germany from late February 2020 to mid January 2021. Using a hierarchical Bayesian framework, we found that the temporal trends of the epidemic in the three countries rapidly reached their peaks and slowly started to decline at the beginning of April and then increased and reached their second maximum in the middle of November. However decline and increase of the temporal trend seems to show different patterns in Spain, Italy and Germany.

Started from Wuhan, the capital of Hubei province, China in December 2019, the outbreak of 2019 novel coronavirus disease has spread rapidly across more than 200 countries, areas or territories in a short period of time with so far (as of January 20, 2021) over 94.96 million confirmed cases and 2.05 million confirmed deaths (World Health Organization 2021) .

The spread of COVID-19 across and within countries has not followed a homogeneous pattern (Giuliani et al. 2020) . The causes of this heterogeneity are not yet clearly identified, but different countries have different levels of national capacity based on their abilities in prevention, detection, response strategies, enabling function, and operational readiness (Kandel et al. 2020) . Besides, different countries have implemented different levels of rigorous quarantine and control measures to prevent and contain the epidemic, which affect the population movement and hence the spread pattern of COVID-19. Given the highly contagious nature of COVID-19, the spatial pattern of the spread of the disease changes rapidly over time. Thus, understanding the spatio-temporal dynamics of the spread of COVID-19 in different countries is undoubtedly critical.

The spatial or geographical distribution of relative location of incidence (new cases) of COVID-19 in a country is important in the analyses of the disease risk across the country. In disease mapping studies, the spatial domain of interest is partitioned into a number of contiguous smaller areas, usually defined by administrative divisions such as provinces, counties, municipalities, towns or census tracts, and the aim of the study is to estimate the relative risk of each area at different times (Lee 2011; Lawson 2018) . Spatio-temporal models are then required to explain and predict the evolution of incidence and risk of the disease in both space and time simultaneously (Anderson and Ryan 2017).

Estimation of area-specific risks over time provides information on the disease burden in specific areas and identifies areas with elevated risk levels (hot spots). In addition, identifying the changes in the spatial patterns of the disease risk over time may result in detecting either regional or global trends, and contributes to make informed and timely public health resource allocation (Wakefield 2007) .

To account for the underlying temporal and spatial autocorrelation structure in the spread of COVID-19, available data on the daily number of new cases and deaths in different countries/regions have already been analyzed in a considerable number of studies. For example, Kang et al. (2020) used Moran's I spatial statistic with various definitions of neighbors and observed a significant spatial association of COVID-19 in daily number of new cases in provinces of mainland China. Gayawan et al. (2020) used a zero-inflated Poisson model for the daily number of new COVID-19 cases in the African continent and found that the pandemic varies geographically across Africa with notable high incidence in neighboring countries. Briz-Redón and Serrano-Aroca (2020) conducted a spatio-temporal analysis for exploring the effect of daily temperature on the accumulated number of COVID-19 cases in the provinces of Spain. They found no evidence suggesting a relationship between the temperature and the prevalence of COVID-19 in Spain. Gross et al. (2020) studied the spatiotemporal spread of COVID-19 in China and compare it to other global regions and concluded that human mobility/ migration from Hubei and the spread of COVID-19 are highly related. Danon et al. (2020) combined 2011 census data to capture population sizes and population movement in England and Wales with parameter estimates from the outbreak in China, and found that the COVID-19 outbreak is going to peak around 4 months after the start of personto-person transmission. Using linear regression, multilayer perceptron and vector autoregression, Sujath et al. (2020) modeled and forecasted the spread of COVID-19 cases in India.

As pointed out in Alamo et al. (2020) , there are many national and international organizations that provide open data on the number of confirmed cases and deaths. However, these data often suffer from incompleteness and inaccuracy, which are considerable limitations for any analyses and modeling conducted on the available data on COVID-19 (Langousis and Carsteanu 2020) . We highlight that we are yet in the center of the pandemic crisis and due to the public health problem, and also to the severe economical situation, we do not have access to all sources of data. Thus reseachers know only a portion of all the elements related to COVID-19. In addition, data on many relevant variables such as population movement and interaction, and the impact of quarantine and social distancing policies are not either available or accurately measured. Combined with the unknown nature of the new COVID-19 virus, any analysis such as the present study only provides an approximate and imprecise description of the underlying spatio-temporal dynamics of the pandemic. Nevertheless, having a vague idea is better than having no idea, and the results should be interpreted with caution.

Currently, a wealth of studies have appeared in the very recent literature. Many of them follow the compartmental models in epidemiology, partition the population into subpopulations (compartments) of susceptible (S), exposed (E), infectious (I) and recovered (R), and fit several variations of the classical deterministic SIR and SEIR epidemiological models (Peng et al. 2020; Roda et al. 2020; Bastos and Cajueiro 2020) . We believe that considering stochastic components is important, if not essential, to explain the complexity and heterogeneity of the spread of COVID-19 over time and space. For this reason, in the present work we propose a spatio-temporal stochastic modeling approach that is able to account for the spatial, temporal and interactions effects, together with possible deterministic covariates.

We acknowledge that the proposed model in its current form requires development and refinements as more information becomes available, but at the stage of the pandemic we are now, it can provide a reasonable modeling framework for the spatio-temporal spread of COVID-19. This is illustrated by modeling the daily number of new confirmed cases in Spain, Italy and Germany from late February 2020 to mid January 2021. The R code for implementing the proposed model can be made available upon request. We also provide a Shiny web application (Chang et al. 2020 ) based on the model discussed in this paper at https://ajalilian.shinyapps.io/shinyapp/.

The structure of the paper is the following. The open data resources used in this study are introduced in Sect. 2. A model for the daily number of regional cases is considered in Sect. 3. As described in Sect. 4, this model explains the spatio-temporal variations in the relative risk of each country in terms of a number of temporal, spatial and spatio-temporal random effects. The results of fitting the considered model to the number of daily confirmed cases in Spain, Italy and Germany are given in Sect. 5. The paper concludes in Sect. 6 with some few last remarks.

Governmental and non-governmental organizations across the world are collecting and reporting regional, national and global data on the daily number of confirmed cases, deaths and recovered patients and provide open data resources. Incompleteness, inconsistency, inaccuracy and ambiguity of these open data are among limitations of any analysis, modeling and forecasting based on the data (Alamo et al. 2020) . Particularly, the number of cases Table 1 summarizes the number of regions, study period and country-wide daily incidence rate of the data for each country.

Data on distribution population of the considered countries are extracted from the Gridded Population of the World, Version 4 (GPWv4), which provides estimates of the number of persons per pixel (1 resolution) for the year 2020 (Center International Earth Science Information Network (CIESIN) Columbia University 2018). These data are consistent with national censuses and population registers.

3 Modeling daily regional counts Suppose that a country, the spatial domain of interest, is partitioned into regions A 1 ; . . .; A m , defined by administrative divisions such as states, provinces, counties, etc (see Table 1 ). Let Y it denote the number of new COVID-19 cases in region A i , i ¼ 1; . . .; m, at time (day) t ¼ 1; . . .; T. Figure 1 shows boxplots of the observed number of new cases Y i1 ; . . .; Y iT for each region i ¼ 1; . . .; m. Clearly, observed values of Y i1 ; . . .; Y iT are heavily right-skewed and directly related to the population of region A i . Boxplots of the observed values of logð1 þ Y 1t Þ; . . .; logð1 þ Y mt Þ for each day t ¼ 1; . . .; m are shown in Fig. 2 . The roughly symmetric boxplots in the logarithmic scale imply right skewness of the regional counts at each day. Figure 2 also depicts temporal trends and variations of the observed number of COVID-19 cases which are due to restrictions and lockdowns imposed by national or regional government authorities and other relevant factors.

Thus, it is reasonable to assume that the expected number of new COVID-19 cases in region A i at time t is given by

where P i is the population of region A i and . it [ 0 is the incidence rate of COVID-19 in region A i at time t. Fig. 1 Boxplots of the daily number of COVID-19 cases y i1 ; . . .; y iT (left) and population P i (right) in each region (first level administrative division), i ¼ 1; . . .; m, of Spain, Italy and Germany

Under the null model of spatial and temporal homogeneity of the incidence rate, we have that

is an estimate of the country-wide homogeneous daily incidence rate (Waller et al. 1997) . The estimated daily incidence rate per million population (10 6 b . 0 ) so far is around 156, 123 and 79 for Spain, Italy and Germany, respectively (see Table 1 ).

3.2 Distribution of daily regional counts Consul and Jain (1973) introduced a generalization of the Poisson distribution, which is a suitable model to most unimodal or reverse J-shaped counting distributions with long right tails. In fact, the generalized Poisson distribution is a mixture of Poisson distributions and, compared to the ordinary Poisson and negative binomial distributions, it is a more flexible model for heavily right-skewed count data (Joe and Rong 2005; Zamani and Ismail 2012) . Given nonnegative random rates K it , i ¼ 1; . . .; m, t ¼ 1; . . .; T, following Zamani and Ismail (2012), we write

where u ! 0 and a 2 R, and assume that Y it 's are independent random variables following the generalized Poisson distribution with

Then,

Thus u is the dispersion parameter and the case u ¼ 0 represents the ordinary Poisson distribution (no dispersion) with

Here, parameter a controls the shape (power) of the relation between the conditional variance of Y it jK it and its conditional mean. For example, the relation between VarðY it jK it Þ and E½Y it jK it is linear if a ¼ 1, quadratic if a ¼ 1:5 and cubic if a ¼ 2 (Zamani and Ismail 2012).

The underlying random rates K it , i ¼ 1; . . .; m; t ¼ 1; . . .; T, account for the extra variability (overdispersion), which may represent unmeasured confounders and model misspecification (Wakefield 2007) . Variations of the random rate K it relative to the expected number of cases E it provide useful information about the spatio-temporal risk of COVID-19 in the whole spatial domain of interest during the study period.

In disease mapping literature, the nonnegative random quantities

are called the area-specific relative risks at time t ( Lawson 2018, Section 5.1.4). Note that the random rate K it and the expected value E it have the same dimension (unit is people per day) as Y it and hence h it is a dimensionless random latent variable. Since

which means that the temporal and spatial correlation structure of the underlying random rates K it determine the spatio-temporal correlations between h it 's. By ignoring these correlations, the standardized inci-

E it provides a naive estimate for the relative risks (Lee 2011). However, in a model-based approach the variations of the relative risks are often related to regional and/or temporal observed covariates and the correlation between h it 's are explained in terms of regional and/or temporal random effects using, for 

In the present study, we consider the log linear model

as the null model, where l is the intercept and d i is the population density of region A i , i.e. the population of A i , P i , divided by the area of A i . The population density is standardized to have mean 0 and variance 1 and b is its regression coefficient. However, Figs. 1 and 2 demonstrate spatio-temporal variations in daily number of COVID-19 cases which are not accounted for in the null model (2). We include zero mean random effects in the null model to account for spatio-temporal variations in relative risks due to temporal and spatial trend and correlation. Among many different possibilities, we consider the following models

where d t represents the temporal trend, e t accounts for temporal correlation and n i and f i explain spatial correlation due to spatial distance and neighborhood relations among regions A 1 ; . . .; A m , respectively (see Table 2 ).

Here, d t , e t , n i and f i are dimensionless random effects that explain temporal and spatial variability in the relative risk h it . The latent (stochastic) temporal trend d t is expected to be a smooth function of t. Since the second order random walk (RW2) model is appropriate for representing smooth curves (Fahrmeir and Kneib 2008) , d ¼ ðd 1 ; . . .; d T Þ is assumed to follow a RW2 model, i.e.,

where 2 ; . . .; TÀ1 are independent and identically distributed (i.i.d.) zero mean Gaussian random variables with variance 1=s d . Here the precision parameter s d [ 0 acts as a smoothing parameter enforcing small or allowing for large variations in d t (Fahrmeir and Kneib 2008) .

To account for temporal correlation, we assume that e t follows a stationary autoregressive model of order 2, AR(2); i.e.,

where À1\w 1 \1 and À1\w 2 \1 are the first and second partial autocorrelations of e t and 0 2 ; . . .; 0 T are i.i.d. zero mean Gaussian random variables with variance 1=s e .

On the other hand, the neighborhood structure of regions A 1 ; . . .; A m may induce spatial correlation among relative risks of regions because neighboring regions often tend to have similar relative risks. To include spatial correlation due to neighborhood structure of regions in the model, we assume that n ¼ ðn 1 ; . . .; n m Þ follows a scaled version of the Besag-York-Mollié (BYM) model (Besag et al. 1991) , i.e., n is a zero mean Gaussian random vector with (Riebler et al. 2016 )

Here Q À denotes the generalized inverse of the m  m spatial precision matrix Q ¼ ½Q ii 0 with entries

where n i is the number of neighbors of region A i and i $ i 0 means that regions A i and A i 0 share a common border. The parameter s n [ 0 represents the marginal precision and 0 / 1 indicates the proportion of the marginal variance explained by the neighborhood structure of regions (Riebler et al. 2016 ).

In addition to neighborhood interconnectivity, to account for spatial correlations due to spatial distance we assume that f ¼ ðf 1 ; . . .; f m Þ follows a Gaussian Markov random field (GMRF). More specifically, we assume that f b Fig. 2 Boxplots (yellow boxes with black whiskers) of logarithm of the regional number of COVID-19 cases logð1 þ y 1t Þ; . . .; logð1 þ y mt Þ in each day t ¼ 1; . . .; m for Spain, Italy and Germany. The red solid line connects the medians of boxplots and depicts the countrywide temporal trends and variations, and the blue dashed line is the overall median of all cases during the study period Small scale temporal trends due to temporal correlation AR(2) n i Spatial dependence due to neighborhood relation between regions BYM f i Spatial dependence due to distance between regions GMRF

is a zero mean Gaussian random vector with the structured covariance matrix

where I m is the m  m identity matrix, 0 x\1 and e max is the largest eigenvalue of the m  m symmetric positive definite matrix C ¼ ½C ii 0 . The entry C ii 0 of matrix C represents to what extend the regions A i and A i 0 are interconnected. For example, C ii 0 can be related to a data on commuting or population movement between regions A i and A i 0 . In absence of most recent and reliable movement data between the regions of Spain, Italy and Germany, we set C ii 0 to be the Euclidean distance between the centroids of A i and A i 0 .

In a Bayesian framework, it is necessary to specify prior distributions for all unknown parameters of the considered model. The Gaussian prior with mean zero and variance 10 6 is considered as a non-informative prior for the dispersion parameter of generalized Poisson distribution, log u, and for the parameters of the log linear model for the relative risks l, b, log s d , log s e , log s f , log s n , log x 1Àx , log / 1À/ , log 1þw 1 1Àw 1 and log 1þw 2 1Àw 2 . The prior distribution for the a parameter of the generalized Poisson distribution is considered to be a Gaussian distribution with mean 1.5 and variance 10 6 . Table 3 summarizes the model parameters and their necessary transformation for imposing the noninformative Gaussian priors.

Since all random effects of the model are Gaussian, the integrated nested Laplace approximation (INLA) method (Rue et al. 2009 ) can be used for deterministic fast approximation of posterior probability distributions of the model parameters and latent random effects (Martins et al. 2013; Lindgren and Rue 2015) . The R-INLA package, an R interface to the INLA program and available at www.rinla.org, is used for the implementation of the Bayesian computations in the present work. The R code can be made available upon request. The initial values for all parameters in the INLA numerical computations are set to be the mean of their corresponding prior distribution. The initial value of a is chosen to be 1.5 (see Table 3 ).

Let # denote the vector of all model parameters and

be the likelihood function for the observed count y it . The deviance information criterion (DIC)

and the Watanabe-Akaike information criterion (WAIC)

are two widely used measures of overall model fit, where b # is the posterior mode (Bayes estimate) of # and E pots and For count data Y it and in a Bayesian framework, a probabilistic forecast is a posterior predictive distribution on Z þ . It is expected to generate values that are consistent with the observations (calibration) and concentrated around their means (sharpness) as much as possible (Czado et al. 2009 ). Following a leave-one-out cross-validation approach, let

be the event of observing all count values except the one for region A i at time t. Dawid (1984) proposed the crossvalidated probability integral transform (PIT)

for calibration checks, where E ÀðitÞ post denotes expectation with respect to the posterior distribution of # based on the leave-one-out data B ÀðitÞ . Thus, PIT it is simply the value that the predictive distribution function of Y it attains at the observation point y it . The conditional predictive ordinate (CPO) (Pettit 1990 ). Moreover, the Bayesian leave-one-out cross-validation u\PIT it ðy it À 1Þ; u À PIT it ðy it À 1Þ PIT it ðy it Þ À PIT it ðy it À 1Þ

; PIT it ðy it À 1Þ u\PIT it ðy it Þ;

which is equivalent to

where 1½ Á is the indicator function, because

The mean PIT

can then be comparing with the standard uniform distribution for calibration. For example, a histogram with heights

and equally spaced bins ðj À 1Þ=J; j=J ½ Þ , j ¼ 1; . . .; J, can be compared with its counterpart from the standard uniform distribution with f j ¼ 1=J. Any departure from uniformity indicates forecast failures and model deficiencies. As mentioned in Czado et al. (2009) , U-shaped (reverse U-shaped) histograms indicate underdispersed (overdispersed) predictive distributions and when central tendencies of the predictive distributions are biased, the histograms are skewed.

The DIC, WAIC and BCV criteria for the fitted nested Table 4 . Compared with the null model M 0 , it can be seen that including the smooth temporal trend d t and spatial dependence due to neighborhood structure n i in the models for all three countries results in considerable improvement in the model fit (smaller DIC and WAIC) and prediction accuracy (larger BCV). Except for DIC criterion for Germany, presence of small scale temporal trends due to temporal correlation e t in the model slightly improves the model fit and prediction accuracy. Finally, including the spatial dependence due to distance between regions f i in the model leads to slightly worse model fit and prediction accuracy. Table 5 presents the Bayesian estimates (posterior means) for every parameter of the full model M 4 fitted to the daily number of new COVID-19 cases in Spain, Italy and Germany. The corresponding 95% credible intervals of the model parameters are also reported in parentheses.

Comparing the estimated parameters among different countries, it can be seen that the dispersion parameter u of the generalized Poisson distribution for Spain is lower than Italy and Germany, but its shape parameter a is around 1.5 for the three countries, which implies that the variance of the daily counts in each region is approximately a quadratic function of their mean. The coefficient of the population density is not significantly different from zero for Spain and Italy, but it is positive for Germany which indicates that regions with higher population density have larger relative risks. The opposite signs of w 1 and w 2 indicate rough oscillations in e t , which explains the small scale oscillatory temporal trends for the three countries, especially Germany.

The contribution of each random effect on the overall variation of relative risks can be quantified by c a ¼ 1=s a 1=s d þ 1=s n þ 1=s e þ 1=s f ; a ¼ d; n; e; f: Using the joint posterior distribution of precision parameters s d ; s n ; s e ; s f , estimates of the posterior means of c d , c n , c e and c f are obtained and summarized in Table 6 . In line with model fit and prediction accuracy criteria in Tables 4 and 6 reveals that more than 98% of the variations in the relative risks are explained by the temporal trend d t . Also, the spatial dependence due to neighborhood relation between regions n i and small scale temporal trends e t have minor contributions in the total variations of the relative risks for the three countries. The random effect term f i representing spatial dependence due to distance between regions has a negligible effect on the relative risks. The Bayesian estimates and 95% credible intervals for the temporal trend d t , t ¼ 1; . . .; T, are shown in Fig. 3 . These plots can be interpreted as a smoothed temporal trend of the relative risk in the whole country. In fact, Fig. 3 suggests that the COVID-19 epidemic in all three countries rapidly reached their first peaks in the middle of March and started to decline at the beginning of April and then increased and reached its maximum in the middle of November. The second wave of the epidemic seems to be differently affecting each of the three countries; it declined in December for Spain and Italy, but seems to be persistent in Germany. In addition, a third peak at the beginning of January is also apparent for Spain, which may be related to end-of-the-year holidays. Figure 4 shows the posterior means of the spatial random effects f i and n i , i ¼ 1; . . .; m, on the corresponding map of each country. The plot illustrates spatial heterogeneity of the relative risk of COVID-19 across regions in Fig. 4 Posterior mean of the spatial random effects n i representing dependence due to neighborhood relation between regions (left) and f i representing dependence due to distance between regions (right) in the fitted full model M 4 under the fitted model M 4 , namely the posterior mean and

Except some discrepancies for Spain and Italy, the observed values are inside the 95% credible intervals and close to the expected values under the fitted model. Figure 5 in addition shows 4-days ahead forecasts of the total daily number of new cases at the end of study period of each country.

Finally, histograms of the normalized PIT values described in Sect. 4.4 are obtained using J ¼ 20 from the fitted full models and plotted in Fig. 6 . The normalized PIT values for the fitted models to data do not show a clear visible pattern and the histograms seem to be close to the standard uniform distribution.

The above results and more details on observed and predicted values from the fitted full model are also provided in an interactive Shiny web application at https:// ajalilian.shinyapps.io/shinyapp/.

There are some limitations in the analyses and modeling of data on the number of new cases of COVID-19, including data incompleteness and inaccuracy, unavailability or inaccuracy of relevant variables such as population movement and interaction, as well as the unknown nature of the new COVID-19 virus. Nevertheless, understanding the underlying spatial and temporal dynamics of the spread of COVID-19 can result in detecting regional or global Observed value (solid line), predicted value (dashed line) and 95% prediction interval (grey area) for the daily number of new COVID-19 cases in the whole country, based on the posterior mean and 95% credible interval of the spatially accumulated relative risks of the fitted model trends and to further make informed and timely public health policies such as resource allocation.

In this study, we used a spatio-temporal model to explain the spatial and temporal variations of the relative risk of the disease in Spain, Italy and Germany. Despite data limitations and the complexity and uncertainty in the spread of COVID-19, the model was able to grasp the temporal and spatial trends in the data.

Obliviously, there are many relevant information and covariates that can be considered in our modeling framework and improve the model's predictive capabilities. One good possibility would be considering most recent and accurate human mobility amongst regions and replace the naive distance matrix C in the covariance matrix of the random spatial effect f i with a matrix constructed based on such mobility data. We would expect our model would benefit from this information, which right now can not be accessed. Moreover, the considered spatio-temporal model in this paper is one instance among many possibilities. For example, one possibility is to include a random effect term in the model that represents variations due to joint spatiotemporal correlations; e.g., a separable sptaio-temporal covariance structure. However, the considered nested models were adequate and no joint spatio-temporal random effect term was considered to avoid increasing the model's complexity.

We focused here on a stochastic spatio-temporal model as a good alternative to existing deterministic compartmental models in epidemiology to explain the spatio-temporal dynamics in the spread of COVID-19. However, it should be emphasized that one step forward would be considering a combination of a deterministic compartmental model in terms of differential equations for the number of susceptible, exposed, infectious and recovered cases with our sort of stochastic modeling approach. This is a clear novelty and a direction for future research.

Data Availibility Statement The datasets analysed during the current study are available on the GitHub repository and were derived from: https://github.com/datadista/datasets/tree/master/COVIDhttps:// github.com/pcm-dpc/COVID-19 (Italy) and https://github.com/ jgehrcke/covid-19-germany-gae. The R code for implementing the proposed model can be made available upon request. We also provide a Shiny web application (Chang et al. 2020 ) based on the model discussed in this paper at https://ajalilian.shinyapps.io/shinyapp/. A preprint of the manuscript reporting this study is available at https:// arxiv.org/abs/2009.13577.

Conflicts of interest The authors declare that they have no conflict of interest.

Open data resources for fighting Covid-19

A comparison of spatio-temporal disease mapping approaches including an application to ischaemic heart disease in new south wales, australia

Modeling and forecasting the Covid-19 pandemic in Brazil

Bayesian image restoration, with two applications in spatial statistics

A spatio-temporal analysis for exploring the effect of temperature on COVID-19 early evolution in Spain

Gridded population of the world, version 4 (gpwv4): population count, revision 11

Shiny: Web application framework for R

A generalization of the Poisson distribution

Predictive model assessment for count data

A spatial model of CoVID-19 transmission in England and wales: early spread and peak timing

Present position and potential developments: some personal views statistical theory the prequential approach

On the identification of trend and correlation in temporal and spatial regression

The spatio-temporal epidemic dynamics of COVID-19 outbreak in Africa. medRxiv Gelman A, Hwang J, Vehtari A (2014) Understanding predictive information criteria for Bayesian models

Modelling and predicting the spatio-temporal spread of coronavirus disease 2019 (COVID-19) in Italy

Spatio-temporal propagation of COVID-19 pandemics

Generalized Poisson distribution: the property of mixture of Poisson and comparison with negative binomial distribution

Health security capacities in the context of COVID-19 outbreak: an analysis of International Health Regulations annual report data from 182 countries

Spatial epidemic dynamics of the COVID-19 outbreak in China

Undersampling in action and at scale: application to the COVID-19 pandemic

A comparison of conditional autoregressive models used in Bayesian disease mapping

Bayesian spatial modelling with R-INLA

Bayesian computing with INLA: new features

Epidemic analysis of COVID-19 in china by dynamical modeling

The conditional predictive ordinate for the normal distribution

An intuitive Bayesian spatial model for disease mapping that accounts for scaling

Why is it difficult to accurately predict the COVID-19 epidemic?

Approximate Bayesian inference for latent Qaussian models by using integrated nested Laplace approximations

A machine learning forecasting model for COVID-19 pandemic in India

Disease mapping and spatial regression with count data

Hierarchical spatiotemporal mapping of disease rates

World Health Organization (2021) WHO coronavirus disease (COVID-19) dashboard

Functional form for the generalized Poisson regression model

Publisher's Note Springer Nature remains neutral with regard to jurisdictional claims in published maps and institutional affiliations

Acknowledgements The authors are thankful to two anonymous reviewers for their helpful comments. J. Mateu has been partially funded by research projects UJI-B2018-04 (Universitat Jaume I), AICO/2019/198 (Generalitat Valenciana) and MTM2016-78917-R (Ministry of Science).