- Overview
- Warning
- Limitations
- Parameter naming convention
- Quick example
- More complete example
- Further
possibilies
- Fixed effects
- Priors
- Data augmentation
- Categorical covariates
- Random effects other than species
- Nested random effects
- Species and site random effect on detection
- Quadratic covariate effects
- WAIC
- Species richness
- Spatial predictions
- Percentage of area occupied (by species)
- bayesplot
- Model diagnostics
- Royle-Nichols model
- Single-species occupancy models

Multi-species occupancy models are a powerful tool for combining information from multiple species to estimate both individual and community-level responses to environmental variables. These models are typically fit using Bayesian frameworks such as BUGS, JAGS, or Nimble. However, preparing the necessary input from camera trapping data, writing the model code, setting parameters to monitor, and providing initial values can be a daunting and time-consuming task, even for experienced users.

To simplify this process, a new extension of camtrapR automatically generates custom community occupancy models for analysis in JAGS or Nimble. This flexible workflow supports:

- standard and Royle-Nichols occupancy models
- continuous and categorical site covariates on detection and occupancy probability
- continuous and categorical observation-level covariates on detection probability
- fixed and random effects of all these covariates
- independent effects of continuous site covariates between the different species
- nested random effects
- fixed, independent or random intercepts of detection and occupancy probability
- data augmentation (fully open community or known maximum species richness)
- species richness estimates
- models in JAGS and Nimble
- derived quantities, e.g. Bayesian p-values (overall and by species) and species richness estimates

The main function is `communityModel()`

. It automates:

- writing model code
- setting parameters to monitor
- setting initial values
- bundling input data

The output is a `commOccu`

object. The workflow provides
methods for:

- fitting models
- plotting marginal effects (response curves)
- plotting effect sizes
- spatial predictions of species occupancy probabilities
- spatial predictions of species richness

The output of fit() is an mcmc.list from the coda package and can be uses as such, so the summary(), plot(), etc. methods are available as well as the functions from bayesplot.

Unlike packages such as unmarked or ubms, the
`communityModel`

function does not currently support R
formula syntax. Instead, covariates and their effect types are specified
through dedicated function arguments. While this non-standard approach
may change in the future to support formula syntax, it may be less
confusing than the lme4 syntax for random effects.

A word of caution, please always check the model code and make sure you understand the code and agree with it. While the workflow is designed for convenience and does not require you to examine the model code, it is highly recommended that you do so. The model code contains many comments to facilitate understanding. This workflow still requires a proper understanding of Bayesian modelling frameworks and the concepts behind multi-species occupancy models.

Please note that the workflow is still experimental and it is important to carefully review and critically evaluate the results. If you encounter any issues, do not hesitate to report them.

Currently, there are some limitations to the workflow:

- priors can not be customized
- models with categorical observation-level covariates cannot be compiled in Nimble and must be fit in JAGS instead
- nested random effects and random effects of site covariates are still experimental, and the plot_* methods don’t work on these yet
- likewise, the output of
`predict`

may not be correct for these - input checks are a bit limited and may not always lead to meaningful error messages. This will improve over time with user feedback.
- Nimble does not yet seem to have an advantage in processing speed. That might require further customization
- independent effects are only supported in continuous site covariates

Over time these will likely be worked on.

The parameter names are assembled from standardized building blocks. The nomenclature is as follows:

Name | Refers to | Description |
---|---|---|

alpha | Submodel | detection submodel |

beta | Submodel | occupancy submodel |

0 | Intercept | denotes the intercepts (alpha0, beta0) |

fixed | Effect type | fixed effects (constant across species) |

indep | Effect type | independent effects (separate and independent between species) |

ranef | Effect type | random effects (of species and/or other categorical covariates) |

cont | Covariate type | continuous covariates |

categ | Covariate type | categorical covariates |

mean | Hyperparameter | mean of random effect |

sigma | Hyperparameter | standard deviation of random effect |

tau | Hyperparameter | precision of random effect (used internally, not returned) |

For example, a fixed intercept of occupancy (constant across species)
is *beta0*, and a fixed intercept of detection probability is
*alpha0*.

An occupancy probability intercept with a random effect of species is:

**beta0.mean** community mean of the occupancy
probability intercept

**beta0.sigma** standard deviation of the community
occupancy probability intercept.

**beta0[1]** occupancy probability intercept of species
1 (likewise for other species).

For effects of site covariates, the pattern is:

*submodel.effectType.covariateType.CovariateName.hyperparameter*

For example:

