In dslramsey/eradicate: Functions to support the "eradication tools" shiny app

knitr::opts_chunk$set(
  collapse = TRUE,
  comment = "#>"
)

Eradication tools user guide

Dave Ramsey (22/9/2020)

The eradicate package contains some functions to undertake basic populaton estimation from data that would typically be collected as part of a pest eradication program. The package borrows heavily from ideas and code used in the unmarked package, and additionally includes some custom models not included in that package. It additionally uses the simpler S3 classes for object orientation and method dispatch, which are less formal that the S4 classes used in unmarked.

Here we illustrate the basic use of the models in the eradicate package by running through a theoretical program to eradicate feral cats on San Nicolas Island, off the coast of California. We start by undertaking a preliminary assessment of the population (assessment phase) using camera trap sampling before attempting to eradicate the population using removal trapping with leg-hold traps (removal phase). During the removal phase we will undertake periodic assessments of the population to track our progress towards eradication.

Since this is a simulated dataset, we know exactly what the population size was before the eradication program began (86 cats). Before we start we load various the packages necessary for our analysis ...

library(eradicate)
library(sf)
library(terra)

... and then read in the shapefile of the region and a raster of habitat values. The raster has a resolution (cell size) of 100 m indicating the presence of grassland or forest.

region<- read_sf(system.file("extdata", "shape/san_nic_region.shp", package="eradicate"))
habitat<- rast(system.file("extdata", "san_nic_habitat.tif", package="eradicate"))

plot(habitat, axes=F)
plot(st_geometry(region), add=T)

Asessment phase

We undertake an initial assessment of the population using a grid of 55 camera traps spread systematically across the island at 1 km spacing. The example dataset san_nic_pre, contains the camera locations as well as simulated counts of cat encounters at these cameras over 30 nights sampling.

traps<- san_nic_pre$traps
counts<- san_nic_pre$counts
dim(counts)

plot(habitat, axes = FALSE)
plot(st_geometry(region), add=TRUE)
points(traps, pch=16)

It was thought that the cats on the island had a tendency to favour the grassland over the forest habitat so this information may be useful for modelling the variation in cat abundance across the island. We therefore extract habitat information from habitat for use in further modelling. Cells in the habitat raster took values of either 1 (grassland) or 0 (forest). Rather than just extract the habitat value from the focal cell containing each trap location, we get the average of all cells within 500 m of the focal cell

traps_sf<- st_as_sf(traps, coords=c(1,2), crs=st_crs(region))
traps_buff<- st_buffer(traps_sf, dist=500)
pgrass<- terra::extract(habitat, vect(traps_buff), fun=mean, na.rm=TRUE)
names(pgrass)<- c("id","pgrass")
site.data<- cbind(pgrass, traps)

We are now ready to fit some models! There are two types of model suitable for the basic count data collected from these cameras. Since the data are counts, and assuming that the counts from each camera are independent (i.e. cats recorded on one camera are not recorded on another camera), we can choose to fit an occupancy type model to estimate the proportion of sites occupied by cats or we can estimate the abundance of cats on the sampled sites. For occupancy analysis, eradicate can fit the basic Mackenzie type occupancy model using the function occuM(). To estimate abundance, we have two types of model, either the Nmixture model of Royle (2004) nmix() or the occupancy model of Royle and Nichols (2003) occuRN(). Although the latter is strictly an occupancy model, it also estimates abundance by positing a relationship between occupancy and abundance.

The first step in the analysis is to collate the count data and the associated covariate information (i.e. habitat) into a custom data class using eFrame(). See ?eFrame for more info on organizing the data necessary to use the models in eradicate.

emf<- eFrame(counts, siteCovs = site.data)
summary(emf)

# Occupancy model
m1<- occM(~pgrass, ~1, data=emf)
summary(m1)

# Royle/Nichols model
m2<- occRN(~pgrass, ~1, K=50, data=emf)
summary(m2)

m3<- nmix(~pgrass, ~1, K=50, data=emf)
summary(m3)

These three models have a common modelling interface. The first argument specifies the model for the abundance or occupancy component, while the second argment specifies the model for the detection component. Unlike unmarked, which can include covariates at the observation level, eradicate currently only uses site-level covariates to model the detection process. For occuRN() and nmix() we also need to specify the parameter K, which is an integer specifying the possible upper value for abundance at each site. Note that all three models give similar parameter estimates.

In order to get estimates of either the proportion of sites occupied or total abundance, we can use the function calcN(), which has an additional argument newdata where the user can supply covariate values for prediction sites. Since we wish to predict only for the sampled sites, we can simply use calcN() without any newdata argument.

nhat1<- calcN(m1)
nhat2<- calcN(m2)
nhat3<- calcN(m3)

nhat1$Occ
nhat2$Nhat
nhat3$Nhat

The return values for calcN contain the predicted occupancy or abundance for each site as well as for the total population over all sites. Confidence intervals for abundance are calculated using a log-normal transformation, which are guaranteed not to include negative values in the interval.

