From b7393d34f2ccb86a1330133d90b5e411bdf59412 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?J=C3=BCrgen=20Niedballa?= Date: Mon, 26 Feb 2024 05:05:04 +0800 Subject: [PATCH] html_vignette instead of html_document --- vignettes/camtrapr5.html | 3484 ++++++++++++++------------------------ 1 file changed, 1298 insertions(+), 2186 deletions(-) diff --git a/vignettes/camtrapr5.html b/vignettes/camtrapr5.html index d9a5aa3..104c3ab 100644 --- a/vignettes/camtrapr5.html +++ b/vignettes/camtrapr5.html @@ -8,1294 +8,37 @@ + - - - - -5. Multi-species occupancy models - - - - - - - - - - - - - - - + - - - - - - - - - + + - - - - - - - - + -
- - - -
-
-
-
-
- -
- - - - - + -
library(camtrapR)
-library(purrr)
-library(DT)
-library(knitr)
-library(ggplot2)
-library(sf)
+
Code +
library(camtrapR)
+library(purrr)
+library(DT)
+library(knitr)
+library(ggplot2)
+library(sf)
+

Overview

Multi-species occupancy models are a powerful tool for combining @@ -4157,53 +3087,61 @@

Quick example

Data preparation

First load the camera trap table and create a camera operation matrix:

-
 data("camtraps")
-
- camop_no_problem <- cameraOperation(CTtable      = camtraps,
-                                     stationCol   = "Station",
-                                     setupCol     = "Setup_date",
-                                     retrievalCol = "Retrieval_date",
-                                     hasProblems  = FALSE,
-                                     dateFormat   = "dmy"
- )
+
Code +
 data("camtraps")
+
+ camop_no_problem <- cameraOperation(CTtable      = camtraps,
+                                     stationCol   = "Station",
+                                     setupCol     = "Setup_date",
+                                     retrievalCol = "Retrieval_date",
+                                     hasProblems  = FALSE,
+                                     dateFormat   = "dmy"
+ )
+

First, we load the record table and create a list of detection histories. While all species are included in this example, you may choose to subset the record table beforehand. It is important to note the includeEffort = TRUE argument, which means that the output for each species will be a list containing both the detection history and effort matrix.

-
 data("recordTableSample")
-
- # list of detection histories
- DetHist_list <- lapply(unique(recordTableSample$Species), FUN = function(x) {
-   detectionHistory(
-     recordTable         = recordTableSample,
-     camOp                = camop_no_problem,
-     stationCol           = "Station",
-     speciesCol           = "Species",
-     recordDateTimeCol    = "DateTimeOriginal",
-     species              = x,     # this gets modifies by lapply
-     occasionLength       = 7,
-     day1                 = "station",
-     datesAsOccasionNames = FALSE,
-     includeEffort        = TRUE,
-     scaleEffort          = FALSE,
-     timeZone             = "Asia/Kuala_Lumpur"
-   )}
- )
-
- # assign species names to the list items
- names(DetHist_list) <- unique(recordTableSample$Species)
-
- # note, DetHist_list is a list containing a list for each species
+
Code +
 data("recordTableSample")
+
+ # list of detection histories
+ DetHist_list <- lapply(unique(recordTableSample$Species), FUN = function(x) {
+   detectionHistory(
+     recordTable         = recordTableSample,
+     camOp                = camop_no_problem,
+     stationCol           = "Station",
+     speciesCol           = "Species",
+     recordDateTimeCol    = "DateTimeOriginal",
+     species              = x,     # this gets modifies by lapply
+     occasionLength       = 7,
+     day1                 = "station",
+     datesAsOccasionNames = FALSE,
+     includeEffort        = TRUE,
+     scaleEffort          = FALSE,
+     timeZone             = "Asia/Kuala_Lumpur"
+   )}
+ )
+
+ # assign species names to the list items
+ names(DetHist_list) <- unique(recordTableSample$Species)
+
+ # note, DetHist_list is a list containing a list for each species
+

Get the detection history of each species and put into a new list (thereby removing the effort matrix).

-
 ylist <- lapply(DetHist_list, FUN = function(x) x$detection_history)
+
Code +
 ylist <- lapply(DetHist_list, FUN = function(x) x$detection_history)
+

Create some made-up covariates (only for demonstration):

-
 sitecovs <- camtraps[, c(1:3)]
- sitecovs$elevation <- c(300, 500, 600)
- sitecovs[, c(2:4)] <- scale(sitecovs[,-1])   # scale numeric covariates
- sitecovs$fact <- factor(c("A", "A", "B"))    # categorical covariate
+
Code +
 sitecovs <- camtraps[, c(1:3)]
+ sitecovs$elevation <- c(300, 500, 600)
+ sitecovs[, c(2:4)] <- scale(sitecovs[,-1])   # scale numeric covariates
+ sitecovs$fact <- factor(c("A", "A", "B"))    # categorical covariate
+

Now we bundle the necessary data for communityModel in a list with 3 items:

    @@ -4212,9 +3150,11 @@

    Data preparation

  • obsCovs: (named) list with observation level covariates. Requires at least the effort matrix.
-
 data_list <- list(ylist    = ylist,
-                   siteCovs = sitecovs,
-                   obsCovs  = list(effort = DetHist_list[[1]]$effort))  # is identical for all species
+
Code +
 data_list <- list(ylist    = ylist,
+                   siteCovs = sitecovs,
+                   obsCovs  = list(effort = DetHist_list[[1]]$effort))  # is identical for all species
+

Now data_list is a list containing the detection histories, site covariates and occasion level covariates.

Please note that it is always necessary to provide an effort matrix. @@ -4238,15 +3178,17 @@

Fit model in JAGS

information for running a community occupancy model. Here we will create a JAGS model (the default). For a model in Nimble, see further below.

-
# text file to save the model
- modelfile1 <- tempfile(fileext = ".txt")
-
- mod.jags <- communityModel(data_list,
-                            occuCovs = list(fixed = "utm_y", ranef = "elevation"),
-                            detCovsObservation = list(fixed = "effort"),
-                            intercepts = list(det = "ranef", occu = "ranef"),
-                            modelFile = modelfile1)
-
## Wrote model to C:\Users\NIEDBA~1\AppData\Local\Temp\Rtmpqkffp1\file1fb862774ec.txt
+
Code +
# text file to save the model
+ modelfile1 <- tempfile(fileext = ".txt")
+
+ mod.jags <- communityModel(data_list,
+                            occuCovs = list(fixed = "utm_y", ranef = "elevation"),
+                            detCovsObservation = list(fixed = "effort"),
+                            intercepts = list(det = "ranef", occu = "ranef"),
+                            modelFile = modelfile1)
+
+
## Wrote model to C:\Users\NIEDBA~1\AppData\Local\Temp\Rtmp4UNIvN\file4c5879e733bc.txt

This call defined a model with a fixed effect of utm_y on occupancy probability of all species, and a species-specific (random) effect of elevation on occupancy probability. Furthermore, we specified effort as @@ -4259,7 +3201,9 @@

Fit model in JAGS

One can also see the model code using (not shown here for brevity):

There is a very basic summary method for commOccu objects:

-
 summary(mod.jags)
+
Code +
 summary(mod.jags)
+
## commOccu object for community occupancy model (in JAGS)
 ## 
 ## 5 species,  3 stations,  7 occasions
@@ -4275,54 +3219,70 @@ 
Fit model in JAGS

 ## Available site-occasion covariates:
 ##  effort

The model is fit with fit().

-
 fit.jags <- fit(mod.jags,
-                 n.iter = 5000,
-                 n.burnin = 2500,
-                 chains = 3)
+
Code +
 fit.jags <- fit(mod.jags,
+                 n.iter = 5000,
+                 n.burnin = 2500,
+                 chains = 3)
+

The output of the fit() method is an mcmc.list from the coda package.

Summarize posterior estimates:

-
fit_summary <- summary(fit.jags)
+
Code +
fit_summary <- summary(fit.jags)
+

Statistics of parameter estimates:

-
# Note, colors may not render correctly in dark themes in RStudio.
-DT::datatable(round(fit_summary$statistics, 3))
-
- +
Code +
# Note, colors may not render correctly in dark themes in RStudio.
+DT::datatable(round(fit_summary$statistics, 3))
+
+
+

Quantiles of parameter estimates:

-
DT::datatable(round(fit_summary$quantiles, 3))
-
- +
Code +
DT::datatable(round(fit_summary$quantiles, 3))
+
+
+

Marginal effect plots (response curves) can be plotted with plot_effects(). Argument submodel defines whether the output is for the detection or occupancy part of the model:

-
 plot_effects(mod.jags,
-              fit.jags,
-              submodel = "state")
+
Code +
 plot_effects(mod.jags,
+              fit.jags,
+              submodel = "state")
+
## $utm_y
-

+

## 
 ## $elevation
-

-
 plot_effects(mod.jags,
-              fit.jags,
-              submodel = "det")
+

+
Code +
 plot_effects(mod.jags,
+              fit.jags,
+              submodel = "det")
+
## $effort
-

+

Likewise, we can plot effect sizes for easier comparison between species and for easily checking significance with plot_coef():

-
 plot_coef(mod.jags,
-           fit.jags,
-           submodel = "state",
-           combine = T)
+
Code +
 plot_coef(mod.jags,
+           fit.jags,
+           submodel = "state",
+           combine = T)
+
## 'combine' and 'ordered' can't both be TRUE. Setting 'ordered = FALSE'
-

-
 plot_coef(mod.jags,
-           fit.jags,
-           submodel = "det")
+

+
Code +
 plot_coef(mod.jags,
+           fit.jags,
+           submodel = "det")
+
## $effort
-

+

There are further arguments for the significance levels, sorting species by parameter estimates, and for combining multiple plots.

