Nothing
Get started with jSDM
In jSDM: Joint Species Distribution Models
knitr::opts_chunk$set(
fig.align = "center",
fig.retina = 1,
fig.width = 6, fig.height = 6,
cache = FALSE,
collapse = TRUE,
comment = "#>",
highlight = TRUE
)
jSDM package
jSDM is an R package for fitting joint species distribution models (JSDM) in a hierarchical Bayesian framework.
The Gibbs sampler is written in C++. It uses Rcpp, Armadillo and GSL to maximize computation efficiency.
|Package: |jSDM |
|:------|:-----|
| Type: |Package|
|Version: |0.2.1|
|Date: |2019-01-11|
|License:|GPL-3|
|LazyLoad:|yes|
The package includes the following functions to fit various species distribution models :
| function | data type | data format |
|:---------|:-------------|:------------:|
| jSDM_binomial_logit() | presence-absence | wide |
| jSDM_binomial_probit() | presence-absence | wide |
| jSDM_binomial_probit_sp_constrained() | presence-absence | wide |
| jSDM_binomial_probit_long_format() | presence-absence | long |
| jSDM_poisson_log() | abundance | wide |
| jSDM_gaussian() | continuous | wide |
jSDM_binomial_probit() :
Ecological process:
$$y_{ij} \sim \mathcal{B}ernoulli(\theta_{ij}),$$
where
| conditions | model specification |
|:---------------------------------------------|:-------------------------------------------------|
| if n_latent=0 and site_effect="none" | probit$(\theta_{ij}) = \beta_{0j} + X_i \beta_j$ |
| if n_latent>0 and site_effect="none" | probit$(\theta_{ij}) = \beta_{0j} + X_i \beta_j + W_i \lambda_j$ |
| if n_latent=0 and site_effect="random" | probit$(\theta_{ij}) = \beta_{0j} + X_i \beta_j + \alpha_i$ and $\alpha_i \sim \mathcal{N}(0,V_\alpha)$ |
| if n_latent>0 and site_effect="fixed" | probit$(\theta_{ij}) = \beta_{0j} + X_i \beta_j + W_i \lambda_j + \alpha_i$ |
| if n_latent=0 and site_effect="fixed" | probit$(\theta_{ij}) = \beta_{0j} + X_i \beta_j + \alpha_i$ |
| if n_latent>0 and site_effect="random" | probit$(\theta_{ij}) = \beta_{0j} + X_i \beta_j + W_i \lambda_j + \alpha_i$ and $\alpha_i \sim \mathcal{N}(0,V_\alpha)$ |
jSDM_binomial_probit_sp_constrained() :
This function allows to fit the same models than the function jSDM_binomial_probit except for models not including latent variables, indeed n_latent must be greater than zero in this function.
At first, the function fit a JSDM with the constrained species arbitrarily chosen as the first ones in the presence-absence data-set.
Then, the function evaluates the convergence of MCMC $\lambda$ chains using the Gelman-Rubin convergence diagnostic [@Gelman1992] ($\hat{R}$). It identifies the species ($\hat{j}l$) having the higher $\hat{R}$ for $\lambda{\hat{j}_l}$. These species drive the structure of the latent axis $l$. The $\lambda$ corresponding to this species are constrained to be positive and placed on the diagonal of the $\Lambda$ matrix for fitting a second model.
This usually improves the convergence of the latent variables and factor loadings. The function returns the parameter posterior distributions for this second model.
jSDM_binomial_logit() :
Ecological process :
$$y_{ij} \sim \mathcal{B}inomial(\theta_{ij},t_i),$$
where
| conditions | model specification |
|:----------------------------------------|:-----------------------------------------------------|
| if n_latent=0 and site_effect="none" | logit$(\theta_{ij}) = \beta_{0j} + X_i \beta_j$ |
| if n_latent>0 and site_effect="none" | logit$(\theta_{ij}) = \beta_{0j} + X_i \beta_j + W_i \lambda_j$ |
| if n_latent=0 and site_effect="fixed" | logit$(\theta_{ij}) = \beta_{0j} + X_i \beta_j + \alpha_i$ |
| if n_latent>0 and site_effect="fixed" | logit$(\theta_{ij}) = \beta_{0j} + X_i \beta_j + W_i \lambda_j + \alpha_i$ |
| if n_latent=0 and site_effect="random" | logit$(\theta_{ij}) = \beta_{0j} + X_i \beta_j + \alpha_i$ and $\alpha_i \sim \mathcal{N}(0,V_{\alpha})$|
| if n_latent>0 and site_effect="random" | logit$(\theta_{ij}) = \beta_{0j} + X_i \beta_j + W_i \lambda_j + \alpha_i$ and $\alpha_i \sim \mathcal{N}(0,V_\alpha)$ |
jSDM_poisson_log() :
Ecological process :
$$y_{ij} \sim \mathcal{P}oisson(\theta_{ij}),$$
where
| conditions | model specification |
|:-----------------------------------------|:-----------------------------------------------------|
| if n_latent=0 and site_effect="none" | log$(\theta_{ij}) = \beta_{0j} + X_i \beta_j$ |
| if n_latent>0 and site_effect="none" | log$(\theta_{ij}) = \beta_{0j} + X_i \beta_j + W_i \lambda_j$ |
| if n_latent=0 and site_effect="fixed"| log$(\theta_{ij}) = \beta_{0j} + X_i \beta_j + \alpha_i$ |
| if n_latent>0 and site_effect="fixed"| log$(\theta_{ij}) = \beta_{0j} + X_i \beta_j + W_i \lambda_j + \alpha_i$ |
| if n_latent=0 and site_effect="random"| log$(\theta_{ij}) = \beta_{0j} + X_i \beta_j + \alpha_i$ and $\alpha_i \sim \mathcal{N}(0,V_\alpha)$ |
| if n_latent>0 and site_effect="random"| log$(\theta_{ij}) = \beta_{0j} + X_i \beta_j + W_i \lambda_j + \alpha_i$ and $\alpha_i \sim \mathcal{N}(0,V_\alpha)$ |
jSDM_binomial_probit_long_format() :
Ecological process:
$$y_{n} \sim \mathcal{B}ernoulli(\theta_n),$$
such as $species_n=j$ and $site_n=i$,
where
| conditions | model specification |
|:---------------------------------------------|:-------------------------------------------------|
| if n_latent=0 and site_effect="none" | probit$(\theta_n) = D_n \gamma + X_n \beta_j$ |
| if n_latent>0 and site_effect="none" | probit$(\theta_n) = D_n \gamma + X_n \beta_j + W_i \lambda_j$ |
| if n_latent=0 and site_effect="random" | probit$(\theta_n) = D_n \gamma + X_n \beta_j + \alpha_i$ and $\alpha_i \sim \mathcal{N}(0,V_\alpha)$ |
| if n_latent>0 and site_effect="fixed" | probit$(\theta_n) = D_n \gamma + X_n \beta_j + W_i \lambda_j + \alpha_i$ |
| if n_latent=0 and site_effect="fixed" | probit$(\theta_n) = D_n \gamma + X_n \beta_j + \alpha_i$ |
| if n_latent>0 and site_effect="random" | probit$(\theta_n) = D_n \gamma + X_n \beta_j + W_i \lambda_j + \alpha_i$ and $\alpha_i \sim \mathcal{N}(0,V_\alpha)$ |
jSDM_gaussian() :
Ecological process:
$$y_{ij} \sim \mathcal{N}(\theta_{ij}, V),$$
where
| conditions | model specification |
|:---------------------------------------------|:-------------------------------------------------|
| if n_latent=0 and site_effect="none" | $\theta_{ij} = \beta_{0j} + X_i \beta_j$ |
| if n_latent>0 and site_effect="none" | $\theta_{ij} = \beta_{0j} + X_i \beta_j + W_i \lambda_j$ |
| if n_latent=0 and site_effect="random" | $\theta_{ij} = \beta_{0j} + X_i \beta_j + \alpha_i$ and $\alpha_i \sim \mathcal{N}(0,V_\alpha)$ |
| if n_latent>0 and site_effect="fixed" | $\theta_{ij} = \beta_{0j} + X_i \beta_j + W_i \lambda_j + \alpha_i$ |
| if n_latent=0 and site_effect="fixed" | $\theta_{ij} = \beta_{0j} + X_i \beta_j + \alpha_i$ |
| if n_latent>0 and site_effect="random" | $\theta_{ij} = \beta_{0j} + X_i \beta_j + W_i \lambda_j + \alpha_i$ and $\alpha_i \sim \mathcal{N}(0,V_\alpha)$ |
Joint Species distribution models (jSDM) are useful tools to explain or predict species range and abundance from various environmental factors and species correlations [@Warton2015]. jSDM is becoming an increasingly popular statistical method in conservation biology.
In this vignette, we illustrate the use of the jSDM R package which aims at providing user-friendly statistical functions using field observations (occurrence or abundance data) to fit jSDMs models.
Package's functions are developed in a hierarchical Bayesian framework and use adaptive rejection Metropolis sampling algorithms or conjugate priors within Gibbs sampling to estimate model's parameters. Using compiled C++ code for the Gibbs sampler reduce drastically the computation time. By making these new statistical tools available to the scientific community, we hope to democratize the use of more complex, but more realistic, statistical models for increasing knowledge in ecology and conserving biodiversity.
Directed Acyclic Graph (DAG)
(ref:cap-DAG) A graphical summary of the jSDM-package statistical framework. In this Directed Acyclic Graph (DAG), the orange boxes refer to data, the blue ellipses to parameters to be estimated, and the arrows to functional relationships described with the help of statistical distributions.
knitr::include_graphics("figures/DAG-jSDM-rand.png")
Model types available in jSDM R package are not limited to those described in this example. jSDM includes various model types for occurrence and abundance data, you can find more examples of use on the jSDM website.
Load librairies
We first load the jSDM library.
# Load libraries
library(jSDM)
Bernoulli probit regression
Below, we show an example of the use of jSDM-package for fitting species distribution model to occurence data for 9 frog's species.
Definition of the model
Referring to the models used in the articles @Warton2015 and @Albert1993, we define the following model :
$$ \mathrm{probit}(\theta_{ij}) =\alpha_i + \beta_{0j}+X_i.\beta_j+ W_i.\lambda_j $$
-
Link function probit: $\mathrm{probit}: q \rightarrow \Phi^{-1}(q)$ where $\Phi$ correspond to the repartition function of the reduced centred normal distribution.
-
Response variable: $Y=(y_{ij})^{i=1,\ldots,nsite}_{j=1,\ldots,nsp}$ with:
$$y_{ij}=\begin{cases}
0 & \text{ if species $j$ is absent on the site $i$}\
1 & \text{ if species $j$ is present on the site $i$}.
\end{cases}$$
- Latent variable $z_{ij} = \alpha_i + \beta_{0j} + X_i.\beta_j + W_i.\lambda_j + \epsilon_{i,j}$, with $\forall (i,j) \ \epsilon_{ij} \sim \mathcal{N}(0,1)$ and such that:
$$y_{ij}=\begin{cases}
1 & \text{if} \ z_{ij} > 0 \
0 & \text{otherwise.}
\end{cases}$$
It can be easily shown that: $y_{ij} \sim \mathcal{B}ernoulli(\theta_{ij})$.
-
Latent variables: $W_i=(W_i^1,\ldots,W_i^q)$ where $q$ is the number of latent variables considered, which has to be fixed by the user (by default q=2).
We assume that $W_i \sim \mathcal{N}(0,I_q)$ and we define the associated coefficients: $\lambda_j=(\lambda_j^1,\ldots, \lambda_j^q)'$. We use a prior distribution $\mathcal{N}(0,1)$ for all lambdas not concerned by constraints to $0$ on upper diagonal and to strictly positive values on diagonal.
-
Explanatory variables: bioclimatic data about each site. $X=(X_i)_{i=1,\ldots,nsite}$ with $X_i=(x_i^1,\ldots,x_i^p)\in \mathbb{R}^p$ where $p$ is the number of bioclimatic variables considered.
The corresponding regression coefficients for each species $j$ are noted : $\beta_j=(\beta_j^1,\ldots,\beta_j^p)'$.
-
$\beta_{0j}$ correspond to the intercept for species $j$ which is assume to be a fixed effect. We use a prior distribution $\mathcal{N}(0,1)$ for all betas.
-
$\alpha_i$ represents the random effect of site $i$ such as $\alpha_i \sim \mathcal{N}(0,V_{\alpha})$ and we assumed that $V_{\alpha} \sim \mathcal {IG}(\text{shape}=0.1, \text{rate}=0.1)$ as prior distribution by default.
Occurrence data-set
(ref:cap-frog) Litoria ewingii [@Wilkinson2019].
knitr::include_graphics("figures/Litoria_ewingii.jpg")
This data-set is available in jSDM-package. It can be loaded with the data() command. The frogs dataset is in "wide" format: each line is a site and the occurrence data (from Species_1 to Species_9) are in columns. A site is characterized by its x-y geographical coordinates, one discrete covariate and two other continuous covariates.
# frogs data
data(frogs, package="jSDM")
head(frogs)
We rearrange the data in two data-sets: a first one for the presence-absence observations for each species (columns) at each site (rows), and a second one for the site characteristics.
We also normalize the continuous explanatory variables to facilitate MCMC convergence.
# data.obs
PA_frogs <- frogs[,4:12]
# Normalized continuous variables
Env_frogs <- cbind(scale(frogs[,1]),frogs[,2],scale(frogs[,3]))
colnames(Env_frogs) <- colnames(frogs[,1:3])
Parameter inference
We use the jSDM_binomial_probit() function to fit the jSDM (increase the number of iterations to achieve convergence).
mod_frogs_jSDM_probit <- jSDM_binomial_probit(
# Chains
burnin=1000, mcmc=1000, thin=1,
# Response variable
presence_data = PA_frogs,
# Explanatory variables
site_formula = ~.,
site_data = Env_frogs,
# Model specification
n_latent=2, site_effect="random",
# Starting values
alpha_start=0, beta_start=0,
lambda_start=0, W_start=0,
V_alpha=1,
# Priors
shape_Valpha=0.1,
rate_Valpha=0.1,
mu_beta=0, V_beta=1,
mu_lambda=0, V_lambda=1,
# Various
seed=1234, verbose=1)
Analysis of the results
We visually evaluate the convergence of MCMCs by representing the trace and density a posteriori of some estimated parameters.
np <- nrow(mod_frogs_jSDM_probit$model_spec$beta_start)
oldpar <- par(no.readonly = TRUE)
## beta_j of the first two species
par(mfrow=c(2,2))
for (j in 1:2) {
for (p in 1:np) {
coda::traceplot(coda::as.mcmc(mod_frogs_jSDM_probit$mcmc.sp[[j]][,p]))
coda::densplot(coda::as.mcmc(mod_frogs_jSDM_probit$mcmc.sp[[j]][,p]),
main = paste(colnames(mod_frogs_jSDM_probit$mcmc.sp[[j]])[p],
", species : ",j), cex.main=0.9)
}
}
## lambda_j of the first two species
n_latent <- mod_frogs_jSDM_probit$model_spec$n_latent
par(mfrow=c(2,2))
for (j in 1:2) {
for (l in 1:n_latent) {
coda::traceplot(coda::as.mcmc(mod_frogs_jSDM_probit$mcmc.sp[[j]][,np+l]))
coda::densplot(coda::as.mcmc(mod_frogs_jSDM_probit$mcmc.sp[[j]][,np+l]),
main = paste(colnames(mod_frogs_jSDM_probit$mcmc.sp[[j]])
[np+l],", species : ",j), cex.main=0.9)
}
}
## Latent variables W_i for the first two sites
par(mfrow=c(2,2))
for (l in 1:n_latent) {
for (i in 1:2) {
coda::traceplot(mod_frogs_jSDM_probit$mcmc.latent[[paste0("lv_",l)]][,i],
main = paste0("Latent variable W_", l, ", site ", i),
cex.main=0.9)
coda::densplot(mod_frogs_jSDM_probit$mcmc.latent[[paste0("lv_",l)]][,i],
main = paste0("Latent variable W_", l, ", site ", i),
cex.main=0.9)
}
}
## alpha_i of the first two sites
plot(coda::as.mcmc(mod_frogs_jSDM_probit$mcmc.alpha[,1:2]))
## V_alpha
par(mfrow=c(2,2))
coda::traceplot(mod_frogs_jSDM_probit$mcmc.V_alpha)
coda::densplot(mod_frogs_jSDM_probit$mcmc.V_alpha)
## Deviance
coda::traceplot(mod_frogs_jSDM_probit$mcmc.Deviance)
coda::densplot(mod_frogs_jSDM_probit$mcmc.Deviance)
## probit_theta
par (mfrow=c(1,2))
hist(mod_frogs_jSDM_probit$probit_theta_latent,
main = "Predicted probit theta",
xlab ="predicted probit theta")
hist(mod_frogs_jSDM_probit$theta_latent,
main = "Predicted theta",
xlab ="predicted theta")
par(oldpar)
Overall, the traces and the densities of the parameters indicate the convergence of the algorithm. Indeed, we observe on the traces that the values oscillate around averages without showing an upward or downward trend and we see that the densities are quite smooth and for the most part of Gaussian form.
Matrice of correlations
After fitting the jSDM with latent variables, the full species residual correlation matrix $R=(R_{ij})^{i=1,\ldots, n_{species}}{j=1,\ldots, n{species}}$ can be derived from the covariance in the latent variables such as :
$$\Sigma_{ij} = \lambda_i^T .\lambda_j $$, then we compute correlations from covariances :
$$R_{i,j} = \frac{\Sigma_{ij}}{\sqrt{\Sigma {ii}\Sigma {jj}}}$$.
We use the function plot_residual_cor() to compute and display the residual correlation matrix between species :
plot_residual_cor(mod_frogs_jSDM_probit)
Predictions
We use the predict.jSDM() S3 method on the mod_frogs_jSDM_probit object of class jSDM to compute the mean (or expectation) of the posterior distributions obtained and get the expected values of model's parameters.
# Sites and species concerned by predictions :
## 50 sites among the 104
Id_sites <- sample.int(nrow(PA_frogs), 50)
## All species
Id_species <- colnames(PA_frogs)
# Simulate new observations of covariates on those sites
simdata <- matrix(nrow=50, ncol = ncol(mod_frogs_jSDM_probit$model_spec$site_data))
colnames(simdata) <- colnames(mod_frogs_jSDM_probit$model_spec$site_data)
rownames(simdata) <- Id_sites
simdata <- as.data.frame(simdata)
simdata$Covariate_1 <- rnorm(50)
simdata$Covariate_3 <- rnorm(50)
simdata$Covariate_2 <- rbinom(50,1,0.5)
# Predictions
theta_pred <- predict(mod_frogs_jSDM_probit, newdata=simdata, Id_species=Id_species,
Id_sites=Id_sites, type="mean")
hist(theta_pred, main="Predicted theta with simulated data", xlab="predicted theta")
References
Try the jSDM package in your browser
Any scripts or data that you put into this service are public.
jSDM documentation built on July 26, 2023, 5:21 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
knitr::opts_chunk$set( fig.align = "center", fig.retina = 1, fig.width = 6, fig.height = 6, cache = FALSE, collapse = TRUE, comment = "#>", highlight = TRUE )
jSDM package
jSDM is an R package for fitting joint species distribution models (JSDM) in a hierarchical Bayesian framework.
The Gibbs sampler is written in C++. It uses Rcpp, Armadillo and GSL to maximize computation efficiency.
|Package: |jSDM | |:------|:-----| | Type: |Package| |Version: |0.2.1| |Date: |2019-01-11| |License:|GPL-3| |LazyLoad:|yes|
The package includes the following functions to fit various species distribution models :
| function | data type | data format |
|:---------|:-------------|:------------:|
| jSDM_binomial_logit() | presence-absence | wide |
| jSDM_binomial_probit() | presence-absence | wide |
| jSDM_binomial_probit_sp_constrained() | presence-absence | wide |
| jSDM_binomial_probit_long_format() | presence-absence | long |
| jSDM_poisson_log() | abundance | wide |
| jSDM_gaussian() | continuous | wide |
jSDM_binomial_probit():
Ecological process:
$$y_{ij} \sim \mathcal{B}ernoulli(\theta_{ij}),$$ where
| conditions | model specification |
|:---------------------------------------------|:-------------------------------------------------|
| if n_latent=0 and site_effect="none" | probit$(\theta_{ij}) = \beta_{0j} + X_i \beta_j$ |
| if n_latent>0 and site_effect="none" | probit$(\theta_{ij}) = \beta_{0j} + X_i \beta_j + W_i \lambda_j$ |
| if n_latent=0 and site_effect="random" | probit$(\theta_{ij}) = \beta_{0j} + X_i \beta_j + \alpha_i$ and $\alpha_i \sim \mathcal{N}(0,V_\alpha)$ |
| if n_latent>0 and site_effect="fixed" | probit$(\theta_{ij}) = \beta_{0j} + X_i \beta_j + W_i \lambda_j + \alpha_i$ |
| if n_latent=0 and site_effect="fixed" | probit$(\theta_{ij}) = \beta_{0j} + X_i \beta_j + \alpha_i$ |
| if n_latent>0 and site_effect="random" | probit$(\theta_{ij}) = \beta_{0j} + X_i \beta_j + W_i \lambda_j + \alpha_i$ and $\alpha_i \sim \mathcal{N}(0,V_\alpha)$ |
jSDM_binomial_probit_sp_constrained():
This function allows to fit the same models than the function jSDM_binomial_probit except for models not including latent variables, indeed n_latent must be greater than zero in this function.
At first, the function fit a JSDM with the constrained species arbitrarily chosen as the first ones in the presence-absence data-set.
Then, the function evaluates the convergence of MCMC $\lambda$ chains using the Gelman-Rubin convergence diagnostic [@Gelman1992] ($\hat{R}$). It identifies the species ($\hat{j}l$) having the higher $\hat{R}$ for $\lambda{\hat{j}_l}$. These species drive the structure of the latent axis $l$. The $\lambda$ corresponding to this species are constrained to be positive and placed on the diagonal of the $\Lambda$ matrix for fitting a second model.
This usually improves the convergence of the latent variables and factor loadings. The function returns the parameter posterior distributions for this second model.
jSDM_binomial_logit():
Ecological process :
$$y_{ij} \sim \mathcal{B}inomial(\theta_{ij},t_i),$$ where
| conditions | model specification |
|:----------------------------------------|:-----------------------------------------------------|
| if n_latent=0 and site_effect="none" | logit$(\theta_{ij}) = \beta_{0j} + X_i \beta_j$ |
| if n_latent>0 and site_effect="none" | logit$(\theta_{ij}) = \beta_{0j} + X_i \beta_j + W_i \lambda_j$ |
| if n_latent=0 and site_effect="fixed" | logit$(\theta_{ij}) = \beta_{0j} + X_i \beta_j + \alpha_i$ |
| if n_latent>0 and site_effect="fixed" | logit$(\theta_{ij}) = \beta_{0j} + X_i \beta_j + W_i \lambda_j + \alpha_i$ |
| if n_latent=0 and site_effect="random" | logit$(\theta_{ij}) = \beta_{0j} + X_i \beta_j + \alpha_i$ and $\alpha_i \sim \mathcal{N}(0,V_{\alpha})$|
| if n_latent>0 and site_effect="random" | logit$(\theta_{ij}) = \beta_{0j} + X_i \beta_j + W_i \lambda_j + \alpha_i$ and $\alpha_i \sim \mathcal{N}(0,V_\alpha)$ |
jSDM_poisson_log():
Ecological process :
$$y_{ij} \sim \mathcal{P}oisson(\theta_{ij}),$$
where
| conditions | model specification |
|:-----------------------------------------|:-----------------------------------------------------|
| if n_latent=0 and site_effect="none" | log$(\theta_{ij}) = \beta_{0j} + X_i \beta_j$ |
| if n_latent>0 and site_effect="none" | log$(\theta_{ij}) = \beta_{0j} + X_i \beta_j + W_i \lambda_j$ |
| if n_latent=0 and site_effect="fixed"| log$(\theta_{ij}) = \beta_{0j} + X_i \beta_j + \alpha_i$ |
| if n_latent>0 and site_effect="fixed"| log$(\theta_{ij}) = \beta_{0j} + X_i \beta_j + W_i \lambda_j + \alpha_i$ |
| if n_latent=0 and site_effect="random"| log$(\theta_{ij}) = \beta_{0j} + X_i \beta_j + \alpha_i$ and $\alpha_i \sim \mathcal{N}(0,V_\alpha)$ |
| if n_latent>0 and site_effect="random"| log$(\theta_{ij}) = \beta_{0j} + X_i \beta_j + W_i \lambda_j + \alpha_i$ and $\alpha_i \sim \mathcal{N}(0,V_\alpha)$ |
jSDM_binomial_probit_long_format():
Ecological process:
$$y_{n} \sim \mathcal{B}ernoulli(\theta_n),$$ such as $species_n=j$ and $site_n=i$, where
| conditions | model specification |
|:---------------------------------------------|:-------------------------------------------------|
| if n_latent=0 and site_effect="none" | probit$(\theta_n) = D_n \gamma + X_n \beta_j$ |
| if n_latent>0 and site_effect="none" | probit$(\theta_n) = D_n \gamma + X_n \beta_j + W_i \lambda_j$ |
| if n_latent=0 and site_effect="random" | probit$(\theta_n) = D_n \gamma + X_n \beta_j + \alpha_i$ and $\alpha_i \sim \mathcal{N}(0,V_\alpha)$ |
| if n_latent>0 and site_effect="fixed" | probit$(\theta_n) = D_n \gamma + X_n \beta_j + W_i \lambda_j + \alpha_i$ |
| if n_latent=0 and site_effect="fixed" | probit$(\theta_n) = D_n \gamma + X_n \beta_j + \alpha_i$ |
| if n_latent>0 and site_effect="random" | probit$(\theta_n) = D_n \gamma + X_n \beta_j + W_i \lambda_j + \alpha_i$ and $\alpha_i \sim \mathcal{N}(0,V_\alpha)$ |
jSDM_gaussian():
Ecological process:
$$y_{ij} \sim \mathcal{N}(\theta_{ij}, V),$$ where
| conditions | model specification |
|:---------------------------------------------|:-------------------------------------------------|
| if n_latent=0 and site_effect="none" | $\theta_{ij} = \beta_{0j} + X_i \beta_j$ |
| if n_latent>0 and site_effect="none" | $\theta_{ij} = \beta_{0j} + X_i \beta_j + W_i \lambda_j$ |
| if n_latent=0 and site_effect="random" | $\theta_{ij} = \beta_{0j} + X_i \beta_j + \alpha_i$ and $\alpha_i \sim \mathcal{N}(0,V_\alpha)$ |
| if n_latent>0 and site_effect="fixed" | $\theta_{ij} = \beta_{0j} + X_i \beta_j + W_i \lambda_j + \alpha_i$ |
| if n_latent=0 and site_effect="fixed" | $\theta_{ij} = \beta_{0j} + X_i \beta_j + \alpha_i$ |
| if n_latent>0 and site_effect="random" | $\theta_{ij} = \beta_{0j} + X_i \beta_j + W_i \lambda_j + \alpha_i$ and $\alpha_i \sim \mathcal{N}(0,V_\alpha)$ |
Joint Species distribution models (jSDM) are useful tools to explain or predict species range and abundance from various environmental factors and species correlations [@Warton2015]. jSDM is becoming an increasingly popular statistical method in conservation biology.
In this vignette, we illustrate the use of the jSDM R package which aims at providing user-friendly statistical functions using field observations (occurrence or abundance data) to fit jSDMs models.
Package's functions are developed in a hierarchical Bayesian framework and use adaptive rejection Metropolis sampling algorithms or conjugate priors within Gibbs sampling to estimate model's parameters. Using compiled C++ code for the Gibbs sampler reduce drastically the computation time. By making these new statistical tools available to the scientific community, we hope to democratize the use of more complex, but more realistic, statistical models for increasing knowledge in ecology and conserving biodiversity.
(ref:cap-DAG) A graphical summary of the jSDM-package statistical framework. In this Directed Acyclic Graph (DAG), the orange boxes refer to data, the blue ellipses to parameters to be estimated, and the arrows to functional relationships described with the help of statistical distributions.
knitr::include_graphics("figures/DAG-jSDM-rand.png")
Model types available in jSDM R package are not limited to those described in this example. jSDM includes various model types for occurrence and abundance data, you can find more examples of use on the jSDM website.
Load librairies
We first load the jSDM library.
# Load libraries library(jSDM)
Bernoulli probit regression
Below, we show an example of the use of jSDM-package for fitting species distribution model to occurence data for 9 frog's species.
Definition of the model
Referring to the models used in the articles @Warton2015 and @Albert1993, we define the following model :
$$ \mathrm{probit}(\theta_{ij}) =\alpha_i + \beta_{0j}+X_i.\beta_j+ W_i.\lambda_j $$
-
Link function probit: $\mathrm{probit}: q \rightarrow \Phi^{-1}(q)$ where $\Phi$ correspond to the repartition function of the reduced centred normal distribution.
-
Response variable: $Y=(y_{ij})^{i=1,\ldots,nsite}_{j=1,\ldots,nsp}$ with:
$$y_{ij}=\begin{cases} 0 & \text{ if species $j$ is absent on the site $i$}\ 1 & \text{ if species $j$ is present on the site $i$}. \end{cases}$$
- Latent variable $z_{ij} = \alpha_i + \beta_{0j} + X_i.\beta_j + W_i.\lambda_j + \epsilon_{i,j}$, with $\forall (i,j) \ \epsilon_{ij} \sim \mathcal{N}(0,1)$ and such that:
$$y_{ij}=\begin{cases} 1 & \text{if} \ z_{ij} > 0 \ 0 & \text{otherwise.} \end{cases}$$
It can be easily shown that: $y_{ij} \sim \mathcal{B}ernoulli(\theta_{ij})$.
-
Latent variables: $W_i=(W_i^1,\ldots,W_i^q)$ where $q$ is the number of latent variables considered, which has to be fixed by the user (by default q=2). We assume that $W_i \sim \mathcal{N}(0,I_q)$ and we define the associated coefficients: $\lambda_j=(\lambda_j^1,\ldots, \lambda_j^q)'$. We use a prior distribution $\mathcal{N}(0,1)$ for all lambdas not concerned by constraints to $0$ on upper diagonal and to strictly positive values on diagonal.
-
Explanatory variables: bioclimatic data about each site. $X=(X_i)_{i=1,\ldots,nsite}$ with $X_i=(x_i^1,\ldots,x_i^p)\in \mathbb{R}^p$ where $p$ is the number of bioclimatic variables considered. The corresponding regression coefficients for each species $j$ are noted : $\beta_j=(\beta_j^1,\ldots,\beta_j^p)'$.
-
$\beta_{0j}$ correspond to the intercept for species $j$ which is assume to be a fixed effect. We use a prior distribution $\mathcal{N}(0,1)$ for all betas.
-
$\alpha_i$ represents the random effect of site $i$ such as $\alpha_i \sim \mathcal{N}(0,V_{\alpha})$ and we assumed that $V_{\alpha} \sim \mathcal {IG}(\text{shape}=0.1, \text{rate}=0.1)$ as prior distribution by default.
Occurrence data-set
(ref:cap-frog) Litoria ewingii [@Wilkinson2019].
knitr::include_graphics("figures/Litoria_ewingii.jpg")
This data-set is available in jSDM-package. It can be loaded with the data() command. The frogs dataset is in "wide" format: each line is a site and the occurrence data (from Species_1 to Species_9) are in columns. A site is characterized by its x-y geographical coordinates, one discrete covariate and two other continuous covariates.
# frogs data data(frogs, package="jSDM") head(frogs)
We rearrange the data in two data-sets: a first one for the presence-absence observations for each species (columns) at each site (rows), and a second one for the site characteristics.
We also normalize the continuous explanatory variables to facilitate MCMC convergence.
# data.obs PA_frogs <- frogs[,4:12] # Normalized continuous variables Env_frogs <- cbind(scale(frogs[,1]),frogs[,2],scale(frogs[,3])) colnames(Env_frogs) <- colnames(frogs[,1:3])
Parameter inference
We use the jSDM_binomial_probit() function to fit the jSDM (increase the number of iterations to achieve convergence).
mod_frogs_jSDM_probit <- jSDM_binomial_probit( # Chains burnin=1000, mcmc=1000, thin=1, # Response variable presence_data = PA_frogs, # Explanatory variables site_formula = ~., site_data = Env_frogs, # Model specification n_latent=2, site_effect="random", # Starting values alpha_start=0, beta_start=0, lambda_start=0, W_start=0, V_alpha=1, # Priors shape_Valpha=0.1, rate_Valpha=0.1, mu_beta=0, V_beta=1, mu_lambda=0, V_lambda=1, # Various seed=1234, verbose=1)
Analysis of the results
We visually evaluate the convergence of MCMCs by representing the trace and density a posteriori of some estimated parameters.
np <- nrow(mod_frogs_jSDM_probit$model_spec$beta_start) oldpar <- par(no.readonly = TRUE) ## beta_j of the first two species par(mfrow=c(2,2)) for (j in 1:2) { for (p in 1:np) { coda::traceplot(coda::as.mcmc(mod_frogs_jSDM_probit$mcmc.sp[[j]][,p])) coda::densplot(coda::as.mcmc(mod_frogs_jSDM_probit$mcmc.sp[[j]][,p]), main = paste(colnames(mod_frogs_jSDM_probit$mcmc.sp[[j]])[p], ", species : ",j), cex.main=0.9) } } ## lambda_j of the first two species n_latent <- mod_frogs_jSDM_probit$model_spec$n_latent par(mfrow=c(2,2)) for (j in 1:2) { for (l in 1:n_latent) { coda::traceplot(coda::as.mcmc(mod_frogs_jSDM_probit$mcmc.sp[[j]][,np+l])) coda::densplot(coda::as.mcmc(mod_frogs_jSDM_probit$mcmc.sp[[j]][,np+l]), main = paste(colnames(mod_frogs_jSDM_probit$mcmc.sp[[j]]) [np+l],", species : ",j), cex.main=0.9) } } ## Latent variables W_i for the first two sites par(mfrow=c(2,2)) for (l in 1:n_latent) { for (i in 1:2) { coda::traceplot(mod_frogs_jSDM_probit$mcmc.latent[[paste0("lv_",l)]][,i], main = paste0("Latent variable W_", l, ", site ", i), cex.main=0.9) coda::densplot(mod_frogs_jSDM_probit$mcmc.latent[[paste0("lv_",l)]][,i], main = paste0("Latent variable W_", l, ", site ", i), cex.main=0.9) } } ## alpha_i of the first two sites plot(coda::as.mcmc(mod_frogs_jSDM_probit$mcmc.alpha[,1:2])) ## V_alpha par(mfrow=c(2,2)) coda::traceplot(mod_frogs_jSDM_probit$mcmc.V_alpha) coda::densplot(mod_frogs_jSDM_probit$mcmc.V_alpha) ## Deviance coda::traceplot(mod_frogs_jSDM_probit$mcmc.Deviance) coda::densplot(mod_frogs_jSDM_probit$mcmc.Deviance) ## probit_theta par (mfrow=c(1,2)) hist(mod_frogs_jSDM_probit$probit_theta_latent, main = "Predicted probit theta", xlab ="predicted probit theta") hist(mod_frogs_jSDM_probit$theta_latent, main = "Predicted theta", xlab ="predicted theta") par(oldpar)
Overall, the traces and the densities of the parameters indicate the convergence of the algorithm. Indeed, we observe on the traces that the values oscillate around averages without showing an upward or downward trend and we see that the densities are quite smooth and for the most part of Gaussian form.
Matrice of correlations
After fitting the jSDM with latent variables, the full species residual correlation matrix $R=(R_{ij})^{i=1,\ldots, n_{species}}{j=1,\ldots, n{species}}$ can be derived from the covariance in the latent variables such as : $$\Sigma_{ij} = \lambda_i^T .\lambda_j $$, then we compute correlations from covariances : $$R_{i,j} = \frac{\Sigma_{ij}}{\sqrt{\Sigma {ii}\Sigma {jj}}}$$.
We use the function plot_residual_cor() to compute and display the residual correlation matrix between species :
plot_residual_cor(mod_frogs_jSDM_probit)
Predictions
We use the predict.jSDM() S3 method on the mod_frogs_jSDM_probit object of class jSDM to compute the mean (or expectation) of the posterior distributions obtained and get the expected values of model's parameters.
# Sites and species concerned by predictions : ## 50 sites among the 104 Id_sites <- sample.int(nrow(PA_frogs), 50) ## All species Id_species <- colnames(PA_frogs) # Simulate new observations of covariates on those sites simdata <- matrix(nrow=50, ncol = ncol(mod_frogs_jSDM_probit$model_spec$site_data)) colnames(simdata) <- colnames(mod_frogs_jSDM_probit$model_spec$site_data) rownames(simdata) <- Id_sites simdata <- as.data.frame(simdata) simdata$Covariate_1 <- rnorm(50) simdata$Covariate_3 <- rnorm(50) simdata$Covariate_2 <- rbinom(50,1,0.5) # Predictions theta_pred <- predict(mod_frogs_jSDM_probit, newdata=simdata, Id_species=Id_species, Id_sites=Id_sites, type="mean") hist(theta_pred, main="Predicted theta with simulated data", xlab="predicted theta")
References
Try the jSDM package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.