Estimating brown bear abundance and harvest rate on the southern Alaska Peninsula

Abundance estimation of hunted brown bear populations should occur on the same geographic scale as harvest data analyses for estimation of harvest rate. Estimated harvest rates are an important statistic for managing hunted bear populations. In Alaska, harvest data is collected over large geographic units, called Game Management Units (GMUs) and sub-GMUs. These sub GMUs often exceed 10,000 km2. In the spring of 2002, we conducted an aerial survey of GMU 9D (12,600 km2) and GMU 10 (4,070 km2) using distance sampling with mark-resight data. We used a mark-resight distance sampling method with a two-piece normal detection function to estimate brown bear abundance as 1,682.9 (SE = 174.29) and 316.9 (SE = 48.25) for GMU 9D and GMU 10, respectively. We used reported hunter harvest to estimate harvest rates of 4.35% (SE = 0.45%) and 3.06% (SE = 0.47%) for GMU 9D and GMU 10, respectively. Management objective for these units support sustained, high quality hunting opportunity which harvest data indicate are met with an annual harvest rate of approximately 5–6% or less.


Introduction
Brown bears (Ursus arctos) are a long-lived species of low reproductive rate that are managed conservatively in most coastal areas of Alaska, but often without the aid of estimates of population abundance or harvest rate. The Alaska Peninsula, which is Game Management Unit (GMU) 9 managed by the Alaska Department of Fish and Game (ADFG), produces large brown bears that are highly prized by hunters and professional bear guides. To manage this resource more actively for sustained yield and trophy characteristics, accurate and precise estimates of population abundance can be obtained to determine and potentially modify the harvest rates of these populations. ADFG monitors bear harvest on a subunit basis, of which there are 5 in GMU 9 (e.g., GMU 9D). Brown bears are managed for high quality hunting opportunity in GMU 9 and GMU 10 on Unimak Island. Trophy-related management objectives include harvesting a minimum number of large male bears � 8 year of age, males with skull size (length plus width) of � 71.1 cm (28 inches), and a low proportion of females (<40%) in the harvest. The ability to estimate harvest rates allows wildlife managers to determine if a1111111111 a1111111111 a1111111111 a1111111111 a1111111111 detection function [2]. A two-piece normal detection function, with or without covariates, has one apex at the distribution mode and can be used in MRDS models that assume point independence [2].
We implemented an aerial distance sampling survey in the spring of 2002 to obtain population estimates of brown bears in GMU 9D and GMU 10 (Unimak Island). Our goals for the survey were to obtain unbiased and precise population abundance and harvest rate estimates for each area. In addition, we wanted to examine if the data from these 2 areas could be combined by treating the 2 areas as separate strata [13] and if so, evaluate any gains in precision.

Study area
Our study area consisted of 2 areas for which bear abundance estimates were desired. The first was the mainland (non-island) portion of GMU 9D and the second was the portion of GMU 10 which consisted of Unimak Island (Figs 1 and 2). GMU 9D was a 12,599.7 km 2 area at the end of the Alaska Peninsula. Elevations ranged from 0 to 2518 m; diverse habitats included: cobble and sand beaches, tidal estuaries and lagoons, sand dunes, the Bristol Bay coastal plain of rolling tundra, mountainsides with alternating bare ground of sand, ash and cinder (referred to as cinder blows), alder (Alnus spp.), willow (Salix spp.) and grass dominated habitats, bare rocky mountain tops and slopes, and snow covered mountains and volcanoes at the highest elevations. Pavlof Volcano on the Pacific side was active and erupts frequently. Three small towns (1,100 people total) with minimal roads are within the GMU 9D study area. Unimak Island was the eastern most Aleutian Island, was 4,069.9 km 2 large, and the only one inhabited by brown bears. Elevations range from 0 to 2857 m; diverse habitats include: cobble and sand beaches, old lava flows, cinder blows, rolling tundra, mountainsides with alternating alder and grass dominated habitats, bare rocky mountain tops and slopes, and snow-covered volcanoes. Volcanoes are prominent features on the island; Shishaldin Volcano is frequently active. The island is roadless and had about 42 inhabitants in the fishing village of False Pass.