@@ -4339,8 +3299,10 @@

Introduction

faster in some situations.

Using Nimble for community occupancy models requires the packages “nimble” and “nimbleEcology”:

-
library(nimble)
-library(nimbleEcology)
+
Code +
library(nimble)
+library(nimbleEcology)
+

Compilation and Rtools

@@ -4352,68 +3314,88 @@

Compilation and Rtools

see something like “C:/rtools40/usr/bin” or “C:/RBuildTools/4.2” somewhere in the output of the following line (for me it’s usually the first or second entry).

-
Sys.getenv("PATH")
+
Code +
Sys.getenv("PATH")
+

Fit models in Nimble

We can fit the same model in Nimble using the exact same functions as for the JAGS model above. Only set nimble = TRUE.

-
 modelfile2 <- tempfile(fileext = ".txt")
-
- mod.nimble <- communityModel(data_list,
-                              occuCovs = list(fixed = "utm_x", ranef = "utm_y"),
-                              detCovsObservation = list(fixed = "effort"),
-                              intercepts = list(det = "ranef", occu = "ranef"),
-                              modelFile = modelfile2,
-                              nimble = TRUE)      # set nimble = TRUE
-
## Wrote model to C:\Users\NIEDBA~1\AppData\Local\Temp\Rtmpqkffp1\file1fb87f3c1ab3.txt
+
Code +
 modelfile2 <- tempfile(fileext = ".txt")
+
+ mod.nimble <- communityModel(data_list,
+                              occuCovs = list(fixed = "utm_x", ranef = "utm_y"),
+                              detCovsObservation = list(fixed = "effort"),
+                              intercepts = list(det = "ranef", occu = "ranef"),
+                              modelFile = modelfile2,
+                              nimble = TRUE)      # set nimble = TRUE
+
+
## Wrote model to C:\Users\NIEDBA~1\AppData\Local\Temp\Rtmp4UNIvN\file4c5826ec4cd4.txt

It is possible to fit uncompiled models, but it is usually extremely slow and should only be done for checking model structure and expected outputs, using very few iterations.

-
 fit.nimble.uncomp <- fit(mod.nimble,
-                          n.iter = 10,
-                          chains = 1)
-
-# the notes and the error message above are harmless
+
Code +
 fit.nimble.uncomp <- fit(mod.nimble,
+                          n.iter = 10,
+                          chains = 1)
+
+# the notes and the error message above are harmless
+

To fit a compiled model, set compile = TRUE:

-
 fit.nimble.comp <- fit(mod.nimble,
-                        n.iter = 20000,
-                        n.burnin = 10000,
-                        chains = 3,
-                        compile = TRUE, 
-)
+
Code +
 fit.nimble.comp <- fit(mod.nimble,
+                        n.iter = 20000,
+                        n.burnin = 10000,
+                        chains = 3,
+                        compile = TRUE, 
+)
+

The output as an mcmc.list, just like for the JAGS model above, and can be worked with just the same:

-
 # parameter summary statistics (not shown)
- summary(fit.nimble.comp)
+
Code +
 # parameter summary statistics (not shown)
+ summary(fit.nimble.comp)
+

Response curves (= marginal effect plots):

-
 plot_effects(object = mod.nimble,
-              mcmc.list = fit.nimble.comp,
-              submodel = "state")
+
Code +
 plot_effects(object = mod.nimble,
+              mcmc.list = fit.nimble.comp,
+              submodel = "state")
+
## $utm_x
-

+

## 
 ## $utm_y
-

-
 plot_effects(object = mod.nimble,
-              mcmc.list = fit.nimble.comp,
-              submodel = "det")
+

+
Code +
 plot_effects(object = mod.nimble,
+              mcmc.list = fit.nimble.comp,
+              submodel = "det")
+
## $effort
-

+

Effect size plots:

-
 plot_coef(mod.nimble,
-           fit.nimble.comp,
-           submodel = "state",
-           combine = TRUE)
+
Code +
 plot_coef(mod.nimble,
+           fit.nimble.comp,
+           submodel = "state",
+           combine = TRUE)
+
## 'combine' and 'ordered' can't both be TRUE. Setting 'ordered = FALSE'
-

-
 plot_coef(mod.nimble,
-           fit.nimble.comp,
-           submodel = "det")
+

+
Code +
 plot_coef(mod.nimble,
+           fit.nimble.comp,
+           submodel = "det")
+
## $effort
-

+

Traceplots (not shown)

-
 plot(fit.nimble.comp)
+
Code +
 plot(fit.nimble.comp)
+

The estimates and marginal effects are not very pretty, and the model likely requires more iterations to achieve smooth estimates.

@@ -4430,77 +3412,94 @@

More complete example

Prepare data

First, install the AHMbook package if necessary:

-
install.packages("AHMbook")
-
library(AHMbook)
+
Code +
install.packages("AHMbook")
+
+
Code +
library(AHMbook)
+

Simulate community with simComm(). Here we simulate a community with 20 species, 30 camera trap stations and 5 occasions.

-
set.seed(10)
-com <- simComm(nspecies = 20, 
-                nsites = 30,
-                nreps = 5,
-                mean.psi = 0.25,     # community mean of psi over all species in community (probability scale)  (0.25 = -1.1 on logit scale)
-                sig.lpsi = 1,        # SD around intercept of logit(psi)
-                mean.p = 0.25,       # community mean of detection probability over all species (0.25 = -1.1 on logit scale)
-                sig.lp = 1,          # SD around intercept of logit(p)
-                mu.beta.lpsi = 0.5,  # mean (community) effect of habitat on logit(psi)
-                sig.beta.lpsi = 0.5, # SD around mean (community) effect of habitat on logit(psi)
-                mu.beta.lp = -0.5,   # mean (community) effect of wind on logit(p)
-                sig.beta.lp = 1)     # SD around mean (community) effect of wind on logit(p)
-

+
Code +
set.seed(10)
+com <- simComm(nspecies = 20, 
+                nsites = 30,
+                nreps = 5,
+                mean.psi = 0.25,     # community mean of psi over all species in community (probability scale)  (0.25 = -1.1 on logit scale)
+                sig.lpsi = 1,        # SD around intercept of logit(psi)
+                mean.p = 0.25,       # community mean of detection probability over all species (0.25 = -1.1 on logit scale)
+                sig.lp = 1,          # SD around intercept of logit(p)
+                mu.beta.lpsi = 0.5,  # mean (community) effect of habitat on logit(psi)
+                sig.beta.lpsi = 0.5, # SD around mean (community) effect of habitat on logit(psi)
+                mu.beta.lp = -0.5,   # mean (community) effect of wind on logit(p)
+                sig.beta.lp = 1)     # SD around mean (community) effect of wind on logit(p)
+
+

+
## Warning: Plotting of output failed
+##    the plotting window is too small.
+##    Try calling 'dev.new()' before the 'sim*' function.

We need to wrangle the data a little to convert the data to input for communityModel(). This is only needed here because the data are generated by AHMbook::simComm().

-
# get continuous site covariate "habitat"
-habitat <- com$habitat
-
-# create (fake) categorical covariate (3 levels)
-# setting "medium" as reference level
-categ <- factor(ifelse(habitat < -0.5, "low", ifelse(habitat < 0.5, "medium", "high")), levels = c("medium", "low", "high"))
-
-# get observation covariate
-wind <- com$wind
-
-# list of detection histories
-ylist <- purrr::array_branch(com$y.all, 3)
-
-# remove unobserved species  (we leave that to data augmentation)
-ylist <- ylist[sapply(ylist, sum) >= 1]
-
-# here, effort was constant at all stations and occasions. Effort is thus a matrix of all 1
-effort <- matrix(1, nrow = nrow(ylist[[1]]), ncol =  ncol(ylist[[1]]))
-
-# bundle data
-input_AHM <- list(ylist = ylist, 
-                  siteCovs = data.frame(habitat = habitat,
-                                        categ   = categ),
-                  obsCovs = list(effort = effort,
-                                 wind = wind))
+
Code +
# get continuous site covariate "habitat"
+habitat <- com$habitat
+
+# create (fake) categorical covariate (3 levels)
+# setting "medium" as reference level
+categ <- factor(ifelse(habitat < -0.5, "low", ifelse(habitat < 0.5, "medium", "high")), levels = c("medium", "low", "high"))
+
+# get observation covariate
+wind <- com$wind
+
+# list of detection histories
+ylist <- purrr::array_branch(com$y.all, 3)
+
+# remove unobserved species  (we leave that to data augmentation)
+ylist <- ylist[sapply(ylist, sum) >= 1]
+
+# here, effort was constant at all stations and occasions. Effort is thus a matrix of all 1
+effort <- matrix(1, nrow = nrow(ylist[[1]]), ncol =  ncol(ylist[[1]]))
+
+# bundle data
+input_AHM <- list(ylist = ylist, 
+                  siteCovs = data.frame(habitat = habitat,
+                                        categ   = categ),
+                  obsCovs = list(effort = effort,
+                                 wind = wind))
+

JAGS Model

Here’s a model in JAGS, with a species-specific effect of habitat on occupancy and wind on detection, and species-specific intercepts.

