Spatial Modeling in Epidemiology

Abstract

The objective of this chapter is to present the methodology of some of the models used in the area of epidemiology, which are used to study, understand, model and predict diseases (infectious and non-infectious) occurring in a given region. These models, which belong to the area of geostatistics, are usually composed of a fixed part and a random part. The fixed part includes the explanatory variables of the model and the random part includes, in addition to the error term, a random term that generally has a multivariate Gaussian distribution. Based on the random effect, the spatial correlation (or covariance) structure of the data will be explained. In this way, the spatial variability of the data in the region of interest is accounted for, thus avoiding that this information is added to the model error term. The chapter begins by introducing Gaussian processes, and then looks at their inclusion in generalized spatial linear models, spatial survival analysis and finally in the generalized extreme value distribution for spatial data. The review also mentions some of the main packages that exist in the R statistical software and that help with the implementation of the mentioned spatial models.

1. Introduction

The term spatial statistics is used to describe a wide range of statistical models and methods for the analysis of geo-referenced data [1]. Its rapid use has been increasing in various fields of science, such as biology, image processing, environmental and earth sciences, ecology, epidemiology, agronomy, forestry, among others [2]. In epidemiology, spatial statistics are used to study the occurrence of health-disease events or deaths in a region of interest. It is now known that several public health problems tend to exhibit spatial dependence (spatial autocorrelation, spatial variability), and that sometimes these problems are related to climatic factors that are generally of a spatially continuous nature or with factors specific to the study region. The use of classical statistical techniques to model spatial data generally leads to an overestimation of model parameters [1]; and although they may eventually help, these models, lacking adequate structure, will not be able to model the spatial variability of the data; valuable information that will be sent to model error and cannot be used to explain the nature of the phenomenon under study.

Recent studies have shown that spatial models can help identify spatial patterns in infectious and non-infectious diseases. These models also help determine the factors that favor them, such as sociodemographic, environmental, etc.; as well as generate maps to visualize the distribution of morbidity or mortality of infectious and non-infectious diseases, and identify critical points in the spatial distribution [3, 4].

Generalized linear spatial models (GLSM), which are a particular class of multilevel or hierarchical models, have been used for the study of certain diseases (infectious and non-infectious). The estimation of GLSM parameters can be done under the frequentist or Bayesian approach [1], some examples are given below. A spatial Poisson regression model, where parameter estimation was performed under the frequentist approach, was used to study esophageal cancer incidence rates [5] and the sociodemographic risk factors for diabetes [6]. Under the Bayesian approach, these models have been used to study the relationship between Visceral Leishmaniasis incidence rates and climatological variables [7], as well as to identify risk factors associated with nontuberculous mycobacterial infections [8]. Spatial Binomial regression models, under the Bayesian approach, have been used to describe patterns of occurrence of dengue and chikungunya [9], and filariasis [10]. Under the classical approach, spatial binomial regression models have been used to investigate environmental and sociodemographic factors associated with leptoserosis disease [11]; are also used to study risk factors associated with HIV infection among drug users [12].

On the other hand, survival analysis under the spatial approach has also received great attention in recent years, because geographic location can play a relevant role in predicting disease survival [13]. Fragility models (spatial survival models) can be an option to analyze the heterogeneity of the data when it cannot be explained by the covariates in a classical survival model. In spatial survival models, in addition to covariates, a random effect known as frailty is added, which modifies the hazard function of an individual, or of spatially correlated individuals [14]. Generally, the random factor, which is assigned a multivariate normal distribution, plays an important role in modeling survival times; since in this term the differences that exist in the socioeconomic level, access to medical care, population density, weather conditions, among others, can be taken into account. It is worth mentioning that spatial survival models have been applied in studies such as: recovery time in patients with COVID-19 [15], hospitalization time in dengue patients [16], HIV/AIDS survival [17] and breast cancer [18] to name a few. In all these works, the estimation of the model parameters was under the Bayesian approach.

Extreme events in public health (for example, the saturation of hospitals) are generally analyzed through measures of central tendency or time series, however, these approaches are not the most appropriate to understand extreme events (unusual events); that when they occur they strongly impact the health care network, thus often collapsing the system [19]. The extreme value theory (EVT) aims to study the probability of occurrence of extreme events (values) of a phenomenon of interest over time, generally these values only occur when they exceed a threshold. Although the applications of EVT in public health are scarce, if they exist at all; an application was presented when predicting extreme events of annual seasonal influenza mortality and the number of emergency department visits in a network of hospitals [20], another application was presented when modeling elevated cholesterol levels using the spikes-over-threshold model [21]. In both cases, the parameters were estimated under the frequentist approach. Given the advantages they have with the application of a spatial model, it would be convenient to study the extreme events of the health sector in space, for which there is already a methodology known as spatial modeling of extreme values [22].

