# Additive multivariate Gaussian processes for joint species distribution
modeling with heterogeneous data
Accepted for publication in Bayesian Analysis

Jarno<EMAIL_ADDRESS>[ Marcelo
<EMAIL_ADDRESS>[ Lari Veneranta [ Department of
Mathematics and Statistics and Organismal and Evolutionary Biology Research
Program, University of Helsinki,  Department of Mathematics and Statistics and
Department of Computer Science, University of Helsinki,  Natural Resources
Institute Finland, Finland,

(2019)

###### Abstract

Species distribution models (SDM) are a key tool in ecology, conservation and
management of natural resources. Two key components of the state-of-the-art
SDMs are the description for species distribution response along environmental
covariates and the spatial random effect that captures deviations from the
distribution patterns explained by environmental covariates. Joint species
distribution models (JSDMs) additionally include interspecific correlations
which have been shown to improve their descriptive and predictive performance
compared to single species models. However, current JSDMs are restricted to
hierarchical generalized linear modeling framework. Their limitation is that
parametric models have trouble in explaining changes in abundance due, for
example, highly non-linear physical tolerance limits which is particularly
important when predicting species distribution in new areas or under scenarios
of environmental change. On the other hand, semi-parametric response functions
have been shown to improve the predictive performance of SDMs in these tasks
in single species models.

Here, we propose JSDMs where the responses to environmental covariates are
modeled with additive multivariate Gaussian processes coded as linear models
of coregionalization. These allow inference for wide range of functional forms
and interspecific correlations between the responses. We propose also an
efficient approach for inference with Laplace approximation and
parameterization of the interspecific covariance matrices on the euclidean
space. We demonstrate the benefits of our model with two small scale examples
and one real world case study. We use cross-validation to compare the proposed
model to analogous semi-parametric single species models and parametric single
and joint species models in interpolation and extrapolation tasks. The
proposed model outperforms the alternative models in all cases. We also show
that the proposed model can be seen as an extension of the current state-of-
the-art JSDMs to semi-parametric models.

60G15,

60K35,

62P12,

linear model of coregionalization,

hierarchical model,

heterogeneous data,

spatial prediction,

model comparison,

Laplace approximation,

covariance transformation,

###### doi:

10.1214/19-BA1158

###### keywords:

[class=MSC]

###### keywords:

††volume: Authors’ accepted version

T1First available in Project Euclid: 3 June 2019: Available in Project Euclid:
https://projecteuclid.org/euclid.ba/1559548823

, and<EMAIL_ADDRESS>
##  1 Introduction

Species distribution models (SDMs) are key tools in the ecologists’ toolbox.
They have been widely used, among other applications, to study species habitat
preferences (Latimer et al., 2006), to improve identification and management
of conservation areas and natural resources (Kallasvuo et al., 2017) and to
evaluate species responses to environmental filtering under climate change
scenarios (Clark et al., 2014; Kotta et al., 2019). The main goals of
statistical inference in these contexts are to use species observations and
information on the associated environment to infer the relationship between
these two attributes and to predict over regions of unsampled locations to
build thematic species distribution maps (Gelfand et al., 2006; Elith and
Leathwich, 2009).

Species distribution modeling is, thus, directly related to inferring species’
responses to environmental factors (Latimer et al., 2006). This is
traditionally done using generalized linear or additive models as well as an
array of machine learning approaches such as regression trees or maximum
entropy modeling (Elith and Leathwich, 2009). These approaches model each
species separately and cannot account for species interactions nor shared
responses to the environment. However, species interaction with other species
is potentially as important factor as its response to environment. Moreover,
in many practical situations, data from species can be patchy or scarce in
which case sharing information between species can significantly improve
models’ predictive performance (Ovaskainen and Soininen, 2011; Clark et al.,
2014; Hui et al., 2013; Thorson et al., 2015; Hartmann et al., 2017). For
these reasons, joint species distribution models (JSDM) have gained increasing
attention in recent years (Warton et al., 2015).

Dunstan et al. (2013) and Hui et al. (2013) introduced species archetype
modeling where species’ responses to the environment are clustered into few
archetype models. Pollock et al. (2014) proposed to use the multivariate
probit regression model (Chib and Greenberg, 1998) to describe geographical
co-occurrence patterns between frogs and eucalyptus trees in Australia and
Clark et al. (2014) built a JSDM to infer richness and loss of species under
climate change scenarios. Thorson et al. (2015) introduced a spatial latent
factor model to predict spatial distribution of breeding birds and rock fish
communities and recently, Ovaskainen et al. (2017) introduced the hierarchical
modeling of species communities (HMSC) framework which includes detailed
description of interspecific correlations in covariate responses and spatial
random effects.

Current JSDMs rely on the classical framework of hierarchical generalized
linear models (GLMs). Even though this approach allows flexible modeling by
describing the randomness of response variables with different probabilistic
models (Nelder and Wedderburn, 1972), it is still limited by its parametric
assumptions. Hence, it may fail to accurately describe a species’ response to
environmental conditions (Vanhatalo et al., 2012; Kotta et al., 2019). Here,
we propose a semi-parametric JSDM model represented with multivariate Gaussian
processes (GPs). GPs are flexible semi-parametric regression models where the
regression function is estimated without restrictive parametric assumptions
about its form (O’Hagan, 1978; Rasmussen and Williams, 2006). Our model
integrates the main strengths of semi-parametric single species models
(Vanhatalo et al., 2012; Golding and Purse, 2016) and generalized linear model
based JSDMs.

First, we model the species responses to environmental covariates with
additive multivariate GPs. This allows us to capture wide range of nonlinear
responses and share information about these responses between species. Second,
due to the hierarchical model structure, the model can simultaneously
accommodate several kinds of outcome variables (observations/measurements) by
assigning different types of probabilistic models for them. This is important,
since it allows us to exploit different types of measurements which commonly
arise in real-case scenarios of multiple species surveys. Our presentation
focuses in the mostly used probabilistic models in practice, the Bernoulli
(Binomial) and Poisson with over-dispersion (Negative-Binomial). We also
propose an efficient computational approach build around Laplace method
(Golding and Purse, 2016; Vanhatalo et al., 2010). Lastly, we present a
structured cross-validation scheme to validate and compare models’ performance
in different kinds of prediction tasks.

This paper is organized as follows. In Section 2 we introduce a motivating
case study. In section 3, we introduce the additive multivariate GPs and in
Section 4 we discuss its predictive properties and introduce the computational
methods for inference. In Section 5, we illustrate the basic properties of the
model through two simple examples and introduce the specific case study model.
In Section 6 we present the case study results and we end by discussion and
conclusions in Section 7.

##  2 Motivating case study: coastal species distributions in the Gulf of
Bothnia

Figure 1: The study region, Gulf of Bothnia, and locations for species data.
The region is divided into five subregions (I-V) to be used in cross-
validation. The environmental conditions (described by the environmental
covariates) are heterogeneous between these regions (Veneranta et al., 2013)
which allows to test models’ extrapolation performance.

As a motivating example, we study spring distribution of four fish species and
three types of algae or macro-vegetation on the coastal region of the Gulf of
Bothnia in the northern Baltic Sea. The Gulf of Bothnia is a brackish water
basin between Sweden and Finland covering an area of approximately $600\times
120\ \mathrm{km}$. Its coastal areas play a central role in the ecosystem and
many Baltic fish stocks are dependent on the coastal regions for their
reproduction (Veneranta et al., 2013; Kallasvuo et al., 2017). Coastal zones
are also the most sensitive regions of the Baltic sea to both natural
variation and anthropogenic pressures (Reusch et al., 2018). Hence, these
areas are of central importance for conservation and there is need for
detailed knowledge on the distribution of coastal species and predictions
concerning their response to environmental changes.

###  2.1 Case study species and data

The studied fish species are the juvenile or adult three-spined stickleback
(_Gasterosteus aculeatus_) and nine-spined stickleback (_Pungitius pungitius_)
as well as larvae of whitefish (_Coregonus lavaretus_) and vendace (_Coregonus
albula_). Both whitefish and vendace are commercially important species and
the sticklebacks are one key fish fauna in the Gulf of Bothnia ecosystem
(Bergström et al., 2015). We treat the studied vegetation and algae in
functional group level comprising of diatomous algae, filamentous and macro-
vegetation. These are the dominating groups of benthic vegetation in shores
where larvae of Coregonids (whitefish and vendace) occur at early
developmental stages. In the scale of Gulf of Bothnia, their occurrence also
reflects the salinity, nutrient balance and wind exposure of studied area
(Veneranta et al., 2013).

The case study species reflect the large scale environmental gradients and the
changes in the Baltic Sea environment in last decades. Sticklebacks in the
Baltic Sea use the shallow coastal zone for reproduction (Bergström et al.,
2015) and high abundances of sticklebacks in the Baltic Sea have been
positively correlated with high occurrence of filamentous algae (Eriksson et
al., 2009, 2012). Coregonids in coastal areas prefer more oligotrophic waters
(Veneranta et al., 2013). Whitefish and vendace reproduce in stony reefs and
islets of Baltic Sea in late autumn and the larvae hatch at ice break-up in
spring. In sheltered areas the overwintering reeds (_Phragmites australis_)
dominate the macro-vegetation in spring. Diatomous algae consist of several
epiphytic species that have an early spring bloom at ice break up and settle
rapidly over bottom when light and water temperature increase (Busse and
Snoeijs, 2002). Filamentous algae in this study consist mostly of _Pilayella
sp_. It is a fast growing annual species that dominates the exposed shores at
Baltic Sea in early spring (Rönnberg and Bonsdorff, 2004).

The whole study area was divided into sampling subareas (Figure 1). Within
each sub-area, data were collected at several sampling sites distributed so
that they covered the range of most important environmental covariates in that
sub-area. Sampling was conducted in 2009-2011. At each site sampling was done
approximately one week after the ice break. In the scale of Gulf of Bothnia,
the ice break-up happens in a period of approximately one month. The species
data used in this work comprise of 160 distinct sampling sites. In 70 sites
(2010 data) all species were sampled and the rest of the sites comprise
samples for the fishes only (Figure 1). Fish samples comprise the number of
caught fish together with information on the sampling effort at each site. The
effort was measured as the volume of sampled water (m3). A more detailed
description of the data collection is provided by Veneranta et al. (2013) For
plant data, the bottom was photographed at 5-13 locations at depth of 30 cm
within a distance of 2 m and parallel to shore line. The occurrence of a
species was recorded at 16 uniformly distributed points in all these
photographs. The case study plants can grow over each others so that presence
of one plant at a spatial location does not exclude the presence of other
plants. The whitefish and vendace data were previously analyzed by Veneranta
et al. (2013) and Vanhatalo et al. (2012) but the data on other species are
unpublished.

###  2.2 Environmental covariates

Variable | Description | Resolution [m]
---|---|---
Openness | The average openness and exposure to winds | 300
Distance to deep | Distance to 20 m depth curve | 200
Sandy bottom index | Area weighted distance to the sandy shores | 90
Ice breakup week | The end date (weeks) of ice cover in 2009 | 1852
Chlorophyll-a | Chlorophyll-a concentration | 1852
River | Distance to the nearest river mouth | 150
Bottom class | Bottom classification to 6 categories | 200
Winter salinity | Length of ice winter | 10,000
Table 1:  Environmental covariates used in the case study and their original
resolution. For this study, the resolution of all the environmental covariates
was scaled to 300m.

Gulf of Bothnia hosts rich variety of environmental conditions. Coastal areas
are affected by inflows from land as well as shallow and complex topography.
In the scale of Gulf of Bothnia, there is a gradient in river influence,
salinity, temperature and length of ice cover period from north to south. We
used seven real-valued and one categorical abiotic environmental covariates
that were available throughout the study area as raster maps. Each species
observation was combined with covariates from the raster cell where the sample
location fell in. These are summarized in table 1 and described and motivated
in detail by Veneranta et al. (2013). Due its fresh water origin vendace is _a
priori_ sensitive to even small changes in salinity levels experienced in the
Gulf of Bothnia whereas the other species respond to salinity only in the
Baltic Sea scale. Hence, we used all covariates only for vendace but excluded
the winter salinity from other species.

##  3 Hierarchical multivariate species distribution model

We start our model building following the generic hierarchical structure
similar to the one presented by, e.g., Wikle (2003), Cressie and Wikle (2011)
and Banerjee et al. (2015)

$\displaystyle\mathrm{[data|process,parameters]:}\hskip 5.69046pt$
$\displaystyle\pi_{\operatorname{\boldsymbol{Y}}}(\operatorname{\boldsymbol{y}}(\operatorname{\mathbf{x}},\operatorname{\mathbf{s}})|\operatorname{\mathbf{f}}(\operatorname{\mathbf{x}},\operatorname{\mathbf{s}}),\operatorname{\boldsymbol{\eta}})$
$\displaystyle\mathrm{[process|parameters]:}\hskip 5.69046pt$
$\displaystyle\pi_{\operatorname{\boldsymbol{f}}}(\operatorname{\mathbf{f}}(\operatorname{\mathbf{x}},\operatorname{\mathbf{s}})|\operatorname{\boldsymbol{\theta}})$
$\displaystyle\mathrm{[parameters]:}\hskip 5.69046pt$
$\displaystyle\pi(\operatorname{\boldsymbol{\eta}},\operatorname{\boldsymbol{\theta}})$
(3.1)

where the first layer of hierarchy is the probabilistic model which defines
the conditional distribution for multivariate data
$\operatorname{\boldsymbol{Y}}(\operatorname{\mathbf{x}},\operatorname{\mathbf{s}})$,
at spatial location $\operatorname{\mathbf{s}}$ with associated covariates
$\operatorname{\mathbf{x}}$, given the model parameters
$\operatorname{\boldsymbol{\eta}}$ and the multivariate latent process
$\operatorname{\boldsymbol{f}}(\operatorname{\mathbf{x}},\operatorname{\mathbf{s}})$.
The second layer defines the prior distribution for the latent process given
the process parameters $\operatorname{\boldsymbol{\theta}}$, and the third
layer defines the prior distribution for all unknown parameters.

Let $j\in\mathcal{J}=\\{1,\cdots,J\\}$ be the species index set and $J$ the
total number of species in the study. Denote by
$\operatorname{\boldsymbol{Y}}(\operatorname{\mathbf{x}},\operatorname{\mathbf{s}})^{T}=[Y_{1}\cdots
Y_{J}]$ the $J$-variate random vector with components
$Y_{j}=Y_{j}(\operatorname{\mathbf{x}}_{j},\operatorname{\mathbf{s}})$ related
to the $j^{th}$ species at spatial location
$\operatorname{\mathbf{s}}^{T}=[s_{1}\ s_{2}]$ (coordinates) under the
influence of environmental covariates
$\operatorname{\mathbf{x}}_{j}^{T}=[x_{j,1}\cdots x_{j,c_{j}}]$ where $c_{j}$
is the number of environmental covariates for the species $j$. We denote by
$\operatorname{\mathbf{x}}$ the full vector of covariates and $\mathcal{X}$
the covariate space. The species specific covariates
$\operatorname{\mathbf{x}}_{j}$ are sub-vectors of $\operatorname{\mathbf{x}}$
such that $\operatorname{\mathbf{x}}_{j}\in\mathcal{X}_{j}$ where
$\mathcal{X}_{j}\subset\mathcal{X}$ and $\mathcal{X}_{j}$ is
$c_{j}$-dimensional. Here, we assume that given a multivariate latent process
$\operatorname{\boldsymbol{f}}(\operatorname{\mathbf{x}},\operatorname{\mathbf{s}})=[f_{1}(\operatorname{\mathbf{x}}_{1},\operatorname{\mathbf{s}})\cdots
f_{J}(\operatorname{\mathbf{x}}_{J},\operatorname{\mathbf{s}})]^{T}:\mathcal{X}\times\mathbb{R}^{2}\rightarrow\mathbb{R}^{J}$,
the probabilistic model for
$\operatorname{\boldsymbol{Y}}=\operatorname{\boldsymbol{Y}}(\operatorname{\mathbf{x}},\operatorname{\mathbf{s}})$
factorizes as follows

$\pi_{\operatorname{\boldsymbol{Y}}}(\operatorname{\boldsymbol{y}}(\operatorname{\mathbf{x}},\operatorname{\mathbf{s}})|\operatorname{\mathbf{f}}(\operatorname{\mathbf{x}},\operatorname{\mathbf{s}}),\operatorname{\boldsymbol{\eta}})=\prod\limits_{j=1}^{J}\pi_{Y_{j}}(y_{j}|\operatorname{\mathrm{f}}_{j},\eta_{j})$
(3.2)

where
$\operatorname{\mathrm{f}}_{j}=\operatorname{\mathrm{f}}_{j}(\operatorname{\mathbf{x}}_{j},\operatorname{\mathbf{s}})$
is fixed but an unknown realization of the $j$’th process with covariates
$\operatorname{\mathbf{x}}_{j}$ at the spatial location
$\operatorname{\mathbf{s}}$. $\eta_{j}$ is the vector of parameters of the
probabilistic model $\pi_{Y_{j}}$ for the species $j$ and
$\operatorname{\boldsymbol{\eta}}=[\eta_{1}\cdots\eta_{J}]^{T}$. In general
any probabilistic model for data could be used and the observation models
should be chosen according to the assumed sampling process. Here, we consider
in detail the Binomial and the Negative-Binomial (over-dispersed Poisson)
models used in our case study (Section 2). These models can be seen as
observation models resulting from inhomogeneous point process model for
species abundance (Gelfand et al., 2006; Warton and Shepherd, 2010).

In practice, each sampling site consists of a small area where the
observations are made. In our case study, the area covered by a sampling site
is so small compared to the whole study region that the latent function is
practically constant in each sampling site. The model for the count of species
presences at uniformly distributed points in a sampling site (plants in our
case study) is then Binomial given by (Gelfand et al., 2006)

$\pi_{Y_{B}}(y|f(\operatorname{\mathbf{x}},\operatorname{\mathbf{s}}),z_{B}(\operatorname{\mathbf{s}}))={{z_{B}(\operatorname{\mathbf{s}})}\choose{y}}p(f(\operatorname{\mathbf{x}},\operatorname{\mathbf{s}}))^{y}[1-p(f(\operatorname{\mathbf{x}},\operatorname{\mathbf{s}}))]^{z_{B}(\operatorname{\mathbf{s}})-y}I_{\\{0,\ldots,z_{B}(\operatorname{\mathbf{s}})\\}}(y)$
(3.3)

where $z_{B}(\operatorname{\mathbf{s}})$ is the total number of observation
points in the sampling site at location $\operatorname{\mathbf{s}}$ and
$p(f(\operatorname{\mathbf{x}},\operatorname{\mathbf{s}}))$ is the success
probability, which will be modeled through the logit here. The natural model
for the number of individuals in a closed area or volume in a sampling site
(fish in our case study) is Poisson which we extend to over-dispersed Poisson
using the Negative-Binomial given by (see Liu and Vanhatalo, 2018, Section 7.1
for details)

$\pi_{Y_{N}}(y|f(\operatorname{\mathbf{x}},\operatorname{\mathbf{s}}),z_{N}(\operatorname{\mathbf{s}}),r)=\dfrac{\Gamma(r+y)}{y!\Gamma(r)}\left(\dfrac{r}{r+m(\operatorname{\mathbf{x}},\operatorname{\mathbf{s}})}\right)^{r}\left(\dfrac{m(\operatorname{\mathbf{x}},\operatorname{\mathbf{s}})}{r+m(\operatorname{\mathbf{x}},\operatorname{\mathbf{s}})}\right)^{y}I_{\mathbb{N}_{0}}(y)$
(3.4)

where $m(\operatorname{\mathbf{x}},\operatorname{\mathbf{s}})$ $=$
$z_{N}(\operatorname{\mathbf{s}})\exp[f(\operatorname{\mathbf{x}},\operatorname{\mathbf{s}})]=E(Y_{N}|f(\operatorname{\mathbf{x}},\operatorname{\mathbf{s}}),z_{N}(\operatorname{\mathbf{x}}),r)$
and $r$ is the over-dispersion parameter. The latent process $f(\cdot)$
corresponds to the logarithm of a species (relative) density and
$z_{{}_{N}}(\operatorname{\mathbf{s}})$ is the sampled volume of water in the
sampling site at location $\operatorname{\mathbf{s}}$. For this
parameterization
$V(Y_{N}|f(\operatorname{\mathbf{x}},\operatorname{\mathbf{s}}),z_{N}(\operatorname{\mathbf{x}}),r)$
$=$ $m(\operatorname{\mathbf{x}},\operatorname{\mathbf{s}})$ $+$
$m(\operatorname{\mathbf{x}},\operatorname{\mathbf{s}})^{2}/r$. Increasing $r$
decreases the variance and when $r\rightarrow\infty$, the model approaches the
Poisson distribution.

###  3.1 Univariate additive latent Gaussian process

First, we assume that marginally for any species $j$, the process model is
given by

$f_{j}(\operatorname{\mathbf{x}}_{j},\operatorname{\mathbf{s}})=\beta_{j,0}+h_{j}(\operatorname{\mathbf{x}}_{j})+\varepsilon_{j}(\operatorname{\mathbf{s}})$
(3.5)

where $\beta_{j,0}$ is the offset weight with distribution
$\beta_{j,0}|\sigma^{2}_{j,0}\sim\mathcal{N}(0,\sigma^{2}_{j,0})$ and
$h_{j}:\mathcal{X}_{j}\rightarrow\mathbb{R}$ is a _predictor function_ of
environmental covariates. The predictor function is assumed to be additive
over the covariates,
$h_{j}(\operatorname{\mathbf{x}}_{j})=\sum_{r=1}^{c_{j}}h_{j,r}(x_{j,r})$, and
each additive term is given an independent zero mean GP prior, so that

$h_{j}(\operatorname{\mathbf{x}}_{j})|\operatorname{\boldsymbol{\theta}}_{j}\sim\mathcal{GP}\left(0,\sum_{r=1}^{c_{j}}k_{h_{j,r}}(x_{j,r},x^{\prime}_{j,r};\operatorname{\boldsymbol{\theta}}_{j,r})\right)$
(3.6)

where
$k_{h_{j,r}}(x_{j,r},x^{\prime}_{j,r};\operatorname{\boldsymbol{\theta}}_{j,r})=\operatorname{\mathrm{Cov}}(h_{j,r}(x_{j,r}),h_{j,r}(x^{\prime}_{j,r})|\operatorname{\boldsymbol{\theta}}_{j,r})$
is a covariance function with parameter
$\operatorname{\boldsymbol{\theta}}_{j,r}$ and
$\operatorname{\boldsymbol{\theta}}_{j}$ $=$
$[\operatorname{\boldsymbol{\theta}}_{j,1}^{T}\cdots\operatorname{\boldsymbol{\theta}}_{j,c_{j}}^{T}]^{T}$.
For example, with $k_{h_{j,r}}(x_{j,r},x^{\prime}_{j,r};\theta_{j,r})$ $=$
$x_{j,r}$ $x^{\prime}_{j,r}\sigma^{2}_{j,r}$ and $\theta_{j,r}$ $=$
$\sigma^{2}_{j,r}$, the predictor function corresponds to linear model
$h_{j,r}(x_{j,r})=x_{j,r}\beta_{j,r}$ where $\beta_{j,r}\sim
N(0,\sigma^{2}_{j,r})$ (Rasmussen and Williams, 2006). With other choices of
covariance functions we can model non-linear responses in which case the model
can be seen as an alternative to the traditional generalized additive models
(GAMs, Hastie and Tibishirani, 1986). The general form of an additive GP prior
would include also joint effects of covariates (Duvenaud et al., 2011) but
this is not considered here.

The term $\varepsilon_{j}(\operatorname{\mathbf{s}})$ is a spatial GP,

$\varepsilon_{j}(\operatorname{\mathbf{s}})|\sigma^{2}_{j},\operatorname{\boldsymbol{\ell}}_{j}\sim\mathcal{GP}\big{(}0,k_{\epsilon_{j}}(\operatorname{\mathbf{s}},\operatorname{\mathbf{s}}^{\prime};\operatorname{\boldsymbol{\ell}}_{j},\sigma^{2}_{j})\big{)}$
(3.7)

where
$k_{\epsilon_{j}}(\operatorname{\mathbf{s}},\operatorname{\mathbf{s}}^{\prime};\operatorname{\boldsymbol{\ell}}_{j},\sigma^{2}_{j})$
is a spatial covariance function with variance $\sigma^{2}_{j}$ and range
parameter $\operatorname{\boldsymbol{\ell}}_{j}$. When we consider independent
models for each species, that is the traditional single species SDMs (SSDMs),
the processes $h_{j,r}$ and $\epsilon_{j}$ are mutually independent among all
species. The model outlined this far is similar to the GP based species
distribution models used by, e.g., Vanhatalo et al. (2012), Kallasvuo et al.
(2017) and Golding and Purse (2016). Next, we introduce models which consider
dependency.

###  3.2 Additive multivariate Gaussian process priors

In order to account for possible species interdependence, we first include
spatial species interaction into the model through the linear model of
coregionalization (LMC) (Mardia and Goodall, 1993; Gelfand et al., 2004).
Write
$\operatorname{\boldsymbol{\varepsilon}}(\operatorname{\mathbf{s}})^{T}=[\varepsilon_{1}(\operatorname{\mathbf{s}})\cdots\varepsilon_{J}(\operatorname{\mathbf{s}})]$
and assume that
$\operatorname{\boldsymbol{\varepsilon}}(\operatorname{\mathbf{s}})$ has LMC
covariance structure with species specific correlation functions
$\tilde{k}_{\epsilon_{j}}(\cdot,\cdot\
;\operatorname{\boldsymbol{\ell}}_{j})$. The covariance structure of the LMC
with interspecific spatial dependence is then,

$\operatorname{\mathrm{Cov}}\big{(}\varepsilon_{j}(\operatorname{\mathbf{s}}),\varepsilon_{j^{\prime}}(\operatorname{\mathbf{s}}^{\prime})|\Sigma_{\epsilon},\operatorname{\boldsymbol{\ell}}\big{)}=\sum_{l=1}^{J}u_{\epsilon,l}(j,j^{\prime})\tilde{k}_{\epsilon_{l}}(\operatorname{\mathbf{s}},\operatorname{\mathbf{s}}^{\prime};\operatorname{\boldsymbol{\ell}}_{l})$
(3.8)

with
$\operatorname{\boldsymbol{\ell}}^{T}=[\operatorname{\boldsymbol{\ell}}_{1}^{T}\cdots\operatorname{\boldsymbol{\ell}}_{J}^{T}]$
and $u_{\epsilon,l}(j,j^{\prime})$ the entry $(j,j^{\prime})$ of
$U_{\epsilon,l}=L_{\epsilon,l}L^{T}_{\epsilon,l}$ where $L_{\epsilon,l}$ is
the $l^{th}$ column of the Cholesky decomposition $L_{\epsilon}$ of the
coregionalization matrix $\Sigma_{\epsilon}$, that is matrix of interspecific
correlations between spatial GP. Hence, the vector-valued latent process
$\operatorname{\boldsymbol{f}}(\operatorname{\mathbf{x}},\operatorname{\mathbf{s}})^{T}=[f_{1}(\operatorname{\mathbf{x}}_{1},\operatorname{\mathbf{s}})\cdots
f_{J}(\operatorname{\mathbf{x}}_{J},\operatorname{\mathbf{s}})]$ unconditional
on $\beta_{1,0}$,$\ldots$,$\beta_{J,0}$ follows a multivariate GP which we
denote as

$\operatorname{\boldsymbol{f}}(\operatorname{\mathbf{x}},\operatorname{\mathbf{s}})|\Lambda_{1}\sim\mathcal{MGP}\bigg{(}\operatorname{\boldsymbol{0}},\Sigma_{0}+\sum_{r=1}^{c}k_{h_{r}}(\operatorname{\mathbf{x}}_{r},\operatorname{\mathbf{x}}^{\prime}_{r};\operatorname{\boldsymbol{\theta}})+\sum_{j=1}^{J}\tilde{k}_{\epsilon_{j}}(\operatorname{\mathbf{s}},\operatorname{\mathbf{s}}^{\prime};\operatorname{\boldsymbol{\ell}}_{j})U_{\epsilon,j}\bigg{)}$
(3.9)

where
$\Lambda_{1}=(\Sigma_{0},\Sigma_{\epsilon},\operatorname{\boldsymbol{\theta}},\operatorname{\boldsymbol{\ell}})$,
$\Sigma_{0}=\operatorname{\mathrm{diag}}(\sigma^{2}_{1,0},\ldots,\sigma^{2}_{J,0})$,
$\operatorname{\boldsymbol{\theta}}^{T}=[\operatorname{\boldsymbol{\theta}}_{1}^{T}\cdots\operatorname{\boldsymbol{\theta}}_{J}^{T}]$
and
$k_{h_{r}}(\operatorname{\mathbf{x}}_{r},\operatorname{\mathbf{x}}^{\prime}_{r};\operatorname{\boldsymbol{\theta}})=\operatorname{\mathrm{diag}}\big{(}k_{h_{1,r}}(x_{1,r},x^{\prime}_{1,r};\operatorname{\boldsymbol{\theta}}_{1,r}),\cdots,k_{h_{J,r}}(x_{J,r},x^{\prime}_{J,r};\operatorname{\boldsymbol{\theta}}_{J,r})\big{)}$.
If the predictor functions, $h_{j,r}$, were linear, the prior (3.9) would
correspond to traditional multivariate spatial model with independent linear
effects and spatial LMC (Gelfand et al., 2004).

We extend (3.9) to an _additive multivariate GP prior_ where the dependence
between the species specific additive predictor functions is modeled with LMC.
We demonstrate this with a model where the set of covariates is equal for all
species, that is, $\mathrm{dim}(\mathcal{X}_{j})=c$ and
$\operatorname{\mathbf{x}}_{j}=\operatorname{\mathbf{x}}$ for all $j$. Then
the model is written as

$\displaystyle\operatorname{\boldsymbol{f}}(\operatorname{\mathbf{x}},\operatorname{\mathbf{s}})|\Lambda_{2}$
$\displaystyle\sim\mathcal{MGP}\bigg{(}\operatorname{\boldsymbol{0}},\Sigma_{0}+\sum_{r=1}^{c}\sum_{j=1}^{J}\tilde{k}_{h_{{}_{j,r}}}(x_{r},x^{\prime}_{r};\operatorname{\boldsymbol{\theta}}_{j,r})U_{h_{{}_{r}},j}+\sum_{j=1}^{J}\tilde{k}_{\epsilon_{j}}(\operatorname{\mathbf{s}},\operatorname{\mathbf{s}}^{\prime};\operatorname{\boldsymbol{\ell}}_{j})U_{\epsilon,j}\bigg{)}$
(3.10)

where $\Lambda_{2}=(\Lambda_{1},\Sigma_{h_{1}},\dots,\Sigma_{h_{c}})$ and
$\Sigma_{h_{r}}$ is the interspecific covariance matrix between the species
specific predictor functions of $r$’th covariate,
$\tilde{k}_{h_{j,r}}(\cdot,\cdot;\operatorname{\boldsymbol{\theta}}_{j,r})$ is
a correlation function related to the predictor function $h_{j,r}$,
$U_{{h_{r}},j}=L_{h_{r},j}L^{T}_{h_{r},j}$ and $L_{{h_{r}},j}$ is the $j^{th}$
column of the Cholesky decomposition $L_{h_{r}}$ of $\Sigma_{h_{r}}$. When the
set of covariates is not the same for all species, covariate specific
predictive functions are correlated only between species that share those same
covariates.

A JSDM with multivariate GP prior (3.10) can be interpreted as extension of GP
based SSDMs (Vanhatalo et al., 2012; Golding and Purse, 2016) to joint species
modeling similarly as done with the generalized linear model based JSDMs
(Pollock et al., 2014; Ovaskainen et al., 2017). The enhanced inferential
ability of the multivariate additive GP compared to univariate GP models lies
in its capability to infer similarity (dissimilarity) in species specific
responses to covariates and in the spatial random effect. The
similarity/dissimilarity of responses of two species $j$ and $j^{\prime}$
along a covariate $r$ is indicated by a positive/negative value for
correlation and hence, examining the LMC covariance matrices,
$\Sigma_{h_{r}}$, can provide new insight to species to species associations.

###  3.3 Marginally uniform priors for correlation parameters

Here, we define the prior for coregionalization covariance matrices but define
the covariance functions and rest of the priors in Section 5. A standard
choice for prior for correlation matrices is the inverse Wishart distribution.
However, if there is no prior information on interspecific correlations, we
follow Hartmann et al. (2017) and suggest to use marginally uniform priors for
the LMC covariance matrices $\Sigma_{\epsilon}$ and $\Sigma_{h_{r}}$. These
can be achieved by the distribution of Barnard et al. (2000) and Tokuda et al.
(2012). Let $\operatorname{\mathcal{P}}$ be a correlation matrix with elements
$\rho_{j,j^{\prime}}$. A prior distribution that is marginally noninformative
for the correlations, that is, the marginal distribution for every
$\rho_{j,j^{\prime}}$ is uniform over $(-1,1)$, is achieved with the
distribution

$\pi(\operatorname{\mathcal{P}}|v)=\dfrac{\Gamma(\frac{v}{2})^{J}}{\Gamma_{J}(\frac{v}{2})}|\operatorname{\mathcal{P}}|^{\frac{1}{2}(v-1)(J-1)-1}\prod_{j=1}^{J}|\operatorname{\mathcal{P}}_{jj}|^{-\frac{v}{2}}\mathds{1}_{(0,\infty)}(\det\operatorname{\mathcal{P}})$
(3.11)

when $v=J-1$. Here, $\Gamma_{J}(\cdot)$ is the multivariate gamma function and
matrix $\operatorname{\mathcal{P}}_{jj}$ is obtained by removing the $j^{th}$
column and $j^{th}$ row of $\operatorname{\mathcal{P}}$. When $v$ increases,
the probability density (3.11) becomes increasingly concentrated around the
origin.

##  4 Posterior inference and predictive properties

Given a set of species observations at locations
$\operatorname{\mathbf{s}}_{i_{j}},i_{j}=1,\ldots,n_{j}$, respectively for
each species $j=1,\ldots,J$, the likelihood can be written as

$\pi(\operatorname{\mathbf{y}}|\operatorname{\mathbf{f}},\operatorname{\boldsymbol{\eta}})=\prod\limits_{j=1}^{J}\displaystyle\prod\limits_{i_{j}=1}^{n_{j}}\pi_{j}(y_{j,i_{j}}|\operatorname{\mathrm{f}}_{j,i_{j}},\eta_{j})$
(4.1)

where $y_{j,i_{j}}$ is the $i_{j}$’th observation of species $j$ at the
$i_{j}$’th spatial location $\operatorname{\mathbf{s}}^{T}_{i_{j}}$ associated
with a set of covariates
$\operatorname{\mathbf{x}}_{i_{j}}\in\mathcal{X}_{j}$. The observed vector
$\operatorname{\mathbf{y}}=[\operatorname{\mathbf{y}}_{1}^{T}\cdots\operatorname{\mathbf{y}}_{J}^{T}]^{T}$
with $\operatorname{\mathbf{y}}_{j}^{T}=[y_{j,1}\cdots y_{j,n_{j}}]$ collects
all the species observations and
$\operatorname{\mathbf{f}}=[\operatorname{\mathbf{f}}_{1}^{T}\ldots\operatorname{\mathbf{f}}_{J}^{T}]^{T}$,
where
$\operatorname{\mathbf{f}}^{T}_{j}=[\operatorname{\mathrm{f}}_{j,1}\cdots\operatorname{\mathrm{f}}_{j,n_{j}}]$,
collects the corresponding latent variables. Hence, the likelihood factorizes
over the latent variables and the inference can be done similarly as with
univariate GP models. Markov chain Monte Carlo (MCMC) would provide exact
answers in the limit of large sample size but deterministic approximations
such as Laplace approximation or expectation propagation algorithm have also
been shown to provide accurate approximate inference for univariate GP models
with much lower computational time (Nickisch and Rasmussen, 2008; Rue and
Marino, 2009; Vanhatalo et al., 2010). In order to study the properties of the
model, we will examine the Laplace approximation for the posterior of latent
variables conditional on the hyperparameters.

###  4.1 Posterior predictive inference conditional on hyperparameters

####  4.1.1 Posterior of latent variables

The Laplace’s method approximates the conditional posterior of the latent
function values
$\operatorname{\mathbf{f}}_{\ast}=\operatorname{\mathbf{f}}(\operatorname{\mathbf{s}}_{\ast},\operatorname{\mathbf{x}}_{\ast})$
at the spatial location $\operatorname{\mathbf{s}}_{\ast}$ with covariates
$\operatorname{\mathbf{x}}_{\ast}$ as111In case of Gaussian observation model,
this equals to the true posterior density of
$\operatorname{\mathbf{f}}|\operatorname{\mathbf{y}},\operatorname{\boldsymbol{\eta}},\Lambda$.

$\displaystyle\pi(\operatorname{\mathbf{f}}_{\ast}|\operatorname{\mathbf{y}},\operatorname{\boldsymbol{\eta}},\Lambda)\approx\mathcal{N}\big{(}\operatorname{\mathbf{f}}_{\ast}|C_{\ast,\operatorname{\mathbf{f}}}C^{-1}\hat{\operatorname{\mathbf{f}}},C_{\ast}-C_{\ast,\operatorname{\mathbf{f}}}(C^{-1}+\operatorname{\mathrm{W}})^{-1}C_{\operatorname{\mathbf{f}},\ast}\big{)}$
(4.2)

where $\hat{\operatorname{\mathbf{f}}}$ is the (conditional) maximum a
posterior (MAP) estimate of latent variables,

$\displaystyle\hat{\operatorname{\mathbf{f}}}=$
$\displaystyle\operatorname*{arg\,max}_{\operatorname{\mathbf{f}}\in\mathbb{R}^{\sum_{j}n_{j}}}\
\sum\limits_{j=1}^{J}\sum\limits_{i_{j}=1}^{n_{j}}\log\pi_{j}(y_{j,i_{j}}|\operatorname{\mathrm{f}}_{j,i_{j}},\eta_{j})+\log\mathcal{N}(\operatorname{\mathbf{f}}|\operatorname{\boldsymbol{0}},C)$
(4.3)

and $\operatorname{\mathrm{W}}$ is the Hessian matrix of the negative log-
likelihood function evaluated at $\hat{\operatorname{\mathbf{f}}}$. Here, $C$
is the prior covariance matrix of $\operatorname{\mathbf{f}}$,
$C_{\ast,\operatorname{\mathbf{f}}}$ is the (prior) covariance matrix between
elements of $\operatorname{\mathbf{f}}_{\ast}$ and
$\operatorname{\mathbf{f}}$. $C_{\ast}$ is the prior covariance of
$\operatorname{\mathbf{f}}_{\ast}$. In case of multivariate additive GP (3.10)
the prior covariance matrix is given by

$\displaystyle C=$
$\displaystyle\left[\begin{array}[]{ccc}\sigma_{0}^{2}\mathrm{J}_{n_{1}}&\cdots&0\\\
\vdots&\ddots&\vdots\\\
0&\cdots&\sigma_{0}^{2}\mathrm{J}_{n_{J}}\end{array}\right]+\sum_{r=1}^{c}\sum_{j=1}^{J}\left[\begin{array}[]{ccc}u_{h_{{}_{r}},j}(1,1)\\{\tilde{K}_{h_{j,r}}\\}_{1,1}&\cdots&u_{h_{{}_{r}},j}(1,J)\\{\tilde{K}_{h_{j,r}}\\}_{1,J}\\\
\vdots&\ddots&\vdots\\\
u_{h_{{}_{r}},j}(J,1)\\{\tilde{K}_{h_{j,r}}\\}_{J,1}&\cdots&u_{h_{{}_{r}},j}(J,J)\\{\tilde{K}_{h_{j,r}}\\}_{J,J}\end{array}\right]$
(4.10)
$\displaystyle+\sum_{j=1}^{J}\left[\begin{array}[]{ccc}u_{\epsilon,j}(1,1)\\{\tilde{K}_{\epsilon_{j}}\\}_{1,1}&\cdots&u_{\epsilon,j}(1,J)\\{\tilde{K}_{\epsilon_{j}}\\}_{1,J}\\\
\vdots&\ddots&\vdots\\\
u_{\epsilon,j}(J,1)\\{\tilde{K}_{\epsilon_{j}}\\}_{J,1}&\cdots&u_{\epsilon,j}(J,J)\\{\tilde{K}_{\epsilon_{j}}\\}_{J,J}\end{array}\right]$
(4.14)

where $\mathrm{J}_{n}$ is an $n\times n$ matrix full of ones
$u_{\epsilon,l}(j,j^{\prime})$ and $u_{h_{r},l}(j,j^{\prime})$ are the entry
$(j,j^{\prime})$ of $U_{\epsilon,l}$ and $U_{h_{r},l}$ (see (3.8)), and
$\\{\tilde{K}_{\epsilon_{l}}\\}_{j,j^{\prime}}$ and
$\\{\tilde{K}_{h_{l,r}}\\}_{j,j^{\prime}}$ are correlation matrices with
elements
$[\\{\tilde{K}_{\epsilon_{l}}\\}_{j,j^{\prime}}]_{i_{j},i_{j^{\prime}}}=\tilde{k}_{\epsilon,l}(\operatorname{\mathbf{s}}_{i_{j}},\operatorname{\mathbf{s}}_{i_{j^{\prime}}}\
;\operatorname{\boldsymbol{\ell}}_{l})$ and
$[\\{\tilde{K}_{h_{l,r}}\\}_{j,j^{\prime}}]_{i_{j},i_{j^{\prime}}}=\tilde{k}_{h_{l,r}}(x_{i_{j},r},x_{i_{j^{\prime}},r}|\operatorname{\boldsymbol{\theta}}_{l,r})$.
The other correlation matrices are constructed similarly. If data on all
species are available at all sampling locations, the covariance matrix reduces
to Kronecker product similarly as in LMC models by Mardia and Goodall (1993)
and Gelfand et al. (2004), so that $C=\Sigma_{0}\otimes
J_{n}+\sum_{r=1}^{c}\sum_{j=1}^{J}U_{h_{{}_{j,r}},j}\otimes\tilde{K}_{h_{j,r}}+\sum_{j=1}^{J}U_{\epsilon,j}\otimes\tilde{K}_{\epsilon_{j}}$
where $n_{j}=n$ for all $j=1,\ldots,J$. Since, in general, the second and
third term of $C$ are full matrices, it can be seen from (4.2) and (4.3) that
the posterior of $\operatorname{\mathbf{f}}$ and the posterior predictive
distribution of $\operatorname{\mathbf{f}}_{\ast}$ are affected by
observations of all species at any spatial locations.

####  4.1.2 The marginal posterior of additive terms

Typically, we are also interested in species’ responses along covariates
encoded by the additive predictor functions. Their marginal posteriors,
conditional on hyperparameters, are analogous to (4.2). The matrices
$C_{\ast}$ and $C_{\ast,\operatorname{\mathbf{f}}}$ are only constructed using
the correlation functions corresponding to the latent function of interest.
For example, to study the response of species $j$ along covariate $x_{r}$ we
evaluate the posterior predictive distribution for $h_{j,r}(x_{r})$ using
(4.2) so that we replace $C_{\ast,\operatorname{\mathbf{f}}}$

$C_{h_{j,r},\operatorname{\mathbf{f}}}=\sum_{l=1}^{J}\left[u_{h_{r},l}(j,1)\\{\tilde{K}_{h_{l,r}}\\}_{j,1},\cdots,u_{h_{{}_{r}},l}(j,J)\\{\tilde{K}_{h_{l,r}}\\}_{j,J}\right],$
(4.15)

and similarly for $C_{\ast}$. However, in case of the traditional LMC (3.9)
only the spatial random effects are correlated between species whereas the
predictor functions are not. Hence, $\Sigma_{h_{r}}$ is diagonal for all $r$
so that $u_{h_{r},l}(j,j^{\prime})=\sigma_{l}^{2}$ if $l=j=j^{\prime}$ and
zero otherwise. The second terms of (4.10) is then block diagonal and only the
$j$’th block column in $\eqref{eq:pred_Chf}$ is non-zero for which reason
there is no information exchange between species specific predictor functions.
In case of univariate GP prior all the terms in (4.10) are block diagonal and
(4.2) reduces to $J$ independent Gaussian distributions with no information
exchange among species at all. Hence, the predictive performance of the
univariate GP and the two multivariate GP models are very different. This is
illustrated in section 5.

####  4.1.3 The (marginal) posterior expectation and variance for new
observations

When predicting species occurrence probability or abundance, we need to
marginalize over the posterior of the latent variables. The posterior
expectation and variance for a new outcome for species $j$ at a location
$(\operatorname{\mathbf{x}}_{*},\operatorname{\mathbf{s}}_{*})$ conditional on
hyperparameters are

$\displaystyle\operatorname{\mathbb{E}}[Y_{j}(\operatorname{\mathbf{x}}_{*},\operatorname{\mathbf{s}}_{*})|\operatorname{\mathbf{y}},\operatorname{\boldsymbol{\eta}},\Lambda]=$
$\displaystyle\operatorname{\mathbb{E}}[\operatorname{\mathbb{E}}(Y_{j}(\operatorname{\mathbf{x}}_{*},\operatorname{\mathbf{s}}_{*})|f_{j}(\operatorname{\mathbf{x}}_{*},\operatorname{\mathbf{s}}_{*}),\operatorname{\mathbf{y}},\operatorname{\boldsymbol{\eta}},\Lambda)]$
(4.16)
$\displaystyle\operatorname{\mathbb{V}}[Y_{j}(\operatorname{\mathbf{x}}_{*},\operatorname{\mathbf{s}}_{*})|\operatorname{\mathbf{y}},\operatorname{\boldsymbol{\eta}},\Lambda]=$
$\displaystyle\operatorname{\mathbb{E}}[\operatorname{\mathbb{V}}(Y_{j}(\operatorname{\mathbf{x}}_{*},\operatorname{\mathbf{s}}_{*})|f_{j}(\operatorname{\mathbf{x}}_{*},\operatorname{\mathbf{s}}_{*}),\operatorname{\mathbf{y}},\operatorname{\boldsymbol{\eta}},\Lambda)]\
+$ $\displaystyle\hskip
85.35826pt\operatorname{\mathbb{V}}[\operatorname{\mathbb{E}}(Y_{j}(\operatorname{\mathbf{x}}_{*},\operatorname{\mathbf{s}}_{*})|f_{j}(\operatorname{\mathbf{x}}_{*},\operatorname{\mathbf{s}}_{*}),\operatorname{\mathbf{y}},\operatorname{\boldsymbol{\eta}},\Lambda)]$

When the probabilistic model for species $j$ is assumed to be the Negative-
Binomial or Binomial model with logistic link function (3.3) we can find
either an exact result or an analytical approximation for these expectations
and variances. These are given in the 7.3. These solutions speed up the
predictive calculations considerably compared to numerically integrating $f$
over $\mathbb{R}$.

####  4.1.4 Conditional scenario predictions

In some applications, one might be interested in scenario predictions
conditional on changes in species composition. For example, one might be
interested in how removing from or introducing specific species into an area
would affect other species. In our case study setting, we could be interested
in, for example, effects of management actions that would clean filamentous
algae from the shoreline. Such scenario predictions would be naturally tackled
with predictive causal inference (Lindley, 2002; Pearl, 2009). In predictive
causal inference the parameters of the model are inferred with the available
data so far and predictions made considering alternative scenarios.

To illustrate this lets first introduce a short notation
$\operatorname{\mathbf{y}}_{\mathcal{J}_{1},*}$ $=$
$\operatorname{\mathbf{y}}_{\mathcal{J}_{1}}(\operatorname{\mathbf{x}}_{*},\operatorname{\mathbf{s}}_{*})$
and $\operatorname{\boldsymbol{f}}_{\mathcal{J}_{1},*}$ $=$
$\operatorname{\boldsymbol{f}}_{\mathcal{J}_{1}}(\operatorname{\mathbf{x}}_{*},\operatorname{\mathbf{s}}_{*})$
where $\mathcal{J}_{1}$ denotes the set of species to be predicted and
$\mathcal{J}_{2}$ the _scenario species_ assumed to be “managed”. For brevity,
we omit the conditioning on hyperparameters and data. The conditional
distribution of
$\operatorname{\boldsymbol{Y}}_{\mathcal{J}_{1}}(\operatorname{\mathbf{x}}_{*},\operatorname{\mathbf{s}}_{*})|\operatorname{\mathbf{y}}_{\mathcal{J}_{2}}$
is

$\displaystyle\pi(\operatorname{\mathbf{y}}_{\mathcal{J}_{1},*}|\operatorname{\mathbf{y}}_{\mathcal{J}_{2}})$
$\displaystyle=\displaystyle\int\pi(\operatorname{\mathbf{y}}_{\mathcal{J}_{1},*}|\operatorname{\mathbf{y}}_{\mathcal{J}_{2}},\operatorname{\mathbf{f}}_{\mathcal{J}_{1},*})\pi(\operatorname{\mathbf{y}}_{\mathcal{J}_{2}},\operatorname{\mathbf{f}}_{\mathcal{J}_{1},*})\mathrm{d}\operatorname{\mathbf{f}}_{\mathcal{J}_{1},*}/\pi(\operatorname{\mathbf{y}}_{\mathcal{J}_{2}})$
(4.17)

where
$Y_{\mathcal{J}_{1}}(\operatorname{\mathbf{x}}_{*},\operatorname{\mathbf{s}}_{*})|\operatorname{\mathbf{y}}_{\mathcal{J}_{2}},\operatorname{\mathbf{f}}_{\mathcal{J}_{1},*}$
only depends on $\operatorname{\mathbf{f}}_{\mathcal{J}_{1},*}$. The Laplace
approximation for this conditional predictive distribution is shown in 7.3.

###  4.2 Parameter inference

We used Laplace approximation (Vanhatalo et al., 2010) to approximate the
conditional posterior of latent variables
$\pi(\operatorname{\mathbf{f}}|\operatorname{\mathbf{y}},\operatorname{\boldsymbol{\eta}},\Lambda)$
and the marginal likelihood of the hyperparameters
$\pi(\operatorname{\mathbf{y}}|\operatorname{\boldsymbol{\eta}},\Lambda)=\int\pi(\operatorname{\mathbf{y}}|\operatorname{\mathbf{f}},\operatorname{\boldsymbol{\eta}})\pi(\operatorname{\mathbf{f}}|\Lambda)d\operatorname{\mathbf{f}}$.
We searched for the (approximate) maximum a posterior (MAP) estimate for
hyperparameters given by

$(\hat{\operatorname{\boldsymbol{\eta}}},\hat{\Lambda})=\underset{\operatorname{\boldsymbol{\eta}},\Lambda}{\operatorname*{arg\,max}}\log
q(\operatorname{\mathbf{y}}|\operatorname{\boldsymbol{\eta}},\Lambda)+\log\pi(\operatorname{\boldsymbol{\eta}},\Lambda).$
(4.18)

where $\log
q(\operatorname{\mathbf{y}}|\operatorname{\boldsymbol{\eta}},\Lambda)$ is the
Laplace approximation for the log marginal likelihood for parameters. This
approach produces also good approximation for the posterior predictive
distribution for latent variables (Vanhatalo et al., 2010; Vehtari et al.,
2016), which are the main interest in this work. Hence, we fixed
hyperparameters to their MAP estimate.

In order to avoid constrained optimization, all the parameters were
transferred to unconstrained space. We used log transformation for covariance
function parameters, and for the interspecific correlation matrices we used
the transformation presented by Kurowicka and Cooke (2003) and Lewandowski et
al. (2009), which is a bijective mapping between the space of correlation
matrices and $\mathbb{R}^{J\choose 2}$. This is summarized in 7.3. We used
scaled conjugate gradient optimization for locating the MAP and checked
carefully that a (local) mode had really been found by verifying that
gradients along all hyperparameters were zero. The required gradients of $\log
q(\operatorname{\mathbf{y}}|\operatorname{\boldsymbol{\eta}},\Lambda)$ were
solved analytically as described by Rasmussen and Williams (2006) and
Vanhatalo et al. (2010). The 7.3 summarizes the derivatives w.r.t. to the
covariance parameters in $\Sigma_{\epsilon}$, $\Sigma_{h_{r}}$, $r$ $=$
$1,\ldots,c$.

##  5 Experiments

Here, we first illustrate the properties of the hierarchical multivariate GP
models with two simple examples. These examples highlight particular
properties of the proposed model. After this we introduce the model and
analysis for our case study (Section 2). We implemented all the models using
the GPstuff package (Vanhatalo et al., 2013).

###  5.1 Demonstration with simulated spatial data

Figure 2 presents simulated data and posterior predictions for spatial
modelling of two species. The model follows spatial LMC; that is model (3.9),
where the covariate terms are dropped out. The spatial correlation function
used in this demo is the Gaussian
$k(\operatorname{\mathbf{s}},\operatorname{\mathbf{s}}^{\prime})=e^{-||\operatorname{\mathbf{s}}-\operatorname{\mathbf{s}}^{\prime}||^{2}/2l^{2}}$,
where $||\operatorname{\mathbf{s}}-\operatorname{\mathbf{s}}^{\prime}||$ is
the Euclidean distance and $l$ is the length-scale. The first row of subplots
shows the posterior predictive probability of observing species 1 ($E[Y_{1}]$,
$Y_{1}(\operatorname{\mathbf{s}})\sim\mathrm{Binomial}$) and the second row
shows the expected number of species 2 ($E[Y_{2}]$,
$Y_{2}(\operatorname{\mathbf{s}})\sim$ Negative-Binomial). Plots (a) and (e)
show the predictions when the species observations are from the same region
but not from exactly the same locations. In this case, the prediction of
multivariate GP is similar to predictions by univariate GPs. Plots (b) and (f)
show the predictions when data are available on both species from the lower
left corner and additionally data on species 1 is available from upper right
corner. There is positive correlation between species, which has been inferred
from the data in the lower left region, so the prediction for species 2 in
upper right corner is informed by data on species 1 in that region. The last
two columns illustrate the conditional scenario predictions (4.17). In the
third column, the training data is the same as in the first column and the
expected values in plots (c) and (g) show joint scenario prediction in new
region of same size and shape. In this scenario, species in the new region
were observed so that species 1 is known to be in the top right, and species 2
in the bottom left. The fourth column shows the conditional scenario
prediction for both of the species separately so that the grey marks show the
locations where the other species would be introduced in these scenarios.

Figure 2: Illustration of spatial multivariate GP prior for JSDM with two
species and spatial only component. Crosses and dots are simulated data
locations of species 1 (Binomial data) and species 2 (Negative-Binomial data).
See Section 5.1 for discussion.

###  5.2 Demonstration with time series of species abundances

Next, we consider the classical predator-prey system containing annual
population counts of hare and lynx in the northern Canada from 1845 to 1935
(Odum, 1953). The data is not spatially explicit since the observations are
total counts over the region so we model it as time series. This illustrates
the behavior of individual additive covariate effect in the multivariate GP
model (3.10) with time being the covariate. Let’s denote by $Y_{1,t}$ and
$Y_{2,t}$ the population sizes of hare and lynx at time $t$, and assume that
$Y_{1,t}|f_{1}(t)\sim\mathrm{Poisson}(\exp(f_{1}(t))$ and
$Y_{2,t}|f_{2}(t)\sim\mathrm{Poisson}(\exp(f_{2}(t))$. We compared two
alternative priors for the latent functions, one with independent GP priors
and another with joint bivariate GP prior with the LMC structure. In order to
compare models’ predictive performances when some species have not been
observed, we removed parts of the original data from the training set, the
period of 1870-1900 for hare and 1850-1870 for lynx. Figure 3 displays the
result of the data analysis. The model with bivariate GP prior clearly
outperforms the independent model since it predicts better the test data
points in periods where training data was removed. Moreover, its predictions
have also smaller uncertainty than the predictions of independent GPs in those
periods of time.

Figure 3: The results from the analysis of hare and lynx interaction
(predator-prey system) analyzed with single and joint SDM. Plots a) and c)
show the posterior mean and 95% credible interval of the latent functions and
plots b) and d) show the training data, test data and posterior expectation of
observations. See Section 5.2 for explanation.

###  5.3 The case study on coastal species distribution

####  5.3.1 Case study models

The plant data are modeled with the Binomial observation model (3.3) and the
fish data are modeled with the Negative-Binomial observation model (3.4). In
order to test the effect of different model components and the effect of semi-
parametric response functions versus standard quadratic response functions we
compare the following models:

  * 1)