-
modelFile_jags <- tempfile(fileext = ".txt")
-
-model1_jags <- communityModel(
-  occuCovs = list(ranef = c("habitat")), 
-  detCovsObservation = list(ranef = "wind"),
-  intercepts = list(det = "ranef", occu = "ranef"),
-  data_list = input_AHM,
-  modelFile = modelFile_jags)
+
Code +
modelFile_jags <- tempfile(fileext = ".txt")
+
+model1_jags <- communityModel(
+  occuCovs = list(ranef = c("habitat")), 
+  detCovsObservation = list(ranef = "wind"),
+  intercepts = list(det = "ranef", occu = "ranef"),
+  data_list = input_AHM,
+  modelFile = modelFile_jags)
+
## Site covariates may not be scaled. Scaling is strongly recommended for numeric covariates.
-
## Wrote model to C:\Users\NIEDBA~1\AppData\Local\Temp\Rtmpqkffp1\file1fb821176852.txt
-
out_ahm_jags <- fit(model1_jags, 
-                    n.iter = 10000, 
-                    n.burnin = 5000,
-                    thin = 5,
-                    chains = 3,
-                    quiet = T
-)
+
## Wrote model to C:\Users\NIEDBA~1\AppData\Local\Temp\Rtmp4UNIvN\file4c5822331a9c.txt
+
Code +
out_ahm_jags <- fit(model1_jags, 
+                    n.iter = 10000, 
+                    n.burnin = 5000,
+                    thin = 5,
+                    chains = 3,
+                    quiet = T
+)
+

Summary statistics of model parameters:

-
DT::datatable(round(summary(out_ahm_jags)$statistics, 3))
-
- +
Code +
DT::datatable(round(summary(out_ahm_jags)$statistics, 3))
+
+
+

Since we used simulated data we can compare model estimates with the truth. For example, the true effect of habitat on occupancy was mean = 0.5 (SD = 0.5). See beta.ranef.cont.habitat.mean and @@ -4518,65 +3517,79 @@

JAGS Model

indicate no lack of fit for any species. This is not surprising because the model uses simulated data.

Plot response curves (marginal effects)

-
plot_resp_jags_occu <- plot_effects(model1_jags, 
-                                    out_ahm_jags)
-plot_resp_jags_occu
+
Code +
plot_resp_jags_occu <- plot_effects(model1_jags, 
+                                    out_ahm_jags)
+plot_resp_jags_occu
+
## $habitat
-

-
plot_resp_jags_det <- plot_effects(model1_jags, 
-                                   out_ahm_jags,
-                                   submodel = "det")
-plot_resp_jags_det
+

+
Code +
plot_resp_jags_det <- plot_effects(model1_jags, 
+                                   out_ahm_jags,
+                                   submodel = "det")
+plot_resp_jags_det
+
## $wind
-

+

Plot coefficient estimates:

-
plot_eff_jags_occu <- plot_coef(model1_jags, 
-                                out_ahm_jags)
-plot_eff_jags_occu
+
Code +
plot_eff_jags_occu <- plot_coef(model1_jags, 
+                                out_ahm_jags)
+plot_eff_jags_occu
+
## $habitat
-

-
plot_eff_jags_det <- plot_coef(model1_jags, 
-                               out_ahm_jags,
-                               submodel = "det")
-plot_eff_jags_det
+

+
Code +
plot_eff_jags_det <- plot_coef(model1_jags, 
+                               out_ahm_jags,
+                               submodel = "det")
+plot_eff_jags_det
+
## $wind
-

+

By default, species are sorted by effect size. Alternatively, you can have them sorted by species names. Set ordered = FALSE:

-
plot_eff_jags_occu2 <- plot_coef(model1_jags, 
-                                out_ahm_jags,
-                                ordered = FALSE)
-plot_eff_jags_occu2
+
Code +
plot_eff_jags_occu2 <- plot_coef(model1_jags, 
+                                out_ahm_jags,
+                                ordered = FALSE)
+plot_eff_jags_occu2
+
## $habitat
-

+

Nimble Model

Here’s the same model in Nimble:

-
modelFile_nimble <- tempfile(fileext = ".txt")
-
-model1_nimble <- communityModel(
-  occuCovs = list(ranef = c("habitat")),
-  detCovsObservation = list(ranef = "wind"),
-  intercepts = list(det = "ranef", occu = "ranef"),
-  data_list = input_AHM,
-  modelFile = modelFile_nimble,
-  nimble = TRUE)
+
Code +
modelFile_nimble <- tempfile(fileext = ".txt")
+
+model1_nimble <- communityModel(
+  occuCovs = list(ranef = c("habitat")),
+  detCovsObservation = list(ranef = "wind"),
+  intercepts = list(det = "ranef", occu = "ranef"),
+  data_list = input_AHM,
+  modelFile = modelFile_nimble,
+  nimble = TRUE)
+
## Site covariates may not be scaled. Scaling is strongly recommended for numeric covariates.
-
## Wrote model to C:\Users\NIEDBA~1\AppData\Local\Temp\Rtmpqkffp1\file1fb860083fb2.txt
+
## Wrote model to C:\Users\NIEDBA~1\AppData\Local\Temp\Rtmp4UNIvN\file4c5859875d92.txt

We use some more iterations, since these Nimble models tend to require more iterations than equivalent JAGS models to achieve good estimates (numbers are still too low for publication quality though).

-
# compiled model run
-out_ahm_nimble2 <- fit(model1_nimble, 
-                       n.iter =20000,
-                       n.burnin = 10000,
-                       thin = 5,
-                       chains = 3,
-                       compile = TRUE,
-                       quiet = T
-)
+
Code +
# compiled model run
+out_ahm_nimble2 <- fit(model1_nimble, 
+                       n.iter =20000,
+                       n.burnin = 10000,
+                       thin = 5,
+                       chains = 3,
+                       compile = TRUE,
+                       quiet = T
+)
+
## Defining model
##   [Note] Using 'y' (given within 'constants') as data.
## Building model
@@ -4596,31 +3609,41 @@

Nimble Model

## running chain 1...
## running chain 2...
## running chain 3...
-
DT::datatable(round(summary(out_ahm_nimble2)$statistics, 3))
+
Code +
DT::datatable(round(summary(out_ahm_nimble2)$statistics, 3))
+

Plot responses (marginal effects):

-
plot_resp_nimble_occu <- plot_effects(model1_nimble, 
-                                      out_ahm_nimble2)
-plot_resp_nimble_occu
+
Code +
plot_resp_nimble_occu <- plot_effects(model1_nimble, 
+                                      out_ahm_nimble2)
+plot_resp_nimble_occu
+
## $habitat
-

-
plot_resp_nimble_det <- plot_effects(model1_nimble, 
-                                     out_ahm_nimble2,
-                                     submodel = "det")
-plot_resp_nimble_det
+

+
Code +
plot_resp_nimble_det <- plot_effects(model1_nimble, 
+                                     out_ahm_nimble2,
+                                     submodel = "det")
+plot_resp_nimble_det
+
## $wind
-

+

Plot coefficient estimates:

-
plot_eff_nimble_occu <- plot_coef(model1_nimble, 
-                                  out_ahm_nimble2)
-plot_eff_nimble_occu
+
Code +
plot_eff_nimble_occu <- plot_coef(model1_nimble, 
+                                  out_ahm_nimble2)
+plot_eff_nimble_occu
+
## $habitat
-

-
plot_eff_nimble_det <- plot_coef(model1_nimble, 
-                                 out_ahm_nimble2,
-                                 submodel = "det")
-plot_eff_nimble_det
+

+
Code +
plot_eff_nimble_det <- plot_coef(model1_nimble, 
+                                 out_ahm_nimble2,
+                                 submodel = "det")
+plot_eff_nimble_det
+
## $wind
-

+

Missing values in detection histories

@@ -4632,69 +3655,81 @@

Missing values in detection histories

multi-species models (at least when fitting with Nimble). Let’s modify the above data set and introduce some NAs in the detection history and the effort matrix:

-
# introduce a few missing values in detection history of all species
-index_na <-  sample(1:length(ylist[[1]]), 10) #
-ylist_with_na <- lapply(ylist, FUN = function(x) {x[index_na] <- NA; x})
-
-
-# set the cells with NA in detection history to NA in effort too
-effort_with_na <- effort
-effort_with_na[index_na] <- NA
-
-# bundle model input
-input_AHM_with_NA <- list(ylist = ylist_with_na, 
-                          siteCovs = data.frame(habitat = habitat,
-                                                categ   = categ),
-                          obsCovs = list(effort = effort_with_na,
-                                         wind = wind))
+
Code +
# introduce a few missing values in detection history of all species
+index_na <-  sample(1:length(ylist[[1]]), 10) #
+ylist_with_na <- lapply(ylist, FUN = function(x) {x[index_na] <- NA; x})
+
+
+# set the cells with NA in detection history to NA in effort too
+effort_with_na <- effort
+effort_with_na[index_na] <- NA
+
+# bundle model input
+input_AHM_with_NA <- list(ylist = ylist_with_na, 
+                          siteCovs = data.frame(habitat = habitat,
+                                                categ   = categ),
+                          obsCovs = list(effort = effort_with_na,
+                                         wind = wind))
+

And fit models in JAGS and Nimble:

-
modelFile_jags_na <- tempfile(fileext = ".txt")
-
-model1_jags_na <- communityModel(
-  occuCovs = list(ranef = c("habitat")), 
-  detCovsObservation = list(ranef = "wind"),
-  intercepts = list(det = "ranef", occu = "ranef"),
-  data_list = input_AHM_with_NA,
-  modelFile = modelFile_jags_na)
+
Code +
modelFile_jags_na <- tempfile(fileext = ".txt")
+
+model1_jags_na <- communityModel(
+  occuCovs = list(ranef = c("habitat")), 
+  detCovsObservation = list(ranef = "wind"),
+  intercepts = list(det = "ranef", occu = "ranef"),
+  data_list = input_AHM_with_NA,
+  modelFile = modelFile_jags_na)
+
## Site covariates may not be scaled. Scaling is strongly recommended for numeric covariates.
-
## Wrote model to C:\Users\NIEDBA~1\AppData\Local\Temp\Rtmpqkffp1\file1fb843216970.txt
-
out_ahm_jags_na <- fit(model1_jags_na, 
-                    n.iter = 5000, 
-                    n.burnin = 2500,
-                    thin = 5,
-                    chains = 3,
-                    quiet = T
-)
+
## Wrote model to C:\Users\NIEDBA~1\AppData\Local\Temp\Rtmp4UNIvN\file4c5836a72214.txt
+
Code +
out_ahm_jags_na <- fit(model1_jags_na, 
+                    n.iter = 5000, 
+                    n.burnin = 2500,
+                    thin = 5,
+                    chains = 3,
+                    quiet = T
+)
+