The objective of this work is to provide a general review of the theoretical framework of spatial statistical models developed in the area of geostatistics, which have been used in the area of epidemiology to analyze, model and predict the phenomena of interest. Some of the packages that exist in the statistical software R [23] to carry out said spatial analyzes are also mentioned.

2. Gaussian processes

A stochastic process $W\left(t\right),t\in T$ is a collection of random variables. That is, for each $t\in T$, $W\left(t\right)$ is a random variable [24]; if the stochastic process is indexed by a coordinate space $s\in A\subset {\mathbb{R}}^d$, then the stochastic process is called a random field [25]. A realization of the random field, $W\left(s\right),s\in A$, is given by $\left(W\left(s_1\right)=y\left(s_1\right),\dots ,W\left(s_n\right)=y\left(s_n\right)\right)$. Generally from the sample $y\left(s_1\right),\dots ,y\left(s_n\right)$ one tries to know the characteristics of the process $W$ in $s_i$, $i=1,\dots ,n$; and with this information to make inference of the process $W\left(s\right)$ on all $A\subset {\mathbb{R}}^d$, a convex set where $s$ varies continuously. To the geo-referenced data $y\left(s_1\right),\dots ,y\left(s_n\right)$ is often referred to as geocoded, geostatistical data or point-referenced data. The study of this type of data is known as geostatistics, which is a part of spatial statistics that studies phenomena with continuous variation in space, a convex region denoted $A$ [26].

A process $W$ is second order stationary if it has finite variance, constant mean and its covariance function depends only on distance. Having second-order stationarity in a stochastic process implies having intrinsic stationarity, i.e., second-order stationarity is stronger than intrinsic stationarity. On the other hand, weak stationarity and second-order stationarity are equivalent in the space [27]. The following defines what is known as a Gaussian process (field).

Definition 1. A stochastic process $\left\{W\left(s\right):s\in A\subset {\mathbb{R}}^2\right\}$, where $s$ varies continuously on a fixed subset $A$ content in ${\mathbb{R}}^2$, is a Gaussian process if for any collection of locations $s_1,\dots ,s_n$ with $s_i\in A$, the joint distribution of $\left(W\left(s_1\right),\dots ,W\left(s_n\right)\right)$ is multivariate Gaussian [1].

What is known as a stationary Gaussian process is defined below.

Definition 2. A Gaussian process $\left\{W\left(s\right):s\in A\subset {\mathbb{R}}^2\right\}$, is stationary if $\forall s\in A$:

$$E\left(W\left(s\right)\right)=0,$$E1

$$\mathit{Var}\left(W\left(s\right)\right)={\sigma }^2,$$E2

and its correlation function depends only on the distance, i.e.

$$\mathit{\text{Corr}}\left(W\left(s\right),W\left(s^{\prime }\right)\right)=\rho \left(h\right),$$E3

where $h=‖s-s^{\prime }‖$ is the Euclidean distance that exists between $s$ and $s^{\prime }$.

That is, the mean and variance of $W\left(s\right)$ are constant and its correlation function only depends on the distance, so that

$$W\sim \mathcal{N}\left(0,{\sigma }^2\rho \left(h\right)\right)$$E4

Given $\mathbf{W}=\left(W_1,\dots ,W_n\right)$, where $W_i=W\left(s_i\right)$, the distribution of $\mathbf{W}$ is normal multivariate $\left(\mathcal{N}\mathcal{M}\right),$ i.e.

$$\mathbf{W}\sim \mathcal{N}\mathcal{M}\left(\mathbf{0},{\sigma }^2\mathbf{R}\right),$$E5

where the $\left(i,j\right)$ element of $\mathbf{R}$ is given by ${\left(\mathbf{R}\right)}_{\mathit{ij}}=\mathit{\text{Corr}}\left(W\left(s_i\right),W\left(s_j\right)\right)=\rho \left(h_{\mathit{ij}}\right)$, $h_{\mathit{ij}}=‖s_i-s_j‖$ is the Euclidean distance between $s_i$ and $s_j$. Note that the covariance of the Gaussian process is given by $\mathit{Cov}\left(\mathbf{W}\right)={\sigma }^2\mathbf{R}$.

In this way, the correlation structure of a stationary Gaussian process can be studied through the $\rho \left(h\right)$ function. Several parametric expressions for this function are shown in the Table 1. In these correlation functions, $\varphi >0$ is a range parameter controlling the spatial decay over distance; $h=‖s-s^{\prime }‖$ is the Euclidean distance between $s$ and $s^{\prime }$ and $h\ge 0$; $\mathrm{\Gamma }\left(\cdot \right)$ denotes the gamma function. $\kappa >0$, in theory of spatial extremes $J_{\kappa }\left(\cdot \right)$ and $K_{\kappa }\left(\cdot \right)$ are the Bessel and modified Bessel function of the third kind of order $\kappa$ [28], while in the spatial survival analysis and generalized linear models $K_{\kappa }\left(\cdot \right)$ is the modified Bessel function of the second kind of order $\kappa$ [29]; $\kappa$ is a shape parameter that determines the analytic smoothness of the underlying process $W$ [1]. In the powered exponential correlation function $0<\kappa \le 2$ and in the Bessel correlation function $\kappa \ge 0$.