In addition to these three models, eradicate additionally includes another model that is exclusively applied to camera data. Introduced by Nakashima et al. (2018), the REST model (short for Random Encounter and Staying Time model) uses data on the number of encounters per camera as well as additional information on the time taken for each detected individual to cross the cameras field of view (the staying time). Under the assumption that detection in the field of view (or some known proportion of it) is perfect, the REST model can estimate animal density. Further information on this model is available in Nakashima et al (2018) and in the help file.

We fit the REST model to the example data found in the san_nic_rest dataset. A more involved example that uses additional distance sampling data from camera traps to estimate the effective sampling area is provided in the vignette REST_example.

counts<- san_nic_rest$y
stay<- san_nic_rest$stay
cens<- san_nic_rest$cens
A<- san_nic_rest$area # in km2
active<- san_nic_rest$active_hours

emf<- eFrameREST(counts, stay, cens, A, active, siteCovs = site.data)

m4<- REST(~pgrass, data=emf)
summary(m4)

calcN(m4)$Nhat

Eradication phase

Now that we have some idea of the size of the cat population on San Nicolas (i.e. around 70-90 individuals, which compares well with the true population size of 86), we now turn to the eradication phase. The eradication will be attempted using 100 leg-hold traps located across the island. At a proportion of these trap locations, additional monitoring will also be undertaken using camera traps. This two-technique approach is useful as the cameras could detect a cat that approaches a trap, but does not get captured. Along with the data on cat removals, information from these extra detections can be used to update the estimate of the residual population size.

The dataset san_nic_rem contains example data from an eradication phase conducted over 10 separate (primary) periods. During each period, the leg-hold traps and cameras were operated for 20 consecutive nights (secondary periods). As well as the leg-hold trap locations, the dataset also contains information on the IDs of the traps where the cameras were placed. Again, since this is a simulated dataset, we know exactly how many cats remained following the eradication effort (9 cats).

The eradicate package contains a several functions for use with removal type data. The function remPois() can be applied to removal data for the particular case when the latent abundance is assumed to follow a Poisson distribution. remPois also assumes that the population is closed over the primary periods. There are two other functions that employ generalized removal models, remGR() and remGRM() where the former is applied to removal only data while the latter is applied to the joint removal and additional detection data. These models are more flexible than remPois() as abundance can be modeled as a negative binomial distribution as well as the Poisson. In addition, the models also allow for a specific form of lack of closure over the sampling period by allowing for temporary emigration. This allows for the estimation of the proportion of the population that is available for detection during each sampling period (availability). We start by collating the data and extracting habitat info for use in the models as before.

rem<- san_nic_rem$rem
ym<- san_nic_rem$ym # detections from additional monitoring
nights<- san_nic_rem$nights
traps<- san_nic_rem$traps

traps_sf<- st_as_sf(traps, coords=c(1,2), crs=st_crs(region))
traps_buff<- st_buffer(traps_sf, dist=500)
pgrass<- terra::extract(habitat, vect(traps_buff), fun=mean, na.rm=TRUE)
names(pgrass)<- c("id","pgrass")
site.data<- cbind(pgrass, traps)

# Poisson abundance 
emf<- eFrameR(rem, siteCovs = site.data)
r1<- remMN(~pgrass, ~1, data=emf)
summary(r1)
nr1<- calcN(r1)
nr1$Nhat
nr1$Nresid

# Generalized Poisson
emf<- eFrameGR(rem, siteCovs = site.data)
r2<- remGR(~pgrass, ~1, data=emf)
summary(r2)
nr2<- calcN(r2)
nr2$Nhat
nr2$Nresid

# Including additional camera monitoring data in ym 
emf<- eFrameGRM(rem, ym, siteCovs = site.data)
r3<- remGRM(~pgrass, ~1, ~1, data=emf)
summary(r3)
nr3<- calcN(r3)
nr3$Nhat
nr3$Nresid

Using data on extra detections has resulted in a higher estimate of the initial population size compared with the model without extra detections. This is due to the fact that cats were detected in a camera in the periods following the last capture of a cat in the leg-hold traps. The calcN() function also produces an estimate of the residual population remaining, which is simply the difference between the estimated initial population and the number removed. This indicates 18 cats are likely still extant using inference from the removal + extra detections model compared with nominal eradication using the removal only model. The estimate of 18 is somewhat higher than the true number of cats remaining (9), but this is better than the inference of complete eradication using the removal only models. Again, this highlights the advantages of using extra monitoring data for inference.

Finally, eradicate also contains a function for estimating the pre-removal and residual population size for data collected from a single site. This might be applicable when the region of interest has no obvious spatial heterogeneity that can be modeled with covariate information using the other functions (i.e. toad eradication from a small lake or pond). This model remGP(), employs the catch-effort maximum likelihood model of Gould & Pollock (1997). It can also make use of additional monitoring data for joint inference. We can run the remGP() model with the current data by pooling the captures and trapnight data from all traps.

catch<- apply(rem,2,sum)
effort<- rep(nrow(rem), length(catch))
# extra monitoring/effort data
index<- apply(ym,2,sum)
ieffort<- rep(nrow(ym), length(index))

emf<- eFrameGP(catch, effort, index=index, ieffort=ieffort)
mod<- remGP(emf)
summary(mod)
nce<- calcN(mod)
nce$Nhat
nce$Nresid