A Review of Spatiotemporal Models for Count Data in R Packages. A Case Study of COVID-19 Data

03/08/2021 ∙ by Maria Victoria Ibáñez, et al. ∙ 0

Spatio-temporal models for count data are required in a wide range of scientific fields and they have become particularly crucial nowadays because of their ability to analyse COVID-19-related data. Models for count data are needed when the variable of interest take only non-negative integer values and these integers arise from counting occurrences. Several R-packages are currently available to deal with spatiotemporal areal count data. Each package focuses on different models and/or statistical methodologies. Unfortunately, the results generated by these models are rarely comparable due to differences in notation and methods. The main objective of this paper is to present a review describing the most important approaches that can be used to model and analyse count data when questions of scientific interest concern both their spatial and their temporal behaviour and we monitor their performance under the same data set. For this review, we focus on the three R-packages that can be used for this purpose and the different models assessed are representative of the two most widespread methodologies used to analyse spatiotemporal count data: the classical approach (based on Penalised Likelihood or Estimating Equations) and the Bayesian point of view. A case study is analysed as an illustration of these different methodologies. In this case study, these packages are used to model and predict daily hospitalisations from COVID-19 in 24 health regions within the Valencian Community (Spain), with data corresponding to the period from 28 June to 13 December 2020. Because of the current urgent need for monitoring and predicting data in the COVID-19 pandemic, this case study is, in itself, of particular importance and can be considered the secondary objective of this work. Satisfactory and promising results have been obtained in this second goal.



There are no comments yet.


page 8

page 9

page 20

page 23

This week in AI

Get the week's most popular data science and artificial intelligence research sent straight to your inbox every Saturday.

1 Introduction

Spatial and temporal models for count data are needed in a variety of settings: agricultural production [Besag and Higdon(1999)], fishing catches [Paradinas et al.(2017)Paradinas, Conesa, López-Quílez, and Bellido], volcano eruptions [Gusev(2008)], crime counts [Law et al.(2014)Law, Quick, and Chan], cases of pulmonary disease [Choi et al.(2011)Choi, Lawson, Cai, and Hossain], etc.

In all these examples, the variable of interest take non-negative integer values and these integers arise from counting occurrences of an event in a geographic areal unit in a certain time unit. The observations refer to a set of contiguous non-overlapping areal units for consecutive time periods. Additionally, a series of covariates are measured for each unit of area and time, these covariates could be common in time or area and they can be discrete, continuous or even factors.

The reasons for modelling these data are diverse and can range from estimating the effect of a risk factor to a response, identifying clusters of contiguous areal units or forecasting future observations. Different modelling strategies have been proposed to deal with this array of scenarios comprising spatiotemporal data. The strategies become more complex when the aim is to build multivariate spatiotemporal models for the joint analysis of different variables that include specific and shared spatial and temporal effects [Gómez-Rubio et al.(2019)Gómez-Rubio, Palmí-Perales, López-Abente, Ramis-Prieto, and Fernández-Navarro].

This kind of data presents two major challenges with respect to classical linear regression models. Firstly, it is well known that the normal assumption is not appropriate for count data modelling and generalised linear models with Poisson, binomial or negative binomial distributions must be used

[Nelder and Wedderburn(1972)]. Secondly, spatiotemporal autocorrelation, i.e. that observations from geographically close areal units and temporally close time periods tend to have more similar values than units and time periods that are further apart, result in complicated correlation structures, and as a result, parameter estimation is not straightforward and different approaches have been developed for this purpose [Anselin(1995), Hardisty and Klippel(2010)].

Due to the diversity of applications, data types and conceptual approaches, there is a broad range of literature on spatiotemporal modelling. Two excellent books that provide a gradual entry to the methodological aspects of spatiotemporal statistics and outline some of the standard techniques used in this area are [Cressie and Wikle(2015)] and [Wikle et al.(2019)Wikle, Zammit-Mangion, and Cressie]. An overview of different spatiotemporal modelling approaches can also be found in [Anderson and Ryan(2017)].

Since December 2019, when the first cases of the illness caused by the coronavirus SARS-CoV-2 were reported in Wuhan, China, the SARS-CoV-2 virus has spread world-wide. According to the website of the World Head Organization (WHO) [WHO()], the virus has caused more than 127 million infections and 2.78 million deaths around the world as of 30 March 2021 and it is currently impossible to predict how many people will be affected by it. The illness caused by the SARS-CoV-2 was given the official name ”COVID-19” by the WHO on 11 February 2020 [(WHO)(2020)].

An effective way to control the spread of the infection is to understand and predict key epidemiological data. Epidemic models have provided powerful insight to study data about the coronavirus pandemic, including the number of new cases and deaths in a given area over time. Epidemic models to study the spread of infectious diseases date back to the beginning of the twentieth century. The susceptible-infected-recovered (SIR) models developed by Kermack and McKendrick [Kermack and McKendrick(1927)] were the first mathematical models developed to study the transmission dynamics of infectious diseases. The SIR and SEIR (Susceptible, Exposed, Infectious and Removed) models have been improved and used for analysing and characterising the COVID-19 epidemic [Fang et al.(2020)Fang, Nie, and Penny, Kucharski et al.(2020)Kucharski, Russell, Diamond, Liu, Edmunds, Funk, Eggo, Sun, Jit, Munday, et al., Tang et al.(2020)Tang, Wang, Li, Bragazzi, Tang, Xiao, and Wu], also in Spain [Guirao(2020), López and Rodo(2020)]. Predicting the future course of epidemics has been another great challenge over time and has become particularly challenging with the rise of new infectious diseases [Held and Meyer(2019)].

In addition to the crucial role played by the above described epidemic models, other models are also being adapted to examine different aspects of the COVID-19 pandemic [Fronterre et al.(2020)Fronterre, Read, Rowlingson, Bridgen, Alderton, Diggle, and Jewell, Dunbar and Held(2020)]. In particular, mathematical modelling of patient hospitalisation is essential, as it may help raise awareness of a possible collapse of the health-care systems due to an increase in the number of patients needing hospitalisation. Robust prediction models are therefore vital to support decisions on population and community-level interventions to control the spread of the virus and to prevent the collapse of health services.

Models for this type of data are rather different from epidemic models (SIR and SEIR) because the prediction of hospitalisations requires previously obtained COVID-19 data such as the number of people tested and/or infected or the population at risk. Additionally, models relating to the number of hospitalisations have to manage observations over time in several geographical areas, such as health departments. Each temporal observation relates to an areal unit and, in this case, refers to count measures for the unit: number of COVID-19 hospitalisations in an areal unit per day. Moreover, given the contribution of people’s mobility to the spread of the virus, these models should be able to adjust for people’s mobility between neighbouring health areas.

The purpose of the current study was two-fold. Firstly, to present a review describing three different approaches that can be used to model and analyse count data when the questions of scientific interest concern both their spatial and their temporal behaviour.

In particular, we revisit three different type of models that are representative of the two most widespread methodologies used to analyse this type of data. The first two models are formulated following the classical statistical paradigm and the last one follows the Bayesian point of view (see [Berger(2013)] for a survey of the hierarchical Bayesian approach and [Lawson(2018)] for Bayesian disease mapping). Of the two classical approaches one is based on Penalised Likelihood and the other on Estimating Equations [Liang and Zeger(1986)].

After an exhaustive search among R packages that implement models and procedures to work with spatiotemporal data, we found just three packages that allow us to work with count areal data and we are going to review them. These three packages are: the surveillance package [Meyer et al.(2016)Meyer, Held, and Höhle], the model implemented in this package is based on a likelihood model, working on the classical methodology; the Mcglm package [Bonat(2018)], whose model is based on estimating equations, working on the classical methodology; and finally, the CARBayesST package [Lee et al.(2018)Lee, Rushworth, and Napier], based on the Bayesian methodology.

The main properties and characteristics of each of them will be discussed. This can be useful as a guide for scientists in different experimental fields.

We do not claim to provide a comprehensive coverage of all existing methods to deal with count data but to describe and compare the three relevant approaches that are implemented in R packages [R Core Team(2020)] to be easily used by researchers. Other approaches include Dynamic Spatial Panel Data models [Elhorst(2012), Liesenfeld et al.(2017)Liesenfeld, Richard, and Vogler]

that are more usual in the econometric literature; Machine-Learning techniques such as Classification and Regression Trees, Support Vector Machine and Multilayer perceptron Neural Network

