# Analysis of camera trapping data

Example analysis with Ivory Coast Maxwell’s duiker.

Eric Howe and Eric Rexstad http://distancesampling.org (CREEM, Univ of St Andrews)https://creem.st-andrews.ac.uk
2022-11-03

# Analysis of camera trapping data using distance sampling

A distance sampling approach to the analysis of camera trapping data offers the potential advantage that individual animal identification is not required. However, accurate animal-to-camera detection distances are required. This requires calibration prior to the survey with images of objects of known size taken at known distances from the camera. See details in Howe, Buckland, Després-Einspenner, & Kühl (2017) for description of the field work and data analysis. Here we present analysis of data from Howe et al. (2017) using the R package Distance Miller, Rexstad, Thomas, Marshall, & Laake (2019).

## Estimating temporal availability for detection

Heat- and motion-sensitive camera traps detect only moving animals within the range of the sensor and the field of view of the camera. Animals are therefore unavailable for detection by camera traps when they are stationary, and when they are above (e.g., semi-arboreal species) or below (e.g., semi-fossorial species) the range of the sensor or the camera, regardless of their distance from the camera in two dimensions. This temporally limited availability for detection must be accounted for to avoid negative bias in estimated densities. When data are abundant, researchers may choose to include only data from times when 100% of the population can be assumed to be active within the vertical range of camera traps Howe et al. (2017). However, for rarely-detected species or surveys with lower effort, it might be necessary to include most or all observations of distance. In these situations, survey duration ($$T_k$$) might be 12- or 24-hours per day, and it becomes necessary to estimate the proportion of time included in $$T_k$$ when animals were available for detection. Methods for estimating this proportion directly from CT data have been described Rowcliffe, Kays, Kranstauber, Carbone, & Jansen (2014), and it can be included in analyses to estimate density Bessone et al. (2020), for example as another multiplier, potentially with an associated standard errors.

## Data input

Times of independent camera triggering events for the period 28 June 21 September 2014 at 23 cameras are recorded in a file described in the data repository Howe, Buckland, Després-Einspenner, Kühl, & Buckland (2018).

trigger.events <- read.table(file="https://datadryad.org/stash/downloads/file_stream/73223",
header=TRUE)

The format of the trigger.events data frame is adjusted to create a datetime field for use in the activity package Rowcliffe (2021)

trigger.events$date <- paste("2014", sprintf("%02i", trigger.events$month),
sprintf("%02i", trigger.events$day), sep="/") trigger.events$time <- paste(sprintf("%02i", trigger.events$hour), sprintf("%02i", trigger.events$minute),
sep=":")
trigger.events$datetime <- paste(trigger.events$date, trigger.events$time) ## Functions in the activity package We will employ two functions from the activity package. First, convert the time of day of a camera triggering event into the fraction of the 24hr cycle when the event took place, measured in radians. In other words, an event occurring at midday is recorded as $$\pi$$ and an event occurring at midnight is recorded as 2$$\pi$$. library(activity) trigger.events$rtime <- gettime(trigger.events$datetime, "%Y/%m/%d %H:%M") With the radian conversion of the camera triggering times, the distribution of the triggering events times is smoothed, using a kernel smoother by the function fitact. The function estimates the proportion of time (in a 24hr day) animals were active. In addition, the triggering time data can be resampled to provide a measure of uncertainty in the point estimate of activity proportion. act_result <- fitact(trigger.events$rtime, sample="data", reps=100)

A plot of the histogram of triggering times, along with the fitted smooth is provided by a plot function applied to the object returned by fitact.

plot(act_result)

The value computed by the smooth through the activity histogram can be extracted from the object created by fitact. The extraction reaches into the object to look at the slot called act. The uncertainty around the point estimate is derived from resampling that takes place within fitact. The slot will display the point estimates, standard error and confidence interval bounds.

print(act_result@act)
       act         se   lcl.2.5%  ucl.97.5%
0.33463831 0.02266198 0.28627336 0.38240833 

The output above would be used to adjust density estimates for temporal activity if the cameras were in operation 24hrs per day. However, in this study, cameras were only active for 11.5 hours per day (0630-1800).

We use the temporal availability information to create a multiplier. Our multiplier must be defined as > proportion of the camera operation time animals were available to be detected

This is not equivalent to the value produced by the fitact function; that value is the proportion of 24hr animals were available to be detected. The availability multiplier must be adjusted based on the daily camera operation period. Uncertainty in this proportion is also included in our computations.

The point estimate and standard error are pulled from the fitact object, adjusted for daily camera operation time and placed into a data frame named creation in a named list, specifically in the fashion shown.

camera.operation.per.day <- 11.5
prop.camera.time <- camera.operation.per.day / 24
avail <- list(creation=data.frame(rate = act_result@act[1]/prop.camera.time,
SE   = act_result@act[2]/prop.camera.time))