additive univariate GP predictor functions and univ. spatial random effects
(3.5),

  * 2)

additive univariate GP predictor functions and LMC spatial random effect
(3.9),

  * 3)

additive multivariate GP predictor function and LMC spatial random effect
(3.10),

  * 4)

additive univariate GP predictor functions only (model (3.5) without
$\epsilon_{j}(\operatorname{\mathbf{s}})$),

  * 5)

univariate spatial random effects only (model (3.5) without
$h_{j}(\operatorname{\mathbf{x}}_{j})$),

  * 6)

additive univariate quadratic predictors and univariate spatial random effects
(model (3.5) with
$h_{j}(\operatorname{\mathbf{x}}_{j})=\operatorname{\mathbf{x}}_{j}\beta_{j}$
and $\operatorname{\mathbf{x}}_{j}$ includes the covariates and their
squares),

  * 7)

additive univariate quadratic predictors and LMC spatial random effect (model
(3.9) with
$h_{j}(\operatorname{\mathbf{x}}_{j})=\operatorname{\mathbf{x}}_{j}\beta_{j}$
and $\operatorname{\mathbf{x}}_{j}$ includes the covariates and their
squares),

  * 8)

additive multivariate quadratic predictors and LMC spatial random effects
(model (3.10) with
$h_{j}(\operatorname{\mathbf{x}}_{j})=\operatorname{\mathbf{x}}_{j}\beta_{j}$
and $\operatorname{\mathbf{x}}_{j}$ includes the covariates and their
squares).

  * 9)