[Martín et al.(2020)Martín, Onrubia, González-Arias, and Vicente-Vírseda]. Generalized Additive Models have also been used in applied real problems with spatiotemporal data [Augustin et al.(1998)Augustin, Borchers, Clarke, Buckland, and Walsh, Beare and Reid(2002), Smith et al.(2019)Smith, Hofner, Lamb, Osenkowski, Allison, Sadoti, McWilliams, and Paton]. With respect to Bayesian models also different types of spatial, temporal and spatiotemporal random effects, not included in the CARBayesST package, can be used such as non-parametric estimation of trends [Knorr-Held(2000)] or splines [Ugarte et al.(2010)Ugarte, Goicoa, and Militino, Bauer et al.(2016)Bauer, Wakefield, Rue, Self, Feng, and Wang].

The second goal of this paper is to apply the reviewed models and compare their performance in the prediction of the number of COVID-19 hospitalisations given the number of infected people in the 24 health departments of the Valencian Community, Spain. The Valencian community is the fourth most populous autonomous community of Spain. It is a rich region with very high residential density along the coast and a lot of tourism and significant exports. Concern about the evolution of the pandemic in this community, led the regional presidency to ask the scientific community for advice about some decision-making regarding the pandemic. Our results will be very useful for this aim.

The article is organised as follows: in Section 2 a descriptive analysis of the data is performed, followed by the mathematical details of the three models. Then, the three models are appplied to COVID-19 data in Section 3. Our results are discussed in Section 4 and finally the conclusions are stated in Section 5.

2 Methodology

2.1 Dataset

The data set comprises the number of daily new positive cases of COVID-19 (tested via PCR (polymerase chain reaction) or the antigen test), and the daily number of hospital admissions due to COVID-19 (daily hospitalisations) in the Valencian Community, a large area of Eastern Spain, over 533 days (from 28 June 2020 to 13 December 2020). This area is organised into different health departments, and the data set comprises both temporal series (hospitalised/new positive cases) for all departments. Figure 1 shows the spatial location of the Valencian Community within Spain and its division into 24 health departments and the distribution of the population ( 100,000 people) in these 24 health departments. As can be observed, the population in the health departments is very heterogeneous.

Figure 1: Location of the Valencian Community within Spain and number of habitants (per 100,000 people) of its 24 health departments.

The number of new daily positive cases by health department is published regularly on an open data platform of the Generalitat Valenciana (Valencian regional goverment) [GVA()]

and the number of daily hospitalisations (by health department) has been provided by the ”Data Science for COVID-19 TaskForce” group of the Valencian Community

[DSA()], with the commitment not to show detailed maps or identifiable information about this variable, which is public only at the aggregated level of the entire Valencian Community. As an illustration, eight of the 24 health departments have been anonymised and will be used to illustrate all the steps in the different analyses. However, the models are fitted using the data of the 24 health departments and the estimations of the parameters and the goodness of fit measures shown in the paper are related to all of them.

Although these data may be subject to temporal biases due to changing testing regimes, among other problems, the mean spatial incidence (number of new cases divided by population size) for three different weeks (2020/07/05-2020/07/11, 2020/09/20-2020/09/26 and 2020/12/04-2020/12/11) plotted in Figure 2, shows strong variation across the different health departments over time.

(a) (b)
Figure 2: Distribution of mean daily incidence (per 100,000 people) in the 24 health departments of the Valencian Community. The mean daily incidence is computed for one week periods: a)2020/07/05-2020/07/11, b)2020/09/20-2020/09/26 and c)2020/12/06-2021/12/11.

Figure 3 a) and b) show time series plots for both variables of interest: the number of daily new positive cases and the number of daily hospitalisations per health department. From now on, daily hospitalisations we will mean the total number of people admitted to hospital due to COVID-19 each day. Figure 3 a) shows a peak of new cases in late August, and another in mid-November, two and a half months apart, and the same temporal pattern appears in the series of hospitalisations.

As different health department have different population sizes, Figures 3 c) and d) show the relative data, it is said, the number of new positive cases and daily hospitalisations corrected by the population size of each health department (incidence values). As can be seen, temporal patters remain unchanged. Figure 3 e) shows the total number of new positive cases and hospitalisations per day (adding the values of the 24 health departments). From Figure 3 e) it can be seen that there is a time lag between both time series.

(a) (b)
(c) (d)
Figure 3:

Temporal trend of COVID-19 for (a) new daily positive cases; (b) daily hospitalisations, per health department, where each point represents the data of one of the 24 health departments;(c) incidence of new daily positive cases, positives divided by region population and multiplied by 100,000; (d)  incidence of daily hospitalisations; e) temporal trend of COVID-19 for the daily sum of the 24 health departments of both time series: daily positive cases (red) and daily hospitalisations (blue) with a mean confidence interval of 95% (grey).

To estimate the delay between both time series a cross-correlogram (see Figure 4) has been used. It plots a measure of correlation of both time series (in this case Pearson correlation), as a function of the displacement (days) of daily positives relative to the daily hospitalisations. Figure 4 shows that there are two maximums in the correlation at the time lag of 9 and 5 days.

Figure 4: Cross-correlogram, between daily new positives and daily hospitalisations. It shows the Pearson correlation between both series as a function of the displacement (days) of daily positives relative to the daily hospitalisations.
a) b)
Figure 5: Data dependency with respect to the day of the week: a)Daily positives vs day of the week; b) Daily hospitalisations vs day of the week; c) Mean of daily positives in the last 4 days vs day of the week.

These data may be subject to temporal biases due to under-reporting at weekends and/or on non-working days. Figure 5 a) shows a great variability in the number of positives depending on the day of the week. As can be seen in Figure 5 b), this effect does not hold for the number of daily hospitalisations. Therefore, this effect is an artefact, due to when the official data is reported rather than a real effect of the virus. To minimise this effect, henceforth, we will take the mean of the last 4 days as the daily hospitalisations. By doing this smoothing, we reduce this effect, as can be seen from Figure 5 c).

Finally, both effects are corrected. Figure 6 shows the smoothed number of positive daily cases (in red) together with the number of people hospitalised due to COVID-19 (in blue) for the eight illustrative health departments, with a time delay of 9 days between both series.

The spatial adjacency matrix between these departments is shown in table 1. In this table, values equal to 1 signal neighbouring regions, i.e. those that share a geographical border.

A 0 1 0 0 0 0 0 0
B 1 0 1 1 0 0 0 0
C 0 1 0 1 1 1 1 0
D 0 1 1 0 1 1 1 1
E 0 0 1 1 0 1 1 1
F 0 0 1 1 1 0 1 1
G 0 0 1 1 1 1 0 0
H 0 0 0 1 1 1 0 0
Table 1: Spatial adjacency matrix of the 8 health departments used for illustration purposes. If two health departments,, are neighbours, the matrix value is one, otherwise 0.
Figure 6: Temporal trend of COVID-19 for daily positives (red line) and daily hospitalisations (blue line) in the eight health departments used as an illustration, with a time lag of 9 days in the daily positive cases, and a smoothing of 4 days.

To conclude, when looking for the best model, we explored the relationship between the mean and the variance of the hospitalisation data collected at each instant of time. Figure 

7 shows that there is a potential relationship between mean and variance. This relationship is crucial for the modelling, as will be seen in the following sections.

Figure 7: Mean vs variance of the number of hospitalisations per day of the 24 health departments for each of the 533 days studied. The data are represented on a log-log scale (black dots). The red line, indicates a smoothing of the data tendency.

2.2 Models

Throughout this section will denote the observation taken in the th areal unit at time , for and . Then, , (), will be a spatiotemporal count series, i.e. count data recorded in the areal units for consecutive discrete time periods. We assume that we also have space-time varying covariates recorded at the same times and locations. Our main objective will be to predict future observations of the spatiotemporal time series , by taking into account the spatiotemporal covariates and the spatial and temporal relationships between the observations.

2.2.1 Endemic-epidemic models. R package surveillance [Meyer et al.(2016)Meyer, Held, and Höhle]

Endemic-epidemic (EE) models are a class of statistical time series models for multivariate surveillance counts proposed by [Held et al.(2005)Held, Höhle, and Hofmann] and extended in [Paul and Held(2011)] and [Meyer et al.(2014)Meyer, Held, et al., Meyer et al.(2016)Meyer, Held, and Höhle]).

In its current formulation and implementation in the R package surveillance [Meyer et al.(2016)Meyer, Held, and Höhle], the EE framework uses incidence from the preceding week, , to explain the incidence in week . So, the counts, , are assumed to be Poisson or Negative Binomial distributed with the conditional mean:


and overdispersion parameter, in the Negative Binomial case, .

