Package 'jSDM'

joint species distribution models

Description

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

Details

The package includes the following functions to fit various species distribution models :

function	data-type
jSDM_binomial_logit	presence-absence
jSDM_binomial_probit	presence-absence
jSDM_binomial_probit_sp_constrained	presence-absence
jSDM_binomial_probit_long_format	presence-absence
jSDM_poisson_log	abundance

• jSDM_binomial_probit :
  
  Ecological process:

  $y_{ij} \sim \mathcal{B}ernoulli(\theta_{ij})$

  where

  if n_latent=0 and site_effect="none"	probit$(\theta_{ij}) = X_i \beta_j$
  if n_latent>0 and site_effect="none"	probit$(\theta_{ij}) = X_i \beta_j + W_i \lambda_j$
  if n_latent=0 and site_effect="fixed"	probit$(\theta_{ij}) = 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}) = X_i \beta_j + W_i \lambda_j + \alpha_i$
  if n_latent=0 and site_effect="random"	probit$(\theta_{ij}) = X_i \beta_j + \alpha_i$
  if n_latent>0 and site_effect="random"	probit$(\theta_{ij}) = 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 ($\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

  if n_latent=0 and site_effect="none"	logit$(\theta_{ij}) = X_i \beta_j$
  if n_latent>0 and site_effect="none"	logit$(\theta_{ij}) = X_i \beta_j + W_i \lambda_j$
  if n_latent=0 and site_effect="fixed"	logit$(\theta_{ij}) = X_i \beta_j + \alpha_i$
  if n_latent>0 and site_effect="fixed"	logit$(\theta_{ij}) = X_i \beta_j + W_i \lambda_j + \alpha_i$
  if n_latent=0 and site_effect="random"	logit$(\theta_{ij}) = 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}) = 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

  if n_latent=0 and site_effect="none"	log$(\theta_{ij}) = X_i \beta_j$
  if n_latent>0 and site_effect="none"	log$(\theta_{ij}) = X_i \beta_j + W_i \lambda_j$
  if n_latent=0 and site_effect="fixed"	log$(\theta_{ij}) = X_i \beta_j + \alpha_i$
  if n_latent>0 and site_effect="fixed"	log$(\theta_{ij}) = X_i \beta_j + W_i \lambda_j + \alpha_i$
  if n_latent=0 and site_effect="random"	log$(\theta_{ij}) = 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}) = 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

  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="fixed"	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="random"	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)$

Author(s)

Ghislain Vieilledent <[email protected]>

Jeanne Clément <[email protected]>

Frédéric Gosselin <[email protected]>

References

Chib, S. and Greenberg, E. (1998) Analysis of multivariate probit models. Biometrika, 85, 347-361.

Warton, D. I.; Blanchet, F. G.; O'Hara, R. B.; O'Hara, R. B.; Ovaskainen, O.; Taskinen, S.; Walker, S. C. and Hui, F. K. C. (2015) So Many Variables: Joint Modeling in Community Ecology. Trends in Ecology & Evolution, 30, 766-779.

Ovaskainen, O., Tikhonov, G., Norberg, A., Blanchet, F. G., 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, 561-576.

---

Distribution of Alpine plants in Aravo (Valloire, France)

Description

This dataset describe the distribution of 82 species of Alpine plants in 75 sites. Species traits and environmental variables are also measured.

Usage

data("aravo")

Format

aravo is a list containing the following objects :

spe

is a data.frame with the abundance values of 82 species (columns) in 75 sites (rows).

env

is a data.frame with the measurements of 6 environmental variables for the sites.

traits

is data.frame with the measurements of 8 traits for the species.

spe.names

is a vector with full species names.

Details

The environmental variables are :

Aspect	Relative south aspect (opposite of the sine of aspect with flat coded 0)
Slope	Slope inclination (degrees)
Form	Microtopographic landform index: 1 (convexity); 2 (convex slope); 3 (right slope);
	4 (concave slope); 5 (concavity)
Snow	Mean snowmelt date (Julian day) averaged over 1997-1999
PhysD	Physical disturbance, i.e., percentage of unvegetated soil due to physical processes
ZoogD	Zoogenic disturbance, i.e., quantity of unvegetated soil due to marmot activity: no; some; high

The species traits for the plants are:

Height	Vegetative height (cm)
Spread	Maximum lateral spread of clonal plants (cm)
Angle	Leaf elevation angle estimated at the middle of the lamina
Area	Area of a single leaf
Thick	Maximum thickness of a leaf cross section (avoiding the midrib)
SLA	Specific leaf area
Nmass	Mass-based leaf nitrogen content
Seed	Seed mass

Source

Choler, P. (2005) Consistent shifts in Alpine plant traits along a mesotopographical gradient. Arctic, Antarctic, and Alpine Research 37,444-453.

Dray S, Dufour A (2007). The ade4 Package: Implementing the Duality Diagram for Ecologists. Journal of Statistical Software, 22(4), 1-20. doi:10.18637/jss.v022.i04.

Examples

data(aravo, package="jSDM")
summary(aravo)

---

birds dataset

Description

The Swiss breeding bird survey ("Monitoring Häufige Brutvögel" MHB) has monitored the populations of 158 common species since 1999.

The MHB sample from data(MHB2014, package="AHMbook") consists of 267 1-km squares that are laid out as a grid across Switzerland. Fieldwork is conducted by about 200 skilled birdwatchers, most of them volunteers. Avian populations are monitored using a simplified territory mapping protocol, where each square is surveyed up to three times during the breeding season (only twice above the tree line).

Surveys are conducted along a transect that does not change over the years. The birds dataset has the data for 2014, except one quadrat not surveyed in 2014.

Usage

data("birds")

Format

A data frame with 266 observations on the following 166 variables.

158 bird species named in latin and whose occurrences are indicated as the number of visits to each site during which the species was observed, including 13 species not recorded in the year 2014 and 5 covariates collected on the 266 1x1 km quadrat as well as their identifiers and coordinates :

siteID	an alphanumeric site identifier
coordx	a numeric vector indicating the x coordinate of the centre of the quadrat.
	The coordinate reference system is not specified intentionally.
coordy	a numeric vector indicating the y coordinate of the centre of the quadrat.
elev	a numeric vector indicating the mean elevation of the quadrat (m).
rlength	the length of the route walked in the quadrat (km).
nsurvey	a numeric vector indicating the number of replicate surveys planned in the quadrat;
	above the tree-line 2, otherwise 3.
forest	a numeric vector indicating the percentage of forest cover in the quadrat.
obs14	a categorical vector indicating the identifying number of the observer.

Details

Only the Latin names of bird species are given in this dataset but you can find the corresponding English names in the original dataset : data(MHB2014, package="AHMbook").

Source

Swiss Ornithological Institute

References

Kéry and Royle (2016) Applied Hierarachical Modeling in Ecology Section 11.3

Examples

data(birds, package="jSDM")
head(birds)
# find species not recorded in 2014
which(colSums(birds[,1:158])==0)

---

eucalypts dataset

Description

The Eucalyptus data set includes 12 taxa recorded in 458 plots spanning elevation gradients in the Grampians National Park, Victoria, which is known for high species diversity and endemism. The park has three mountain ranges interspersed with alluvial valleys and sand sheet and has a semi-Mediterranean climate with warm, dry summers and cool, wet winters.

This dataset records presence or absence at 458 sites of 12 eucalypts species, 7 covariates collected at these sites as well as their longitude and latitude.

Usage

data("eucalypts")

Format

A data frame with 458 observations on the following 21 variables.

