Incorporating Biotic Information in Species Distribution Models: A Coregionalized Approach

Abstract

In this work, we discuss the use of a methodological approach for modelling spatial relationships among species by means of a Bayesian spatial coregionalized model. Inference and prediction is performed using the integrated nested Laplace approximation methodology to reduce the computational burden. We illustrate the performance of the coregionalized model in species interaction scenarios using both simulated and real data. The simulation demonstrates the better predictive performance of the coregionalized model with respect to the univariate models. The case study focus on the spatial distribution of a prey species, the European anchovy (Engraulis encrasicolus), and one of its predator species, the European hake (Merluccius merluccius), in the Mediterranean sea. The results indicate that European hake and anchovy are positively associated, resulting in improved model predictions using the coregionalized model.

1. Introduction

Gaining knowledge about where the species are present has become nearly mandatory in many areas of research such as Ecology and Epidemiology. The most popular way of modelling the distribution of a species are the so-called Species Distribution Models (SDMs) (see, e.g., [1,2,3], for detailed revisions on SDMs). These models basically link species observations with environmental data with the final purpose of predicting where the species are in unsampled locations or time periods.

The inherent complexity of nature and the wide variety of sampling approaches encountered in ecology have prompted researchers to design SDMs that carefully consider different sources of bias and error. Martínez-Minaya et al. [4] have reviewed some of them (the presence of a temporal effect, preferential sampling, spatial misalignment, non-stationarity, imperfect detection, and the excess of zeros) while presenting some guidance about how to deal with the inherent complexity.

In this work, we deal with another important source of complexity that appears when analyzing species: the fact that species interact among them. So, when predicting where a species is likely to be present, we should try to include as many biotic relationships as possible in the modelling [5]. Species could compete, leading them to only partially occupy their potential habitat [6,7], but they could also be related as predator-prey [8,9,10]. Other relationships such as parasitism [11,12] and mutualism [13,14] could also explain the possible presence/absence or the abundance of a species. Despite its enormous relevance, these factors are completely ignored in most classical SDM-based studies.

There are two ways of incorporating these relationships in a model, the first one being through the covariates [15,16]. The second one consists of considering a multivariate response. In the case of SDMs, this leads us to a multivariate geostatistical approach (see for example, [17,18]). Although there have been various ways for analysing multivariate geostatistics, all of them have faced the problem of predicting in non-sampled locations. A good option is to use a hierarchical Bayesian conditional coregionalized linear model for multivariate data based on intrinsic correlation [19,20], that allows us to simultaneously predict the values of related species (in non-sampled locations) based on environmental covariates and common residual spatial patterns. The Bayesian paradigm has proved to be a good alternative to tackle spatial hierarchical models. The reason comes from the ability to treat observed data and model parameters as random variables, which in general provides a proper assessment of uncertainty in the spatial context [20,21].

Despite the enormous potential of hierarchical Bayesian conditional coregionalized linear models, still there are not too many applications of it [22,23,24]. The contribution of this work is to fill this gap presenting both a simulation study of its behavior and a practical example about fisheries.

