Lecture 9b

class: center, middle, inverse, title-slide

.title[
# Lecture 9b
]
.subtitle[
## Static Occupancy Models
]
.author[
### Spring 2025
]

---

# Readings: 
 
> Kéry & Schaub 383-411

> Occupancy Estimation and Modeling : Inferring Patterns and Dynamics of Species Occurrence by Darryl I. MacKenzie, James D. Nichols, J. Andrew Royle, Kenneth H. Pollock, Larissa Bailey, and James E. Hines

---
## Estimating abundance?

#### For the past few weeks, we've been modeling abundance:

`$$\Large N_t \sim Poisson(\lambda)$$`

--
#### Occupancy is the probability a site is occupied

- Occupancy comes free from abundance
    + If `$\large N_i > 0$`, `$\large z_i = 1$`

+ If `$\large N_i = 0$`, `$\large z_i = 0$`

--
#### So even when occupancy is the state-variable, we are still modeling something related to abundance

---
## Estimating abundance?

#### Typical state model for occupancy

`$$\Large z_i \sim Bernoulli(\psi)$$`

--
#### If the expected abundance is `$\large \lambda$`, what is the probability `$\large N = 0$`?

`$$\Large Pr(N=0|\lambda)=\frac{\lambda^0e^-\lambda}{0!}=e^{-\lambda}$$`

--
#### If the expected abundance is `$\large \lambda$`, what is the probability `$\large N > 0$`?

`$$\Large 1 - Pr(N=0|\lambda)=1 -e^{-\lambda}$$`

--
#### So (if our assumptions are valid):

`$$\Large \psi = 1 - e^{-\lambda}$$`

---
## Why estimate occupancy?

#### Less effort

--
#### Historical data sets

--
#### More reliable when `$\large N$` very small

--
#### Occupancy = abundance (e.g., breeding territory)

--
#### Metapopulation dynamics

--
#### Distribution/range size

--
#### Disease dynamics

---
## Why not just use observed presence/absence?

#### As in all ecological studies, we rarely (if ever) observe the state process perfectly

--
#### In the case of occupancy, some sites will be occupied but we will fail to detect the species

- i.e., `$\large p < 1$`

--
#### Also possible (though generally more rare) that we record the species when it's not present (false positive)

- see Royle & Link 2006

--
#### Similar to N-mixture models, estimating `$\large p$` requires temporal replication

---
## Single-season (static) occupancy model

#### State-space formulation

- State-model

`$$\Large z_i \sim Bernoulli(\psi_i)$$`

`$$\Large logit(\psi_i) = \alpha_0 + \alpha_1x_i$$`

- Observation-model

`$$\Large y_{ik} \sim Bernoulli(z_i \times p_{ik})$$`

`$$\Large logit(p_{ik}) = \beta_0 + \beta_1x_{ik}$$`

---
## Assumptions of closed occupancy models

1) Binomial distribution is a true description of state/observation processes

2) Occupancy at each site is random and independent of occupancy at all other sites

3) Population is closed between surveys

4) Observers do not produce false positives

---
## Estimating `$\LARGE p$`

#### Imagine a single site surveyed 3 times:

- Assume site is closed across samples

- Assume constant `$\large p$`

`$$\LARGE y_i = [111]$$`

--
#### What is the likelihood of this observation?

--
P(occupied) x P(observed 3 times)

`$$\LARGE \psi p^3$$`

---
## Estimating `$\LARGE p$`

#### What about?

`$$\LARGE y_i = [011]$$`

--
P(occupied) x P(not observed 1 time) x P(observed 2 times)

`$$\LARGE \psi (1-p)p^2$$`
---
## Estimating `$\LARGE p$`

#### Is this scenario any different?

`$$\LARGE y_i = [101]$$`

--
No, the order of detection/non-detection does not matter.

P(occupied) x P(not observed 1 time) x P(observed 2 times)

`$$\LARGE \psi (1-p)p^2$$`

---
## Estimating `$\LARGE p$`

#### What about?

`$$\LARGE y_i = [000]$$`

