Modelling density surfaces of intraspecific classes using camera trap distance sampling

Spatially explicit densities of wildlife are important for understanding environmental drivers of populations, and density surfaces of intraspecific classes allow exploration of links between demographic ratios and environmental conditions. Although spatially explicit densities and class densities are valuable, conventional design‐based estimators remain prevalent when using camera‐trapping methods for unmarked populations. We developed a density surface model that utilized camera trap distance sampling data within a hierarchical generalized additive modelling framework. We estimated density surfaces of intraspecific classes of a common ungulate, white‐tailed deer Odocoileus virginianus, across three large management regions in Indiana, United States. We then extended simple statistical theory to test for differences in two ratios of density. Deer density was influenced by landscape fragmentation, wetlands and anthropogenic development. We documented class‐specific responses of density to availability of concealment cover, and found strong evidence that increased recruitment of young was tied to increased resource availability from anthropogenic agricultural land use. The coefficients of variation of the total density estimates within the three regions we surveyed were 0.11, 0.10 and 0.06. Synthesis and applications. Our strategy extends camera trap distance sampling and enables managers to use camera traps to better understand spatial predictors of density. Our density estimates were more precise than previous estimates from camera trap distance sampling. Population managers can use our methods to detect finer spatiotemporal changes in density or ratios of intraspecific‐class densities. Such changes in density can be linked to land use, or to management regimes on habitat and harvest limits of game species.