Methods
We conducted aerial distance sampling, out of Cold Bay, Alaska, in GMU 9D (degrees latitude 52.20, longitude 162.71), from 29 May to 11 June 2002 using 20 km long transects. We used 5-Piper super cub aircraft, a tandem aircraft with a pilot and backseat observer to fly daily surveys. To ensure crew safety, we split the study area up into 5 zones and assigned each plane their own zone to survey the transects located within their zone. If a transect straddled 2 zones, radio-communication between the 2 aircraft ensured it was safely surveyed by one of the 2 planes. Periodically, we shifted zone boundaries to allocate more resources to areas that needed additional resources to complete the survey. No permits were required for this project because surveys were conducted in the air and no animals were pursued or captured. Survey conditions ranged from good to poor; poor conditions included high winds, turbulence, or fog. We instructed crews not to survey transects under poor conditions, instead go to other transects with acceptable conditions and survey them. If conditions for all transects were unacceptable, the aircraft either returned to Cold Bay or landed on the beach to wait for improved conditions. Acceptable conditions included wind and turbulence conditions that allowed the safe and low level (100 m above the ground) flight operation of the aircraft, and sufficient lighting to see bears. The quality of the lighting was variable and included soft warm light, harsher direct sunlight, and flat light. Lighting is hard to quantify; however, crews were instructed to discontinue surveying if they felt lighting conditions were making it too difficult to see bears.
Survey methods to obtain the distance data and mark-resight data and meet MRDS model assumptions was described by Becker and Quang [12]. In the mountainous terrain, contour transects [14] were used. To sample the narrow irregular portions of GMU 9D, straight-line transects were modified to contain a hinge. For each randomly selected midpoint 2 angles were randomly selected to delineate the legs of the transect, if the angle created between the 2 legs was less than 30 degrees, the angle selection process was repeated. If the transect did not fit within the study area, the angle selection process was repeated. Hinged transects in the unexpected mountainous terrain on the western end of GMU 9D, were converted to contour transects by randomly selecting an end point and flying the elevational contour of that point towards the other end of the transect for 20 km. In the study area, brown bear habitat is constrained by elevation, but the upper elevations of bear habitat in these 2 GMUs was unknown. For GMU 10 (latitude 54.82, longitude 164.14), transect selection was restricted to 610 m or lower, and for GMU 9D to 762 m or lower.
Mark-resight data was collected during the distance sampling survey. A screen separated the pilot and backseat observer, and a light system was used to indicate a potential bear sighting. When a potential bear was observed by the pilot or backseat observer, they turned their light on but did not communicate the sighting. When the potential bear was 5 seconds past the trailing wing edge, the sighting was announced, and the lights checked. To get credit for sighting the bear your light must be on. The mark-resight data will be pilot only, backseat observer only or both; essentially each crewmember was marking the bear for the other. After a sighting was announced the plane was flown off transect; if the object was identified as bear(s), a GPS location was collected. The flight path used to get the GPS location was flown parallel to the transect to minimize measurement error on calculating the closest distance to the transect from this location. A location of how far out the crew was searching when the bear was spotted is determined and overflown to get a search distance location. Distance to bear and search distance were calculated using GIS software (ESRI ArcGIS Desktop. Redlands, CA: Environmental Systems Research Institute) as the closest distance from the transect to that location. Additional details are available in Becker and Quang [12].
Covariates are important for modeling detection probabilities [6,15]. To obtain covariates for each bear group observation, the following data was recorded: bear group size, bear group activity, bear group type, plane speed, percent cover, percent snow, pilot name, backseat observer name, transect type, GMU, and search distance. Bear group activity was recorded as bedded, feeding, standing, walking, or running. Bear group type was recorded as subadult, adult (sex unknown), boar, sow, breeding pair, sow with cubs, sow with yearlings, or sow with 2-year old offspring. Plane speed was the average speed of the aircraft in the 0.5 km transect segment prior to spotting the bear group, it was calculated using the aircraft location and time data. Percent cover included all obscurations and obstructions, with vegetation being the primary cause, within 10 m of the bear group; the crew referred to a reference sheet that depicted percent cover to the nearest 10% (Fig 3). Transect type was recorded as straight or contour. GMU was recorded as 9 for GMU 9D and 10 for Unimak Island.
Some covariates are very difficult to quantify and may even interact with other covariates. Additionally, measuring covariates on the correct scale can be problematic; a bear in 50% vegetative cover within 10 m on an open hillside is easier to detect than one in the same vegetative cover on a hillside with hundreds of patches of 50% vegetative cover. The quality of lighting is hard to quantify but can dramatically affect detection. To account for macro scale and other unmeasured covariates, such as light conditions, we used the search distance covariate. The idea is that if conditions are "good" you can generally search farther than when conditions are average.
Group size was a potential covariate to guard against size-biased sampling [11,16]. Exploratory data techniques, such as stem and leaf plots, LOESS plots, scatter plots, and group summary statistics [17] to identify potentially useful covariates for MCDS and MR modeling. These covariates included modifying the collected covariate data into new covariates. For example, the bear group activity level, bedded bears may have shorter detection distances than other bear activity levels and a new covariate would be created to reflect this. Potential covariates were rescaled to have values which ranged from -5 to 5 to facilitate numerical optimization (J. Laake pers. comm.).
For the MCDS modeling exploratory data analysis identified group size, Pcvr5, PcvrBin3, TranFlat, MCDS2GrpPilot, MCDS2GrpObsvr, Bin2SD1000, and GMU as potential covariates. Pcvr5 was calculated as (percent cover/20). Pcvr5Bin3 was coded 1 for 0% cover, 2 for [10% -30%] cover, and 3 for cover greater than 30%. TransFlat was a covariate which recoded transect type data into a binary variable which was 1 for straight and 0 for contour transects. Exploratory data analysis revealed differences in the distribution of bear detection distances among both pilots and observers. For pilots (n = 5) we created the variable MCDS2GrpPilot which was coded 0 if pilot group = A (n = 1) and 1 if pilot group = B (n = 4). For observers (n = 6) we created MCDS2GrpObsvr which was coded 0 if observer group = A (n = 5) and 1 if observer group = B (n = 1). Exploratory data analysis suggested search distance would be a good potential covariate. In order the met the assumption of independence between covariates and distance for MCDS models [15], we binned search distance into 2 bins and created the covariate Bin2SD1000. We coded Bin2SD1000 as 0 if search distance was 1000 m (w) or less, and 1 if it was greater than 1000 m this covariate was denoted as Bin2SD1000. To consider an interaction between IGrMew and Bin2SD1000 we created the covariate ImewBinSD1000 taking the product of the 2 covariates (IGrMew x Bin2SD1000).
For the MR modeling exploratory data analysis identified observer, distance, bed, Pcvr5, MR2GrpPilot, and MR2GrpObsevr. Conditional probability plots by bear activity suggested bedded bears were harder to detect. We created the covariate bed, which was coded 1 for bedded bears and 0 otherwise. Scatterplots of conditional probabilities of detection by pilots and observers suggested differences among each group. For the pilots we created the covariate