--
 
Three 0's can arise from two different scenarios - the site is not occupied, or the site is occupied but we missed it 3 times.

P(not occupied) + P(occupied) x P(not observed 3 times)

`$$\LARGE (1 - \psi) + \psi (1-p)^3$$`

---

## In Class Exercise

Simulate 15 years of an occupancy system using the following information. Occupancy can change between years, but this change is not modeled:
 
- `set.seed(55)`

- `$\psi$` = .32

- nsites = 40

- Binomial detection process with 3 site visits per year.

-  `$E(p)$` = .4 + `$\eta_t$`

-  `$\eta_t \sim Normal(0, 1)$`

---

## In Class Exercise

``` r
set.seed(55)
nyears <- 15
psi <- .32
nsites <- 40
nvisits <- 3
mean_p <- .4
p <- plogis(mean_p + rnorm(nyears, 0, 1))

z <- y <-array(NA, c(nsites, nyears))
for(t in 1:nyears){
 for(i in 1:nsites){
 z[i,t] <- rbinom(1, 1, psi)
 y[i,t] <- rbinom(1, z[i,t]*nvisits, p[t])
 } #end i
} #end t
```

---

## In Class Exercise

What proportion of sites were occupied in each year?

``` r
colSums(z)/nsites
```

```
##  [1] 0.550 0.350 0.250 0.425 0.400 0.150 0.300 0.350 0.300 0.425 0.350 0.275
## [13] 0.325 0.175 0.350
```

How about observed occupancy?

``` r
colSums(y > 0)/nsites
```

```
##  [1] 0.525 0.050 0.200 0.275 0.400 0.150 0.275 0.350 0.300 0.375 0.325 0.175
## [13] 0.100 0.150 0.200
```

---
# Expected vs Realized Occupancy

---
## Nimble Model

``` r
occ.mod <- nimbleCode({
for(t in 1:nyears){
 for(i in 1:nsites){
 z[i,t] ~dbern(psi)
 y[i,t] ~dbinom(p = p[t], size = z[i,t]*nvisits)
 } #end i 
 
 logit(p[t]) <- p0+p1[t]
 p1[t] ~ dnorm(0, 1)
} #end t
 
psi ~ dbeta(1, 1)
p0 ~ dnorm(0, 1)
})
```

Quiz- is `p1` being treated as a random or fixed effect?

---
## Nimble Model

``` r
params <- c('p', 'p0', 'p1', 'psi')
n.c <- list(nsites = nsites, nvisits = nvisits, nyears = nyears)
n.d <- list(y =y)
n.i <- function()list(p0 = rnorm(1), p1 = rnorm(nyears), psi = runif(1))
```

Initial values for `$z$` can be tricky - we need to make sure we don't suggest a site is unoccupied if we had observations there 
 
--
One option is to start anywhere with y = 0 as unoccupied

``` r
z.init <- (y >0)*1
n.i <- function()list(p0 = rnorm(1), p1 = rnorm(nyears), psi = runif(1), z =z.init)
```

--
Alternatively, flip a coin for any sites where y = 0

``` r
z.init <- matrix(rbinom(nsites*nyears, 1, .5), ncol = nyears)
z.init[which(y > 0)] <- 1 
n.i <- function()list(p0 = rnorm(1), p1 = rnorm(nyears), psi = runif(1), z =z.init)
```

---
## Run the model in Nimble

``` r
occ_out <- nimbleMCMC(code = occ.mod,
 data = n.d,
 constants = n.c,
 inits = n.i(),
 monitors = params,
 thin = 1,
 niter = 5000,
 nburnin = 2500,
 nchains = 3,
 samplesAsCodaMCMC = TRUE
 )
```

---
## Check the model
<img src="9b_ExtraOcc_files/figure-html/unnamed-chunk-11-1.png" width="576" style="display: block; margin: auto;" />

---
## Compare with simulation
<img src="9b_ExtraOcc_files/figure-html/unnamed-chunk-12-1.png" width="576" style="display: block; margin: auto;" />