The first component of the summation is called the endemic component and captures information not directly linked to observed cases from the previous day. This component can cover exogenous factors such as temporal trends, seasonality, sociodemographics, and/or population. As an example, in spatial applications, can refer to the fraction of the population living in region at time . The remaining terms in Eq. 1 constitute the epidemic component and describe how the incidence in region is linked to previous cases in the same and adjacent regions. The two terms of this epidemic component are usually denoted as ’autoregressive’ and ’spatiotemporal’ component, respectively.

The parameters , and are constrained to be non-negative and can be modelled by allowing for log-linear predictors in all three components, as sine-cosine terms to account for seasonality [Held and Paul(2012)], long-term temporal trends or/and covariates [Bauer and Wakefield(2018), Cheng et al.(2016)Cheng, Lu, Wu, Liu, and Huang].


This form allows for fixed intercepts , region-specific intercepts and exogenous covariates in each model compartment. Population fraction, population density, border effects, etc. can be used as covariates. The region-specific intercepts , can be treated as fixed effects or as random effects accounting for heterogeneity between the regions. When they are treated as random effects, they are assumed to be independent and identically distributed across

, but can be correlated across the model components, following a Gaussian distribution:

We will see this part in more detail in Section 3.1.

Maximum likelihood (ML) estimates are obtained using penalised likelihood approaches.

This basic model has been extended to cover other different aspects of disease modelling (see [Dunbar and Held(2020)] for references). Recent extensions include methodology to adjust for under-reporting [Dunbar and Held(2020), Bracher(2019a)] or to allow different lags in the auto-regressive part of the model (package hhh4addon [Bracher(2019b), Bracher and Held(2017)]), modelling the conditional mean as:


where is the maximum lag considered.

In [Meyer and Held(2017)], the authors extend the basic endemic-epidemic spatiotemporal model to fit multivariate time series of counts stratified by (age) groups in addition to spatial regions. They therefore define a contact matrix , where quantifies the average number of contacts of an individual of group with individuals of group , and the spatiotemporal model is now modelled as:


where both the endemic and epidemic predictors may gain group-specific effects. This model is implemented in the R-package hhh4contacts [Bracher(2019)].


The surveillance package uses the function hhh4 to fit the models and implements the oneStepAhead() function, which computes successive one-step-ahead predictions for the fitted model, also providing confident intervals for the predictions and plot methods. The associated scores-method computes a number of (strictly) proper scoring rules based on such one-step-ahead predictions; see [Paul and Held(2011)] for details.

A discussion of suitable measures to evaluate the quality of a point forecast can be found in [Gneiting(2011)] and several scoring rules based on the one-step-ahead predictions [Paul and Held(2011)] are implemented in the function scores, although we will consider the root mean squared error of the predictions (RMSEp). Another function implemented in the package related to the oneStepAhead() function is the calibrationTest function, which implements calibration tests for Poisson or Negative Binomial predictions of count data based on proper scoring rules; it is described in detail in [Wei and Held(2014)].

Long-term predictions do not have much sense in our context because we do not know the long-term evolution of the covariates.

2.2.2 Multivariate covariance generalised linear models. R package Mcglm

Under the same previous assumption of predicting in terms of spatiotemporal correlations and covariates, we can use the multivariate covariance generalised linear model (McGLM) introduced in [Bonat and Jørgensen(2016)]. This model is a general and flexible statistical model to deal with multivariate count data that explicitly models the marginal covariance matrix combining a covariance link function and a matrix linear predictor composed of known matrices.

Let be the outcome matrix and let denote the corresponding matrix of expected values.

Given the

covariance matrix within the response variable

for and the correlation matrix whose components denote the correlation between outcomes, the McGLM as proposed by [Bonat and Jørgensen(2016)] is given by


where are monotonic differentiable link functions, denotes an design matrix,

is a regression parameter vector to be estimated, and

is the generalised Kronecker product [Martinez-Beneito(2013)]. The matrix denotes the lower triangular matrix of the Cholesky decomposition of . The operator denotes a block diagonal matrix, and denotes a identity matrix.

A key point for specifying McGLMs for mixed types of outcomes is the specification of the covariance matrix within outcomes . Following [Bonat and Jørgensen(2016)], we define the covariance within outcomes by


For modelling count outcomes, they propose to adopt the Poisson-Tweedie dispersion function [Jørgensen and Kokonendji(2016)] so that


is a diagonal matrix whose main entries are given by the power variance function. The Poisson-Tweedie family of distributions provides a rich class of models to deal with count outcomes, since many important distributions appear as special cases; examples include the Hermite ( = 0), Neyman Type A ( = 1), Negative Binomial (p = 2) and Poisson-inverse Gaussian ( = 3).

The dispersion matrix describes the part of the covariance within outcomes that does not depend on the mean structure. Jorgensen et al. [Jørgensen and Kokonendji(2016)], among others, propose to model the dispersion matrix using a matrix linear predictor combined with a covariance link function, i.e.


where is the covariance link function, with are known matrices reflecting the covariance structure within the response variable , and is a vector of dispersion parameters.

McGLMs are fitted based on the estimating function approach described in detail by [Bonat and Jørgensen(2016)]and [Jørgensen and Knudsen(2004)]. A general overview of the algorithm and the asymptotic distribution of the estimating function estimators can be found in [Bonat(2018)]

. As a method for selecting the components of the matrix linear predictor (variable selection), the score information criterion (SIC) is proposed. This is an important tool to assist with the selection of the linear and matrix linear predictor components, but unfortunately it is less useful for comparing models fitted using different link, variance or covariance functions.


Unfortunately, the mcglm package does not have any function implemented to predict future observations. If we do not know to implement a more sophisticated routine, once the model has been estimated, the mc link function can be used. This function returns the inverse of the link function applied to the linear predictor i.e. , as an approximation of the predictions sought.

2.2.3 Bayesian hierarchical generalised linear models. CARst package

A great variety of spatiotemporal models for count data using generalised linear models (GLM) can be found in the Bayesian literature. To model these data a hierarchical model with spatiotemporal structured prior distributions is used. The spatiotemporal structure is modelled via sets of autocorrelated random effects with conditional autoregressive and its spatiotemporal extensions priors. An excellent review can be found in [Lee et al.(2018)Lee, Rushworth, and Napier].

These methods have had a remarkable development, especially in disease mapping, thanks to the availability of estimation methods based on Monte Carlo Markov Chain. With respect to the software packages for implementing these models, although a great quantity of software packages can be found for implementing purely spatial models, such as BUGS

[Lunn et al.(2009)Lunn, Spiegelhalter, Thomas, and Best] and R-INLA [Rue et al.(2009)Rue, Martino, and Chopin], software for spatiotemporal modelling is much less well developed and mainly focuses on geostatistical data. This was the motivation for developing the CARBayesST R package [Lee et al.(2018)Lee, Rushworth, and Napier]. This package can fit several models for count data with different spatiotemporal structures. A useful tutorial is provided by [Lee(2020)]. CARBayesST package has been recently used to study the case-fatality risk by COVID-19 in Colombia [Polo et al.(2020)Polo, Acosta, Soler-Tovar, Villamil, Palencia, Penagos, Martinez, Bobadilla, Martin, Durán, et al.].

The general Bayesian hierarchical model for spatiotemporal count data is as following:


The probability function

is in the exponential family (not necessarily a Gaussian distribution), is the vector of covariate regression parameters, and a multivariate Gaussian prior is assumed. can be any monotonic differentiable link function, and is a latent component for areal unit and time period encompassing one or more sets of spatiotemporally autocorrelated random effects, we denote .

In this paper, we are just focusing on the models implemented in the CARBayesST package. In this package, binomial, Gaussian and Poisson data models can be used for the first level of the model, in Eq. 10, and different spatiotemporal structures for in Eq.11 are given.

All models implemented in this package use random effects to introduce spatial autocorrelation into the response variable. For this purpose, CAR-type prior distributions and their space-time extensions are used. Spatial autocorrelation is induced via a non-negative symmetric matrix of adjacency , where represents the spatial closeness between units . Larger valued elements represent spatial closeness between the two areas in question and spatially autocorrelated random effects, whereas zero values correspond to areas that are not spatially close and conditionally independent random effects given the remaining.

The models are outlined in Table 2 [Lee et al.(2018)Lee, Rushworth, and Napier]. In all cases, inference is based on Markov chain Monte Carlo (MCMC) simulation.