A more robust way of incorporating uncertainty in the temporal availability estimate will be described later.

# Detection data analysis

Detection distances for the full daytime data set is also available on Howe et al. (2018) and is read in the code chunk below:

DuikerCameraTraps <- read.csv(file="https://datadryad.org/stash/downloads/file_stream/73221",
DuikerCameraTraps$Area <- DuikerCameraTraps$Area / (1000*1000)
DuikerCameraTraps$object <- NA DuikerCameraTraps$object[!is.na(DuikerCameraTraps$distance)] <- 1:sum(!is.na(DuikerCameraTraps$distance))

Data file recorded study area size in square meters; second line above converts this to area in square kilometers; the remaining lines create an object field, which uniquely identify each observation.

## Exploratory Data Analysis

A quick summary of the data set including: How many camera stations and how many detections in total.

sum(!is.na(DuikerCameraTraps$distance)) [1] 11180 table(DuikerCameraTraps$Sample.Label)

A1   A2   A3   A4   B1   B2   B3   B4   C1   C2   C3   C4   C5   C6
388   66  988  420    3 1951   73  208   52  195  767  153   41 2682
D3   D4   D5   E3   E4   E5   E6
342  193  524  518    1  375 1241 

Note, three sampling stations (B1, C5, E4) had no detections. The one record for each of those stations has distance recorded as NA, but the record is important because it contains effort information.

## Distance recording

A quick examination of the distribution of detection distances; note the bespoke cutpoints causing distance bins to be narrow out to 8m, then increasing in width to the maximum detection distance of 21m.

breakpoints <- c(seq(0,8,1), 10, 12, 15, 21)
hist(DuikerCameraTraps$distance, breaks=breakpoints, main="Peak activity data set", xlab="Radial distance (m)") ## Truncation decisions As described by Howe et al. (2017): a paucity of observations between 1 and 2 m but not between 2 and 3 m, so we left-truncated at 2 m. Fitted detection functions and probability density functions were heavy-tailed when distances >15 m were included, so we right truncated at 15 m. ## Detection function fits The conversion factor must be included both in the call to ds() and the call to bootdht(). Candidate models considered here differ from the candidate set presented in Howe et al. (2017). This set includes • uniform key with 1, 2 and 3 cosine adjustments, • half normal key with 0, 1 and 2 cosine adjustment and • hazard rate key with 0, 1 simple polynomial adjustments. The maximum number of parameters in models within the candidate model set is no more than 3. library(Distance) trunc.list <- list(left=2, right=15) mybreaks <- c(seq(2,8,1), 10, 12, 15) conversion <- convert_units("meter", NULL, "square kilometer") uni1 <- ds(DuikerCameraTraps, transect = "point", key="unif", adjustment = "cos", nadj=1, convert_units = conversion, cutpoints = mybreaks, truncation = trunc.list) uni2 <- ds(DuikerCameraTraps, transect = "point", key="unif", adjustment = "cos", nadj=2, convert_units = conversion, cutpoints = mybreaks, truncation = trunc.list) uni3 <- ds(DuikerCameraTraps, transect = "point", key="unif", adjustment = "cos", nadj=3, convert_units = conversion, cutpoints = mybreaks, truncation = trunc.list) hn0 <- ds(DuikerCameraTraps, transect = "point", key="hn", adjustment = NULL, convert_units = conversion, cutpoints = mybreaks, truncation = trunc.list) hn1 <- ds(DuikerCameraTraps, transect = "point", key="hn", adjustment = "cos", nadj=1, convert_units = conversion, cutpoints = mybreaks, truncation = trunc.list) hn2 <- ds(DuikerCameraTraps, transect = "point", key="hn", adjustment = "cos", nadj=2, convert_units = conversion, cutpoints = mybreaks, truncation = trunc.list) hr0 <- ds(DuikerCameraTraps, transect = "point", key="hr", adjustment = NULL, convert_units = conversion, cutpoints = mybreaks, truncation = trunc.list) hr1 <- ds(DuikerCameraTraps, transect = "point", key="hr", adjustment = "poly", nadj=1, convert_units = conversion, cutpoints = mybreaks, truncation = trunc.list) We do not present the density estimates produced from the fitted detection function models because a) we have not chosen a preferred model and b) the density estimates have not been adjusted for viewing angle and temporal availability. ## Model selection adjustments from overdispersion Overdispersion causes AIC to select overly-complex models, so analysts should specify the number/order of adjustment terms manually when fitting distance sampling models to data from camera traps, rather than allowing automated selection using AIC. Howe, Buckland, Després-Einspenner, & Kühl (2019) describes two methods for performing model selection of distance sampling models in the face of overdispersion. Here we provide R functions to perform the first of these methods. The first method of Howe et al. (2019) employs a two-step process. First, an overdisersion factor $$(\hat{c})$$ is computed for each key function family from the most complex model in each family. The $$\hat{c}$$ is derived from the $$\chi^2$$ goodness of fit test statistic divided by its degrees of freedom. This results in an adjusted AIC score for each model in the key function family: $QAIC = -2 \left \{ \frac{log(\mathcal{L}(\hat{\theta}))}{\hat{c}} \right \} + 2K$ Code to perform this QAIC computation is found in the function QAIC in the Distance package: Tables of QAIC values for each key function family are shown below (code for kable() calls suppressed for easier readability of results). Table 1: QAIC values for uniform key models. df QAIC uni1 1 2819.103 uni2 2 2820.611 uni3 3 2817.732 Table 2: QAIC values for half normal key models. df QAIC hn0 1 2294.191 hn1 2 2290.255 hn2 3 2292.093 Table 3: QAIC values for hazard rate key models. df QAIC hr0 2 2465.119 hr1 3 2466.395 From this first pass of model selection based on QAIC values, we find the model with the uniform key function preferred by QAIC has three cosine adjustment terms. The preferred model from the half normal key function family has one cosine adjustment term. Finally, the preferable model from the hazard rate key function family has no adjustment terms. The second step of model selection ranks the models by their $$\hat{c}$$ values. chats <- chi2_select(uni3, hn1, hr0)$criteria
modnames <- unlist(lapply(list(uni3, hn1, hr0), function(x) x$ddf$name.message))
results <- data.frame(modnames, chats)
Dhat = ests$individuals$D$Estimate)) } ## Multiplier function This rather complex list makes use of make_activity_fn that exists in the Distance package used to call the fitact function from the activity package. For the user, your responsibility is to provide three arguments to this function: • vector containing the detection times in radians (computed in earlier section), • the manner in which precision of the temporal availability estimate is produced and • the number of hours per day the cameras are in operation mult <- list(availability= make_activity_fn(trigger.events$rtime, sample="data",
detector_daily_duration=camera.operation.per.day))