PLOS ONE
Estimating brown bear abundance and harvest in Alaska MR2GrpPilot, which was coded 0 for pilot group a) (n = 1), and 1 for pilot group b) (n = 4). For observers we created the covariate MR2GrpObsevr, which was coded 0 for observer group a) (n = 5) and 1 for observer group b) (n = 1). The covariate observer, coded 0 for pilots and 1 for backseat observers was always used in conjunction with either MR2GrpPilot or MR2GrpObsevr in order to avoid the assumption that pilot group a) and backseat observer group a) had the same conditional probability of detection. Sample sizes of categorical covariates for the MRDS models are listed in Table 1.
Natural breakpoints around the 95 th percentile of the distance data was used to truncate the dataset and determine the truncation distance, w [4]. The two-piece normal detection function was used to fit the MCDS model with parameter μ denoting the mode of the distribution and σ used as a shape parameter that also incorporated covariates into the model [2]. The covariate IGrMew was used to model the asymmetry of the two-piece normal distribution and was coded 1 if the detection distance exceededm, 0 otherwise. The effect of search distance may occur only to the right of the detection apex, in which case an interaction between IGrMew and search distance would be required.
A chi-square goodness of fit test [4], a Kolmogorov-Smirnov goodness of fit test, and a q-q plot [18] were used to assess the fit of the MCDS model to the distance data. A chi-square goodness of fit test was used to assess the fit of the MR model to the MR data [6]. Marques and Buckland [19] recommend avoiding models with estimated detection probabilities less than 0.1 or if more than 5% of these probabilities are less than 0.2. On all models with nonsignificant (α>0.05) goodness of fit statistics and which had estimated detection probabilities that met the recommendations of Marques and Buckland [19], AIC [20] was used to select the most supported MCDS and MR models for GMU 9D, GMU 10 and the combined data. A forward selection process was used and change AIC < -2 to add an additional covariate [20]. AIC was used to determine if the combined model fit better than the 2 separate models.
Population abundance estimates were obtained using MRDS formulae [6,11] that use a 2-piece-normal detection function [2] and assume point independence [6]. The S1 File was used to perform the analysis. The functions used in the analysis are contained in S2 File, which uses a 2-piece-normal detection function in an MRDS model and allows for stratified sampling to obtain the population estimate and standard errors that include variation in the encounter rate (ER, [number of bears detected/number km transect flown]) [21]. For the combined GMU 9D and 10 analysis, we treated the 2 GMUs as strata to obtain individual GMU estimates. For the separate analyses, the whole dataset is partitioned into 2 subsets (GMU 9D and GMU 10), that was both exhaustive and mutually exclusive of the whole dataset; which allows the AIC statistics for GMU 9D and GMU 10 models to be added and compared to the combined model [20] for both the MCDS and MR models. Brown bear hunting during the survey period was by biennial registration permit in GMU 9; seasons were open every other regulatory year (RY) (RY01 = 1 July 2001 through 30 June 2002). This GMU 9 restriction for bears, a unique harvest scheme in Alaska, was established nearly 50 years ago to substantially limit harvest. Consequently, annual harvest has been reported as the average for a 2-year hunt period: (closed year + open year)/2. Fall and spring seasons were limited to 2-3 weeks in length. Bear hunting on Unimak Island (GMU 10) was by limited-entry drawing permit with 8 permits available annually. The bag limit for GMU's 9 and 10 was 1 bear every 4 years. Hunters were required to report their hunting activity including location of kill on a drainage-specific basis. Skulls and hides must be salvaged and presented to ADFG for data collection and tagging with identification. This sealing process was required for transportation of hides and skulls out of state and was therefore a reliable method for monitoring brown bear harvest. Legal (reported) harvest rate was calculated as: (number of legally harvested bears)/(estimated number of bears + number of legally harvested bears). Bear harvest occurred prior to the survey and must be accounted for in the post-hunt population estimate. A 95% harvest rate confidence interval was constructed by dividing the reported harvest by the 95% confidence limits of the population abundance estimate. Harvest rate standard error was calculated by applying the coefficient of variation to the estimated harvest rate. Illegal (unreported harvest) was unknown and not accounted for in our estimates.

