The simplest Species Distribution Model in OpenBUGS & R

September 27, 2012
By

(This article was first published on Are you cereal? » R, and kindly contributed to R-bloggers)

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")

To leave a comment for the author, please follow the link and comment on his blog: Are you cereal? » R.

R-bloggers.com offers daily e-mail updates about R news and tutorials on topics such as: visualization (ggplot2, Boxplots, maps, animation), programming (RStudio, Sweave, LaTeX, SQL, Eclipse, git, hadoop, Web Scraping) statistics (regression, PCA, time series, trading) and more...



If you got this far, why not subscribe for updates from the site? Choose your flavor: e-mail, twitter, RSS, or facebook...

Comments are closed.