key: cord-168710-a5pst4gf
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: 2020-09-28
journal: nan
DOI: nan
sha: 
doc_id: 168710
cord_uid: a5pst4gf

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 to mid September 2020. Using a hierarchical Bayesian framework, we found that the temporal trend 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 August. However decline and increase of the temporal trend seems to be sharper in Spain and smoother in Germany. The spatial heterogeneity of the relative risk of COVID-19 in Spain is also more pronounced than 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 over 4.4 million confirmed cases and 296 thousand confirmed deaths (World Health Organization, 2020) .

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 spatio-temporal 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 person-to-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 dynamic 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 to mid August 2020. 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 Section 2. A model for the daily number of regional cases is considered in Section 3. As described in Section 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 Section 5. The paper concludes in Section 6 with a 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 mainly consist of cases confirmed by a laboratory test and do not include infected asymptomatic cases and infected symptomatic cases without a positive laboratory test.

In this study, we focus on the daily number of confirmed cases in Spain, Italy and Germany and used the following open data resources. Italy: Data on the daily accumulated number of confirmed cases in the 20 regions of Italy are reported by the Civil Protection Department (Dipartimento della Protezione Civile), a national organization in Italy that deals with the prediction, prevention and management of emergency events. These data are available at the GitHub repository https://github.com/ pcm-dpc/COVID-19 and are being constantly updated every day at 18:00.

Germany: The Robert Koch Institute, a federal government agency and research institute responsible for disease control and prevention, collects data and publishes official daily situation reports on COVID-19 in Germany. Data on the daily accumulated number of confirmed cases in the 16 federal states of Germany extracted from the situation reports of the Robert Koch Institute are available at the GitHub repository https: //github.com/jgehrcke/covid-19-germany-gae and are being updated on a daily basis. 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 degree 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

The expected number of new cases is given by E it = P i it , where P i is the population of region A i and it is the incidence rate of COVID-19 in region A i at time t. Under the null model of spatial and temporal homogeneity of the incidence rate it = 0 and

provides an estimate for E it , where

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 0 ) so far is around 68, 46 and 29 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. Given nonnegative random rates Λ it , i = 1, . . . , m, t = 1, . . . , T , we assume that Y it 's are independent random variables following the generalized Poisson distribution with

and the parameterization (Zamani and Ismail, 2012 )

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

Here, parameter α controls the shape (power) of the relation between the conditional variance of Y it |Λ it and its conditional mean. For example, the relation between (Zamani and Ismail, 2012) .

The underlying random rates Λ 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 Λ 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).

Obviously Eθ it = 1 and

which means that the temporal and spatial correlation structure of the underlying random rates Λ it determine the spatio-temporal correlations between θ it 's. temporal correlation AR(2) ζ i spatial correlation due to distance between regions GMRF ξ i spatial correlation due to neighborhood relation between regions BYM By ignoring these correlations, the standardized incidence ratio θ it = Y it / 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 θ it 's are explained in terms of regional and/or temporal random effects using, for example, a log linear model (Wakefield, 2007; Lee, 2011; Lawson, 2018) .

In the present study, we consider the log linear model

where µ 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 β is its regression coefficient. Moreover, η it is a zero mean random effect which represents spatio-temporal variations in relative risks due to temporal and spatial trend and correlation. Among many different possibilities, we assume that η it takes the additive form

where δ i represents the temporal trend, ε t accounts for temporal correlation and ζ i and ξ i explain spatial correlation due to spatial distance and neighborhood relations among regions A 1 , . . . , A m , respectively (see Table 2 ). The latent (stochastic) temporal trend δ 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) , δ = (δ 1 , . . . , δ 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/τ δ . Here the precision parameter τ δ > 0 acts as a smoothing parameter enforcing small or allowing for large variations in δ t (Fahrmeir and Kneib, 2008) .

To account for temporal correlation, we assume that ε t follows a stationary autoregressive model of order 2, AR(2); i.e., ε t = ψ 1 (1 − ψ 2 )ε t−1 + ψ 2 ε t−2 + t , t = 2, . . . , T, where −1 < ψ 1 < 1 and −1 < ψ 2 < 1 are the first and second partial autocorrelations of ε t and 2 , . . . , T are i.i.d. zero mean Gaussian random variables with variance 1/τ ε .

On the other hand, to account for spatial correlation, we assume that ζ = (ζ 1 , . . . , ζ m ) follows a Gaussian Markov random field (GMRF). More specifically, we assume that ζ is a zero mean Gaussian random vector with the structured covariance matrix

where I m is the m × m identity matrix, 0 ≤ ω < 1 and e max is the largest eigenvalue of the m × m symmetric positive definite matrix C = [C ii ]. The entry C ii of matrix C represents to what extend the regions A i and A i are interconnected. For example, C ii can be related to a data on commuting or population movement between regions A i and A i . In absence of most recent and reliable movement data between the regions of Spain, Italy and Germany, we set C ii to be the Euclidean distance between the centroids of A i and A i . In addition to interconnectivity and correlations due to spatial distance, the neighbourhood structure of regions A 1 , . . . , A m may induce spatial correlation among relative risks of regions because neighbouring regions often tend to have similar relative risks. To include spatial correlation due to neighborhood structure of regions in the model, we assume that ξ = (ξ 1 , . . . , ξ m ) follows a scaled version of the Besag-York-Mollié (BYM) model (Besag et al., 1991) ; i.e., ξ 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 ] with entries