Family	Correlation function
Exponential	$\rho \left(h,\varphi \right)=exp\left(-\frac{h}{\varphi }\right)$
Gaussian	$\rho \left(h,\varphi \right)=exp\left(-\frac{h^2}{{\varphi }^2}\right)$
Spherical	$\rho \left(h,\varphi \right)=1-1.5\left(\frac{h}{\varphi }\right)+0.5{\left(\frac{h}{\varphi }\right)}^3$
Circular	$\rho \left(h,\varphi \right)=1-\frac{2}{\pi }\left(a\sqrt{1-a^2}+{sin}^{-1}\sqrt{a}\right)$
Cubic	$\rho \left(h,\varphi \right)=1-\left[7{\left(\frac{h}{\varphi }\right)}^2-\frac{35}{4}{\left(\frac{h}{\varphi }\right)}^3+\frac{7}{2}{\left(\frac{h}{\varphi }\right)}^5-\frac{3}{4}{\left(\frac{h}{\varphi }\right)}^7\right]$
Wave	$\rho \left(h,\varphi \right)=\frac{\varphi }{h}sin\left(\frac{h}{\varphi }\right)$
Matérn	$\rho \left(h,\varphi ,\kappa \right)=\frac{1}{2^{\kappa -1}\mathrm{\Gamma }\left(\kappa \right)}{\left(\frac{h}{\varphi }\right)}^{\kappa }{K}_{\kappa }\left(\frac{h}{\varphi }\right)$
Powered exponential	$\rho \left(h,\varphi ,\kappa \right)=exp\left[-{\left(\frac{h}{\varphi }\right)}^{\kappa }\right]$
Cauchy	$\rho \left(h,\varphi ,\kappa \right)={\left[1+{\left(\frac{h}{\varphi }\right)}^2\right]}^{-\kappa }$
Stable	$\rho \left(h,\varphi ,\kappa \right)=exp\left(-\frac{h}{\varphi }\right)$
Bessel	$\rho \left(h,\varphi ,\kappa \right)={\left(\frac{2\varphi }{h}\right)}^{\kappa }\mathrm{\Gamma }\left(\kappa +1\right)J_{\kappa }\left(\frac{h}{\varphi }\right)$

Table 1.

Models for the spatial correlation structure of a spatial process.

3. Gaussian spatial model

Generally from process $\left\{W\left(s\right):s\in A\subset {\mathbb{R}}^2\right\}$, there is a noisy version, i.e., a set of observation $y\left(s_1\right),\dots ,y\left(s_n\right)$ of the random variables $Y\left(s_1\right),\dots ,Y\left(s_n\right)$, $s_i\in A$. In this way $Y\left(s\right)$ is a measurement process of $W\left(s\right)$, $s\in A$ [1, 26].

The Gaussian geostatistical model, in the absence of independent variables, is given by

$$Y\left(s\right)=\mu +W\left(s\right)+Z\left(s\right),s\in A,$$E6

where $\mu$ is a constant mean effect, $W\left(s\right)$ is a stationary Gaussian process (1) and $Z\left(s\right)$ is the error term in the model with $Z\left(s\right)\sim \mathcal{N}\left(0,{\tau }^2\right)$; ${\tau }^2$ is the nugget effect variance. $Z\left(s\right)$ is known as measurement error, micro-scale variation or a non-identifiable combination of the two [22, 26].

Thus for a realization of a stationary Gaussian spatial process, $\mathbf{Y}\left(s\right)=\left(Y\left(s_1\right),\dots ,Y\left(s_n\right)\right)$, $s_i\in A$ and $i=1,\dots ,n$ with

$$Y\left(s_i\right)=\mu +W\left(s_i\right)+Z\left(s_i\right),$$E7

where

• $W\left(s_i\right)\sim \mathcal{N}\left(0,{\sigma }^2\right)$.

• $Z\left(s_i\right)$ are mutually independent and identically distributed, $Z\left(s_i\right)\sim \mathcal{N}\left(0,{\tau }^2\right),i=1,\dots ,n$.

• $Z\left(s_i\right)$ are independent of the process $W\left(\cdot \right)$ [26].

• Conditional on $W\left(\cdot \right)$, random variables $Y\left(s_i\right)$, $i=1,\dots ,n$, are mutually independent with normal distribution,

$$Y\left(s_i\right)\mid W\left(\cdot \right)\sim \mathcal{N}\left(\mu +W\left(s_i\right),{\tau }^2\right).$$E8

The joint distribution of $\mathbf{Y}\left(s\right)$ is normal multivariate given by