additive multivariate quadratic predictors only (model (3.10) without
$\epsilon(\operatorname{\mathbf{s}})$ with
$h_{j}(\operatorname{\mathbf{x}}_{j})=\operatorname{\mathbf{x}}_{j}\beta_{j}$
and $\operatorname{\mathbf{x}}_{j}$ includes the covariates and their
squares).

In each model, the spatial random effects were given the Mátern covariance
function with $\nu=3/2$ degrees of freedom
$k_{\epsilon_{j}}(\operatorname{\mathbf{s}},\operatorname{\mathbf{s}}^{\prime};\sigma^{2}_{j},\operatorname{\boldsymbol{\ell}}_{j})=\sigma^{2}_{j}\big{(}1+\sqrt{3}\
r_{j}(\operatorname{\mathbf{s}},\operatorname{\mathbf{s}}^{\prime})\big{)}e^{-\sqrt{3}r_{j}(\operatorname{\mathbf{s}},\operatorname{\mathbf{s}}^{\prime})}$
where $\sigma_{j}^{2}$ is the variance parameter,
$\operatorname{\boldsymbol{\ell}}_{j}^{T}=[\ell_{j,1}\ \ell_{j,2}]$ is the
vector of length-scales and
$r_{j}(\operatorname{\mathbf{s}},\operatorname{\mathbf{s}}^{\prime})=[(\operatorname{\mathbf{s}}-\operatorname{\mathbf{s}}^{\prime})^{T}(\operatorname{\mathrm{diag}}(\operatorname{\boldsymbol{\ell}}_{j})^{2})^{-1}(\operatorname{\mathbf{s}}-\operatorname{\mathbf{s}}^{\prime})]^{1/2}$.
The continuous covariate effects in the additive GP models (models 1-4) were
given the Gaussian correlation function
$\tilde{k}_{h_{j},r}(x_{r},x_{r}^{\prime})=e^{-||x_{r}-x_{r}^{\prime}||^{2}/2l_{j,r}^{2}}$.
The additive linear models (models 6-8) were coded as additive GPs with
$\tilde{k}_{h_{j,r}}(x_{r},x^{\prime}_{r})=x_{r}x^{\prime}_{r}\sigma^{2}_{j,r}$.
Even though this is not a proper correlation function, when used in LMC it
implies a generalized linear model with interspecific dependencies between
weights, $\beta_{j}$, that are the key components in state-of-the-art
parametric JSDMs. See Discussion for details. For the categorical covariate
(Bottom class, table 1) we used a delta function so that
$\tilde{k}_{h_{j,r}}(x_{r},x_{r}^{\prime};\operatorname{\boldsymbol{\theta}})=\sigma_{\delta_{j}}^{2}\delta_{x_{r}}(x_{r}^{\prime})$,
where $\delta_{x_{r}}(x_{r}^{\prime})=1$ if $x_{r}=x_{r}^{\prime}$ and zero
otherwise. This corresponds to having an own intercept for each category. We
used the marginally uniform priors (section 4.1.3) for the between species
correlations and weakly informative priors for the rest of the parameters. The
variance parameters were given
$\mathrm{Half}\mbox{-}\mathrm{Student}\mbox{-}\it{t}$($\mu=0$, $\sigma^{2}=4$,
$\nu=4$) priors and the length-scale parameters were given
$\mathrm{Half}\mbox{-}\mathrm{Inverse\mbox{-}Student}\mbox{-}\it{t}$($\mu=0$,
$\sigma^{2}=1$, $\nu=4$). Hence, _a priori_ more weight is given for smooth
functions with small variability.