Summary statistics of model parameters:

-
DT::datatable(round(summary(out_ahm_jags_na)$statistics, 3))
-
- +
Code +
DT::datatable(round(summary(out_ahm_jags_na)$statistics, 3))
+
+
+

As can be seen, the JAGS model handles NAs well and parameter estimates are very similar to the model without missing data above.

Nimble models can technically also be fit with NAs present, but there will be numerous warnings and residuals as well as Bayesian p-values cannot be computed.

-
modelFile_nimble_na <- tempfile(fileext = ".txt")
-
-model1_nimble_na <- communityModel(
-  occuCovs = list(ranef = c("habitat")), 
-  detCovsObservation = list(ranef = "wind"),
-  intercepts = list(det = "ranef", occu = "ranef"),
-  data_list = input_AHM_with_NA,
-  modelFile = modelFile_nimble_na,
-  nimble = T)
+
Code +
modelFile_nimble_na <- tempfile(fileext = ".txt")
+
+model1_nimble_na <- communityModel(
+  occuCovs = list(ranef = c("habitat")), 
+  detCovsObservation = list(ranef = "wind"),
+  intercepts = list(det = "ranef", occu = "ranef"),
+  data_list = input_AHM_with_NA,
+  modelFile = modelFile_nimble_na,
+  nimble = T)
+
## Site covariates may not be scaled. Scaling is strongly recommended for numeric covariates.
-
## Wrote model to C:\Users\NIEDBA~1\AppData\Local\Temp\Rtmpqkffp1\file1fb838e4c60.txt
+
## Wrote model to C:\Users\NIEDBA~1\AppData\Local\Temp\Rtmp4UNIvN\file4c582b2b1f46.txt

When fitting the model there will be numerous warnings about values of data nodes y[,,] and determinisitc nodes res[,], R2[] and R3 being NA (not shown here).

-
out_ahm_nimble_na <- fit(model1_nimble_na, 
-                    n.iter = 5000, 
-                    n.burnin = 2500,
-                    thin = 5,
-                    chains = 3,
-                    quiet = T
-)
+
Code +
out_ahm_nimble_na <- fit(model1_nimble_na, 
+                    n.iter = 5000, 
+                    n.burnin = 2500,
+                    thin = 5,
+                    chains = 3,
+                    quiet = T
+)
+
## Defining model
##   [Note] Using 'y' (given within 'constants') as data.
## Building model
@@ -4716,14 +3751,18 @@

Missing values in detection histories

## running chain 3...

Parameter summaries cannot be computed and will return an error (not shown here).

-
DT::datatable(round(summary(out_ahm_nimble2)$statistics, 3))
+
Code +
DT::datatable(round(summary(out_ahm_nimble2)$statistics, 3))
+

To alleviate this, it is necessary (for Nimble models) to convert NAs in the detection history to 0, e.g. with:

-
# replace NAs with 0 in the ylist component of input_AHM (for all species at once)
-input_AHM_with_NA$ylist <- lapply(input_AHM_with_NA$ylist, FUN = function(x) {
-  x[is.na(x)] <- 0
-  x
-})
+
Code +
# replace NAs with 0 in the ylist component of input_AHM (for all species at once)
+input_AHM_with_NA$ylist <- lapply(input_AHM_with_NA$ylist, FUN = function(x) {
+  x[is.na(x)] <- 0
+  x
+})
+

After converting NAs to 0, Nimble models can be run without warnings and summaries can be calculated.

To prevent the additional 0s from affecting the estimation of @@ -4772,52 +3811,58 @@

Data augmentation

We know there are 20 species in the simulated community example. We can inform the model that the community is known to consist of 20 species via: augmentation = c(maxknown = 20):

-
modelFile_jags_aug1 <- tempfile(fileext = ".txt")
-
-model1_jags_aug1 <- communityModel(
-  occuCovs           = list(ranef = c("habitat")), 
-  detCovsObservation = list(ranef = "wind"),
-  intercepts         = list(det = "ranef", occu = "ranef"),
-  data_list          = input_AHM,
-  augmentation       = c(maxknown = 20),
-  modelFile          = modelFile_jags_aug1)
-
-# short model run for demonstration
-out_ahm_jags_aug1 <- fit(model1_jags_aug1, 
-                         n.iter = 500, 
-                         n.burnin = 250,
-                         thin = 5,
-                         chains = 3,
-                         quiet = T
-)
+
Code +
modelFile_jags_aug1 <- tempfile(fileext = ".txt")
+
+model1_jags_aug1 <- communityModel(
+  occuCovs           = list(ranef = c("habitat")), 
+  detCovsObservation = list(ranef = "wind"),
+  intercepts         = list(det = "ranef", occu = "ranef"),
+  data_list          = input_AHM,
+  augmentation       = c(maxknown = 20),
+  modelFile          = modelFile_jags_aug1)
+
+# short model run for demonstration
+out_ahm_jags_aug1 <- fit(model1_jags_aug1, 
+                         n.iter = 500, 
+                         n.burnin = 250,
+                         thin = 5,
+                         chains = 3,
+                         quiet = T
+)
+

Here is another example with fully open data augmentation (metacommunity size is unknown). The model will estimate the probability for each species to be part of the metacommunity. Fully open data augmentation is specified with augmentation = c(full = 30) (for a maximum potential community size of 30 species in this example).

-
modelFile_jags_aug2 <- tempfile(fileext = ".txt")
-
-model1_jags_aug2 <- communityModel(
-  occuCovs = list(ranef = c("habitat")), 
-  detCovsObservation = list(ranef = "wind"),
-  intercepts = list(det = "ranef", occu = "ranef"),
-  data_list = input_AHM,
-  augmentation = c(full = 30),
-  modelFile = modelFile_jags_aug2)
-
-# short model run for demonstration
-out_ahm_jags_aug2 <- fit(model1_jags_aug2, 
-                         n.iter = 500, 
-                         n.burnin = 250,
-                         thin = 5,
-                         chains = 3,
-                         quiet = T
-)
+
Code +
modelFile_jags_aug2 <- tempfile(fileext = ".txt")
+
+model1_jags_aug2 <- communityModel(
+  occuCovs = list(ranef = c("habitat")), 
+  detCovsObservation = list(ranef = "wind"),
+  intercepts = list(det = "ranef", occu = "ranef"),
+  data_list = input_AHM,
+  augmentation = c(full = 30),
+  modelFile = modelFile_jags_aug2)
+
+# short model run for demonstration
+out_ahm_jags_aug2 <- fit(model1_jags_aug2, 
+                         n.iter = 500, 
+                         n.burnin = 250,
+                         thin = 5,
+                         chains = 3,
+                         quiet = T
+)
+

See parameters “omega” and “w[1]” to “w[30]” in the model output:

-
DT::datatable(round(summary(out_ahm_jags_aug2)$statistics, 3))
-
- +
Code +
DT::datatable(round(summary(out_ahm_jags_aug2)$statistics, 3))
+
+
+

Higher numbers of species in data augmentation increase model run time.

@@ -4833,32 +3878,38 @@

Site covariates

of species.

Here is an example with a categorical site covariate “categ” on occupancy:

-
modelFile_jags_categ1 <- tempfile(fileext = ".txt")
-
-model_jags_categ1 <- communityModel(
-  occuCovs = list(ranef = c("categ")),
-  detCovsObservation = list(ranef = "wind"),
-  intercepts = list(det = "ranef", occu = "ranef"),
-  data_list = input_AHM,
-  modelFile = modelFile_jags_categ1)
-
-# short model run for demonstration
-out_ahm_jags_categ1 <- fit(model_jags_categ1, 
-                         n.iter = 500, 
-                         n.burnin = 250,
-                         thin = 5,
-                         chains = 3,
-                         quiet = T
-)
+
Code +
modelFile_jags_categ1 <- tempfile(fileext = ".txt")
+
+model_jags_categ1 <- communityModel(
+  occuCovs = list(ranef = c("categ")),
+  detCovsObservation = list(ranef = "wind"),
+  intercepts = list(det = "ranef", occu = "ranef"),
+  data_list = input_AHM,
+  modelFile = modelFile_jags_categ1)
+
+# short model run for demonstration
+out_ahm_jags_categ1 <- fit(model_jags_categ1, 
+                         n.iter = 500, 
+                         n.burnin = 250,
+                         thin = 5,
+                         chains = 3,
+                         quiet = T
+)
+

Plot the results:

-
plot_effects(object = model_jags_categ1,
-             mcmc.list = out_ahm_jags_categ1)
+
Code +
plot_effects(object = model_jags_categ1,
+             mcmc.list = out_ahm_jags_categ1)
+
## $categ
-

-
plot_coef(object = model_jags_categ1,
-          mcmc.list = out_ahm_jags_categ1)
+

+
Code +
plot_coef(object = model_jags_categ1,
+          mcmc.list = out_ahm_jags_categ1)
+
## $categ
-

+

In the table below, see the arguments beta.ranef.categ.categ. The species effects have two indices: [species, factorLevel]. Accordingly, there are three estimates @@ -4866,9 +3917,11 @@

Site covariates

for the first factor level are 0, since the intercept refers to the first factor level (i.e. reference level, which we defined as “medium” above).

-
DT::datatable(round(summary(out_ahm_jags_categ1)$statistics, 3))
-
- +
Code +
DT::datatable(round(summary(out_ahm_jags_categ1)$statistics, 3))
+
+
+