**beta.ranef.cont.habitat.mean** is the mean community
effect of the continuous site covariate ‘habitat’ on occupancy
probability.

**beta.ranef.cont.habitat[1]** is the effect of
continuous site covariate ‘habitat’ on occupancy probability of species
1.

Site-occasion covariates are denoted by *.obs* after the
submodel, e.g.:

**alpha.obs.fixed.cont.effort** is the fixed effect of
the continuous observation-level covariate ‘effort’ on detection
probability

Below is a brief example using the sample data included in camtrapR. This data set contains only five species at three fictional camera trap stations and uses made-up covariates. As a result, the parameter estimates will not have any ecological significance. This example serves only as a technical demonstration of how community models fit into the camtrapR workflow.

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

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
```

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

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

Now we bundle the necessary data for `communityModel`

in a
list with 3 items:

- ylist: the list of detection histories
- siteCovs: data frame with site covariates
- obsCovs: (named) list with observation level covariates. Requires at least the effort matrix.

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. This ensures that the detection probability is set to 0 when there was no effort. This is achieved by creating a binary effort matrix based on the user-provided effort matrix. In the binary effort matrix, all cells with effort = NA are set to 0, while all cells with values are set to 1, regardless of the value of effort (as long as it is not NA). As a result, all occasions with effort = NA in the provided effort matrix will have a detection probability of 0.

This does **not** automatically include effort as a
covariate on detection probability. To include effort as a covariate on
detection probability, include “effort” in argument
`detCovsObservation`

(you can specify a fixed or random
effect of effort on detection probability).

Now the interesting part. The function communityModel() creates a new class, “commOccu”. It is an object that contains all relevant information for running a community occupancy model. Here we will create a JAGS model (the default). For a model in Nimble, see further below.

`## 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 a covariate on detection probability (constant across species), and there are random effects of species on the intercepts of detection and occupancy probabilities.

Note that throughout this vignette we use a temporary file to save the model text. In real examples you’d probably want to use a persistent file.

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

There is a very basic summary method for commOccu objects:

```
## commOccu object for community occupancy model (in JAGS)
##
## 5 species, 3 stations, 7 occasions
## 21 occasions with effort
## Number of detections (by species): 2 - 11
##
## Available site covariates:
## Station, utm_y, utm_x, elevation, fact
##
## Used site covariates:
## utm_y, elevation
##
## Available site-occasion covariates:
## effort
```

The model is fit with `fit()`

.

The output of the `fit()`

method is an
`mcmc.list`

from the coda package.

Summarize posterior estimates:

Statistics of parameter estimates:

Quantiles of parameter estimates:

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:

`## $utm_y`

```
##
## $elevation
```

`## $effort`

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

:

`## 'combine' and 'ordered' can't both be TRUE. Setting 'ordered = FALSE'`

`## $effort`

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

NIMBLE is a system for building and sharing analysis methods for statistical models, especially for hierarchical models and computationally-intensive methods. NIMBLE is built in R but compiles your models and algorithms using C++ for speed.

It allows you to run the same models as BUGS or JAGS, but can be much faster in some situations.

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

Using nimble only makes sense if one can compile the models (which runs them as fast C++ code). Compiling Nimble code requires Rtools.

You can obtain Rtools from: https://cran.r-project.org/bin/windows/Rtools/.

It then needs to be put in PATH, see here for how to do it: https://cran.r-project.org/bin/windows/Rtools/#putting-rtools-on-the-path.

Depending on the version and where you installed Rtools, you should 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).

We can fit the same model in Nimble using the exact same functions as
for the JAGS model above. Only 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.

To fit a compiled model, set compile = TRUE:

The output as an `mcmc.list`

, just like for the JAGS model
above, and can be worked with just the same:

Response curves (= marginal effect plots):

`## $utm_x`

```
##
## $utm_y
```

`## $effort`

Effect size plots:

`## 'combine' and 'ordered' can't both be TRUE. Setting 'ordered = FALSE'`

`## $effort`

Traceplots (not shown)

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

The example above was very small and unrealistic. We can’t provide a
large camera trapping data set in camtrapR, but we can simulate one
using the `simComm`

function from the `AHMbook`

package. Using a simulated data set also means we can compare model
estimates to the “true” parameters (which were used to create the
simulated data).

