A hierarchical spatio-temporal model to analyze relative risk variations of COVID-19: a focus on Spain, Italy and Germany.

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.

Introduction

Started from Wuhan, the capital of Hubei province, China in December 2019, the outbreak of 2019 novel coronavirus disease (COVID-19) 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 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 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.

Data on the daily number of COVID-19 cases

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.

DATADISTA, a Spanish digital communication medium that extracts the data on confirmed cases (PCR and antibody test) registered by the Autonomous Communities of Spain and published by the Ministry of Health and the Carlos III Health Institute. DATADISTA makes the data available in an accessible format at the GitHub repository https://github.com/datadista/datasets/tree/master/COVID The daily accumulated number of total confirmed cases registered in the 19 Autonomous Communities of Spain are updated by DATADISTA on a daily bases.

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.

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.

Table 1

Considered countries with their corresponding number of regions (m), length of study period and the estimated country-wide daily incidence rate ()

Country	Number of regions	Study period	Incidence rate
Spain	18a Autonomous communities	2020-02-25 to 2021-01-16	155.7 × 10−6
Italy	20 Regions	2020-02-25 to 2021-01-16	122.6 × 10−6
Germany	16 Federal states	2020-03-03 to 2021-01-16	79.1 × 10−6

aThe Autonomous Communities Ceuta and Melilla are merged into “Ceuta y Melilla”

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.

Modeling daily regional counts

Suppose that a country, the spatial domain of interest, is partitioned into regions , defined by administrative divisions such as states, provinces, counties, etc (see Table 1). Let denote the number of new COVID-19 cases in region , , at time (day) .

Figure 1 shows boxplots of the observed number of new cases for each region . Clearly, observed values of are heavily right-skewed and directly related to the population of region . Boxplots of the observed values of for each day 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.

Fig. 1

Boxplots of the daily number of COVID-19 cases (left) and population (right) in each region (first level administrative division), , of Spain, Italy and Germany

Fig. 2

Boxplots (yellow boxes with black whiskers) of logarithm of the regional number of COVID-19 cases in each day for Spain, Italy and Germany. The red solid line connects the medians of boxplots and depicts the country-wide temporal trends and variations, and the blue dashed line is the overall median of all cases during the study period

Thus, it is reasonable to assume that the expected number of new COVID-19 cases in region at time t is given bywhere is the population of region and is the incidence rate of COVID-19 in region at time t.

The null model of homogeneous incidence rates

Under the null model of spatial and temporal homogeneity of the incidence rate, we have that , andprovides an estimate for , whereis an estimate of the country-wide homogeneous daily incidence rate (Waller et al. 1997). The estimated daily incidence rate per million population () so far is around 156, 123 and 79 for Spain, Italy and Germany, respectively (see Table 1).

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 , , , following Zamani and Ismail (2012), we writewhere and , and assume that ’s are independent random variables following the generalized Poisson distribution withThen,andThus is the dispersion parameter and the case represents the ordinary Poisson distribution (no dispersion) withHere, parameter controls the shape (power) of the relation between the conditional variance of and its conditional mean. For example, the relation between and is linear if , quadratic if and cubic if (Zamani and Ismail 2012).

Modeling relative risks

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

Relative risks

In disease mapping literature, the nonnegative random quantitiesare called the area-specific relative risks at time t ( Lawson 2018, Section 5.1.4). Note that the random rate and the expected value have the same dimension (unit is people per day) as and hence is a dimensionless random latent variable. Sinceit follows that andwhich means that the temporal and spatial correlation structure of the underlying random rates determine the spatio-temporal correlations between ’s.

By ignoring these correlations, the standardized incidence ratio 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 ’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).

A model for relative risks

In the present study, we consider the log linear modelas the null model, where is the intercept and is the population density of region , i.e. the population of , , divided by the area of . The population density is standardized to have mean 0 and variance 1 and 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 modelswhere represents the temporal trend, accounts for temporal correlation and and explain spatial correlation due to spatial distance and neighborhood relations among regions , respectively (see Table 2). Here, , , and are dimensionless random effects that explain temporal and spatial variability in the relative risk .

Table 2

Considered terms in the additive model for the spatio-temporal random effect of the log-linear model for relative risks

Term	Description	Model
δ_t	Smooth large scale temporal trend	RW2
ε_t	Small scale temporal trends due to temporal correlation	AR(2)
ξ_i	Spatial dependence due to neighborhood relation between regions	BYM
ζ_i	Spatial dependence due to distance between regions	GMRF

The latent (stochastic) temporal trend 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), is assumed to follow a RW2 model, i.e.,where are independent and identically distributed (i.i.d.) zero mean Gaussian random variables with variance . Here the precision parameter acts as a smoothing parameter enforcing small or allowing for large variations in (Fahrmeir and Kneib 2008).

To account for temporal correlation, we assume that follows a stationary autoregressive model of order 2, AR(2); i.e.,where and are the first and second partial autocorrelations of and are i.i.d. zero mean Gaussian random variables with variance .