ST.CARlinear [Bernardinelli et al.(1995)Bernardinelli, Clayton, Pascutto, Montomoli, Ghislandi, and Songini] Spatially varying linear time trends
ST.CARanova [Knorr-Held(2000)] Spatial and temporal autoregressive main
effects and independent interaction model
ST.CARsepspatial [Napier et al.(2016)Napier, Lee, Robertson, Lawson, and Pollock] Common temporal trend but varying spatial
surfaces model
ST.CARar [Rushworth et al.(2014)Rushworth, Lee, and Mitchell] Spatially autocorrelated autoregressive of
order 1 time series model
ST.CARadaptive [Rushworth et al.(2017)Rushworth, Lee, and Sarran] Spatially adaptive smoothing model for
localised spatial smoothing
ST.CARlocalised [Lee and Lawson(2016)] spatiotemporal clustering model
Table 2: Summary of the models available in the CARBayesST package.


In order to predict future observations of the response variable, simulations from the posterior predictive distribution (the distribution of possible unobserved values conditional on the observed values) have to be obtained.

This predictive density can be approximated by Monte Carlo integration. If we denote the vector with all the parameters of the model by :

If a representative value is wanted, the mean of the predictive density can be obtained by taking into account the property that and approximating the mean again with respect to :

Unfortunately, the CARBayesST package does not have any function implemented to predict future observations. for =1 it can be easily implemented using the samples of the posterior distributions of the parameters. For the simulation of the posterior distribution of the random effects should be implemented. In this case, for users who are not experts in R programming, just an approximation can be obtained, approximating the value of the random effect by the that of the previous prediction. This approximation will be used in Section 3.3.

3 Application

In this section, all the models reviewed are applied to the COVID-19 data described in Section 2.1. Continuing with the notation introduced, refers to the observed hospitalisation data at time in the health departments of the Valencian Community, refers to the positive cases at time , where is the previously determined and is the population at risk in the -th health department.

Since our primary goal is forecasting, the mean squared errors of the predictions (RMSEp) up to a five-day horizon are calculated. This is a classical approach to measure their performance. Therefore, data from 28 of June 2020 to 8 of December 2020 are used to fit the different models. The root mean square error comparing observed and fitted values (RMSEf) are computed in all cases to describe the goodness the fits. Data from 9 to 13 of December 2020, jointly with the five-days horizon forecasts of each model, are used to compute the root mean square error of the predictions (RMSEp). If we want to use real data of positive cases, the horizon of prediction is limited, but taking a horizon of five is considered enough, within our possibilities, to prevent the collapse of hospitals.

Each model is adjusted using the corresponding statistical methodology included in its corresponding R package. Additional specific measures of goodness of fit are given; these other measures will be useful to compare different models within the same R package.

3.1 Endemic-epidemic models

As stated in Section 2.2.1, let us consider that the counts of daily hospitalisations in the th health department, on the th day,

, follows a Poisson distribution with mean as in Eq.


If denotes the population of the th health department, we assume in Eq. 1 known population fractions

and that , i.e, if both health areas have a common geographic border (assuming the epidemic only arrives from adjacent health areas) and 0 otherwise. Weights are normalised and restricted to be positive.

Covariates such as number of positive cases can be added to the model in different ways [Herzog et al.(2011)Herzog, Paul, and Held, Meyer et al.(2016)Meyer, Held, and Höhle]. The simplest way is to include the covariates in the formulation in the endemic part of the model, for example considering:

  • Model 1:


We are going to consider two possibilities regarding the covariates. Case 1: consider the smoothed number of new positive cases at a lag 9 as a covariate. Case 2: consider the smoothed number of positive new cases at a lag 9 and the smoothed number of new positive cases at a lag 5 as covariates. Many works include a seasonal effect in the model of the parameters and/or , but we would expect this seasonal effect to be included in the covariates (time series of positive cases), so we include it only in the model of and .

All unknown parameters are estimated directly by maximising the corresponding log-likelihoods using numerical optimisation routines (see [Paul et al.(2008)Paul, Held, and Toschke]). The estimates obtained and several goodness of fit measures for this model are shown in Table  3.

Having estimated the parameters of the model, the fitted mean can be compared with the observed counts in order to check the goodness of fit, but additionally, we can see the contribution to this fitted mean of the endemic and autoregressive components. The average of the proportions of the mean explained by the different components are also shown in Table 3. Note that the proportion explained by the epidemic component is around 97 in both cases, being by far the component with the greatest influence on the value of the total fit. So, there is a high influence of the within-health area autoregressive component, with very little contribution of adjacent areas and a rather small endemic incidence.

Model 1.1 Model 1.2
Estimate Std. Error Estimate Std. Error
0.985 0.006 0.985 0.006
0.002 0.0008 0.002 0.0008
4.061 0.709 4.085 0.7117
1.019 0.003 1.005 0.019
- - 1.014 0.0153
Log-likelihood: -8559.39 -8558.74
AIC: 17126.78 17127.48
BIC: 17151.64 17158.56
RMSEf 0.54 0.53
RMSEp 6.13 6.18
endemic 1.39 1.47
epi.own 97.28 97.21
epi.neigbours 1.33 1.32
Table 3: Estimations, goodness of fit measures and contribution of the endemic and autoregressive components to the global fit. RMSEf is the RMSE of fitted values and RMSEp the RMSE of predictions.

We have assumed a Poisson distribution to model the observations but, as seen in Figure 7, there is a clear overdispersion in the data set. To account for this overdispersion, the Poisson distribution may be replaced by two alternatives included in the hhh4 function: ’NegBin1’, that is a Negative Binomial model with a common overdispersion parameter for all areas and ’NegBinM’ that has different overdispersion parameters () for the different health areas, but these distributions do not improve the fit (see Table 4).

Poisson NegBin1 NegBinM
Log-likelihood: -8539
AIC: 17177.99 31554.37 31309.46
BIC: 17488.74 31566.80 31464.84
Table 4: Goodness of fit comparison between Poisson and two Negative Binomial distributions.

This basic model can be refined. Paul et al. [Paul and Held(2011)] introduced random effects for endemic-epidemic models, which are useful if the areas exhibit heterogeneous incidence levels not explained by observed covariates allowing for area-specific intercepts in the endemic or epidemic components. Due to the heterogeneity shown in the different health departments, Model 2 considers an area-specific baseline incidence, ; the population fraction has been included as a multiplicative offset and reflects the mean spatial force of influence of the neighbouring health areas.

Both and have been modelled as fixed effects.

  • Model 2:


So, considering a Poisson distribution for the observations, the goodness of fit measures of Model 2 and the contribution of the three components of the mean to the global fit are shown in Table 5. Once again, results are shown considering the smoothed number of positive new cases at a lag of 9 (Model 2.1) and the smoothed number of new positive cases at lags of 9 and 5 (Model 2.2) as covariates. Individual estimates for the parameters are not shown here, because of the quantity of parameters to estimate.

Model 2.1 Model 2.2
Log-likelihood: -8539 -8538.49
AIC: 17177.99 17178.98
BIC: 17488.74 17495.94
RMSEf 0.71 0.69
RMSEp 5.76 5.78
endemic 1.89 2.04
epi.own 94.73 94.60
epi.neigbours 3.38 3.36
Table 5: Goodness of fit measures and contribution of the endemic and autoregressive components to the global fit. RMSEf is the RMSE of fitted values and RMSEp the RMSE of predictions.

As can be seen in Tables 3 and 5, there is not much difference between the four models. In all of them, the largest portion of the fitted mean results from the within-area autoregressive component (between 94 and 97), with very little contribution of cases from adjacent areas and a rather small endemic incidence. There are also no great differences in the goodness of fit parameters provided by the Likelihood inference (Log-likelihood, Akaike Information Criteria (AIC), Bayesian Information Criteria (BIC)) or in the values of the RMSE.

Figure 8 shows up to 5-day predictions obtained with Model 2.2, together with the fitted and the true observed values.

Figure 8: Observed values (in blue) jointly with fitted values (in brown) and predictions (in green).

3.2 Multivariate Covariance Generalised Linear Models

In this section, we apply the MCGLM approach to analyse the multivariate count data set that was presented in Section 2.1, Eq. 5.

As seen in Eq. 5, MCGLM takes non-normality into account, defining a variance function and modelling the mean structure by means of a link function and a linear predictor. In this application, we have daily observations from 24 health departments () on consecutive days.

Despite the previous model, we are no longer dealing with an autoregressive model, but this autoregressive effect can be included in the linear predictor expression together with other covariates


with being, in this case, the log-link function and the offset.

A key point for specifying McGLMs for mixed types of outcomes is the specification of the covariance matrix within outcomes . Following [Bonat and Jørgensen(2016)], we define the covariance within outcomes by


The matrix linear predictor is defined as: , where denotes the total number of observations in the data set and is the identity matrix. is the intercept of the covariance linear model. If denotes the number of observations in time in each spatial region, , where denotes the Kronecker product and with if and otherwise. measures the effect of the ’time’. Finally, with being the spatial adjacency matrix between the 24 health departments.

