#==============================================
# jSDM_poisson_log()
# 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
#= Ecological process (suitability)
set.seed(seed)
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,1)
l.other <- runif(nsp*2-3,-1,1)
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,-1,1), byrow=TRUE, nrow=nsp)
V_alpha.target <- 0.5
alpha.target <- rnorm(nsite,0,sqrt(V_alpha.target))
log.theta <- X %*% t(beta.target) + W %*% t(lambda.target) + alpha.target
theta <- exp(log.theta)
Y <- apply(theta, 2, rpois, n=nsite)
#= Site-occupancy model
# Increase number of iterations (burnin and mcmc) to get convergence
mod <- jSDM_poisson_log(# Chains
burnin=200,
mcmc=200,
thin=1,
# Response variable
count_data=Y,
# 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
oldpar <- par(no.readonly = TRUE)
#= Parameter estimates
## beta_j
mean_beta <- matrix(0,nsp,np)
pdf(file=file.path(tempdir(), "Posteriors_beta_jSDM_log.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_log.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(log.theta, mod$log_theta_latent,
main="log(theta)",
xlab="obs", ylab="fitted")
abline(a=0 ,b=1, col="red")
plot(theta, mod$theta_latent,
main="Expected abundance theta",
xlab="obs", ylab="fitted")
abline(a=0 ,b=1, col="red")
par(oldpar)
Run the code above in your browser using DataLab