First, install the AHMbook package if necessary:

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)
```

```
## 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))
```

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

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

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

Summary statistics of model parameters:

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
`beta.ranef.cont.habitat.sigma`

in the table above for
estimates.

The mean effect of wind on detection probability was -0.5 (SD = 1).
See `alpha.obs.ranef.cont.wind.mean`

and
`alpha.obs.ranef.cont.wind.sigma`

in the table above for
estimates.

Bayesian p-values can be inspected to check for lack of fit (overall or by species). Lack of fit at significance level = 0.05 is indicated by Bayesian p-values below 0.025 or greater than 0.975. In the example above, the overall Bayesian p-value (Bpvalue) indicates no problems with lack of fit. Likewise, species-level Bayesian p-values (Bpvalue_species) indicate no lack of fit for any species. This is not surprising because the model uses simulated data.

Plot response curves (marginal effects)

`## $habitat`

`## $wind`

Plot coefficient estimates:

`## $habitat`

`## $wind`

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

:

`## $habitat`

Here’s the same model in Nimble:

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

`## 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).

`## Defining model`

`## [Note] Using 'y' (given within 'constants') as data.`

`## Building model`

`## Setting data and initial values`

```
## Running calculate on model
## [Note] Any error reports that follow may simply reflect missing values in model variables.
```

`## Checking model sizes and dimensions`

```
## [Note] This model is not fully initialized. This is not an error.
## To see which variables are not initialized, use model$initializeInfo().
## For more information on model initialization, see help(modelInitialization).
```

```
## Compiling
## [Note] This may take a minute.
## [Note] Use 'showCompilerOutput = TRUE' to see C++ compilation details.
## Compiling
## [Note] This may take a minute.
## [Note] Use 'showCompilerOutput = TRUE' to see C++ compilation details.
```

`## running chain 1...`

`## running chain 2...`

`## running chain 3...`

Plot responses (marginal effects):

`## $habitat`

`## $wind`

Plot coefficient estimates:

`## $habitat`

`## $wind`

The examples so far did not have any missing data in the detection histories. Real-world data commonly have missing data though, e.g. due to camera malfunction or specific study designs. In single-species occupancy models (e.g. in unmarked or ubms) this is handled automatically, but requires a little more attention in these 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))
```

And fit models in JAGS and Nimble:

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

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

Summary statistics of model parameters:

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.

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

`## 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).

`## Defining model`

`## [Note] Using 'y' (given within 'constants') as data.`

`## Building model`

`## Setting data and initial values`

```
## Running calculate on model
## [Note] Any error reports that follow may simply reflect missing values in model variables.
```

`## Checking model sizes and dimensions`

```
## [Note] This model is not fully initialized. This is not an error.
## To see which variables are not initialized, use model$initializeInfo().
## For more information on model initialization, see help(modelInitialization).
```

```
## Compiling
## [Note] This may take a minute.
## [Note] Use 'showCompilerOutput = TRUE' to see C++ compilation details.
## Compiling
## [Note] This may take a minute.
## [Note] Use 'showCompilerOutput = TRUE' to see C++ compilation details.
```

`## running chain 1...`

`## running chain 2...`

`## running chain 3...`

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

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

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
detection probabilities, information about occasions without effort is
provided via the mandatory argument `effortCov`

in the
`communityModel`

function. Use the `$effort`

part
of the output of `detectionHistory()`

to provide the
necessary information, and of course ensure that
`detectionHistory`

is run with
`includeEffort = T`

. All occasions with NA in the effort
matrix will have their detection probabilities fixed at 0. This ensures
that those occasions do not affect the model fit. Be aware that this
does not include effort as a covariate on detection. Including effort as
a covariate on detection would need to be done via
`detCovsObservation`

.

This section demonstrates additional possibilities. The models here are (mostly) fit in JAGS, but can also be fit with Nimble.

One can define fixed effects of covariates if the responses are thought to be constant across all species. This might make sense for the effect of sampling effort in camera trapping studies.

It is not yet possible to customize priors. The current defaults are:

- dnorm(0, 0.05) : random effect means
- dgamma(0.1, 0.1): random effect precision (tau)

Future versions may allow more flexibility. Please contact me if this is relevant for you. In the meantime feel free to adjust the model code in the txt file as needed.

`communityModel`

has an argument
`augmentation`

, which can implement fully open data
augmentation or data augmentation up to a known maximum number of
species.

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
)
```

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
)
```

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

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

Site covariates will be interpreted as categorical if they are factors in the input data frame. If a random effect of a categorical covariate is defined, there will be independent species-level random effects within each factor level of the categorical covariate. In other words, each factor level will have its own and independent random effect 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
)
```

Plot the results:

`## $categ`

`## $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
for `mean`

and `SD`

, respectively. All estimates
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).

Currently categorical observation-level covariates are only fully 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.

`## 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
)
```

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`

.