Results
Weather did not allow surveys to be conducted on 1, 9, and 10 June; survey effort was reduced to 3-4 planes on 4, 8, and 11 June due to weather, pilot illness and 2 pilots leaving for other job commitments. Daily local weather within the 5 zones of the study area curtailed survey efforts on an infrequent basis. Not accounting for ferry time to and from the study area, our survey was conducted with approximately 38 plane days of effort. A total of 609 20-km long transects were flown during the survey, with 370 in GMU 9D and 239 in GMU 10. For GMU 9D, no brown bears were detected above 610 m; and for GMU 10 none were detected above 305 m, which we used to define spring brown bear habitat. GMU 9D brown bear habitat consisted of 8,648.2 km 2 and was sampled with 363 transects with a total length of 7,436.4 km from which 253 bear groups containing 458 bears were detected (Fig 1). GMU 10 brown bear habitat was 2,868.7 km 2 and was sampled with 219 transects with a total length of 4,393.3 km from which 106 bear groups containing 174 bears were detected (Fig 2).
A stem and leaf plot of the distance data [17] indicated 1000 m was a natural breakpoint in the vicinity of the 95 th percentile of the distance data. We used 1000 m for the truncation distance (w) and 22 m for the edge of the blind spot underneath the aircraft (w b = 22, [1]) for the GMU 9D, GMU 10, and the combined dataset; this truncated 5.2% (n = 14) of the GMU 9D data, 5.4% (n = 6) of the GMU 10 data, and 5.3% (n = 20) of the combined GMU 9D and 10 dataset.
Two-hundred and fifty-three bear groups containing 458 bears were detected within the transect searched area (22-1000 m strip next to transect) in GMU 9D, with 204 groups detected by the pilot and 184 detected by the backseat observer, and 135 of these groups were detected by both. For GMU 9D, the most supported MCDS model (AIC = 3262.41) contained an intercept parameter, IGrMew, IMewBin2SD1000, and TranFlat (Table 2). Based on QQplots and goodness of fit tests (χ 2 = 10.942, df = 10, p = 0.362; Kolmogorov-Smirnov Statistic p = 0.902) this model fit the data well. The distance of apex of detection (m) was estimated to be 110.6 m. Model fit at the apex of detection is very important [3]. A good fit was obtained with this model (0.860 contribution to the χ 2 statistic of 10.942) for the bin containing the apex. The most supported MR model (AIC = 502.96) contained parameters for the intercept, bed, and Pcvr5 (Table 3). This model fit the MR data well (χ 2 = 23.558, df = 21, p = 0.315) and estimated detection at 110.6 m (p � ðmÞ) to be 0.872 (SE = 0.024). The combined MRDS model (Fig 4)  One-hundred and six bear groups containing 174 bears were detected within the transect searched area in GMU 10, with 88 groups detected by the pilot and 80 detected by the backseat observer, and 62 of these groups were detected by both. For GMU 10, the most supported MCDS model (AIC = 1388.85) contained an intercept parameter, IGrMew, Bin2SD1000, and TranFlat (Table 2). Based on QQ-plots and goodness of fit tests (χ 2 = 2.340, df = 4, p = 0.673; Kolmogorov-Smirnov Statistic p = 0.275) this model fit the data well. The distance of apex of detection (m) was estimated to be 111.8 m. A good fit at the detection apex was obtained with this model (0.112 contribution to the χ 2 statistic of 2.340) for the bin containing the apex. Our MR data indicated the most supported MR model (AIC = 197.18) contained parameters for an intercept, distance, and bed (Table 3). This model fit the MR data well (χ 2 = 14.518, df = 13, p = 0.338) and estimatedp � ðmÞ at 111.8 m to be 0.952 (SE 0.026). The combined MRDS model (Fig 4) estimated 309.9 brown bears (SE = 53.2), CV = 17.2%, or 108.03 brown bears/1000 km 2 (SE = 18.53) in GMU 10.
We fit an MRDS model to the combined GMU 9D and GMU 10 dataset. The most supported MCDS model AIC = 4644.75) contained an intercept parameter, IGrMew, Bin2SD1000, and TranFlat (Table 2, Figs 5 and 6). Based on QQ-plots (Fig 7), and goodness of fit tests (χ 2 = 11.566, df = 13, p = 0.563; Kolmogorov-Smirnov Statistic p = 0.645) this model fit the data well. The distance of apex of detection (m) was estimated to be 110.5 m. A good apex fit was obtained with this model (0.437 contribution to the χ 2 statistic of 11.566) for the bin containing the apex. The most supported MR model (AIC = 692.14) contained parameters for an intercept, bed, observer, MRPilotGrp2, distance, Pcvr5, and MR2GrpObsevr (Table 3) (Table 4) than models based on the individual GMU data. The precision gain was slight for GMU 10 (change CV = 0.0192, 1.92%) and miniscule for GMU 9D (change CV = 0.0024, 0.24%). The precision of the combined population Table 3 Calculating bear harvest rate for GMU 9 is complicated by hunts that are primarily open every other regulatory year as a means to limit hunter harvest. The reported hunter harvest for GMU 9D was 153 bears (69 in fall and 84 in spring) for RY01 and zero in RY02 (Table 5). Average annual harvest for these 2 years was (153 + 0)/2 = 76.5 bears. The resulting harvest rate was 4.35% (se = 0.45%, CV = 10.4%, Table 5) of the estimated population size. The estimated harvest rate in GMU 10 for RY01 was 3.06% (se = 0.47%, CV = 15.2%, Table 5).

