A cost efficient spatially balanced hierarchical sampling design for monitoring boreal birds incorporating access costs and habitat stratification

Predicting and mitigating impacts of climate change and development within the boreal biome requires a sound understanding of factors influencing the abundance, distribution, and population dynamics of species inhabiting this vast biome. Unfortunately, the limited accessibility of the boreal biome has resulted in sparse and spatially biased sampling, and thus our understanding of boreal bird population dynamics is limited. To implement effective conservation of boreal birds, a cost-effective approach to sampling the boreal biome will be needed. Our objective was to devise a sampling scheme for monitoring boreal birds that would improve our ability to model species-habitat relationships and monitor changes in population size and distribution. A statistically rigorous design to achieve these objectives would have to be spatially balanced and hierarchically structured with respect to ecozones, ecoregions and political jurisdictions. Therefore, we developed a multi-stage hierarchically structured sampling design known as the Boreal Optimal Sampling Strategy (BOSS) that included cost constraints, habitat stratification, and optimization to provide a cost-effective alternative to other common monitoring designs. Our design provided similar habitat and spatial representation to habitat stratification and equal-probability spatially balanced designs, respectively. Not only was our design able to achieve the desired habitat representation and spatial balance necessary to meet our objectives, it was also significantly less expensive (1.3−2.6 times less) than the alternative designs we considered. To further balance trade-offs between cost and representativeness prior to field implementation, we ran multiple iterations of the BOSS design and selected the one which minimized predicted costs while maximizing a multi-criteria evaluation of representativeness. Field implementation of the design in three vastly different regions over three field seasons showed that the approach can be implemented in a wide variety of logistical scenarios and ecological conditions. We provide worked examples and scripts to allow our approach to be implemented or adapted elsewhere. We also provide recommendations for possible future refinements to our approach, but recommend that our design now be implemented to provide unbiased information to assess the status of boreal birds and inform conservation and management actions.

Predicting and mitigating impacts of climate change and development within the boreal biome requires a sound understanding of factors influencing the abundance, distribution, and population dynamics of species inhabiting this vast biome. Unfortunately, the limited accessibility of the boreal biome has resulted in sparse and spatially biased sampling, and thus our understanding of boreal bird population dynamics is limited. To implement effective conservation of boreal birds, a cost effective approach to sampling the boreal will be needed. Our objective was to devise a sampling scheme for monitoring boreal birds that would improve our ability to model species-habitat relationships and monitor changes in population size and distribution. A statistically rigorous design to achieve these objectives would have to be spatially balanced and hierarchically structured with respect to ecozones, ecoregions and political jurisdictions. Therefore, we developed a multi-stage hierarchically-structured sampling design known as the Boreal Optimal Sampling Strategy (BOSS) that included cost constraints, habitat stratification, and optimization to provide a cost-effective alternative to other common monitoring designs. Our design provided similar habitat and spatial representation to habitat stratification and equal-probability spatially balanced designs, respectively. Not only was our design able to achieve the desired habitat representation and spatial balance necessary to meet our objectives, it was also significantly less expensive (1.3−2.6 times less) than alternative designs we considered. Field implementation of the design in three vastly different regions over three field seasons showed that the approach can be implemented in a wide variety of logistical scenarios and ecological conditions. We provide worked examples and scripts to allow our approach to be implemented or adapted elsewhere. We also provide recommendations for possible future refinements to our approach, but recommend that our design be implemented to provide unbiased information to assess the status of boreal birds and inform conservation and management within the boreal biome.

Abstract
Predicting and mitigating impacts of climate change and development within the boreal biome requires a sound understanding of factors influencing the abundance, distribution, and population dynamics of species inhabiting this vast biome.
Unfortunately, the limited accessibility of the boreal biome has resulted in sparse and spatially biased sampling, and thus our understanding of boreal bird population dynamics is limited. To implement effective conservation of boreal birds, a cost effective approach to sampling the boreal will be needed. Our objective was to devise a sampling scheme for monitoring boreal birds that would improve our ability to model specieshabitat relationships and monitor changes in population size and distribution. A statistically rigorous design to achieve these objectives would have to be spatially balanced and hierarchically structured with respect to ecozones, ecoregions and political jurisdictions. Therefore, we developed a multi-stage hierarchically-structured sampling design known as the Boreal Optimal Sampling Strategy (BOSS) that included cost constraints, habitat stratification, and optimization to provide a cost-effective alternative to other common monitoring designs. Our design provided similar habitat and spatial representation to habitat stratification and equal-probability spatially balanced designs, respectively. Not only was our design able to achieve the desired habitat representation and spatial balance necessary to meet our objectives, it was also significantly less expensive (1.3−2.6 times less) than alternative designs we considered.
Field implementation of the design in three vastly different regions over three field seasons showed that the approach can be implemented in a wide variety of logistical