| INTRODUC TI ON
Ecologists and wildlife managers require estimates of population density for more comprehensive ecological knowledge and effective application of management (Williams et al., 2002). For many species, density estimates specific to age, sex or other intraspecific classes can provide additional information to ecologists and managers. For instance, class-specific densities may aid management of game species when classes are unequally harvested to achieve population goals (Keehner et al., 2015), or of species for which recruitment is associated with food availability (Costello et al., 2003).
Camera traps are commonly used to estimate wildlife density (Delisle et al., 2021). When monitoring intraspecific classes, camera traps are especially advantageous because they allow viewers of images to carefully consider class membership of individuals.
Other methods requiring the physical presence of surveyors in the field may not facilitate such careful determination of class membership. Thus, inaccuracies may arise when surveyors must determine class membership in challenging field conditions (Gerrodette & Forcada, 2005). Difficulties may be especially great when class membership must be assigned quickly before sampled individuals are out of sight, as often happens when sampling extremely mobile species, species that flee in response to surveyors or when surveyors are aboard a moving platform.
Several methods exist for estimating density of unmarked populations using camera traps (Gilbert et al., 2021). Although each method has unique advantages, camera trap distance sampling remains appealing because it rests upon a well-established and intuitive statistical foundation (Buckland et al., 2001;Howe et al., 2017).
However, one limitation of conventional distance sampling is that a single design-based density, estimated from an average count across sampled locations, is inferred across entire study areas (see Buckland et al., 2016 for more on design-based estimation of density). This approach can obscure local fluctuations in density along environmental gradients (e.g. topography, food proximity). Accordingly, ecologists are frequently interested in relationships between density and local covariates.
Density surface modelling uses distance sampling data to estimate spatially explicit density in two steps (Miller et al., 2013).
Firstly, detectability is modelled with a detection function fit to distance sampling data, which accounts for decreased detection probability associated with increasing distance from the surveyor.
Secondly, counts on sampling locations are modelled as a function of spatially explicit predictors. Generalized additive modelling is often used within the second step, which facilitates nonlinear relationships that are common in ecology (Wood, 2017). The fitted model is then used to predict density in unsampled areas, which yields more accurate density estimates across study areas, especially for species whose density relates to spatial fluctuations of habitat characteristics (Miller et al., 2013).
Within the second step of density surface modelling, hierarchical generalized additive models are uniquely appealing for researchers interested in intraspecific classes. Specifically, they can simultaneously model multiple biologically determined subsets of data while maintaining global relationships exhibited by the entire dataset (Pedersen et al., 2019). Such fitted relationships, termed factor-smooth interactions, can enable unique spatial modelling of class-specific densities in a single model. Class specificity is especially useful when class density exhibits unique relationships with spatially explicit predictors. Considering class-specific relationships may yield more accurate estimates of overall density and betweenclass ratios of density, more robust understanding of ecological relationships and better fit of models.
Although the advantages of camera trap distance sampling, density surface modelling and hierarchical generalized additive models are evident, to our knowledge, no study has combined these approaches to model the spatially explicit density of classes. We combine these tools to model the class density of a common North American game ungulate, white-tailed deer Odocoileus virginianus, across three large regional areas over which management is implemented in Indiana, United States. We then use density estimates of classes to compare recruitment rates and adult sex ratios in these three areas. First, because recruitment of white-tailed deer is directly linked to food availability (McCullough, 1979) and deer readily forage on agricultural crops (Rouleau et al., 2002), we hypothesized that recruitment rates of deer would be higher in areas with more agricultural land use. Specifically, we predicted larger ratios of young deer:adult female deer (fawns and does, respectively) densities in areas where the fraction of land being used for agriculture is larger.
Second, annual deer harvest reports from the years and areas we sampled (Boggess & Vaught, 2021;Caudell & Vaught, 2019 showed that ratios of harvested adult male (buck) deer:harvested does were similar in two of the areas we sampled, but higher in the third area we sampled that had the highest relative amount of land devoted to agriculture. Therefore, we predicted that the ratios of buck:doe densities would be lowest in the most highly agricultural region.
Within RMUs, we sampled in randomly selected 41.44-km 2 cells from the deer-harvest reporting grid that the Indiana Department of Natural Resources uses to collect spatially explicit harvest data. We did not preferentially sample 41.44-km 2 cells that we believed would have more or less deer. Specifically, we sampled inside 10.36-km 2 areas (subareas) that we placed inside the larger cell. We attempted to place subareas within larger cells to ensure property permission was homogeneously distributed, and habitat composition within subareas was representative of the larger cell. In some instances, private property access limited our ability to place subareas in perfectly representative locations. During each year, we sampled 20 subareas across all three RMUs. Thus, in 2 years, we sampled seven subareas in each RMU, and in 1 year, we sampled six subareas. We repeatedly sampled two subareas in each RMU during each year to assess interannual variation uncoupled from spatial variation, resulting in sampling of 48 unique subareas. We acquired permission to sample on all properties before we conducted fieldwork.

| Data collection and analysis
We randomly deployed Browning Strike Force HD or BTC-5HDE motion-triggered camera traps in forest, grassland, wetland and agricultural fields within subareas in each RMU using ArcMap 10.7. We placed the two models of camera traps in cover types randomly (i.e. one model was not placed more often in a particular cover type than another model). In forests, we affixed cameras to trees at 1-m height and oriented cameras north to avoid sun glare at dawn and dusk.
In non-forested areas, we affixed cameras to metal posts (that we hammered into the ground) at the same height and orientation. We used a height of 1 m for two reasons: (1) 1 m was low enough to capture even smaller deer at very close distances (Swihart et al., 1998); differed slightly each year due to logistical constraints during camera deployment. Before each sampling window, cameras were deployed for 7 days minimum to allow for deer habituation to cameras.
In rare instances (5%), cameras did not sample the entire 2 weeks because of unpredictable circumstances in the field (e.g. human moving camera). When triggered, cameras captured a three-photo burst with a 0.3-s delay between photos and a minimum 1-s delay before subsequent bursts. In some instances (~25%), settings were slightly altered accidentally or because BTC-5HDE cameras did not permit minimum delays <5 s. We recorded the spatiotemporal effort of each camera as in Howe et al. (2017). When cameras did not sample the entire 2-week window, or when photo delays were set incorrectly, we adjusted the spatiotemporal sampling effort accordingly. If we deployed cameras in subareas that we sampled each year, we treated repeatedly sampled points as a single spatial replicate and combined the spatiotemporal sampling effort across all years.
We did not adjust the sampling effort of cameras when distances to reactive deer were removed because nonreactive deer were frequently detected in the same photo as reactive deer.
We used camera trap distance sampling to model the observation process (Howe et al., 2017). We recorded reference videos as in Howe et al. (2017), which we used to measure the distances between cameras and deer within each image from 1 to 15 m integer F I G U R E 1 Land cover types within regional management units 3 (westcentral), 4 (southern) and 9 (north-eastern) in Indiana, United States. We deployed cameras within 10.36-km 2 subareas. distances. We did not measure distances to bedded deer, or deer that became interested in the camera. Each deer was recorded as doe, buck, fawn or unknown. We binned distances between cameras and deer at 0, 1, 2, 3, 4, 5, 6, 7, 8, 10, 12 and 15 m and estimated several detection functions with the Distance package in R (Miller, 2022; R Core Team, 2022). We did not consider distances >15 m. Candidate detection functions included half-normal key functions with 0, 1 and 2 Hermite polynomial adjustments; uniform key functions with 1 and 2 cosine adjustments; and hazard-rate key functions with 0, 1 and 2 cosine adjustments. We also considered half-normal and hazard-rate key functions with the following factor covariates: (1) night or day (determined by camera flash), (2) microhabitat surrounding the camera (corn or soybean field; deciduous, mixed or evergreen forest; woody or herbaceous wetland; herbaceous grassland; and pasture/hay), (3) RMU, (4) RMU and night or day, (5) night or day and microhabitat, (6) RMU and microhabitat and (7)  Density surface modelling assumes that all individuals are available for sampling. Because deer are not available for camera sampling when bedded, we estimated the activity level of deer to account for this. We recorded the times individual deer were first detected by cameras upon entering the field of view. Because deer are crepuscular (Beier & McCullough, 1990), we double-anchored detection times with the average sunset and sunrise times during our sampling (Vazquez et al., 2019). We estimated the activity level of deer in each RMU by fitting circular kernel densities to doubleanchored detection times, and standard errors of activity levels with nonparametric bootstrapping using the activity package in R (Rowcliffe, 2021;Rowcliffe et al., 2014).
After selecting detection functions and estimating deer activity levels, we fit density surface models using the Dsm package in R (Miller, Rexstad, et al., 2022) in the form: where n i = the count of deer at camera i, 0 = the intercept, f m = the smooth functions of spatially explicit predictors x im , fac = any factor variables considered and v i = the product of detection probability and activity level of deer at camera i used as an offset. We modelled smooths with thin plate regression splines (Wood, 2003). In preliminary analysis, overfitting was apparent (e.g. extremely wiggly relationships). To prevent overfitting, we specified the gamma parameter at 2 (Wood, 2017 section 4.6.1). We found goodness of fit (via observed vs. expected counts) of the quasi-Poisson distribution to outperform other count distributions (e.g. Poisson, negative binomial, Tweedie). Therefore, we used this distribution for all models.
We tested several spatially explicit predictors of density at each camera including metrics of distance and landscape composition or structure within buffers around cameras. For distance metrics, we tested the distance to wetland, which we calculated in R using the 2019 National Land Cover Database land cover raster (Dewitz & Geological Survey, 2021). We used the lanDscapemetrics package (Hesselbarth et al., 2019) in R to test several landscape composition or structure indices within buffers, including the contagion index, coefficient of variation of the core area of patches and total area of concealment cover (defined as forest or wetland) and wetland (McGarigal & Marks, 1995). We used the Indiana primary and secondary roads state-based shapefile (US Census Bureau, Department of Commerce, 2018) and R to calculate the total road length within buffers. Other metrics were tested and eliminated from consideration due to concurvity (a measure similar to collinearity for smooth models; Wood, 2017) with these better predictors. We evaluated indices in buffers with radii of 250, 750, 1425, 4000 and 8000 m (see Supporting Information for buffer choices). We chose a final buffer size for each metric based on the strength of the relationship with deer counts. We did not use the same metrics with multiple buffer radii due to concurvity.
For each predictor, we used observed versus expected counts to choose between the following types of class-smooth interactions for each predictor: (1) single common smoother across all classes (i.e. global smoother; no factor-smooth interaction); (2) global smoother and class-specific smoothers with identical wiggliness; (3) global smoother and class-specific smoothers with differing wiggliness; (4) class-specific smoothers with identical wiggliness but no global smoother; and (5) class-specific smoothers with differing wiggliness but no global smoother (Pedersen et al., 2019). Additionally, we considered habitat type (forest, wetland, grassland and agricultural field) and an interaction between class and RMU as factor variables. Upon deciding which smooth type to fit for each predictor, we fit a global model containing all the best factor-smooth types for each predictor. We then used F-tests to identify and remove factor-smooths in the global model with weak (p > 0.05) relationships with deer counts.
To predict spatially explicit class density, we created a grid over each RMU. We specified the resolution of grid cells to be 30 × 30 m because we did not expect deer density to change perceptibly over this area. Within each grid cell, we calculated each predictor of deer density that we parameterized in our final density surface model, and used the final fitted model to predict buck, doe, fawn and unknown density. We did not predict density in developed, barren (rock, sand or clay) or scrub/shrub habitats because we did not sample these habitats. Additionally, we used the Dsmextra package in R (Bouchet et al., 2020) to calculate Euclidean and Mahalanobis distances (Mesgaran et al., 2014) and identify cells with covariates exhibiting univariate or combinatorial extrapolation outside the range of the covariates we sampled. We estimated densities that both considered and did not consider extrapolated cells (extrapolated and non-extrapolated densities respectively).
We used a modified formula from Gerrodette and Forcada (2005) to prorate unknown-specific density in each cell to the known classes: where i refers to known class i (bucks, does or fawns), u refers to the unknown class, k refers to all other classes ≠ i, D i (pr) = the prorated density of class i, D i = the density of class i before proration, D u = the density of the unknown class and w i = the fraction of class i that is at risk of being unidentifiable. For does and fawns, w i = 1. However, because only a fraction of bucks had shed their antlers during our sampling, we specified w i for bucks as the fraction of bucks detected that shed both antlers in the RMU the cell was within. We considered bucks with antlers to be wholly identifiable, as antlers are easily distinguishable in images. To estimate the final density of intraspecific classes in each RMU and habitat-specific densities, we averaged the predicted class density across the entire grid of each RMU and across habitats within each RMU respectively. We estimated total deer densities (all classes combined) by summing the predicted density of all classes in each cell, and averaging the sum across cells.
Because of the size of our prediction grids (~4 to ~16 million cells), we used posterior simulation in combination with Welford's online algorithm (Knuth, 2014;Welford, 1962) to estimate the vari- 2. Use sampled parameters to generate class-specific predictions of density across the prediction grid.
3. Save necessary summary statistics needed to calculate the variance per Welford's algorithm.
We repeated the above algorithm 1000 times and prorated the simulated unknown densities using equation 2. We estimated variances of all density estimates from the 1000 simulated densities (using Welford's algorithm) and propagated variance from the activity level via the delta method. We approximated 95% confidence intervals using the percentile method (Efron, 1981). Cov D (f),D(d) = the covariance of D (f) and D (d) estimated from the 1000 simulated densities. We repeated this process for buck:doe ratios. To infer differences in ratios between RMUs, we used a Z test: where rat 1 and rat 2 = the two ratios being compared, and Var rat 1 − rat 2 = Var rat 1 + Var rat 2 . We implemented a one-tailed test in instances when our hypotheses on ratios of density between RMUs were directional.

| RE SULTS
We deployed 1295 cameras but removed 73 cameras from analysis because of faulty placement. After accounting for repeatedly sampled points across years, we used data from 1018 independent locations and measured 83,824 distances ( Table 1).

We observed a spike in detections at close distances in open
areas. Therefore, we removed all hazard rate models from consideration in open areas, because this model can fit unnatural spikes at close distances. In this instance, the hazard rate key function can model unnaturally abrupt declines in detection probability, which underestimates detectability and overestimates density (Buckland et al., 2001). After fitting each type of factor-smooth interaction for each predictor, goodness of fit via observed versus expected counts suggested that class-specific smoothers with the same wiggliness without a global smoother were best for the coefficient of variation Var of the core area of patches <4000 m, and a global smoother in addition to class-level smoothers with differing wiggliness was best for the contagion index <4000 m, distance to wetland, area of wetland <8000 m, area of concealment cover <1425 m and total length of road <8000 m. F-tests on a global model containing all the best factor-smooth interactions suggested that factor-smooth interactions were unnecessary for the coefficient of variation of the core area of patches <4000 m, contagion index <4000 m, distance to wetland, area of wetland <8000 m and total length of road <8000 m. Therefore, we replaced these factor-smooth interactions with global smooths, but kept the factor-smooth interaction for the area of concealment cover <1425 m within our final model (deviance explained = 48.7%).
We found strong evidence to include factor terms for habitat Extrapolated densities were much more variable (Supporting Information). Henceforth, we only present non-extrapolated densities. Within each RMU, density was highest in wetlands and lowest in agricultural fields ( Table 2). The CV of total densities in RMU 3, 4 and 9 was 0.11, 0.10 and 0.06 respectively. The average classand habitat-specific CV was 0.15 (SE = 0.02) and 0.14 (SE = 0.03) respectively. Buck:doe ratios between RMUs and fawn:doe ratios in RMU 3 and 9 (observed difference = 0.10, Z = 0.83, p = 0.203) did not exhibit strong differences (Table 3). However, the fawn:doe ratio in RMU 3 was significantly larger than RMU 4 (observed difference = 0.19, Z = 1.70, p = 0.045) and RMU 9 was marginally larger than RMU 4 (observed difference = 0.09, Z = 1.38, p = 0.084).

| DISCUSS ION
We present a novel strategy to estimate and test for differences between density surfaces of classes within a population. We first implemented a density surface model with camera trap distance sampling data, which linked variation in local density to

F I G U R E 2 Probability density functions of distances (m) between white-tailed deer and camera traps in Indiana, United
States. In open areas, the uniform key function with 1 cosine adjustment term is presented. In concealed areas, the uniform key function with 2 cosine adjustment terms is presented.

F I G U R E 3
Activity probability distributions of white-tailed deer in regional management units (RMU) 3, 4 and 9 in Indiana, United States. Solid lines represent the fit kernel density, and dotted lines represent 95% confidence intervals.
environmental predictors across large spatial expanses with widely differing landscape characteristics. Secondly, the framework of hierarchical generalized additive modelling allowed fitting and assessment of various factor-smooth interactions, thereby illuminating differential class-specific responses to external factors within a single statistical model. Thirdly, we extended foundational statistical theory to tractably test for differences in ratios of density, which further elucidated differences in ratios of classes related to landscape characteristics.
We found metrics associated with landscape fragmentation, anthropogenic development and concealment cover to predict local density. These predictors can be tied to naturally or anthropogenically sourced increases in food availability. Density was largest in areas with moderately high contagion values, which indicates moderately low amounts of landscape fragmentation. Similarly, density responded positively to increased road lengths, which can increase edges that contain greater concentrations of plants that deer forage (Ford et al., 1993). We extended existing statistical theory to provide a tractable method to compare ratios of density. Using this method, we found F I G U R E 4 Partial effects plot of our final density surface model fit to camera trap distance sampling data from white-tailed deer within Indiana, United States. Predictors include the contagion index of habitat <4000 m of the camera (a), coefficient of variation of the core area of patches <4000 m of the camera (b), distance (m) from the camera to the nearest wetland (c), amount (m 2 ) of wetland <8000 m of the camera (d), amount (m 2 ) of concealment cover <1425 m of the camera (e) and length (m) of road <8000 m of the camera (f). For the amount of concealment cover <1425 m of the camera, we used a factor-smooth interaction that implemented a global smoother and class-level smoothers with differing wiggliness.
F I G U R E 5 Predicted densities (deer/ km 2 ) of white-tailed deer across regional management units (RMU) 3, 4 and 9 of Indiana, United States, from a density surface model fit to camera trap distance sampling data. Cells exhibiting univariate or combinatorial extrapolation are not depicted.
F I G U R E 6 Coefficient of variation (CV) of predicted densities of white-tailed deer across regional management units (RMU) 3, 4 and 9 of Indiana, United States, from a density surface model fit to camera trap distance sampling data.
statistically higher recruitment rates in areas with more agricultural use, which supported our hypothesis. Although agriculture is one of the greatest causes for biodiversity loss (MEA, 2005), many adaptable species, such as deer, consume crops for nutritional gains (Putman & Moore, 1998). In these species, supplemental caloric intake from agriculture may support unnaturally high recruitment rates for dense populations, even when natural food resources are limited.
Extrapolation beyond the range of covariates sampled in the field can lead to dubious inference in ecological studies, as it assumes that the form of the fitted model remains appropriate (Jones et al., 2022).
We exemplified this by using Euclidean and Mahalanobis distances to identify cells in our prediction grids that exhibited univariate or combinatorial extrapolation (Mesgaran et al., 2014), and computed densities that included and did not include these cells. Although extrapolation did not appear to affect the point estimates of density in two of the RMUs, the variances of extrapolated densities were considerably higher, and larger spatially explicit CVs usually corresponded to extrapolation ( Figure 6). Conversely, point estimates of density were affected in RMU 9 by including extrapolated cells. Therefore, we encourage readers to interpret extrapolated and non-extrapolated densities accordingly; although the total nonextrapolated density in RMU 9 was smaller than RMU 4, cells exhibiting extrapolation in RMU 9 could alter this conclusion to an unknown degree. More thorough sampling across the range of environmental predictors of density that we used in our density surface model would be needed to minimize such sampling-based discrepancies and should be a design consideration for future studies that intend to use the methods we describe.
Wildlife research and management benefits from precise density estimates (Williams et al., 2002). The CV typically is used to assess relative precision of estimates (Skalski et al., 2010). In conventional camera trap distance sampling, the encounter rate variance is   Cappelle et al., 2021). Density surface modelling avoids design-based estimation of the encounter rate variance by modelling encounter rates at cameras as a function of environmental predictors. This difference in modelling approaches may be partially responsible for the lower CVs of our density estimates, as CVs from density estimates using conventional camera trap distance sampling on the same data were 1.34 (SE = 0.15) times larger than those from density surface modelling. Therefore, we believe future managers using camera traps to sample wildlife will benefit from using spatially explicit methods such as ours, as this will facilitate the detection of meaningful changes in density and provide confidence in single estimates.

CO N FLI C T O F I NTE R E S T S TATE M E NT
The authors have no conflict of interests to disclose.

PEER R E V I E W
The peer review history for this article is available at https:// w w w.web of scien ce.com/api/g atew ay/wos/p e er-revie w/ 10.1111/2041-210X.14093.

DATA AVA I L A B I L I T Y S TAT E M E N T
The data, files and code used to conduct our analysis as well as estimate deer density using conventional camera trap distance sampling have been archived at the Purdue University Research Repository   Table S1. Habitat-specific, class-specific, and total density estimates (deer/km 2 ) of white-tailed deer (Odocoileus virginianus) in Regional Management Units (RMU) 3, 4, and 9 of Indiana, USA. Densities were estimated with a density surface model fitted to camera trap distance sampling data collected in the winter of 2019, 2020, and 2021. Habitat-specific densities are presented for agricultural fields (Ag), grasslands and pastures (Grass), wetlands (Wet), and forests (For). Densities were estimated across all areas, including areas that exhibited combinatorial or univariate extrapolation past the range of environmental covariates that we sampled. UCI = upper 95% confidence interval. LCI = lower 95% confidence interval.