In this case, we have obtained the best fits when using the exponential covariance link function as .

To define , we will adopt the Poisson-Tweedie dispersion function [Jørgensen and Kokonendji(2016)] so that


in accordance with the plot in Figure 7.

We employed a step-wise procedure for selecting the components of the linear predictor. As in the previous modelling, we are going to consider the number of new positive cases at a lag of 9, the number of new positive cases at a lag of 5 and the number of hospitalisations at lag of 1 as potential covariates in Eq. 14, and an additional categorical covariate ’Health Department’ to allow different intercepts ( instead of ) in Eq. 14. The population of each health area will also be used as an offset. The SIC using penalty and the Wald test were used in the forward and backward steps, respectively. We defined a stopping criterion for the selection procedure as SIC 0, since the penalty is larger than the score statistics in that case.

The SIC is an important tool to assist with the selection of the linear and matrix linear predictor components, but it is less useful for comparing models fitted using different link, variance or covariance functions. The mcglm package implements the SIC to select the linear and matrix linear predictor components, with the mc_sic and mc_sic_covariance functions. In this case, the SIC values indicate that all components, except the value of the daily hospitalisations with a lag of 1 in neighbouring regions, should be included in the model (SIC 0).

To compare the goodness of non-nested models, the gof function provides the pseudo Akaike information criterion (pAIC), the pseudo Bayesian information criterion (pBIC) and the pseudo Kullback-Leibler information criterion (pKLIC). We are going to use it to select the best structure for the matrix linear predictor. The results can be seen in Table 6. We will fit the model considering .

plogLik -13540.43 -10448.25 -10444.49
df 28 29 30
pAIC 27136.86 20954.5 20948.98
pKLIC 27138.46 20958.28 20954.1
pBIC 27311.06 21134.92 21135.62
RMSEf 13.72 13.62 11.39
RMSEp 16.46 11.77 11.6
Table 6: Goodness of non-nested models, defined from different matrix linear predictors .

Figure 9 shows the up to 5-day predictions obtained with the resulting model, together with the fitted and the true observed values. In this case, we have obtained an RMSEp equal to 11.6, quite higher than the obtained with the other models/packages. Fig. 9 shows that in this case there are health departments where the model provides neither a good fit of the observations nor precise predictions, while the fits and predictions of other health departments are very accurate. With the information available, we have not found a model with better fits and predictions for all the health departments.

Figure 9: Observed values (in blue) together with fitted values (in brown) and predictions (in green).

3.3 Bayesian spatiotemporal models

In this section, we apply different models included in the CARBayesST package to the COVID-19 data given in Section 2.1. For the reasons explained in this section, we use the Poisson log-linear model for , Eq. 10. Overdispersion cannot be controlled with this package. Although that could be regarded as a handicap, as we will see in the results sections, the adjusted and forecasted values are quite good.

In order to induce spatial smoothness between the random effects, we use the binary adjacency matrix used in previous sections. Element =1 if areas and share a common border and =0 otherwise, whereas =0 for all .

Taking into account again the characteristics of our data, the plots shown in Section 2.1 and the fact that our main objective is the prediction of future observations, just two of the spatiotemporal correlation structures included in this package (see Table 2) have been fitted for in Eq. 11: the CARar and CAR adaptative structures. Neither spatially varying linear time trend models nor spatiotemporal clustering models are appropriate for our data and CARanova and CARsepspatial structures assume a symmetric temporal correlation that does not allow us to obtain future predictions.

In both cases, the spatiotemporal structure is modelled with a multivariate first-order autoregressive process with a spatially correlated precision matrix:


is the vector of random effects for time period , the precision corresponds to the CAR models proposed in [Leroux et al.(2000)Leroux, Lei, and Breslow] and has the expression:

is the vector of ones and the identity matrix. respectively control the levels of spatial and temporal autocorrelation, with values of 0 corresponding to independence while a value of 1 corresponds to strong autocorrelation.

The random effects from CARar have a single level of spatial dependence that is controlled by the parameter . That means that all pairs of adjacent areal units will have the same degree of autocorrelation: strongly if is close to one, while no spatial dependence will exist if is close to zero.

The CARadaptative model allows for localised spatial autocorrelation, that is, it allows it to be stronger in some parts of the study region. This could be adequate for our data because it would be possible for spatial autocorrelation between adjacent health departments to be correlated or conditionally independent, depending, for example, on whether these departments are in the same big city or according to the socioeconomic characteristics of their inhabitants.

The CARadaptative model allows this spatial autocorrelation heterogeneity by allowing spatially neighbouring random effects, which is achieved by modelling the non-zero elements of the neighbourhood matrix as unknown parameters rather than assuming they are fixed constants.

With respect to the covariates to be included in the mean, for the same reasons explained in the previous sections, we again tried two possibilities. Case 1: the smoothed number of new positive cases at a lag of 9. Case 2: the smoothed number of new positive cases at a lag of 9 plus the smoothed number of new positive cases at a lag of 5.

Additionally, as a basic area-specific measure of disease incidence, the population fraction has been also included as an offset.

Inference for all models is based on thinning (by 10) 60,000 posterior samples, including a burn-in period of a further 1000 samples. Convergence plots assured that it was reached in all cases.

The Table 7 displays the overall fit of each model by presenting the deviance information criterion (DIC) and the effective number of parameters (pd). It shows that the adaptive model fits the data better than the pure AR model, with reductions in the DIC in both cases. As a complementary measure of goodness of fit and in order to compare across models, the mean square error of adjusted data (RMSEf) has been also calculated and it is showed in Table 7. In this case better results are obtained with the CARar model but as it can be seen the differences are very slight.

case 1 case 2 case 1 case 2
DIC 18730.92 18741.31 18671.93 18678.58
pd 793.5918 787.9688 780.8599 770.0481
RMSEf 1.939859 2.027429 1.970366 2.059419
RMSEp 5.610704 5.854628 5.861029 5.402237
Table 7: Deviance information criterion (DIC), effective number of parameters (pd), RMSE of fitted values (RMSEf) and RMSE of predictions (RMSEp) for each model and scenario.

The medians of the posterior distribution of each parameter and its credible intervals are displayed in Table 8.

As can be seen, the estimated parameters of spatial and temporal correlations show a strong spatial and temporal correlation in all cases. Regarding the covariates, both are significative.

Finally, with respect to the number of step-changes between two spatially adjacent areas detected in the CARadaptative models, only one step change is detected in case 1, while no changes are detected in case 2.

Since our main objective is prediction, we calculate the mean square error of the prediction for up to 5 days for the four cases. These values can be found in the last row of Table 7, and again no great differences are observed. In this case, the best result is obtained with the CARadaptative model with both covariates: positive cases at lags of 9 and 5. Adjusted, observed and predicted values of these models can be seen in Figure 10.

CARar case 1 CARar case 2 CARadaptative case 1 CARadaptative case 2
Median 2.5 % 97.5% Median 2.5 % 97.5% Median 2.5 % 97.5% Median 2.5 % 97.5%
Intercept 5.7011 5.5108 5.7307 5.6233 5.2602 5.6621 5.6955 5.6223 5.7247 5.6107 5.3673 5.6532
Posit9 0.0022 0.0011 0.0108 0.0027 0.0018 0.0100 0.0022 0.0012 0.0054 0.0028 0.0018 0.0081
Posit5 - - - 0.0025 0.0014 0.0097 - - - 0.0026 0.0017 0.0075
tau2 0.0225 0.0197 0.0313 0.0222 0.0193 0.0313 0.0136 0.0099 0.0206 0.0128 0.0094 0.0198
rho.S 0.9432 0.8767 0.9591 0.9666 0.9349 0.9784 0.9652 0.9225 0.9774 0.9671 0.9206 0.9792
rho.T 0.9895 0.9811 0.9949 0.9376 0.8420 0.9548 0.9898 0.9836 0.9952 0.9890 0.9824 0.9946
tau2.w - - - - - - 171.0751 98.3797 289.3351 177.6793 104.1783 293.4227
Table 8: Medians of the posterior distribution of each parameter and 95 % credible intervals for each model and scenario
Figure 10: Fitted and predicted values in brown and green together with the observed counts in blue.

4 Discussion

Having reviewed and applied the different models, it has been seen that the second model (mcglm package) is the most flexible for modelling the variable of interest, , allowing any type of function of time in the expression of the mean and a huge variety of spatiotemporal neighbourhoods to model both the mean and the variance-covariance matrix (when necessary). The first model (surveillance package) can also use any function of time in the expression of the mean, but the temporal relationship with the neighbours only allows an autoregressive (AR) structure of order 1. This model allows us to choose between a Poisson distribution for modelling the observations or a Negative Binomial model when the data show overdispersion. Regarding modelling the mean as a function of , the third model (CARBayesST package) only offers the possibility of linear relationships. In this third model, spatiotemporal relationships between the observations can be introduced into the model of the mean using random effects with different correlation structures.