12 eucalypts species which presence on sites is indicated by a 1 and absence by a 0 :

ALA

a binary vector indicating the occurrence of the species E. alaticaulis

ARE

a binary vector indicating the occurrence of the species E. arenacea

BAX

a binary vector indicating the occurrence of the species E. baxteri

CAM

a binary vector indicating the occurrence of the species E. camaldulensis

GON

a binary vector indicating the occurrence of the species E. goniocalyx

MEL

a binary vector indicating the occurrence of the species E. melliodora

OBL

a binary vector indicating the occurrence of the species E. oblique

OVA

a binary vector indicating the occurrence of the species E. ovata

WIL

a binary vector indicating the occurrence of the species E. willisii subsp. Falciformis

ALP

a binary vector indicating the occurrence of the species E. serraensis, E. verrucata and E. victoriana

VIM

a binary vector indicating the occurrence of the species E. viminalis subsp. Viminalis and Cygnetensis

ARO.SAB

a binary vector indicating the occurrence of the species E. aromaphloia and E. sabulosa

7 covariates collected on the 458 sites and their coordinates :

Rockiness

a numeric vector taking values from 0 to 95 corresponding to the rock cover of the site in percent estimated in 5 % increments in field plots

Sandiness

a binary vector indicating if soil texture categorie is sandiness based on soil texture classes from field plots and according to relative amounts of sand, silt, and clay particles

VallyBotFlat

a numeric vector taking values from 0 to 6 corresponding to the valley bottom flatness GIS-derived variable defining flat areas relative to surroundings likely to accumulate sediment (units correspond to the percentage of slope e.g. 0.5 = 16 %slope, 4.5 = 1 %slope, 5.5 = 0.5 %slope)

PPTann

a numeric vector taking values from 555 to 1348 corresponding to annual precipitation in millimeters measured as the sum of monthly precipitation estimated using BIOCLIM based on 20m grid cell Digital Elevation Model

Loaminess

a binary vector indicating if soil texture categorie is loaminess based on soil texture classes from field plots and according to relative amounts of sand, silt, and clay particles

cvTemp

a numeric vector taking values from 136 to 158 corresponding to coefficient of variation of temperature seasonality in percent measured as the standard deviation of weekly mean temperatures as a percentage of the annual mean temperature from BIOCLIM

T0

a numeric vector corresponding to solar radiation in $WH/m^2$ measured as the amount of incident solar energy based on the visible sky and the sun's position. Derived from Digital Elevation Model in ArcGIS 9.2 Spatial Analyst for the summer solstice (December 22)

latitude

a numeric vector indicating the latitude of the studied site

longitude

a numeric vector indicating the longitude of the studied site

Source

Wilkinson, D. P.; Golding, N.; Guillera-Arroita, G.; Tingley, R. and McCarthy, M. A. (2018) A comparison of joint species distribution models for presence-absence data. Methods in Ecology and Evolution.

Examples

data(eucalypts, package="jSDM")
head(eucalypts)

---

frogs dataset

Description

Presence or absence of 9 species of frogs at 104 sites, 3 covariates collected at these sites and their coordinates.

Format

frogs is a data frame with 104 observations on the following 14 variables.

Species_

1 to 9 indicate by a 0 the absence of the species on one site and by a 1 its presence

Covariates_

1 and 3 continuous variables

Covariates_

2 discrete variables

x

a numeric vector of first coordinates corresponding to each site

y

a numeric vector of second coordinates corresponding to each site

Source

Wilkinson, D. P.; Golding, N.; Guillera-Arroita, G.; Tingley, R. and McCarthy, M. A. (2018) A comparison of joint species distribution models for presence-absence data. Methods in Ecology and Evolution.

Examples

data(frogs, package="jSDM")
head(frogs)

---

fungi dataset

Description

Presence or absence of 11 species of fungi on dead-wood objects at 800 sites and 12 covariates collected at these sites.

Usage

data("fungi")

Format

A data frame with 800 observations on the following 23 variables :

11 fungi species which presence on sites is indicated by a 1 and absence by a 0 :

antser

a binary vector

antsin

a binary vector

astfer

a binary vector

fompin

a binary vector

hetpar

a binary vector

junlut

a binary vector

phefer

a binary vector

phenig

a binary vector

phevit

a binary vector

poscae

a binary vector

triabi

a binary vector

12 covariates collected on the 800 sites :

diam

a numeric vector indicating the diameter of dead-wood object

dc1

a binary vector indicating if the decay class is 1 measured in the scale 1, 2, 3, 4, 5 (from freshly decayed to almost completely decayed)

dc2

a binary vector indicating if the decay class is 2

dc3

a binary vector indicating if the decay class is 3

dc4

a binary vector indicating if the decay class is 4

dc5

a binary vector indicating if the decay class is 5

quality3

a binary vector indicating if the quality is level 3

quality4

a binary vector indicating if the quality is level 4

ground3

a binary vector indicating if the ground contact is level 3 as 2 = no ground contact, 3 = less than half of the log in ground contact and 4 = more than half of the log in ground contact

ground4

a binary vector a binary vector indicating if the ground contact is level 4

epi

a numeric vector indicating the epiphyte cover

bark

a numeric vector indicating the bark cover

Source

Wilkinson, D. P.; Golding, N.; Guillera-Arroita, G.; Tingley, R. and McCarthy, M. A. (2018) A comparison of joint species distribution models for presence-absence data. Methods in Ecology and Evolution.

Examples

data(fungi, package="jSDM")
head(fungi)

---

Extract covariances and correlations due to shared environmental responses

Description

Calculates the correlation between columns of the response matrix, due to similarities in the response to explanatory variables i.e., shared environmental response.

Usage

get_enviro_cor(mod, type = "mean", prob = 0.95)

Arguments

mod	An object of class "jSDM"
type	A choice of either the posterior median (type = "median") or posterior mean (type = "mean"), which are then treated as estimates and the fitted values are calculated from. Default is posterior mean.
prob	A numeric scalar in the interval $(0,1)$ giving the target probability coverage of the intervals, by which to determine whether the correlations are "significant". Defaults to 0.95.

Details

In both independent response and correlated response models, where each of the columns of the response matrix $Y$ are fitted to a set of explanatory variables given by $X$, the covariance between two columns $j$ and $j'$, due to similarities in their response to the model matrix, is thus calculated based on the linear predictors $X \beta_j$ and $X \beta_j'$, where $\beta_j$ are species effects relating to the explanatory variables. Such correlation matrices are discussed and found in Ovaskainen et al., (2010), Pollock et al., (2014).

Value

results, a list including :

cor, cor.lower, cor.upper	A set of $np \times np$ correlation matrices, containing either the posterior median or mean estimate over the MCMC samples plus lower and upper limits of the corresponding 95 % highest posterior interval.
cor.sig	A $np \times np$ correlation matrix containing only the “significant" correlations whose 95 % highest posterior density (HPD) interval does not contain zero. All non-significant correlations are set to zero.
cov	Average over the MCMC samples of the $np \times np$ covariance matrix.

Author(s)

Ghislain Vieilledent <[email protected]>

Jeanne Clément <[email protected]>

References

Hui FKC (2016). “boral: Bayesian Ordination and Regression Analysis of Multivariate Abundance Data in R.” Methods in Ecology and Evolution, 7, 744–750.

Ovaskainen et al. (2010). Modeling species co-occurrence by multivariate logistic regression generates new hypotheses on fungal interactions. Ecology, 91, 2514-2521.