Introduction
Tackling ongoing [1][2][3] and projected [4,5] global biodiversity losses will require difficult decisions about resource allocation and the need to consider when, where, and how to focus conservation efforts. Species conservation is often most successful and least expensive when implemented early, before species require dedicated recovery efforts [6]. Early identification of species declines should improve the likelihood of successful conservation. In addition, developing effective conservation strategies requires identifying the locations and factors contributing to species declines [7]. Sparse or biased data may inaccurately identify which species require conservation actions or incorrectly identify key drivers, resulting in misdirected or ineffective conservation actions [8]. Thus, well-designed ecological monitoring is necessary to effectively prioritize conservation activities [9][10][11].
Despite the need for monitoring data, even comparatively well monitored taxa such as terrestrial birds have significant gaps in species coverage that hinder effective status and trend assessments and associated conservation actions [12]. In North America, many species of terrestrial birds are primarily monitored using the North American Breeding Bird Survey (hereafter BBS; [12]). However, this survey, which is based on roadside surveys, mainly run by volunteers, has very limited coverage in remote northern locations such as the boreal forest [13,14]. Furthermore, much of the coverage for species that breed in the boreal forest is limited to the southern edges of their range, where population change may be quite different from elsewhere. The lack of data from most of the boreal forest may present significant risks to conservation given expanding resource development [15][16][17] and projected climate change impacts on boreal birds [18].
Monitoring data from northern ecoregions are needed to inform conservation actions (e.g. prioritize selection for protected areas), contribute to management decisions (e.g. adaptive management), detect range shifts, assess threats to species (e.g. species response to developments), and meet legislative requirements. In Canada, monitoring data are required to support conservation under the Migratory Birds Convention Act (1994, c. 22), which includes informing listing decisions under Canada's Species at Risk Act (hereafter SARA; [19]), impact assessments required by the Impact Assessment Act [20], or assessments required by provinces or territories. Key criteria used by the Committee on the Status of Endangered Wildlife in Canada (hereafter COSEWIC) to recommend listing a species under SARA include decline in the total number of mature individuals, population size, and extent of species occurrence [21]. In addition, habitats and region specific estimates of species' abundance or density are required baseline data in environmental impact assessments [20]. To support these diverse requirements, multi-species bird monitoring should be capable of accurately capturing not only trends in species' population sizes, but also changes in both spatial and temporal variation in species status and distribution.
Given the breadth of these objectives and the large and remote areas to be monitored a cost effective sampling approach is required. As most of these locations have limited road networks [22] with few human settlements, there is little potential for volunteer surveys [23]. Similarly, a simple, design-based sampling strategy (e.g., a randomized design in which sample units are selected with the same probabilities) is logistically and financially impractical [9,24]. Instead, we propose a design-based sample design that integrates multiple approaches from the statistical sampling-theory literature including spatially balanced sampling [25,26], stratified sampling [27][28][29][30], costconstrained sampling [29] and optimization [31]. The integration of these approaches optimizes the spatial balance and habitat representation at minimal program costs, for the entire boreal landbird community.
Here, we describe a cost-effective hierarchically structured unequal probability sampling design known as the Boreal Optimal Sampling Strategy (BOSS). We compared our design against three competing spatially balanced sampling designs: a purely spatial (equal-probability) design, a habitat stratification design, and a design that minimized access costs. We hypothesized that the BOSS design should provide the second best habitat representation after the habitat stratification design since we incorporated habitat in the inclusion probabilities. We also hypothesized that our BOSS design should rank as the second-least expensive alternative after the cost design owing to incorporation of costs within the inclusion probabilities. Specifically, our objectives are to: (1) outline the development of the BOSS sampling design; (2) describe BOSS design implementation in three ecologically distinct regions in Canada (Newfoundland and Labrador, Saskatchewan, Yukon Territory); (3) evaluate how the BOSS design meets spatial balance, habitat representation, and cost objectives in these three distinct regions; and 4) evaluate implementation costs in field trials.