Observation-level covariates

@@ -4876,47 +3929,57 @@

Observation-level covariates

supported in JAGS, not Nimble. Observation-level covariates are matrices and thus can’t be factors. Thus they need to be defined as character matrices.

-
# Create example categorical observation covariate
-input_AHM$obsCovs$wind_categ <- ifelse(input_AHM$obsCovs$wind > 0, "windy", "calm")
-
-# result is a character matrix:
-str(input_AHM$obsCovs$wind_categ)
+
Code +
# Create example categorical observation covariate
+input_AHM$obsCovs$wind_categ <- ifelse(input_AHM$obsCovs$wind > 0, "windy", "calm")
+
+# result is a character matrix:
+str(input_AHM$obsCovs$wind_categ)
+
##  chr [1:30, 1:5] "calm" "calm" "windy" "windy" "calm" "calm" "windy" "calm" "calm" "calm" "windy" "calm" "calm" "windy" "calm" "calm" "windy" "calm" "calm" ...

Model example:

-
modelFile_jags_categ2 <- tempfile(fileext = ".txt")
-
-model_jags_categ2 <- communityModel(
-  occuCovs = list(fixed = c("habitat")), 
-  detCovsObservation = list(ranef = c("wind_categ")),     # "wind_categ" now
-  intercepts = list(det = "ranef", occu = "ranef"),
-  data_list = input_AHM,
-  modelFile = modelFile_jags_categ2)
-
-out_ahm_jags_categ2 <- fit(model_jags_categ2, 
-                         n.iter = 500, 
-                         n.burnin = 250,
-                         thin = 5,
-                         chains = 3,
-                         quiet = T,
-                         compile = FALSE
-)
+
Code +
modelFile_jags_categ2 <- tempfile(fileext = ".txt")
+
+model_jags_categ2 <- communityModel(
+  occuCovs = list(fixed = c("habitat")), 
+  detCovsObservation = list(ranef = c("wind_categ")),     # "wind_categ" now
+  intercepts = list(det = "ranef", occu = "ranef"),
+  data_list = input_AHM,
+  modelFile = modelFile_jags_categ2)
+
+out_ahm_jags_categ2 <- fit(model_jags_categ2, 
+                         n.iter = 500, 
+                         n.burnin = 250,
+                         thin = 5,
+                         chains = 3,
+                         quiet = T,
+                         compile = FALSE
+)
+

Model estimates. As above, the estimates for the first factor level are 0. The relevant parameters here are called alpha.obs.ranef.categ.wind_categ.

-
DT::datatable(round(summary(out_ahm_jags_categ2)$statistics, 3))
-
- +
Code +
DT::datatable(round(summary(out_ahm_jags_categ2)$statistics, 3))
+
+
+

Plot the results:

-
plot_effects(object = model_jags_categ2, 
-             mcmc.list = out_ahm_jags_categ2,
-             submodel = "det")
+
Code +
plot_effects(object = model_jags_categ2, 
+             mcmc.list = out_ahm_jags_categ2,
+             submodel = "det")
+
## $wind_categ
-

-
plot_coef(object = model_jags_categ2, 
-          mcmc.list = out_ahm_jags_categ2,
-          submodel = "det")
+

+
Code +
plot_coef(object = model_jags_categ2, 
+          mcmc.list = out_ahm_jags_categ2,
+          submodel = "det")
+
## $wind_categ
-

+

@@ -4931,30 +3994,34 @@

Random effects other than species

detCovs(ranef = c(“covariate1|covariate2”))

Here is an example, where the effect of “habitat” differs between levels of “categ”.

-
modelFile_jags_ranef <- tempfile(fileext = ".txt")
-
-model1_jags_ranef <- communityModel(
-  occuCovs = list(ranef = c("habitat|categ")), 
-  detCovsObservation = list(ranef = "wind"),
-  intercepts = list(det = "ranef", occu = "ranef"),
-  data_list = input_AHM,
-  modelFile = modelFile_jags_ranef)
-
-# short model run for demonstration
-out_ahm_jags_ranef <- fit(model1_jags_ranef, 
-                         n.iter = 500, 
-                         n.burnin = 250,
-                         thin = 5,
-                         chains = 3,
-                         quiet = T
-)
+
Code +
modelFile_jags_ranef <- tempfile(fileext = ".txt")
+
+model1_jags_ranef <- communityModel(
+  occuCovs = list(ranef = c("habitat|categ")), 
+  detCovsObservation = list(ranef = "wind"),
+  intercepts = list(det = "ranef", occu = "ranef"),
+  data_list = input_AHM,
+  modelFile = modelFile_jags_ranef)
+
+# short model run for demonstration
+out_ahm_jags_ranef <- fit(model1_jags_ranef, 
+                         n.iter = 500, 
+                         n.burnin = 250,
+                         thin = 5,
+                         chains = 3,
+                         quiet = T
+)
+

Search for the term “habitat_categ” in the table below. There are 2 estimates (for stations with “A” and “B”), plus the mean and sigma of the random effect (this is only for demonstration, 2 factor levels is of course not enough for applying random effects).

-
DT::datatable(round(summary(out_ahm_jags_ranef)$statistics, 3))
-
- +
Code +
DT::datatable(round(summary(out_ahm_jags_ranef)$statistics, 3))
+
+
+

Currently random effects other than species can not be plotted with plot_effects and plot_coef.

@@ -4966,29 +4033,33 @@

Nested random effects

occuCovs(ranef = c(“covariate1|covariate2+Species”))

Here is an example, with independent random effect of “categ” on “habitat” for each species:

-
modelFile_jags_nested <- tempfile(fileext = ".txt")
-
-model1_jags_nested <- communityModel(
-  occuCovs = list(ranef = c("habitat|categ+Species")), 
-  detCovsObservation = list(ranef = "wind"),
-  intercepts = list(det = "ranef", occu = "ranef"),
-  data_list = input_AHM,
-  modelFile = modelFile_jags_nested)
-
-# short model run for demonstration
-out_ahm_jags_nested <- fit(model1_jags_nested, 
-                         n.iter = 500, 
-                         n.burnin = 250,
-                         thin = 5,
-                         chains = 3,
-                         quiet = T
-)
+
Code +
modelFile_jags_nested <- tempfile(fileext = ".txt")
+
+model1_jags_nested <- communityModel(
+  occuCovs = list(ranef = c("habitat|categ+Species")), 
+  detCovsObservation = list(ranef = "wind"),
+  intercepts = list(det = "ranef", occu = "ranef"),
+  data_list = input_AHM,
+  modelFile = modelFile_jags_nested)
+
+# short model run for demonstration
+out_ahm_jags_nested <- fit(model1_jags_nested, 
+                         n.iter = 500, 
+                         n.burnin = 250,
+                         thin = 5,
+                         chains = 3,
+                         quiet = T
+)
+

Search for the term “habitat_categ_Species” in the table below. Depending on the number of species and factor levels The number of parameter estimates is very high when using nested random effects.

-
DT::datatable(round(summary(out_ahm_jags_nested)$statistics, 3))
-
- +
Code +
DT::datatable(round(summary(out_ahm_jags_nested)$statistics, 3))
+
+
+

Currently nested random effects can not be plotted with plot_effects and plot_coef.

@@ -5008,30 +4079,34 @@

Species and site random effect on detection

Site is defined as the sampling location (e.g. camera trap station in this context).

Below is an example model:

-
modelFile_jags_species_station_ran <- tempfile(fileext = ".txt")
-
-model1_jags_species_station_ran <- communityModel(
-  occuCovs = list(ranef = c("habitat")), 
-  detCovsObservation = list(ranef = "wind"),
-  intercepts = list(det = "ranef", occu = "ranef"),
-  speciesSiteRandomEffect = list(det = TRUE),
-  data_list = input_AHM,
-  modelFile = modelFile_jags_species_station_ran)
-
-# short model run for demonstration
-out_ahm_jags_species_station_ran <- fit(model1_jags_species_station_ran, 
-                                        n.iter = 500, 
-                                        n.burnin = 250,
-                                        thin = 5,
-                                        chains = 3,
-                                        quiet = T
-)
+
Code +
modelFile_jags_species_station_ran <- tempfile(fileext = ".txt")
+
+model1_jags_species_station_ran <- communityModel(
+  occuCovs = list(ranef = c("habitat")), 
+  detCovsObservation = list(ranef = "wind"),
+  intercepts = list(det = "ranef", occu = "ranef"),
+  speciesSiteRandomEffect = list(det = TRUE),
+  data_list = input_AHM,
+  modelFile = modelFile_jags_species_station_ran)
+
+# short model run for demonstration
+out_ahm_jags_species_station_ran <- fit(model1_jags_species_station_ran, 
+                                        n.iter = 500, 
+                                        n.burnin = 250,
+                                        thin = 5,
+                                        chains = 3,
+                                        quiet = T
+)
+

The model returns the standard deviation of this random effect as alpha.speciessite.ranef.sigma, which can be found in the table below.

-
DT::datatable(round(summary(out_ahm_jags_species_station_ran)$statistics, 3))
-
- +
Code +
DT::datatable(round(summary(out_ahm_jags_species_station_ran)$statistics, 3))
+
+
+

Quadratic covariate effects

@@ -5042,37 +4117,43 @@

Quadratic covariate effects

Covariates should be scaled to prior to squaring (“habitat” in this examples is already scaled).

Create quadratic covariate:

-
input_AHM$siteCovs$habitat_squared <- input_AHM$siteCovs$habitat ^2
+
Code +
input_AHM$siteCovs$habitat_squared <- input_AHM$siteCovs$habitat ^2
+

Create and fit model:

