Comparing methods to estimate the proportion of turbine-induced bird and bat mortality in the search area under a road and pad search protocol

Estimating bird and bat mortality at wind facilities typically involves searching for carcasses on the ground near turbines. Some fraction of carcasses inevitably lie outside the search plots, and accurate mortality estimation requires accounting for those carcasses using models to extrapolate from searched to unsearched areas. Such models should account for variation in carcass density with distance, and ideally also for variation with direction (anisotropy). We compare five methods of accounting for carcasses that land outside the searched area (ratio, weighted distribution, non-parametric, and two generalized linear models (glm)) by simulating spatial arrival patterns and the detection process to mimic observations which result from surveying only, or primarily, roads and pads (R&P) and applying the five methods. Simulations vary R&P configurations, spatial carcass distributions (isotropic and anisotropic), and per turbine fatality rates. Our results suggest that the ratio method is less accurate with higher variation relative to the other four methods which all perform similarly under isotropy. All methods were biased under anisotropy; however, including direction covariates in the glm method substantially reduced bias. In addition to comparing methods of accounting for unsearched areas, we suggest a semiparametric bootstrap to produce confidence-based bounds for the proportion of carcasses that land in the searched area.


Introduction
Impacts to birds (particularly raptors) from collisions with wind turbines have been of concern in the U.S. since the early 1990s (Orloff and Flannery 1992) and impacts to bats since the early-2000s (Kerns and Kerlinger 2004). The magnitude of the impact is typically assessed through post-construction monitoring (PCM) (Howell and DiDonato 1991;Kunz et al. 2007;Arnett et al. 2008;Strickland et al. 2011). PCM to estimate mortality-the total number of turbine-induced bird or bat fatalities occurring over a specified period of time at individual wind projects-typically involves conducting carcass surveys, where human observers (or well-trained human-dog teams) search for carcasses within plots centered at randomly selected turbines (US Fish and Wildlife Service 2012). Under ideal circumstances, all carcasses from the sampled turbines land and remain in the designated search plots, and the entire area within each plot is searchable. In addition, under ideal conditions, all carcasses are detected by searchers with probability = 1.
In practice, search conditions tend to be far from this ideal and counts of observed turbine-caused wildlife fatalities will necessarily underestimate total mortality (Huso 2011), because not all turbine-caused fatalities are observed in the search process and each carcass may have a unique probability of detection depending on when or where it falls and its specific characteristics (e.g., size or color). Carcasses can be removed by scavengers very soon after they arrive (Bernardino et al. 2011;Bispo et al. 2013;Flint et al. 2010), and humans and even the best scent-dogs fail to detect some carcasses (Homan et al. 2001;Arnett 2006;Mathews et al. 2013;Beebe et al. 2016;Reyes et al. 2016). Trials designed to estimate the probability that a carcass persists to the next search must be flexible enough to allow for non-constant removal and must account for potential differences due to carcass size, season, vegetation density in which the carcass lands, and other carcass-specific factors (Bernardino et al. 2011;Bispo et al. 2013;Turner et al. 2017). Likewise, searcher efficiency-the ability of a searcher to find a carcass that is present-may change in time (Wolpert 2013). The number of trial carcasses necessary to adequately estimate the effects of each of these factor combinations can be expensive. For example, placing 50 trial carcasses, as recommended by Strickland et al. (2011), in each of 4 size classes, 4 visibility/vegetation classes, and 4 seasons would result in 3200 total carcasses during a year and could potentially attract scavengers to the site.
The maximum distance at which a carcass can land from the turbine is unknown (Hull and Muir 2010), but empirical evidence suggests that species as small as bats can land ≥ 90 m from the turbine (Erickson et al. 2003;Chatfield et al. 2009) (Fig. 1). Extending the most commonly implemented plot shapes, square or circular (US Fish and Wildlife Service 2012), to or past 90 m is expensive and inefficient, because area increases with the square of the distance searched, but even with constant numbers within each distance class, carcass density (carcasses per unit area) decreases (Fig. 1). Observed (dots) carcass density (number of carcasses per meter 2 ) for each 10-m annulus and estimated carcass density (dashed line). Poisson regression was used to obtain the fitted curve assuming the log mean carcass density follows a quadratic function of distance Furthermore, even if the designated search plot encompasses the entire population of carcasses, thick vegetation, crops, dangerous terrain, and rivers or lakes within the plot boundaries can make surveying the entire plot impractical. Accurate estimation of mortality requires accounting for unsearched areas by estimating the proportion of total carcasses that landed within the searched area (Jain 2005;Huso and Dalthorp 2014;Huso et al. 2017), not the proportion of the designated plot area that was searched. A method to accurately estimate the probability density of carcasses at any distance and direction from the turbine would enable accurate estimation of the proportion of carcasses falling inside the searched area of plots of varying shapes and sizes, accommodating the site-specific conditions that might prevent full searches of standard plots.