Plot the results:

`## $wind_categ`

`## $wind_categ`

By default random effects are by species. Alternatively, you can specify random effects of a categorical site covariate2 on continuous covariate1 using the pattern:

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

The same syntax can be used in

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

and

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
)
```

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).

Currently random effects other than species can not be plotted with
`plot_effects`

and `plot_coef`

.

Consider a numeric covariate1 and categorical covariate2.

Use the keyword “+Species” to specify a species-specific random effect within each factor level of covariate2:

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
)
```

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.

Currently nested random effects can not be plotted with
`plot_effects`

and `plot_coef`

.

Species and site random effect on detection probability can be
helpful if Bayesian p-values indicate lack of fit (overall and / or for
individual species). Such an effect can be added via argument
`speciesSiteRandomEffect`

. Setting it to TRUE adds a random
effect “ran” to detection probability in the form:

logit.p[i,j,k] <- alpha0[i] + … + **ran[i, j]**

where i,j are species and station and occasion indices.
`ran`

is defined as:

ran[i,j] ~ dnorm(0, tau.ran)

with `tau.ran`

being the precision of the normal
distribution (precision = 1 / standard deviation).

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
)
```

The model returns the standard deviation of this random effect as
`alpha.speciessite.ranef.sigma`

, which can be found in the
table below.

The easiest way to include quadratic effects is by squaring the covariate and including it as a new covariate. The plot_* and predict methods will be able to identify the quadratic term via a required keyword (the suffix “_squared” by default).

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

Create quadratic covariate:

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
)
```

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”).

`## $habitat`

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

Nimble models can return WAIC. Set argument
`fit(..., 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:

WAIC:

```
## nimbleList object of type waicNimbleList
## Field "WAIC":
## [1] 1071.5
## Field "lppd":
## [1] -490.7924
## Field "pWAIC":
## [1] 44.95773
```

Global species richness estimates is returned in parameter “Nspecies” (see below).

Species richness by station is returned in parameter “Nspecies_station” (see below).

This is only possible in JAGS models. Argument “richnessCategories” 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
)
```

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”, “low”, and “high” respectively.
- “Nspecies_by_station[1]” to “…[30]” are species richness estimates for each sampling station.

You can create spatial predictions based on the occupancy models
created with `communityModel`

, their model fits and covariate
rasters.

This requires the terra package.

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
)
```

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

**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).

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

Species occupancy maps:

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

Species richness maps:

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 RAM. Check RAM usage during computations. Parameter “batch” can help reduce memory load.

The `predict`

function can be used to calculate Percentage
of Area Occupied (PAO) for all species, optionally for a given area of
interest (AOI).

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.

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.

Now the prediction of PAO:

```
## 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.