On the other hand, the neighborhood structure of regions 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 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 denotes the generalized inverse of the spatial precision matrix with entrieswhere is the number of neighbors of region and means that regions and share a common border. The parameter represents the marginal precision and 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 follows a Gaussian Markov random field (GMRF). More specifically, we assume that is a zero mean Gaussian random vector with the structured covariance matrixwhere is the identity matrix, and is the largest eigenvalue of the symmetric positive definite matrix . The entry of matrix represents to what extend the regions and are interconnected. For example, can be related to a data on commuting or population movement between regions and . In absence of most recent and reliable movement data between the regions of Spain, Italy and Germany, we set to be the Euclidean distance between the centroids of and .

Prior specification and implementation

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 is considered as a non-informative prior for the dispersion parameter of generalized Poisson distribution, , and for the parameters of the log linear model for the relative risks , , , , , , , , and . The prior distribution for the parameter of the generalized Poisson distribution is considered to be a Gaussian distribution with mean 1.5 and variance . Table 3 summarizes the model parameters and their necessary transformation for imposing the non-informative Gaussian priors.

Table 3

Parameters of the considered models for the daily number of new cases, their transformation for non-informative uniform (flat) priors and initial values

Parameter	Notation	Transformation	Initial value	Models
Dispersion parameter of generalized Poisson	φ	logφ	0	M_0,M_1,…,M_4
Shape parameter of generalized Poisson	α	α	1.5	M_0,M_1,…,M_4
intercept	μ	μ	0	M_0,M_1,…,M_4
Coefficient of population density	β	β	0	M_0,M_1,…,M_4
Precision (smoothness) of the temporal trend	τ_δ	logτ_δ	0	M_1,…,M_4
Precision of ξ_i	τ_ξ	logτ_ξ	0	M_2,M_3,M_4
Contribution of Q^- in the variance of ξ	ϕ	logϕ1-ϕ	0	M_2,M_3,M_4
Precision of ε_t	τ_ε	logτ_ε	0	M_3,M_4
First partial autocorrelation of ε_t	ψ_1	log1+ψ_11-ψ_1	0	M_3,M_4
Second partial autocorrelation of ε_t	ψ_2	log1+ψ_21-ψ_2	0	M_3,M_4
Precision of ζ_i	τ_ζ	logτ_ζ	0	M_4
Contribution of C in the variance of ζ	ω	logω1-ω	0	M_4

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.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 1.5 (see Table 3).

Bayesian model posterior predictive checks

Let denote the vector of all model parameters andbe the likelihood function for the observed count . The deviance information criterion (DIC)and the Watanabe–Akaike information criterion (WAIC)are two widely used measures of overall model fit, where is the posterior mode (Bayes estimate) of and and denote expectation and variance with respect to the posterior distribution of (Gelman et al. 2014).

For count data and in a Bayesian framework, a probabilistic forecast is a posterior predictive distribution on . 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, letbe the event of observing all count values except the one for region at time t. Dawid (1984) proposed the cross-validated probability integral transform (PIT)for calibration checks, where denotes expectation with respect to the posterior distribution of based on the leave-one-out data . Thus, is simply the value that the predictive distribution function of attains at the observation point . The conditional predictive ordinate (CPO)is another Bayesian model diagnostic. Small values of indicate possible outliers, high-leverage and influential observations (Pettit 1990). Moreover, the Bayesian leave-one-out cross-validationcomputes predictive accuracy (Gelman et al. 2014).

For count data, Czado et al. (2009) suggested a nonrandomized yet uniform version of the PIT withwhich is equivalent towhere is the indicator function, becauseThe mean PITcan then be comparing with the standard uniform distribution for calibration. For example, a histogram with heightsand equally spaced bins , , can be compared with its counterpart from the standard uniform distribution with . 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.

Results

The DIC, WAIC and BCV criteria for the fitted nested models are presented in Table 4. Compared with the null model , it can be seen that including the smooth temporal trend and spatial dependence due to neighborhood structure 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 in the model slightly improves the model fit and prediction accuracy. Finally, including the spatial dependence due to distance between regions in the model leads to slightly worse model fit and prediction accuracy.

Table 4

Deviance information criterion (DIC), Watanabe–Akaike information criterion (WAIC) and Bayesian leave-one-out cross-validation (BCV) for the considered models fitted to the daily number of new COVID-19 cases in Spain, Italy and Germany

Country	Model	DIC	WAIC	BCV
Spain	M_0	76999.0 (–)	76999.9 (–)	− 38499.9 (–)
	M_1	66296.1 (13.90%)	66310.9 (13.88%)	− 33155.7 (13.88%)
	M_2	63372.9 (17.70%)	63388.5 (17.68%)	− 31694.3 (17.68%)
	M_3	63298.8 (17.79%)	63321.4 (17.76%)	− 31660.8 (17.76%)
	M_4	63333.4 (17.75%)	63344.1 (17.73%)	− 31672.1 (17.73%)