Sampling frame
We defined the target population as all adult individuals of terrestrial bird populations occupying an ecoregion, province/territory (hereafter jurisdiction), and Bird Conservation Region (BCR) within the boreal region of Canada [32] during the breeding season for a particular year. We defined the Primary Sampling Unit (PSU) as a 5 km diameter hexagonal grid cell. We generated a hexagonal grid of PSUs covering the entire sampling frame: the boreal region of Canada, plus a 100 km northern buffer ( Fig   1). We used hexagonal geometry because it introduces fewer edge effects and performs well in nearest neighbour analyses since the distances between centroids are the same in all directions [33]. We included a northern buffer to allow the design to capture possible breeding range shifts associated with climatic changes in the north [34]. Each PSU was assigned a unique number/identifier. With the exclusion of PSUs that have cumulative inclusion probabilities of zero (e.g. exclusively open water), all PSUs in the sampling frame were available to be sampled and as a result, the BOSS sampling design can make valid inferences to entire populations of species inhabiting the boreal [35].

Stratification
We developed a hierarchical stratification scheme to estimate population variables at multiple spatial scales. At the first level, we stratified the sampling frame by the intersection of BCRs and jurisdictions (e.g., area of Yukon Territory within BCR 4-Northwestern Interior Forest). At the second level, we further stratified each BCR by jurisdiction combination by ecoregion (e.g., the area of Selwyn Mountains ecoregion within Yukon Territory within BCR 4), where ecoregions are a subdivision of an ecozone characterized by similarity in surface forms, flora, fauna, hydrology, soil, and macro climate [36]. Regional monitoring programs as part of the BOSS design are defined by the area of each jurisdiction within a BCR. We set target sample sizes for the number of PSUs to select within an ecoregion using the process outlined below (see sample size allocation).

Sample size allocation
Optimal stratified sampling accounts for both the size of the stratum (e.g., landarea for geographic regions), and the variance within the stratum [37]. Greater sampling effort is allocated using Neyman allocation to strata that are larger and more variable to optimize the precision of the overall estimate [27,29]. Specifically, we calculated the sample size of PSUs (n) for a given jurisdiction using the following equation to allocate sampling effort between all (k) ecoregions within the jurisdiction: Where ne = is the target sample size for a given ecoregion within a given jurisdiction, N = total target sample size to sample within the jurisdiction, Areae = area of ecoregion e, σe= is a measure of observed or predicted variation within ecoregion e for the variable(s) of interest, and re = described species richness of the bird community in ecoregion e (i.e. the number of all bird species with ranges overlapping the ecoregion derived from overlays of species range maps with the ecodistricts of Canada [36] by GeoInsight Corporation [38]). Since our goal is to monitor multiple species, employing a univariate measure of variance would provide an inadequate representation of variance [29]. We included species richness as a weighting factor on within-stratum variation (σe), to increase sampling in ecoregions where more species occur and improve the (1) likelihood of adequately monitoring a greater number of species, and (2) ability to monitor and predict species distributions.
Many of the ecoregions in our sampling frame have little or no historic sampling from which we could derive estimates of spatial or temporal variance in the bird community. As an alternative, we used proxy variables that correspond with the spatial and/or temporal variance in bird abundance within an ecoregion to set relative PSU sample size targets for the ecoregion. We analyzed the proxy variables [29] representing spatial and temporal variation in bird abundance and distribution by calculating a measure of multivariate dispersion (variance estimate) for each ecoregion following [39,40]. We used several proxy variables to calculate multivariate dispersion given well described links between species abundance and distribution and spatial and temporal climatic variation [41], vegetation [41,42] and variation in forest age [42,43].
Specifically, we used the following proxy variables: (1) (6) percent area of open water since riparian edges significantly increase avian species richness [45]. We calculated variance in May-July (1) temperatures and (2) precipitation from monthly gridded (0.5 x 0.5 degrees) values from the climate datasets of Harris et al. [46]. We calculated the between year standard deviation of burn rates from data obtained from [47]. Since burn rates from [47] were calculated in sample units of 10,000

Sampling design
We developed a two-stage sampling design in which the first stage of the design involved selection of PSUs from ecoregions within the sampling frame and the second stage involved selection of Secondary Sampling Units (SSU) within the selected PSUs ( Fig 2). The SSU represent on-the-ground locations where point count surveys [48] will be conducted using trained observers and/or Autonomous Recording Units [49].