####  5.3.2 Model validation and comparison

We aim to provide models that give reliable posterior predictions. Hence, it
is natural to compare models with the goodness of their posterior predictive
distributions. This can be done efficiently with cross-validation (CV) using
the log predictive density diagnostics (Vehtari and Ojanen, 2012). Let
$\mathcal{D}$ denote the full data-set. Fix $K$ disjoint sub-sets of
$\mathcal{D}$, say $\mathcal{D}_{1},\ldots,\mathcal{D}_{K}$, such that their
union is $\mathcal{D}$. The K-fold CV log point-wise marginal predictive
density is then
$L_{K}(\mathcal{D})=\frac{1}{n}\sum_{i=1}^{n}\log\pi(y_{i}|\mathcal{D}_{\sim
k(i)})$ where $\mathcal{D}_{\sim k(i)}=\\{\cup_{r=1}^{K}\mathcal{D}_{r}:r\neq
k(i)\\}$ and $k(i)$ is such that $y_{i}\in\mathcal{D}_{k(i)}$. We compare the
models with the leave-one-out CV (LOO-CV, $K=n$) and structured 5-fold CV. We
used the MAP estimate for the hyperparameters. In 5-fold CV we found the MAP
for each training set separately. The LOO-CV was conducted at the MAP found
with full data and we used Laplace approximation to approximate the LOO-CV log
predictive densities (Vehtari et al., 2016). Laplace approximation for LOO-CV
is shown to work well in GP models and, since our data is rather large,
leaving only one data point out of the training set has only negligible effect
on the posterior of the hyperparameters (Vehtari et al., 2016).

