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:
The main function is communityModel()
. It automates:
The output is a commOccu
object. The workflow provides
methods for:
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:
predict
may not be correct for
theseOver 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:
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:
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.
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.