-
modelFile_jags_quadratic <- tempfile(fileext = ".txt")
-
-
-model1_jags_quadratic <- communityModel(
-  occuCovs = list(ranef = c("habitat", "habitat_squared")), 
-  detCovsObservation = list(ranef = "wind"),
-  intercepts = list(det = "ranef", occu = "ranef"),
-  data_list = input_AHM,
-  modelFile = modelFile_jags_quadratic,
-  keyword_quadratic = "_squared")   # this is the default, added for clarity
-
-# short model run for demonstration
-out_ahm_jags_quadratic<- fit(model1_jags_quadratic, 
-                         n.iter = 10000, 
-                         n.burnin = 5000,
-                         thin = 5,
-                         chains = 3,
-                         quiet = T
-)
+
Code +
modelFile_jags_quadratic <- tempfile(fileext = ".txt")
+
+
+model1_jags_quadratic <- communityModel(
+  occuCovs = list(ranef = c("habitat", "habitat_squared")), 
+  detCovsObservation = list(ranef = "wind"),
+  intercepts = list(det = "ranef", occu = "ranef"),
+  data_list = input_AHM,
+  modelFile = modelFile_jags_quadratic,
+  keyword_quadratic = "_squared")   # this is the default, added for clarity
+
+# short model run for demonstration
+out_ahm_jags_quadratic<- fit(model1_jags_quadratic, 
+                         n.iter = 10000, 
+                         n.burnin = 5000,
+                         thin = 5,
+                         chains = 3,
+                         quiet = T
+)
+

Thanks to to keyword “_squared” in the covariate name, plot_effects knows that “habitat” and “habitat_squared” are the same covariate and will combine the effects of “habitat” and “habitat_squared” in response curves. If a different keyword was used, set the parameter “keyword_quadratic” accordingly (it defaults to “_squared”).

-
plot_effects(model1_jags_quadratic, 
-             out_ahm_jags_quadratic)
+
Code +
plot_effects(model1_jags_quadratic, 
+             out_ahm_jags_quadratic)
+
## $habitat
-

+

The simulated data have no quadratic relationship with covariates, hence the quadratic effects are rather weak in this example.

@@ -5080,25 +4161,31 @@

Quadratic covariate effects

WAIC

Nimble models can return WAIC. Set argument fit(..., WAIC = TRUE).

-
out_ahm_nimble3 <- fit(model1_nimble, 
-                       n.iter = 500,
-                       n.burnin = 250,
-                       thin = 5,
-                       chains = 3,
-                       compile = TRUE,
-                       quiet = F, 
-                       WAIC = TRUE
-)
+
Code +
out_ahm_nimble3 <- fit(model1_nimble, 
+                       n.iter = 500,
+                       n.burnin = 250,
+                       thin = 5,
+                       chains = 3,
+                       compile = TRUE,
+                       quiet = F, 
+                       WAIC = TRUE
+)
+

The message about “individual pWAIC values that are greater than 0.4” in the example is likely due to the short model run.

If WAIC = TRUE, the output will be a list with 2 elements. The first (“samples”) is the familiar mcmc.list, the second (“WAIC”) is a waicList.

mcmc.list is accessible with $samples:

-
# not printed here
-summary(out_ahm_nimble3$samples)$statistics
+
Code +
# not printed here
+summary(out_ahm_nimble3$samples)$statistics
+

WAIC:

-
out_ahm_nimble3$WAIC
+
Code +
out_ahm_nimble3$WAIC
+
## nimbleList object of type waicNimbleList
 ## Field "WAIC":
 ## [1] 1071.5
@@ -5125,36 +4212,40 @@ 
Richness estimates for categorical covariates

 can be used to calculate separate species richness estimates for
 different levels of categorical site covariates. If it is not defined,
 there will only be global and by-station species richness estimates.

-
modelFile_jags_richness <- tempfile(fileext = ".txt")
-
-model1_jags_richness <- communityModel(
-  occuCovs = list(ranef = c("habitat")), 
-  detCovsObservation = list(ranef = "wind"),
-  intercepts = list(det = "ranef", occu = "ranef"),
-  data_list = input_AHM,
-  richnessCategories = "categ",
-  modelFile = modelFile_jags_richness)
-
-# short model run for demonstration
-out_ahm_jags_rich <- fit(model1_jags_richness, 
-                         n.iter = 500, 
-                         n.burnin = 250,
-                         thin = 5,
-                         chains = 3,
-                         quiet = T
-)

+
Code
+
modelFile_jags_richness <- tempfile(fileext = ".txt")
+
+model1_jags_richness <- communityModel(
+  occuCovs = list(ranef = c("habitat")), 
+  detCovsObservation = list(ranef = "wind"),
+  intercepts = list(det = "ranef", occu = "ranef"),
+  data_list = input_AHM,
+  richnessCategories = "categ",
+  modelFile = modelFile_jags_richness)
+
+# short model run for demonstration
+out_ahm_jags_rich <- fit(model1_jags_richness, 
+                         n.iter = 500, 
+                         n.burnin = 250,
+                         thin = 5,
+                         chains = 3,
+                         quiet = T
+)

+

 
Now search for “Nspecies” in the table below. You will see entries
 for three different parameters.

 
    
 
  • “Nspecies” is overall species richness,
    • 
 
  • “Nspecies_by_category[1]”, “…[2]” and “…[3]” are for the stations in
-category “medium” and “low”, “high” respectively.
    • 
+category “medium”, “low”, and “high” respectively.
 
  • “Nspecies_by_station[1]” to “…[30]” are species richness estimates
 for each sampling station.
    • 
 

-
DT::datatable(round(summary(out_ahm_jags_rich)$statistics, 3))

-

-
+
Code
+
DT::datatable(round(summary(out_ahm_jags_rich)$statistics, 3))

+

+

+
 
 
 

@@ -5163,72 +4254,94 @@ 
Spatial predictions

 created with communityModel, their model fits and covariate
 rasters.

 
This requires the terra package.

-
  library(terra)

+
Code
+
  library(terra)

+

 
For demonstration we will use a model that has both a continuous and
 a categorical covariate.

-
modelFile_jags_categ3 <- tempfile(fileext = ".txt")
-
-model_jags_categ3 <- communityModel(
-  occuCovs = list(ranef = c("habitat", 
-                            "categ")),
-  detCovsObservation = list(ranef = "wind"),
-  intercepts = list(det = "ranef", occu = "ranef"),
-  data_list = input_AHM,
-  modelFile = modelFile_jags_categ3)
-
-# short model run for demonstration
-out_ahm_jags_categ3 <- fit(model_jags_categ3, 
-                         n.iter = 500, 
-                         n.burnin = 250,
-                         thin = 5,
-                         chains = 3,
-                         quiet = T
-)

+
Code
+
modelFile_jags_categ3 <- tempfile(fileext = ".txt")
+
+model_jags_categ3 <- communityModel(
+  occuCovs = list(ranef = c("habitat", 
+                            "categ")),
+  detCovsObservation = list(ranef = "wind"),
+  intercepts = list(det = "ranef", occu = "ranef"),
+  data_list = input_AHM,
+  modelFile = modelFile_jags_categ3)
+
+# short model run for demonstration
+out_ahm_jags_categ3 <- fit(model_jags_categ3, 
+                         n.iter = 500, 
+                         n.burnin = 250,
+                         thin = 5,
+                         chains = 3,
+                         quiet = T
+)

+

 
Create a habitat raster. We’ll use the volcano
 raster for demonstration. The coordinates are not geographic though
 since this is only for demonstration.

-
data("volcano")
-r_volcano <- scale(rast(volcano))

+
Code
+
data("volcano")
+r_volcano <- scale(rast(volcano))

+

 
Create a categorical raster. Since it is not
 possible to have categorical raster (with data type “factor”), the
 raster must consist of integer numbers that refer to the factor levels
 (as defined in the model input).

-
levels(input_AHM$siteCovs$categ)    # show the order of factor levels, first one is reference

+
Code
+
levels(input_AHM$siteCovs$categ)    # show the order of factor levels, first one is reference

+

 
## [1] "medium" "low"    "high"

-
raster_categ <- rast(x = matrix(rep(c(1,2, 3), each = ncell(r_volcano) / 3), 
-                                  ncol = ncol(r_volcano),
-                                  nrow = nrow(r_volcano)))

-
# add raster to list
-stack_prediction <- rast(list(habitat = r_volcano,
-                                 categ = raster_categ))
-
-plot(stack_prediction, asp = 1)

-

-
# raster categ, from west to east: 1 = medium, 2 = low, 3 = high

+
Code
+
raster_categ <- rast(x = matrix(rep(c(1,2, 3), each = ncell(r_volcano) / 3), 
+                                  ncol = ncol(r_volcano),
+                                  nrow = nrow(r_volcano)))

+

+
Code
+
# add raster to list
+stack_prediction <- rast(list(habitat = r_volcano,
+                                 categ = raster_categ))
+
+plot(stack_prediction, asp = 1)

+

+

+
Code
+
# raster categ, from west to east: 1 = medium, 2 = low, 3 = high

+

 
Species occupancy maps:

-
# species occupancy estimates
-  predictions_psi <- predict(object    = model_jags_categ3, 
-                             mcmc.list = out_ahm_jags_categ3,
-                             x         = stack_prediction,
-                             type      = "psi",
-                             draws     = 1000)

+
Code
+
# species occupancy estimates
+  predictions_psi <- predict(object    = model_jags_categ3, 
+                             mcmc.list = out_ahm_jags_categ3,
+                             x         = stack_prediction,
+                             type      = "psi",
+                             draws     = 1000)

+

 
## draws (1000) is greater than the number of available samples (150).

-
  plot(predictions_psi$mean, zlim = c(0,1), 
-       col = hcl.colors(100), 
-       maxnl = 9,   # plotting only the first 9 species for space reasons
-       asp = 1)  