The rationale for calculating both LOO-CV and 5-fold CV comparison is the
following. Since multiple species were sampled at every sampling site and each
of the sampling sites has other sites spatially nearby it, the LOO-CV log
predictive densities are affected significantly by the spatial random effects.
The LOO-CV, thus, measures models’ interpolation performance which can be good
even if the models were not able to represent well responses along covariates
(predictor functions) (Vanhatalo et al., 2012; Veneranta et al., 2013). For
this reason, we structured the 5-fold CV by dividing the data into five
spatially distinct groups corresponding to regions I-V in Figure 1. The
sampling sites in different groups are spatially so far from each others that
the spatial random effects do not affect the posterior predictive
distributions for test groups. Hence, the 5-fold CV tests mostly the
extrapolation performance of a model, which is governed by the goodness of the
predictor functions, whereas the LOO-CV tests the interpolation performance,
which is governed also by the spatial random effects.

GP models | LOO Cross Validation | 5-fold Cross validation
---|---|---
1) Univariate GP $h(\operatorname{\mathbf{x}})$, univ. $\epsilon(\operatorname{\mathbf{s}})$ | -2.230 (0.082) | -2.465 (0.094)
2) Univariate GP $h(\operatorname{\mathbf{x}})$, multiv. $\epsilon(\operatorname{\mathbf{s}})$ | -2.081 (0.080) | -2.480 (0.100)
3) Multiv. GP $h(\operatorname{\mathbf{x}})$, multiv. $\epsilon(\operatorname{\mathbf{s}})$ | -1.669 (0.080) |  -2.316 (0.086)
4) Multiv. GP $h(\operatorname{\mathbf{x}})$ only | -2.228(0.089) | -2.590(0.012)
5) Multiv. $\epsilon(\operatorname{\mathbf{s}})$ only | -2.351(0.087) | -2.500(0.086)
6) Univ. quadr. $h(\operatorname{\mathbf{x}})$, univ. $\epsilon(\operatorname{\mathbf{s}})$ | -2.262 (0.083) | -2.517 (0.095)
7) Univ. quadr. $h(\operatorname{\mathbf{x}})$, multiv. $\epsilon(\operatorname{\mathbf{s}})$ | -2.117 (0.084) | -2.484 (0.109)
8) Multiv. quadr. $h(\operatorname{\mathbf{x}})$, multiv. $\epsilon(\operatorname{\mathbf{s}})$ | -2.105 (0.084) | -2.416 (0.099)
9) Multiv. quadr. $h(\operatorname{\mathbf{x}})$ only | -10.246(1.297) | -9.305(1.008)
Table 2:  Model comparison with leave-one-out (LOO) and 5-fold cross
validation using the mean point wise log marginal predictive density
statistics (and its standard error) over all species. see Section 5.3.1 for
model descriptions.

##  6 Results

###  6.1 Predictive performance of models

Table 2 summarizes the CV comparisons over all species and tables 1-2 in 7.3
show the models’ predictive performance for each species separately. Tables
3-4 in 7.3 show the pairwise comparison of the CV point wise log predictive
densities between the best GP/parametric model (models 3 and 8) and the other
GP/parametric models with both the covariate effects and spatial random effect
(models 1-2 and 6-7). The results show that model 3, which includes
multivariate GP predictors and multivariate spatial random effects (3.10), is
the best in both LOO-CV and 5-fold CV with a difference of $10^{-2}$ or more
to the other models. Since the mean log predictive density is an average of
$n=850$ point-wise predictions the difference of $10^{-2}$ corresponds to a
difference of 8.5 in log (point-wise) joint predictive densities. Analogously
to Bayes factors, which compare log prior predictive distribution, this can be
considered a significant difference between two models (Kass and Raftery,
1995).

Hence, the multivariate GP (model 3) has significantly better overall
predictive performance over all species than the other models. According to
posterior predictive checks (Gelman et al., 2013) there are no serious
discrepancies between its predictions and observed data. In extrapolation
(5-fold CV) model 3 performs best also for all species separately. However, in
interpolation (LOO-CV) it is not the best for all species separately (tables
1-2 in 7.3). Moreover, the semi-parametric GP models (models 1-4) work better
than the corresponding parametric models (quadratic environmental response,
models 6-9) in both interpolation and extrapolation in general. Dropping
either spatial random effect or covariate effect out from the model decreases
its performance clearly. All models work better in interpolation than in
extrapolation and compared to univariate models including interspecific
correlations either in spatial random effects or also in environmental
predictors improves models’ performance in both of these tasks. Moreover,
including interspecific correlations between environmental responses,
$h_{j,r}(\operatorname{\mathbf{x}}_{r})$, improves the extrapolation
performance relatively more than including interspecific correlations only for
spatial random effects. Hence, multivariate spatial random effect improved
more interpolation whereas multivariate predictors has larger effect for
extrapolation as would intuitively be expected.

Figure 4: The maximum a posterior estimate of the correlation matrices in the
multivariate additive GP model with Gaussian covariance function (model 3).

###  6.2 Effects of environmental covariates

Figure 5: Marginal latent responses along covariates, $h_{j,r}$. Grey
corresponds to single species GP model (model 1) and black to multivariate
additive GP model (model 3).

The interspecific correlations (Figure 4) show that the responses to
environment are similar among Coregonids and among sticklebacks whereas there
are clear differences among these groups. These fish groups respond
differently to sandy bottom and distance to deep. Moreover, there is strong
positive spatial correlation among Coregonids and among sticklebacks but not
between these groups. All species have negative spatial correlation with
filamentous algae whereas there is either weak positive or negligible spatial
correlation between fish species and macrovegetation and diatomous algae.

Figure 5 shows the marginal effects of the predictor functions for both
univariate (SSDM) and multivariate (JSDM) GP models calculated as described in
section 4.1.2. In general, the responses in JSDM and SSDM are similar. In most
of the cases the uncertainty in response function is smaller (narrower
credible regions) in JSDM than in SSDM although the differences are not large.
Most of the responses are smooth but show patterns that would be hard to
capture with a quadratic function. For example, three and nine spine
stickleback, macrovegetation and filamentous algae show logistic style
response to either distance to deep or openness. Moreover, the response of
vendace on ice break up week is first constant but has very steep increase
after week 16. The responses to ice break up or chlorophyll-a show non-linear
and non-quadratic responses also for white fish, macrovegetation and
filamentous algae. The MAP estimates of other hyper-parameters are summarized
in 7.3.

###  6.3 Spatial predictions