Pollock et al. (2014). Understanding co-occurrence by modelling species simultaneously with a Joint Species Distribution Model (JSDM). Methods in Ecology and Evolution, 5, 397-406.

See Also

cov2cor get_residual_cor jSDM-package jSDM_binomial_probit
jSDM_binomial_logit jSDM_poisson_log

Examples

library(jSDM)
# frogs data
 data(frogs, package="jSDM")
 # Arranging data
 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])
 Env_frogs <- as.data.frame(Env_frogs)
 # Parameter inference
# Increase the number of iterations to reach MCMC convergence
mod <- jSDM_binomial_probit(# Response variable
                             presence_data=PA_frogs,
                             # Explanatory variables
                             site_formula = ~.,
                             site_data = Env_frogs,
                             n_latent=0,
                             site_effect="random",
                             # Chains
                             burnin=100,
                             mcmc=100,
                             thin=1,
                             # Starting values
                             alpha_start=0,
                             beta_start=0,
                             V_alpha=1,
                             # Priors
                             shape=0.5, rate=0.0005,
                             mu_beta=0, V_beta=10,
                             # Various
                             seed=1234, verbose=1)
# Calcul of residual correlation between species 
 enviro.cors <- get_enviro_cor(mod)

---

Calculate the residual correlation matrix from a latent variable model (LVM)

Description

This function use coefficients $(\lambda_{jl}$ with $j=1,\dots,n_{species}$ and $l=1,\dots,n_{latent})$, corresponding to latent variables fitted using jSDM package, to calculate the variance-covariance matrix which controls correlation between species.

Usage

get_residual_cor(mod, prob = 0.95, type = "mean")

Arguments

mod	An object of class "jSDM"
prob	A numeric scalar in the interval $(0,1)$ giving the target probability coverage of the highest posterior density (HPD) intervals, by which to determine whether the correlations are "significant". Defaults to 0.95.
type	A choice of either the posterior median (type = "median") or posterior mean (type = "mean"), which are then treated as estimates and the fitted values are calculated from. Default is posterior mean.

Details

After fitting the jSDM with latent variables, the fullspecies residual correlation matrix : $R=(R_{ij})$ with $i=1,\ldots, n_{species}$ and $j=1,\ldots, n_{species}$ can be derived from the covariance in the latent variables such as : $\Sigma_{ij}=\lambda_i' .\lambda_j$, in the case of a regression with probit, logit or poisson link function and

$\Sigma_{ij}$	$= \lambda_i' .\lambda_j + V$	if i=j
	$= \lambda_i' .\lambda_j$	else,

, in the case of a linear regression with a response variable such as

$y_{ij} \sim \mathcal{N}(\theta_{ij}, V)$

. Then we compute correlations from covariances :

$R_{ij} = \frac{\Sigma_{ij}}{\sqrt{\Sigma_ii\Sigma _jj}}$

.

Value

results A list including :

cov.mean	Average over the MCMC samples of the variance-covariance matrix, if type = "mean".
cov.median	Median over the MCMC samples of the variance-covariance matrix, if type = "median".
cov.lower	A $n_{species} \times n_{species}$ matrix containing the lower limits of the ($100 \times prob \%$) HPD interval of variance-covariance matrices over the MCMC samples.
cov.upper	A $n_{species} \times n_{species}$ matrix containing the upper limits of the ($100 \times prob \%$) HPD interval of variance-covariance matrices over the MCMC samples.
cov.sig	A $n_{species} \times n_{species}$ matrix containing the value 1 corresponding to the “significant" co-variances and the value 0 corresponding to "non-significant" co-variances, whose ($100 \times prob \%$) HPD interval over the MCMC samples contain zero.
cor.mean	Average over the MCMC samples of the residual correlation matrix, if type = "mean".
cor.median	Median over the MCMC samples of the residual correlation matrix, if type = "median".
cor.lower	A $n_{species} \times n_{species}$ matrix containing the lower limits of the ($100 \times prob \%$) HPD interval of correlation matrices over the MCMC samples.
cor.upper	A $n_{species} \times n_{species}$ matrix containing the upper limits of the ($100 \times prob \%$) HPD interval of correlation matrices over the MCMC samples.
cor.sig	A $n_{species} \times n_{species}$ matrix containing the value $1$ corresponding to the “significant" correlations and the value $0$ corresponding to "non-significant" correlations, whose ($100 \times prob \%$) HPD interval over the MCMC samples contain zero.

Author(s)

Ghislain Vieilledent <[email protected]>

Jeanne Clément <[email protected]>

References

Hui FKC (2016). boral: Bayesian Ordination and Regression Analysis of Multivariate Abundance Data in R. Methods in Ecology and Evolution, 7, 744–750.

Ovaskainen and al. (2016). Using latent variable models to identify large networks of species-to-species associations at different spatial scales. Methods in Ecology and Evolution, 7, 549-555.

Pollock and al. (2014). Understanding co-occurrence by modelling species simultaneously with a Joint Species Distribution Model (JSDM). Methods in Ecology and Evolution, 5, 397-406.

See Also

get_enviro_cor cov2cor jSDM-package jSDM_binomial_probit
jSDM_binomial_logit jSDM_poisson_log

Examples

library(jSDM)
# frogs data
 data(frogs, package="jSDM")
 # Arranging data
 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])
 Env_frogs <- as.data.frame(Env_frogs)
 # Parameter inference
# Increase the number of iterations to reach MCMC convergence
 mod <- jSDM_binomial_probit(# Response variable
                             presence_data=PA_frogs,
                             # Explanatory variables
                             site_formula = ~.,
                             site_data = Env_frogs,
                             n_latent=2,
                             site_effect="random",
                             # Chains
                             burnin=100,
                             mcmc=100,
                             thin=1,
                             # Starting values
                             alpha_start=0,
                             beta_start=0,
                             lambda_start=0,
                             W_start=0,
                             V_alpha=1,
                             # Priors
                             shape=0.5, rate=0.0005,
                             mu_beta=0, V_beta=10,
                             mu_lambda=0, V_lambda=10,
                             # Various
                             seed=1234, verbose=1)
# Calcul of residual correlation between species 
result <- get_residual_cor(mod, prob=0.95, type="mean")
# Residual variance-covariance matrix
result$cov.mean
## All non-significant co-variances are set to zero.
result$cov.mean * result$cov.sig
# Residual correlation matrix
result$cor.mean
## All non-significant correlations are set to zero.
result$cor.mean * result$cor.sig

---

Generalized inverse logit function

Description

Compute generalized inverse logit function.

Usage

inv_logit(x, min = 0, max = 1)

Arguments

x	value(s) to be transformed
min	Lower end of logit interval
max	Upper end of logit interval

Details

The generalized inverse logit function takes values on [-Inf,Inf] and transforms them to span [min, max] :

$y = p' (max-min) + min$

$where$

$p =\frac{exp(x)}{(1+exp(x))}$

Value

y Transformed value(s).

Author(s)

Gregory R. Warnes <[email protected]>

Examples

x <- seq(0,10, by=0.25)
xt <- jSDM::logit(x, min=0, max=10)
cbind(x,xt)
y <- jSDM::inv_logit(xt, min=0, max=10)
cbind(x,xt,y)

---

Binomial logistic regression

Description

The jSDM_binomial_logit function performs a Binomial logistic regression in a Bayesian framework. The function calls a Gibbs sampler written in 'C++' code which uses an adaptive Metropolis algorithm to estimate the conditional posterior distribution of model's parameters.

Usage