Having fitted the models, the surveillance package provides the oneStepAhead function to compute successive one-step-ahead predictions for the fitted model, in addition to confident intervals for the predictions and plot methods. However, the other two packages do not have any function implemented to predict future observations from the fitted models, making this process more difficult for non-experts in R programming.

Regarding the estimates of the parameters of the models, significant positive parameters for the covariates are obtained in all cases and the parameters that indicate temporal correlation show high values too for all models. Regarding spatial autocorrelation, there is a difference between the models depending on how the spatial neighbourhood has been included in the model. In general, it can be considered that in order to predict the number of hospitalisations per day, both the number of hospitalisations from the previous day and the number of new cases in the region of interest and in adjacent areas are needed.

Because each type of model has a different methodology that provides different measures of goodness of fit, to compare the performance across different approaches two common and habitual measures are calculated; the RMSE of the predictions up to 5 days, in order to compare prediction performance between different models and the RMSE of adjusted data as a measure of goodness of fit.

With respect to the results, the RMSE of the predictions up to five days of the best model in each package ranges from 5.4 (using the CARBayesST package) to 13.72 (using the mcglm package). There are no large differences between the fitted models using CARBayesST and Surveillance packages, ranging in this case from 5.4 to 5.78, both values are acceptable in clinical practice. RMSE of fitted data are too excellent using CARBayesST and Surveillance (from 0.53 to 1.93), obtaining in this case the smallest error with the Endemic-Epidemic model 3.1 (case 2).

As was previously explained using Multivariate Covariance generalised Linear Models there are health departments where the model provides neither a good fit of the observations nor precise predictions, while the fits and predictions of other health departments are very accurate. With the information available, we have not found a model with better results for all the health departments. But, as has also been said before, this is the most flexible model for modelling spatiotemporal count data response variables and, probably, good results could have been obtained if there had been more information or in other types of applications.

These models can truly be used in the current situation. As far as we know, the surveillance package has been used to address several problems related to the COVD-19 pandemic. In fact, the endemic-epidemic models included in the surveillance package, have been applied to a multitude of infectious diseases (see [Dunbar and Held(2020)] for references) such as Influenza [Paul et al.(2008)Paul, Held, and Toschke], Norovirus [Held et al.(2017)Held, Meyer, and Bracher] and COVID-19 [Fronterre et al.(2020)Fronterre, Read, Rowlingson, Bridgen, Alderton, Diggle, and Jewell, Giuliani et al.(2020)Giuliani, Dickson, Espa, and Santi]. In fact, a regularly updated table of use cases is maintained by S. Meyer at https://github.com/rforge/surveillance/blob/master/www/applications_EE.csv. The CARBayesST package has also been used recently for the study of case-fatality risk due to COVID-19 in Colombia [Polo et al.(2020)Polo, Acosta, Soler-Tovar, Villamil, Palencia, Penagos, Martinez, Bobadilla, Martin, Durán, et al.].

From the beginning of the 2020, other works have been also focused on modeling the number of COVID-19 related hospitalizations, but as far as we know, all them have objectives, covariates and methodologies different from those seen in our work. Ferstad et al. [Ferstad et al.(2020)Ferstad, Gu, Lee, Thapa, Shin, Salomon, Glynn, Shah, Milstein, Schulman, et al.] model the number of people in each county in the United States who are likely to require hospitalization as a result of COVID-19 given the age distribution of the county per the US Census. G. Perone [Perone(2020)] compares several time series forecasting methods to predict the number of patients hospitalized with mild symptoms, and in intensive care units (ICU) in Italy, over the period after October 13, 2020, getting RMSE values greater than ours. Reno et al. [Reno et al.(2020)Reno, Lenzi, Navarra, Barelli, Gori, Lanza, Valentini, Tang, and Fantini] model the spread of COVID-19 and its burden on hospital care under different conditions of social distancing in Lombardy and Emilia-Romagna, the two regions of Italy most affected by the epidemic, using a Susceptible-Exposed-Infectious-Recovered (SEIR) deterministic model, which encompasses compartments relevant to public health interventions such as quarantine. Goic et al. [Goic et al.(2021)Goic, Bozanic-Leal, Badal, and Basso] combine autoregressive, machine learning and epidemiological models to provide a short-term forecast of ICU utilization at the regional level in Chile.

Most of them do not provide a goodness of fit measurement that can be used to compare with our results. Just [Ferstad et al.(2020)Ferstad, Gu, Lee, Thapa, Shin, Salomon, Glynn, Shah, Milstein, Schulman, et al.] and [Goic et al.(2021)Goic, Bozanic-Leal, Badal, and Basso] give values of RMSE. These values are in general greater than the obtained in our work, but they are not directly comparable, because their objective is to predict the number of patients hospitalized with mild symptoms, and/or in intensive care units (ICU) separately.

Within the framework of the government support group of the Generalitat Valenciana our models are intended to be used and updated weekly, helping the government make public health decisions such as, the possible need to open new COVID-19 wards in hospitals in the most affected regions.

5 Conclusions

The aim of this work is to review three different spatiotemporal models for count data, implemented in R packages, and to test their performance on an actual case study using three completely different approaches.

Due to these different statistical methodologies, the different packages provide different goodness of fit measures, there being no measure in common between them. Therefore, as the final aim of our case study has been the short-term prediction of the evolution of hospitalisations in the different spatial areas, the root mean squared prediction error (RMSE) has been obtained in all cases. We can achieve very satsifactory results using each of the packages reviewed and there is not much difference between them. These results are very promising in the particular case of the Valencian Community, but they are also very valuable because they can be applied to any region and the use of these models can be promoted to help in short-term government decision making regarding preventive measures against the collapse of hospitals. Additionally, they can provide tools to know in advance whether it is necessary to expand hospital capacity in terms of beds and/or workers.

Our objective in this work has not been to say that one package or one type of model is better than others, but to show possibilities that can be used in practice to analyse this type of data. The choice of the type of model to use will depend on the application at hand.