Italy	M_0	78704.8 (–)	78705.9 (–)	− 39352.94 (–)
	M_1	66019.4 (16.12%)	66035.4 (16.10%)	− 33017.7 (16.10%)
	M_2	64103.4 (18.55%)	64128.3 (18.52%)	− 32064.2 (18.52%)
	M_3	63796.2 (18.94%)	63838.9 (18.89%)	− 31920.2 (18.89%)
	M_4	63815.5 (18.92%)	63849.5 (18.88%)	− 31924.9 (18.88%)
Germany	M_0	63630.2 (–)	63630.6 (–)	− 31815.3 (–)
	M_1	53579.0 (15.80%)	53583.1 (15.79%)	− 26792.0 (15.79%)
	M_2	51755.7 (18.66%)	51818.8 (18.56%)	− 25910.2 (18.56%)
	M_3	51786.7 (18.61%)	51803.2 (18.59%)	− 25902.2 (18.59%)
	M_4	51809.5 (18.58%)	51812.7 (18.57%)	− 25907.1 (18.57%)

The values in parentheses represent improvement with respect to the null model for each criterion

Numbers represented in bold font show best performing models in temrs of the DIC, WAIC and BCV criteria (smallest values)

Table 5 presents the Bayesian estimates (posterior means) for every parameter of the full 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.

Table 5

Posterior mean and 95% credible interval (in parentheses) for every parameter of the considered model fitted to the daily number of new COVID-19 cases in Spain, Italy and Germany

parameter	Spain	Italy	Germany
φ	0.522 (0.518, 0.526)	1.439 (1.407, 1.468)	1.451 (1.411, 1.486)
α	1.499 (1.495, 1.503)	1.328 (1.322, 1.336)	1.261 (1.255, 1.267)
μ	− 0.859 (− 1.094, − 0.65)	− 1.335 (− 1.41, − 1.271)	− 1.15 (− 1.237, − 1.06)
β	− 0.191 (− 0.461, 0.076)	0.052 (− 0.041, 0.147)	0.236 (0.071, 0.395)
τ_δ	0.058 (0.047, 0.072)	0.112 (0.056, 0.175)	0.102 (0.061, 0.151)
τ_ξ	3.788 (2.23, 6.088)	12.767 (7.011, 22.564)	15.579 (6.69, 30.182)
ϕ	0.490 (0.156, 0.824)	0.818 (0.386, 0.988)	0.738 (0.160, 0.997)
τ_ε	89.296 (74.854, 103.413)	52.822 (34.857, 78.715)	9.96 (7.361, 14.027)
ψ_1	0.668 (0.649, 0.687)	0.592 (0.56, 0.635)	0.587 (0.542, 0.645)
ψ_2	− 0.863 (− 0.899, − 0.801)	− 0.875 (− 0.937, − 0.806)	− 0.756 (− 0.836, − 0.668)
τ_ζ	17372.2 (1199.2, 66901.7)	20531.5 (1084.0, 72227.4)	20100.1 (1249.8, 70790.0)
ω	0.320 (0.001, 0.959)	0.314 (0.001, 0.958)	0.318 (0.001, 0.958)

Comparing the estimated parameters among different countries, it can be seen that the dispersion parameter of the generalized Poisson distribution for Spain is lower than Italy 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 it is positive for Germany which indicates that regions with higher population density have larger relative risks. The opposite signs of and indicate rough oscillations in , 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 byUsing the joint posterior distribution of precision parameters , estimates of the posterior means of , , and 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 . Also, the spatial dependence due to neighborhood relation between regions and small scale temporal trends have minor contributions in the total variations of the relative risks for the three countries. The random effect term representing spatial dependence due to distance between regions has a negligible effect on the relative risks.

Table 6

Contribution (in percent) of each random effect term in the full model on the overall variation of relative risks

country	γ_δ (%)	γ_ξ (%)	γ_ε (%)	γ_ζ (%)
Spain	98.39	1.54	0.06	<0.01
Italy	98.93	0.87	0.20	<0.01
Germany	98.31	0.71	0.98	<0.01

The Bayesian estimates and 95% credible intervals for the temporal trend , , 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.

Fig. 3

Smoothed temporal trend of the relative risks of COVID-19, obtained from posterior mean and 95% credible interval of the structured temporal random effect of the fitted full model

Figure 4 shows the posterior means of the spatial random effects and , , on the corresponding map of each country. The plot illustrates spatial heterogeneity of the relative risk of COVID-19 across regions in each country. Once again, the negligibility of can be seen in Fig. 4. Regions with positive (negative) values are expected to have elevated (lower) relative risks than the baseline country-wide risk during the study period.

Fig. 4

Posterior mean of the spatial random effects representing dependence due to neighborhood relation between regions (left) and representing dependence due to distance between regions (right) in the fitted full model

In order to see how the estimated relative risks under the fitted full model are in agreement with the observed data, Fig. 5 shows the spatially accumulated daily number of cases , , and their expected values under the fitted model , namely the posterior mean and 95% credible interval of , . 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.

Fig. 5

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

Finally, histograms of the normalized PIT values described in Sect. 4.4 are obtained using 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.

Fig. 6

Histograms of normalized PIT values obtained from the fitted full model to check for uniformity

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/.

Concluding remarks

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.

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 in the covariance matrix of the random spatial effect 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 spatio-temporal 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.