Here we will go through a simple hierarchical model of species occupancy that accounts for imperfect detectability of a species.

1 The Data

We will work with Marc Kery’s data from chapter 20 of his Introduction to WinBUGS for ecologists. The data describe detections of Gentianella germanica at 150 sites, each visited 3 times. We wish to estimate the relationship between \(P\) of occurence, and a humidity index. However, we also know that the detectability of the species depends on humidity itself.

Loading the data from the web:

  gentiana <- read.csv("http://www.petrkeil.com/wp-content/uploads/2014/02/gentiana.csv")
  gentiana <- gentiana[,-1]

Explore the data a little bit:

head(gentiana)
##   humidity visit1 visit2 visit3 pres.abs
## 1    -0.99      0      0      0        0
## 2    -0.98      0      0      0        0
## 3    -0.96      0      0      0        0
## 4    -0.93      0      0      0        0
## 5    -0.89      0      0      0        0
## 6    -0.88      0      0      0        0
nrow(gentiana)
## [1] 150

1.1 The “naive” analysis by GLM:

This is the classical LOGISTIC REGRESSION for presence-absence data:

  naive.m1 <- glm(pres.abs~humidity, data=gentiana, family="binomial")

Or the response can be number of successes out of 3 visits:

  succ.fail <- cbind(rowSums(gentiana[,2:4]), rowSums(gentiana[,2:4]))
  naive.m2 <- glm(succ.fail~humidity, data=gentiana, family="binomial")

Plotting the data and the predictions of the two naive models:

  plot(gentiana$humidity, gentiana$pres.abs,
       xlab="Humidity index", ylab="Gentiana detections")
  lines(gentiana$humidity, predict(naive.m1, type="response"), col="blue")
  lines(gentiana$humidity, predict(naive.m2, type="response"), col="red")
  legend("right", lwd=c(2,2), col=c("blue","red"),
         legend=c("model1","model2"))