Researchers and wildlife managers strive for low bias and high precision (i.e. high accuracy) when estimating animal population sizes. Distance sampling is currently one of the most widely used monitoring methods. However, it relies on strict sampling designs and modeling assumptions that can be difficult to meet in the field. Here, we use data from two sub-populations of non-migratory wild Svalbard reindeer Rangifer tarandus platyrhynchus inhabiting flat, open and isolated coastal tundra plains, to demonstrate some challenges related to the distance sampling methodology. To achieve this, we compared distance sampling line transect estimates with repeated total population counts and combined available software tools (R packages unmarked, Distance and dsm) to fulfill the analytical requirements of small study sites in which large areas are surveyed relative to the study area size. Based on low variation among repeated total counts (CV = 0.02 - 0.06) and the virtual absence of false negatives and positives of marked animals, the total counts could be used as reference population sizes. Distance sampling estimates were not statistically different from the total count estimates. Our relatively large sample size of 143 observations enabled precise distance sampling abundance estimates (CV = 0.16 - 0.26) compared with other studies in the wild. However, capturing the processes shaping population dynamics would likely require even higher sampling effort or other, more resource demanding monitoring tools, such as total counts or mark-recapture. In this type of ecosystem, distance sampling nevertheless represents a cost-effective tool suitable for ‘population state’ assessment and studies of large-scale spatial distribution patterns. Our study stresses the importance of choosing the appropriate analytical tools and estimating the accuracy of the monitoring methods that are used to achieve specific scientific, management or conservation goals.
Population size estimates with low bias and high precision (i.e. accuracy, Williams et al. 2002, p. 45) are important to understand spatiotemporal patterns of wildlife populations, thereby informing management and conservation decisions. While there are numerous challenges associated with obtaining unbiased and precise estimates of population size and demographic rates in space and time (Williams et al. 2002), such estimates are essential to understand population fluctuations and their causes (Abadi et al. 2010, Zipkin et al. 2014). Bias and imprecision originate from sources of errors that can occur at multiple levels in the measurement of population size and vital rates (Cressie et al. 2009, Lebreton and Gimenez 2013). These uncertainties are related to process variation (i.e. demographic and environmental stochasticity) and observational error (Clark and Bjørnstad 2004, Buckland et al. 2007, Sæther et al. 2007). Because observational error is not part of the process variation, but inherent to the sampling methodology used, it is important to identify its different sources (Ahrestani et al. 2013).
Worldwide, the most common method to estimate population abundance of wild animals is distance sampling (hereafter referred to as DS) (Buckland et al. 2004). Surveys are conducted along transect lines or at transect points where the detection probability is a function of the perpendicular distance from the line or radial distance from the point to the object of interest (reviewed by Buckland et al. 2015). DS relies on four key assumptions related to study design and statistical analysis of the data (Buckland et al. 2015): 1) animals are distributed independently of the transects; 2) objects on or close to the transects are always detected; 3) distances are measured without error; and 4) objects are detected at their original position. However, the degree to which the assumptions of DS are violated in wild populations is largely unknown (Morellet et al. 2007, 2011). Comparisons of bias and precision have been made between DS and other sampling methods. Some examples are comparisons with capture—mark—recapture, total counts (hereafter referred to as TC) and strip transects, in different environments (e.g. forest: Focardi et al. 2005, Wegge and Storaas 2009, Amos et al. 2014, grassland: Kruger et al. 2008, Amos et al. 2014 and steppe: Seddon et al. 2003, Bårdsen and Fox 2006). Even for similar types of DS, bias can vary greatly depending on how the field sampling is conducted (e.g. foot versus road line transects; Wegge and Storaas 2009, Marques et al. 2013). Overall, few studies have assessed the performance of DS by comparing estimates to populations of known size (but see Wegge and Storaas 2009, Porteus et al. 2011 for ungulates and Glass et al. 2015 for kangaroos).
Recent spatial modeling developments have been incorporated into the DISTANCE interface (Thomas et al. 2010), whereby a two-stage approach analyzes detection and density separately. This is particularly suitable for large-scale study areas. By contrast, in small study areas, the transect width will often cover most or all of the study area. Thus, there could be information about the spatial distribution of animals in the observed distances, as well as their detectability (Miller et al. 2013). Therefore, in small study areas, the detection and density functions should be estimated simultaneously (i.e. a one-stage approach, Miller et al. 2013; see Glass et al. 2015 for a counter-example of a two-stage approach where the search area of all transects combined included the entire 76 ha study area). Hence, in line transect studies, different analytical approaches may be required depending on the type of study area (Miller et al. 2013).
In this study, we use monitoring data from a high-Arctic wild subspecies of reindeer, the Svalbard reindeer Rangifer tarandus platyrhynchus, which inhabits tundra with sparse vegetation of low stature. The distinctive landscape characteristics of Svalbard are highly suitable to evaluate the precision and sources of error in two methods of estimating animal abundance; DS and TC. In this open tundra landscape, the detection of reindeer should in principle only vary with distance from the transect line as visibility is good. Additionally, many coastal sub-populations are isolated in small areas by glaciers, steep mountains and the sea and so are possible to census on foot by TC. Although never assessed quantitatively, it is assumed that TCs of Svalbard reindeer give precise and unbiased population size estimates, partly due to the open habitat and partly because of the restricted ranging and solitary behavior of the reindeer (Aanes et al. 2000, Kohler and Aanes 2004, Tyler et al. 2008, Hansen et al. 2011). Recent studies suggest that the rapidly changing climate in Polar regions (Larsen et al. 2014, Nordli et al. 2014) will strongly impact ungulate population dynamics (Rennert et al. 2009, Hansen et al. 2011, 2013). This underlines the need for robust estimates of population abundance. Here we take advantage of the characteristics of this simple, high-Arctic model system to compare reindeer abundance estimates made using DS and TC methodologies. In particular, we show that TC estimates are accurate and therefore usable as reference points (i.e. ‘known’ population sizes). We are then able to assess whether estimates based on DS analysis, using a combination of R packages (unmarked, Distance and dsm), are different from TCs. We further evaluate some sources of error and imprecision.
Figure 1.
Map of the Svalbard archipelago including the two study sites (Sarsoyra and Kaffiøyra) where total counts and distance sampling surveys of Svalbard reindeer were conducted. The vegetation map from Johansen et al. (2012) was plotted as background of the inset map (the darker the grey, the more of the area is vegetated).
Methods
Study system
The study was conducted in two sites, the Sarsøyra (40 km2) and Kaffiøyra (35 km2) peninsulas, close to the Ny-Ålesund scientific base (78°55/N, 11°55′E; Fig. 1), on Svalbard. The study sites lie in the northern Arctic tundra zone (Elvebakk 2005) which is characterized by graminoid and dwarf-shrub tundra. Reindeer summer habitat is generally confined to areas below 200 m elevation, excluding large moraines and glaciers. Natural barriers to reindeer movement include tidewater glaciers, steep mountains and more recently, year-round open water fjords. Both peninsulas are dominated by ‘pioneer vegetation’ (41% of the total coverage) and ‘established Dryas tundra’ (39%). The ‘pioneer vegetation’ is characterized by vegetation communities with low vascular plant diversity and mosses, highly affected by erosion and flooding. ‘Established Dryas tundra’ is dominated by Saxifraga oppositifolia and short-growing graminoids on coastal plains (Johansen et al. 2012). With the exception of some graminoids, plant height rarely reaches more than 5 cm above the ground.
The Svalbard reindeer is an endemic resident key herbivore that lives in a predator-free environment (but see Derocher 2000 for predation events by polar bear Ursus maritimus). This sub-species is tame compared to most other wild Rangifer, and individuals in our study sites can typically be approached by humans in summer to within 100 m before their behavior is affected (Hansen and Aanes 2015). Their space use during summer is relatively limited (Tyler and Øritsland 1989), although dispersal events can occur during winter (Hansen et al. 2010). Reindeer were absent from our study sites for a century before the sub-populations reestablished after a major dispersal event from Brogger Peninsula to Sarsoyra in 1994 (Kohler and Aanes 2004), with onward dispersion to Kaffiøyra in 1996 (Fig. 1). Fluctuations in the reindeer population sizes are large and typically associated with ‘rain-on-snow’ events, which cause high mortality and reduced population growth rates (Kohler and Aanes 2004, Hansen et al. 2011).
Reindeer data collection
Total counts and marked individuals
We conducted repeated TCs of Svalbard reindeer in 2009 and 2013 (July to August) in the two study sites (Table 1). The natural barriers, limited ranging behavior (Tyler and Øritsland 1989) and negligible summer mortality rates (Reimers 1983) mean that the assumption of constant population size within the area and field season (i.e. closed population) is likely met. This also means that the same population was counted during each repeated TC. Four predefined routes from south to north, always less than 1 km apart, allowed each of the study sites to be covered in one day by four observers walking simultaneously. Observers were not strictly confined to stay on their route, but expected to deviate from the line to utilize terrain features (e.g. small mounds) to get the best possible overview, and keep visual contact with other observers (e.g. to avoid double counts). Observers scanned the area along their route with binoculars (10 × 42 mm, allowing reindeer detection up to 2–3 km distance) and communicated by VHF-radio to reduce the potential for double counts. Routes were switched between observers to reduce bias related to observer heterogeneity (Field et al. 2005). All reindeer positions (i.e. of single individuals or groups) were marked on a topographic map (1:50 000). Repeated TCs were always separated by a minimum of four days. All TCs had similar weather conditions with good visibility and little wind.
Table 1.
Overview and summary statistics of the repeated total counts (TCs) of reindeer In the two study sites, Sarsøyra and Kaffiøyra, Svalbard, Norway. TC values = reindeer abundance from each TC, TC estimates = mean reindeer abundance and 95% confidence Intervals, SE = standard error and CV = coefficient of variation.
In addition, we used data from TC and marked animals in previous years (1999 and 2000) to evaluate 1) the closure assumption based on data from regularly (every 2–3 days) tracked VHF-equipped females (19 in 1999 and 23 in 2000, Hansen et al. 2009); 2) the probability of missing animals in the TC based on the number of VHF-collared females known to be present in the study site immediately before the TC (27 individuals in 2000, which includes the 23 females followed every second-third day); and 3) the probability of double counting based on re-sightings of the VHF-marked animals (27 individuals) as well as other marked animals (26 individuals) during the TC in 2000.
Distance sampling
One single observer conducted line transect DS twice in each of the two study sites in 2013 (12 and 19 July in Sarsøyra and 26 and 27 July in Kaffiøyra). We chose one random latitude for each of the two DS surveys and placed additional parallel transect lines systematically 3 km away either north or south from this latitude so as to avoid overlapping reindeer observations and violation of the assumption of independence (Hammond et al. 2014, Buckland et al. 2015). Lines were orientated east/west from the sea-shore to the mountain foothills. We chose this systematic orientation to reduce any potential bias from parallel animal density gradients along the line (e.g. due to plant phenology gradients; Marques et al. 2013, Barabesi and Fattorini 2013). In total, 11 transect lines were walked in the study sites (Sarsøyra, total length 19 029 m, three transect lines in first survey and two in second survey; Kaffiøyra, 14 937 m, four transect lines first survey and two in second survey). Because of the random placements of lines in the small study locations and because two transect lines (one in each study site) were not completed due to bad weather, the surveys had different total numbers of transects.
The line transects were walked by the observer at a constant speed (2–3 km h-1) without stops, except during measurements. A handheld GPS was used to keep the line direction. Reindeer were detected on both sides of the transect line with the naked eye. When a reindeer or reindeer group was spotted, the observer looked only in its direction until measurements were taken. Each observation is hereafter referred to as a cluster, regardless of whether it is an individual or a group to meet the assumption of independent detection between observations (Buckland et al. 2001, Guillera-Arroita et al. 2012). The geographic position of the observer was recorded. We used laser binocular and compass to measure the distance and angle from the observer to the reindeer. For practical reasons using the laser, measurements were taken to the largest reindeer (e.g. a mother rather than her calf), or the left-most individual of a group of adults. We acknowledge this as sub-optimal to measuring the center animal but considered the associated potential positioning bias as negligible (reindeer belonging to the same group are close to each other and mean group size is <2) and evened out by following the same procedure on both sides of the transect line. If a reindeer individual or cluster was beyond the distance that the laser could measure with confidence (sometimes down to ∼500 m), positions were marked on a topographic map (1:50 000) and the exact perpendicular distance was subsequently read from an electronic map (< www.toposvalbard.no>). Even in this relatively flat terrain, we consider based on qualitative inference that the mapping method is correct due to the fine map resolution (e.g. every little creek is shown) and the supporting use of a GPS. In order to reduce error in our estimation of mean cluster size, the observer counted the number of individuals in a cluster with binoculars. If a cluster that had not been observed initially became apparent while measuring the distance or cluster size (i.e. because of using binoculars), it was not included in the observation data. We only conducted DS surveys when conditions were adequate in terms of good visibility and little wind.
Data analysis
Total counts
To estimate reindeer abundance from TC data and its uncertainties we considered two types of errors related to the observers; 1) an animal could be counted twice with a probability p or 2) an animal could be undetected with a probability q. We expected the proportion of these two error types to be constant across years and similar in the two study sites, which share the same flat and open tundra landscape and therefore detectability of reindeer. Let X denote the total number of reindeer counted twice, Y the number of undetected animals, and N and n the estimated and true population size, respectively. Animal abundance from TCs can be expressed as N = n + X-Y. Given that p and q are small, N is approximately unbiased for the true population size n. We have several N (i.e. repeated TCs) of the same true abundance n at the different sites and years (Table 1). Note also that the number of individuals counted twice, never or once, i.e. X,Y and n-X-Y respectively, is multinomially distributed with parameters p,q, 1-p-q, if individuals are counted independently. From variance and covariance formulas for the multinomial distribution, a straightforward calculation then shows that the variance of a single count, Var (N) = Var (n + X-Y), is proportional to the mean E(N) = E(n + X–Y) with a scale parameter
corresponding to the dispersion parameter of a quasi-Poisson likelihood model (McCullagh and Neider 1989). Estimates, standard errors and profile likelihood 95%-confidence intervals for the true abundances n at each study site and year were obtained by fitting a generalized linear model from the quasi-Poisson family, with no intercept and with a four level fixed effect factor representing the two study sites and the two summers (i.e. 2009 and 2013). We obtained the coefficient of variation (CV) for n by dividing the standard error by its mean abundance estimate.
Distance sampling
Prior to estimating reindeer abundance from the DS line transect data, we divided the 11 line transects into 33 segments, i.e. three equal segments per transect (Sarsøyra: 15 segments 728–1544 m, Kaffiøyra: 18 segments 107-1047 m). This allowed us to fit density models that included covariates measured at the segment level (Royle et al. 2004, Miller et al. 2013). Transect lengths were relatively short compared with their half-widths due to the small scale of our study sites and the long detection distances of reindeer. Note that segment lengths could not be twice the half-transect width, because this would give only one to two segments per transect and reduce our ability to detect effects of habitat heterogeneity. We calculated the proportion of vegetated area within each segment area (segment length × truncation distance, ranging from 0.43 to 2.89 km2) from a digital vegetation map (Johansen et al. 2012). This was the ratio of pixels (spatial resolution of 30 X 30 m) classified as vegetated (corresponding to classes 8 to 18 in Johansen et al. 2012) to the total number of pixels. We also right-truncated the DS data, as suggested by Buckland et al. (2001), by removing 5% of the reindeer clusters that were most distant from the line. The furthest observation after truncation was 953 m from the line and twice this distance defined the width of the surveyed area along the transect lines. Following Buckland et al. (2001, p. 109), we calculated along the segments of total length L, the encounter rate w/L, the encounter rate variance (n/L) (also corresponding to Fewster et al. 2009 ‘R3’ estimate) and evaluated the homogeneity of clusters position using the ratio between the expected number of objects detected E (n) and the sampling variance var (n). A ratio close to 1 gives no evidence against a Poisson distribution. We used Pearsons correlation tests to investigate the correlation between the segment-based encounter rate or vegetation cover of adjacent segments on the same transect (22 total possibilities).
We estimated reindeer abundance by combining statistical tools available in R ver. 3.2.2 (< www.r-project.org>) and the packages unmarked ver. 0.10-2 (Fiske and Chandler 2011), Distance ver. 0.9.4 (Miller 2014) and dsm ver. 2.2.4 (Miller et al. 2013). Model selection was done in unmarked because it uses a one-stage model selection procedure (i.e. full likelihood approach), which is required for spatially confined study sites (Miller et al. 2013), as illustrated in earlier studies (Royle et al. 2004, Royle and Dorazio 2008, Johnson et al. 2010, Miller et al. 2013). In unmarked, the data have to be pooled into distance intervals. These were set to 1 m, simulating a continuous fit, because the cluster positions were precisely measured and we wanted to fit similar models in Distance/dsm. A hierarchical DS model implemented in the function ‘distsamp’ uses the multinomial-Poisson mixture (Royle et al. 2004, Fiske and Chandler 2011) to compute detection and density parameters simultaneously. The covariates Zi of transect i were related to the detection parameter σi and mean density parameter λi¡ with the log link function and α and β as their respective parameters estimates (Royle et al. 2004, Fiske and Chandler 2011, Siliert et al. 2012) as follows:
Different model combinations were analyzed using a half normal, hazard rate, exponential or uniform key for the detection function. No adjustment terms to the key functions could be implemented in unmarked. For each key detection function we included either the study site, the vegetation cover or none of these as covariates, while for the density estimation we included either the study site, the vegetation cover, both (i.e. additive effect) or no covariate. To avoid over-parameterization, no more than six parameters were allowed per model. We retained the best model ranked by Akaike information criteron (AIC, Burnham and Anderson 2002) from the unmarked model selection and implemented it in Distance/dsm using the same distribution, key functions and generalized linear models as presented in Eq. 2–3 (Marques et al. 2007, Miller et al. 2013). Only models with ΔAIC < 2.00 were presented in Table 2. A Freeman–Tukey goodness-of-fit statistics with 500 bootstrap iterations (h2 , Brooks et al. 2000, Siliert et al. 2012) evaluated the model fit to the data. We verified whether group size was a potential covariate influencing detection at the observation level in Distance, since it was not possible in unmarked, using the top models from Table 2. A linear regression investigated whether group size increased with increasing vegetation cover (the vegetation cover around each DS cluster was extracted with a 200 m buffer; see details below).
We compared three different approaches to estimate animal abundances from the DS analyses. Firstly, we estimated cluster density using the top ranked model (termed ‘estimate 1’). Secondly, we estimated density using model averaging, including all models with ΔAIC < 2.00 (‘estimate 2’). In both estimate 1 and 2, computed in unmarked, we used the function ‘predict’ (Fiske and Chandler 2011) to extract mean cluster density and its standard error for the proportion of vegetation found across the study sites. Then, we transformed estimated cluster densities to cluster abundances by multiplying by the respective area of each study site. Thirdly, we estimated cluster density using spatial density modeling methodology (Miller et al. 2013) in dsm. The top ranked model (Table 2) was used in the function dsm with segment length as an offset and the family set as quasi-Poisson with a logarithmic link function. We verified that the quasi-Poisson scale parameter was close to 1, indicating an absence of over-dispersion. The mean cluster abundance and standard error were extracted from the variance propagation function ‘dsm.var.prop’ (‘estimate 3’). This standard error obtained from dsm do not account for the variance related to the detection function. Following La Morgia et al. (2014) we thus summed the coefficient of variation related to the detection function (Distance) with the coefficient of variation related to the spatial density modeling (dsm) using the delta method. From these calculations, we obtained the mean cluster standard error. Density from equation 3, using the function ‘predict’, was estimated at the pixel level (spatial resolution of 30 X 30 m) and afterward summed over the study site to obtain the site-specific cluster abundance estimate. Detection and density covariates (i.e. vegetation proportion) were thus required for each pixel. Because vegetation information from the vegetation map is binary (vegetated or not) at the cell level, we assigned each pixel a value calculated from a neighborhood consisting of a circular buffer zone (radius of 200 m, corresponding to an average of 63 pixels) around each mid-point of the pixel.
Table 2.
Parameter estimates of the seven top ranked models (ΔAIC<2.00) for estimating Svalbard reindeer abundance using distance sampling analyses computed with the R package unmarked. These models used a continuous fit of the observed distances (1 m distance Intervals), and the hazard rate (hz), half normal (hn) or exponential (exp) detection key function. Proportion of vegetation In transect line segments (veg) and/or study site (Sarsøyra or Kaffiøyra) were Included as detection covarlates (σ) and/or density covarlates (λ). Estimated abundance = mean abundance estimate with 95% confidence Intervals (corresponding to estimate 1 in Methods); CV = coefficient of variation.
Estimated cluster abundances 
and their corresponding standard errors 
from estimates 1, 2 and 3 were combined with mean expected clusters sizes 
and their standard errors 
for each site i, into estimates of the abundance of individuals, where Assuming independence between 
and 
and using an exact formula for the variance of products (Goodman 1960) rather than relying on the approximate delta method, the standard errors of abundance of individuals are
Upper and lower 95% confidence intervals of individual abundance estimates were obtained using a normal approximation.
Total counts and distance sampling comparison
We calculated the difference between TC and DS reindeer abundance estimates with their 95% confidence intervals obtained using a normal approximation of mean 
and standard error 
. Thereafter for each site i in 2013, the probability that TC and DS estimates were significantly different from each other (in a two-sided test) used the absolute value of the normally distributed Z-statistic:
We compared cluster density with vegetation cover intervals (every 10% from 0% to 100%) from both DS modeling and TC mapped cluster positions. In DS modeling we used the relationship from the top ranked model and predicted the cluster density using the function ‘predict’ from unmarked (corresponding to estimate 1). For each vegetation cover interval, we also summed the pixel-wise (30 × 30 m resolution) density estimated from spatial modeling (corresponding to estimate 3). For the TC, the vegetation cover around each cluster was extracted with a 200 m buffer and number of clusters were summed in intervals (every 10% from 0% to 100%).
Results
Total counts
The TC precision was high in both study sites and years (CV ranging between 0.02 and 0.06, Table 1). The high re-sighting rate of the VHF-marked female reindeer (100% in 2000) and the low rate of double counts (1.9%, i.e. one out of 53 reindeer) suggested a low bias of the TCs. Virtually all the VHF-collared female reindeer that were closely tracked stayed within their respective study site throughout the summer (97.5%, only 1 out of 42 changed study site; 19 in 1999 and 23 in 2000). This largely confirmed our assumption of closed populations. We expected the variance of the animal abundance to increase proportionally to the mean (Eq. 1) and estimated the dispersion parameter Φ to be 0.45 [0.10:0.92] (point estimate and 95% confidence interval derived via the approximate χ2-distribution of 
, Venables and Ripley 2002, p. 210). This under-dispersion relative to a Poisson variance is not surprising since the variability of the total counts originate from small probabilities of individuals being undetected or double counted.
Distance sampling
In total, we observed 143 reindeer clusters (88 in Sarsøyra and 55 in Kaffiøyra) during the DS surveys. There was a slight increase in the number of clusters detected between approximately 90 m and 170 m from the transect lines (Fig. 2). The mean cluster size was 1.68 ± 0.10 (mean ± SE) and 1.58 ± 0.11 in Sarsøyra and Kaffiøyra, respectively. The largest group had five animals and the correlation between group size and detected distances, if considered independent, was at the critical level to be statistically significant (cor = 1.16, t = 1.97, df = 141, p = 0.05). Nonetheless, we found no evidence in our data that cluster size influenced detection probabilities (ΔAIC = 4.17). Similarly, we found no evidence that vegetation cover influenced reindeer cluster size (linear regression: intercept = 1.47 ± 0.26, slope = 0.22 ± 0.32, p = 0.50). The encounter rate had a low segment-based variance, 4.62 ± 0.59 (n/ L + SE[n/L], Buckland et al. 2001, p. 109) and 3.68 ± 0.61 observations per km in Sarsoyra and Kaffiøyra, respectively. Adjacent segments on a transect were not found to be correlated in either their segment-based encounter rate (cor = —0.16, t = -.72, df = 20, p = 0.48) or their vegetation proportion (cor = 0.20, t = 0.91, df = 20, p = 0.37). The low ratios between the expected number of objects detected and sampling variance (Sarsøyra = 1.42 and Kaffiøyra = 1.53) indicated no evidence against a homogenous Poisson distribution of reindeer clusters along the segments, i.e. they could occur at any position. Using a higher truncation percentage did not decrease this ratio. Accordingly, no over-dispersion was found in dsm (quasi-Poisson scale parameter of 1.04).
Figure 2.
Detection probability function based on line transect distance sampling of Svalbard reindeer. The best model was fitted at a continuous scale for observed distances and included a hazard rate key detection function with no covariate. Observations of reindeer clusters are illustrated by dots along the curve. The rescaled histogram was plotted at 30 m distance intervals for ease of the visual illustration.
All models with ΔAIC < 2.00 included vegetation cover (proportion of vegetated area in each transect segment, ranging from 0.13 to 0.86) as a covariate positively influencing the cluster density function (Table 2). None of these models showed a significant lack of fit to the data (Freeman— Tukey goodness-of-fit statistics; 260.26 ± 19.77 < h2< 263.87 ± 22.26 where 0.40 < p < 0.53, the top-ranked model having the best fit). The top-ranked model used the hazard rate key detection function, with no covariate or factor influencing detection probability. As expected, transferring the best model from unmarked (with 1 m distance intervals) to Distance gave a similar detection probability (α0 from equation 2 is 5.77 ± 0.27 in unmarked and 5.77 ± 0.36 in Distance) and density parameter estimate (β1 from equation 2 is 2.02 ± 0.53 in unmarked and 2.02 ± 0.54 in Distance). The vegetation cover along the transect lines (69.6% in Sarsøyra and 54.1% in Kaffløyra) was comparable to the total vegetation cover of the study sites (66.1% in Sarsøyra and 50.5% in Kaffløyra). The precisions of the abundance estimates from all three DS estimates were ranging between CV = 0.16–0.26 (Table 3).
Total counts and distance sampling comparison
The mean reindeer abundance estimated from the TCs in both study sites were always within the DS 95% confidence intervals from the top ranked models (ΔAIC<2.00; Table 2). Accordingly, from the Z-statistics (Eq. 5, Table 3), none of the three DS abundance estimates differed significantly from the TC estimates in either of the sites. However, the precision of DS estimates was considerably lower than from the TC estimates, with the upper limit of the confidence interval close to twice the mean of TCs abundance estimates, Table 1, 3).
The importance of vegetation cover for spatial cluster density modeling was also supported by reindeer cluster positions from repeated TCs (Fig. 3). Nonetheless, based on the TC cluster positions this relationship appeared non-linear, with a sharp increase in reindeer abundance when more than about half of the ground was vegetated. Overall, the estimated density in the DS modeling tracked this apparent non-linearity well, yet slightly over- and underestimated reindeer abundance at low and high vegetation cover, respectively (Fig. 3).
Table 3.
Svalbard reindeer abundance In the two study sites (Sarsøyra and Kaffløyra) estimated using distance sampling. The estimated reindeer abundances are shown according to the three estimation methods (see Methods for details). Estimate 1 = the estimated abundance using the top ranked model, and corresponding proportional cover of vegetation across the study site (unmarked). Estimate 2 = the estimated abundance using model averaging (models with ΔAIC<2.00) and corresponding to the proportion of vegetation across the study site (unmarked). Estimate 3 = the estimated abundance using the top ranked model and corresponding sum of densities projected for each pixel across the study site (30×30 m resolution, pixels have a vegetation proportion value, see Methods, Distance/dsm). CV = coefficient of variation, Z = the Z-statistic (see Methods), p = p-value from the Z-statistic. Difference = the difference between total counts and the respective distance sampling abundance estimates. Numbers In brackets are 95% confidence Intervals.
Figure 3.
Histogram of the frequency of Svalbard reindeer clusters observed during total counts (left y-axis), pooled across study sites, plotted in relation to the proportion of vegetation covered area within a 200 m buffer zone around each reindeer cluster. Reindeer cluster observations were mapped during four repeated total counts (TC) in Sarsøyra and one total count in Kaffiøyra in 2013 (see Table 1). Estimated reindeer cluster density (right y-axis, solid curve) and associated standard error (dashed curves) in relation to vegetation cover, based on the selected model obtained from the R package unmarked. Red lines on the x-axis show the vegetation cover measured of each segment.
Discussion
In the present study, we used total count (TC) and distance sampling (DS) line transect data from two isolated sub-populations of high-arctic Svalbard reindeer to compare methodologies of abundance estimation. We evaluated some of the potential sources of error and imprecision in both methodologies. We obtained unbiased and precise TC abundance estimates (Table 1) from the two reindeer sub-populations that inhabit a flat and open tundra landscape. These provided ‘known’ population sizes and were used as a reference for DS abundance estimates. We found no statistically significant differences between the DS and the TC abundance estimates, but DS estimates were considerably less precise (Table 3). The vegetation cover proved an important covariate for estimating reindeer density spatially (Fig. 3, 4).
Biased abundance estimates can result from various reasons, including violation of the major assumptions of the DS method, selection of statistical models for density estimation, or software limitations. In this study, the reindeer DS abundance estimates were not statistically different from the respective TC estimates. However, although not significant, our DS estimates tended to be consistently larger than TC estimates at both study sites. If this tendency is reflecting a true difference, the lack of a significant difference may be due to the large variance of DS estimates (Eq. 5). Porteus et al. (2011) showed that attraction of sheep toward the observer inflated detected clusters around 100 m, resulting in underestimation of the detection probability and positively biased density estimation. Similarly, our histogram of detected distances close to the line indicated a hump between 90m and 170 m (Fig. 2), and the hump's effect on the hazard rate function could have led to an underestimated detection probability which further could overestimate reindeer abundance slightly (Fewster et al. 2008, Marques et al. 2010). Nonetheless, reindeer movement towards the line is not likely because of lack of Svalbard reindeer reaction towards human presence (Reimers et al. 2011, Hansen and Aanes 2015). We suggest this hump could have occurred by chance (see also a study on robin in Buckland et al. 2015, p. 71). Another explanation for the tendency of the DS estimates to be larger than the TCs could come from bias in the TCs themselves, but we consider this very unlikely, given the information from marked animals on false positives and negatives.
Figure 4.
Map of the estimated reindeer abundance for each 30 × 30 m pixel (900 m2, estimate 3; see Methods) in the two study sites, Sarsøyra (40 km2, upper figure) and Kaffiøyra (35 km2, lower figure), obtained from the R package Distanceldsm. Red horizontal lines show the distance sampling transect lines, the dotted lines show the 5% data truncation distance at 953 m and the black dots the positions of reindeer clusters from the first survey in both study sites.
Other studies comparing TC and DS line transect performance have been conducted in ecosystems with more environmental complexity, for instance at sea (Williams and Thomas 2009 killer whales Orcinus orca), in forests (White et al. 1989 mule deer Odocoileus hemionus; Wegge and Storaas 2009 Nepal's ungulates) or in human disturbed landscapes with high animal density (>100 animals kmsup>-2 in Porteus et al. 2011 for sheep and Glass et al. 2015 for kangaroo Macropus giganteus). The precision of our three DS abundance estimates was high (CV= 16-26%, Table 3) relative to the precision levels in other studies from the wild. However, the upper confidence interval of around 400 reindeer in Sarsoyra represents a density of ∼10 reindeer km2, which is much higher than reported densities from more productive areas in Svalbard (∼6 reindeer km2 in Adventdalen, Tyler et al. 2008). A larger sample size (i.e. more transect lines and a larger number of observations) would likely increase model quality by reducing the confidence intervals and thus excluded such biologically unreasonable abundances (La Morgia et al. 2015).
Habitat structure and heterogeneity are also important potential sources of uncertainty (Pedersen et al. 2012, Sillett et al. 2012), yet the habitat structure was well captured overall by our DS study design, with only a negligible 3.5% difference in vegetation cover within the transects areas versus the total study site areas. Thereafter, a design-based method extrapolating the density from the transects areas to the whole study area would give comparable results to estimate 1 (predicting density to the vegetation cover of the total study site). While Miller et al. (2013) recommend a larger grid cell size to match the segment scale, our smaller grid cell size (i.e. 200 m circular buffer) captures changes in reindeer density in response to covariate values. Both estimated DS cluster density (from unmarked, using the segment scale) and the cluster number from TC (using the 200 m circular buffer scale) clearly responded to the changes in covariate values by a sharp increase at about half of the ground being vegetated (Fig. 3). The precision of the three different reindeer density estimates (i.e. using unmarked or Distance) was fairly similar, while estimating reindeer density at a small spatial scale (estimate 3, spatial resolution of 30 × 30 m in Distance) tended to give the largest and more imprecise abundance estimates (Table 3). One reason for this could be that although overall similar, the estimated vegetation cover effect on density deviated slightly from the actual relationship (Fig. 3) at low and high vegetation cover, as visualized by the animal positions from TC. Greater flexibility in the model building within a single analytical tool, could select a more precise model. A one-step modeling approach tool that enabled spatial density modeling (available in e.g. R package DSpat, Johnson et al. 2010, but no transect overlap is possible) with adjustment parameters for the detection key functions and covariates at the individual observation level (only available in Distance) would be ideal for small study areas.
Clearly, even with our relatively large sample sizes (recommended minimum of 60–80 observations, Buckland et al. 2001) and the associated DS uncertainty level, we may not be able to relate our DS abundance estimates to environmental drivers of annual population size fluctuations. By contrast, using TC time-series from these sub-populations, drivers like ‘rain-on-snow’ events and density dependence, have previously determined population growth rates (Kohler and Aanes 2004, Hansen et al. 2011). Hence, a mechanistic understanding of the system relating the state variables’ responses to the environment requires a year to year unbiased and precise estimation of abundances (Yoccoz et al. 2001) or estimation of ecological indicators combining abundance indices, animal performance (i.e. body condition, reproductive and survival rate) and habitat quality (Morellet et al. 2007, Alb on et al. 2016). In practice, mark—recapture is resource demanding, and TCs are not possible across large study areas, such as the much larger coastal plains found elsewhere in Svalbard, due to high demands for resources and logistics, as well as the challenges associated with false positives and negatives. In such cases, DS represents a promising alternative monitoring method to estimate animal abundance. With wider applications than analyses of annual population dynamics per se, DS may provide information on long-term population trends (if repeated annually) and may be used as a state evaluation tool across wider spatial scales (Buckland et al. 2004, Harrison et al. 2016). This could relate to major changes in animal spatial distribution and abundance declines after extreme climate events or overharvest like the local extinction and following re-colonization of Svalbard reindeer after protection in 1925 (Reimers 1983). Indeed, it is likely that a relatively low-cost DS monitoring program, repeated at regular intervals, would have detected the dramatic population crash (from 360 to 80 animals) on Brøgger Peninsula during the extremely icy winter of 1993–1994 (Aanes et al. 2000, Kohler and Aanes 2004), as well as the associated re-establishment of the Sarsøyra and the Kaffiøyra sub-populations. Nevertheless, a population abundance monitoring method without an investigation of the accuracy and estimation of detection probabilities is never reliable enough to draw strong inferences about the monitored system (Yoccoz et al. 2001).
Management implications
Accurate population size estimates improve our ability to detect the effects of environmental drivers or anthropogenic stressors on population dynamics in space and time (Yoccoz et al. 2001, Clark and Bjørnstad 2004, Knape et al. 2013). Because biased or imprecise estimates may lead to erroneous conclusions (Buckland et al. 2007), we encourage researchers and managers to apply the most appropriate statistical modeling tools to investigate the uncertainties and sources of error in their wildlife monitoring programs, and then adjust their management aims and decisions accordingly. TC is a resource demanding and logistically difficult field method for large study areas, especially when the detection probability is low due to environmental characteristics such as rugged terrain. The alternative DS methodology is widely used worldwide, often by non-experts. Nonetheless, a target sample size and precision of around 80 observations and a CV of 0.15 (Buckland et al. 2001, Porteus et al. 2011), may be insufficient to understand and predict dynamic processes related to e.g. climate change or harvest. Glass et al. (2015) stated that “considerably more than the recommended 60–80 observations (Buckland et al. 2001)” should be aimed for to estimate abundances. If the resources limit larger sample sizes and a higher precision, we advise the use of DS methodology for wide spatial scale population state assessments. For a mechanistic understanding of the system, DS should preferably be combined with other measures like highly precise local abundance estimates or animal performance (Morellet et al. 2007). When it comes to field work costs, the monitoring of our study area is approximately twice as expensive using TC versus DS methods, mainly due to the difference in number of people involved. This may be an important argument when choosing methods for long-term population monitoring at large spatial scales (La Morgia et al. 2015).
When choosing DS methodology, animal observation sample size and relevant habitat covariate information should always be maximized to obtain unbiased detection curves and accurate habitat structure effects (Miller et al. 2013, Buckland et al. 2015). Underestimation of detection probability and overestimation of animal density, even for a onetime state assessment, can have fundamental implications in management and conservation of wildlife populations (Thompson et al. 1998). To prevent this, we recommend the further development of easily accessible tools, preferably within a single R package. Such developments should simultaneously allow for: 1) a one-step modeling approach for relatively small study areas (i.e. the area covered is large relative to the entire study area); 2) modeling with data on a continuous scale; 3) adjustment parameters for the detection key functions; 4) detection and density covariates at the individual observation level; and 5) spatial density modeling tools that permit transect lines to overlap.
Acknowledgements
Thanks to the Sverdrup Research Station staff for logistical field help. We also thank B.-J. Bårdsen, S. T. Killengreen, S. Hamel, M. L. Bender, J. M. Milner and L. Little who provided comments on an early draft of the paper and V. Grøtan and I. Herfindal who gave useful advice in the analysis of the data.
Funding — This project was funded by the Research Council of Norway (POLARPROG project 216051, the Arctic Field Grant project 227456, NORKLIMA project 178561, SFF-III project 223257), the Norwegian Polar Institute and Centre for Biodiversity Dynamics at the Norwegian University of Science and Technology.