## Remaining arguments to bootdht

Just as with dht2 there are arguments for the model, flatfile, sample_fraction, convert.units and multipliers (although for bootdht multipliers uses a function rather than a single value). The only novel arguments to dht2 are resample_transects indicating camera stations are to be resampled with replacement, nboot for the number of bootstrap replicates and cores which can take advantage of multiple cores in your computer to speed computation.

daytime.boot.uni <- bootdht(model=uni3, flatfile=DuikerCameraTraps,
resample_transects = TRUE, nboot=600, cores = 11,
summary_fun=mysummary, sample_fraction = samfrac,
convert_units = conversion, multipliers=mult)

## Confidence limits computed via the percentile method of the bootstrap.

print(summary(daytime.boot.uni))
Bootstrap results

Boostraps          : 600
Successes          : 577
Failures           : 23

median  mean   se  lcl   ucl  cv
Dhat  17.71 18.66 7.08 7.85 34.99 0.4
hist(daytime.boot.uni$Dhat, breaks = 20, xlab="Estimated density", main="D-hat estimates bootstraps") abline(v=quantile(daytime.boot.uni$Dhat, probs = c(0.025,0.975), na.rm=TRUE), lwd=2, lty=3)

The confidence interval derived from the bootstrap is considerably wider than the confidence interval derived from analytical methods.

Bessone, M., Kühl, H. S., Hohmann, G., Herbinger, I., N’Goran, K. P., Asanzi, P., … Fruth, B. (2020). Drawn out of the shadows: Surveying secretive forest species with camera trap distance sampling. Journal of Applied Ecology, 57(5), 963–974. https://doi.org/10.1111/1365-2664.13602
Howe, E. J., Buckland, S. T., Després-Einspenner, M.-L., & Kühl, H. S. (2017). Distance sampling with camera traps. Methods in Ecology and Evolution, 8(11), 1558–1565. https://doi.org/10.1111/2041-210X.12790
Howe, E. J., Buckland, S. T., Després-Einspenner, M.-L., & Kühl, H. S. (2019). Model selection with overdispersed distance sampling data. Methods in Ecology and Evolution, 10(1), 38–47. https://doi.org/10.1111/2041-210X.13082
Howe, E. J., Buckland, S. T., Després-Einspenner, M.-L., Kühl, H. S., & Buckland, S. T. (2018). Data from: Distance sampling with camera traps. https://doi.org/https://doi.org/10.5061/dryad.b4c70
Miller, D., Rexstad, E., Thomas, L., Marshall, L., & Laake, J. (2019). Distance sampling in r. Journal of Statistical Software, Articles, 89(1), 1–28. https://doi.org/10.18637/jss.v089.i01