Site occupancy model using nimbleEcology

occ_nimbleEcology.Rmd

---
title: "R Notebook"
output: html_notebook
---

# Site occupancy model analyzed using nimbleEcology

```{r setup}
library(AHMbook)
library(nimbleEcology)
```

## Generation of simulated data

Simulated data are generated using the `simOcc` function of the `AHMbook` package.

- Number of sites: 300
- Number of observation for each site: 2
- Average occupancy probability: 0.6
- Average detection probability: 0.3
- Effecf of wind on detection : -3
- Standard deviation of site random effect: 0.5

```{r gen_site_occ_data}
set.seed(20220220)
occ_data <- AHMbook::simOcc(
  M = 300,                # Number of sites
  J = 2,                  # Number of temporal replicates
  mean.occupancy = 0.6,   # Mean occupancy prob.
  beta1 = -2,             # Main effect of elev. on occ.
  beta2 = 2,              # Main effect of forest on occ.
  beta3 = 0,              # Interaction of elev. and forest
  mean.detection = 0.3,   # Mean detection prob.
  time.effects = c(0, 0), # Time effect on det.
  alpha1 = 0,             # Main effect of elev. on det.
  alpha2 = -3,            # Main effect of wind speed on det.
  alpha3 = 0,             # Interaction of elev and wind
  sd.lp = 0.5,            # S.D. of random site effects
  b = 0,                  # Behavioral response
  show.plot = FALSE)
```

## View data

True number of occupied sites.

```{r view_site_occ_data1}
sum(occ_data$z)
```

Observed number of occupied sites.

```{r view_site_occ_data2}
sum(apply(occ_data$y, 1, max))
```

## Fitting to the occupancy model using nimbleEcology

### Model

```{r occ_nimble_model1}
mc <- nimbleCode({
  for (m in 1:M) {
    logit(psi[m]) <- beta[1] + beta[2] * elev[m] + beta[3] * forest[m]
    for (j in 1:J) {
      logit(p[m, j]) <- alpha[1] + alpha[2] * wind[m, j]
    }
    y[m, 1:J] ~ dOcc_v(psi[m], p[m, 1:J], len = J)
  }
  # priors
  for (i in 1:2) {
    alpha[i] ~ dnorm(0, 1.0e-2)
  }
  for (i in 1:3) {
    beta[i] ~ dnorm(0, 1.0e-2)
  }
  # derived values
  for (m in 1:M) {
    pr[m] <- step(sum(y[m, 1:J]) - 1) +
             equals(sum(y[m, 1:J]), 0) *
               (psi[m] * prod(1 - p[m, 1:J]) /
                  (psi[m] * prod(1 - p[m, 1:J]) + (1 - psi[m])))
    z[m] ~ dbinom(size = 1, prob = pr[m])
  }
  N <- sum(z[1:M])
})
```

### Fitting

```{r}
nimble_const <- list(M = occ_data$M, J = occ_data$J)
nimble_data <- list(y = occ_data$y,
                    elev = occ_data$elev, forest = occ_data$forest,
                    wind = occ_data$wind)
nimble_init <- function() {
  list(alpha = runif(2, -2, 2),
       beta = runif(3, -2, 2),
       z = apply(occ_data$y, 1, max))}
nimble_fit <- nimbleMCMC(code = mc,
                         constants = nimble_const,
                         data = nimble_data,
                         inits = nimble_init,
                         monitors = c("beta", "alpha", "N"),
                         niter = 8000,
                         nburnin = 4000,
                         nchains = 3,
                         setSeed = 1:3,
                         summary = TRUE,
                         samples = TRUE)
```

### Results

```{r}
nimble_fit$summary
```