Selection of primary sampling units
To select PSU within each jurisdiction (Stage 1) we first needed to develop an approach to combine the following three design components: access costs, habitat representation, and spatial balance. We followed the steps outlined below in each of our three jurisdictions (Newfoundland and Labrador, Saskatchewan, and Yukon Territory) to integrate all components of our sampling design in areas spread widely across our sampling frame (Fig 3). First, we created separate cost surfaces for logistically feasible access methods within each of our jurisdictions. Specifically, we developed spatially explicit models predicting the cost to access a given PSU assuming access from roads (trucks), helicopters, float planes, snowmobile and/or all terrain vehicles as appropriate for a given jurisdiction (see examples in Supporting Information S2). The resulting cost models were represented as pixel based (250 x 250m) estimates of cost, from which we subsequently calculated the average cost of access for all pixels within a PSU separately for each access type (trucks, helicopters, float planes). We assumed that the lowest cost access method will generally be used to access a given PSU and calculated the minimum cost across all access cost surfaces. We subsequently used this value as the estimated access cost for each PSU.
Second, we calculated inclusion probabilities based on landcover classes (hereafter habitat classes). We calculated habitat-based inclusion probabilities in an attempt to achieve balanced representation of habitat classes within an ecoregion.
Inclusion probabilities for habitat classes used a weighted sampling approach to balance sampling across common (high occurrence/high proportion) and uncommon (low occurrence/low proportion) habitats within an ecoregion. Weighted habitat inclusion probabilities were calculated as follows: Where the inclusion probability for a given pixel p of habitat class i within a given ecoregion (phabitat ice ) was a function of the inverse of the number of habitat classes within the ecoregion divided by the area of habitat i within ecoregion e. Our weighted habitat inclusion probability thus results in the cumulative probability of sampling a given habitat within an ecoregion being equal among all habitat classes despite differences in habitat area.
Third, we combined access costs with weighted habitat inclusion probabilities to derive our final summed value for inclusion probabilities for each PSU as follows: Where PSU k = the sum of all inclusion probabilities within a given primary sample unit (PSU) k, Phabitat pie = area weighted selection probabilities for habitat class i within pixel p in ecoregion e, and cost for a given PSU was the minimum access cost for sample unit k (as described above). Cost was included as the inverse square root of costs to preferentially (all else being equal) favour sampling at sites that are less expensive to sample in accordance with optimal allocation theory [50].
Fourth, we used unequal-probability sampling based on the calculated inclusion probabilities (PSU k ) to draw a spatially balanced sample within each ecoregion using Generalized Random Tessellation Stratified sampling (GRTS; [25]). The random draw to select PSU within ecoregions in each jurisdiction uses PSU sample size targets for each ecoregion (where ecoregion boundary was clipped to the jurisdiction boundary).
We included a 20% oversample of primary sampling units for each ecoregion to allow substitution of inaccessible or unsafe PSUs (e.g. cliffs or other features not visible in the available input layers). All PSUs were labeled with their draw order and oversample PSUs were additionally labeled as oversamples. Should a PSU be inaccessible or unsafe to survey, the oversample and draw order can be used to select a replacement PSU that retains the spatial balance of the sample [25].
Finally, we derived an approximately optimal draw of PSU by a combination of repeated random sampling and multi-criteria evaluation. Specifically, we propose that multiple random unequal-probability GRTS samples can be obtained for which the predicted sampling costs and metrics of habitat representation and spatial balance can be derived. The iteration resulting in the highest combined representativeness of habitat sampling combined with spatial balance for the lowest predicted sampling cost should provide an approximately optimal design. As a metric of habitat representation, we calculated the sum of squared differences (SSD) between the area of each habitat class within all PSUs in a randomized draw versus the proportional area of each habitat class.
As a metric of spatial balance, we use Pielou evenness (PE in eq. 4 below) to represent the even spread of PSU across the sample frame [28]. In order to evaluate the optimal draw across multiple (1…j) iterations, we propose evaluating which iteration minimizes predicted costs while maximizing the weighted sum of the two metrics: Where SSD for sample draw i is relativized to the minimum SSD across draws 1 through j, and subtract the resulting value from 1 so that the metric scales from 0-1 and thus a value of 1 would represent perfect representation of habitat composition within the draw. We applied equal weights to both the habitat (wh) and spatial balance (ws) components as we did not want either component to be more heavily valued than the other in our final design (i.e. wh = ws = 0.5). This approach can be generalized by allowing weights to vary between habitat representation and spatial representation if the desire was to emphasize habitat representation (e.g. environmental impact assessments) versus distribution modeling or vice versa.