-

+
Code
+
  plot(predictions_psi$mean, zlim = c(0,1), 
+       col = hcl.colors(100), 
+       maxnl = 9,   # plotting only the first 9 species for space reasons
+       asp = 1)  

+

+

 
Species richness maps:

-
# species richness estimates
-  predictions_rich <- predict(object   = model_jags_categ3, 
-                             mcmc.list = out_ahm_jags_categ3,
-                             x         = stack_prediction,
-                             type      = "richness")
-  
-  plot(predictions_rich, col = hcl.colors(100), asp = 1)

-

-
  # mean = mean species richness estimate
-  # sd = standard deviation of richness estimate

+
Code
+
# species richness estimates
+  predictions_rich <- predict(object   = model_jags_categ3, 
+                             mcmc.list = out_ahm_jags_categ3,
+                             x         = stack_prediction,
+                             type      = "richness")
+  
+  plot(predictions_rich, col = hcl.colors(100), asp = 1)

+

+

+
Code
+
  # mean = mean species richness estimate
+  # sd = standard deviation of richness estimate

+

 
In these examples we did not return confidence intervals. Doing so
 would be easy by setting intervals =  "confidence".

 
Computation can take long on larger rasters and can require lots of
@@ -5243,74 +4356,94 @@ 
Percentage of area occupied (by species)

 
For demonstration, let’s create a AOI as a simple feature (from the
 sf package) Normally you’d just load a shapefile or other spatial vector
 data set.

-
Sr1 = st_polygon(list(cbind(c(10, 10, 30, 30, 10),
-                    c(40, 70, 80, 40, 40))))
-aoi <- data.frame(name = "My AOI")
-st_geometry(aoi) <- st_geometry(Sr1)

-
plot(r_volcano)
-plot(vect(aoi), add = T)   # easier plotting via SpatVect object from terra 

-

+
Code
+
Sr1 = st_polygon(list(cbind(c(10, 10, 30, 30, 10),
+                    c(40, 70, 80, 40, 40))))
+aoi <- data.frame(name = "My AOI")
+st_geometry(aoi) <- st_geometry(Sr1)

+

+
Code
+
plot(r_volcano)
+plot(vect(aoi), add = T)   # easier plotting via SpatVect object from terra 

+

+

 
To use the vector layer as an AOI, we need to rasterize it. PAO will
 refer to all cells in the AOI raster that are not NA.

-
r_aoi <- rasterize(aoi, r_volcano)
-plot(r_aoi)

-

+
Code
+
r_aoi <- rasterize(aoi, r_volcano)
+plot(r_aoi)

+

+

 
Now the prediction of PAO:

-
  predictions_pao <- predict(object    = model_jags_categ3, 
-                             mcmc.list = out_ahm_jags_categ3,
-                             x         = stack_prediction,
-                             type      = "pao", 
-                             aoi       = r_aoi)

+
Code
+
  predictions_pao <- predict(object    = model_jags_categ3, 
+                             mcmc.list = out_ahm_jags_categ3,
+                             x         = stack_prediction,
+                             type      = "pao", 
+                             aoi       = r_aoi)

+

 
## No id variables; using all as measure variables
 ## No id variables; using all as measure variables

 
Output is a list with several components:

 
$pao_summary is a summary table with quantiles of predicted PAO for
 all species.

-
predictions_pao$pao_summary

+
Code
+
predictions_pao$pao_summary

+

 
##       Min. lower.ci.0.025 1st Qu. Median  Mean 3rd Qu. upper.ci.0.975  Max.
-## sp1  0.017          0.038   0.183  0.371 0.393   0.607          0.819 0.941
-## sp2  0.014          0.047   0.181  0.359 0.372   0.538          0.765 0.846
-## sp3  0.027          0.062   0.178  0.297 0.350   0.510          0.821 0.883
-## sp4  0.226          0.321   0.655  0.771 0.735   0.862          0.936 0.966
-## sp5  0.040          0.065   0.176  0.301 0.376   0.561          0.887 0.907
-## sp6  0.240          0.385   0.760  0.879 0.821   0.923          0.969 0.979
-## sp7  0.046          0.105   0.326  0.551 0.522   0.721          0.881 0.903
-## sp8  0.014          0.023   0.089  0.175 0.207   0.286          0.566 0.660
-## sp9  0.019          0.103   0.290  0.506 0.479   0.657          0.849 0.913
-## sp10 0.019          0.033   0.088  0.161 0.228   0.354          0.605 0.744
-## sp11 0.146          0.189   0.569  0.744 0.673   0.825          0.889 0.919
-## sp12 0.071          0.151   0.476  0.644 0.597   0.754          0.892 0.913
-## sp13 0.010          0.020   0.087  0.184 0.229   0.319          0.701 0.850
-## sp14 0.116          0.190   0.471  0.660 0.623   0.786          0.915 0.940
-## sp15 0.030          0.060   0.227  0.450 0.432   0.619          0.835 0.889
-## sp16 0.070          0.128   0.306  0.488 0.484   0.671          0.800 0.893
-## sp18 0.010          0.040   0.188  0.386 0.384   0.565          0.801 0.867
-## sp19 0.026          0.058   0.243  0.461 0.434   0.597          0.829 0.859
-## sp20 0.014          0.088   0.285  0.479 0.462   0.617          0.835 0.923

+## sp1  0.031          0.050   0.130  0.229 0.296   0.430          0.745 0.820
+## sp2  0.007          0.028   0.104  0.267 0.297   0.440          0.734 0.814
+## sp3  0.007          0.030   0.120  0.232 0.309   0.490          0.775 0.886
+## sp4  0.084          0.151   0.617  0.761 0.697   0.865          0.928 0.956
+## sp5  0.041          0.056   0.132  0.223 0.284   0.400          0.765 0.833
+## sp6  0.247          0.332   0.691  0.850 0.779   0.911          0.968 0.990
+## sp7  0.017          0.056   0.228  0.408 0.440   0.659          0.871 0.933
+## sp8  0.016          0.033   0.072  0.114 0.166   0.219          0.484 0.657
+## sp9  0.030          0.062   0.243  0.383 0.407   0.580          0.820 0.857
+## sp10 0.014          0.024   0.076  0.136 0.175   0.198          0.605 0.736
+## sp11 0.040          0.141   0.555  0.739 0.663   0.814          0.911 0.936
+## sp12 0.051          0.108   0.387  0.603 0.569   0.767          0.928 0.937
+## sp13 0.003          0.015   0.068  0.134 0.200   0.299          0.665 0.701
+## sp14 0.074          0.141   0.399  0.613 0.586   0.789          0.912 0.961
+## sp15 0.024          0.047   0.179  0.335 0.371   0.557          0.835 0.870
+## sp16 0.049          0.111   0.289  0.526 0.493   0.702          0.820 0.846
+## sp18 0.027          0.039   0.172  0.387 0.402   0.619          0.821 0.899
+## sp19 0.003          0.032   0.139  0.307 0.363   0.560          0.839 0.869
+## sp20 0.006          0.028   0.216  0.446 0.447   0.667          0.873 0.909

$pao_matrix contains the PAO estimates for all species from all posterior draws (not shown here)

-
predictions_pao$pao_matrix
+
Code +
predictions_pao$pao_matrix
+

$pao_df is similar to $pao_matrix, but in a data frame in long format, ready for use in ggplot2.

-
  ggplot2::ggplot(predictions_pao$pao_df, aes(x = Species , y = PAO)) +
-  geom_boxplot(outlier.shape = NA) +
-  theme(axis.text.x = element_text(angle = 45, hjust = 1))
-

+
Code +
  ggplot2::ggplot(predictions_pao$pao_df, aes(x = Species , y = PAO)) +
+  geom_boxplot(outlier.shape = NA) +
+  theme(axis.text.x = element_text(angle = 45, hjust = 1))
+
+

$aoi is a raster layer (SpatRaster) showing the AOI (now with cells removed that are NA in the prediction rasters)

-
plot(predictions_pao$aoi)
-

+
Code +
plot(predictions_pao$aoi)
+
+

bayesplot

We can also use the powerful bayesplot package for visualizing model estimates. One downside is that it won’t automatically insert the species names.

-
library(bayesplot)
+
Code +
library(bayesplot)
+

Here is an example for density plots showing the species-level effects of habitat on occupancy:

-
mcmc_areas(out_ahm_nimble2, regex_pars = "beta.ranef.cont.habitat")
-

+
Code +
mcmc_areas(out_ahm_nimble2, regex_pars = "beta.ranef.cont.habitat")
+
+

bayesplot has many other visualization options beyond this example. See the respective documentation, e.g. https://mc-stan.org/bayesplot/.

@@ -5319,12 +4452,16 @@

Model diagnostics

Output of fit() is a coda::mcmc.list, so plot() provides traceplots (not shown here for space reasons):

-
plot(out_ahm_nimble2)
+
Code +
plot(out_ahm_nimble2)
+

The Gelman-Rubin convergence statistic (potential scale reduction factor) can be calculated with gelman.diag() from coda:

-
gd <- coda::gelman.diag(out_ahm_nimble2,  multivariate = FALSE)
-gd
+
Code +
gd <- coda::gelman.diag(out_ahm_nimble2,  multivariate = FALSE)
+gd
+
## Potential scale reduction factors:
 ## 
 ##                                 Point est. Upper C.I.
@@ -5415,8 +4552,10 @@ 
Model diagnostics

 
One can also use gelman.plot to create
 Gelman-Rubin-Brooks plots (evolution of Gelman-Rubin statistic over
 iterations):

-
# not shown here due to space reasons
-coda::gelman.plot(out_ahm_nimble2,  multivariate = FALSE)

