The simplest Species Distribution Model in OpenBUGS & R

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)

# 4.2 - loading the R2OpenBUGS package
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

res.summary <- read.bugs(res)
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")