**3. Gaussian spatial model**

Generally from process *W s*ð Þ : *<sup>s</sup>* <sup>∈</sup> *<sup>A</sup>* <sup>⊂</sup> <sup>2</sup> � �, there is a noisy version, i.e., a set of observation *y s*ð Þ<sup>1</sup> , … , *y s*ð Þ*<sup>n</sup>* of the random variables *Y s*ð Þ<sup>1</sup> , … , *Y s*ð Þ*<sup>n</sup>* , *si* ∈ *A*. In this way *Y s*ð Þ is a measurement process of *W s*ð Þ, *s*∈ *A* [1, 26].

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

$$Y(\mathfrak{s}) = \mu + \mathcal{W}(\mathfrak{s}) + Z(\mathfrak{s}), \ \mathfrak{s} \in A,\tag{6}$$

where *μ* is a constant mean effect, *W s*ð Þ is a stationary Gaussian process (1) and *Z s*ð Þ is the error term in the model with *Z s*ðÞ� <sup>N</sup> 0, *<sup>τ</sup>*<sup>2</sup> ð Þ; *<sup>τ</sup>*<sup>2</sup> is the nugget effect variance. *Z s*ð Þ 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, **Y**ðÞ¼ *s* ð Þ *Y s*ð Þ<sup>1</sup> , … , *Y s*ð Þ*<sup>n</sup>* , *si* ∈ *A* and *i* ¼ 1, … , *n* with

$$Y(\mathfrak{s}\_i) = \mu + W(\mathfrak{s}\_i) + Z(\mathfrak{s}\_i),\tag{7}$$

where


$$Y(\mathfrak{s}\_i)|\mathcal{W}(\cdot) \sim \mathcal{N}\left(\boldsymbol{\mu} + \mathcal{W}(\mathfrak{s}\_i), \boldsymbol{\tau}^2\right). \tag{8}$$

The joint distribution of **Y**ð Þ*s* is normal multivariate given by

$$\mathbf{Y}(\mathfrak{s}) \sim \mathcal{N} \mathcal{M}(\mathfrak{\mu}\mathbf{1}, \sigma^2 \mathbf{R}(\phi) + \tau^2 \mathbf{I}),\tag{9}$$

where


$$(\mathbf{R}(\phi))\_{\vec{\eta}} = \rho\_Y(h\_{\vec{\eta}}, \phi), \tag{10}$$

where *hij* ¼ *si* � *sj* is the euclidean distance that exists between *si* and *sj* that are in *A*, and *ϕ* is a spatial scale parameter.


When **Y**ð Þ*s* can be explained by a set of covariates that also depend on the location, **<sup>X</sup>**ðÞ¼ *<sup>s</sup> <sup>X</sup>*1ð Þ*<sup>s</sup>* , *:* … ,*Xp*ð Þ*<sup>s</sup>* , then the model is given by

$$\mathbf{Y}(\mathbf{s}) = \mathbf{X}(\mathbf{s})\boldsymbol{\beta} + \mathbf{W}(\mathbf{s}) + \mathbf{Z}(\mathbf{s}), \ \mathbf{s} \in A,\tag{11}$$

with

$$\mathbf{Y}(\mathfrak{s}) \sim \mathcal{N}\mathcal{M}(\mathbf{X}(\mathfrak{s})\boldsymbol{\beta}, \ \sigma^2 \mathbf{R}(\boldsymbol{\phi}) + \tau^2 \mathbf{I}),\tag{12}$$

where *β* ¼ *β*0, … , *β<sup>p</sup>* is a vector of unknown regression parameters; in this case also *Cov*ð Þ¼ **<sup>Y</sup>**ð Þ*<sup>s</sup> <sup>σ</sup>*<sup>2</sup>**R**ð Þþ *<sup>ϕ</sup> <sup>τ</sup>*<sup>2</sup>**I**. The unknown parameters in this model are *<sup>β</sup>*, *<sup>σ</sup>*2, *<sup>τ</sup>*<sup>2</sup> and *ϕ*. 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.
