The simplest Species Distribution Model in OpenBUGS & R
[This article was first published on Are you cereal? » R, and kindly contributed to R-bloggers]. (You can report issue about the content on this page here)
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.
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) # 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 their blog: Are you cereal? » R.
R-bloggers.com offers daily e-mail updates about R news and tutorials about learning R and many other topics. Click here if you're looking to post or find an R/data-science job.
Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.