As usual, the more complex a model, the more difficult is to perform both inferences in the parameters governing it and prediction. In this case, the resulting posterior distributions of the Bayesian conditional coregionalized linear model here used do not have closed expressions. Moreover, the resulting posterior predictive distribution for the unsampled locations do not have either closed expressions, and so, in both cases, numerical approaches are needed to approximate them. We propose the use of the Integrated Nested Laplace Approximation (INLA) methodology [25] (see http://www.r-inla.org accessed on 11 January 2021 for more details). Nowadays, INLA has become a very popular alternative to Markov Chain Monte Carlo (MCMC) methods [26] due to its speed of calculation and the ease with which model comparison and prediction can be handled. Indeed, as stated in [24,27], performing prediction via MCMC and the compositional method can be effective although computationally expensive in those cases with a large number of new locations to make a prediction. Not even the use of parallelization can achieve the computational burden reduction provided by INLA.

The remainder of this article is organized as follows. In Section 2, we describe the hierarchical Bayesian Coregionalisation model and how to perform inference and prediction on it. Section 3 provides a simulation study to show how the use of multivariate conditional models outperforms the univariate case. In Section 4, we apply this methodology in a real setting modeling a prey–predator interaction between the European anchovy (Engraulis encrasicolus) and the European hake (Merluccius merluccius), in the Mediterranean Sea. The final section concludes the paper and presents lines for future research.

2. Coregionalized Models for Multivariate SDMs

For the sake of simplicity and in line with [24], in what follows, we first present a method for modeling two species (three or more would be similar) based on the conditional regionalisation approach [22]. In the second place, we describe how to perform inference and prediction on it.

2.1. The Hierarchical Bayesian Coregionalisation Model

Let $\mathit{Z}=(Z_1,Z_2)$ be a random vector of the abundance of two species and $\mathit{s}=({\mathit{s}}_1,\dots ,{\mathit{s}}_n)$ a subset of n locations in a region. Then, $\mathit{Z}\left(\mathit{s}\right)={(\mathit{Z}{\left({\mathit{s}}_1\right)}^T,\mathit{Z}{\left({\mathit{s}}_2\right)}^T,\dots ,\mathit{Z}{\left({\mathit{s}}_n\right)}^T)}^T$ is a multivariate vector of $2n$ length that represents the abundances at all locations. Let also $\mathit{X}\left({\mathit{s}}_i\right)$ represent the linear predictor (which can be formed with environmental or any other covariates of interest) for location i. The multivariate spatial regression model is then:

$$\mathit{Z}\left({\mathit{s}}_i\right)=\left[{\mathit{I}}_2\otimes {\mathit{X}}^T\left({\mathit{s}}_i\right)\right]\mathit{\beta }+\mathit{W}\left({\mathit{s}}_i\right)+\mathit{\epsilon }\left({\mathit{s}}_i\right),$$

(1)

where $\mathit{\beta }=({\mathit{\beta }}_1,{\mathit{\beta }}_2)$ is the vector of parameters (each ${\mathit{\beta }}_j$ having its own dimension depending of the linear predictor dimension); $\mathit{W}\left(\mathit{s}\right)$ is the Gaussian distributed spatial effect with cross-covariance function ${\mathit{C}}_{\mathit{W}}$, $\mathit{W}\left(\mathit{s}\right)\sim {\mathcal{N}}_{2n}(\mathbf{0},{\mathit{C}}_{\mathit{W}}(·,·))$; and, finally $\mathit{\epsilon }\left(\mathit{s}\right)$ corresponds to the uncorrelated error vector, also Gaussian distributed, ${\epsilon }_j\left(s_i\right)\sim \mathcal{N}(0,{\psi }_j^2)$, that explains the small-scale variability.

The above-mentioned cross-covariance function is a 2 × 2 matrix-valued function, which it is defined for any pairs of locations as

$${\mathit{C}}_{\mathit{W}}({\mathit{s}}_i,{\mathit{s}}_k)=\left[Cov\left({\mathit{W}}_j\left({\mathit{s}}_i\right),{\mathit{W}}_l\left({\mathit{s}}_k\right)\right)\right]\mathrm{with}j,l=1,2.$$

(2)

There are several cross covariance functions for multivariate geostatistics (see [28], for a detailed review of functions). In our case, and in line with [29] that propose a conditional approach to multivariate spatial covariance models, we define it by means of a correlation function $\rho$ and a positive definite matrix $\mathit{D}$,

$${\mathit{C}}_{\mathit{W}}({\mathit{s}}_i,{\mathit{s}}_k)=\rho ({\mathit{s}}_i,{\mathit{s}}_k)\mathit{D}.$$

(3)

The resulting covariance matrix of $\mathit{W}$ turns out to be ${\mathbf{\Sigma }}_{\mathit{W}}=\mathit{H}\otimes \mathit{D}$, being ${\mathit{H}}_{ik}=\rho ({\mathit{s}}_i,{\mathit{s}}_k)$.

In order to model bivariate species abundance distribution, we focus on the regionalization idea, in particular, the one proposed by [22], which takes advantage of the Bayesian approach. The basic idea is to express the spatial model in (1) hierarchically [30], that is:

$$\begin{array}{cc}\hfill \left(I\right)& \mathit{Z}\left(\mathit{s}\right)|\mathit{\beta },\mathit{W},\mathit{\psi }\sim {\mathcal{N}}_{2n}\left({\mathit{X}}^T\left(\mathit{s}\right)\mathit{\beta }+\mathit{W}\left(\mathit{s}\right),\mathit{\psi }\otimes {\mathit{I}}_n\right)\\ \hfill \left(II\right)& \mathit{W}\left(\mathit{s}\right)|\mathit{\theta }\sim {\mathcal{N}}_{2n}\left(\mathbf{0},{\mathbf{\Sigma }}_{\mathit{Z}}\left(\mathit{s}\right)\right)\\ \hfill \left(III\right)& p(\mathit{\beta },\mathit{\psi },\mathit{\theta }),\end{array}$$

(4)

where $\mathit{\psi }=\mathrm{diag}({\mathit{\psi }}_1^2,{\mathit{\psi }}_2^2)$ and $\mathit{\theta }$ stands for the parametric vectors of the correlation functions of $\mathit{W}\left(\mathit{s}\right)$. As it can be noticed, the last level requires the elicitation of the priors (of the parameters), and the hyperpriors (of the hyperparameters) of the model.

In line with [20], the model in (4) can be rewritten keeping interpretability but gaining flexibility and, more importantly, computational tractability. More precisely, the joint distribution of the abundances of both species can be expressed as the product of the following conditional distributions

$$p\left[{\mathit{Z}}_1,{\mathit{Z}}_2\right]=p\left[{\mathit{Z}}_2|{\mathit{Z}}_1\right]p\left[{\mathit{Z}}_1\right],$$

(5)

where from the sake of simplicity, $\mathit{Z}\left(\mathit{s}\right)$, $\mathit{X}\left(\mathit{s}\right)$ and $\mathit{W}\left(\mathit{s}\right)$ are, from now on, denoted by $\mathit{Z}$, $\mathit{X}$ and $\mathit{W}$, respectively. Note that this assumption implies that the cross-covariance function in (3) is separable [28]. Based on this, expression (4) becomes

$$\begin{array}{cc}\hfill \left(I\right)& \left\{\begin{array}{ccc}\hfill {\mathit{Z}}_1& \sim & {\mathcal{N}}_n\left({\mathit{X}}^T{\mathit{\beta }}_1+{\mathit{W}}_1,{\psi }_1^2{\mathit{I}}_n\right)\hfill \\ \hfill {\mathit{Z}}_2|{\mathit{Z}}_1& \sim & {\mathcal{N}}_n\left({\mathit{X}}^T{\mathit{\beta }}_2+\alpha {\mathit{W}}_1+{\mathit{W}}_2,{\psi }_2^2{\mathit{I}}_n\right)\hfill \end{array}\right.\\ \hfill \left(II\right)& \left\{\begin{array}{ccc}\hfill {\mathit{W}}_1& \sim & {\mathcal{N}}_n\left(\mathbf{0},{\sigma }_1^2{\mathit{H}}_1\left({\mathit{\theta }}_1\right)\right)\hfill \\ \hfill {\mathit{W}}_2& \sim & {\mathcal{N}}_n\left(\mathbf{0},{\mathit{\sigma }}_2^2{\mathit{H}}_2\left({\mathit{\theta }}_2\right)\right)\hfill \end{array}\right.\\ \hfill \left(III\right)& p(\mathit{\beta },\alpha ,{\mathit{\sigma }}^2,{\mathit{\psi }}^2,\mathit{\theta }),\end{array}$$

(6)

where $\alpha$ is an unknown parameter that relates both geographical patterns, ${\mathit{\sigma }}^{\mathbf{2}}=({\sigma }_1^2,{\sigma }_2^2)$ are the spatially structured species variabilities, ${\mathit{\psi }}^{\mathbf{2}}=({\mathit{\psi }}_1^2,{\mathit{\psi }}_2^2)$ are the non-structured spatial variabilities and $\mathit{\theta }=({\mathit{\theta }}_1,{\mathit{\theta }}_2)$ are the parameters of the correlation functions. Note that the last level requires again the elicitation of the priors and hyperpriors of the model.

With respect to the variance-covariance matrices for the abundance of the two species, we assume that

$$\begin{array}{c}{\mathbf{\Sigma }}_{{\mathit{Z}}_1}={\sigma }_1^2{\mathit{H}}_1\left({\mathit{\theta }}_1\right)+{\psi }_1^2{\mathit{I}}_n;\hfill \\ {\mathbf{\Sigma }}_{{\mathit{Z}}_2}={\sigma }_2^2{\mathit{H}}_2\left({\mathit{\theta }}_2\right)+\alpha {\sigma }_1^2{\mathit{H}}_1\left({\mathit{\theta }}_1\right)+{\psi }_2^2{\mathit{I}}_n,\hfill & \hfill \end{array}$$

(7)

where ${\mathit{H}}_j\left({\mathit{\theta }}_j\right)$ represents the correlation matrix between locations.

Among other possible ones, we propose here the use of the family of correlation matrices due to Matérn [31], which turns out to be a very flexible family that generalizes many of the widely used covariance models in spatial statistics. Its expression, giving the covariance between the values of a random field at locations separated by a distance $h>0$, can be parameterized as

$$C\left(h\right)=\frac{{\sigma }^2}{2^{\nu -1}\Gamma \left(\nu \right)}{\left(\kappa h\right)}^{\nu }{K}_{\nu }\left(\kappa h\right),$$

where $\kappa >0$ stands for a scaling parameter, ${\sigma }^2$ represents the marginal variance and finally $K_{\nu }$ is the modified Bessel function of the second kind and order $\nu >0$ ([32] section 9.6).

2.2. Inference and Prediction within INLA

Making inference and obtain predictions in unsampled locations using this kind of geostatistical models where an indexed continuous Gaussian Field (GF) is part of the model can be rather complicated. To alleviate this, Lindgren et al. [33] proposed an explicit link between latent Gaussian Markov Random Fields (GMRF) and GF with a Matérn covariance structure via a weak solution to a Stochastic Partial Differential Equation (SPDE). In this case, the spatial term can be reparametrised as $\mathit{W}\sim \mathcal{N}(0,{\mathit{Q}}^{-1}(\kappa ,\tau ))$, still a Gaussian distribution but depending now on two different parameters, $\kappa$ and $\tau$. The relationship between these new ones and the previous ones is that the effective range, the distance at which the correlation between two points is close to 0.1, is approximately $\varphi =\frac{\sqrt{8\nu }}{\kappa }$, while the variance is ${\sigma }_{\mathit{w}}^2=\frac{1}{4\pi {\kappa }^2{\tau }^2}$ [33]. Nevertheless, Krainski et al. [34] recommend a more intuitive parametrisation, in particular, using the marginal standard deviation and the range.

The final aim of converting the GF into a GRMF is that the latter ones have sparse precision matrices and so inference and prediction can be then performed via the Integrated Nested Laplace Approximation (INLA) [25], a computational alternative to the well-known Monte Carlo Markov Chain methods. Nowadays, INLA is considered as a well-established tool for Bayesian inference in several research fields [35]. The R-INLA is a R package to easily use INLA. For more details on INLA and its relationship with the SPDE approach and, more generally, for more details on spatial and spatio–temporal models within INLA, we refer the reader respectively Krainski et al. [34] and to Blangiardo and Cameletti [36] where practical examples and code guidelines are also provided.

As mentioned above, the final step in a hierarchical Bayesian approach is to elicit the priors and hyperpriors of the corresponding parameters and hyperparameters. A good option to do so is to use previous information based on available expert opinion. For those situations where there is no expert opinion, the best choice is to use non-informative prior distributions. In these latter case, as usual in the INLA context, we propose the use of normal vague priors with mean and precision for the regression coefficients and, in line with Fuglstad et al. [37], the use of Penalized Complexity priors (PC-priors) [38] for the ranges and the standard deviations of the spatial fields.

3. Simulation Study

In this section, we discuss a simulation study, the final aim being to illustrate the effectiveness of a coregionalized model and to emphasize the difference we would obtain by using a univariate model for each species. For the sake of simplicity, simulations do not include covariates.

3.1. Generation of the Simulated Dataset

In particular, we generate 40 realisations of 999 points from the following two related (coregionalized) Gaussian Random Fields, ${\mathit{Z}}_1$ and ${\mathit{Z}}_2$, over a regular grid (10 × 10) using the RandomFields package [39]:

$$\begin{array}{ccc}{\mathit{Z}}_1\hfill & =& -5+{\mathit{W}}_1+\mathcal{N}(0,0.1^2)\\ {\mathit{Z}}_2\hfill & =& 3+1.5\times {\mathit{W}}_1+{\mathit{W}}_2+\mathcal{N}(0,0.3^2)\end{array}$$

(8)

where ${\mathit{W}}_1$ and ${\mathit{W}}_2$ are realisations of two stationary isotropic covariance models belonging to the Matérn family, both with smoothness parameters $\nu =1$, with corresponding random field marginal variances of $0.5$ and 1 respectively; and also both with the same range of 4. Note that the non-structured error follows a Gaussian process with 0 mean and standard deviation $0.1$ for ${\mathit{Z}}_1$ and with standard deviation of $0.3$ for ${\mathit{Z}}_2$. Note that the intercepts have been set to $-1.5$ for ${\mathit{Z}}_1$ and 3 for ${\mathit{Z}}_2$, while the spatial relationship between the two random fields is produced by setting $\alpha =1.5$.

Once we have the 40 realizations, we randomly divide them into two subsamples, a train sample with 60% of the 999 simulated points, and a prediction sample with the remaining 40%. The train samples will be used to fit the models and estimate the effects, while the prediction sample will be used to evaluate the prediction accuracy of the models.

3.2. Fitting Univariate and Coregionalized Models

Fitting an univariate modelling consists of erroneously considering no interaction between the two random fields

$$\begin{array}{ccc}{\mathit{Z}}_1\hfill & =& {\beta }_{10}+{\mathit{W}}_1+{ϵ}_1\\ {\mathit{Z}}_2\hfill & =& {\beta }_{20}+{\mathit{W}}_2+{ϵ}_2,\end{array}$$

(9)

while fitting a coregionalized modelling considers the interaction between the two random fields

$$\begin{array}{ccc}{\mathit{Z}}_1\hfill & =& {\beta }_{10}+{\mathit{W}}_1+{ϵ}_1\\ {\mathit{Z}}_2\hfill & =& {\beta }_{20}+\alpha {\mathit{W}}_1+{\mathit{W}}_2+{ϵ}_2.\end{array}$$

(10)

Note that in both modellings, ${\mathit{W}}_1$ and ${\mathit{W}}_2$ are spatial random effects as in (6) and (7).

As above mentioned, inference on the parameters governing both modelings is performed via INLA and the SPDE approach, the main interest being the hyperparameters of both spatial random effects (range and spatial variance), and the coefficient that relates both effects. In line with Section 2, we use a PC-prior for the standard deviation of the latent effect, defined by P$(\sigma >2)=0.01$; a PC-prior for the range defined by P$(\varphi <0.5)=0.01$; and the default (non-informative) priors for the remaining parameters. These selections are based on the fact that we know the real value in advance (see [37], for hints about how to express the lack of previous information).

3.3. Results

The good fit of the coregionalized modeling can be appreciated in Figure 1, where we have represented a box-plot of the simulated posterior means of the parameter $\alpha$ (one for each realization). Note that the estimated values of the most important parameter in the coregionalization modeling are close to the real one.

Figure 2 compares the predictive capacity of both modelings through the predictive Mean Absolute Error (MAE) and Root Mean Square Error (RMSE) obtained for the univariate and coregionalized models in each realization. In general, coregionalized models show lower values than univariate models, both in terms of MAE and RMSE, indicating a better predictive capacity. Indeed, the median of the MAEs corresponding to the 40 univariate realizations (0.2140) is lower than the median for the coregionaliszed approach (0.2083), while the median of the RMSE corresponding the 40 univariate realizations (0.2705) is lower than the one in the coregionalized approach (0.2620).

The R-code of this simulated study is an open-source code that can be found at the following web page https://github.com/StMoEcGroup/Incorporating-biotic-information-in-SDMs-a-coregionalised-approach accessed on 11 January 2021.

4. Describing the Prey-Predator Interaction between European Anchovy and Hake

In this section, we present a practical application in the fishery ecology context. In particular, we consider data on a prey-predator example between the European anchovy (Engraulis encrasicolus) and one of its predator species, the European hake (Merluccius merluccius), both collected during trawling surveys in the Mediterranean sea (Figure 3).

4.1. Data Collection

The data came from the EU-funded survey project named MEDIterranean Trawl Survey (MEDITS) [40]. This survey is performed every spring (April to June) using a stratified sampling design based on 5 bathymetric strata (10–50, 51–100, 101–200, 201–500 and 501–700 m). Data from 2000 to 2016 were provided by the Instituto Español de Oceanografía (IEO, Spanish Oceanographic Institute). For each station, all individuals were identified and their weight (kg) was recorded. Estimates of biomass (kg/km${}^2$) were calculated with the standard swept area method for both the European hake and anchovy species. In order to ensure the normality theoretical assumptions of these two response variables, a logarithmic transformation was applied.

Information about the bathymetry was obtained from the European Marine Observation and Data Network (EMODnet, http://www.emodnet.eu/ accessed on 11 January 2021). The spatial resolution was 0.002 × 0.002 decimal degrees (20 m).

4.2. Coregionalized Model

The coregionalized model that tries to reflect the relationship between the European anchovy and hake species is given by:

$$\begin{array}{ccc}log\left(\mathrm{Hake}\right)\hfill & =& {\beta }_{10}+{\beta }_{11}\mathrm{Bathymetry}+{\mathit{W}}_1+{ϵ}_1\\ log\left(\mathrm{Anchovy}\right)\hfill & =& {\beta }_{20}+{\beta }_{21}\mathrm{Bathymetry}+\alpha {\mathit{W}}_1+{\mathit{W}}_2+{ϵ}_2\end{array}$$

(11)

where log(Hake) and log(Anchovy) represent the logarithmic transformation of the mean biomass of both species in each sampling station. The linear predictors, which indicate how parameters behave in space, are formed by: ${\beta }_{10}$ and ${\beta }_{20}$, i.e., the terms representing the intercepts of each variable; ${\beta }_{11}$ and ${\beta }_{21}$, that represent the regression coefficients of a covariate introduced into the model, the bathymetry of the seabed for each sampling station; and ${\mathit{W}}_1$ and ${\mathit{W}}_2$, that are the spatially structured effects. No interacting univariate models, i.e., $\alpha =0$, were also implemented.

In order to facilitate the interpretation, we used a linear relationship between the fish species and the bathymetry variable, but a non-linear relationship could also have been fitted by incorporating second-order random walk latent models [41].

As no other information was available, we used non-informative Gaussian prior distributions with zero mean and standard deviation of 100 for all the fixed effects. PC-priors were used to describe prior knowledge of hyperparameters of the spatial terms. These priors were set as follows: if the prior probability of the spatial range was smaller than 0.1, it was set at 0.05, if the probability of the spatial variance was larger than 3, it was also set at 0.05.

4.3. Results

The coregionalized model in (11) and its corresponding univariate version were compared with submodels including fewer terms using a model fit measure, in particular WAIC [42]. All submodels resulted in a worse model fit.

Results of the inference process for the univariate and for the coregionalized model are summarised in

Table 1. Point estimates (mean, median and standard deviation) along with 95% credible interval of the posterior distribution of the fixed effects and the hyperparameters of the univariate models for both the European anchovy (Engraulis encrasicolus) and the European hake (Merluccius merluccius).

		Mean	sd	${\mathit{q}}_{0.025}$	${\mathit{q}}_{0.5}$	${\mathit{q}}_{0.975}$
Hake
	${\beta }_{10}$	6.5634	0.2711	6.0059	6.5727	7.0698
	${\beta }_{11}$	−0.0039	0.0007	−0.0053	−0.0039	−0.0025
	${\varphi }_1$	0.154	0.027	0.107	0.152	0.214
	${\sigma }_1$	2.147	0.154	1.863	2.141	2.466
Anchovies
	${\beta }_{20}$	4.4756	0.4160	3.58831	4.4972	5.2391
	${\beta }_{21}$	−0.0025	0.0007	−0.0039	−0.0025	−0.0010
	${\varphi }_2$	0.614	0.166	0.371	0.585	1.015
	${\sigma }_2$	1.623	0.220	1.251	1.602	2.113

and 

Table 2. Point estimates (mean, median and standard deviation) along with 95% credible interval of the posterior distribution of the fixed effects and the hyperparameters of the coregionalized approach for both the European anchovy (Engraulis encrasicolus) and the European hake (Merluccius merluccius).

	Mean	sd	${\mathit{q}}_{0.025}$	${\mathit{q}}_{0.5}$	${\mathit{q}}_{0.975}$
${\beta }_{10}$	6.5258	0.2759	5.9600	6.5343	7.0434
${\beta }_{20}$	4.4426	0.3934	3.6091	4.4609	5.1692
${\beta }_{11}$	−0.0038	0.0007	−0.0052	−0.0038	−0.0024
${\beta }_{21}$	−0.0021	0.0008	−0.0036	−0.0021	−0.0006
${\varphi }_1$	0.161	0.028	0.113	0.158	0.223
${\sigma }_1$	2.148	0.152	1.864	2.143	2.463
${\varphi }_2$	0.575	0.154	0.350	0.548	0.948
${\sigma }_2$	1.544	0.194	1.212	1.527	1.971
$\alpha$	0.143	0.073	−0.003	0.144	0.285

, respectively. The association between European hake and anchovy is positive (i.e., mean $\alpha$ = 0.143), indicating that the geographical patterns of abundance of both species are related, and absorbs part of anchovies residual pattern (i.e., smaller ${\sigma }_2$ in the coregionalized model). This positive value is in line with [43], who concluded that predator–prey interactions may display positive or negative associations. The negative bathymetric effect was consistent across approaches, and suggests high abundances in shallow waters.

Figure 4 shows the mean of the spatial effect in the univariate (left map) and coregionalized models (middle map), as well as their difference (right map). As it can be seen, differences manifest in the transition zones between high and low abundances, providing a subtle improvement in the prediction of anchovy.

This improvement is also observed when analyzing the predictive capacity of both modelings through the predictive Mean Absolute Error (MAE) and Root Mean Square Error (RMSE) obtained for the univariate and coregionalized models (see

Table 3. Predictive capacity of Univariate and coregionalized models obtained by cross-validation of real 2016 data with the previous adjusted model.

Model	Measure	Value
Univariate	MAE	4212.35
Coregionalisation	MAE	3842.35
Univariate	RMSE	8422.54
Coregionalisation	RMSE	7697.88

). In particular, the model was fitted using data from 2000–2015, and the prediction was made for 2016 and compared with the observed values. Again, the coregionalized model shows lower values than the univariate model, both in terms of MAE and RMSE, indicating a better predictive capacity.

5. Concluding Remarks

There are still many issues to be tackled in species distribution modeling. Among them, this work has focused on the modeling of more than one species. We have revisited a spatial hierarchical Bayesian conditional coregionalized model to infer biotic associations in space. The main contribution has been to validate its behavior by means of a simulated example and a real application in the fisheries context.

Our results with simulated data indicate a better predictive capacity of the coregionalized approach. The coregionalized model successfully captured the association between species and the median MAE and RMSE scores were lower than those for the univariate model scores.

The practical application assessing the spatial distribution of prey, the European anchovy (Engraulis encrasicolus), and one of its predator species, European hake (Merluccius merluccius), in the Mediterranean sea has also shown a better behavior in the coregionalized model that includes biotic information. Main differences appear in the transition zones between high and low abundances, providing a subtle improvement in the prediction of anchovy. Moreover, the coregionalized model is better in terms of predictive capacity according to MAE and RMSE scores. Consequently, we conclude that this modeling could provide a better understanding of both species’ mesoscale ecology.

It is worth also noting that all this can be done thanks to the integrated nested Laplace approximation methodology and software that, jointly with the SPDE approach can help to minimize the computational burden while constituting a flexible tool in order to fit complex geostatistical models [44].

Finally, we would like to mention that this modeling could also include a temporal component that would allow us to analyze any possible temporal effect, in a similar way as in [23]. Moreover, it could also include other parametric or semiparametric constructions to reflect non-linear, autoregressive or more complex behaviors that would help to better describe the distribution of a species.