jSDM_binomial_logit(
  burnin = 5000,
  mcmc = 10000,
  thin = 10,
  presence_data,
  site_formula,
  trait_data = NULL,
  trait_formula = NULL,
  site_data,
  trials = NULL,
  n_latent = 0,
  site_effect = "none",
  beta_start = 0,
  gamma_start = 0,
  lambda_start = 0,
  W_start = 0,
  alpha_start = 0,
  V_alpha = 1,
  shape_Valpha = 0.5,
  rate_Valpha = 5e-04,
  mu_beta = 0,
  V_beta = 10,
  mu_gamma = 0,
  V_gamma = 10,
  mu_lambda = 0,
  V_lambda = 10,
  ropt = 0.44,
  seed = 1234,
  verbose = 1
)

Arguments

burnin	The number of burnin iterations for the sampler.
mcmc	The number of Gibbs iterations for the sampler. Total number of Gibbs iterations is equal to burnin+mcmc. burnin+mcmc must be divisible by 10 and superior or equal to 100 so that the progress bar can be displayed.
thin	The thinning interval used in the simulation. The number of mcmc iterations must be divisible by this value.
presence_data	A matrix $n_{site} \times n_{species}$ indicating the number of successes (or presences) and the absence by a zero for each species at the studied sites.
site_formula	A one-sided formula of the form '~x1+...+xp' specifying the explanatory variables for the suitability process of the model, used to form the design matrix $X$ of size $n_{site} \times np$.
trait_data	A data frame containing the species traits which can be included as part of the model. Default to NULL to fit a model without species traits.
trait_formula	A one-sided formula of the form '~ t1 + ... + tk + x1:t1 + ... + xp:tk' specifying the interactions between the environmental variables and the species traits to be considered in the model, used to form the trait design matrix $Tr$ of size $n_{species} \times nt$ and to set to $0$ the $\gamma$ parameters corresponding to interactions not taken into account according to the formula. Default to NULL to fit a model with all possible interactions between species traits found in trait_data and environmental variables defined by site_formula.
site_data	A data frame containing the model's explanatory variables by site.
trials	A vector indicating the number of trials for each site. $n_i$ should be superior or equal to $y_{ij}$, the number of successes for observation $n$. If $n_i=0$, then $y_{ij}=0$. The default is one visit by site.
n_latent	An integer which specifies the number of latent variables to generate. Defaults to $0$.
site_effect	A string indicating whether row effects are included as fixed effects ("fixed"), as random effects ("random"), or not included ("none") in the model. If fixed effects, then for parameter identifiability the first row effect is set to zero, which analogous to acting as a reference level when dummy variables are used. If random effects, they are drawn from a normal distribution with mean zero and unknown variance, analogous to a random intercept in mixed models. Defaults to "none".
beta_start	Starting values for $\beta$ parameters of the suitability process for each species must be either a scalar or a $np \times n_{species}$ matrix. If beta_start takes a scalar value, then that value will serve for all of the $\beta$ parameters.
gamma_start	Starting values for $\gamma$ parameters that represent the influence of species-specific traits on species' responses $\beta$, gamma_start must be either a scalar, a vector of length $nt$, $np$ or $nt.np$ or a $nt \times np$ matrix. If gamma_start takes a scalar value, then that value will serve for all of the $\gamma$ parameters. If gamma_start is a vector of length $nt$ or $nt.np$ the resulting $nt \times np$ matrix is filled by column with specified values, if a $np$-length vector is given, the matrix is filled by row.
lambda_start	Starting values for $\lambda$ parameters corresponding to the latent variables for each species must be either a scalar or a $n_{latent} \times n_{species}$ upper triangular matrix with strictly positive values on the diagonal, ignored if n_latent=0. If lambda_start takes a scalar value, then that value will serve for all of the $\lambda$ parameters except those concerned by the constraints explained above.
W_start	Starting values for latent variables must be either a scalar or a $nsite \times n_latent$ matrix, ignored if n_latent=0. If W_start takes a scalar value, then that value will serve for all of the $W_{il}$ with $i=1,\ldots,n_{site}$ and $l=1,\ldots,n_{latent}$.
alpha_start	Starting values for random site effect parameters must be either a scalar or a $n_{site}$-length vector, ignored if site_effect="none". If alpha_start takes a scalar value, then that value will serve for all of the $\alpha$ parameters.
V_alpha	Starting value for variance of random site effect if site_effect="random" or constant variance of the Gaussian prior distribution for the fixed site effect if site_effect="fixed". Must be a strictly positive scalar, ignored if site_effect="none".
shape_Valpha	Shape parameter of the Inverse-Gamma prior for the random site effect variance V_alpha, ignored if site_effect="none" or site_effect="fixed". Must be a strictly positive scalar. Default to 0.5 for weak informative prior.
rate_Valpha	Rate parameter of the Inverse-Gamma prior for the random site effect variance V_alpha, ignored if site_effect="none" or site_effect="fixed" Must be a strictly positive scalar. Default to 0.0005 for weak informative prior.
mu_beta	Means of the Normal priors for the $\beta$ parameters of the suitability process. mu_beta must be either a scalar or a $np$-length vector. If mu_beta takes a scalar value, then that value will serve as the prior mean for all of the $\beta$ parameters. The default value is set to 0 for an uninformative prior, ignored if trait_data is specified.
V_beta	Variances of the Normal priors for the $\beta$ parameters of the suitability process. V_beta must be either a scalar or a $np \times np$ symmetric positive semi-definite square matrix. If V_beta takes a scalar value, then that value will serve as the prior variance for all of the $\beta$ parameters, so the variance covariance matrix used in this case is diagonal with the specified value on the diagonal. The default variance is large and set to 10 for an uninformative flat prior.
mu_gamma	Means of the Normal priors for the $\gamma$ parameters. mu_gamma must be either a scalar, a vector of length $nt$, $np$ or $nt.np$ or a $nt \times np$ matrix. If mu_gamma takes a scalar value, then that value will serve as the prior mean for all of the $\gamma$ parameters. If mu_gamma is a vector of length $nt$ or $nt.np$ the resulting $nt \times np$ matrix is filled by column with specified values, if a $np$-length vector is given, the matrix is filled by row. The default value is set to 0 for an uninformative prior, ignored if trait_data=NULL.
V_gamma	Variances of the Normal priors for the $\gamma$ parameters. V_gamma must be either a scalar, a vector of length $nt$, $np$ or $nt.np$ or a $nt \times np$ positive matrix. If V_gamma takes a scalar value, then that value will serve as the prior variance for all of the $\gamma$ parameters. If V_gamma is a vector of length $nt$ or $nt.np$ the resulting $nt \times np$ matrix is filled by column with specified values, if a $np$-length vector is given, the matrix is filled by row. The default variance is large and set to 10 for an uninformative flat prior, ignored if trait_data=NULL.
mu_lambda	Means of the Normal priors for the $\lambda$ parameters corresponding to the latent variables. mu_lambda must be either a scalar or a $n_{latent}$-length vector. If mu_lambda takes a scalar value, then that value will serve as the prior mean for all of the $\lambda$ parameters. The default value is set to 0 for an uninformative prior.
V_lambda	Variances of the Normal priors for the $\lambda$ parameters corresponding to the latent variables. V_lambda must be either a scalar or a $n_{latent} \times n_{latent}$ symmetric positive semi-definite square matrix. If V_lambda takes a scalar value, then that value will serve as the prior variance for all of $\lambda$ parameters, so the variance covariance matrix used in this case is diagonal with the specified value on the diagonal. The default variance is large and set to 10 for an uninformative flat prior.
ropt	Target acceptance rate for the adaptive Metropolis algorithm. Default to 0.44.
seed	The seed for the random number generator. Default to 1234.
verbose	A switch (0,1) which determines whether or not the progress of the sampler is printed to the screen. Default is 1: a progress bar is printed, indicating the step (in %) reached by the Gibbs sampler.

