#======================================
# jSDM_binomial_probit_sp_constrained()
# Example with simulated data
#====================================
#=================
#== Load libraries
library(jSDM)
#==================
#== Data simulation
#= Number of sites
nsite <- 30
#= Set seed for repeatability
seed <- 1234
set.seed(seed)
#= Number of species
nsp <- 10
#= 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_sp_constrained(# Iteration
burnin=100,
mcmc=100,
thin=1,
# parallel MCMCs
nchains=2, ncores=2,
# 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=c(123,1234), verbose=1)
# ===================================================
# Result analysis
# ===================================================
#==========
#== Outputs
oldpar <- par(no.readonly = TRUE)
burnin <- mod[[1]]$model_spec$burnin
ngibbs <- burnin + mod[[1]]$model_spec$mcmc
thin <- mod[[1]]$model_spec$thin
require(coda)
arr2mcmc <- function(x) {
return(mcmc(as.data.frame(x),
start=burnin+1 , end=ngibbs, thin=thin))
}
mcmc_list_Deviance <- mcmc.list(lapply(lapply(mod,"[[","mcmc.Deviance"), arr2mcmc))
mcmc_list_alpha <- mcmc.list(lapply(lapply(mod,"[[","mcmc.alpha"), arr2mcmc))
mcmc_list_V_alpha <- mcmc.list(lapply(lapply(mod,"[[","mcmc.V_alpha"), arr2mcmc))
mcmc_list_param <- mcmc.list(lapply(lapply(mod,"[[","mcmc.sp"), arr2mcmc))
mcmc_list_lv <- mcmc.list(lapply(lapply(mod,"[[","mcmc.latent"), arr2mcmc))
mcmc_list_beta <- mcmc_list_param[,grep("beta",colnames(mcmc_list_param[[1]]))]
mcmc_list_beta0 <- mcmc_list_beta[,grep("Intercept", colnames(mcmc_list_beta[[1]]))]
mcmc_list_lambda <- mcmc.list(
lapply(mcmc_list_param[, grep("lambda", colnames(mcmc_list_param[[1]]))],
arr2mcmc))
# Compute Rhat
psrf_alpha <- mean(gelman.diag(mcmc_list_alpha, multivariate=FALSE)$psrf[,2])
psrf_V_alpha <- gelman.diag(mcmc_list_V_alpha)$psrf[,2]
psrf_beta0 <- mean(gelman.diag(mcmc_list_beta0)$psrf[,2])
psrf_beta <- mean(gelman.diag(mcmc_list_beta[,grep("Intercept",
colnames(mcmc_list_beta[[1]]),
invert=TRUE)])$psrf[,2])
psrf_lambda <- mean(gelman.diag(mcmc_list_lambda,
multivariate=FALSE)$psrf[,2],
na.rm=TRUE)
psrf_lv <- mean(gelman.diag(mcmc_list_lv, multivariate=FALSE)$psrf[,2])
Rhat <- data.frame(Rhat=c(psrf_alpha, psrf_V_alpha, psrf_beta0, psrf_beta,
psrf_lambda, psrf_lv),
Variable=c("alpha", "Valpha", "beta0", "beta",
"lambda", "W"))
# Barplot
library(ggplot2)
ggplot2::ggplot(Rhat, aes(x=Variable, y=Rhat)) +
ggtitle("Averages of Rhat obtained with jSDM for each type of parameter") +
theme(plot.title = element_text(hjust = 0.5, size=15)) +
geom_bar(fill="skyblue", stat = "identity") +
geom_text(aes(label=round(Rhat,2)), vjust=0, hjust=-0.1) +
geom_hline(yintercept=1, color='red') +
coord_flip()
#= Parameter estimates
nchains <- length(mod)
## beta_j
pdf(file=file.path(tempdir(), "Posteriors_species_effect_jSDM_probit.pdf"))
plot(mcmc_list_param)
dev.off()
## Latent variables
pdf(file=file.path(tempdir(), "Posteriors_latent_variables_jSDM_probit.pdf"))
plot(mcmc_list_lv)
dev.off()
## Random site effect and its variance
pdf(file=file.path(tempdir(), "Posteriors_site_effect_jSDM_probit.pdf"))
plot(mcmc_list_V_alpha)
plot(mcmc_list_alpha)
dev.off()
## Predictive posterior mean for each observation
# Species effects beta and factor loadings lambda
param <- list()
for (i in 1:nchains){
param[[i]] <- matrix(unlist(lapply(mod[[i]]$mcmc.sp,colMeans)),
nrow=nsp, byrow=TRUE)
}
par(mfrow=c(1,1))
for (i in 1:nchains){
if(i==1){
plot(t(beta.target), param[[i]][,1:np],
main="species effect beta",
xlab ="obs", ylab ="fitted")
abline(a=0,b=1, col='red')
}
else{
points(t(beta.target), param[[i]][,1:np], col=i)
}
}
par(mfrow=c(1,2))
for(l in 1:n_latent){
for (i in 1:nchains){
if (i==1){
plot(t(lambda.target)[,l],
param[[i]][,np+l],
main=paste0("factor loadings lambda_", l),
xlab ="obs", ylab ="fitted")
abline(a=0,b=1, col='red')
} else {
points(t(lambda.target)[,l],
param[[i]][,np+l],
col=i)
}
}
}
## W latent variables
mean_W <- list()
for (i in 1:nchains){
mean_W[[i]] <- sapply(mod[[i]]$mcmc.latent,colMeans)
}
par(mfrow=c(1,2))
for (l in 1:n_latent) {
for (i in 1:nchains){
if (i==1){
plot(W[,l], mean_W[[i]][,l],
main = paste0("Latent variable W_", l),
xlab ="obs", ylab ="fitted")
abline(a=0,b=1, col='red')
}
else{
points(W[,l], mean_W[[i]][,l], col=i)
}
}
}
#= W.lambda
par(mfrow=c(1,2))
for (i in 1:nchains){
if (i==1){
plot(W%*%lambda.target,
mean_W[[i]]%*%t(param[[i]][,(np+1):(np+n_latent)]),
main = "W.lambda",
xlab ="obs", ylab ="fitted")
abline(a=0,b=1, col='red')
}
else{
points(W%*%lambda.target,
mean_W[[i]]%*%t(param[[i]][,(np+1):(np+n_latent)]),
col=i)
}
}
# Site effect alpha et V_alpha
plot(alpha.target,colMeans(mod[[1]]$mcmc.alpha),
xlab="obs", ylab="fitted",
main="Random site effect alpha")
abline(a=0,b=1, col='red')
points(V_alpha.target, mean(mod[[1]]$mcmc.V_alpha),
pch=18, cex=2)
legend("bottomright", legend=c("V_alpha"), pch =18, pt.cex=1.5)
for (i in 2:nchains){
points(alpha.target, colMeans(mod[[i]]$mcmc.alpha), col=i)
points(V_alpha.target, mean(mod[[i]]$mcmc.V_alpha),
pch=18, col=i, cex=2)
}
#= Predictions
## Occurence probabilities
plot(pnorm(probit_theta), mod[[1]]$theta_latent,
main="theta",xlab="obs",ylab="fitted")
for (i in 2:nchains){
points(pnorm(probit_theta), mod[[i]]$theta_latent, col=i)
}
abline(a=0,b=1, col='red')
## probit(theta)
plot(probit_theta, mod[[1]]$probit_theta_latent,
main="probit(theta)",xlab="obs",ylab="fitted")
for (i in 2:nchains){
points(probit_theta, mod[[i]]$probit_theta_latent, col=i)
}
abline(a=0,b=1, col='red')
## Deviance
plot(mcmc_list_Deviance)
par(oldpar)
Run the code above in your browser using DataLab