---
title: "Get started with jSDM"
output:
  bookdown::html_document2:
    #base_format: rmarkdown::html_vignette
    #highlight: tango
    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.width = 6, fig.height = 6,
	cache = FALSE,
	collapse = TRUE,
	comment = "#>",
	highlight = TRUE
)
```

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 wich aims at providing user-friendly statistical functions using field observations (occurence or abundance data) to fit jSDMs models. Package's functions are developped in a hierarchical Bayesian framework and use conjugate priors 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. 

Below, we show an example of the use of `jSDM` for fitting species distribution model to abundance data for 9 frog's species. Model types available in `jSDM` are not limited to those described in this example. `jSDM` includes various model types for occurrence data: 

# 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}ernouilli(\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,10)$ for all lambdas not concerned by constraints. 

- Explanatory variables: bioclimatic data about each site. $X=(X_i)_{i=1,\ldots,nsite}$ avec $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,10^6)$ 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.5, \text{rate}=0.005)$ as prior distribution by default. 


# Data-set

(ref:cap-frog) **_Litoria ewingii_** [@Wilkinson2019].

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

We first load the `jSDM` library.

```{r libraries}
# Load libraries
library(jSDM)
```

This data-set is available in the `jSDM` R package. It can be loaded with the `data` command. The data 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 explicative variables to facilitate MCMC convergence. 

```{r arranging-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_probit_block()` function.

```{r jSDM, cache=TRUE}
mod_jSDM_block_frogs <- jSDM_probit_block (
  # Response variable 
  presence_site_sp = as.matrix(PA_frogs), 
  # Explanatory variables 
  site_suitability = ~.,   
  site_data = as.data.frame(Env_frogs), n_latent=2,
  # Chains
  burnin=20000, mcmc=5000, thin=5,
  # Starting values
  alpha_start=0, beta_start=0,
  lambda_start=0, W_start=0,
  V_alpha_start=1, 
  # Priors
  shape=0.5, rate=0.0005,
  mu_beta=0, V_beta=1.0E6,
  mu_lambda=0, V_lambda=10,
  # Various 
  seed=1234, verbose=1)
```

# Analysis of the results

```{r plot-results}
## alpha_i of the first two sites
plot(coda::as.mcmc(mod_jSDM_block_frogs$mcmc.alpha[,1:2]))

## Valpha
par(mfrow=c(1,2))
coda::traceplot(mod_jSDM_block_frogs$mcmc.Valpha, main="V_alpha")
coda::densplot(mod_jSDM_block_frogs$mcmc.Valpha, main="V_alpha")

np <- nrow(mod_jSDM_block_frogs$model_spec$beta_start)

## beta_j of the first two species
par(mfrow=c(np,2))
for (j in 1:2) {
  for (p in 1:np) {
      coda::traceplot(coda::as.mcmc(mod_jSDM_block_frogs$mcmc.sp[[paste0("sp_",j)]][,p]))
      coda::densplot(coda::as.mcmc(mod_jSDM_block_frogs$mcmc.sp[[paste0("sp_",j)]][,p]), 
main = paste(colnames(mod_jSDM_block_frogs$mcmc.sp[[paste0("sp_",j)]])[p],
", species : ",j))
  }
}

## lambda_j of the first two species
n_latent <- mod_jSDM_block_frogs$model_spec$n_latent
par(mfrow=c(n_latent*2,2))
for (j in 1:2) {
  for (l in 1:n_latent) {
      coda::traceplot(coda::as.mcmc(mod_jSDM_block_frogs$mcmc.sp[[paste0("sp_",j)]][,np+l]))
      coda::densplot(coda::as.mcmc(mod_jSDM_block_frogs$mcmc.sp[[paste0("sp_",j)]][,np+l]), 
      main = paste(colnames(mod_jSDM_block_frogs$mcmc.sp[[paste0("sp_",j)]])[np+l],
      ", species : ",j))
    }
  }

## Latent variables W_i for the first two sites
par(mfrow=c(1,2))
for (l in 1:n_latent) {
  plot(mod_jSDM_block_frogs$mcmc.latent[[paste0("lv_",l)]][,1:2], main = paste0("Latent variable W_", l))
}

## probit_theta
par (mfrow=c(1,1))
hist(mod_jSDM_block_frogs$probit_theta_pred, main = "Predicted probit theta", xlab ="predicted probit theta")

## Deviance
plot(mod_jSDM_block_frogs$mcmc.Deviance,main = "Deviance")

```

# Matrice of correlations 

We use the `plot_residual_cor()` function.

```{r correlation-matrix}
plot_residual_cor(mod_jSDM_block_frogs)
```

# Predictions 

We use the `predict_jSDM()` function.

```{r predictions}
# 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_jSDM_block_frogs$model_spec$site_data))
colnames(simdata) <- colnames(mod_jSDM_block_frogs$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.jSDM(mod_jSDM_block_frogs, 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