This research is supported by Ayudas Fundación BBVA a Equipos de Investigación Científica SARS-CoV-2 y COVID-19 and FONDO SUPERA COVID19 by Banco Santander. We would like to thank to Dr. Nuria Oliver (chief head of Generalitat Valenciana Task force) and Generalitat Valenciana for their support providing the data set.


  • [1]
  • [Besag and Higdon(1999)] Besag, J.; Higdon, D. Bayesian analysis of agricultural field experiments. Journal of the Royal Statistical Society: Series B (Statistical Methodology) 1999, 61, 691–746.
  • [Paradinas et al.(2017)Paradinas, Conesa, López-Quílez, and Bellido] Paradinas, I.; Conesa, D.; López-Quílez, A.; Bellido, J.M. Spatio-temporal model structures with shared components for semi-continuous species distribution modelling. Spatial Statistics 2017, 22, 434–450.
  • [Gusev(2008)] Gusev, A. Temporal structure of the global sequence of volcanic eruptions: Order clustering and intermittent discharge rate. Physics of the Earth and planetary interiors 2008, 166, 203–218.
  • [Law et al.(2014)Law, Quick, and Chan] Law, J.; Quick, M.; Chan, P. Bayesian spatio-temporal modeling for analysing local patterns of crime over time at the small-area level. Journal of quantitative criminology 2014, 30, 57–78.
  • [Choi et al.(2011)Choi, Lawson, Cai, and Hossain] Choi, J.; Lawson, A.B.; Cai, B.; Hossain, M.M. Evaluation of Bayesian spatiotemporal latent models in small area health data. Environmetrics 2011, 22, 1008–1022.
  • [Gómez-Rubio et al.(2019)Gómez-Rubio, Palmí-Perales, López-Abente, Ramis-Prieto, and Fernández-Navarro] Gómez-Rubio, V.; Palmí-Perales, F.; López-Abente, G.; Ramis-Prieto, R.; Fernández-Navarro, P. Bayesian joint spatio-temporal analysis of multiple diseases. SORT-Statistics and Operations Research Transactions 2019, pp. 51–74.
  • [Nelder and Wedderburn(1972)] Nelder, J.A.; Wedderburn, R.W. Generalized linear models. Journal of the Royal Statistical Society: Series A (General) 1972, 135, 370–384.
  • [Anselin(1995)] Anselin, L. Local indicators of spatial association—LISA. Geographical analysis 1995, 27, 93–115.
  • [Hardisty and Klippel(2010)] Hardisty, F.; Klippel, A. Analysing spatio-temporal autocorrelation with LISTA-Viz. International Journal of Geographical Information Science 2010, 24, 1515–1526.
  • [Cressie and Wikle(2015)] Cressie, N.; Wikle, C.K. Statistics for spatio-temporal data; John Wiley & Sons, 2015.
  • [Wikle et al.(2019)Wikle, Zammit-Mangion, and Cressie] Wikle, C.K.; Zammit-Mangion, A.; Cressie, N. Spatio-temporal Statistics with R; CRC Press, 2019.
  • [Anderson and Ryan(2017)] Anderson, C.; Ryan, L.M. A comparison of spatio-temporal disease mapping approaches including an application to ischaemic heart disease in New South Wales, Australia. International journal of environmental research and public health 2017, 14, 146.
  • [WHO()] World Health Organization Coronavirus (COVID-19) Dashboard. https://covid19.who.int/, year=2021,.
  • [(WHO)(2020)] (WHO), W.H.O. Naming the coronavirus disease (COVID-19) and the virus that causes it. https://www.who.int/emergencies/diseases/novel-coronavirus-2019/technical-guidance/naming-the-coronavirus-disease-(covid-2019)-and-the-virus-that-causes-it, 2020.
  • [Kermack and McKendrick(1927)] Kermack, W.O.; McKendrick, A.G. A contribution to the mathematical theory of epidemics. Proceedings of the royal society of london. Series A, Containing papers of a mathematical and physical character 1927, 115, 700–721.
  • [Fang et al.(2020)Fang, Nie, and Penny] Fang, Y.; Nie, Y.; Penny, M. Transmission dynamics of the COVID-19 outbreak and effectiveness of government interventions: A data-driven analysis. Journal of medical virology 2020, 92, 645–659.
  • [Kucharski et al.(2020)Kucharski, Russell, Diamond, Liu, Edmunds, Funk, Eggo, Sun, Jit, Munday, et al.] Kucharski, A.J.; Russell, T.W.; Diamond, C.; Liu, Y.; Edmunds, J.; Funk, S.; Eggo, R.M.; Sun, F.; Jit, M.; Munday, J.D.; others. Early dynamics of transmission and control of COVID-19: a mathematical modelling study. The lancet infectious diseases 2020.
  • [Tang et al.(2020)Tang, Wang, Li, Bragazzi, Tang, Xiao, and Wu] Tang, B.; Wang, X.; Li, Q.; Bragazzi, N.L.; Tang, S.; Xiao, Y.; Wu, J. Estimation of the transmission risk of the 2019-nCoV and its implication for public health interventions. Journal of clinical medicine 2020, 9, 462.
  • [Guirao(2020)] Guirao, A. The Covid-19 outbreak in Spain. A simple dynamics model, some lessons, and a theoretical framework for control response. Infectious Disease Modelling 2020, 5, 652–669.
  • [López and Rodo(2020)] López, L.; Rodo, X. A modified SEIR model to predict the COVID-19 outbreak in Spain and Italy: simulating control scenarios and multi-scale epidemics. Available at SSRN 3576802 2020.
  • [Held and Meyer(2019)] Held, L.; Meyer, S. Forecasting based on surveillance data. Handbook of Infectious Disease Data Analysis. Chapman & Hall/CRC 2019.
  • [Fronterre et al.(2020)Fronterre, Read, Rowlingson, Bridgen, Alderton, Diggle, and Jewell] Fronterre, C.; Read, J.M.; Rowlingson, B.; Bridgen, J.; Alderton, S.; Diggle, P.J.; Jewell, C.P. COVID-19 in England: spatial patterns and regional outbreaks. medRxiv 2020.
  • [Dunbar and Held(2020)] Dunbar, M.B.N.; Held, L. Endemic-Epidemic framework used in COVID-19 modelling. REVSTAT–Statistical Journal 2020, 18, 565–574.
  • [Berger(2013)] Berger, J.O. Statistical decision theory and Bayesian analysis; Springer Science & Business Media, 2013.
  • [Lawson(2018)] Lawson, A.B. Bayesian disease mapping: hierarchical modeling in spatial epidemiology; CRC press, 2018.
  • [Liang and Zeger(1986)] Liang, K.Y.; Zeger, S.L. Longitudinal data analysis using generalized linear models. Biometrika 1986, 73, 13–22.
  • [Meyer et al.(2016)Meyer, Held, and Höhle] Meyer, S.; Held, L.; Höhle, M. hhh4: Endemic-epidemic modeling of areal count time series. Journal of Statistical Software 2016.
  • [Bonat(2018)] Bonat, W.H. Multiple response variables regression models in R: The mcglm package. Journal of Statistical Software 2018, 84.
  • [Lee et al.(2018)Lee, Rushworth, and Napier] Lee, D.; Rushworth, A.; Napier, G. Spatio-temporal areal unit modelling in R with conditional autoregressive priors using the CARBayesST package. Journal of Statistical Software 2018, 84.
  • [R Core Team(2020)] R Core Team. R: A Language and Environment for Statistical Computing. R Foundation for Statistical Computing, Vienna, Austria, 2020.
  • [Elhorst(2012)] Elhorst, J.P. Dynamic spatial panels: models, methods, and inferences. Journal of geographical systems 2012, 14, 5–28.
  • [Liesenfeld et al.(2017)Liesenfeld, Richard, and Vogler] Liesenfeld, R.; Richard, J.F.; Vogler, J. Likelihood-Based Inference and Prediction in Spatio-Temporal Panel Count Models for Urban Crimes. Journal of applied econometrics 2017, 32, 600–620.
  • [Martín et al.(2020)Martín, Onrubia, González-Arias, and Vicente-Vírseda] Martín, B.; Onrubia, A.; González-Arias, J.; Vicente-Vírseda, J.A. Citizen science for predicting spatio-temporal patterns in seabird abundance during migration. PloS one 2020, 15, e0236631.
  • [Augustin et al.(1998)Augustin, Borchers, Clarke, Buckland, and Walsh] Augustin, N.H.; Borchers, D.L.; Clarke, E.D.; Buckland, S.T.; Walsh, M. Spatiotemporal modelling for the annual egg production method of stock assessment using generalized additive models. Canadian Journal of Fisheries and Aquatic Sciences 1998, 55, 2608–2621.
  • [Beare and Reid(2002)] Beare, D.J.; Reid, D.G. Investigating spatio-temporal change in spawning activity by Atlantic mackerel between 1977 and 1998 using generalized additive models. ICES Journal of Marine Science 2002, 59, 711–724.
  • [Smith et al.(2019)Smith, Hofner, Lamb, Osenkowski, Allison, Sadoti, McWilliams, and Paton] Smith, A.; Hofner, B.; Lamb, J.S.; Osenkowski, J.; Allison, T.; Sadoti, G.; McWilliams, S.R.; Paton, P. Modeling spatiotemporal abundance of mobile wildlife in highly variable environments using boosted GAMLSS hurdle models. Ecology and evolution 2019, 9, 2346–2364.
  • [Knorr-Held(2000)] Knorr-Held, L. Bayesian modelling of inseparable space-time variation in disease risk. Statistics in medicine 2000, 19, 2555–2567.
  • [Ugarte et al.(2010)Ugarte, Goicoa, and Militino] Ugarte, M.; Goicoa, T.; Militino, A. Spatio-temporal modeling of mortality risks using penalized splines. Environmetrics: The official journal of the International Environmetrics Society 2010, 21, 270–289.
  • [Bauer et al.(2016)Bauer, Wakefield, Rue, Self, Feng, and Wang] Bauer, C.; Wakefield, J.; Rue, H.; Self, S.; Feng, Z.; Wang, Y. Bayesian penalized spline models for the analysis of spatio-temporal count data. Statistics in medicine 2016, 35, 1848–1865.
  • [GVA()] GVA. Portal de dades obertes. https://dadesobertes.gva.es/es/dataset.
  • [DSA()] Data Science against COVID-19 Taskforce. https://ellisalicante.org/ai4covid19/en, year=2020,.
  • [Held et al.(2005)Held, Höhle, and Hofmann] Held, L.; Höhle, M.; Hofmann, M. A statistical framework for the analysis of multivariate infectious disease surveillance counts. Statistical modelling 2005, 5, 187–199.
  • [Paul and Held(2011)] Paul, M.; Held, L. Predictive assessment of a non-linear random effects model for multivariate time series of infectious disease counts. Statistics in medicine 2011, 30, 1118–1136.
  • [Meyer et al.(2014)Meyer, Held, et al.] Meyer, S.; Held, L.; others. Power-law models for infectious disease spread. The Annals of Applied Statistics 2014, 8, 1612–1639.
  • [Held and Paul(2012)] Held, L.; Paul, M. Modeling seasonality in space-time infectious disease surveillance data. Biometrical Journal 2012, 54, 824–843.
  • [Bauer and Wakefield(2018)] Bauer, C.; Wakefield, J. Stratified space–time infectious disease modelling, with an application to hand, foot and mouth disease in China. Journal of the Royal Statistical Society: Series C (Applied Statistics) 2018, 67, 1379–1398.
  • [Cheng et al.(2016)Cheng, Lu, Wu, Liu, and Huang] Cheng, Q.; Lu, X.; Wu, J.T.; Liu, Z.; Huang, J. Analysis of heterogeneous dengue transmission in Guangdong in 2014 with multivariate time series model. Scientific reports 2016, 6, 33755.
  • [Bracher(2019a)] Bracher, J. hhh4underreporting. R package. https://github.com/jbracher/hhh4underreporting, 2019.
  • [Bracher(2019b)] Bracher, J. hhh4addon: extending the funcionality of surveillance: hhh4. R package. https://github.com/jbracher/hhh4addon, 2019.
  • [Bracher and Held(2017)] Bracher, J.; Held, L. Periodically stationary multivariate autoregressive models. arXiv preprint arXiv:1707.04635 2017.
  • [Meyer and Held(2017)] Meyer, S.; Held, L. Incorporating social contact data in spatio-temporal models for infectious disease spread. Biostatistics 2017, 18, 338–351.
  • [Bracher(2019)] Bracher, J. hhh4contacts: Age-structured spatio-temporal models for infectious disease counts. R package. https://github.com/cran/hhh4contacts, 2019.
  • [Gneiting(2011)] Gneiting, T. Making and evaluating point forecasts. Journal of the American Statistical Association 2011, 106, 746–762.
  • [Wei and Held(2014)] Wei, W.; Held, L. Calibration tests for count data. Test 2014, 23, 787–805.
  • [Bonat and Jørgensen(2016)] Bonat, W.; Jørgensen, B. Multivariate covariance generalized linear models. Journal of the Royal Statistical Society, Series C (applied Statistics) 2016, 65, 649–675.
  • [Martinez-Beneito(2013)] Martinez-Beneito, M.A. A general modelling framework for multivariate disease mapping. Biometrika 2013, 100, 539–553.
  • [Jørgensen and Kokonendji(2016)] Jørgensen, B.; Kokonendji, C.C. Discrete dispersion models and their Tweedie asymptotics. AStA Advances in Statistical Analysis 2016, 100, 43–78.
  • [Jørgensen and Knudsen(2004)] Jørgensen, B.; Knudsen, S.J. Parameter orthogonality and bias adjustment for estimating functions. Scandinavian Journal of Statistics 2004, 31, 93–114.
  • [Lunn et al.(2009)Lunn, Spiegelhalter, Thomas, and Best] Lunn, D.; Spiegelhalter, D.; Thomas, A.; Best, N. The BUGS project: Evolution, critique and future directions. Statistics in medicine 2009, 28, 3049–3067.
  • [Rue et al.(2009)Rue, Martino, and Chopin] Rue, H.; Martino, S.; Chopin, N. Approximate Bayesian inference for latent Gaussian models by using integrated nested Laplace approximations. Journal of the royal statistical society: Series b (statistical methodology) 2009, 71, 319–392.
  • [Lee(2020)] Lee, D. A tutorial on spatio-temporal disease risk modelling in R using Markov chain Monte Carlo simulation and the CARBayesST package. Spatial and Spatio-temporal Epidemiology 2020, 34, 100353.
  • [Polo et al.(2020)Polo, Acosta, Soler-Tovar, Villamil, Palencia, Penagos, Martinez, Bobadilla, Martin, Durán, et al.] Polo, G.; Acosta, C.M.; Soler-Tovar, D.; Villamil, J.F.P.; Palencia, N.P.; Penagos, M.; Martinez, J.M.; Bobadilla, J.N.; Martin, L.V.; Durán, S.; others. Bayesian Spatio-Temporal Modeling of COVID-19: Inequalities on Case-Fatality Risk. medRxiv 2020.
  • [Bernardinelli et al.(1995)Bernardinelli, Clayton, Pascutto, Montomoli, Ghislandi, and Songini] Bernardinelli, L.; Clayton, D.; Pascutto, C.; Montomoli, C.; Ghislandi, M.; Songini, M. Bayesian analysis of space–time variation in disease risk. Statistics in medicine 1995, 14, 2433–2443.
  • [Napier et al.(2016)Napier, Lee, Robertson, Lawson, and Pollock] Napier, G.; Lee, D.; Robertson, C.; Lawson, A.; Pollock, K.G. A model to estimate the impact of changes in MMR vaccine uptake on inequalities in measles susceptibility in Scotland. Statistical methods in medical research 2016, 25, 1185–1200.
  • [Rushworth et al.(2014)Rushworth, Lee, and Mitchell] Rushworth, A.; Lee, D.; Mitchell, R. A spatio-temporal model for estimating the long-term effects of air pollution on respiratory hospital admissions in Greater London. Spatial and spatio-temporal epidemiology 2014, 10, 29–38.
  • [Rushworth et al.(2017)Rushworth, Lee, and Sarran] Rushworth, A.; Lee, D.; Sarran, C. An adaptive spatiotemporal smoothing model for estimating trends and step changes in disease risk. Journal of the Royal Statistical Society: Series C (Applied Statistics) 2017, 66, 141–157.
  • [Lee and Lawson(2016)] Lee, D.; Lawson, A. Quantifying the spatial inequality and temporal trends in maternal smoking rates in Glasgow. The annals of applied statistics 2016, 10, 1427.
  • [Herzog et al.(2011)Herzog, Paul, and Held] Herzog, S.; Paul, M.; Held, L. Heterogeneity in vaccination coverage explains the size and occurrence of measles epidemics in German surveillance data. Epidemiology & Infection 2011, 139, 505–515.
  • [Paul et al.(2008)Paul, Held, and Toschke] Paul, M.; Held, L.; Toschke, A.M. Multivariate modelling of infectious disease surveillance data. Statistics in medicine 2008, 27, 6250–6267.
  • [Leroux et al.(2000)Leroux, Lei, and Breslow] Leroux, B.G.; Lei, X.; Breslow, N. Estimation of disease rates in small areas: a new mixed model for spatial dependence. In Statistical models in epidemiology, the environment, and clinical trials; Springer, 2000; pp. 179–191.
  • [Held et al.(2017)Held, Meyer, and Bracher] Held, L.; Meyer, S.; Bracher, J. Probabilistic forecasting in infectious disease epidemiology: the 13th Armitage lecture. Statistics in medicine 2017, 36, 3443–3460.
  • [Giuliani et al.(2020)Giuliani, Dickson, Espa, and Santi] Giuliani, D.; Dickson, M.M.; Espa, G.; Santi, F. Modelling and predicting the spatio-temporal spread of Coronavirus disease 2019 (COVID-19) in Italy. Available at SSRN 3559569 2020.
  • [Ferstad et al.(2020)Ferstad, Gu, Lee, Thapa, Shin, Salomon, Glynn, Shah, Milstein, Schulman, et al.] Ferstad, J.O.; Gu, A.J.; Lee, R.Y.; Thapa, I.; Shin, A.Y.; Salomon, J.A.; Glynn, P.; Shah, N.H.; Milstein, A.; Schulman, K.; others. A model to forecast regional demand for COVID-19 related hospital beds. medRxiv 2020.
  • [Perone(2020)] Perone, G. Comparison of ARIMA, ETS, NNAR and hybrid models to forecast the second wave of COVID-19 hospitalizations in Italy. arXiv preprint arXiv:2010.11617 2020.
  • [Reno et al.(2020)Reno, Lenzi, Navarra, Barelli, Gori, Lanza, Valentini, Tang, and Fantini] Reno, C.; Lenzi, J.; Navarra, A.; Barelli, E.; Gori, D.; Lanza, A.; Valentini, R.; Tang, B.; Fantini, M.P. Forecasting COVID-19-associated hospitalizations under different levels of social distancing in Lombardy and Emilia-Romagna, Northern Italy: results from an extended SEIR compartmental model. Journal of clinical medicine 2020, 9, 1492.
  • [Goic et al.(2021)Goic, Bozanic-Leal, Badal, and Basso] Goic, M.; Bozanic-Leal, M.S.; Badal, M.; Basso, L.J. COVID-19: Short-term forecast of ICU beds in times of crisis. PloS one 2021, 16, e0245272.