Details

We model an ecological process where the presence or absence of species $j$ on site $i$ is explained by habitat suitability.

Ecological process :

$y_{ij} \sim \mathcal{B}inomial(\theta_{ij},n_i)$

where

if n_latent=0 and site_effect="none"	logit$(\theta_{ij}) = X_i \beta_j$
if n_latent>0 and site_effect="none"	logit$(\theta_{ij}) = X_i \beta_j + W_i \lambda_j$
if n_latent=0 and site_effect="fixed"	logit$(\theta_{ij}) = X_i \beta_j + \alpha_i$
if n_latent>0 and site_effect="fixed"	logit$(\theta_{ij}) = X_i \beta_j + W_i \lambda_j + \alpha_i$
if n_latent=0 and site_effect="random"	logit$(\theta_{ij}) = 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}) = X_i \beta_j + W_i \lambda_j + \alpha_i$ and $\alpha_i \sim \mathcal{N}(0,V_\alpha)$

In the absence of data on species traits (trait_data=NULL), the effect of species $j$: $\beta_j$; follows the same a priori Gaussian distribution such that $\beta_j \sim \mathcal{N}_{np}(\mu_{\beta},V_{\beta})$, for each species.

If species traits data are provided, the effect of species $j$: $\beta_j$; follows an a priori Gaussian distribution such that $\beta_j \sim \mathcal{N}_{np}(\mu_{\beta_j},V_{\beta})$, where $\mu_{\beta_jp} = \sum_{k=1}^{nt} t_{jk}.\gamma_{kp}$, takes different values for each species.

We assume that $\gamma_{kp} \sim \mathcal{N}(\mu_{\gamma_{kp}},V_{\gamma_{kp}})$ as prior distribution.

We define the matrix $\gamma=(\gamma_{kp})_{k=1,...,nt}^{p=1,...,np}$ such as :

	$x_0$	$x_1$	...	$x_p$	...	$x_{np}$
	__________	________	________	________	________	________
