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

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`

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