where n i is the number of neighbors of region A i and i ∼ i means that regions A i and A i share a common border. The parameter τ ξ > 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 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 ϕ, and for the parameters of the log linear model for the relative risks µ, β, log τ δ , log τ ε , log τ ζ , log τ ξ , log ω 1−ω , 

1−ψ1 and log 1+ψ2 1−ψ2 . The prior distribution for the α 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 non-informative 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 et al., 2015) . The R-INLA package, an R interface to the INLA program and available at www.r-inla.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 α is chosen to be one (see Table 3 ).

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 cross-validated probability integral transform (PIT)

for calibration checks. 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)

is another Bayesian model diagnostic. Small values of CPO it (y it ) indicate possible outliers, high-leverage and influential observations (Pettit, 1990) .

For count data, Czado et al. (2009) suggested a nonrandomized yet uniform version of the PIT with

which is equivalent to

The mean PIT

can then be comparing with the standard uniform distribution for calibration. For example, a histogram with heights Table 4 presents the Bayesian estimates (posterior means) for every parameter of the considered model 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 ϕ of the generalized Poisson distribution for Italy is higher than Spain and Germany, but its shape parameter α 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 is positive for Germany which indicates that regions with higher population density have larger relative risks.

The precision parameters of the temporal random effects imply that the temporal trend δ t has at least 35 times larger contribution (smaller precision) than ε t which represents temporal correlation. The opposite signs of ψ 1 and ψ 2 indicate rough oscillations in ε t . The spatial random effect ζ i has larger contribution (smaller precision) than ξ i in the total variations of the relative risks only in Spain, while for Italy and Germany it is the opposite. This could be a result of large Euclidean distance between Spain continental European territory from two archipelagos territories, which is affecting the considered covariance structure of ζ i .

In summary, the higher contribution (lower precision) in the total variations of the relative risks for Spain, Italy and Germany is due to the temporal trend, spatial correlation and finally temporal correlation, respectively. This may hint that spatial correlations have a greater impact on the relative risks of COVID-19 than temporal correlations.

The Bayesian estimates and 95% credible intervals for the temporal trend δ t , t = 1, . . . , T , are shown in Figure 1 . These plots can be interpreted as a smoothed temporal trend of the relative risk in the whole country. In fact, Figure 1 suggests that the COVID-19 epidemic in all three countries rapidly reached their peaks and slowly started to decline at the beginning of April and then increased and reached its maximum in August. In addition, the second wave of the epidemic seems to be stronger in Spain and Germany shows a more smoother trend during the study period. Figure 2 shows the the posterior means of the spatial random effects ζ i and ξ 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 each country, particularly in Spain. Regions with positive (negative) ζ i + ξ i values are expected to have elevated (lower) relative risks than the the baseline country-wide risk during the study period.

In order to see how the estimated relative risks under the fitted model are in agreement with the observed data, Figure 3 shows the spatially accumulated daily number of cases m i=1 Y it , t = 1, . . . , T , and their expected values under the fitted model, namely the posterior mean and 95% credible interval of m i=1 E it θ it , t = 1, . . . , T . 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 3 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 Section 4.4 are obtained using J = 20 from the fitted models and plotted in Figure 4 . The normalized PIT values for the fitted models to data do not show a clear visible pattern and the histograms seems to be close to the standard uniform distribution.

The above results and more details on observed and predicated values from the fitted 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 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. However, the posterior predictive checks using the normalized probability integral transform (PIT) showed that there is room for the model improvements.

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. 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 spatio-temporal correlations; e.g., a separable sptaio-temporal covariance structure. However, the considered model was 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.

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. Science of The Total Environment

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

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

Health security capacities in the context of covid-19 outbreak: an analysis of international health regulations annual report data from 182 countries. The Lancet

Spatial epidemic dynamics of the covid-19 outbreak in china

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

Bayesian Disease Mapping: Hierarchical Modeling in Spatial Epidemiology

A comparison of conditional autoregressive models used in bayesian disease mapping. Spatial and spatio-temporal epidemiology

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. Statistical methods in medical research

Approximate bayesian inference for latent gaussian 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 spatio-temporal mapping of disease rates

Coronavirus disease (covid-19) outbreak

Functional form for the generalized poisson regression model

The authors declare that they have no conflict of interest.