$$\mathbf{Y}\left(s\right)\sim \mathcal{N}\mathcal{M}\left(\mu \mathbf{1},{\sigma }^2\mathbf{R}\left(\varphi \right)+{\tau }^2\mathbf{I}\right),$$E9

where

• $\mu$ is the mean of the Gaussian process $W\left(\cdot \right)$ and $\mathbf{1}$ is a vector of dimension $n\times 1$_.

• ${\sigma }^2$ is the variance of the process $W\left(\cdot \right)$_.

• $\mathbf{R}\left(\varphi \right)$ is a matrix of correlations of dimension $n\times n$, whose elements given by

$${\left(\mathbf{R}\left(\varphi \right)\right)}_{\mathit{ij}}={\rho }_Y\left(h_{\mathit{ij}},\varphi \right),$$E10

where $h_{\mathit{ij}}=‖s_i-s_j‖$ is the euclidean distance that exists between $s_i$ and $s_j$ that are in $A$, and $\varphi$ is a spatial scale parameter.

• ${\tau }^2$ is the variance of $Z$_and$\mathbf{I}$ is the identity matrix of dimension $n\times n$_.

• Note that the covariance of the $\mathbf{Y}\left(s\right)$ is given by $\mathit{Cov}\left(\mathbf{Y}\left(s\right)\right)={\sigma }^2\mathbf{R}\left(\varphi \right)+{\tau }^2\mathbf{I}$_.

When $\mathbf{Y}\left(s\right)$ can be explained by a set of covariates that also depend on the location, $\mathbf{X}\left(s\right)=\left(X_1\left(s\right),.\dots ,X_p\left(s\right)\right)$, then the model is given by

$$\mathbf{Y}\left(s\right)=\mathbf{X}\left(\mathbf{s}\right)\beta +\mathbf{W}\left(s\right)+\mathbf{Z}\left(s\right),s\in A,$$E11

with

$$\mathbf{Y}\left(s\right)\sim \mathcal{N}\mathcal{M}\left(\mathbf{X}\left(\mathbf{s}\right)\beta ,{\sigma }^2\mathbf{R}\left(\varphi \right)+{\tau }^2\mathbf{I}\right),$$E12

where $\beta =\left({\beta }_0,\dots ,{\beta }_p\right)$ is a vector of unknown regression parameters; in this case also $\mathit{Cov}\left(\mathbf{Y}\left(s\right)\right)={\sigma }^2\mathbf{R}\left(\varphi \right)+{\tau }^2\mathbf{I}$. The unknown parameters in this model are $\beta$, ${\sigma }^2$, ${\tau }^2$ and $\varphi$. The parameters of the Models (4) y (5) can be estimated under the classical approach (maximum likelihood or maximum restricted likelihood) and under the Bayesian statistical approach [1, 30]. Among the most important points in geostatistics is the modeling of the covariance structure of the spatial process and the identification of the interpolation method that will be used to perform the prediction of the process in the non sampled points in $A$. Regarding the last point, [31] made a compilation of the most used criteria for assessing the performance of the spatial interpolation method.

The geoR package contains the likfit function that allows to estimate, under Maximum likelihood (ML) or restricted maximum likelihood (REML), the parameters of a Gaussian process [32]. The function likfit estimates the coefficients of the models (4) y (5).

The function krige.cov of the same package helps to perform the spatial prediction of a Gaussian process using simple kriging (SK), ordinary kriging (OK), external trend kriging (KTE) and universal kriging (UK) [33]. The package glmmfields allows to fit Gaussian models [34] under the Bayesian approach.

On the other hand, with the function glmmfields of the package glmmfields, the coefficients of the models (4) and (5) can be estimated under the Bayesian approach. The function glmmfields reports the posterior median of the parameters with their respective $95\%$ credible intervals; this function, also reports the values of the Gelman and Rubin statistic [35], where values less than $1.20$ would indicate convergence of the chain.

4. Generalized linear spatial models

Generalized linear models (GLM) [36, 37] are very useful when the response variable does not follow a normal distribution. The assumptions of GLMs are

1. 1. $Y_i$, $i=1,\dots ,n$ are mutually independent with expectations ${\mu }_i$.

2. The ${\mu }_i$ are specified by $g\left({\mu }_i\right)={\eta }_i$, where $g\left(\cdot \right)$ is a known link function.

3. The linear predictor is given by ${\eta }_i=x_i^{\prime }\beta$, where $x_i$ is a vector of explanatory variables associated with the response $Y_i$, and $\beta$ is a vector of unknown parameters.

The $Y_i$ follow a common distributional family, indexed by their expectations ${\mu }_i$, and possibly by additional parameters common to all $n$ responses.

