Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

This post demonstrates the simplest Species Distribution Model based on logistic regression for presence/absence data. I heavily simplified the example from Kéry (2010): Introduction to WinBUGS for Ecologists, Chapter 20.

```
# -------------------------------------------------------
# 1. SOME USEFUL FUNCTIONS

# logit function
logit <- function(x) log(x/(1-x))

# "unlogit" function
unlogit <- function(x) exp(x)/(1+exp(x))

# -------------------------------------------------------
# 2. GENERATING ARTIFICIAL DATA

n.site <- 150  # 150 sites visited
humidity <- sort(runif(n = n.site, min = -1, max =1))

alpha.occ <- 0 # Logit-linear intercept
beta.occ  <- 3 # Logit-linear slope

occ.prob <- unlogit(alpha.occ + beta.occ*humidity)
true.presence <- rbinom(n=n.site, size=1, prob=occ.prob)

plot(true.presence ~ humidity)
lines(occ.prob ~ humidity)

# -------------------------------------------------------
# 3. ORDINARY GLM ANALYSIS
m1 <- glm(true.presence ~ humidity, family="binomial")
summary(m1)
lines(unlogit(predict(m1))~humidity, col="red")

# -------------------------------------------------------
# 4. BAYESIAN GLM ANALYSIS

# 4.1 - putting the data in the OpenBugs-friendly format
my.data <- list(humidity=humidity,
true.presence=true.presence,
n.site=n.site)

library(R2OpenBUGS)

# 4.3 - sinking the MODEL into a file
sink("sdm.txt")
cat("
model
{
# priors
alpha.occ ~ dnorm(0,0.001)
beta.occ ~ dnorm(0,0.001)

# likelihood
for(i in 1:n.site)
{
logit(p[i]) <- alpha.occ + beta.occ*humidity[i]
true.presence[i] ~ dbern(p[i])
}
}
")
sink()

# 4.4 - setting the PARAMETERS to be monitored
params <- c("alpha.occ", "beta.occ")

# 4.5 - specifying the INITIAL MCMC VALUES
# Note what happens if you specify too broad initial distributions.
inits <- function()
{
list (alpha.occ=rnorm(1,0,10),
beta.occ=rnorm(1,0,10))
}

# 4.6 - calling OpenBugs to sample from the posterior distributions
# (you may need to change the OpenBUGS.pgm directory)
res <- bugs(data=my.data,
inits=inits,
parameters.to.save=params,
n.iter=2000,
n.chains=3,
n.burnin=1000,
model.file="sdm.txt",
debug=TRUE,
codaPkg=TRUE,
OpenBUGS.pgm="C:/Program Files (x86)/OpenBUGS/OpenBUGS321/OpenBUGS.exe",
working.directory=getwd())

# -------------------------------------------------------
# 4. GETTING OpenBUGS RESULTS BACK TO R

plot(res.summary)
quant <- summary(res.summary)\$quantiles
alpha.est <- quant[1,3]
beta.est <- quant[2,3]

# probabilities predicted by the Bayesian model
est.prob <- unlogit(alpha.est + beta.est*humidity)

# plotting everything in one graph
plot(true.presence ~ humidity)
lines(occ.prob ~ humidity)
lines(unlogit(predict(m1))~humidity, col="red")
lines(est.prob ~ humidity, col="green")

```