Discussion
The goal of this survey was to estimate brown bear population abundance on a sub-GMU basis to calculate harvest rates and confidence intervals. To our knowledge, this was the first time in Alaska that brown bear abundance was estimated on the same geographic scale as was harvest data, which allowed for harvest rates and confidence intervals to be calculated. The point estimates for GMU 9D and 10 were below 5% indicating a that harvest level supporting trophy-  nonresident hunters in spring 2020. Estimating abundance of bears in GMU 9D again using MRDS while harvest is historically low would be of benefit in monitoring the bear population and evaluating this technique.
A secondary goal of this work was to determine if data from the 2 surveys could be combined and if so, the resulting precision gain. The MRDS model with combined data from GMUs 9D and 10 were more supported (AIC) than individual models for each GMU. The combined model for these data produced precise (CV = 10%) results for GMU 9D (CV = 10.4%) and less precise estimates for GMU 10 (CV = 15.2%). Fewster et al. [21] states "the dominant source of variance in distance sampling is usually the encounter rate". Our results support this observation since the estimated CV of the encounter rates for GMUs 9D and 10 (8.9% and 13.5% respectively) were only slightly smaller than the CV of the population estimate, which incorporates model variance and encounter rate variance. This would explain the small gains in precision we obtained from using a combined model versus single GMU based models; since the estimates of encounter rate and variance are calculated within GMU, and only model variance would differ between the two approaches. This survey had the same