+
Code
+
# not shown here due to space reasons
+coda::gelman.plot(out_ahm_nimble2,  multivariate = FALSE)

+

 
 

 
Royle-Nichols model

@@ -5427,104 +4566,130 @@ 
Royle-Nichols model

 is a derived parameter: the probability that at least one individual is
 available for detection at the site. It is achieve by setting the
 parameter “model” to “RN”:

-
modelFile_jags_RN <- tempfile(fileext = ".txt")
-
-model1_jags_RN <- communityModel(
-  model = "RN",
-  occuCovs = list(ranef = c("habitat")),
-  detCovsObservation = list(ranef = "wind"),
-  intercepts = list(det = "ranef", occu = "ranef"),
-  data_list = input_AHM,
-  modelFile = modelFile_jags_RN)

+
Code
+
modelFile_jags_RN <- tempfile(fileext = ".txt")
+
+model1_jags_RN <- communityModel(
+  model = "RN",
+  occuCovs = list(ranef = c("habitat")),
+  detCovsObservation = list(ranef = "wind"),
+  intercepts = list(det = "ranef", occu = "ranef"),
+  data_list = input_AHM,
+  modelFile = modelFile_jags_RN)

+

 
## Site covariates may not be scaled. Scaling is strongly recommended for numeric covariates.

-
## Wrote model to C:\Users\NIEDBA~1\AppData\Local\Temp\Rtmpqkffp1\file1fb846ea4aff.txt

-
out_ahm_jags_RN <- fit(model1_jags_RN, 
-                    n.iter = 5000, 
-                    n.burnin = 2500,
-                    thin = 5,
-                    chains = 3,
-                    quiet = T
-)

+
## Wrote model to C:\Users\NIEDBA~1\AppData\Local\Temp\Rtmp4UNIvN\file4c5845247c3f.txt

+
Code
+
out_ahm_jags_RN <- fit(model1_jags_RN, 
+                    n.iter = 5000, 
+                    n.burnin = 2500,
+                    thin = 5,
+                    chains = 3,
+                    quiet = T
+)

+

 
The models can be processed as shown before, e.g.:

 
Summary statistics of model parameters:

-
DT::datatable(round(summary(out_ahm_jags_RN)$statistics, 3))

-

-
+
Code
+
DT::datatable(round(summary(out_ahm_jags_RN)$statistics, 3))

+

+

+
 
Response curves (marginal effect plots).

-
# In the occupancy (state) submodel, one can show abundance (the latent quantity modelled) 
-plot_effects(model1_jags_RN, 
-             out_ahm_jags_RN,
-             submodel = "state",
-             response = "abundance")

+
Code
+
# In the occupancy (state) submodel, one can show abundance (the latent quantity modelled) 
+plot_effects(model1_jags_RN, 
+             out_ahm_jags_RN,
+             submodel = "state",
+             response = "abundance")

+

 
## $habitat

-

-
# In the occupancy (state) submodel, one can also show occupancy (which is derived from abundance):
-plot_effects(model1_jags_RN, 
-             out_ahm_jags_RN,
-             submodel = "state",
-             response = "occupancy")

+

+
Code
+
# In the occupancy (state) submodel, one can also show occupancy (which is derived from abundance):
+plot_effects(model1_jags_RN, 
+             out_ahm_jags_RN,
+             submodel = "state",
+             response = "occupancy")

+

 
## $habitat

-

-
plot_effects(model1_jags_RN, 
-             out_ahm_jags_RN,
-             submodel = "det")

+

+
Code
+
plot_effects(model1_jags_RN, 
+             out_ahm_jags_RN,
+             submodel = "det")

+

 
## $wind

-

+

 
Coefficient plots:

-
plot_coef(model1_jags_RN,
-          out_ahm_jags_RN,
-          ordered = FALSE)

+
Code
+
plot_coef(model1_jags_RN,
+          out_ahm_jags_RN,
+          ordered = FALSE)

+

 
## $habitat

-

-
plot_coef(model1_jags_RN,
-          out_ahm_jags_RN,
-          ordered = FALSE,
-          submodel = "det")

+

+
Code
+
plot_coef(model1_jags_RN,
+          out_ahm_jags_RN,
+          ordered = FALSE,
+          submodel = "det")

+

 
## $wind

-

+

 
Predictions work as before, but a Royle-Nichols model allows for
 returning expected species abundance (lambda) via type = “abundance” and
 type = “lambda_array”. Occupancy predictions internally convert
 abundance to occupancy.

 
Occupancy predictions (type = “psi”)

-
predictions_psi_RN <- predict(object    = model1_jags_RN, 
-                             mcmc.list = out_ahm_jags_RN,
-                             x         = stack_prediction,
-                             type      = "psi",
-                             draws     = 200)
-plot(predictions_psi_RN$mean, 
-     col = hcl.colors(100),
-       maxnl = 9)  # first 9 species only

-

+
Code
+
predictions_psi_RN <- predict(object    = model1_jags_RN, 
+                             mcmc.list = out_ahm_jags_RN,
+                             x         = stack_prediction,
+                             type      = "psi",
+                             draws     = 200)
+plot(predictions_psi_RN$mean, 
+     col = hcl.colors(100),
+       maxnl = 9)  # first 9 species only

+

+

 
Abundance predictions (type = “abundance”):

-
predictions_abundance_RN <- predict(object    = model1_jags_RN, 
-                             mcmc.list = out_ahm_jags_RN,
-                             x         = stack_prediction,
-                             type      = "abundance",
-                             draws     = 200)
-plot(predictions_abundance_RN$mean, 
-     col = hcl.colors(100),
-     maxnl = 9)   # first 9 species only

-

+
Code
+
predictions_abundance_RN <- predict(object    = model1_jags_RN, 
+                             mcmc.list = out_ahm_jags_RN,
+                             x         = stack_prediction,
+                             type      = "abundance",
+                             draws     = 200)
+plot(predictions_abundance_RN$mean, 
+     col = hcl.colors(100),
+     maxnl = 9)   # first 9 species only

+

+

 
Richness predictions (type = “richness”)

-
  predictions_richness_RN <- predict(object    = model1_jags_RN, 
-                             mcmc.list = out_ahm_jags_RN,
-                             x         = stack_prediction,
-                             type      = "richness",
-                             draws     = 200)
-plot(predictions_richness_RN$mean, col = hcl.colors(100))

-

+
Code
+
  predictions_richness_RN <- predict(object    = model1_jags_RN, 
+                             mcmc.list = out_ahm_jags_RN,
+                             x         = stack_prediction,
+                             type      = "richness",
+                             draws     = 200)
+plot(predictions_richness_RN$mean, col = hcl.colors(100))

+

+

 
PAO predictions (type = “pao”):

-
  predictions_pao_RN <- predict(object    = model1_jags_RN, 
-                             mcmc.list = out_ahm_jags_RN,
-                             x         = stack_prediction,
-                             type      = "pao",
-                             draws     = 200)

+
Code
+
  predictions_pao_RN <- predict(object    = model1_jags_RN, 
+                             mcmc.list = out_ahm_jags_RN,
+                             x         = stack_prediction,
+                             type      = "pao",
+                             draws     = 200)

+

 
## No id variables; using all as measure variables
 ## No id variables; using all as measure variables

-
DT::datatable(predictions_pao_RN$pao_summary)

-

-
+
Code
+
DT::datatable(predictions_pao_RN$pao_summary)

+

+

+
 

 

 
Single-species occupancy models

@@ -5541,24 +4706,32 @@ 
Single-species occupancy models

 covariates as in the example below.

 
In this example we use species 1 from the community created with
 AHMbook::simComm.

-
input_AHM_sp1 <- list(ylist = ylist[1],    # ylist[1] is still a list, with 1 element
-                      siteCovs = input_AHM$siteCovs,
-                      obsCovs = input_AHM$obsCovs)

+
Code
+
input_AHM_sp1 <- list(ylist = ylist[1],    # ylist[1] is still a list, with 1 element
+                      siteCovs = input_AHM$siteCovs,
+                      obsCovs = input_AHM$obsCovs)

+

 
Create model:

-
modelFile_jags_sp1 <- tempfile(fileext = ".txt")
-
-model_sp1_jags_ranef <- communityModel(
-  occuCovs = list(ranef = c("habitat|categ")),    # response to habitat differs between categ via random effect
-  detCovsObservation = list(fixed = "wind"),
-  intercepts = list(det = "fixed", occu = "fixed"),
-  data_list = input_AHM_sp1,
-  modelFile = modelFile_jags_sp1)

+
Code
+
modelFile_jags_sp1 <- tempfile(fileext = ".txt")
+
+model_sp1_jags_ranef <- communityModel(
+  occuCovs = list(ranef = c("habitat|categ")),    # response to habitat differs between categ via random effect
+  detCovsObservation = list(fixed = "wind"),
+  intercepts = list(det = "fixed", occu = "fixed"),
+  data_list = input_AHM_sp1,
+  modelFile = modelFile_jags_sp1)

+

 
Fit model:

-
fit_jags_sp1 <- fit(model_sp1_jags_ranef, n.iter = 1000, n.burnin = 500, thin = 5)

+
Code
+
fit_jags_sp1 <- fit(model_sp1_jags_ranef, n.iter = 1000, n.burnin = 500, thin = 5)

+

 
Print summary statistics:

-
DT::datatable(round(summary(fit_jags_sp1)$statistics, 3))

-

-
+
Code
+
DT::datatable(round(summary(fit_jags_sp1)$statistics, 3))

+

+

+
 
The plot_effects() and plot_coef()
 functions currently don’t work for this case, since single-species
 models are mostly a side-effect of the intended community occupancy
@@ -5570,69 +4743,8 @@ 
Single-species occupancy models

 
 
 
-

-
-
-
-
-
-
-
-
-
-
 
 
-