---
title: "Get started with jSDM"
output:
  bookdown::html_document2:
    #base_format: rmarkdown::html_vignette
    number_sections: true
    toc: true
    toc_float: true
    fig_caption: yes
link-citations: yes
bibliography: bib/biblio-jSDM.bib
biblio-style: bib/jae.bst
csl: bib/journal-of-applied-ecology.csl
pkgdown:
  as_is: true
vignette: >
 %\VignetteEngine{knitr::rmarkdown}
 %\VignetteIndexEntry{Get started with jSDM}
 %\VignetteEncoding{UTF-8}
---
 
```{r setup, include=FALSE}
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. 

<center> __Directed Acyclic Graph (DAG)__ </center>

(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.

```{r plot-DAG, echo=FALSE, out.width=500, fig.retina=1, fig.cap="(ref:cap-DAG)", eval=TRUE}
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](https://ecology.ghislainv.fr/jSDM/). 

# Load librairies 

We first load the `jSDM` library.

```{r libraries}
# 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].

```{r frog-picture, echo=FALSE, out.width=400, out.height=300, fig.retina=1, fig.cap="(ref:cap-frog)"}
knitr::include_graphics("figures/Litoria_ewingii.jpg")
```

This data-set is available in [`jSDM-package`](https://ecology.ghislainv.fr/jSDM/reference/jSDM-package.html). It can be loaded with the `data()` command. The [`frogs`](https://ecology.ghislainv.fr/jSDM/reference/frogs.html) 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. 

```{r frogs-data}
# 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. 

```{r arranging-frogs-data}
# 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). 

```{r jSDM-probit}
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.

```{r plot-results-probit}
np <- nrow(mod_frogs_jSDM_probit$model_spec$beta_start)

## 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")
```

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 :

```{r correlation-matrix-probit}
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.

```{r predictions-probit}
# 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