$t_0$ |	$\gamma_{0,0}$	$\gamma_{0,1}$	...	$\gamma_{0,p}$	...	$\gamma_{0,np}$	{ effect of
|	intercept						environmental
|							variables
$t_1$ |	$\gamma_{1,0}$	$\gamma_{1,1}$	...	$\gamma_{1,p}$	...	$\gamma_{1,np}$
... |	...	...	...	...	...	...
$t_k$ |	$\gamma_{k,0}$	$\gamma_{k,1}$	...	$\gamma_{k,p}$	...	$\gamma_{k,np}$
... |	...	...	...	...	...	...
$t_{nt}$ |	$\gamma_{nt,0}$	$\gamma_{nt,1}$	...	$\gamma_{nt,p}$	...	$\gamma_{nt,np}$
	average
	trait effect		interaction	traits	environment

Value

An object of class "jSDM" acting like a list including :

mcmc.alpha	An mcmc object that contains the posterior samples for site effects $\alpha_i$, not returned if site_effect="none".
mcmc.V_alpha	An mcmc object that contains the posterior samples for variance of random site effect, not returned if site_effect="none" or site_effect="fixed".
mcmc.latent	A list by latent variable of mcmc objects that contains the posterior samples for latent variables $W_l$ with $l=1,\ldots,n_{latent}$, not returned if n_latent=0.
mcmc.sp	A list by species of mcmc objects that contains the posterior samples for species effects $\beta_j$ and $\lambda_j$ if n_latent>0.
mcmc.gamma	A list by covariates of mcmc objects that contains the posterior samples for $\gamma_p$ parameters with $p=1,\ldots,np$ if trait_data is specified.
mcmc.Deviance	The posterior sample of the deviance ($D$) is also provided, with $D$ defined as : $D=-2\log(\prod_{ij} P(y_{ij}|\beta_j,\lambda_j, \alpha_i, W_i))$.
logit_theta_latent	Predictive posterior mean of the probability to each species to be present on each site, transformed by logit link function.
theta_latent	Predictive posterior mean of the probability associated to the suitability process for each observation.
model_spec	Various attributes of the model fitted, including the response and model matrix used, distributional assumptions as link function, family and number of latent variables, hyperparameters used in the Bayesian estimation and mcmc, burnin and thin.

The mcmc. objects can be summarized by functions provided by the coda package.

Author(s)

Ghislain Vieilledent <[email protected]>

Jeanne Clément <[email protected]>

References

Gelfand, A. E.; Schmidt, A. M.; Wu, S.; Silander, J. A.; Latimer, A. and Rebelo, A. G. (2005) Modelling species diversity through species level hierarchical modelling. Applied Statistics, 54, 1-20.

Latimer, A. M.; Wu, S. S.; Gelfand, A. E. and Silander, J. A. (2006) Building statistical models to analyze species distributions. Ecological Applications, 16, 33-50.

Ovaskainen, O., Tikhonov, G., Norberg, A., Blanchet, F. G., 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, 561-576.

See Also

plot.mcmc, summary.mcmc jSDM_binomial_probit jSDM_poisson_log

Examples

#==============================================
# jSDM_binomial_logit()
# Example with simulated data
#==============================================

#=================
#== Load libraries
library(jSDM)

#==================
#== Data simulation

#= Number of sites
nsite <- 50
#= Number of species
nsp <- 10
#= Set seed for repeatability
seed <- 1234

#= Number of visits associated to each site
set.seed(seed)
visits <- rpois(nsite,3)
visits[visits==0] <- 1

#= Ecological process (suitability)
x1 <- rnorm(nsite,0,1)
set.seed(2*seed)
x2 <- rnorm(nsite,0,1)
X <- cbind(rep(1,nsite),x1,x2)
np <- ncol(X)
set.seed(3*seed)
W <- cbind(rnorm(nsite,0,1),rnorm(nsite,0,1))
n_latent <- ncol(W)
l.zero <- 0
l.diag <- runif(2,0,2)
l.other <- runif(nsp*2-3,-2,2)
lambda.target <- matrix(c(l.diag[1],l.zero,l.other[1],
                          l.diag[2],l.other[-1]),
                        byrow=TRUE, nrow=nsp)
beta.target <- matrix(runif(nsp*np,-2,2), byrow=TRUE, nrow=nsp)
V_alpha.target <- 0.5
alpha.target <- rnorm(nsite,0,sqrt(V_alpha.target))
logit.theta <- X %*% t(beta.target) + W %*% t(lambda.target) + alpha.target
theta <- inv_logit(logit.theta)
set.seed(seed)
Y <- apply(theta, 2, rbinom, n=nsite, size=visits)

#= Site-occupancy model
# Increase number of iterations (burnin and mcmc) to get convergence
mod <- jSDM_binomial_logit(# Chains
  burnin=150,
  mcmc=150,
  thin=1,
  # Response variable
  presence_data=Y,
  trials=visits,
  # Explanatory variables
  site_formula=~x1+x2,
  site_data=X,
  n_latent=n_latent,
  site_effect="random",
  # Starting values
  beta_start=0,
  lambda_start=0,
  W_start=0,
  alpha_start=0,
  V_alpha=1,
  # Priors
  shape_Valpha=0.5,
  rate_Valpha=0.0005,
  mu_beta=0,
  V_beta=10,
  mu_lambda=0,
  V_lambda=10,
  # Various
  seed=1234,
  ropt=0.44,
  verbose=1)
#==========
#== Outputs

#= Parameter estimates
oldpar <- par(no.readonly = TRUE)
## beta_j
# summary(mod$mcmc.sp$sp_1[,1:ncol(X)])
mean_beta <- matrix(0,nsp,np)
pdf(file=file.path(tempdir(), "Posteriors_beta_jSDM_logit.pdf"))
par(mfrow=c(ncol(X),2))
for (j in 1:nsp) {
  mean_beta[j,] <- apply(mod$mcmc.sp[[j]][,1:ncol(X)],
                         2, mean)
  for (p in 1:ncol(X)) {
    coda::traceplot(mod$mcmc.sp[[j]][,p])
    coda::densplot(mod$mcmc.sp[[j]][,p],
      main = paste(colnames(
        mod$mcmc.sp[[j]])[p],
        ", species : ",j))
    abline(v=beta.target[j,p],col='red')
  }
}
dev.off()

## lambda_j
mean_lambda <- matrix(0,nsp,n_latent)
pdf(file=file.path(tempdir(), "Posteriors_lambda_jSDM_logit.pdf"))
par(mfrow=c(n_latent*2,2))
for (j in 1:nsp) {
  mean_lambda[j,] <- apply(mod$mcmc.sp[[j]]
                           [,(ncol(X)+1):(ncol(X)+n_latent)], 2, mean)
  for (l in 1:n_latent) {
    coda::traceplot(mod$mcmc.sp[[j]][,ncol(X)+l])
    coda::densplot(mod$mcmc.sp[[j]][,ncol(X)+l],
                   main=paste(colnames(mod$mcmc.sp[[j]])
                              [ncol(X)+l],", species : ",j))
    abline(v=lambda.target[j,l],col='red')
  }
}
dev.off()

# Species effects beta and factor loadings lambda
par(mfrow=c(1,2))
plot(beta.target, mean_beta,
     main="species effect beta",
     xlab ="obs", ylab ="fitted")
abline(a=0,b=1,col='red')
plot(lambda.target, mean_lambda,
     main="factor loadings lambda",
     xlab ="obs", ylab ="fitted")
abline(a=0,b=1,col='red')

## W latent variables
par(mfrow=c(1,2))
for (l in 1:n_latent) {
  plot(W[,l],
       summary(mod$mcmc.latent[[paste0("lv_",l)]])[[1]][,"Mean"],
       main = paste0("Latent variable W_", l),
       xlab ="obs", ylab ="fitted")
  abline(a=0,b=1,col='red')
}

## alpha
par(mfrow=c(1,3))
plot(alpha.target, summary(mod$mcmc.alpha)[[1]][,"Mean"],
     xlab ="obs", ylab ="fitted", main="site effect alpha")
abline(a=0,b=1,col='red')
## Valpha
coda::traceplot(mod$mcmc.V_alpha)
coda::densplot(mod$mcmc.V_alpha)
abline(v=V_alpha.target,col='red')

## Deviance
summary(mod$mcmc.Deviance)
plot(mod$mcmc.Deviance)

#= Predictions
par(mfrow=c(1,2))
plot(logit.theta, mod$logit_theta_latent,
     main="logit(theta)",
     xlab="obs", ylab="fitted")
abline(a=0 ,b=1, col="red")
plot(theta, mod$theta_latent,
     main="Probabilities of occurence theta",
     xlab="obs", ylab="fitted")
abline(a=0 ,b=1, col="red")
par(oldpar)

---

Binomial probit regression

Description

The jSDM_binomial_probit function performs a Binomial probit regression in a Bayesian framework. The function calls a Gibbs sampler written in 'C++' code which uses conjugate priors to estimate the conditional posterior distribution of model's parameters.

Usage

jSDM_binomial_probit(
  burnin = 5000,
  mcmc = 10000,
  thin = 10,
  presence_data,
  site_formula,
  trait_data = NULL,
  trait_formula = NULL,
  site_data,
  n_latent = 0,
  site_effect = "none",
  lambda_start = 0,
  W_start = 0,
  beta_start = 0,
  alpha_start = 0,
  gamma_start = 0,
  V_alpha = 1,
  mu_beta = 0,
  V_beta = 10,
  mu_lambda = 0,
  V_lambda = 10,
  mu_gamma = 0,
  V_gamma = 10,
  shape_Valpha = 0.5,
  rate_Valpha = 5e-04,
  seed = 1234,
  verbose = 1
)

Arguments

burnin	The number of burnin iterations for the sampler.
mcmc	The number of Gibbs iterations for the sampler. Total number of Gibbs iterations is equal to burnin+mcmc. burnin+mcmc must be divisible by 10 and superior or equal to 100 so that the progress bar can be displayed.
thin	The thinning interval used in the simulation. The number of mcmc iterations must be divisible by this value.
presence_data	A matrix $n_{site} \times n_{species}$ indicating the presence by a 1 (or the absence by a 0) of each species on each site.
site_formula	A one-sided formula of the form '~x1+...+xp' specifying the explanatory variables for the suitability process of the model, used to form the design matrix $X$ of size $n_{site} \times np$.
trait_data	A data frame containing the species traits which can be included as part of the model. Default to NULL to fit a model without species traits.
trait_formula	A one-sided formula of the form '~ t1 + ... + tk + x1:t1 + ... + xp:tk' specifying the interactions between the environmental variables and the species traits to be considered in the model, used to form the trait design matrix $Tr$ of size $n_{species} \times nt$ and to set to $0$ the $\gamma$ parameters corresponding to interactions not taken into account according to the formula. Default to NULL to fit a model with all possible interactions between species traits found in trait_data and environmental variables defined by site_formula.
site_data	A data frame containing the model's explanatory variables by site.
n_latent	An integer which specifies the number of latent variables to generate. Defaults to $0$.
site_effect	A string indicating whether row effects are included as fixed effects ("fixed"), as random effects ("random"), or not included ("none") in the model. If fixed effects, then for parameter identifiability the first row effect is set to zero, which analogous to acting as a reference level when dummy variables are used. If random effects, they are drawn from a normal distribution with mean zero and unknown variance, analogous to a random intercept in mixed models. Defaults to "none".
lambda_start	Starting values for $\lambda$ parameters corresponding to the latent variables for each species must be either a scalar or a $n_{latent} \times n_{species}$ upper triangular matrix with strictly positive values on the diagonal, ignored if n_latent=0. If lambda_start takes a scalar value, then that value will serve for all of the $\lambda$ parameters except those concerned by the constraints explained above.
W_start	Starting values for latent variables must be either a scalar or a $nsite \times n_latent$ matrix, ignored if n_latent=0. If W_start takes a scalar value, then that value will serve for all of the $W_{il}$ with $i=1,\ldots,n_{site}$ and $l=1,\ldots,n_{latent}$.
beta_start	Starting values for $\beta$ parameters of the suitability process for each species must be either a scalar or a $np \times n_{species}$ matrix. If beta_start takes a scalar value, then that value will serve for all of the $\beta$ parameters.
alpha_start	Starting values for random site effect parameters must be either a scalar or a $n_{site}$-length vector, ignored if site_effect="none". If alpha_start takes a scalar value, then that value will serve for all of the $\alpha$ parameters.
gamma_start	Starting values for $\gamma$ parameters that represent the influence of species-specific traits on species' responses $\beta$, gamma_start must be either a scalar, a vector of length $nt$, $np$ or $nt.np$ or a $nt \times np$ matrix. If gamma_start takes a scalar value, then that value will serve for all of the $\gamma$ parameters. If gamma_start is a vector of length $nt$ or $nt.np$ the resulting $nt \times np$ matrix is filled by column with specified values, if a $np$-length vector is given, the matrix is filled by row.
V_alpha	Starting value for variance of random site effect if site_effect="random" or constant variance of the Gaussian prior distribution for the fixed site effect if site_effect="fixed". Must be a strictly positive scalar, ignored if site_effect="none".
mu_beta	Means of the Normal priors for the $\beta$ parameters of the suitability process. mu_beta must be either a scalar or a $np$-length vector. If mu_beta takes a scalar value, then that value will serve as the prior mean for all of the $\beta$ parameters. The default value is set to 0 for an uninformative prior, ignored if trait_data is specified.
V_beta	Variances of the Normal priors for the $\beta$ parameters of the suitability process. V_beta must be either a scalar or a $np \times np$ symmetric positive semi-definite square matrix. If V_beta takes a scalar value, then that value will serve as the prior variance for all of the $\beta$ parameters, so the variance covariance matrix used in this case is diagonal with the specified value on the diagonal. The default variance is large and set to 10 for an uninformative flat prior.
mu_lambda	Means of the Normal priors for the $\lambda$ parameters corresponding to the latent variables. mu_lambda must be either a scalar or a $n_{latent}$-length vector. If mu_lambda takes a scalar value, then that value will serve as the prior mean for all of the $\lambda$ parameters. The default value is set to 0 for an uninformative prior.
V_lambda	Variances of the Normal priors for the $\lambda$ parameters corresponding to the latent variables. V_lambda must be either a scalar or a $n_{latent} \times n_{latent}$ symmetric positive semi-definite square matrix. If V_lambda takes a scalar value, then that value will serve as the prior variance for all of $\lambda$ parameters, so the variance covariance matrix used in this case is diagonal with the specified value on the diagonal. The default variance is large and set to 10 for an uninformative flat prior.
mu_gamma	Means of the Normal priors for the $\gamma$ parameters. mu_gamma must be either a scalar, a vector of length $nt$, $np$ or $nt.np$ or a $nt \times np$ matrix. If mu_gamma takes a scalar value, then that value will serve as the prior mean for all of the $\gamma$ parameters. If mu_gamma is a vector of length $nt$ or $nt.np$ the resulting $nt \times np$ matrix is filled by column with specified values, if a $np$-length vector is given, the matrix is filled by row. The default value is set to 0 for an uninformative prior, ignored if trait_data=NULL.
V_gamma	Variances of the Normal priors for the $\gamma$ parameters. V_gamma must be either a scalar, a vector of length $nt$, $np$ or $nt.np$ or a $nt \times np$ positive matrix. If V_gamma takes a scalar value, then that value will serve as the prior variance for all of the $\gamma$ parameters. If V_gamma is a vector of length $nt$ or $nt.np$ the resulting $nt \times np$ matrix is filled by column with specified values, if a $np$-length vector is given, the matrix is filled by row. The default variance is large and set to 10 for an uninformative flat prior, ignored if trait_data=NULL.
shape_Valpha	Shape parameter of the Inverse-Gamma prior for the random site effect variance V_alpha, ignored if site_effect="none" or site_effect="fixed". Must be a strictly positive scalar. Default to 0.5 for weak informative prior.
rate_Valpha	Rate parameter of the Inverse-Gamma prior for the random site effect variance V_alpha, ignored if site_effect="none" or site_effect="fixed" Must be a strictly positive scalar. Default to 0.0005 for weak informative prior.
seed	The seed for the random number generator. Default to 1234.
verbose	A switch (0,1) which determines whether or not the progress of the sampler is printed to the screen. Default is 1: a progress bar is printed, indicating the step (in %) reached by the Gibbs sampler.

Details

We model an ecological process where the presence or absence of species $j$ on site $i$ is explained by habitat suitability.

Ecological process:

$y_{ij} \sim \mathcal{B}ernoulli(\theta_{ij})$

where

if n_latent=0 and site_effect="none"	probit$(\theta_{ij}) = X_i \beta_j$
if n_latent>0 and site_effect="none"	probit$(\theta_{ij}) = X_i \beta_j + W_i \lambda_j$
if n_latent=0 and site_effect="random"	probit$(\theta_{ij}) = 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}) = X_i \beta_j + W_i \lambda_j + \alpha_i$
if n_latent=0 and site_effect="fixed"	probit$(\theta_{ij}) = X_i \beta_j + \alpha_i$
if n_latent>0 and site_effect="random"	probit$(\theta_{ij}) = X_i \beta_j + W_i \lambda_j + \alpha_i$ and $\alpha_i \sim \mathcal{N}(0,V_\alpha)$