Figure 6 shows the distribution maps as posterior median of intensity and
expected count of individuals in replicate sampling produced by SSDM (model 1)
and JSDM (model 3) for vendace. The distribution maps for other species
together with maps on posterior predictive variance are shown in 7.3. In broad
scale the posterior median of the intensity looks similar with SSDM and JSDM
whereas SSDM predicts larger species counts than JSDM for all species
throughout the study region and the uncertainty related to SSDM predictions is
larger than that of JSDM predictions. Both SSDM and JSDM predict that vendace
is distributed mostly in the northern Gulf of Bothnia. However, SSDM predicts
somewhat higher median intensity than JSDM also in the southern Gulf of
Bothnia. In relation to median intensity similar pattern that JSDM predicts
more restricted distribution range than the single species model is seen also
in the predictions concerning sticklebacks and whitefish. In case of diatomous
and filamentous algae SSDM and JSDM predict the distribution pattern similarly
whereas for macrovegetation JSDM predicts slightly larger distribution ranges.
The posterior distributions of spatial lenght-scale parameters were
concentrated near one kilometer or less (see table 5 in 7.3). Hence, the
differences in distribution predictions cannot be explained by the spatial
random effects over large areas but the spatial random effect explains local
deviations from the predictions based on covariates.

Figure 6: Posterior predictive median of intensity and the expected count of
individuals in replicate sampling for vendace predicted by SSDM (model 1) and
JSDM (model 3).

##  7 Discussion and concluding remarks

###  7.1 Case study results

The GP based SDMs had better predictive performance than the parametric SDMs
and JSDMs had better predictive performance than the corresponding SSDMs. The
differences between SSDM and JSDM are most apparent in the distribution maps
and predictor functions. The JSDM predicted in general smaller distribution
ranges than SSDMs (macrovegetation being the only exception) and the posterior
uncertainty in their predictions were smaller than in SSDMs.

When interpreting the results, it should be remembered that the sampling was
targeted to larvae of sea-spawning Coregonids (whitefish and vendace) and
planned to cover their plausible distribution range in terms of spatial
(coastal areas of Gulf of Bothnia) and temporal extent (spring). Other species
were sampled alongside Coregonids and the sampling area covers only limited
portion of their full distribution range in the spring. Moreover, the
abundance and distribution of all the studied species vary annually. Fish
change their distribution areas seasonally and their larval stages last only
few weeks. Vegetation and algal cover vary due changes in temperature and ice
effects. For these reasons the results are most representative for larvae of
whitefish and vendace. For other species the results describe their spring
distribution in the shallow coastal regions only.

Our results correspond rather well to the earlier knowledge on the studied
species. The responses to the environmental covariates and the interspecific
correlations (Figures 5 and 4) indicate that whitefish and vendace larvae are,
in general, distributed in different areas than sticklebacks. Most of the
literature on Baltic Sea sticklebacks focus on three-spined sticklebacks
whereas the nine-spined stickleback is not well studied. In early spring,
stickleback abundances have been found to be highest in sheltered archipelago
areas, where part of the population overwinter. High abundance of sticklebacks
are typically thought to indicate structural complexity on the bottom; such as
the presence of stones and boulders as well as reeds and other macrovegetation
that function as shelter and provide food (Peltonen et al., 2004; Sieben et
al., 2011). On contrary, highest densities of whitefish and vendace larvae are
observed in open sandy shores near deep areas and in structurally simple
shores, without macrovegetation, boulders and stones (Hudd et al., 1988;
Leskelä et al., 1991; Veneranta et al., 2013). The distribution of vendace is
highly influenced by ice break up week and salinity so that long ice cover
period and low salinity increase their abundance. Ice winter is longer and
salinity is lower in the northern than in the southern Gulf of Bothnia and
thus it correlates positively to Coregonid presence (Vanhatalo et al., 2012).
In general, optimal habitats for Coregonid larvae are mainly located in the
northernmost Gulf of Bothnia.

Sticklebacks are known to prey mainly on mesozooplankton, but also on grazers
(Peltonen et al., 2004; Sieben et al., 2011). Sticklebacks feed also on fish
larvae if available (Byström et al., 2015), but there are no studies on
potential predation risk to Coregonids. Based on our results, in the scale of
coastal region of Gulf of Bothnia, sticklebacks are not thread to Coregonid
larvae since their high density areas do not overlap. Moreover, there was no
spatial correlation between sticklebacks and coregonids (Figure 4) which could
indicate interspecific interaction of any kind. The presence of sticklebacks
has been connected to higher eutrophication status and stickleback
reproduction and growth benefit from increasing temperature and eutrophication
(Lefébure et al., 2011; Candolin et al., 2008; Meier et al., 2012) On
contrary, these environmental characteristics are assumed to affect negatively
the Coregonid reproduction (M.Cingi et al., 2010; Müller, 1992; Veneranta et
al., 2013; Vanhatalo et al., 2012). In our results whitefish and vendace
larvae have positive response to Chlorophyl-a, which is a strong indicator for
increasing eutrophication. However, in the scale of Gulf of Bothnia
Chlorophyl-a concentration is typically higher near river mouths where
salinity is lower and nutrient inflow high. Hence the result more likely
reflects the high river influence through low salinity than preference for
eutrophicated water.

The vegetation and algae distribution in our results reflect the nutrient
status during winter and early spring as well as effect of ice cover and ice
scraping in wind exposed shallow areas. Both the SSDMs and JSDMs indicate that
filamentous algae occur in all coastal areas in high densities, except in some
sheltered inner coastal areas. Macrovegetation and diatomous algae had highest
presence probabilities at areas with lower filament presence probability. This
pattern agrees well with the general understanding of these species groups.
Filamentous algae are typically distributed in wind and wave exposed shores,
that epiphytic diatoms and reeds that require shelter cannot tolerate.
Filamentous algae can, however, grow over macrovegetation and it is possible
that macrovegetation distribution is underestimated if filamentous algae
growth over macrovegetation has hided macrovegetation from the sampling
pictures. The higher nutrient levels in sheltered areas have been found to
affect positively reed belt growth, and high abundances of macrovegetation
species have been found from archipelago areas, lagoons, bays and river inlets
(Pitkänen et al., 2013; Altartouri et al., 2014). The JSDM model reflects
these smaller scale occurrences better than SSDM. In our results, diatomous
algae are more common in estuaries and northern parts of the study area (see
7.3). This is likely explained by longer ice winter towards northern Gulf of
Bothnia and longer distance to deep water.

###  7.2 Multivariate additive Gaussian process

Next we discuss some of the similarities and differences between our model and
the existing JSDMs and highlight the novel methodological contributions in
this paper.

####  7.2.1 Interspecific correlations

From the predictive point of view, the interspecific correlations in
predictive functions are attractive for many reasons. In many applications of
SDMs the aim is to predict species distribution over regions that include
locations spatially far from the data (Record et al., 2013) or to conduct
scenario predictions related to, for example, climate change or land use
(Guisan et al., 2013). In these applications, predictions are based on the
responses to environmental covariates. The inclusion of interspecific
dependence allows information flow between species which improves the
estimates for predictive functions especially for species with only scarce
data (Thorson et al., 2015; Ovaskainen et al., 2017; Hui et al., 2013; Clark
et al., 2014).

We used the marginally uniform priors of (Barnard et al., 2000; Tokuda et al.,
2012) for the interspecific correlations. This is justified by prior ignorance
on correlations. Prior information on interspecific correlations could,
however, be added into model with, e.g., informative inverse Wishart prior.
With many species inferring the full covariance matrix is hard in which case
we could use spatially dependent latent factors which induce for the spatial
random effects a covariance structure
$\mathrm{Cov}\big{(}\epsilon_{j}(\operatorname{\mathbf{s}}),\epsilon_{j}^{\prime}(\operatorname{\mathbf{s}}^{\prime})\big{)}=\sum_{q=1}^{M}\lambda_{qj}\lambda_{qj^{\prime}}k_{q}(\operatorname{\mathbf{s}},\operatorname{\mathbf{s}}^{\prime})$,
where $k_{q}$ is the $q$’th spatial covariance function and $\lambda_{qj}$ the
corresponding species specific factor loading. Here, $M<J$ so that this
covariance structure is a low rank representation of LMC (Lopez, 2000; Lopez
et al., 2008). Latent factor representation was first introduced to SDMs by
Thorson et al. (2015) and it is used also in HMSC Ovaskainen et al. (2017).
Recently Taylor-Rodríguez et al. (2017) proposed a clustering scheme where the
species are clustered to less than $J$ factor loading vectors.

Interspecific correlation between response functions has received less
interest. Typically the response functions are assumed to be independent among
species. Exception are species archetype models (Dunstan et al., 2013) and the
HMSC framework of Ovaskainen et al. (2017). In the latter, the response
functions are defined as
$h_{j}(\operatorname{\mathbf{x}}_{j})=\operatorname{\mathbf{x}}_{j}\theta_{j\cdot}$
where $\theta=[\theta_{1\cdot}^{T},\dots,\theta_{J\cdot}^{T}]^{T}$ is a
$J\times p$ matrix of regression weights with hierarchical prior. The species
specific weights are given Gaussian prior $\theta_{j\cdot}\sim
N(\mu_{j\cdot},V_{j})$ where $\mu_{j\cdot}$ is the expected response of a
species that can be common for all species, common within groups of species or
modelled through species traits $\mu_{jr}=\tau^{T}_{j}\gamma_{r}$ where
$\tau_{j}$ is a vector of trait values of species $j$ and $\gamma_{r}\sim
N(0,V_{\gamma})$ are the effects of traits to response along covariate $r$.
With $V=\sigma^{2}I_{J\times J}$ and common prior mean $\mu_{j\cdot}=\mu\sim
N(0,\sigma_{\mu}^{2})$ for all $j\in 1,\dots,J$, the induced covariance
between species specific additive response functions is

$\displaystyle\operatorname{\mathrm{Cov}}\left(h_{j,r}(x_{r}),h_{j^{\prime},r}(x^{\prime}_{r})\right)=$
$\displaystyle\mathbb{E}[\operatorname{\mathrm{Cov}}\left(x_{r}\theta_{jr},x^{\prime}_{r}\theta_{j^{\prime}r}\right)]+\operatorname{\mathrm{Cov}}\left(\mathbb{E}[x_{r}\theta_{jr}],\mathbb{E}[x^{\prime}_{r}\theta_{j^{\prime}r}]\right)$
$\displaystyle=$
$\displaystyle(\sigma_{\mu}^{2}+\delta_{j}(j^{\prime})\sigma^{2})x_{r}x^{\prime}_{r}$

and with trait dependent prior mean
$\operatorname{\mathrm{Cov}}\left(h_{j,r}(x_{r}),h_{j^{\prime},r}(x^{\prime}_{r})\right)=(\tau_{j}^{T}V_{\gamma}\tau_{j^{\prime}}+\delta_{j}(j^{\prime})\sigma^{2})x_{r}x^{\prime}_{r}$.
Hence, these alternatives have the same covariance structure as an additive
multivariate GP (3.10) with $k_{h_{j},r}=x_{r}x^{\prime}_{r}$ for all
$j=1,\dots,J$, and interspecific covariances
$\Sigma_{h_{r}}=\sigma_{\mu}^{2}+\delta_{j}(j^{\prime})\sigma^{2}$ and
$\Sigma_{h_{r}}=\tau_{j}^{T}V_{\gamma}\tau_{j^{\prime}}+\delta_{j}(j^{\prime})\sigma^{2}$.
Hence, the hierarchical prior structure of HMSC could easily be extended to
semiparametric models as well if we had trait information from our species.
The benefit would be that these hierarchical model structures contain
ecologically relevant prior information and, hence, using the induced
interspecific covariance structure in the additive multivariate GP model would
make it ecologically more realistic and potentially improve its predictive
performance. The HMSC framework has many other features as well, such as
phylogenetic relationships Ovaskainen et al. (2017), which could in principle
be incorporated with our approach as well. The structure of interspecific
correlations is likely to be especially important in scenario based predictive
analyses, such as climate change predictions.

Since our model does not make any mechanistic assumptions on species
interdependencies the results concerning interspecific correlations have to be
interpreted with care. For example, a positive correlation between spatial
random effects could arize due mutualism, predator prey association or similar
response to unobserved environmental covariates. Hence, interpreting these
correlations should be done in light of more general ecological knowledge on
the species. This is, however, a common challenge with all current JSDMs
(Thorson et al., 2015; Dunstan et al., 2013; Ovaskainen et al., 2017; Taylor-
Rodríguez et al., 2017). Moreover, the property of our model, as well as most
other JSDM, is that if the training data shows positive or negative spatial
correlation between two species this positive correlation is assumed to hold
everywhere in the study domain. In many cases this might be an unrealistic
assumption for which reason one interesting future development would be to
extend our models to allow species-to-species associations to depend on the
environmental context(Fox and Dunson, 2015; Tikhonov et al., 2017). However,
the colocalization patterns are explained by both environmental covariates and
spatial random effects and regional differences in colocalization patterns can
be explained by different environmental covariates in different regions.

The interspecific correlation in our approach are modelled on the latent
variable level and, hence, are not directly comparable to interspecific
correlations in data. Some authors have argued that this is a conceptual
weakness of latent variable models since it makes interpretating their results
hard and thus less transparent. An alternative would be to model the
correlations in the observation error model (see, e.g., Ovaskainen et al.,
2010; Pollock et al., 2014). Clark et al. (2016) proposed generalized joint
attribute model (GJAM) to fix this interpretation challenge by aiming to model
the correlations between species on the data scale. However, for example, in
case of presence absence observations GJAM corresponds to multivariate probit
model that is the marginal likelihood of our model in case of Bernoulli
observation model. It is also unclear how to disentangle the process
underlying the species occurrence and abundance from the observation process
leading to the actual data in GJAM approach. Hence, despite the challenge of
interpretating the correlations, we prefer the hierarchical latent variable
modeling since it provides a coherent approach to separate these two
processes.

####  7.2.2 Response functions

The response functions along environmental covariates provide basis to
understand how environment affects species distribution and to predict the
species distribution in new areas. Linear and other parametric models are
efficient in finding overall trends in responses but have trouble in locating
abrupt changes and change points due, for example, physical tolerance limits.
Such limits, for example, along temperature and salinity are typical for large
variety of taxa in aquatic domains (MacKenzie et al., 2007; Kotta et al.,
2019). With SSDMs Vanhatalo et al. (2012), Shelton et al. (2014), Golding and
Purse (2016) and Kallasvuo et al. (2017) have demonstrated the benefits of
semiparametric models in such situations. The GP approach for modeling
environmental responses is similar to generalized additive models (Guisan et
al., 2002) but the latter have not been implemented as JSDMs. The flexibility
of semiparametric models provides also challenges compared to the more
restricted parametric models. The prior for covariance function parameters has
to be set with care (Golding and Purse, 2016) and we may also need to pose
monotonicity constraints to the response functions Kotta et al. (2019) in
order to prevent overfitting. Inferring the responses reliably requires
typically more data compared to parametric models. The multivariate additive
GP helps in these challenges as well since the interspecific correlations
increase the effectively amount of data to be used to infer each response
function and hence decreases the uncertainty in them.

