5. Multi-species occupancy models
Juergen Niedballa (camtrapr@gmail.com)
-2024-02-24
+2024-02-25
--
+
- 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 +
+
library(camtrapR)
-library(purrr)
-library(DT)
-library(knitr)
-library(ggplot2)
-library(sf)Code
+ +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
+ +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 speciesCode
+ 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 speciesGet 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)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 covariateCode
+ +Now we bundle the necessary data for communityModel in a
list with 3 items:
-
@@ -4212,9 +3150,11 @@
- obsCovs: (named) list with observation level covariates. Requires at least the effort matrix.
Data preparation
data_list <- list(ylist = ylist,
- siteCovs = sitecovs,
- obsCovs = list(effort = DetHist_list[[1]]$effort)) # is identical for all speciesCode
+ +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.txtCode
+ +## Wrote model to C:\Users\NIEDBA~1\AppData\Local\Temp\Rtmp4UNIvN\file4c5879e733bc.txtThis 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)## 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:
## effortThe model is fit with fit().
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)Statistics of parameter estimates:
-# Note, colors may not render correctly in dark themes in RStudio.
-DT::datatable(round(fit_summary$statistics, 3))Code
+ +Quantiles of parameter estimates:
-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")## $utm_y
##
## $elevation
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)## 'combine' and 'ordered' can't both be TRUE. Setting 'ordered = FALSE'
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.
Introduction
faster in some situations.Using Nimble for community occupancy models requires the packages “nimble” and “nimbleEcology”:
-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") 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.txtCode
+ +## Wrote model to C:\Users\NIEDBA~1\AppData\Local\Temp\Rtmp4UNIvN\file4c5826ec4cd4.txtIt 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 harmlessCode
+ +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
+ +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)Response curves (= marginal effect plots):
- 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")
## $effort
Effect size plots:
- 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")
## $effort
Traceplots (not shown)
- plot(fit.nimble.comp)The estimates and marginal effects are not very pretty, and the model likely requires more iterations to achieve smooth estimates.
More complete example
Prepare data
First, install the AHMbook package if necessary:
-install.packages("AHMbook")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
+ +## 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.txtout_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.txtCode
+ +Summary statistics of model parameters:
-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## $habitat
plot_resp_jags_det <- plot_effects(model1_jags,
- out_ahm_jags,
- submodel = "det")
-plot_resp_jags_det
Code
+ +## $wind
Plot coefficient estimates:
-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
+ +## $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_occu2Code
+ +## $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
+ +## 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.txtWe 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
+ +## Defining model## [Note] Using 'y' (given within 'constants') as data.## Building modelNimble Model
## running chain 1...## running chain 2...## running chain 3...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_occuCode
+ +## $habitat
plot_resp_nimble_det <- plot_effects(model1_nimble,
- out_ahm_nimble2,
- submodel = "det")
-plot_resp_nimble_det
Code
+ +## $wind
Plot coefficient estimates:
-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
+ +## $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
+ +## 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.txtout_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.txtCode
+ +Summary statistics of model parameters:
-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.txtWhen 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
+ +## Defining model## [Note] Using 'y' (given within 'constants') as data.## Building modelMissing 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))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
+ +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))Higher numbers of species in data augmentation increase model run time.
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)## $categ
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))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
+ +## 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))Plot the results:
-plot_effects(object = model_jags_categ2,
- mcmc.list = out_ahm_jags_categ2,
- submodel = "det")Code
+ +## $wind_categ
plot_coef(object = model_jags_categ2,
- mcmc.list = out_ahm_jags_categ2,
- submodel = "det")
Code
+ +## $wind_categ