Meters
Estimators of mortality, M, the total number of turbine-caused fatalities at a wind project over a specified period of time, have primarily been Horvitz-Thompson type estimators where observed carcass counts are divided by an estimated probability of detecting a carcass (Shoenfeld 2004;Huso 2011;Huso et al. 2012;Korner-Nievergelt et al. 2013;Wolpert 2015;Dalthorp et al. 2018), but see Etterson (2013) and  for alternative approaches. Although detection probability must account for all the major sources of imperfect detection outlined above, in this paper we focus on estimating one aspect: the proportion of carcasses landing within the searched area, a. Because the magnitude and variance of a can greatly influence the precision with which M is estimated, in this study we look at accuracy and precision of estimating both a and M. To isolate the influence of the magnitude and uncertainty ofâ we evaluate its properties under idealized conditions where other detection factors are assumed known, but realistic in their magnitudes.
We consider two primary factors that could influence estimates of a: (i) the mortality rate at each turbine, which in turn influences the number of carcasses discovered with which to fit models of a; and (ii) the spatial distribution of carcasses around each turbine relative to the searched area that in turn can influence the number of carcasses discovered with which to fit models of a, particularly if prevailing winds or other factors cause anisotropic spatial distributions.
In practice, resources available to conduct PCM are limited, and trade-offs must be made between the size of the search plot beneath each turbine and the fraction of turbines that can be visited : the larger the search plot, the fewer turbines that can be sampled for a fixed monitoring budget. In addition, larger search plots will often comprise several visibility/vegetation classes with different probabilities of carcass persistence or searcher efficiency within each, incurring additional expense to conduct trials to estimate these probabilities within each class. In this context, it would be logical to consider limiting searching to areas in a single visibility/vegetation class (where carcass density or probability of detection is highest and consequently, variance of estimated mortality within the searched area is likely lowest) and rely on models of the proportion of carcasses landing within the searched area to extrapolate to the unsearched areas at the site. Such an approach would enable researchers to design carcass surveys that deviate from traditional square or circular plots and are conducive to increasing the maximum distance searched from the turbine without sacrificing efficiency. Road and pad (R&P) surveys represent a convenient form of implementation of this methodology.
In principle, an R&P search protocol has many potential advantages. The critical question is how reliably total mortality M can be estimated from data primarily (or entirely) collected on roads and pads. The general approach is to first estimate the number of carcasses arriving on R&P and then divide by the estimated proportion of carcasses falling on R&P,â. Ratio estimation (Thompson 2012) has also been utilized to obtainâ (Good et al. 2011) by clearing vegetation from a small number of square or circular plots out to a specified distance, exhaustively searching the entire cleared areas, and calculating the proportion of the total carcasses found on R&P. A potentially more robust and less costly alternative might be to model the distribution of carcasses as a function of distance (and possibly direction) from turbines and integrate over the R&P. Such methods have been implemented by Arnett et al. (2009), Huso and Dalthorp (2014), and Hallingstad et al. (2018. In this paper we compare the performance of five different estimators of a under a range of conditions that might be found in practice. For example, the configurations of roads might be along the cardinal directions in flat agricultural areas, but along ridgelines in more rugged terrain. The reported average number of carcasses per turbine can vary considerably among facilities, which affects the number of carcasses observed on R&P. Spatial carcass distributions can vary in unpredictable ways or may be highly dependent on direction of predominant winds. Comparisons focus on both potential bias and precision of each method for estimating a, and M from it, under varying conditions. In Sect. 2, we describe five methods for estimating the proportion of carcasses that fell within the searched areas: four methods developed from concepts proposed by Jain (2005), Good et al. (2011), Arnett et al. (2009 and Huso and Dalthorp (2014), and one novel method using a parametric weighted distribution. In Sect. 3 we describe our simulation of spatial distributions of carcasses and the search process under a range of conditions. In Sect. 4 we evaluate the statistical properties of these methods. In Sect. 5 we propose a semi-parametric bootstrap method to capture the uncertainty associated with estimating a, the proportion of carcasses falling within searched area boundaries. We conclude with a general discussion of the performance of the methodologies and potential limitations of estimation from R&P surveys in Sect. 6.

Methods of accounting for unsearched area
After colliding with a moving turbine blade, birds and bats fall to the ground in variable patterns that depend on several factors, including location of the strike along the blade, speed of the blade, wind speed, body weight, wing span, and possibly even animal movement after getting struck (for example, injured animals flapping in the air or crawling to relative safety after hitting the ground). The result is a spatial distribution of carcass density that is not constant but generally decreases with increasing distance from the turbine (Fig. 1). Because density is not constant, models that account for non-constant density must be employed to estimate the proportion of carcasses, a, that have landed in searched areas of variable configuration. Conditioned on a given carcass realization, a is a parameter. Several methods have been proposed to estimate a (Jain 2005;Huso and Dalthorp 2014;Good et al. 2011;Studyvin et al. 2019). A fourth method used by, for example, Gruver et al. (2009) and Drake et al. (2012) calculates a as the simple proportion of area searched withing the designated plot. This assumes that the designated plot comprises all areas within which a carcass could fall and the carcass density is independent of distance, "an implausible assumption when turbineinduced mortality occurs" (Huso and Dalthorp 2014). We compared the first three methods, but the method used by Studyvin et al. (2019) had not been published when we completed our analysis and is not included.
The methods described in this section assume data are derived from carcass surveys conducted on plots that extend at most r max m from the center of the turbines selected for surveying. Let i in {1, . . . , C} index each observed carcass with C as the total number of observed carcasses. Each method can be implemented by recording the location of each carcass as a distance, r i , and direction, θ i in [0, 2π), with respect to the nearest turbine. Note that we assume no carcass can land within 2 m of the center of any turbine due to the girth of the tower; therefore, r i is in [2, r max ].
In the description of the methods in this section, we assume that the searched area has been mapped out in a GIS software, so that information such as the proportion or amount of area searched at a specified distance or location relative to nearest turbines is known.
In the description of some methods below, and in the results presented in Sect. 4, we assume the conditional probability, G, that a carcass is observed given it landed on R&P is known. In addition, for the ratio method, we assume the probability, G clr , that a carcass is observed given it landed on non-R&P at a turbine cleared of vegetation and surveyed is also known. In practice, these quantities are unknown, and numerous methods have been proposed for their estimation (Shoenfeld 2004;Huso 2011;Korner-Nievergelt et al. 2015;Wolpert 2015;. Including uncertainty in mortality estimates from the probability of detection is beyond the scope of this research. For the comparison of all methods, except the ratio method, it is acceptable to assume G is known, because a is estimated independently of G, so uncertainty from estimating G would have a similar effect on these methods. In the case of the ratio method, the detection probabilities in the two visibility classes are used to estimate a, so assuming these values are known results in more precise estimates of a than would be observed in practice; however, as discussed in Sect. 4, the performance of this method is generally poor, so this slight advantage does not obscure the conclusions of this research. In this section we introduce five methods for obtaining a point estimate for a. In Sect. 5 we propose a general method that can be used to calculate a standard error of a or confidence interval for a using any of the methods described below. Jain (2005) proposed accounting for unsearched areas by partitioning the area beneath turbines into 10 m concentric annuli, or rings, and calculating mortality independently for each ring. This method intrinsically assumes that carcass density is homogeneous within an annulus, which under isotropic carcass spatial distributions, becomes more reasonable as the buffer-distance between the outer and inner radii of a ring-goes to zero. In practice, the lower bound of the buffer will be determined by the resolution of the equipment used to record carcass distance with respect to the center of the nearest turbine. The Cake method refers to the methodology proposed by Jain (2005) with a buffer < 10 m. We introduce the notation of the Cake method in this section, and examine the relative performance of this methodology in the results presented in Sect. 4 utilizing a 1 m buffer.

Cake method
To implement the Cake method, let j in {2, . . . , (r max − 1)} index the rings, C j = C i=1 I(r i = j) be the number of observed carcasses j m from the center of the surveyed turbines, and A j be the known proportion of area searched j m from the surveyed turbines. For an individual ring, the total number of fatalities, M j , is estimated as,M This method directly estimates M with M cake = (r max −1) j=2M j . Given M cake , the estimate of a is calculated aŝ We refer to this as the "Cake" method, because annuli closest to the turbine tend to have the highest density of carcasses ( Fig. 1), and the assumption of homogeneity within an annulus creates a step function yielding a tiered wedding cake appearance when plotted in three dimensions. Good et al. (2011) used ratio estimation (Thompson 2012) at Fowler Ridge Wind Farm in Indiana, USA to account for unsearched areas. The R&P of 109 out of 355 turbines were searched out to 40 m, and in addition, a square 80 × 80 m plot around 9 of the turbines was cleared of vegetation and searched, albeit with a lower searcher efficiency than on the R&P. Implementation of the ratio method requires each carcass be assigned an indicator, I(clr ) i , to indicate if the carcass was found at a turbine cleared of vegetation (regardless of visibility class), and an indicator, I(R P) i , to indicate if the carcass was observed on R&P. This information can be obtained either using the location data, or within survey records.

Ratio method
To implement this method, define C R P = C i=1 I(clr ) i · I(R P) i as the number of observed carcasses on R&P at turbines cleared of vegetation, and C non R P = C i=1 I(clr ) i · 1 − I(R P) i as the number of observed carcasses on non-R&P at turbines cleared of vegetation; and recall that G and G clr are the known conditional probabilities a carcass is observed given it landed on R&P or non-R&P, respectively, at surveyed turbines. Implementing the ratio method, a is estimated aŝ (2) When applying the ratio method in the simulations described in Sect. 3, we deviated from the protocol used by Good et al. (2011), simulating estimates from a protocol where the cleared plots encompass the entire spatial carcass distribution (out to 120 m), and the R&P are searched out to 120 m at every turbine. Obviously this would produce more precise and accurate estimates relative to cleared plots that do not comprise the full spatial carcass distribution and R&P searches conducted on less than one third of the total number of turbines.

Weighted distribution (WD) method
The weighted distribution method uses maximum likelihood estimation to estimate the univariate carcass distance distribution by viewing the observed carcass distances, {r 1 , . . . , r c }, as a random sample from the observed carcass distribution, f o (r | ).
Assuming that carcass distance (both observed and unobserved) follows a parametric family distribution (e.g. Weibull, Gamma, Log-normal, etc.), indexed by a parameter vector , and the spatial carcass distribution is isotropic-independent of directionthe pdf of the observed carcass distribution is where A(r ) is a known function that yields the fraction of area searched r meters from the center of surveyed the turbines (Patil and Rao 1978). Heuristically, this assumes the probability a carcass is found at a given distance is proportional to the probability a carcass lands at that distance, times the proportion of that distance searched. The maximum likelihood estimate (mle) of ,ˆ mle , can be estimated by maximizing the log-likelihood function, which can be used to estimate a aŝ In practice, this method can accommodate any parametric family distribution, and AIC (Akaike 1998) can be used to compare between fitted models. In obtaining the results discussed Sect. 3, we assumed carcass distance can be modeled with a Weibull distribution, because it produced the lowest AIC of the parametric models considered for modeling the distance of the carcass distances (Fig. 2). We also assumed distance to the nearest turbine could be measured within a meter, and recorded r i as either 2.5, 3.5, . . . , or 119.5. We defined A(r ) as the proportion of area within 0.5 meters of r contained in R&P, which is analogous to A (r −0.5) in the notation introduced for the Cake method. After development, we shared this method with researchers who implemented the weighted distribution method in the context of estimating eagle fatalities at wind projects (Hallingstad et al. 2018).

Generalized linear model method
Originally proposed by Arnett et al. (2009), Huso and Dalthorp (2014) examined the statistical properties of accounting for unsearched areas by partitioning the area beneath surveyed turbines into a grid and modeling the spatial distribution of carcasses with a generalized linear model (glm). Huso and Dalthorp (2014) introduced their method using logistic regression. Warton and Shepherd (2010) demonstrated that a Poisson point process model produces similar coefficient estimates and standard errors as a logistic regression as long as the grid used to fit the models is sufficiently dense. In the simulations described in Sect. 3, we used a 1m dense grid and results for the logistic and Poisson regression models were virtually indistinguishable. The Poisson model is advantageous in this setting as it allows for more than one carcass to occupy a grid cell. It is worth noting that this method can also be established by partitioning the area beneath turbines into rings, similar to the Cake method. We chose to use grid notation to draw attention to the similarity of this method to the Anisotropic glm method introduced below.
The method to account for unsearched areas developed by Huso and Dalthorp (2014) in the context of a Poisson regression model fit on a grid where cells are indexed by distance, r in [2, r max ] m, and direction, θ in [0, 2π) radians, to the nearest turbine, assumes the number of observed carcasses in a grid cell, , follows a Poisson distribution. Its mean, λ r ,θ , modeled as a function of distance from the turbine, where O r ,θ is the offset for the cell indexed by (r , θ), and is equal to the number of turbines at which that cell is searched. The fitted Poisson regression model can be utilized to estimate the spatial carcass distribution,f (r , θ), aŝ The estimate of a is obtained by integratingf (r , θ) over the searched area. If an smeter dense grid is utilized, with S representing the collection of points indexed by (r , θ), then the integral can be calculated numerically aŝ where T is the total number of turbines at the wind facility.

Anisotropic generalized linear model method
The cake, weighted distribution, and glm methods all assume the spatial carcass distribution is isotropic. Under an anisotropic distribution, carcass density changes with direction, in addition to distance, with respect to the nearest turbine. The glm method can readily be modified to account for anisotropic carcass distributions by incorporating direction, θ , into the Poisson regression model. Thus, the Anisotropic glm method, like the glm method also partitions the area beneath surveyed turbines into a grid with cells indexed by r in [2, r max ] m and θ in [0, 2π) radians, and assumes that the number of observed carcasses for a specific grid cell, , follows a Poisson distribution; however, the anisotropic glm allows E(C r ,θ ) = λ r ,θ to differ by direction.
Direction has to be included into the model as a circular variable, which accounts for the cyclic pattern of direction in which, for example, 0 • and 360 • are equivalent and 5 • and 355 • are equidistant from 0 • . There are numerous techniques for incorporating circular variables into regression models (see Pewsey et al. (2013)). Perhaps the simplest approach for modeling a continuous response, Y , as a function of a circular variable, θ , is E(Y ) = β 0 +β 1 cos(θ )+β 2 sin(θ ). This model can be made more flexible by considering the K th Fourier expansion on θ , and including cosine and sine terms for each of kθ , k ∈ {1, . . . , K }. Accounting for an anisotropic spatial carcass distribution in the context of Poisson regression requires allowing the predicted carcass density for a given distance to differ by direction, which can be achieved using interactions between distance terms and direction terms (represented by sine and cosine terms), log(λ r ,θ ) = β 0 + β 1 r + β 2 r 2 + K k=1 β 4k−1 r cos(kθ) + β 4k r sin(kθ) + β 4k+1 r 2 cos(kθ) + β 4k+2 r 2 sin(kθ) The main effects for angle (i.e. sin(θ ) and cos(θ )) are not incorporated into (5), because it seems reasonable to assume a common intercept. In this application, no carcass can land within 2 m of the turbine, but theoretically the generalized linear model used to estimate the spatial carcass distribution produces an estimated carcass density at the origin. If the main effects for angle were included in the model, the estimated density at the origin could change depending on the value for θ that is used, providing multiple density estimates for a single location. To avoid this issue, the main effects are not included in (5).
The value of K represents the order of the Fourier expansion on θ . Higher values of K allow for more flexible models; however, increasing K by 1 requires estimation of four additional parameters. An excess of terms in the model could potentially increase the variation in mortality estimates. In practice, different values of K can be considered, along with nested models, such as those excluding sine terms or terms with r 2 , and the appropriate Poisson regression model can be obtained using any reasonable model selection method. To obtain the Poisson regression model for the anisotropic glm method in the simulation described in Sect. 3, we used forward model selection, increasing the value of K one unit at a time (up to K = 3), and using the likelihood ratio test with α = 0.05 to compare models.
The anisotropic glm method differs from the (isotropic) glm method only in the Poisson regression model used to obtain the estimated bivariate carcass distribution. After obtaining the fitted model, the procedure for estimating a is identical. The fitted values from anisotropic method, log(λ r ,θ ), can be used to model the anisotropic spatial carcass distribution (designated by superscript (a)) aŝ Analogous to (4), if an s-meter dense grid is utilized, with S representing the collection of points indexed by (r , θ), the estimate of a is calculated aŝ

Simulation study
We compared the performance of the five methods described above using simulations designed to (i) evaluate accuracy and precision of estimates for the proportion of carcasses that land in the searched area, i.e., on R&P; (ii) evaluate the influence of various conditions, e.g., per-turbine fatality rate, R&P configuration, spatial density patterns, etc. on the accuracy and precision of each method; and (iii) understand how variability inâ contributes to variation in mortality estimates,M. The simulations considered four spatial carcass distributions, two R&P configurations of actual wind projects, and four per-turbine fatality rates. One thousand trials were performed for every combination of settings. Each trial simulated a data set of observed carcasses, designed to mimic the random process that would be observed in practice.

Distance distributions
Our capacity to accurately estimate the fraction of carcasses that land within our search boundaries depends on whether the spatial distribution is isotropic and if not, whether our search boundaries comprise a disproportionate amount of area with either higher or lower carcass density. Actual spatial distributions of carcasses will likely depend on (i) wind speed at the time of collision that will affect distance the carcass travels downwind, (ii) location of the collision along the blade that will determine the relative speed at which the blade is moving when the carcass is struck, and (iii) position of the blade in its rotational arc that will determine whether the carcass is lifted as it is struck or forced directly downward, as well as turbine height, topography, scavenger behavior, and other factors. We simulated carcass locations using four spatial distributions from factorial combinations of isotropy/anisotropy and empirical/modeled based on 200 bat carcasses observed at a wind project in the Midwestern U.S. from April to November in 2013 (Figs. 1, 2, and 3). Searches were conducted daily on 15 circular plots cleared of vegetation out to 90 m from the turbine. Assuming constant probability of observing a carcass within each cleared plot, the distribution of observed carcasses provides an estimate for the true underlying spatial carcass distribution. Furthermore, searchers were able to determine that the 200 carcasses arrived the night prior to the search, enabling us to assign covariate information, such as average nightly wind direction, which was used to create an anisotropic spatial carcass distribution.
We refer to the two spatial carcass distributions that are independent of direction as the Isotropic Empirical distribution and the Isotropic Gamma distribution (Fig. 3a). The Isotropic Empirical distribution is the empirical distribution of the 200 observed carcass distances to the nearest meter, plus 10 simulated carcass distances. The simulated carcass distances were included because the empirical distribution (and corresponding parametric fits) indicated that carcasses likely landed past 90 m (Fig. 2). Therefore we used the triangular distribution on 90 m to 120 m, to simulate ten additional locations increasing the total by 5% to create a more realistic carcass distribution. The Isotropic Gamma parameters, shape = 2.4 and scale = 19.5, were calculated by fitting a Gamma distribution to the observed carcass distances, and then using this fit to extrapolate to 120 m. The Gamma distribution was utilized, because it produced the third lowest AIC of the parametric methods we considered to model carcass distance (Fig. 2). Although the Weibull distribution and scaled-Poisson regression model had lower AICs, we chose to model the spatial distribution with the Gamma to avoid favoring either the weighted distribution (3), or the glm (4), methods by generating and fitting data with the same model.
We refer to the two spatial carcass distributions that depend on direction as the Anisotropic Empirical distribution and the Anisotropic Gamma distribution. The Anisotropic Empirical distribution was created using the location of the 200 observed carcasses, maintaining their distances but rotating them by an amount equal to the average wind direction on the night they were killed minus 180 o . This simulates the location relative to a prevailing wind originating from the south (Fig. 3b). To these, we added the same 10 simulated distances used for the Isotropic Empirical distribution, and the angle was determined by sampling with replacement from the rotated angles that landed at least 45 m from the turbine. Simulated carcass location under the Anisotropic Empirical distribution was assigned by randomly sampling a location displayed in Fig. 3b and shifting the coordinates uniformly within a 10-m radius of the sampled location (Fig. 3c). The 10-m shift represents a compromise between reproducing the observed data too closely vs. not closely enough. The Anisotropic Gamma distribution (Fig. 3d) was created by distorting the spatial carcass distribution of the Isotropic Gamma distribution, to create a realistic anisotropic distribution, while maintaining the same univariate carcass distance distribution to its isotropic counterpart. The Anisotropic Gamma distribution simulates conditions where the prevailing winds originate in the east; carcasses are twice as likely to land west (downwind) of the turbine, where they are also more dispersed relative to carcasses that land east (upwind) of the turbine. There is a vast variety of possible patterns of anisotropy, and it is well beyond the scope of this study to analyze or characterize them all in any meaningful way. Instead, we chose two contrasting patterns that represent two somewhat extreme cases because the wind directions in each are so constant, as illustrated in the contrast between the observed carcass dispersion pattern (Fig. 1a) with variable winds and the anisotropic dispersion pattern derived from the same data but adjusted to mimic a constant prevailing wind direction (Fig. 3b).
To simulate carcass location for a specific spatial carcass distribution, we assigned a probability mass to each 1 m grid cell, and sampled grid cells proportionally to the probability a carcass lands in that cell. For more detail see the R code in Appendix 1.

R&P configurations
We selected two current wind projects to demonstrate how the methods perform under completely different R&P configurations. WP1 is located amid corn and soy farms in Indiana, has 355 turbines, and road orientation is predominately aligned with cardinal directions. WP2 is located along two mountain ridges in Pennsylvania and has 23 turbines. The roads at WP2 have more twist and turns with no pronounced unifying pattern like WP1. The R&P configurations for both facilities were obtained using Google Earth, mapped in ArcGIS, and imported into R R-statistical software (R Core Team 2015) as polygons (Fig. 4).

Average number of fatalities per turbine
In a given simulation trial, mortality was determined by the number of turbines at the wind project, 355 at WP1 and 23 at WP2, and the average number of carcasses per turbine: 2, 10, 25, or 70. We chose these mortality levels to represent the range of rates reported in the publicly available literature. A rate of 70 carcasses per turbine represents the average of two of the highest estimated mortality per turbine rates we are aware of. BHE Environmental (2011) (Table 1), and the probability that 0 carcasses are observed is less than 0.00015. Averages of 10 and 25 carcasses per turbine were selected to assess the performance of the methods when the rate is between 2 and 70.
While mortality was fixed for a set of conditions, the number of observed carcasses differed in each trial. The expected number of carcasses observed depends on spatial carcass distribution, R&P configuration, and average number of carcasses per turbine (Table 1).

Simulating carcass counts
Simulations were designed to mimic a setting where R&P carcass surveys were conducted at every turbine, to the maximum assignable distance from the turbine: 120 m. To simulate a set of observed carcasses used to estimate mortality, M simulated carcasses were assigned (i) a turbine number using the discrete uniform distribution on the number of turbines at the wind project, (ii) a carcass location using one of the four spatial carcass distributions, (iii) an indicator of whether the carcass landed on an R&P based on its location and the specific R&P configuration of its turbine, and (iv) an indicator of whether the carcass was observed based on a random draw from a Bernoulli with p = 0.85 (Good et al. 2011). In each trial M is a fixed value equal to the number of turbines times the mean fatality per turbine rate. In addition to R&P surveys, the ratio method also requires information from clearedplot surveys. Good et al. (2011) searched 9 out of 355 (2.5%) turbines cleared of vegetation within an 80 m × 80 m plot centered at the turbine to provide information for the ratio estimator. To reflect this method, we simulated carcass surveys on non-R&P areas within 120 m of 9 turbines for WP1, although, this represents a far greater area than was searched in practice. We are not aware of the ratio method having been applied at facilities with only 23 turbines. Clearing the same percent of turbines at WP2 as WP1 corresponds to clearing < 1 turbine; therefore, we arbitrarily chose to conduct simulated surveys on non-R&P at 3 out of 23 (13%) turbines. The turbines selected for clearing differed in each trial, and were selected using a simple random sample. Each carcass on non-R&P was assigned an indicator of whether it was observed based on a random draw from a Bernoulli with p = 0.32 as reported by Good et al. (2011).
Simulations were perfomed using R-statistical software. One thousand trials were generated for every combination of settings (2 R&P configurations × 4 distributions × 4 rates =32 combinations). We applied all five methods to the same data generated in every trial, or simulated data set, to obtain the estimated proportion,â, of carcasses that land on R&P for each combination of settings, except at WP2 where the anisotropic glm method was not applied when the rate was 2 due to the low number of simulated carcasses observed (Table 1).

Results
The objective of this research is to compare several methods used to estimate a; however, in practice, a is a nuisance parameter that must be estimated in order to estimate the true quantity of interest: mortality, M. Therefore, in addition to displaying and discussing each method's performance in estimating a (Figs. 5 and 6), we show how variability inâ contributes to variability in mortality estimates by considering mortality estimate,M where C was the number of observed carcasses on R&P and G is the known probability a carcass is observed given it landed within the searched area (Table 2). In practice G is unknown and must be estimated, in which case the estimatorM would involveĜ. In displaying the results ofâ for each method, we calculatedâ a for each estimate in every trial. We dividedâ by a, because a is a random variable that changes from trial to trial even under the same settings (e.g. underlying carcass distribution, R&P configurations, and fatality per turbine rate). We chose to plot the results on the log scale, because underestimating a leads to greater bias in mortality estimates relative to an equidistant overestimate of a.

Isotropic spatial carcass distributions
In general, each method performed at its relative best at both wind projects under isotropic spatial carcass distributions with high average fatality rates per turbine (Table 2, Figs. 5 and 6). The cake, weighted distribution, glm, and anisotropic glm methods tend to produce similar results for estimates of mortality, (M), with the anisotropic glm method having a slightly larger standard deviation in some cases. The maximum observed bias for these methods under isotropic distributions with an average fatality per turbine rate of 70 is 3.9%. The ratio method produces the most biased results and the largest standard deviation in almost all cases. Under the Isotropic Empirical distribution, when the average rate of carcasses per turbine is 70, the standard deviation of mortality estimates for the ratio method at WP1 and WP2 are 3.5 and 1.5 times larger, respectively, than the standard deviation of the other four methods. To prevent its large variation from obscuring results of other estimators, we do not include results from the ratio estimator in Figs. 5 or 6, but the results are given in Table 2.
As expected, as the rate of carcasses per turbine decreases, the number of observed carcasses decreases and the variation ofâ a increases for all methods. For the ratio method, both bias and the variation ofâ, and consequentlyM, increase as the average rate of carcasses per turbine decreases. For example, under the Isotropic Empirical distribution at WP1 the relative bias of the ratio method increases from 2.4% when the average rate is 70 to 15.5% when the rate is 2 (Table 2). In addition, the ratio method cannot produce a mortality estimate if 0 carcasses are observed on the R&P at cleared turbines. When the mortality rate was low (2/turbine) this occurred in about 10% of trials at WP1 and more than 33% of trials at WP2, for every spatial carcass distribution. Results for which more than 1/3 of the estimates are inestimable are reported as NR.

Anisotropic spatial carcass distributions
The effect that anisotropic spatial carcass distributions have on each of the methods differed by wind project and the two anisotropic distributions. Consider the results from the cake, weighted distribution, and isotropic glm methods. Again, these methods all performed similarly for both anisotropic distributions at both wind project R&P configurations (Table 2, Figs. 5 and 6). It is not surprising that these methods are affected similarly in each setting, because none of the methods take direction into consideration in addition to distance, and therefore intrinsically assume the spatial carcass distribution is isotropic. The results from the simulations suggest that when this assumption is violated the effect depends on the spatial carcass distribution, and R&P configurations. Under the Anisotropic Empirical distribution, these methods overestimate a and consequently underestimate mortality at WP1 and WP2; while under the Anisotropic Gamma distribution, these methods overestimate mortality at WP1 and underestimate mortality at WP2 (Figs. 5 and 6).
The direction of the bias when estimating mortality using a method that assumes isotropy can be explained by examining the anisotropic distribution and a composition plot of the R&P configurations at the wind project of interest (Fig. 7). A composition plot displays the proportion of turbines' R&P that contain a specific point. Points Table 2 Mean (and standard deviation) for the estimated average number of fatalities per turbine, calculated asM divided by the number of turbines (355 at WP1 and 23 at WP2). The anisotropic method was not applied (NA) at WP2 when the rate was two, because of the low number of expected carcasses. The results for the ratio method at WP2 when the rate was two were not reported (NR), because this method failed to produce an estimate in more than 33% of trials The anisotropic method was not applied (NA) at WP2 when the rate was two, because of the low number of expected carcasses. The results for the ratio method at WP2 when the rate was two were not reported (NR), because this method failed to produce an estimate in more than 33% of trials a Based on the 90% of trials where an estimate was obtained b Calculated excluding one or two trials that are darker are contained in the R&P of many turbines, and points that are lighter are contained in the R&P on few turbines. From this plot it becomes apparent that most of the roads at WP1 are aligned with cardinal directions. Furthermore, most roads approach the turbine from the west, and fewest from the north. Now consider the two anisotropic distributions displayed in Fig. 3. Recall, under the Anisotropic Gamma carcass distribution, most carcasses land west of the turbine; and given the R&P configuration at WP1, the most data for any distance greater than 10 m from the turbine is collected west of the turbine. This causes the isotropic methods to overestimate the number of carcasses at most distances, producing an overestimate of mortality. Under the Anisotropic Empirical distribution relatively few carcasses land west of the turbine producing an underestimate of mortality. Using the same logic and Fig. 7(right), the direction of the bias for mortality estimates at WP2 using an isotropic method can be explained. Three trials were excluded when calculating the results for the anisotropic glm method at WP2: two under the Anisotropic Empirical distribution when the rate was 10, and one under the Anisotropic Gamma distribution when the average rate was 25. In each of these trials, the forward selection method identified the Poisson regression model with the third order Fourier expansion of θ as the most appropriate model, and the estimated mortality exceeded 4000 carcasses per turbine. These unreasonable estimates occur when the quadratic term for distance, which is the sum of seven terms for a fixed direction (see Eq. 5), is positive, causing predicted carcass density to increase with distance from the turbine. Researchers can easily identify this situation in practice with a graphical display, and use a lower order Fourier expansion or avoid this anisotropic method altogether; therefore, it is reasonable to exclude these trials when reporting the results. Attention to these trials is included here, so researchers are aware of a potential issue that may arise when including quadratic terms in the generalized linear model. The Anisotropic glm method had success in reducing the bias in mortality estimates under anisotropic spatial carcass distributions. When compared to the isotropic glm method, under the Anisotropic Empirical distribution the bias inâ was reduced by 52% and 75% at WP1 for average rates of 2 and 70 respectively, with negligible increases in the standard deviation. Under the same distribution at WP2 the bias ofâ was reduced by 73% and 94% for average rates 10 and 70 respectively; however, there was a 31% and 16% increase in standard deviation.
Recall, the order of the Fourier expansion, K , was determined using forward selection by increasing K one unit at a time and using the likelihood ratio test with α = 0.05 to compare models. Under the isotropic models, the coefficients for terms with cos and sin are all zero, and expectedly about 95% of the trials select the Poisson regression model with K = 0, which is the model used in the isotropic glm method (Table 3). In general, under the anisotropic distributions, the anisotropic glm method routinely identified Poisson regression models with K > 0.
The bias of the ratio method was relatively unaffected by the anisotropy at WP1 (Table 2). At WP2, there was an increase in the standard deviation of the ratio method under the anisotropic distributions. This is likely attributed to an increase in variation in the proportion of carcasses landing on the R&P at individual turbines. The ratio method assumes the proportion is similar among turbines. Under anisotropy, this assumption Table 3 Proportion of trials in which the anisotropic glm method used the k th order Fourier expansion for estimation using a forward model selection process based on p-values (α = 0.05) from the likelihood ratio test becomes less tenable at wind projects with large variation in R&P configurations at individual turbines.

Accounting for the uncertainty ofâ
When a is unknown and estimated using the methods described above, uncertainty induced byâ must be taken into account when calculating a confidence interval for M. In this context, a confidence interval for a and M can be calculated using a semiparametric bootstrap (Efron and Gong 1983); however in verifying this, we discovered using a standard bootstrap-drawing C locations with replacement from the distribution of observed locations and applying one of (1) through (5) to simulate estimates of a-in conjunction with the percentile method yields intervals below the desired coverage. This occurs because such a method does not account for binomial uncertainty in a.
Recall, a becomes a parameter after conditioning on a specific carcass realization that is generated from an underlying point process. Under this paradigm, F, the number of carcasses that land on R&P, follows a Binomial(M, P a ) distribution, where P a is the probability a carcass lands on R&P for the underlying point process specific to a wind project during the period when searches are conducted. The methods discussed in the previous sections actually provide direct estimates for P a , therein providing a reasonable point estimate for a = F/M, because E(a) = P a . Both the uncertainty in estimating P a and binomial uncertainty in F must be accounted for in calculating a confidence interval for a, which can be accomplished using a semi-parametric bootstrap to simulate the distribution of (â − a).
A single bootstrap is obtained by simulating F boot from a Binomial (M,â) distribution. Given F boot , a simulated number of observed carcasses, C boot , is obtained by drawing from a Binomial(F boot ,G) distribution. Observed distances are simulated by drawing with replacement C boot locations from the set of locations at which fatalities were observed, andâ boot can be obtained by applying the method used to estimate a to the simulated locations. A (1 − α) confidence interval for a is created by adding the α/2 and 1 − α/2 quantile from the distribution of (â boot − a boot ) toâ, where a boot = F boot /M .
As previously discussed, G is unknown in practice and estimated withĜ. Under such circumstances a confidence interval for a can be calculated by replacing G witĥ G in the above description.
Confidence intervals for M can be calculated by simulating a bootstrap distribution forM. For unknown G, this requires generating a bootstrap distribution forĜ, which can be achieved independently of the bootstrap distribution ofâ (Huso et al. 2012). Given C boot ,â boot , andĜ boot , a bootstrap mortatlity estimate is calculated as,

Discussion
The primary objective of this research was to compare five approaches to account for unsearched areas when estimating mortality at wind power facilities. A convenient example comes from defining the searched area as roads leading up to turbines, and pads beneath turbines (R&P). We compared the methods via a simulation study that considered different R&P configurations, spatial carcass distributions, and average carcass per turbine rates. We also proposed a semi-parametric bootstrap method that can be used to estimate uncertainty in estimates of a when using the weighted distribution, cake, glm, or anisotropic glm methods.
In nearly all settings we considered, at least one, if not all, of the weighted distribution, cake, glm, and anisotropic glm methods outperformed the ratio method in terms of accuracy and precision in estimating M. The ratio method displayed some degree of bias in every setting, which was expected because ratio estimation is always biased (Lohr 2009). In the context of turbine-induced bird and bat mortality estimation, the bias becomes increasingly prevalent as the average fatality rate of carcasses per turbine decreases.
In theory, the bias and variance of the ratio method could be reduced by clearing more turbines or increasing the probability of detection in the cleared area; however, exploratory work suggested that an unreasonable number of turbines (more than 25%) would have to be cleared to be competitive with the other methods we considered, and increasing the probably of detection to that of the R&P only produced modest increases in performance. It is interesting to note that in the exploratory research our results suggested that simply treating the data at cleared turbines as a random sample for mortality estimation produced much better results relative to using the data from the remaining R&P in ratio estimation.
In addition to its inferior statistical properties relative to the other methods we considered, the ratio method is also more expensive and not always possible to implement. Under the protocol used in the simulations, all methods require searching the R&P of every turbine out to 120 m. The weighted distribution, cake, glm, and anisotropic glm methods also require that the facility be mapped in a GIS program before estimation. The ratio method, though, requires additional data from some plots that are cleared of vegetation. When the vegetation is dense, terrain is treacherous, or if bodies of water are within the carcass fall zone, obtaining counts for cleared turbines is impractical or even impossible. Avoiding turbines with these conditions may result in a non-representative sample of turbines. When clearing beneath turbines is possible, the cost of machine clearing as well as compensation to the land owner for lost crops add to the overall cost of the monitoring effort. In non-agricultural settings, clearing will damage the surrounding habitat, increasing the ecological impact of the wind facility. The cleared areas may also represent a substantial increase to the total area searched. At WP1 the area comprising the R&P of all 355 turbines is equivalent to 13.1 turbines cleared to 120 m. Therefore, under the protocol we simulated where 9 turbines were cleared and searched out to 120 m, the ratio method requires searching 66% more area than is required for the other 4 methods. At WP2, searching the R&P of all 23 turbines is equivalent to searching 1.8 turbines cleared to 120 m. Therefore, under the protocol we simulated where 3 turbines were cleared and searched out to 120 m, the ratio method requires searching 156% more area than is required for the other 4 methods.
The weighted distribution, cake, and glm methods all performed similarly in every simulation, raising the question of which method to implement in practice if it is reasonable to assume the spatial carcass distribution is isotropic. Each method possesses its own advantages.
The weighted distribution method yields a parametric fit obtained by maximizing a likelihood function, which enables comparison of parametric models, using AIC, for example. This method completely specifies the distribution of carcass distance, and allows convenient estimation of parameters such as the mean, standard deviation, relevant quantiles, and the proportion of carcasses that landed beyond the maximum distance searched. For example, when modeling the data used to generate the underlying spatial carcass distributions in the simulations along with a Weibull distribution, it is estimated that 4.7% of carcasses landed past the maximum distance searched (90 m). Using a Gamma distribution, this estimate increases to 9.1%. Given the ease of extending the maximum distance from the turbine searched along access roads, it is highly preferred that researchers search at distances beyond the expected maximum fall distance of a carcass, if possible, instead of relying on extrapolation beyond the range of the data. If, however, it becomes impractical to extend the plots far enough, due to resources or terrain, the parametric methods (including the weighted distribution and the glm) can be used to estimate the proportion of fatalities falling beyond the maximum distance searched, whereas such extrapolation is not possible with the cake or ratio methods.
The cake method estimates mortality non-parametrically allowing for relatively unbiased estimates under a variety of isotropic spatial distributions when searches are conducted to or past the maximum fall distance of a carcass. It becomes biased, however, if searches do not extend far enough, or the carcass distribution is anisotropic. The cake method proposed in this research partitions the fall zone into concentric annuli, which is reasonable under isotropic settings. It may be possible to create an anisotropic version of the cake method using the same non-parametric approach, by partitioning the fall zone into elements that are internally homogeneous with respect to the anisotropic spatial carcass distribution (see Maurer (2017)); however, determining an appropriate partition requires information about the spatial carcass distribution, which is unknown.
The glm model is advantageous because, unlike the weighted distribution and cake methods, which only use the univariate carcass distances, it fits a model on the observed spatial carcass distribution. The model can easily be modified to form the anisotropic glm by incorporating additional parameters that reflect dependence of carcass density on direction as well as distance. In addition to reducing bias under anisotropic spatial carcass distributions, a statistical hypothesis test can be used to determine if anisotropy is detectable by testing if terms including direction are statistically significant to the model. If there is insufficient evidence of anisotropy, then the isotropic glm method can be selected as the most appropriate model, and results will be similar to the weighted distribution and cake methods.
As discussed in Sect. 1, R&P searches provide a convenient example for implementing and comparing the statistical methods discussed above. In addition, there are some logical advantages to R&P searches relative to full-plot searches conducted by human searchers: (i) the probability of detection tends to be maximized on R&P (Kerns et al. 2005;Good et al. 2012); (ii) searcher efficiency and carcass persistence trials only need be conducted in a single visibility class; (iii) extending the maximum distance searches increases required resources at a linear rate, if at all; and (iv) Searches can be conducted at a larger number of turbines, yielding information on spatial patterns across the entire project.
It was beyond our scope to compare these methods to a protocol that searches the entirety of a randomly selected subset of turbines when holding the total searched area constant; however, we did conduct some exploratory research. Our results suggested that under isotropic carcass distributions, estimation based on R&P was more precise relative to searches conducted on an equivalent amount of area at cleared turbines, when the probability of detection was lower in cleared areas than on R&P. As the probability of detection in the cleared areas approached the probability of detection on R&P, the performance became similar. This suggests that it would behoove researchers to focus on areas with high detection probabilities before searching areas with lower detection probabilities. It should be noted that the areas with higher detection probabilities do not necessarily align perfectly with areas of higher searcher efficiency; increased visibility of carcasses to humans might apply to scavengers, as well, resulting in higher carcass removal rates and potentially lower overall detection rates in those areas. Also, the model assumes that the distribution of carcasses with respect to distance from turbine is the same in searched and unsearched areas. It is not likely that carcass deposition distances vary between searched and unsearched areas. However, it is possible that animal movement after hitting the ground could vary significantly by ground cover, as injured animals move from exposed areas (like roads and pads) to more protected areas before dying. The degree to which this may introduce bias in practice is unknown but may be an interesting topic for future research.
Of course, accurate and precise estimation using R&P searches in conjunction with the methods described in this research has limitations. The weighted distribution, cake, glm, and anisotropic glm methods estimate mortality by estimating the spatial carcass distribution, which requires a reasonable sample size. For example, with only 23 turbines installed at WP2, an average of 2 carcasses per turbine equated to an average of eight total carcasses observed on R&P across all simulations. In a similar setting, the precision of the weighted distribution, cake, and glm methods will likely be poor, and the anisotropic glm inestimable. Therefore mortality estimation based on data solely collected on R&P should be avoided when few carcasses are expected to be observed and precise estimates of mortality are needed. If R&P searches are performed and few carcasses are observed, then it may be possible to use information from a nearby wind project with similar fauna, topography, wind, and turbines to estimate a spatial carcass distribution.
These methods may also fail to produce reasonable estimates if the sample of observed distances is non-representative of the underlying fall distribution. This could happen if anisotropy is present and either (i) major proportions of certain directions go unsearched, or (ii) the anisotropic carcass distribution at individual turbines interacts with the R&P configurations. For example, (i) could occur if every turbine access road came in from the same direction. In this instance there would be no way of testing or accounting for an anisotropic carcass distribution, because of the lack of information in the other directions. Situation (ii) is possible to encounter in practice, because road engineers often take topography into consideration and winds respond to topography. This could be problematic if, for example, all access roads are parallel to prevailing wind direction, although access direction may not be constant. In this situation, the carcass locations obtained from R&P searches would clearly not be representative. Researchers should be able to identify such situations by consulting with the engineers that built the project or using wind direction data collected at individual turbines.
Mortality estimation based entirely on R&P searches may not be applicable in all settings, but when it is, it can provide an efficient way to gather useful estimates at a relatively low cost. The simulations preformed in this research demonstrate that it is possible to obtain accurate and precise mortality estimates from data collected entirely on R&P searches under a wide variety of settings.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.