Selection of secondary sampling units
To select secondary sample units (SSU) within each PSU in each jurisdiction (Stage 2), we first created a uniform (systematic) grid of potential secondary sample unit locations at 300 m spacing over the extent of each ecoregion. Each SSU was separated by 300 m since the effective detection radius for most boreal songbirds are less than 150 meters [48] and therefore double counting should be infrequent for most species.
We used the grid of SSU to query the habitat-based inclusion probabilities. We selected a stratified random sample of four SSU where higher inclusion probability values indicating rare/uncommon habitat classes influenced the selection of SSU. We used the four SSU selected during this draw as plot centroids and the eight SSUs surrounding the plot centroid to construct plots of nine SSU for sampling (see Fig 2). We set a sampling target of two plots of nine SSUs per PSU, but drew a sample of four plots of SSU per PSU to create an oversample in case of inaccessibility (e.g. across a large water barrier) or safety concerns (e.g. hazardous terrain or aggressive wildlife). In addition, we aimed for two plots of nine SSUs so that a team of two staff can complete one PSU per day in most circumstances. Furthermore, having two survey plots allows staff to have a partner(s) within the same PSU to act as a first responder if needed. In addition, our PSU size of 5 km diameter was specifically selected to ensure field staff would be working within a reasonable walking or flight distance if assistance or emergency response is required.

Data collection
We

Ethics statement
We obtained the required permits from provinces/territories, Provincial and

Statistical analyses
Prior to analysis of our multi-variate proxy data, we centered and standardized all variables to zero means and unit variance and then calculated Euclidean distances between all points (ecodistricts) in the dataset. We derived the multivariate estimate of dispersion by calculating mean distances to ecoregion centroids using the 'betadisper' function in the vegan package [51] in the R statistical computing environment [52]. The resulting estimates of dispersion for all ecodistricts within the boreal region as per Brandt [32] are available in Supporting Information S1.
In order to assess our BOSS design, we used the 'spsurvey package' [53] within the R statistical computing environment [52] to draw samples under four alternative spatially balanced sampling designs. Specifically, we drew samples under one equalprobability spatially balanced design and three unequal-probability designs: (a) a habitat stratification design, calculated as per equation (2) above; (b) a cost design in which inclusion probabilities were based solely on access costs as per equation (3); and (c) our BOSS design in which inclusion probabilities included both habitat stratification and access costs as per equation (4) above. For each scenario, we ran 100 iterations in which we drew a random sample of n= 400 PSUs from each of the three jurisdictions considered here. In each case, we drew a spatially balanced sample using the GRTS algorithm [53]. Within each iteration, our script calculated the access costs (CDN $), spatial balance (Pielou evenness [28,53]), and habitat representation as SSD (see above). All implementations are available within scripts in the supplemental material (Supporting Information S2).
We used linear models to test for differences (for each jurisdiction separately) in each of our response variables (cost, SSD, and Pielou evenness) between the alternative designs described above (i.e. spatial, habitat, cost, and BOSS). Study design was included as a factor, and in each model we specified an a priori reference factor level to test our hypotheses. We specified the cost design as the reference factor level for models examining access costs, as we anticipated this design should have the lowest access costs on average. Similarly, we specified the habitat stratification design as the reference level for models examining variation in SSD as we predicted this design should have the lowest SSD.
All four designs we considered employed a GRTS algorithm to generate spatially balanced samples; however, the inclusion of unequal inclusion probabilities based on spatially structured variables could alter the ability of the design to achieve spatial balance. We therefore tested for differences in spatial balance between our designs and we treated the spatial design (equal-probability spatial sampling design) as the reference factor level when comparing spatial balance between sampling designs, but did not have a priori predictions. We calculated robust standard errors ('sandwich' estimators) to overcome problems with heteroscedastic errors using the 'lmtest' [54] and 'sandwich' [55] packages within the R statistical computing environment [52]. We applied a Bonferroni correction and considered results statistically significant when α ≤ 0.05/9.

Sample Size Allocation
Sample size allocation to geographic strata (ecoregions) under the BOSS design resulted in a similar sample size on average as allocating samples based on stratum area as data roughly fall along a 1:1 correspondence line (Fig 4). However, allocation to individual strata differed markedly, with the BOSS design suggesting significantly low sample sizes for some strata than an area-based allocation while other strata were allocated larger sample size than would be assigned based on stratum area (Fig 4).