PLOS ONE
truncation distance, crews, and survey dates the survey was conducted differ among surveys which usually results in it being inadvisable to combine these data.
Our MRDS model fit the combined data and had parameter estimates that make biological sense. Our MCDS model indicated that detection probabilities increased when conditions   allowed for farther distances to be searched more effectively (Bin2SD1000 = 0.6620, Table 2), similarity straight line transects were more efficient (TranFlat = 0.3724, Table 2) for detection than contour transects in mountainous terrain. Our MR model estimated apex detection at 0.920 (0.019, Table 3) indicating the need to estimate this parameter and how erroneous an assumption of perfect detection would be. The need to estimate apex detection has been documented in other distance sampling surveys of bears [22]. The MR model also indicated bedded bears were harder to detect than more upright bears (bed = -1.1194, Table 3), as distance increased more bears were missed (distance = -0.0013, Table 3), and bears were harder to detect as percent cover increased (Pcvr5 = -0.2971, Table 3). By adding the coefficients for observer, MRPilot2Grp and MRObsvr2Grp (Table 3), we obtain the following cumulative detection coefficients: Pilot Group a (0), Backseat Group a (0.2977), Pilot Group b (0.8515), Backseat Group b (0.2977+1.0011 = 1.2988, Table 3). Pilots have a better view than backseat observers; however, one pilot (Pilot Group a, n = 1) was suffering from the flu during the survey, additionally one backseat observer (Backseat Group b, n = 1) did exceptionally well detecting bears. When the data was collected, the plan was to obtain a population estimate with an MRDS model based on a gamma shaped detection function [12]. In the process of analyzing the data we were puzzled how a MCDS model could obtain a bigger abundance estimate than one from our MRDS model. Eventually we read and understood Borchers et al. [6] that indicated the full independence assumption of MRDS models was the cause of this bias. It made no sense to publish an analysis which was biased, but we did not have a detection function that could implement the recommended solution of an MRDS model that assumed point independence [6]. A gamma detection function has an apex for every covariate level [2] while the assumption of point independence allows only 1 apex [6]. At this point there was no existing solution on how to model this data. In 2011, the senior author and S. Buckland (University of St. Andrews) collaborated on an effort to use a Level of independence modeling approach [23] based on a logistic detection function for modeling the bear detection data. Initially, progress was made, but ultimately the logistic function was not flexible enough to model the shape of the bear detection data. In desperation, the senior author used a half normal detection function to fit the data, the cost was an additional left truncation of 11.6% to 27.2% of the data which was deemed very inefficient [22]. The left truncated data was also a half normal distribution with a different scale parameter. This insight lead to the development of the two-piece normal detection model [2]. The r package mrds [24] which performs MRDS modeling could not be easily adapted to the two-piece normal detection function; so major r-code writing was required to fit an MRDS model that uses a two-piece normal detection function, graph model results, and perform goodness of fit tests. The complexity of the r-coding coupled with lack of r-coding expertise by the senior author slowed progress in reporting an MRDS model based on a two-piece detection model 2]. The variance in the two-piece normal MRDS model [2] was based on Borchers et al. [6]. The ER variance work of Fewster et al. [21] was not incorporated into the initial MRDS model based on a two-piece normal detection function [2]. The rcode in this paper corrects that error and can analyze stratified data. To summarize, the 2002 data was collected just prior to a time of rapid advancement in MRDS models which caused delays in obtaining a model that was suitable for these data; the complexity of the r-coding caused additional delays in the analysis, this delayed analysis addresses the issues raised above to obtain scientifically sound abundance and harvest rate estimates. The GMU 9D and GMU 10 survey took 38 plane days to accomplish and detected 359 bear groups that were used in the analysis. Currently survey aircraft with gas charge $375 an hour, assuming 8 hours of flying per day and 38 plane days puts the survey cost at $114,000. The lack of knowledge about brown bear distribution in the study area resulted in a simple random sample (SRS) survey design which is an inefficient design [13].
Survey cost is a major impediment to using distance sampling for brown bear management. The individual GMU 9D and 10 results indicate 253 detections and 106 detections from an inefficient SRS survey resulted in CV's of 10.60% and 17.20% respectively (Table 4). Survey costs can be managed and reduced by: 1) stratifying on encounter rate; 2) reducing the number of detections for the desired sample size for planning a survey; 3) use better survey designs to more uniformly spread transects out within strata, 4) lowering the desired precision, 5) conducting the survey over 2 or 3 years to make the cost more manageable and treat the abundance estimate as an average [12], and 6) sharing costs with the USFWS when surveying National Wildlife Refuges such as GMU 9D and Unimak Island. Successful stratification will require pre-survey information form a previously done survey or a pilot study [25]. We stratify on encounter rate versus density [26] since encounter rate is responsible for 68% to 90% of the CV of abundance from SRS distance sampling surveys (E.Becker unpublished data). Examination of Fig 1 indicate large geographic areas with lower ERs than other areas. ERs by 250 feet elevational bands ranges from 0.050 to 0.092 for GMU 9D and 0.028 to 0.046 for GMU 10. Clearly ER can be stratified for these areas. It is hard to determine the effect of stratified designs on precision which would impact the setting the desired target sample size, The results of SRS design on the GMU 9D and 10 induvial estimates suggest 200 to 250 detections is a good target range for planning a survey and can be refined once more data is available. Spreading transects out in a more uniform coverage pattern would help with survey precision [25], the use of zig-zag transects [26] in U-shaped river valleys of interior Alaska and constricted flat areas like the narrow portions of GMU 9D (Fig 1) would allow for consistent survey effort in these areas. Precision goals of 10-12% CV seem reasonable.
We should be able to reduce sample size which would reduce precision and gain that precision back with a good stratification plan. Stratifying the study area based on encounter rate by elevational band and sampling proportional to the encounter rate we crudely estimated 65% of our 2002 survey effort would yield about 250 group observations to build a detection model. We hypothesis that the modeling variance would increase slightly, and the encounter rate variance would be significantly decreased. If this conjecture is true, a stratified survey would of GMU 9D and 10 would take 24.7 plane days and cost approximately $74,000; we would anticipate a 10% to 12% CV for GMU 9D and GMU10. Incorporating spatial components into the stratification would improve the efficiency of a new survey and possibly farther reduce costs.
The r-package dssd [25] can implement transect study designs with both straight and zigzag transects and also depict transect coverage, currently contour transects are not supported. Density surface models [27,28] can be used to increase the precision of the abundance estimates and estimate abundance in sub-regions.
Managing brown bear populations in large remote areas is extremely difficult without reliable abundance estimates. An unbiased and precise estimate of population abundance on the same geographic scale on which harvest data is collected allows for precise estimates of reported harvest rates. The ability to estimate harvest rates is an invaluable tool for wildlife managers to determine if changes to hunting regulations need to be made relative to over or under harvest. Stratified aerial distance sampling with MR data should be considered for these areas. Care must be taken to ensure these surveys can met assumptions. We anticipate costs can be reduced and efficiency increased in future surveys.
Supporting information S1 File. An r-script file containing the r-code used in the analysis of this manuscript. The r-packages mrds [24] and optimx [29,30] are required to run the analyses. (R)

S2
File. An r-script file containing the functions called by GMU9D10_MRDS_Analysis.r. This code is an improvement of the r-code found in [2]. The most important improvements are the incorporation of encounter rate variance into the population estimate and the ability to use stratified designs. Additional refinements are documented by comments in the -r-code. The r-package mrds [24] is required. (R) S3 File. A csv file containing the brown bear distance data that was analyzed in this manuscript; this file is required by GMU9D10_MRDS_Analysis.r. (CSV) S4 File. A csv file containing the transect data that was analyzed in this manuscript; this file is required by GMU9D10_MRDS_Analysis.r. (CSV) data. We are grateful to J. Stetz, ADF&G, for his helpful comments and suggestions on an early draft of this manuscript. We thank Academic Editor T. Mousseau and 2 anonymous reviewers for their many helpful suggestions. We thank The Centre of Research into Ecological and Environmental Monitoring, University of St. Andrews, and ADF&G Commissioners C. Campbell and D. Vincent-Lang for supporting the senior authors' collaboration efforts with S. Buckland.