An important extension of this basic class of models is the generalized linear mixed model (GLMM) [38], in which $Y_1,\dots ,Y_n$ are mutually independent conditional on the realized values of a set latent random variables (random effects) $U_1,\dots ,U_n$ and the conditional expectations are given by $g\left({\mu }_i\right)=U_i+x_i^{\text{'}}\beta$. A generalized linear spatial model is a GLMM in which the $U_1,\dots ,U_n$ are derived from spatial process. Diggle and Ribeiro in $2007$ [1], refers to these models as generalized linear geostatistical model (GLGM). In accordance with Diggle et al. [39], the assumptions of the generalized linear spatial models are as follows

1. $W$ is a stationary Gaussian process, $W\sim \mathcal{N}\left(0,{\sigma }^2\rho \left(h\right)\right)$, (Eq. (1)).

2. Conditionally an $W$, the random variables $Y_i$, $i=1,\dots ,n$ are mutually independent, with distributions $f_i\left(y|W\left(s_i\right)\right)=f\left(y;M_i\right)$, specified by the values of the conditional expectations $M_i=E\left[Y_i|W\left(s_i\right)\right]$.

3. $g\left(M_i\right)={\mathbf{x}}_i^{\prime }\beta +W\left(s_i\right)$ for some known link function $g$ and explanatory variable ${\mathbf{x}}_i=\mathbf{x}\left(s_i\right)$.

Then $M_i=g^{-1}\left({\mathbf{x}}_i^{\prime }\beta +W\left(s_i\right)\right)$, where the linear predictor would be given by ${\eta }_i={\mathbf{x}}_i^{\prime }\beta +W\left(s_i\right)$.

Taking Diggle and Tawn as a precedent (1998) [39]; Jing and De Oliveira in $2015$ [40] state the GLSM as follows

$$Y_i\mid W_i\sim p\left(\cdot |{\mu }_i\right).$$E13

where

$$\mathbf{W}\sim \mathcal{N}\mathcal{M}\left(\mathbf{X}\beta ,{\sigma }^2\mathbf{R}\right)$$E14

$\mathbf{R}$ is of the same form as the Gaussian process (2)

• $Y\left(s_i\right):i=1,\dots ,n$ are conditionally independent given $\mathbf{W}$ with pdfs or pmfs $p\left(\cdot |{\mu }_i\right)$.

• $E\left(Y_i|W_i\right)={\mu }_i$ and $g\left(\cdot \right)$ is a known one-to-one link function.

• $\mathbf{X}=\left(\mathbf{1},{\mathbf{x}}_1,\dots ,{\mathbf{x}}_p\right)$ is a known $n\times \left(p+1\right)$ design matrix assumed of full-rank, with $\mathbf{1}$ a vector of $n\times 1$ of ones and ${\mathbf{x}}_j={\left(x_j\left(s_1\right),.\dots ,x_j\left(s_n\right)\right)}^{\text{'}}$, where $x_j\left(s_i\right)$ is the value of the $j$-th covariate of the $i$-th sampling location, and $\beta =\left({\beta }_0,{\beta }_1,\dots ,{\beta }_p\right)$ is the vector of unknown regression parameters.

Since $g$ is the link function then $g\left({\mu }_i\right)={\eta }_i$ and ${\mu }_i=g^{-1}\left({\eta }_i\right)$, $i=1,\dots ,n$, where the linear predictor is given by ${\eta }_i=W_i$, then ${\mu }_i=g^{-1}\left(W_i\right)$. The unknown parameters in GLSM are $\beta$, ${\sigma }^2$ and $\varphi$.

The two most widely used GLSM for spatial count data are the Poisson and Binomial spatial models [39, 41].

The geoCount [40] package implements the GLSM; the function runMCMC is used to generate posterior samples of the Gaussian process and the GLSM parameters, with which the parameter estimates and their credibility intervals can be obtained.

In the package geoRglm [42, 43], the functions glsm.krige, pois.krige and binom.krige implement the GLSMs, in this case, parameter estimation is performed under the frequentist approach. While the functions krige.bayes, pois.krige.bayes and binom.krige.bayes, which also implement the GLSMs, estimate the parameters under the Bayesian approach. These functions report estimates of $\beta$, ${\sigma }^2$ and $\varphi$.

The glmmfields package implements the Gamma, Poisson, Negative Binomial, Binomial and Lognormal models using the function glmmfields [34], parameter estimation is performed under the Bayesian approach. The function glmmfields reports the parameter estimates using the posterior median with their respective $95\%$ percentile credible intervals; it also reports the Gelman and Rubin diagnostic values.

5. Spatial survival models

Generally, survival analysis models are specified through their hazard function, $h\left(t\right)$, whose intuitive interpretation is that $h\left(t\right)\mathit{\delta t}$ is the conditional probability that a patient will die in the interval $\left(t,t+\mathit{\delta t}\right)$, given tat they have survived until time $t$. The most widely used approach to modeling $h\left(t\right)$, at least in medical applications, is to use a semi-parametric formulation [44]. In this approach, the hazard for the $i$-th patient is modeled as

$$h\left(t_i\right)=h_0\left(t_i\right)exp\left(x_i^{\prime }\beta \right),$$E15

where $x_i$ is a vector of explanatory variables for patient $i$ and $h_0\left(t\right)$ is an unspecified baseline hazard function. This is known as a proportional hazards (PH) model, because for any two patients $i$ and $j$, $h\left(t_i\right)/h\left(t_j\right)$ does not change over time [1].

Another key idea in survival analysis is frailty, this corresponds to the random effects term used; time-to-event data will be group into strata, such as clinical sites, geographic regions, etc. This gives rise to mixed models, which include a random effect (the frailty) that correspond to a stratum’s overall health status [30]. To illustrate, let $t_{\mathit{ij}}$ be the time to death or censoring for subject $j$ in stratum $i$, $j=1,\dots ,n_i$, $i=1,\dots ,m$. Let ${\mathbf{x}}_{\mathit{ij}}$ be a vector of individual specific covariates, then

$$h\left(t_{\mathit{ij}},x_{\mathit{ij}}\right)=h_0\left(t_{\mathit{ij}}\right)exp\left({\mathbf{x}}_{\mathit{ij}}^{\prime }\beta +W_i\right),$$E16

where $W_i$ es the stratum-specific frailty term, designed to capture differences among strata; strata are typically denoted by $s_i$, $i=1,\dots ,m$, so $s_i$ denotes the location of the $i$-th patient and $W_i=W\left(s_i\right)$. It can be assumed that the $W_i$ are independent identical distribution (iid), i.e.

$$W_i\sim \mathcal{N}\left(0,{\sigma }^2\right).$$E17

But it can also be assumed that $W_i$ arises from a Gaussian process, i.e. if $\mathbf{W}=\left(W_1,\dots ,W_m\right)$, then

$$\mathbf{W}\sim \mathcal{N}\mathcal{M}\left(\mathbf{0},{\sigma }^2\mathbf{R}\left(\varphi \right)\right).$$E18

This way, suppose subjects are observed at $m$ distinct spatial locations $s_1,\dots ,s_m\in A$. Let $t_{\mathit{ij}}$ be a random event time associated with the $j$-th subject in $s_i$, assume the survival time $t_{\mathit{ij}}$ lies in the interval $\left(a_{\mathit{ij}},b_{\mathit{ij}}\right)$, $i=1,\dots ,m$, $j=1,\dots ,n_i$; and ${\mathbf{x}}_{\mathit{ij}}$ be a related p-dimensional vector of covariates, then are defined proportional hazard (PH) frailty models, accelerated failure time (AFT) frailty models and proportional odds (PO) frailty models.

PH frailty models are the extensions of the population hazards model which is best known as the Cox model [44] a widely pursued model in survival analysis. PH frailty models extends the Cox model such that the hazard of an individual depends in addition on an unobserved random variable $W$, then introducing an additive frailty term $W_i$ for each individual in the exponent of the hazard function as follows

$$h\left(t_{\mathit{ij}},{\mathbf{x}}_{\mathit{ij}}\right)=h_0\left(t_{\mathit{ij}}\right)e^{{\mathbf{x}}_{\mathit{ij}}^{\prime }\beta +W_i}.$$E19

The corresponding survival function and the density are given by

$$\begin{array}{c}S\left(t_{\mathit{ij}},{\mathbf{x}}_{\mathit{ij}}\right)=S_0{\left(t_{\mathit{ij}}\right)}^{e^{{\mathbf{x}}_{\mathit{ij}}^{\prime }\beta +W_i}},\\ f\left(t_{\mathit{ij}},{\mathbf{x}}_{\mathit{ij}}\right)=e^{{\mathbf{x}}_{\mathit{ij}}^{\prime }\beta +W_i}{S}_0{\left(t_{\mathit{ij}}\right)}^{e^{{\mathbf{x}}_{\mathit{ij}}^{\prime }\beta +W_i}-1}{f}_0\left(t_{\mathit{ij}}\right),\end{array}$$E20

where $S_0\left(\cdot \right)$, $f_0\left(\cdot \right)$ and $h_0\left(\cdot \right)$ are the baseline survival function, baseline density and baseline hazard function assumed to be unique for all individual in the study population respectively.

Accelerated failure time frailty model extends the AFT model such that the hazard of an individual depends in addition on an unobserved random variable $W$ [45, 46, 47]. Introducing an additive frailty term $W_i$ for each individual in the exponent of the hazard function it becomes:

$$h\left(t_{\mathit{ij}},{\mathbf{x}}_{\mathit{ij}}\right)=h_0\left(e^{{\mathbf{x}}_{\mathit{ij}}^{\prime }\beta +W_i}{t}_{\mathit{ij}}\right)e^{{\mathbf{x}}_{\mathit{ij}}^{\prime }\beta +W_i}.$$E21

The survival function and density are given by

$$\begin{array}{c}S\left(t_{\mathit{ij}},{\mathbf{x}}_{\mathit{ij}}\right)=S_0\left(e^{{\mathbf{x}}_{\mathit{ij}}^{\prime }\beta +W_i}{t}_{\mathit{ij}}\right),\\ f\left(t_{\mathit{ij}},{\mathbf{x}}_{\mathit{ij}}\right)=e^{{\mathbf{x}}_{\mathit{ij}}^{\prime }\beta +W_i}{f}_0\left(e^{{\mathbf{x}}_{\mathit{ij}}^{\prime }\beta +W_i}{t}_{\mathit{ij}}\right),\end{array}$$E22

where $S_0\left(\cdot \right)$, $f_0\left(\cdot \right)$ and $h_0\left(\cdot \right)$ are the baseline survival function, baseline density and baseline hazard function assumed to be unique for all individual in the study population respectively.

Finally, proportional odds frailty model is given by

$$h\left(t_{\mathit{ij}},{\mathbf{x}}_{\mathit{ij}}\right)=h_0\frac{1}{1+\left[e^{-{\mathbf{x}}_{\mathit{ij}}^{\prime }\beta -W_i}-1\right]S_0\left(t_{\mathit{ij}}\right)}.$$E23

The survival function and density are given by

$$\begin{array}{c}S\left(t_{\mathit{ij}},{\mathbf{x}}_{\mathit{ij}}\right)=\frac{S_0\left(t_{\mathit{ij}}\right)e^{-{\mathbf{x}}_{\mathit{ij}}^{\prime }\beta -W_i}}{1+\left(e^{-{\mathbf{x}}_{\mathit{ij}}^{\prime }\beta -W_i}-1\right)S_0\left(t_{\mathit{ij}}\right)},\\ f\left(t_{\mathit{ij}},{\mathbf{x}}_{\mathit{ij}}\right)=\frac{f_0\left(t_{\mathit{ij}}\right)e^{-{\mathbf{x}}_{\mathit{ij}}^{\prime }\beta -W_i}}{{\left[1+\left(e^{-{\mathbf{x}}_{\mathit{ij}}^{\prime }\beta -W_i}-1\right)S_0\left(t_{\mathit{ij}}\right)\right]}^2},\end{array}$$E24

where $S_0\left(\cdot \right)$, $f_0\left(\cdot \right)$ and $h_0\left(\cdot \right)$ are the baseline survival function, baseline density and baseline hazard function assumed to be unique for all individual in the study population respectively.

In the frailty models, it is possible to deal with left, right and interval censoring of the data. Among the packages that exist in the R statistical software to perform spatial survival analysis is the spBayesSurv package [48]; the function survregbayes estimates the parameters of the PH, AFT and PO spatial models under the classical and Bayesian approach; also reports the posterior mean and median of the regression coefficients and of the parameters of the covariance function of the Gaussian process, ${\sigma }^2$ and $\varphi$, with their $95\%$ credible intervals. The spBayesSurv package uses the powered exponential function (Table 0) to model the spatial correlation of the data.

Also in R, there is the spatsurv package [49], which implements the function survspat that fits parametric PH spatial survival models. This function reports the estimates and posterior median of the parameters $\beta$, ${\sigma }^2$ and $\varphi$ with the respective credibility intervals.

6. Spatial generalized extreme value model

According to Coles (2001) [50], given $Y_1,\dots ,Y_n$ a sequence of independent random variables with a common distribution function $F$ with $M_n=max\left\{Y_1,\dots ,Y_n\right\}$, if there a sequences of constants $a_n>0$ and $b_n$ such that

$$P\left(\left(M_n-b_n\right)/a_n\le z\right)\to G\left(z\right),$$E25

when $n\to \infty$, for a non-degenerative distribution function $G$, then $G$ is a member of the generalized extreme value (GEV) distribution family

$$G\left(y;\eta ,\tau ,\xi \right)=exp\left\{-{\left[1+\xi \left(\frac{y-\eta }{\tau }\right)\right]}^{\frac{-1}{\xi }}\right\},$$E26

defined on $\left\{z:1+\xi \left(z-\eta \right)/\tau >0\right\}$, where $-\infty <\eta <\infty$, $\tau >0$ and $-\infty <\xi <\infty$.

Davison et al. in $2012$ [51], describe spatial GVE as follows. For each $s$ in ${\mathbb{R}}^2$, suppose that $Y\left(s\right)$ is GEV distributed whose parameters $\mu \left(s\right)$, $\sigma \left(s\right)$ and $\xi \left(s\right)$ vary smoothly for $s$ in ${\mathbb{R}}^2$ according to a stochastic process $W\left(s\right)$. We assume that the processes for each GEV parameters are mutually independent Gaussian processes [52]. Then

$$\begin{array}{c}\eta \left(s\right)=f_{\eta }\left(s;{\beta }_{\eta }\right)+W_{\eta }\left(s;{\sigma }_{\eta },{\varphi }_{\eta },{\kappa }_{\eta }\right),\\ \tau \left(s\right)=f_{\tau }\left(s;{\beta }_{\tau }\right)+W_{\tau }\left(s;{\sigma }_{\tau },{\varphi }_{\tau },{\kappa }_{\tau }\right),\\ \xi \left(s\right)=f_{\xi }\left(s;{\beta }_{\xi }\right)+W_{\xi }\left(s;{\sigma }_{\xi },{\varphi }_{\xi };{\kappa }_{\xi }\right),\end{array}$$E27

where $f_{\eta }$, $f_{\tau }$ and $f_{\mathit{xi}}$ are deterministic functions depending on a regression parameters ${\beta }_{\eta }$, ${\beta }_{\tau }$ and ${\beta }_{\xi }$ respectively. While $W_{\eta }$, $W_{\tau }$ and $W_{\xi }$ are a zero mean stationary Gaussian process with correlation function $\rho \left(h,{\varphi }_{\eta },{\kappa }_{\eta }\right)$, $\rho \left(h,{\varphi }_{\tau },{\kappa }_{\tau }\right)$ and $\rho \left(h,{\varphi }_{\xi },{\kappa }_{\xi }\right)$ respectively, i.e.

$$\begin{array}{c}{W}_{\eta }\sim \mathcal{N}(0,{\sigma }_{\eta }^2\rho \left(h,{\varphi }_{\eta },{\kappa }_{\eta }\right),\\ W_{\tau }\sim \mathcal{N}(0,{\sigma }_{\tau }^2\rho \left(h,{\varphi }_{\tau },{\kappa }_{\tau }\right),\\ W_{\xi }\sim \mathcal{N}(0,{\sigma }_{\xi }^2\rho \left(h,{\varphi }_{\xi },{\kappa }_{\xi }\right).\end{array}$$E28

Then conditional on the values of the tree Gaussian process at the sites $\left(s_1,\dots ,s_k\right)$, the maxima are assumed to follow GEV distributions

$$Y_{s_i}\mid \eta \left(s_j\right),\tau \left(s_j\right),\xi \left(s_j\right)\sim \mathit{GEV}\left(\eta \left(s_j\right),\tau \left(s_j\right),\xi \left(s_j\right)\right)$$E29

z independently for each location $s_1,\dots ,s_k$, $j=1,\dots ,k$ and $i=1,\dots ,n$.

Davison et al. in $2012$ [51], proposed the construction of Bayesian hierarchical models for spatial extremes.

The SpatialExtremes package [53] allows modeling spatial extremes, through max-stable processes with the function fitmaxstab, which reports the values of the parameter estimates with their respective standard errors.

To implement hierarchical Bayesian models, the function latent is used, this reports the posterior median of the scale, shape and location parameters with their respective credible intervals.

Another package in the literature to model spatial extremes is glmmfields [34], with the function glmmfields, parameter estimation is performed under the Bayesian approach. The function glmmfields also allows modeling spatial extreme events incorporating temporally, that is, time, these models are known as spatio-temporal models.

7. Conclusions

The main characteristic of spatial data is that observations close in space tend to be correlated, and in spatial modeling this correlation is used to understand the behavior of the phenomenon under study in a region of interest.

Omitting the spatial dependence of the data can generate a bias of the information and, consequently, lead to an incorrect inference. Therefore, adequately describing the spatial pattern of an event can provide sufficient elements to elaborate possible hypotheses of its cause. As we have seen, the spatial variability of georeferenced data can be studied with the spatial models developed in geostatistics. The usefulness of these models has been demonstrated in several applications related to the identification of social structures, disease patterns, occupational patterns, as well as in the identification of populations (or subgroups) that are at greater or lesser risk of an event. As we have seen, in statistics, all correctly processed information helps in correct decision making. In this sense, this paper aims to introduce the reader to the use of spatial models in geostatistics.

If the response or variable of interest is the cases (counts) of sick people in a given region, or the new cases of a disease in a given period of time (incidence), then Poisson GLSMs can be useful to know the spread of the disease in the population of interest, predict new cases, and identify the variables that influence the occurrence of the disease. On the other hand, when the response variable is a binary or ratio variable, such as mortality rates or infection rates, then binomial GLSMs can be helpful. These models have been used to study the prevalence of dengue and to identify the variables associated with the event.

Conflict of interest

The authors declare no conflict of interest.

Abbreviations

AFT	Accelerated failure time
EVT	Extreme value theory.
GEV	generalized extreme value.
GLSM	Generalized Linear Spatial Models.
GLM	Generalized linear models.
GLMM	generalized linear mixed model
GLGM	generalized linear geostatistical model
iid	Independent identical distribution
NM	Normal multivariate.
PH	Proportional hazards.
PO	Proportional odds.