hSDM.binomial.RdThe hSDM.binomial 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.
hSDM.binomial(presences, trials, suitability, data,
suitability.pred = NULL, burnin = 5000, mcmc = 10000, thin = 10,
beta.start, mubeta = 0, Vbeta = 1e+06, seed = 1234, verbose = 1, save.p
= 0)A vector indicating the number of successes (or presences) for each observation.
A vector indicating the number of trials for each observation. \(t_n\) should be superior or equal to \(y_n\), the number of successes for observation \(n\). If \(t_n=0\), then \(y_n=0\).
A one-sided formula of the form '~x1+...+xp' with p terms specifying the explicative variables for the suitability process of the model.
A data frame containing the model's explicative variables.
An optional data frame in which to look for variables with which to predict. If NULL, the observations are used.
The number of burnin iterations for the sampler.
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.
The thinning interval used in the simulation. The number of mcmc iterations must be divisible by this value.
Starting values for beta parameters of the
suitability process. If beta.start takes a scalar value, then
that value will serve for all of the betas.
Means of the priors for the \(\beta\) parameters
of the suitability process. mubeta must be either a scalar or a
p-length vector. If mubeta takes a scalar value, then that value will
serve as the prior mean for all of the betas. The default value is set
to 0 for an uninformative prior.
Variances of the Normal priors for the \(\beta\)
parameters of the suitability process. Vbeta must be either a
scalar or a p-length vector. If Vbeta takes a scalar value,
then that value will serve as the prior variance for all of the
betas. The default variance is large and set to 1.0E6 for an
uninformative flat prior.
The seed for the random number generator. Default to 1234.
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.
A switch (0,1) which determines whether or not the
sampled values for predictions are saved. Default is 0: the
posterior mean is computed and returned in the theta.pred
vector. Be careful, setting save.p to 1 might require a large
amount of memory.
An mcmc object that contains the posterior sample. This object can be summarized by functions provided by the coda package. The posterior sample of the deviance \(D\), with \(D=-2\log(\prod_i P(y_i|\beta,t_i))\), is also provided.
If save.p is set to 0 (default),
theta.pred is the predictive posterior mean of the
probability associated to the suitability process for each
prediction. If save.p is set to 1, theta.pred is an
mcmc object with sampled values of the probability associated
to the suitability process for each prediction.
Predictive posterior mean of the probability associated to the suitability process for each observation.
We model an ecological process where the presence or absence of the species is explained by habitat suitability.
Ecological process: $$y_i \sim \mathcal{B}inomial(\theta_i,t_i)$$ $$logit(\theta_i) = X_i \beta$$
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.
if (FALSE) {
#==============================================
# hSDM.binomial()
# Example with simulated data
#==============================================
#=================
#== Load libraries
library(hSDM)
#==================
#== Data simulation
#= Number of sites
nsite <- 200
#= 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)
set.seed(seed)
x1 <- rnorm(nsite,0,1)
set.seed(2*seed)
x2 <- rnorm(nsite,0,1)
X <- cbind(rep(1,nsite),x1,x2)
beta.target <- c(-1,1,-1)
logit.theta <- X %*% beta.target
theta <- inv.logit(logit.theta)
set.seed(seed)
Y <- rbinom(nsite,visits,theta)
#= Data-sets
data.obs <- data.frame(Y,visits,x1,x2)
#==================================
#== Site-occupancy model
mod.hSDM.binomial <- hSDM.binomial(presences=data.obs$Y,
trials=data.obs$visits,
suitability=~x1+x2,
data=data.obs,
suitability.pred=NULL,
burnin=1000, mcmc=1000, thin=1,
beta.start=0,
mubeta=0, Vbeta=1.0E6,
seed=1234, verbose=1,
save.p=0)
#==========
#== Outputs
#= Parameter estimates
summary(mod.hSDM.binomial$mcmc)
pdf(file="Posteriors_hSDM.binomial.pdf")
plot(mod.hSDM.binomial$mcmc)
dev.off()
#== glm resolution to compare
mod.glm <- glm(cbind(Y,visits-Y)~x1+x2,family="binomial",data=data.obs)
summary(mod.glm)
#= Predictions
summary(mod.hSDM.binomial$theta.latent)
summary(mod.hSDM.binomial$theta.pred)
pdf(file="Pred-Init.pdf")
plot(theta,mod.hSDM.binomial$theta.pred)
abline(a=0,b=1,col="red")
dev.off()
}