  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. 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. Current JSDMs are restricted to hierarchical generalized linear
modeling framework. These parametric models have trouble in explaining changes
in abundance due, e.g., 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.