Additive multivariate GP framework opens also new questions related to the
choice of covariance functions and interspecific correlations. A typical
choice for general purpose covariance function is a radial basis function.
However, predictions using these covariance functions revert to prior
predictive distributions when predicting beyond covariate range covered by
data. Hence, in extrapolation tasks stationary covariance functions may not be
the optimal choice (Vanhatalo et al., 2012; Kallasvuo et al., 2017) and
combining the multivariate GP models with functional constraints (e.g. Kotta
et al., 2019) could improve their predictive performance further.

####  7.2.3 Spatial random effects

Our model comparison results show clearly that including spatial random
effects into model improves their predictive performance in both inter- and
extrapolation. This result is well in line with earlier works demonstrating
that environmental conditions alone may not sufficiently explain species
distribution and spatial random effects improve their predictive performance
(Latimer et al., 2006; Vanhatalo et al., 2012; Clark et al., 2014; Thorson et
al., 2015; Kallasvuo et al., 2017). This is reasonable, since the distribution
of species is shaped by the interplay of environmental covariates, stochastic
processes and species interaction (Ovaskainen et al., 2017). Hence, a
justified model should account for random processes and as demonstrated also
by our results also to species interactions in these processes. However,
adding spatial random effects into model may lead to problems in model
identifibiality which we discuss in the next section.

####  7.2.4 Posterior inference and predictive performance

The JSDM built with multivariate additive GP had clearly the best overall
predictive performance. Moreover, the uncertainty in the predictions by JSDMs
were smaller than in the SSDMs which, together with the best log predictive
density performance, indicates that the predictions were also more accurate.
These findings have important implications to practical use of SDMs for
management and other purposes. For example Kallasvuo et al. (2017) used SSDMs
to classify coastal areas of Finland to unsuitable, suitable and highly
productive spawning regions for four commercially important fish species. They
based their estimates on the expected numbers of larvae in the coastal region
which, as shown by Figure 6 is highly sensitive to the uncertainty of the
predictions. Moreover, by using SSDMs we can overestimate the total biomass of
all species (Clark et al., 2014).

An inherent challenge with the type of models considered here is the model
identifibiality. Spatial random effects are known to affect the fixed effects
estimates (Hodges and Reich, 2010) and in some cases the spatial random effect
can actually capture variability otherwise explainable with fixed effects
(Paciorek, 2010; Bose et al., 2018). Moreover, with uniform (uninformative)
priors the length-scale and variance hyperparameters of Matérn class of
covariance functions are non-identifiable (Zhang, 2004) in general. In order
to alleviate these identifibiality challenges we used weakly informative
priors for the variance and length-scale parameters defined through the
principles proposed by Gelman (2006), Simpson et al. (2017) and Fuglstad et
al. (2018). Our priors for the covariance function parameters favor small
variance and large length-scale parameter values. The former penalizes for
large magnitudes and the latter penalizes for wiggly responses along
covariates or space. In case of covariate responses this is ecologically
justified since according to ecological niche theory species’ responses to
environmental covariates are typically either monotonic or have only one mode
(Ovaskainen et al., 2016). Moreover, Fuglstad et al. (2018) show that these
type of priors work well for spatial random effects with Matérn covariance
functions in general. Favoring longer spatial length-scales is justified also
from the model identifibiality point of view. Identifiability between fixed
and spatial random effects is of less concern if spatial random effect varies
with larger spatial range than the environmental covariates (Paciorek, 2010,
Section 3). Since the environmental covariates in Finnish coastal region vary
with relatively small spatial range, our weakly informative prior effectively
favors models that explain this variation with covariate responses instead of
with spatial random effects. Our model comparison suggests that the models
with both covariate effects and spatial random effects lead to most reliable
inference on environmental responses and spatial random effects. A model’s
extrapolation performance is governed mainly by the environmental covariates
so better extrapolation performance indicates more reliable response
functions. Similarly more accurate interpolation performance indicates more
reliable joint effect of environmental covariates and spatial random effect.

An evident challenge with multivariate additive GP is the computation. The
core element in the multivariate additive GP is the covariance matrix induced
by the GP components. With many species and sampling sites the covariance
matrix can become so large that it makes the implementation infeasible. Here
we used Laplace method to speed up the computation by decreasing the number of
costly posterior density calculations compared to full MCMC. However, in order
to scale up the methods for large data sets involving thousands of sampling
sites, the model would need to be implemented even more efficiently. One
option to reduce the computational time could be to exploit the property that
linear LMC can be parameterized through species specific conditional
distributions (Gelfand et al., 2004). We could also implement multivariate GPs
with sparse GP approximations (Vanhatalo et al., 2010; Alvarez et al., 2010)
and replace the full rank LMC model with latent factor model (Thorson et al.,
2015; Ovaskainen et al., 2017). However, we leave these considerations for
future.

###  7.3 Summary

Our JSDM based on multivariate additive GP combines the key ideas from semi-
parametric SSDMs and state-of-the-art parametric JSDMs. In our case study, it
showed superior predictive performance compared to existing GP based SSDMs and
parametric SSDMs and JSDMs in both interpolation and extrapolation. Hence, the
multivariate additive GP can be seen as the first step towards integration of
the semi-parametric SSDMs and hierarchical JSDMs. We propose also an efficient
approach for inference by utilizing Laplace approximation and gradient based
optimization for hyperparameters. The multivariate additive GP model is not
restricted only to species distribution modeling but can be used in wide
variety of other applications as well.

Supplementary material The supplementary material contains additional
mathematical formulation of the methodology proposed in this paper and
additional figures and tables for the case study analysis.
[url]http://www.some-url-address.org/dowload/0000.zip

This work has bee funded by the Academy of Finland (grant 317255) and Research
Funds of the University of Helsinki (decision No. 465/51/2014).