Predicted sampling costs
Predicted costs to access 400 primary sample units varied by the jurisdiction and sampling design considered (Fig 6). Across all three jurisdictions, designs in which the inclusion probabilities derived solely based on access costs (cost design) were the least expensive designs (Fig 6) and all the alternative designs were statistically more expensive regardless of jurisdiction (Table 1). The BOSS design had the next lowest cost (Fig 6), with costs predicted to be between 1.02 (Yukon Territory) and 1.2 (Saskatchewan) times more expensive on average than the cost design for a given jurisdiction (Table 1). In both Newfoundland and Saskatchewan the equal-probability spatial design had the third lowest cost followed by the habitat stratification design, whereas it was the most expensive design for Yukon Territory (Fig 6); it was predicted to be 1.4−2.6 times more expensive on average (~$362,000−562,000; Table 1) than the cost design (maximum of 1.5−3.0 times; Fig 6). The habitat stratification design was predicted to be 1.4−2.0 times more expensive on average (~$122,000−1,020,000; Table   1), or a maximum of 1.6−2.3 times more expensive across draws and jurisdictions ( Fig   6).

Habitat representation
As predicted, the habitat stratification design had the lowest sum of squared differences between area represented within sample draws and the desired equal habitat representation (Fig 7). The BOSS design provided the second most representative sampling of habitat classes on average (Fig 7) with habitat representation in Saskatchewan and Yukon Territory showing no statistical difference, and Newfoundland having statistically lower habitat representation than the habitat design ( Table 2). In contrast, both the equal-probability spatial design and the cost design had large (Fig 7) significantly lower habitat representation relative to the habitat stratification design in all three jurisdictions (Table 2).

Spatial representation
All sampling designs provided relatively balanced spatial sampling of all three jurisdictions (range of Pielou evenness across all designs: 0.96−1.01; Fig 8).
Interestingly, the spatial balance of the cost design provided the closest approximation to the equal-probability spatial design (Fig 8) and was not statistically different compared to the equal-probability spatial design in either Newfoundland or Saskatchewan and small statistical differences in mean spatial balance in Yukon Territory ( Table 3). The BOSS design showed small but statistically significant differences in spatial balance relative to a purely spatial design in both Newfoundland and Yukon Territory (Fig 8, Table 3), but not Saskatchewan ( Table 3). The habitat stratification design had the largest overall differences in spatial balance relative to a pure spatial design, with small but statistically significant differences in all three jurisdictions and the largest difference occurring in Yukon Territory (Fig 8, Table 3).

