#======================================
# 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)
Run the code above in your browser using DataLab