In the absence of data on species traits (trait_data=NULL), the effect of species $j$: $\beta_j$; follows the same a priori Gaussian distribution such that $\beta_j \sim \mathcal{N}_{np}(\mu_{\beta},V_{\beta})$, for each species.

If species traits data are provided, the effect of species $j$: $\beta_j$; follows an a priori Gaussian distribution such that $\beta_j \sim \mathcal{N}_{np}(\mu_{\beta_j},V_{\beta})$, where $\mu_{\beta_jp} = \sum_{k=1}^{nt} t_{jk}.\gamma_{kp}$, takes different values for each species.

We assume that $\gamma_{kp} \sim \mathcal{N}(\mu_{\gamma_{kp}},V_{\gamma_{kp}})$ as prior distribution.

We define the matrix $\gamma=(\gamma_{kp})_{k=1,...,nt}^{p=1,...,np}$ such as :

	$x_0$	$x_1$	...	$x_p$	...	$x_{np}$
	__________	________	________	________	________	________
$t_0$ |	$\gamma_{0,0}$	$\gamma_{0,1}$	...	$\gamma_{0,p}$	...	$\gamma_{0,np}$	{ effect of
|	intercept						environmental
|							variables
$t_1$ |	$\gamma_{1,0}$	$\gamma_{1,1}$	...	$\gamma_{1,p}$	...	$\gamma_{1,np}$
... |	...	...	...	...	...	...
$t_k$ |	$\gamma_{k,0}$	$\gamma_{k,1}$	...	$\gamma_{k,p}$	...	$\gamma_{k,np}$
... |	...	...	...	...	...	...
$t_{nt}$ |	$\gamma_{nt,0}$	$\gamma_{nt,1}$	...	$\gamma_{nt,p}$	...	$\gamma_{nt,np}$
	average
	trait effect		interaction	traits	environment

Value

An object of class "jSDM" acting like a list including :

mcmc.alpha	An mcmc object that contains the posterior samples for site effects $\alpha$, not returned if site_effect="none".
mcmc.V_alpha	An mcmc object that contains the posterior samples for variance of random site effect, not returned if site_effect="none" or site_effect="fixed".
mcmc.latent	A list by latent variable of mcmc objects that contains the posterior samples for latent variables $W_l$ with $l=1,\ldots,n_{latent}$, not returned if n_latent=0.
mcmc.sp	A list by species of mcmc objects that contains the posterior samples for species effects $\beta_j$ and $\lambda_j$ if n_latent>0 with $j=1,\ldots,n_{species}$.
mcmc.gamma	A list by covariates of mcmc objects that contains the posterior samples for $\gamma_p$ parameters with $p=1,\ldots,np$ if trait_data is specified.
mcmc.Deviance	The posterior sample of the deviance ($D$) is also provided, with $D$ defined as : $D=-2\log(\prod_{ij} P(y_{ij}|\beta_j,\lambda_j, \alpha_i, W_i))$.
Z_latent	Predictive posterior mean of the latent variable Z.
probit_theta_latent	Predictive posterior mean of the probability to each species to be present on each site, transformed by probit link function.
theta_latent	Predictive posterior mean of the probability to each species to be present on each site.
model_spec	Various attributes of the model fitted, including the response and model matrix used, distributional assumptions as link function, family and number of latent variables, hyperparameters used in the Bayesian estimation and mcmc, burnin and thin.

The mcmc. objects can be summarized by functions provided by the coda package.

Author(s)

Ghislain Vieilledent <[email protected]>

Jeanne Clément <[email protected]>

References

Chib, S. and Greenberg, E. (1998) Analysis of multivariate probit models. Biometrika, 85, 347-361.

Warton, D. I.; Blanchet, F. G.; O'Hara, R. B.; O'Hara, R. B.; Ovaskainen, O.; Taskinen, S.; Walker, S. C. and Hui, F. K. C. (2015) So Many Variables: Joint Modeling in Community Ecology. Trends in Ecology & Evolution, 30, 766-779.

Ovaskainen, O., Tikhonov, G., Norberg, A., Blanchet, F. G., 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, 561-576.

See Also

plot.mcmc, summary.mcmc jSDM_binomial_logit jSDM_poisson_log jSDM_binomial_probit_sp_constrained

Examples

#======================================
# jSDM_binomial_probit()
# Example with simulated data
#====================================

#=================
#== Load libraries
library(jSDM)

#==================
#' #== Data simulation

#= Number of sites
nsite <- 150

#= Set seed for repeatability
seed <- 1234
set.seed(seed)

#= Number of species
nsp<- 20

#= Number of latent variables
n_latent <- 2

#= Ecological process (suitability)
x1 <- rnorm(nsite,0,1)
x2 <- rnorm(nsite,0,1)
X <- cbind(rep(1,nsite),x1,x2)
np <- ncol(X)
#= Latent variables W
W <- matrix(rnorm(nsite*n_latent,0,1), nsite, n_latent)
#= Fixed species effect beta 
beta.target <- t(matrix(runif(nsp*np,-2,2),
                        byrow=TRUE, nrow=nsp))
#= Factor loading lambda  
lambda.target <- matrix(0, n_latent, nsp)
mat <- t(matrix(runif(nsp*n_latent, -2, 2), byrow=TRUE, nrow=nsp))
lambda.target[upper.tri(mat, diag=TRUE)] <- mat[upper.tri(mat, diag=TRUE)]
diag(lambda.target) <- runif(n_latent, 0, 2)
#= Variance of random site effect 
V_alpha.target <- 0.5
#= Random site effect alpha
alpha.target <- rnorm(nsite,0 , sqrt(V_alpha.target))
# Simulation of response data with probit link
probit_theta <- X%*%beta.target + W%*%lambda.target + alpha.target
theta <- pnorm(probit_theta)
e <- matrix(rnorm(nsp*nsite,0,1),nsite,nsp)
# Latent variable Z 
Z_true <- probit_theta + e
# Presence-absence matrix Y
Y <- matrix (NA, nsite,nsp)
for (i in 1:nsite){
  for (j in 1:nsp){
    if ( Z_true[i,j] > 0) {Y[i,j] <- 1}
    else {Y[i,j] <- 0}
  }
}

#==================================
#== Site-occupancy model

# Increase number of iterations (burnin and mcmc) to get convergence
mod<-jSDM_binomial_probit(# Iteration
                           burnin=200,
                           mcmc=200,
                           thin=1,
                           # Response variable
                           presence_data=Y,
                           # Explanatory variables
                           site_formula=~x1+x2,
                           site_data = X,
                           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.5,
                           rate_Valpha=0.0005,
                           mu_beta=0, V_beta=1,
                           mu_lambda=0, V_lambda=1,
                           seed=1234, verbose=1)
# ===================================================
# Result analysis
# ===================================================

#==========
#== Outputs

oldpar <- par(no.readonly = TRUE) 

#= Parameter estimates

## beta_j
mean_beta <- matrix(0,nsp,ncol(X))
pdf(file=file.path(tempdir(), "Posteriors_beta_jSDM_probit.pdf"))
par(mfrow=c(ncol(X),2))
for (j in 1:nsp) {
  mean_beta[j,] <- apply(mod$mcmc.sp[[j]]
                         [,1:ncol(X)], 2, mean)  
  for (p in 1:ncol(X)){
    coda::traceplot(mod$mcmc.sp[[j]][,p])
    coda::densplot(mod$mcmc.sp[[j]][,p],
      main = paste(colnames(mod$mcmc.sp[[j]])[p],", species : ",j))
    abline(v=beta.target[p,j],col='red')
  }
}
dev.off()

## lambda_j
mean_lambda <- matrix(0,nsp,n_latent)
pdf(file=file.path(tempdir(), "Posteriors_lambda_jSDM_probit.pdf"))
par(mfrow=c(n_latent*2,2))
for (j in 1:nsp) {
  mean_lambda[j,] <- apply(mod$mcmc.sp[[j]]
                           [,(ncol(X)+1):(ncol(X)+n_latent)], 2, mean)  
  for (l in 1:n_latent) {
    coda::traceplot(mod$mcmc.sp[[j]][,ncol(X)+l])
    coda::densplot(mod$mcmc.sp[[j]][,ncol(X)+l],
                   main=paste(colnames(mod$mcmc.sp[[j]])
                              [ncol(X)+l],", species : ",j))
    abline(v=lambda.target[l,j],col='red')
  }
}
dev.off()

# Species effects beta and factor loadings lambda
par(mfrow=c(1,2))
plot(t(beta.target), mean_beta,
     main="species effect beta",
     xlab ="obs", ylab ="fitted")
abline(a=0,b=1,col='red')
plot(t(lambda.target), mean_lambda,
     main="factor loadings lambda",
     xlab ="obs", ylab ="fitted")
abline(a=0,b=1,col='red')

## W latent variables
par(mfrow=c(1,2))
for (l in 1:n_latent) {
  plot(W[,l],
       summary(mod$mcmc.latent[[paste0("lv_",l)]])[[1]][,"Mean"],
       main = paste0("Latent variable W_", l),
       xlab ="obs", ylab ="fitted")
  abline(a=0,b=1,col='red')
}

## alpha
par(mfrow=c(1,3))
plot(alpha.target, summary(mod$mcmc.alpha)[[1]][,"Mean"],
     xlab ="obs", ylab ="fitted", main="site effect alpha")
abline(a=0,b=1,col='red')
## Valpha
coda::traceplot(mod$mcmc.V_alpha)
coda::densplot(mod$mcmc.V_alpha)
abline(v=V_alpha.target,col='red')

## Deviance
summary(mod$mcmc.Deviance)
plot(mod$mcmc.Deviance)

#= Predictions

## Z
par(mfrow=c(1,2))
plot(Z_true,mod$Z_latent,
     main="Z_latent", xlab="obs", ylab="fitted")
abline(a=0,b=1,col='red')

## probit_theta
plot(probit_theta,mod$probit_theta_latent,
     main="probit(theta)",xlab="obs",ylab="fitted")
abline(a=0,b=1,col='red')

## probabilities theta
par(mfrow=c(1,1))
plot(theta,mod$theta_latent,
     main="theta",xlab="obs",ylab="fitted")
abline(a=0,b=1,col='red')

par(oldpar)

---