## References

  * Altartouri et al. (2014) Altartouri, A., Nurminen, L., and Jolma, A. (2014).  “Modeling the role of the close-range effect and environmental variables in the occurrence and spread of Phragmites australis in four sites on the Finnish coast of the Gulf of Finland and the Archipelago Sea.”  Ecology and Evolution, 4: 987–1005.
  * Alvarez et al. (2010) Alvarez, M. A., Luengo, D., Titsias, M. K., and Lawrence., N. D. (2010).  “Efficient multi-output Gaussian processes through variational inducing kernels.”  JMLR Workshop and conference proceedings, 9: 25–32.
  * Banerjee et al. (2015) Banerjee, S., Carlin, B. P., and Gelfand, A. E. (2015).  Hierarchical Modelling and Analysis for Spatial Data.  Chapman Hall/CRC, second edition.
  * Barnard et al. (2000) Barnard, J., McCulloch, R., and Meng, X.-L. (2000).  “Modelling covariance matrices in terms of standard deviations and correlations, with applications to shrinkage.”  Statistical Sinica.
  * Bergström et al. (2015) Bergström, U., Olsson, J., Casini, M., Eriksson, B. K., Fredriksson, R., Wennhage, H., and Appelberg, M. (2015).  “Stickleback increase in the Baltic Sea – A thorny issue for coastal predatory fish.”  Estuarine, Coastal and Shelf Science, 163: 134–142.
  * Bose et al. (2018) Bose, M., Hodges, J. S., and Banerjee, S. (2018).  “Toward a diagnostic toolkit for linear models with Gaussian-process distributed random effects.”  Biometrics, 74(3): 863–873.
  * Busse and Snoeijs (2002) Busse, S. and Snoeijs, P. (2002).  “Gradient responses of diatom communities in the Bothnian Bay, northern Baltic Sea.”  Nova Hedwigia, 3-4: 501–525.
  * Byström et al. (2015) Byström, P., Bergström, U., Hjälten, A., Ståhl, S., Jonsson, D., and Olsson, J. (2015).  “Declining coastal piscivore populations in the Baltic Sea: Where and when do sticklebacks matter?”  Ambio, 44: 462–471.
  * Candolin et al. (2008) Candolin, U., Engström‐Öst, J., and Salesto, T. (2008).  “Human‐induced eutrophication enhances reproductive success through effects on parenting ability in sticklebacks.”  Oikos, 117: 459–465.
  * Chib and Greenberg (1998) Chib, S. and Greenberg, E. (1998).  “Analysis of multivariate probit models.”  Biometrika, 85(2): 347–361.
  * Clark et al. (2014) Clark, J. S., Gelfand, A., Woodall, C. W., and Zhu, K. (2014).  “More than the sum of the parts: forest climate response from joint species distribution models.”  Ecological Applications, 24(5): 990–999.
  * Clark et al. (2016) Clark, J. S., Nemergut, D., Seyednasrollah, B., Turner, P. J., and Zhang, S. (2016).  “Generalized joint attribute modeling for biodiversity analysis: median-zero, multivariate, multifarious data.”  Ecological Monographs, 87(1): 34–56.
  * Cressie and Wikle (2011) Cressie, N. and Wikle, C. K. (2011).  Statistics for Spatial-Temporal Data.  Wiley Series in Probability and Statistics.
  * Dunstan et al. (2013) Dunstan, P. K., Foster, S. D., Hui, F. K. C., and Warton, D. I. (2013).  “Finite Mixture of Regression Modeling for High-Dimensional Count and Biomass Data in Ecology.”  Journal of Agricultural, Biological, and Environmental Statistics, 18(3): 357–375.
  * Duvenaud et al. (2011) Duvenaud, D., Nickisch, H., and Rasmussen, C. (2011).  “Additive Gaussian Processes.”  Neural Information Processing Systems.
  * Elith and Leathwich (2009) Elith, J. and Leathwich, J. R. (2009).  “Species Distributions Models: Ecological Explanation and Predictions Across Space and Time.”  The Annual Review of Ecology, Evolution and Systematics, 40(677-697).
  * Eriksson et al. (2009) Eriksson, B. K., Ljunggren, L., Sandstrom, A., Johansson, G., Mattila, J., Rubach, A. E., Råberg, S., and Snickars, M. (2009).  “Declines in predatory fish promote bloomforming macroalgae.”  Ecological Applications, 19: 1975e1988.
  * Eriksson et al. (2012) Eriksson, B. K., Rubach, A., Batsleer, J., and Hillebrand, H. (2012).  “Cascading predator control interacts with productivity to determine the trophic level of biomass accumulation in a benthic food web.”  Ecological Research, 27: 203e210.
  * Fox and Dunson (2015) Fox, E. B. and Dunson, D. B. (2015).  “Bayesian Nonparametric Covariance Regression.”  Journal of Machine Learning research, 16(1): 2501–2542.
  * Fuglstad et al. (2018) Fuglstad, G.-A., Simpson, D., Lindgren, F., and Rue, H. (2018).  “Constructing priors that penalize the complexity of Gaussian random fields.”  Journal of the American Statistical Association, 1(1): 1–8.
  * Gelfand et al. (2006) Gelfand, A. E., Jr, J. A. S., Wu, S., Latimer, A., Lewis, P. O., Rebelo, A. G., and Holder, M. (2006).  “Explaining Species Distribution Patterns through Hierarchical Modelling.”  Bayesian Analysis, 1(1): 41–92.
  * Gelfand et al. (2004) Gelfand, A. E., Schmidt, A. M., Banerjee, S., and Sirmans, C. F. (2004).  “Nonstationary multivariate process modeling through spatially varying coregionalization.”  Test, 13(2): 263–312.
  * Gelman (2006) Gelman, A. (2006).  “Prior distributions for variance parameters in hierarchical models (comment on article by Browne and Draper).”  Bayesian Analysis, 1(3): 515–534.
  * Gelman et al. (2013) Gelman, A., Carlin, J. B., Stern, H. S., Dunson, D. B., Vehtari, A., and Rubin, D. B. (2013).  Bayesian Data Analysis.  Chapman & Hall/CRC, third edition.
  * Golding and Purse (2016) Golding, N. and Purse, B. V. (2016).  “Fast and flexible Bayesian species distribution modelling using Gaussian processes.”  Methods in Ecology and Evolution, 7: 598–608.
  * Guisan et al. (2002) Guisan, A., Edwards, T. C., and Hastie, T. (2002).  “Generalized linear and generalized additive models in studies of species distributions: setting the scene.”  Ecological Modelling, 157(2-3): 89–100.
  * Guisan et al. (2013) Guisan, A., Tingley, R., Baumgartner, J. B., Naujokaitis-Lewis, I., Sutcliffe, P. R., Tulloch, A. I. T., Regan, T. J., Brotons, L., McDonald-Madden, E., Mantyka-Pringle, C., Martin, T. G., Rhodes, J. R., Maggini, R., Setterfield, S. A., Elith, J., Schwartz, M. W., Wintle, B. A., Broennimann, O., Austin, M., Ferrier, S., Kearney, M. R., Possingham, H. P., and Buckley, Y. M. (2013).  “Predicting species distributions for conservation decisions.”  Ecology Letters, 16(12): 1424–1435.
  * Hartmann et al. (2017) Hartmann, M., Hosack, G. R., Hillary, R. M., and Vanhatalo, J. (2017).  “Gaussian process framework for temporal dependence and discrepancy functions in Ricker-type population growth models.”  Ann. Appl. Statist., 11(3): 1375–1402.
  * Hastie and Tibishirani (1986) Hastie, T. and Tibishirani, R. (1986).  “Generalize additive models.”  Statistical Science, 1(3): 297–318.
  * Hodges and Reich (2010) Hodges, J. S. and Reich, B. J. (2010).  “Adding Spatially-Correlated Errors Can Mess Up the Fixed Effect You Love.”  The American Statistician, 64(4): 325–334.
  * Hudd et al. (1988) Hudd, R., Lehtonen, H., and Kurttila, I. (1988).  “Growth and abundance of fry; factors which influence the year-class strength of whitefish (_Coregonus lavaretus widegreni_) in the southern Bothnian Bay (Baltic).”  Finnish Fisheries Research, 9: 213–220.
  * Hui et al. (2013) Hui, F. K. C., Warton, D. I., Foster, S. D., and Dunstan, P. K. (2013).  “To mix or not to mix: Comparing the predictive performance of mixture models vs. separate species distribution models.”  Ecology, 94(9): 1913–1919.
  * Kallasvuo et al. (2017) Kallasvuo, M., Vanhatalo, J., and Veneranta, L. (2017).  “Modeling the spatial distribution of larval fish abundance provides essential information for management.”  Canadian Journal of Fisheries and Aquatic Sciences, 74: 636–649.
  * Kass and Raftery (1995) Kass, R. E. and Raftery, A. E. (1995).  “Bayes Factors.”  Journal of the American Statistical Association, 90(430): 773–795.
  * Kotta et al. (2019) Kotta, J., Vanhatalo, J., Jänes, H., Orav-Kotta, H., Rugiu, L., Jormalainen, V., Bobsien, I., Viitasalo, M., Virtanen, E., Sandman, A. N., Isaeus, M., Leidenberger, S., Jonsson, P., and Johannesson, K. (2019).  “Integrating experimental and distribution data to predict future species patterns.”  Scientific Reports, 9(1): 1821.
  * Kurowicka and Cooke (2003) Kurowicka, D. and Cooke, R. (2003).  “A parameterization of positive definite matrices in terms of partial correlation vines.”  Linear Algebra and its Applications, 372(Supplement C): 225–251.
  * Latimer et al. (2006) Latimer, A. M., Wu, S., Gelfand, A. E., and Silander, Jr., J. A. (2006).  “Building Statistical Models to Analyze Species Distributions.”  Ecological Applications, 16(1): 33–50.
  * Lefébure et al. (2011) Lefébure, R., Larsson, S., and Byström, P. (2011).  “A temperature dependent growth model for the three‐spined stickleback _Gasterosteus aculeatus_.”  Journal of fish biology, 79: 1815–1827.
  * Leskelä et al. (1991) Leskelä, A., Hudd, R., Lehtonen, H., Huhmarniemi, A., and Sandström, O. (1991).  “Habitats of whitefish (_Coregonus lavaretus (L.) s.l._) larvae in the Gulf of Bothnia.”  Aqua Fennica, 21: 145–151.
  * Lewandowski et al. (2009) Lewandowski, D., Kurowicka, D., and Joe, H. (2009).  “Generating random correlation matrices based on vines and extended onion method.”  Journal of Multivariate Analysis, 100(9): 1989–2001.
  * Lindley (2002) Lindley, D. V. (2002).  “Seeing and Doing: The Concept of Causation.”  International Statistical Review, 70(2): 191–214.
  * Liu and Vanhatalo (2018) Liu, J. and Vanhatalo, J. (2018).  “Bayesian model-based spatiotemporal survey design for log-Gaussian Cox process.”  arXiv e-prints, arXiv:1808.09200.
  * Lopez (2000) Lopez, H. F. (2000).  “Bayesian Analysis in Latent Factor and Longitudinal Models.”  Ph.D. thesis, Institute of Statistics and Decision Sciences, Duke University.
  * Lopez et al. (2008) Lopez, H. F., Salazar, E., and Gamerman, D. (2008).  “Spatial Dynamic Factor Analysis.”  Bayesian Statistics, 3(4): 759–792.
  * MacKenzie et al. (2007) MacKenzie, B. R., Gislason, H., Möllmann, C., and Köster, F. W. (2007).  “Impact of 21st century climate change on the Baltic Sea fish community and fisheries.”  Global Change Biology, 13: 1348–1367.
  * Mardia and Goodall (1993) Mardia, K. V. and Goodall, C. R. (1993).  “Spatio-Temporal analysis of Multivariate Environmental Monitoring Data.”  Multivariate Environmental Statistics, 347–386.
  * M.Cingi et al. (2010) M.Cingi, S., Keinänen, and Vuorinen, P. J. (2010).  “Elevated water temperature impairs fertilization and embryonic development of whitefish _Coregonus lavaretus_.”  Journal of Fish Biology, 76: 502–521.
  * Meier et al. (2012) Meier, H. E. M., Hordoir, R., Andersson, H. C., Dieterich, C., Eilola, K., …, B. G. G., and Schimanke, S. (2012).  “Modeling the combined impact of changing climate and changing nutrient loads on the Baltic Sea environment in an ensemble of transient simulations for 1961–2099.”  Climate Dynamics, 39: 2421–2441.
  * Müller (1992) Müller, R. (1992).  “Trophic state and its implications for natural reproduction of salmonid fish.”  In Dynamics and Use of Lacustrine Ecosystems. Springer, Dordrecht.
  * Nelder and Wedderburn (1972) Nelder, J. A. and Wedderburn, R. W. M. (1972).  “Generalized Linear Models.”  Journal of the Royal Statistical Association, 135(3): 370–384.
  * Nickisch and Rasmussen (2008) Nickisch, H. and Rasmussen, C. E. (2008).  “Approximations for binary Gaussian process classification.”  Journal of Machine learning, 9: 2035–2078.
  * Odum (1953) Odum, E. P. (1953).  Fundamentals of ecology.  Wiley Subscription Services, Inc., A Wiley Company.
  * O’Hagan (1978) O’Hagan, A. (1978).  “Curve Fitting and Optimal Design for Prediction.”  Journal of Royal Statistical Society B, 40(1): 1–42.
  * Ovaskainen et al. (2016) Ovaskainen, O., de Knegt, H. J., and del Mar Delgado, M. (2016).  Quantitative Ecology and Evolutionary Biology -Integrating models with data.  Oxford University Press.
  * Ovaskainen et al. (2010) Ovaskainen, O., Hottola, J., and Siitonen, J. (2010).  “Modeling species co-occurence by multivariate logistic regression generates new hypotheses on fungal interaction.”  Ecology, 91(9): 2514–2521.
  * Ovaskainen and Soininen (2011) Ovaskainen, O. and Soininen, J. (2011).  “Making more out of sparse data: Hierarchical modelling of species communities.”  Ecology, 92(2): 289–295.
  * Ovaskainen et al. (2017) Ovaskainen, O., Tikhonov, G., Norberg, A., Guillaume Blanchet, F., Duan, L., Dunson, D., Roslin, T., and Abrego, N. (2017).  “How to make more out of community data? A conceptual framework and its implementation as models and software.”  Ecology Letters, 20(5): 561–576.
  * Paciorek (2010) Paciorek, C. J. (2010).  “The Importance of Scale for Spatial-Confounding Bias and Precision of Spatial Regression Estimators.”  Statistical Science, 25(1): 107–125.
  * Pearl (2009) Pearl, J. (2009).  “Causal inference in statistics: An overview.”  Statistics Surveys, 3: 96–146.
  * Peltonen et al. (2004) Peltonen, H., Vinni, M., Lappalainen, A., and Pönni, J. (2004).  “Spatial feeding patterns of herring (_Clupea harengus_ L.), sprat (_Sprattus sprattus_ L.), and the three-spined stickleback (emphGasterosteus aculeatus L.) in the Gulf of Finland, Baltic Sea.”  ICES Journal of Marine Science, 61: 966–971.
  * Pitkänen et al. (2013) Pitkänen, H., Peuraniemi, M., Westerbom, M., Kilpi, M., and von Numers, M. (2013).  “Longterm changes in distribution and frequency of aquatic vascular plants and charophytes in an estuary in the Baltic Sea.”  Annales Botanica Fennica, 50: 1e54.
  * Pollock et al. (2014) Pollock, L. J., Tingley, R., Morris, W. K., Golding, N., Hara, R. B. O., Parris, K. M., Vesk, P. A., and McCarthy, M. A. (2014).  “Understanding co-occurence by modelling species simultaneously with joint species distribution model (JSDM).”  Methods in Ecology and Evolution, 5: 397–406.
  * Rasmussen and Williams (2006) Rasmussen, C. E. and Williams, C. K. I. (2006).  Gaussian Processes for Machine Learning.  The MIT Press.
  * Record et al. (2013) Record, S., Fitzpatrick, M. C., Finley, A. O., Veloz, S., and Ellison, A. M. (2013).  “Should species distribution models account for spatial autocorrelation? A test of model projections across eight millennia of climate change.”  Global Ecology and Biogeography, 22(6): 760–771.
  * Reusch et al. (2018) Reusch, T. B. H., Dierking, J., Andersson, H. C., Bonsdorff, E., Carstensen, J., Casini, M., Czajkowski, M., Hasler, B., Hinsby, K., Hyytiäinen, K., Johannesson, K., Jomaa, S., Jormalainen, V., Kuosa, H., Kurland, S., Laikre, L., MacKenzie, B. R., Margonski, P., Melzner, F., Oesterwind, D., Ojaveer, H., Refsgaard, J. C., Sandström, A., Schwarz, G., Tonderski, K., Winder, M., and Zandersen, M. (2018).  “The Baltic Sea as a time machine for the future coastal ocean.”  Science Advances, 4(5).
  * Rönnberg and Bonsdorff (2004) Rönnberg, C. and Bonsdorff, E. (2004).  “Baltic Sea eutrophication: area-specific ecological consequences.”  Hydrobiologia, 514: 227–241.
  * Rue and Marino (2009) Rue, H. and Marino, S. (2009).  “Approximate Bayesian Inference for latent Gaussian models by using integrated Laplace approximantions.”  Journal of the Royal Statistical Society, 71(2): 319–392.
  * Shelton et al. (2014) Shelton, A. O., Thorson, J. T., Ward, E. J., and Feist, B. E. (2014).  “Spatial semiparametric models improve estimates of species abundance and distribution.”  Canadian Journal of Fisheries and Aquatic Sciences, 71(July): 1655–1666.
  * Sieben et al. (2011) Sieben, K., Ljunggren, L., Bergström, U., and Eriksson, B. K. (2011).  “A meso-predator release of stickleback promotes recruitment of macroalgae in the Baltic Sea.”  Journal of Experimental Marine Biology and Ecology, 397: 79e84.
  * Simpson et al. (2017) Simpson, D. P., Rue, H. v., Riebler, A., Martins, T. G., and Sørbye, S. H. (2017).  “Penalising model component complexity: A principled, practical approach to constructing priors.”  Statistical Science, 32(1): 1–28.
  * Taylor-Rodríguez et al. (2017) Taylor-Rodríguez, D., Kaufeld, K., Schliep, E. M., Clark, J. S., and Gelfand, A. E. (2017).  “Joint Species Distribution Modeling: Dimension Reduction Using Dirichlet Processes.”  Bayesian Analysis, 12(4): 939–967.
  * Thorson et al. (2015) Thorson, J. T., Scheuerell, M. D., Shelton, A. O., See, K. E., Skaug, H. J., and Kristensen, K. (2015).  “Spatial factor analysis: a new tool for estimating joint species distributions and correlations in species range.”  Methods in Ecology and Evolution, 6(6): 627–637.
  * Tikhonov et al. (2017) Tikhonov, G., Abrego, N., Dunson, D., and Ovaskainen, O. (2017).  “Using joint species distribution models for evaluating how species-to-species associations depend on the environmental context.”  Methods in Ecology and Evolution, 8(4): 443–452.
  * Tokuda et al. (2012) Tokuda, T., Goodrich, B., Mechelen, I. V., and Gelman, A. (2012).  “Visualizing Distributions of Covariance Matrices.”
  * Vanhatalo et al. (2010) Vanhatalo, J., Pietiläinen, V., and Vehtari, A. (2010).  “Approximate inference for disease mapping with sparse Gaussian processes.”  Statistics in Medicine, 29(15): 1580–1607.
  * Vanhatalo et al. (2013) Vanhatalo, J., Riihimäki, J., Hartikainen, J., Jylänki, P., Tolvanen, V., and Vehtari, A. (2013).  “GPstuff : Bayesian Modeling with Gaussian Processes.”  Journal of Machine Learning Research, 14: 1175–1179.
  * Vanhatalo et al. (2012) Vanhatalo, J., Veneranta, L., and Hudd, R. (2012).  “Species distribution modelling with Gaussian processes: A case study with youngest stages of sea spawning whitefish (Coregonus lavatus L. s.l.) larvae.”  Ecological Modelling, (228): 49–58.
  * Vehtari et al. (2016) Vehtari, A., Mononen, T., Tolvanen, V., Sivula, T., and Winther, O. (2016).  “Bayesian Leave-One-Out Cross-Validation Approximations for Gaussian Latent Variable Models.”  Journal of Machine Learning Research, 17: 1–38.
  * Vehtari and Ojanen (2012) Vehtari, A. and Ojanen, J. (2012).  “A survey of Bayesian predictive methods for model assessment, selection and comparison.”  Statistics Surveys, 6: 141–228.
  * Veneranta et al. (2013) Veneranta, L., Hudd, R., and Vanhatalo, J. (2013).  “Reproduction areas of sea-spawning coregionids reflect the environment in shallow coastal areas.”  Marine Ecology Progress Series, 477: 231–250.
  * Warton et al. (2015) Warton, D. I., Blanchet, F. G., O’Hara, R. B., Ovaskainen, O., Taskinen, S., Walker, S. C., and Hui, F. K. (2015).  “So Many Variables: Joint Modeling in Community Ecology.”  Trends in Ecology and Evolution, 30(12): 766–779.
  * Warton and Shepherd (2010) Warton, D. I. and Shepherd, L. C. (2010).  “Poisson point process models solve the ”pseudo-absence problem” for presence-only data in ecology.”  Annals of Applied Statistics, 4(3): 1383–1402.
  * Wikle (2003) Wikle, C. K. (2003).  “Hierarchical Models in Environmental Science.”  International Statistical Review, 71(2): 181–199.
  * Zhang (2004) Zhang, H. (2004).  “Inconsistent Estimation and Asymptotically Equal Interpolations in Model-Based Geostatistics.”  Journal of the American Statistical Association, 99(465): 250–261.