Optimal design selection
The BOSS design provided the most optimal combination of cost savings and both spatial and habitat representation. Regardless of jurisdiction, the BOSS design had costs competitive with the cost design (Fig 6 and 9), with no loss of habitat or spatial representation. For all three jurisdictions combined, selecting the draw that simultaneously maximized representation while minimizing cost (see filled squares in

Field implementation
Of the 295 PSUs sampled across all three jurisdictions, we sampled 85 PSUs in remote areas that we predicted would require charter flights to access. In Newfoundland, cost models accurately predicted an average cost of $1,721 (SD = $86) per PSU to access 16 PSUs, compared to actual costs which averaged $1,721 (SD = $441) per PSU (Fig 8). In comparison, the cost models for Saskatchewan predicted an average access cost of $3,741 (SD= $1,928) for the 43 remote PSUs sampled, but the actual access costs were substantially lower (mean = $1,573, SD = $1,562). In Yukon, the mean predicted access cost to reach the 26 remote PSUs sampled was $1,574 (SD = $82; Fig 10) per PSU, but the mean cost to access these PSUs was $3,016 (SD = $132). Ordinary least squares regression suggests across all three jurisdictions the access costs were on average 48% of the model predicted cost on average (β = 0.48, SE = 0.10; Fig 10).

Discussion
We have integrated key concepts from sampling theory to develop a hierarchically structured survey design that can provide stratified sampling of predefined strata (e.g. ecoregions) and achieve the desired spatial balance and habitat representation while minimizing access costs. The probability-based BOSS design allows for valid statistical inference from field samples across the sampling frame and reduces operational costs. We found that when compared with equal-probability spatial designs (spatial design) or habitat-stratified unequal-probability spatial designs (habitat stratification design), our integrated approach was significantly less expensive.
Predicted access costs were 1.3−1.7 times more expensive for the habitat stratification design and 1.6−2.7 times more expensive on average for the spatial design compared to using our hierarchical BOSS design. Additionally, habitat representation in the BOSS design was statistically indistinguishable from the habitat stratification design for two out of three jurisdictions and was only marginally different for Newfoundland. While we did find minor statistical differences in spatial balance between the designs considered, spatial balance was always close to a value of one indicating the samples were well spread over jurisdictions irrespective of the type of design. In addition, the distribution of spatial balance metrics across iterations did not overlap one across all jurisdictions for the spatial design, and thus it does not appear that the inclusion of unequal sampling probabilities introduced any systematic biases in spatial representation. Combined with our optimization, we therefore simultaneously achieved habitat and spatial representation goals while minimizing access costs.
Our pilot field seasons suggest that logistic considerations might result in access costs differing from those suggested by our cost models. While cost models for Newfoundland have thus far provided accurate prediction of average access costs per PSU, our cost models slightly underestimated access costs in Yukon. Overall, regression analysis suggested that field implementation was 52% cheaper than suggested by our cost models. Many reasons likely contributed to the cost differences we observed, including external logistical support. One key difference is that our cost models largely portray cost of access to a sample unit as independent of access to other nearby sampling units. In practice, we were able to combine logistics and access methods for many sampling units and thereby substantially lowered the average cost per sample unit. For example, in the summer of 2019 we chartered a Twin Otter to drop off two crews on or near opposite sides of Cree Lake, Saskatchewan (57°23′57″N 106°40′10″W) to access six PSUs via canoe. When the plane returned to retrieve crews, an additional crew was dropped off to canoe down a river back to a roadside pickup location, allowing access to two additional PSUs. We were therefore able to access eight PSUs for the cost of two charter flights, substantially reducing the average cost per PSU. We therefore recommend that programs revise their cost models as field implementation improves knowledge of local logistics. As a result, final program costs will likely be lower than initially predicted by preliminary cost models, Our sampling design has several key differences from other large-scale terrestrial bird monitoring programs that improve its' utility for remote areas such as the boreal forest. Other programs such as the Integrated Monitoring in Bird Conservation Regions (IMBCR) program [9], the UK Breeding Bird Survey [56], the North American BBS [12] and the US National Park Service Vital Signs Monitoring program [57] also have defined sampling frames and randomly selected sample units that facilitate valid statistical inference from field samples [58]. Importantly, our design has fewer restrictions on the sampling frame while minimizing sampling costs, and implicitly accommodates health and safety considerations. While other programs restrict the sampling frame to roadsides or trail networks to facilitate access and alleviate safety concerns, we instead dealt with access by incorporating access costs into sample unit selection, and designed primary and secondary sampling units to have first responders available if needed. Thus, our design increases safety and reduces costs without limiting statistical inference about bird populations to accessible areas [58,59,60]. In addition, our hierarchical stratification with static strata (BCRs, jurisdictions and ecoregions) allows for assessment of the population-level responses of birds to natural and human disturbance and climate change [9,61]. Indeed, Pavlacky et al. [9] recently demonstrated the ability of their hierarchical spatial design to allow estimation of both strata-specific population sizes and species-habitat relationships. Our approach of integrating both cost and weighted habitat-inclusion probabilities within a similar spatially balanced design will allow us to draw similar inferences to Pavlacky et al. [9] at multiple spatial scales, but at reduced costs compared to alternative approaches. The ability to model species-habitat relationships is a key feature of the BOSS design, which will greatly improve our capacity to develop and test models of species response to habitat, which are crucial to predicting responses to industrial development [15,16,62] and climate change [18].
Our design directly incorporated several factors into the stratification and layout of primary and secondary sampling units to address the need to minimize costs and provide for the health and safety of field staff. Field implementation of the optimal designs selected here during pilot field seasons across all three jurisdictions showed that the BOSS design was achievable on the ground. We were able to apply the design across a wide variety of challenging conditions from mountainous terrain (e.g. >2000 m), wetlands (floating bogs, fens and marshes), forests with blow-down (recent burns, tornados and insect outbreaks), as well as steep, rocky, and densely vegetated slopes across all three jurisdictions. Furthermore, drawing an oversample using the GRTS algorithm allowed us to replace unsafe PSUs without altering the representativeness of the sample [25]. Using pairs of observers to conduct surveys within each PSU means that staff can work alone to each survey a plot of nine SSU within a PSU or they can work together to complete one plot of nine SSU before travelling together to complete the second plot. The latter option may be appropriate if a PSU contains dangerous terrain or wildlife or limits communication between staff (radio, satellite transmitter).
One complication imposed by our design is that most analyses assume that sampling is proportional to habitat availability and, by extension, population sizes. We used habitat stratification to ensure representation across the range of available habitat classes and thereby provide reasonable sample sizes to facilitate modeling abundance and distributions of rare species often associated with rare habitats. Habitat stratification should also improve our future capacity to detect trends of species using rare habits. As a result, extrapolating from the sample to estimate stratum-specific population sizes and/or population trends is more complicated because strata boundaries and sizes can change with shifts in vegetation through time [10]. This problem can be addressed analytically via post-hoc stratification [63,64], weighting via the inclusion probabilities [65,66], including random-effects for the strata indexed by the inclusion probabilities [67], or inclusion of interval-specific covariates reflecting habitat supply [66].
While our BOSS design presented here incorporates several key improvements over current large-scale bird monitoring designs, we foresee several areas where further continued improvements could occur. First, we stratified habitat using coarse satellite-derived products due to their availability and consistency across the entire sampling frame. Incorporating higher-resolution products capable of further stratifying habitat classes, and potentially better quantifying the amount and distribution of rare habitats, should improve our capacity to monitor species associated with these habitats.
One possible approach worth investigation would be to create an avian habitat classification (sensu [42]) from higher-resolution data, such as forest inventory data, for the portion of the sample frame for which they may be available. Higher-resolution data could be integrated into the current framework by (1) using the coarse-resolution satellite data to select PSUs and using the high-resolution habitat classification to stratify SSUs within the selected PSUs, or (2) splitting the selection of PSUs between areas with versus without access to these high-resolution data products.
Second, once sufficient on-the-ground data are available, it would be desirable to re-examine the stratification by using the design based sample to directly estimate spatial and temporal variance in bird community composition for Neyman allocation rather than the proxy variables used in the initial sample. Regional differences in timing of implementation, available resources, and historic sampling efforts mean that this data will be available earlier in some regions than others. We anticipate that sampling efficiency will improve by incorporating these direct estimates of variability into the design, but overall costs could be inflated if these data are not incorporated until they are available for all regions. Using data from preliminary rounds of sampling under the sampling design to validate Neyman allocation based on multivariate dispersion of proxy variables will ensure allocation between strata is cost effective.
To avoid inefficiencies, there is a clear need for improvements in cost models or habitat layers to be reflected in revised sample draws; however, revised sample draws should account for existing sampling under the BOSS design to avoid inducing spatial imbalance. Foster et al. [26] recently proposed using a squared loss distance metric to alter the inclusion probabilities around "legacy" monitoring sites (i.e. pre-existing sampling locations are known to be a randomly selected representative sample from a known sample frame [26]). This approach would provide a clear way to improve future sampling under the BOSS design without introducing biases and inefficiencies due to spatial imbalances. In addition to data from legacy sites, a large collection of historic data exist within the Boreal Avian Modeling (hereafter BAM) project database [68] for a large portion of boreal Canada. The difficulty with incorporating data from the BAM database is that these largely represent "iconic sites" (sensu [26]) which are either known to be non-randomly selected and/or there is insufficient documentation to know whether the sampling frame included the full suite of habitat classes. Including iconic sites in a randomized draw could introduce bias into the design and result in biased status and trend estimates that may not match those in the population [26]. Since these historic datasets could potentially add valuable ecological insights into factors affecting populations of boreal birds, it would be fruitful to consider when and how they could be included within a long-term monitoring design.
We have demonstrated that including access costs and habitat within a hierarchically structured sampling scheme was the most cost-efficient design that maintained the desired spatial and habitat representation. treated as legacy sites) into long-term monitoring to provide more responsive conservation efforts in light of widespread declines in biodiversity [1][2][3]. Our hierarchical sampling design should be widely implemented to monitor boreal birds because it can provide an unbiased representation of when and where conservation and management should be targeted [9,70]. Using monitoring data to inform conservation decisions is crucial, given ongoing and projected ecological changes in the boreal biome [18,34] and population declines already observed for boreal birds [3].