Estimating the density of wildlife populations is still a difficult task, especially when you work with spatially open populations and you must relax the assumption of closure, which is the basis of most methods currently used. Further difficulties arise when obtaining density estimates at small spatial scales. Using eight years (1996-2003) to monitor data from a roe deer Capreolus capreolus population that lives in a sub-Mediterranean environment in central Italy, we were able to estimate local density (at a spatial scale of one home range) by using a large sample of radio-marked animals. Local density estimates could be obtained only in zones in which radio-marked deer were available in sufficient numbers. To estimate local density in the whole study area, we developed a calibration model, which allowed us to infer density where radio-marked deer were absent or scarce. To do this, we computed the mark-resight density estimates (using radio-marked animals) and related these estimates to linear and non-linear functions of animal count and surface area of fields, to obtain a set of density estimators. Then, we selected a linear combination of such estimators, whose quality was assessed by cross-validation. Our results show that the method we propose can be effective in investigating small-scale spatial structure of density in a roe deer population. We see several potential applications of this method for both research and management purposes.
Estimation of population abundance is a fast developing issue. Several approaches, such as distance-sampling and capture-mark-recapture (Borchers et al. 2002), are available for closed populations. The estimation of spatially open populations' size is still a difficult task, especially if one is interested in densities rather than abundances.
For spatially closed populations, the estimation of animal numbers and densities is equivalent, since the area size is fixed. However, this is not so for spatially open populations, which are the rule rather than the exception in wildlife management. The use of radio-telemetry may provide a solution to this problem, as shown by Eberhardt (1990) and White & Schenk (2001). Density estimates based on radio-telemetry were calculated successfully for bears Ursus sp. (Miller et al. 1987, 1997), goshawks Accipiter gentilis (Kenward et al. 1981), roe deer Capreolus capreolus (Focardi et al. 2002b, Hewison et al. 2007) and CMR methods were used for the house mouse Mus muscutus (Efford 2004), and distance sampling can also provide population density estimates (Buckland et al. 2001, Focardi et al. 2006).
For spatially open populations such methods are reliable provided that the survey is so fast that one obtains a snapshot estimate (see Kenward et al. 1981 for another approach to relax the assumption of closure). Focardi et al. (2002b) estimated the density of the roe deer population at Tredozio (in central Italy) at a 200-ha scale using mark-resighting and then calculated the average density with a maximum likelihood approach by pooling four count sessions (White 1996). Their assumption was that even if some movement occurred across the border of the area, it was negligible given the small border-area ratio. However, during the counts at Tredozio some radio-marked deer, which were detected outside the area just before the beginning of the survey, were subsequently sighted during the two-hours trial. This problem, probably negligible at a large spatial scale, may become relevant at small scale when the border-area ratios increase.
The importance of investigating ecological patterns at various spatial scales was reviewed by Ray & Hastings (1996). Local density, i.e. the density at a very small scale, is recognised as being more and more important. For roe deer in particular, the importance of local density was pointed out by Pettorelli et al. (2005), who detected an influence of environmental variables on fawn survival on the scale of a single home range. Moreover, local density can be used to estimate large-scale density if necessary, whereas the inverse, of course, does not hold.
Local density estimation is a relatively novel issue with theoretical questions still being unaddressed. One major point is that for local density estimation, the closure assumption does not hold. Each individual, which inhabits the surveyed zone, has a given probability (usually < 1) of staying inside the zone at a given time. Density can be redefined as the average number of individuals staying inside the surveyed zone. Of course, if the closure assumption is valid, i.e. the probability of staying inside the zone is one, this interpretation is completely consistent with the usual definition of density.
Exploiting an idea presented by Kenward et al. (1981), Focardi et al. (2002a) relaxed the assumption of spatial closure and were successful in estimating roe deer density at a very small scale (i.e 10 ha) by combining capture-mark-resighting techniques using radio-telemetry. However, since the size of the count unit in this case was very small, simply by chance many zones contained too few marked animals or none at all. Accordingly, the estimates obtained by Focardi et al. (2002a) were scattered over the study area (24 reliable estimates out of 45 surveyed zones) and were concentrated in the areas with the highest presence of marked animals. This outcome was unfortunate because, due to the lack of local density estimates for the whole study area, we could not analyse the effect of local density on roe deer demography. The aim of our paper is to further develop the work of Focardi et al. (2002a) and to present a method that is amenable in estimating the local density even in areas without marked animals.
Basically, our method originates from sightability models, which have often been used to estimate the size of ungulate populations (see e.g. Poole 2007, Freddy et al. 2004 and White & Schenk 2001 for a summary). In sightability models, the detection probability is computed using a trial survey where detectability of radio-marked animals is modelled as a function of environmental variables in a logistic regression framework. Once sightability is computed, one can derive population size (Borchers et al. 2002).
Central to our work is the idea that individual sightability is influenced by population density, vegetation openness and fragmentation. We suppose that biological and geometrical issues can influence the rate at which deer can be seen at a survey point. An example of a biological reason is that at higher densities, food shortage can force deer to appear in open areas; a geometrical reason could be that the edge of a survey area grows as the square root of its extension. These relations have a non-linear form and we expect other relations of this type to influence deer sightability.
More specifically, we derive a calibration model using local density estimates as dependent variable and linear and non-linear functions of counted animals and survey zone area as regressors.
Material and methods
Study area
Our study area extends over approximately 5 km2 of a hilly landscape in north-central Italy (44°04′3”N, 11°44′30”E) with woods covering 47% of the whole surface. Woods in the area are predominantly coppice (i.e. deciduous oak Quercus spp., hornbeam Ostrya carpinifolia, chestnut Castanea sativa, ash Fraxinus ornus and maple Acer sp.) and coniferous plantations (mostly pine Pinus nigra) with a few cultivated chestnut and old woods. Open vegetation covers 48% of the whole area, consisting mostly of cultivated fields (hay, cereals, sunflower Helianthus annuus and pastures) with a significant portion of abandoned fields.
A stream and an adjoining paved road divide the study area into two subareas. The subareas (Monti on the southeastern side and Collinaccia on the northwestern side) share most of their relevant characteristics for vegetation composition of woods and types of cultivations. Nonetheless, the average open field area is 6 ha for Collinaccia and 11 ha for Monti.
The climate in our study area is favourable for roe deer. The average monthly temperature never drops below 0°C and there was no snow cover during the surveys, as snow cover in March is occasional and short lasting.
Roe deer captures and radio-tracking
Roe deer were captured by driving and the use of falling long-nets during autumn-winter. Captured deer were sexed, aged and fitted with Televilt TXE3 radio-collars. Newborns were captured and fitted with Televilt TXH2 radio-collars during the hiding period (i.e. end of May through 10 June; Raganella-Pelliccioni et al. 2004).
We performed animal handling according to the present regulations related to animal welfare and with constant veterinary assistance. Upon capture all animals were fitted with a soft mask (to reduce the animals' distress) which allowed deer to breathe normally. Deer were positioned on the right side to avoid ruminal meteorism which might reduce respiratory efficiency. A total of 220 roe deer were captured between 1996 and 2003.
Each deer was localised 12 times per month with a time-stratified sampling which was distributed over 24 hours.
Single operators equipped with a receiver (Advanced Telemetry Systems, model R2000 or Wildlife Materials RX 1000 S) collected at least three azimuths per fix. Animal location was computed using software LOCATE (Nams 1990), which estimates the position of the radio source based on a maximum likelihood criterion.
Surveys of roe deer population
We performed surveys at the end of March (i.e. green-up season in the study area) during 1996-2003. The observers were trained operators and were equipped with binoculars and telescopes in order to be able to read the numbers of the collars and, thus, individually identify the marked deer. Counts were performed on four consecutive occasions at dawn and dusk. The observers were located at vantage points (Mayle et al. 1999) such as farmyards, hunting hides and roads. Each observation session lasted for two hours, just after dawn or before dusk. The number of vantage points ranged within 19-29 per year in response to changes in the surveyed area or personnel availability. The operators were located where they could achieve a full visual coverage of all fields. Accordingly, while one field could be monitored by a single operator, sometimes two or more operators could survey the same field and a single operator could cover two or more field portions. Where necessary, the operators communicated with each other using a radio to avoid double counts. With a total of 135 fields, 16-19 fields were surveyed each year.
Statistical methods
Our method consisted of four steps: 1) we computed the mark-resight density estimates to be used as a set of reference densities, 2) we related these estimates to linear and non-linear functions of animal count and surface area of fields to obtain a set of density estimators, 3) we selected a linear combination of such estimators based on a set of statistical criteria and 4) we cross-validated the selected model to test whether the selected model was able to predict densities at three different spatial scales: the entire study area, the subareas and small plots that were the size of a single home range.
Concerning 1), we obtained the reference density estimates at small scale using the method described by Focardi et al. (2002a). The number of deer associated with each field was estimated by mark-resight (White & Schenk 2001) using the software NOREMARK (White 1996). We computed density, correcting for a factor accounting for the fraction of the time that each animal spent inside the field (Kenward et al. 1981) and divided by the area of the field. This way we obtained M reference density estimates representing our sample units, each of which pertained to a given field in a given year.
We obtained subarea densities using the same procedure. The whole area densities were the average of all the subarea densities, weighted for subareas extension.
Concerning 2), we used two variables to account for deer density: the number of counted deer, Oi, and the field area, (Ai), with i = 1,…,M. Oi was computed as the mean number of deer recorded in a given field during the four count sessions. Ai was the area of the horizontal projection of the surveyed field computed using a GIS.
Let us define ρ(O,A) a generic function such that (ρi − Di)2 be minimum, where Di is the estimated density in the ith count zone. There is no a priori reason to select a certain function as ρ. A graphical exploration of the relationships among D, A and O allowed us to select a certain number of candidate models (e.g. linear, exponential and hyperbolic).
The fitting of non-linear functions was made using PROC NLIN of SAS (SAS Institute Inc. 2000). We obtained G different models (ρ1,ρ2,…,ρg,…,ρG) …where each ρg was an estimator of density, D.
All the regressors are some function of Oi and Ai. For a more intuitive understanding of their biological meaning, we put Ei = Oi/Ai and Ui = Ei/Di. Ei is an apparent density, while Ui is the ratio between Ei and the CMR-estimates of density in a given field and, hence, it represents the intensity of use of open fields by roe deer. The regressors we used, the function that was fitted to compute parameter values and the function that we used as regressor are shown in Table 1.
Concerning 3), we then exploited all of these estimators to compute a best linear function, R(ρ1, ρ2, …, ρG). The fitting, computed so that (Ri − Di)2 be minimum, was performed using PROC REG (SAS Institute Inc. 2000). All possible models were evaluated. We obtained one model with G regressors, G models with G-1 regressors, G(G-1)/2 models with G-2 regressors,…, and G models with one regressor. We discarded all of the models that were not full-rank (i.e. at least one regressor was a linear combination of the others), and then we ranked the remaining models according to the Akaike Information Criterion (AIC). Finally, we tested for collinearity among regressors, discarding the models in which at least one regressor had Variance Inflation Factors (VIF) ≥ 10 (Belsey et al. 1980). The model with the best AIC and VIF < 10 for each ρg was selected.
The final model reads: R = s1b1 ρ1+ …+ sgbgρg +…+ sGbGρG + c, where sg = 0 if ρg was discarded, and sg = 1 if ρg entered the model. bg is the coefficient of ρg and c is the intercept.
Concerning 4), we compared Ri with Di to evaluate the performance of the estimator R. We evaluated first the significance of the model used (Fisher F-test). The level of precision attained by the model was described computing the median absolute error, which was reported in percentage of the reference density, D, and its interquartile variation, where D is scaled to 0. Then we tested for the presence of a significant bias of the predictor (Student's t-test), comparing model estimates with reference densities. Finally the Pearson correlation (r) between R and D, is an index of a model's sensitivity or, in other words, a measure of how R is able to detect a trend. Finally we cross-validated the model. We divided our data set into years and subareas, then we recursively calibrated the model excluding one stratum, which was then used to test the performance of the model itself. In such a way, we tested the robustness of the method with respect to different environmental and temporal conditions. Note that R estimates at subarea and at the whole area scales were computed by averaging field estimates.
Results
The model for the whole data set
We do not show here the outcome for each of the parameters of the regressors. The plots of some of the relations that were minimised to compute regressors in order to provide an intuitive look of the fits are shown in Figure 1. In Figure 1A, we plotted apparent density U as a function of field area. Clearly, this relationship is not linear. In this example we have fitted a hyperbole to the data. Note how U decreases as the field area increases. In Figure 1B, we plotted density D against O. Again, the relationship does not look linear. In the figure we minimised a Holling type-2 function. The point is that density is not usually a linear function of the number of observed animals.
Two regressors were included in the selected model (Table 2) which is highly significant (F3,60 = 18.0, P < 0.0001, R2 = 0.47).
How the estimator (R) was good at fitting the reference density D at three different spatial scales is shown in Table 3. The estimates were unbiased and their precision was good at all scales. At large and medium scale the error was well below 20%. As expected the error decreased at an increasing scale. No significant bias was evident at any scale. The correlation indicated that model sensitivity was quite good for small and medium scales, suggesting efficient trend detection. The apparent lack of correlation for the large scale was possibly an effect of the reduced data set characterising this spatial scale.
Cross validation
The models that we used for cross validation per year and their results are reported in Tables 4 and 5, respectively. The results are qualitatively very similar to the results previously described for the complete model (see Table 3) suggesting that one can obtain reliable density estimates for one year using a model developed for different years. The bias is always non-significant at all scales and remarkably low at the small scale. In particular, the correlations relative to small and medium scales appear to be high, whereas the performance is lower at the large spatial scale.
Finally, we investigated whether one can use a model developed in one subarea to estimate density in another subarea (Tables 6 and 7), i.e. how much our model can be exported. These results (when compared to the previous cross validation per year) exhibit a lower precision at the medium scale, which remains below 20%, while the precision at small scale is even better than in the case of cross validation per time (see Table 5). Bias is always not significant. Correlations are good at both scales.
Looking in particular at the small scale, the correlation between the true density and our estimates is always significant (P < 0.01; see Tables 3, 5 and 7), while precision is always around 25%.
Our experience showed that the use of both AIC and VIF for model selection sensibly improved the cross validations between subareas with respect to the use of AIC only.
The cross validation by year (see Table 5) is comparable with the whole data set model (see Table 3) and the difference is partly due to the reduced data set used to calibrate the model. Overall, the results shown in Tables 3 and 5 suggest that the method is very reliable to estimate local density even in years following the model's calibration or in different subareas.
Discussion
The basic message in our work is that it is possible to use calibration models with simple and cheap data as inputs (in our case field area and number of counted roe deer) to estimate density at different scales, provided that reliable density data are available, e.g. obtained by use of mark-resighting. Our results are also encouraging if one wants to predict density using data collected in different years within the same area. Errors are usually < 20%, which is acceptable in most of the applications.
Our modelling approach is purely pragmatic in the sense that the functions used (the ρs) have no specific mechanistic interpretation. In this paper, we used a sort of ‘blind’ calibration of our density estimates. Calibration is usually defined as the reverse process of linear regression, by deducing an unknown value of the dependent variable (density) given known independent variables. However, given the complexity of factors linking sightability, habitat structure and density, we do not expect their relationship to be linear. In fact, even in sightability models non-linear (exponential) functions are used to predict animal abundance. The use of non-linear relationships is in fact commonplace when modelling animal detectability (White & Schenk 2001). A non-linear relationship between population indices and population abundance was also found by Morellet et al. (2007) for the roe deer population of Dourdan, France. In distance sampling non-linear detectability functions are used to account for different detection patterns (Buckland et al. 2001). Thus, our choice of ρs is arbitrary (i.e. not derived by mechanistic considerations), but empirically useful and justified from the literature. The method we propose has the potential to be used for validating population indices such as the kilometric index of abundance (KIA), which are often used as management tools for roe deer populations throughout Europe (Vincent et al. 1991). Recently, the use of KIA was extended to other species of ungulates as well; Maillard et al. (2001) supported its use for detecting population trends in an African ecosystem, and Acevedo et al. (2008) showed that the KIA was well correlated with line-transect density estimates of red deer Cervus elaphus populations in Spain.
We argue that our method was very successful at estimating local density of roe deer. We expected to obtain the best results for the estimates at larger scales than at a smaller scale. While this turned out to be true for the precision of our estimates, which was very good at intermediate and large scales, correlation and accuracy were worse at the larger scale. We assume that this unexpected result was caused by the reduced data set available when the spatial scale increased, but density differences between fields and woods may also have played a role.
Our results confirm the usefulness of our method both for research and management purposes. For research, our method can provide reliable local density estimates where a radio-marked set of roe deer is available. It allows insights in population dynamics and regulation, mainly by providing the researcher with information about the spatial structure of density.
Note also that the information about density in an area can be extended to zones in which animals were not captured, provided that the model is used only for zones similar to those in which the model was calibrated. This result was achieved by testing for collinearity using variance inflation factors (VIF). In a preliminary phase of our work, we saw that estimators obtained only with AIC were biased when used to estimate density between subareas. Interestingly, a single regressor was the best solution only in one case (see Table 6). In all other cases a combination of regressors performed better than any single regressor. This supports our choice of using a multiple regression framework.
In management, the use of simple counts may determine the adoptions of harvesting quotas which are inappropriate to fulfill the aims of the management. Integrating field counts with a zone-specific method calibrated in place, game managers can smooth steep interannual variations of counts, thereby gaining trust among hunters and allowing more stable game returns.
The application of our method relies on appropriate calibrations. Wildlife managers should invest resources (at least for some years) to radio-mark animals. Alternatively, it could be interesting to investigate the use of surface density modelling based on distance sampling survey (Buckland et al. 2004, Focardi et al. 2006) as a basis for the calibration model. Distance sampling surveys could be performed for some years at relatively low costs (Franzetti & Focardi 2006). Although we have not yet tried to develop this approach, we believe that the outcomes could be very fruitful.
Finally, our results will allow us to test the relevance of local density on roe deer population dynamics in the same population that was described in Focardi et al. (2002a,b). We will investigate the scale at which density is relevant to shape life history traits.
Although our method is quite specific for our needs, we believe that the general framework can be adapted to other species and environments.
Acknowledgments
Our field study was supported by a research grants from Regione Emilia Romagna, Provincia di Forlì-Cesena and Federazione Italiana della Caccia. We are especially indebted to S. Diverio (at the Università di Perugia), V. Guberti and coworkers (at ISPRA) for veterinary assistance. M. Amadei, B. Amadesi, F. Baldinelli, L. Carnevali, G. Colombi, P. Genovesi, V. Gervasi, A. Lenuzza, A. Monaco, M. Panzacchi, R. Petrucco, A. Piazzi, G. Roccasalvo, M. Sacchi, S. Santoro, S. Savini, A. Scappi and M. Scremin and many volunteers helped us during the capturing and tracking of roe deer. We thank the AFV Consorzio Iniziative Tredoziesi for logistic support. Comments by M. Hewison improved a first draft of the manuscript. We are grateful to the referees for their suggestions which substantially improved our paper.
References
Figure 1.
Examples of non-linear relations that were minimised to estimate the ρs; in A) the use of fields has an hyperbolic relation with the field area and in B) the relationship between D and O is Holling-type 2. The parameters of the equations were estimated from the whole data set. U is the ratio between apparent and true densities.
Table 1.
List of regressors used in our density estimation model. The functions in the right column (Regressor functions) were obtained by expliciting density D into the equation reported in the left column (Fitted functions) and replacing it with ρ. The abbreviations are: A = field area, O = counted deer, E = O/A, U = E/D, a, b, c, d and k are parameters computed by least squares.
Table 2.
Coefficients of the regressors (ρ1 and ρ7 ) entered into the selected model for the whole data set, standard errors, Student's t-tests for difference from 0 and the variance inflation factors (VIF).
Table 3.
Precision and bias of the estimator calibrated upon the whole data set, and the correlation between the estimator, R, and reference density, D, at the three different spatial scales small (≈ 10 ha), medium (≈ 200 ha) and large (≈ 400 ha). The median error (in %) was computed using the absolute value of the difference between the estimates derived from the model and the reference density. The percentile error was computed considering the percentage deviation from the reference density. The Student's t-tests for the presence of a bias were performed with N-1 degrees of freedom. Significant correlations are indicative of a good trend detection.
Table 4.
Cross-validation per year. The year reported in the left column refers to the data predicted by the procedure of cross-validation, which were hence excluded by the model's computation. For each year of study, we report the selected model, F-test and adjusted R2 for the selected model.
Table 5.
Cross-validation per year (see Table 4). Precision, bias and correlation between the estimator and reference density, at the three different spatial scales (small: ≈ 10 ha, medium: ≈ 200 ha and large: ≈ 400 ha). The median error (in %) was computed using the absolute value of the difference between the estimates derived from the model and the reference density, D. The percentile error was computed considering the percentage deviation from the reference density. The Student's t-tests for the presence of a bias were performed with N-1 degrees of freedom. Significant correlations are indicative of a good trend detection.
Table 6.
Cross-validation per subarea. The selected model and F-test and adjusted R2 for the cross-validation models are reported by subareas. The selection procedure was performed on our data set removing data relative to the subarea reported in the left column.
Table 7.
Cross-validation by subarea (see Table 6). Precision, bias and correlation between the estimator and reference density, at the two different spatial scales small (≈ 10 ha) and medium (≈ 200 ha; note that the large scale cannot be evaluated in this analysis). The median error (in %) was computed using the absolute value of the difference between the estimates derived from the model and the reference density, D. The percentile error was computed considering the percentage deviation from the reference density. The Student's t-tests for the presence of a bias were performed with N-1 degrees of freedom. Significant correlations are indicative of a good trend detection.
