In the last section we looked at a simple single-species, single-season closed occupancy model. We still keep those assumptions, but now we model probability of occupancy as a function of habitat covariates.

### Swiss willow tits

Our example data come from the Swiss Breeding Bird Survey (“Monitoring Häufige Brutvögel” MHB). Surveys are done 2-3 times per year during the breeding season at 267 1-km square quadrats across Switzerland. We will use the data for willow tits (*Poecile montanus*) analysed by Royle and Dorazio (2008, p.87ff)^{1}Royle, J.A. & Dorazio, R.M. (2008) *Hierarchical modeling and inference in ecology: the analysis of data from populations, metapopulations and communities*. Academic Press..

You can download and check the data with the following code:

**# Read in data
wt <- read.csv("http://bcss.org.my/data/willowtits.csv", comment="#")
str(wt)
#'data.frame': 237 obs. of 15 variables:
# $ y.1 : int 0 0 0 0 0 0 0 0 1 0 ...
# $ y.2 : int 0 0 0 0 0 0 0 0 1 0 ...
# $ y.3 : int 0 0 0 0 0 0 0 0 0 1 ...
# $ c.1 : int 0 0 0 0 0 0 0 0 2 0 ...
# $ c.2 : int 0 0 0 0 0 0 0 0 1 0 ...
# $ c.3 : int 0 0 0 0 0 0 0 0 0 1 ...
# $ elev : int 420 450 1050 1110 510 630 590 530 ...
# $ forest: int 3 21 32 35 2 60 5 13 50 57 ...
# $ dur.1 : int 240 160 120 180 210 150 115 155 165 ...
# $ dur.2 : int 240 155 105 170 225 180 105 140 165 ...
# $ dur.3 : int 240 140 85 145 235 195 95 135 180 ...
# $ day.1 : int 29 13 30 23 28 17 16 24 25 21 ...
# $ day.2 : int 58 39 47 44 56 56 37 47 46 38 ...
# $ day.3 : int 73 62 74 71 73 73 76 74 70 50 ...
# $ length: num 6.2 5.1 4.3 5.4 3.6 6.1 5.1 3.7 ...**

Here’s what the variables are:

- y.1 to y.3 : Detection/nondetection (1/0) on each of three visits to each quadrat;
- c.1 to c.3 : Counts of territories on each of the three visits;
- elev : quadrat elevation (m);
- forest : forest cover (%);
- dur.1 to dur.3 : duration of each visit (mins);
- day.1 to day.3 : Julian date of each visit, 1 = 1st April;
- length : length of the transect walked within the quadrat (km).

In this section we will use the detection/nondetection data (y.1 to y.3) and the elevation (elev) and forest covariates.

### Data preparation

Since we are not using the survey-level covariates (dur.1-dur.3 and day.1-day.3), we can simplify the model by aggregating the detection data to produce a vector with the number of visits to each site $n$ and the number of visits when willow tits were detected $y$.

**# Extract and aggregate the detection data:
Y <- wt[, 1:3]
( n <- rowSums(!is.na(Y)) ) # number of visits
#[1] 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3 3...
( y <- rowSums(Y, na.rm=TRUE) ) # number with detection
#[1] 0 0 0 0 0 0 0 0 2 1 0 0 0 0 1 2 0 3 0 1 0 0...**

We will also standardise the elevation and forest data so that our predictors have mean 0 and standard deviation (SD) 1. We use the `standardize`

function in the `wiqid`

package.

**# Standardise the covariates:
elevS <- standardize(wt$elev)
mean(elevS)
# [1] 8.442321e-17
sd(elevS)
# [1] 1
forestS <- standardize(wt$forest)
mean(forestS)
# [1] 8.606644e-17
sd(forestS)
# [1] 1**

I like to add an ‘S’ to the end of the variable to indicate that it’s standardised. The means are not quite zero because the SD is adjusted after centering.

### The model

As with the salamander model, $z_i$ is drawn from a Bernoulli distribution. But now the probability of occupancy depends on the elevation and forest cover, and is different for each site:$$z_i \sim {\rm Bernoulli}(\psi_i)$$

Because $\psi$ is a probability we use a logistic link, with an expression that’s essentially a logistic regression model. Switzerland is mountainous, with elevations ranging from 250 to 2750m; it’s likely that willow tit occupancy is highest at some middle elevation, so we include elevation^2 as a predictor to give a “humpy” curve rather than straight line. $${\rm logit}(\psi_i) = \beta_0 + \beta_{for} forest + \beta_{elev} elev + \beta_{elev2} elev^2$$

We won’t use covariates for the probability of detection, so the observation model is the same as before:$$y_i \sim {\rm Binomial}(n, pz_i)$$

### Priors

### The JAGS code

Work in progress