N-mixure model using nimble with parallel computing

nmix_nimble_parallel.Rmd

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

# N-mixture model

```{r setup}
library(AHMbook)
library(nimble)
library(coda)
library(parallel)
```

## Generation of simulated data

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

- Number of sites: 267
- Number of replicated surveys: 3
- Proportion of suitable sites : 1
- Expected abundance: 2
- Mean detection probability: 0.6
- Coefficient of site covariate 2 in abundance model: 1
- Coefficient of site covariate 3 in detection model: 1
- coefficient of survey covariate in detection model: -1

```{r gen_data}
set.seed(20220222)
nmix_data <- AHMbook::simNmix(
  nsites = 267,
  nvisits = 3,
  mean.lam = 2,
  mean.p = 0.6,
  beta2.lam = 1,
  beta3.p = 1,
  beta.p.survey = -1,
  show.plot = FALSE)
```

## View data

True abundance

```{r view_total}
sum(nmix_data$N)
```

Histogram

```{r view_data1}
barplot(hist(nmix_data$N, breaks = 0:(max(nmix_data$N) + 1),
             right = FALSE, plot = FALSE)$counts,
        names.arg = 0:max(nmix_data$N))
```

Histogram of observed number of individuals

```{r view_data2}
barplot(hist(nmix_data$C, breaks = 0:(max(nmix_data$C) + 1),
             right = FALSE, plot = FALSE)$counts,
        names.arg = 0:max(nmix_data$C))
```

Relationship between covariate 2 and N

```{r view_data3}
plot(nmix_data$site.cov[, 2], nmix_data$N)
```


## Model

```{r model}
mc <- nimbleCode({
  for (m in 1:M) {
    log(lambda[m]) <- beta[1] + beta[2] * X2[m]
    N[m] ~ dpois(lambda[m])
    for (r in 1:R) {
      logit(p[m, r]) <- beta[3] + beta[4] * X3[m] + beta[5] * Xs[m, r]
      y[m, r] ~ dbinom(size = N[m], prob = p[m, r])
    }
  }
  # priors
  for (i in 1:5) {
    beta[i] ~ dnorm(0, 1.0e-2)
  }
  # generated quantities
  Ntotal <- sum(N[1:M])
})
```

## Fitting

computing in parallel

see https://r-nimble.org/html_manual/cha-building-models.html#using-models-in-parallel

```{r fit}
nimble_const <- list(M = nmix_data$nsites,
                     R = nmix_data$nvisits)
nimble_data <- list(y = nmix_data$C,
                    X2 = nmix_data$site.cov[, 2],
                    X3 = nmix_data$site.cov[, 3],
                    Xs = nmix_data$survey.cov)
nimble_init <- function() {
  list(beta = runif(5, -2, 2),
       N = apply(nmix_data$C, 1, max))}
nimble_fit <- mclapply(1:3,
                       function(i) {
                         nimbleMCMC(code = mc,
                                    constants = nimble_const,
                                    data = nimble_data,
                                    inits = nimble_init,
                                    monitors = c("beta", "Ntotal"),
                                    niter = 8000,
                                    nburnin = 4000,
                                    nchains = 1,
                                    setSeed = i,
                                    samples = TRUE,
                                    samplesAsCodaMCMC = TRUE)},
                       mc.cores = 3)
```

### Results

```{r results}
fit <- mcmc.list(nimble_fit)
plot(fit)
```

```{r}
gelman.diag(fit)
summary(fit)
```