```
## Min. lower.ci.0.025 1st Qu. Median Mean 3rd Qu. upper.ci.0.975 Max.
## 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)

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

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

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

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

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

Output of `fit()`

is a coda::mcmc.list, so
`plot()`

provides traceplots (not shown here for space
reasons):

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

from
`coda`

:

```
## Potential scale reduction factors:
##
## Point est. Upper C.I.
## alpha.obs.ranef.cont.wind[1] 1.00 1.00
## alpha.obs.ranef.cont.wind[2] 1.00 1.00
## alpha.obs.ranef.cont.wind[3] 1.00 1.00
## alpha.obs.ranef.cont.wind[4] 1.00 1.00
## alpha.obs.ranef.cont.wind[5] 1.00 1.00
## alpha.obs.ranef.cont.wind[6] 1.00 1.01
## alpha.obs.ranef.cont.wind[7] 1.00 1.00
## alpha.obs.ranef.cont.wind[8] 1.00 1.00
## alpha.obs.ranef.cont.wind[9] 1.00 1.00
## alpha.obs.ranef.cont.wind[10] 1.00 1.00
## alpha.obs.ranef.cont.wind[11] 1.00 1.00
## alpha.obs.ranef.cont.wind[12] 1.00 1.00
## alpha.obs.ranef.cont.wind[13] 1.00 1.00
## alpha.obs.ranef.cont.wind[14] 1.01 1.02
## alpha.obs.ranef.cont.wind[15] 1.00 1.00
## alpha.obs.ranef.cont.wind[16] 1.00 1.00
## alpha.obs.ranef.cont.wind[17] 1.00 1.00
## alpha.obs.ranef.cont.wind[18] 1.00 1.00
## alpha.obs.ranef.cont.wind[19] 1.00 1.01
## alpha.obs.ranef.cont.wind.mean 1.00 1.00
## alpha.obs.ranef.cont.wind.sigma 1.00 1.01
## alpha0[1] 1.00 1.01
## alpha0[2] 1.00 1.01
## alpha0[3] 1.00 1.00
## alpha0[4] 1.00 1.01
## alpha0[5] 1.00 1.00
## alpha0[6] 1.00 1.00
## alpha0[7] 1.00 1.00
## alpha0[8] 1.01 1.02
## alpha0[9] 1.00 1.00
## alpha0[10] 1.00 1.00
## alpha0[11] 1.00 1.01
## alpha0[12] 1.00 1.01
## alpha0[13] 1.00 1.00
## alpha0[14] 1.00 1.00
## alpha0[15] 1.00 1.00
## alpha0[16] 1.00 1.00
## alpha0[17] 1.00 1.01
## alpha0[18] 1.00 1.00
## alpha0[19] 1.00 1.01
## alpha0.mean 1.00 1.01
## alpha0.sigma 1.00 1.01
## beta.ranef.cont.habitat[1] 1.00 1.01
## beta.ranef.cont.habitat[2] 1.00 1.01
## beta.ranef.cont.habitat[3] 1.00 1.00
## beta.ranef.cont.habitat[4] 1.00 1.00
## beta.ranef.cont.habitat[5] 1.01 1.02
## beta.ranef.cont.habitat[6] 1.00 1.01
## beta.ranef.cont.habitat[7] 1.00 1.00
## beta.ranef.cont.habitat[8] 1.00 1.00
## beta.ranef.cont.habitat[9] 1.00 1.01
## beta.ranef.cont.habitat[10] 1.00 1.00
## beta.ranef.cont.habitat[11] 1.00 1.00
## beta.ranef.cont.habitat[12] 1.00 1.01
## beta.ranef.cont.habitat[13] 1.00 1.00
## beta.ranef.cont.habitat[14] 1.00 1.00
## beta.ranef.cont.habitat[15] 1.00 1.01
## beta.ranef.cont.habitat[16] 1.00 1.00
## beta.ranef.cont.habitat[17] 1.00 1.00
## beta.ranef.cont.habitat[18] 1.00 1.00
## beta.ranef.cont.habitat[19] 1.00 1.00
## beta.ranef.cont.habitat.mean 1.00 1.00
## beta.ranef.cont.habitat.sigma 1.00 1.01
## beta0[1] 1.00 1.00
## beta0[2] 1.00 1.00
## beta0[3] 1.01 1.02
## beta0[4] 1.00 1.00
## beta0[5] 1.00 1.00
## beta0[6] 1.00 1.02
## beta0[7] 1.00 1.00
## beta0[8] 1.00 1.02
## beta0[9] 1.00 1.00
## beta0[10] 1.00 1.00
## beta0[11] 1.00 1.01
## beta0[12] 1.00 1.01
## beta0[13] 1.00 1.01
## beta0[14] 1.00 1.00
## beta0[15] 1.00 1.00
## beta0[16] 1.00 1.00
## beta0[17] 1.00 1.01
## beta0[18] 1.00 1.00
## beta0[19] 1.00 1.01
## beta0.mean 1.01 1.02
## beta0.sigma 1.00 1.02
```

One can also use `gelman.plot`

to create
Gelman-Rubin-Brooks plots (evolution of Gelman-Rubin statistic over
iterations):

Instead of a standard occupancy model one can fit the occupancy model of Royle and Nichols (2003), which relates probability of detection of the species to the number of individuals available for detection at each site (instead of their occupancy status alone). Probability of occupancy 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”:

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

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

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

Summary statistics of model parameters:

Response curves (marginal effect plots).

`## $habitat`

`## $habitat`

`## $wind`

Coefficient plots:

`## $habitat`

`## $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”)

Abundance predictions (type = “abundance”):

Richness predictions (type = “richness”)

PAO predictions (type = “pao”):

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

Incidentally, the workflow above can be used for single-species models too. The code will be less streamlined and possibly slower than normal code for single-species models (due to the species loop which is unnecessary for single-species models), but it is possible.

To create a single-species occupancy model, simply pass a single
detection history in the data_list$ylist of the input to
`communityModel`

.

Take care to only define fixed effects for covariates and the intercepts, not random effects of species. It is possible however to specify random effects of categorical site covariates on continuous site covariates as in the example below.

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

.

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)
```

Fit model:

Print summary statistics:

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
workflow and not the main goal. There are other great options for
Bayesian single-species occupancy models, e.g. in the packages ubms or
wiqid, or unmarked for a frequentist implementation.