get_residual_cor function

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

This function use coefficients $(\lambda_jl$ with $j=1,...,n_species$ and $l=1,...,n_latent)$, corresponding to latent variables fitted using jSDM package, to calculate the variance-covariance matrix which controls correlation between species.

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.

Returns

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 x n_species$ matrix containing the lower limits of the ($100 x prob$) HPD interval of variance-covariance matrices over the MCMC samples.

• cov.upper: A $n_species x n_species$ matrix containing the upper limits of the ($100 x prob$) HPD interval of variance-covariance matrices over the MCMC samples.

• cov.sig: A $n_species x 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 x 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 x n_species$ matrix containing the lower limits of the ($100 x prob$) HPD interval of correlation matrices over the MCMC samples.

• cor.upper: A $n_species x n_species$ matrix containing the upper limits of the ($100 x prob$) HPD interval of correlation matrices over the MCMC samples.

• cor.sig: A $n_species x n_species$ matrix containing the value $1$ corresponding to the “significant" correlations and the value $0$ corresponding to "non-significant" correlations, whose ($100 x prob$) HPD interval over the MCMC samples contain zero.

Details

After fitting the jSDM with latent variables, the fullspecies residual correlation matrix : $R=(R_ij)$ with $i=1,..., n_species$ and $j=1,..., 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)y_ij ~ N(\theta_ij, V),$$

. Then we compute correlations from covariances :

$$R_{ij} = \frac{\Sigma_{ij}}{\sqrt{\Sigma_ii\Sigma _jj}}R_ij = Sigma_ij / sqrt(Sigma_ii.Sigma _jj)$$

.

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

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

Author(s)

Ghislain Vieilledent ghislain.vieilledent@cirad.fr

Jeanne Clément jeanne.clement16@laposte.net