West-wide Flammulated Owl Analysis
3/30/2022
Authorship is still being negotiated. Requirements include at least 2, and preferably 3, of the following: 1) Funding, 2) Survey Design and Management, 3) Fieldwork Implementation and Management, 4) Data Analysis, 5) Manuscript Writing (more than just reviewing), and 6) Manuscript Editing. As this program involves 15 - 18 different datasets, these six categories should be considered with the multi-program scope in mind, not a single local program.
Introduction
The Flammulated Owl (Psiloscops flammeolus) is a small insectivorous owl that is common in the mixed conifer forests of the western United States and Mexico (Linkhart and McCallum 2020). It subsists almost entirely on insects, especially moths and beetles, and appears to be highly migratory (Linkhart and McCallum 2020). The species is a secondary cavity nester primarily associated with forests of commercially valuable trees (e.g., ponderosa pine [Pinus ponderosa] and Douglas fir [Pseudotsuga menziesii]), and timber management practices may influence its viability (Linkhart and McCallum 2020).
The Flammulated Owl is ranked as G4 (apparently secure) globally but is considered critically imperiled (state ranking of S1) in Wyoming and South Dakota, vulnerable (S3) in eight western states, Texas, and British Columbia, and apparently secure (S4) in Arizona and Colorado (NatureServe 2020). The current population estimates used for the Partners in Flight ranking, Breeding Bird Surveys, and eBird are not well suited for this species. As a result, there is very low confidence in the global estimated population size of 4,900-5,500 individuals.
In 2009 – 2012, multiple organizations across the West completed surveys for Flammulated Owls in at least six states. Each of these survey efforts utilized similar, but not identical protocols. However, a region-wide analysis of this data was never completed. Since that time other studies have been implemented (north-central Idaho, Wyoming), but each has only been analyzed locally. This current analysis effort was kicked off to integrate as many of these programs as possible into a west-wide view of the habitat use and population of this often elusive species.
Methods
Many organizations implemented programs contributing to these analyses. In most cases, consistent implementation was performed, allowing for the integration of the data used in this report. However, some organizations deviated from the consistent implementation to better meet local objectives. We strived in the development of this analysis and report to address these issues, integrating analyses when possible, but separating when appropriate. Some fidelity may have been lost in this process, so this report should not replace local analyses that were able to implement consistency end-to-end.
Study Areas
Surveys were implemented within the geography including these 9 states in the western and central United States (Figure 1). Surveys were only performed within six of the states.
Figure 1. Geographic locations of contributed datasets. Ownership includes a mixture of primarily National Forest land, with some Bureau of Land Management and state lands included.
Datasets
We assembled datasets from many different organizations that participated in the initial survey effort and those that have implemented projects since that time using compatible point protocols (Table 1).
| Ownership | State | Year | Sampling.Unit | Number.of.Transects | Number.of.Points | Number.of.Visits |
|---|---|---|---|---|---|---|
| Roosevelt NF | CO | 2010 | 4 km2 | 16 | 16 | 1 - 3 |
| Ashley NF | UT | 2010 | 4 km2 | 15 | 16 | 1 - 3 |
| Shoshone BLM | ID | 2010 | 4 km2 | 4 | 16 | 1 - 3 |
| Sawtooth NF | ID | 2010 | 4 km2 | 22 | 16 | 1 - 3 |
| Plumas NF | CA | 2010 | 4 km2 | 19 | 16 | 1 - 3 |
| Wasatch-Cache NF | UT | 2011 | 4 km2 | 16 | 16 | 1 - 2 |
| White River & Routt NF | CO | 2009 | 1 km2 | NA | 4 | 1 |
| Caribou-Targhee NF | ID | 2011 | 1 km2 | 40 | 4 | 1 - 3 |
| USFS Region 2 (9 forest) | CO & WY | 2011 | 1 km2 | NA | 4 | 1 - 3 |
| Medicine Bow NF | WY | 2012 | 1 km2 | NA | 4 | 1 - 3 |
| Nez Perce-Clearwater NF | ID | 2018 | 1 km2 | 39 | 9 | 1 |
| Sawtooth NF | ID | 2010 | Road Transect | 7 | 6 - 12 | 1 |
| Wyoming State-wide | WY | 2019 | Road Transect | 75 | 1 - 13 | 1 |
| Washington State-wide | WA | 2012 | Road Transect | 31 | 10 | 3 |
| Caribou-Targhee NF | ID | 2011 | Road Transect | NA | 10 | 1 - 3 |
| Sequoia NF | CA | 2011 | Indiv. Points | NA | 1 | 1 - 3 |
Figure 2. Geographic locations of contributed dataset (zoomable)s. Ownership includes a mixture of primarily National Forest land, with some Bureau of Land Management and state lands included.
Field Methods
Spatially-balanced Grid Surveys
All programs implementing off-road surveys utilized a Generalized Random Tessellation Stratification (GRTS) sampling procedure (Stevens and Olsen 2004). The survey unit was a square grid cell of either 1 km2 or 4 km2 in area. All programs selected grids within a predefined stratum definition that were all based upon some form of potential habitat. For example, most included some form of potential habitat including ponderosa pine, Douglas-fir, quaking aspen (Populus tremuloides), mixed xeric forest, lodgepole pine (Pinus contorta), subalpine fir (Abies lasiocarpa), Douglas-fir/lodgepole pine, mixed conifer/broadleaf, mixed conifer forest, juniper (Juniperus spp.), coniferous riparian, and broadleaf riparian. However, programs differed in the source of this data and the quantity of potential habitat required to be within the grid cell. As the detail of these stratifications have been lost with time, we are careful in these analyses to maintain a stratum for each potential grid selection procedure and not use the initial stratum definition in the reported results within this report. We instead implemented a post-stratification procedure for all habitat composition reporting.
Our analyses allowed up to three visits to each grid. Some programs attempted to survey all grids three times (e.g., Washington and Wyoming statewide programs), most used the GRTS process to select a subset of grids to be surveyed up to three times, and some program implemented only single visits to survey grids (e.g., Nez Perce-Clearwater National Forest). Our analyses consistently acknowledge the possibility of three repeat visits, even for programs with only a single visit.
Road-based Surveys
Road-based surveys were conducted similarity to grid-based surveys with the exception that points were laid out linearly along road or trail features instead of in randomly selected grid locations. Road-based surveys implemented the same point protocol as grid-based surveys. The road-based surveys were selected using a similar habitat stratification, or in some cases, such as within the Sawtooth National Forest, they were positioned in relation to sampled grid-cells for a comparison between methods.
Point Survey Protocol
We used a 10-minute point survey protocol developed by Partners In Flight - Western Working Group (Fylling et al. 2010) for standardized nocturnal surveys. All surveys started 30 minutes after sunset and lasted no more than six hours. The protocol consisted of five listening periods each separated by 30 seconds of Flammulated Owl call broadcasting. Surveys began with two minutes of silent listening (interval 1) followed by 30 seconds of call playback, followed by alternating periods of one-and-a-half minutes of silence and 30 second playbacks. This design allowed the survey period to be broken up into five equal intervals, which could later be used in occupancy estimation. We used a Foxpro NX3 (or similar) handheld game caller to broadcast the Flammulated Owl call (the same recording as Cilimburg 2007 and Smucker et al. 2008) and we standardized the volume level to 100-110 db @ 1 m using a handheld sound level meter. We programmed the game caller for the whole 10-minute survey period, so surveyors only had to turn on the caller once and were minimally distracted. During each playback session, we broadcasted calls in each of the four cardinal directions.
For each point we recorded survey start and end times, sky condition, wind condition, and noise conditions. When a Flammulated Owl was detected we recorded its estimated distance, bearing, and the interval and minutes in which it was detected. Bearing was determined by using a compass and we estimated distance to the best of our ability. When possible, we tried to calibrate our distance estimations by walking off the distance to a Flammulated Owl and recording it with our GPS unit. Most programs, but not all, utilized two surveyors for safety precautions and observers were allowed to communicate with each other to ensure agreement about detections.
Statistical Methods
Occupancy Rates and Habitat Associations
We used single species Spatial Occupancy Modeling to evaluate occupancy rates and habitat associations (Doser et al. 2021). We included year (as a factor) and stratum (survey program - e.g., national forest or state) as random intercepts in all models. As our point surveys consisted of five intervals, one silent and four broadcast/listen, we included a interval predictor specifying silent (0) and broadcast (1). We included this predictor in all models. Thus, our “null” model included the broadcast predictor, and the two random effects.
Due to the computationally intensive nature of this analysis, we performed model selection in stages. In stage one we focused on evaluating predictor variables possibly influencing the probability of detecting an owl given that one was present at the point. We evaluated julian date (day of year), start time (times after midnight were adjusted by adding 24 hours - 1:15am = 2515), and the noise level at the point. We evaluated both linear and quadratic effects of julian date, start time and noise. We compared models of all combinations of these potential predictors, selecting the best-supported model using the Widely Applicable Information Criterion (WAIC; Watanabe 2010).
For stage two, we chose just the predictors influencing the probability of detection in the best-supported model from stage one, while evaluating potential covariates influencing occupancy on the grid or transect. Once again, we selected the best supported model using WAIC.
Habitat Suitability and Predicted Distribution
For the predicted distribution of Flammulated Owls we used a machine learning approach called Random Forests (Breiman, L. 2001, 2002). This method has been shown to produce more reliable estimates than other machine learning approaches such as MaxEnt or Boosted Regression Trees (Cutler et al. 2007). Within each set of results we summarize the data used in the model, present the model fit using a Receiver Operating Characteristic (ROC) curve and Area Under the Curve (AUC; Fig. 2), present the importance of each variable influencing the model, plot the effect of each variable used in the final model, generate a west-wide prediction for the species, and lastly report the distribution and acres of this prediction.
Figure 3. Example Receiver Operating Characteristic curve with False Positive Rate and True Positive Rate. The diagonal shows the performance of a random classifier. Three example classifiers (blue, orange, green) are shown. Drawn by CMG Lee. Used unmodified under Creative Common Attribution - https://commons.wikimedia.org/wiki/File:Roc_curve.svg
We first prepared the predictor data at a resolution of 800 m by 800 m, using the USGS Digital Elevation Model as the template. We aligned numeric datasets to this template using bilinear interpolation, and aligned categorical datasets using nearest neighbor methods. We removed variables due to multi-colinearity, which occurred primarily within the proposed climate data (Evans and Murphy 2019). We used all presence points for the Flammulated Owls, and randomly selected twice that number of psuedo-absence points from the pool of negative survey data. This selection helped to ensure we met the “Rule-of-Thirds” guideline between positive and negative data (Evans and Murphy 2019). We then created a Random Forest using 1000 trees. We repeated the process nine more times, each time selecting a new pool of psuedo-absence points (Evans and Murphy 2019). These ten forests were averaged together and assessed to choose the variables of interest in the final model (i.e., model selection). We repeated the forest creation with the final set of predictors, once again generating ten random forests, each with 1000 trees. This set of averaged forests is the basis of all reported results.
Statistical Software
TBD…Update to include the new spOccupancy software…
Results
We completed 3360 point surveys, located on 568 transects, on 18 sampling strata, within 6 western states.
For fitting spatial occupancy models we used the following settings:
## Number of chains: 5
## Burn-in samples: 2000
## Chain thinning factor: 4
## Batch size: 400
## Batch length: 25
## Samples per chain (batch size * length) 10000
## Nearest Neigbor fitting (Datta et al. 2016): TRUE
## Number of Neighbors: 8Phase One: Choosing variables influencing the probability of detecting at least one Flammulated Owl, given that one was present.
| Model | WAIC | Delta |
|---|---|---|
| Psi(.)p(Broadcast + Julian + Julian^2 + Noise + Start + Start^2) | 3296.42 | 0.00 |
| Psi(.)p(Broadcast + Julian + Julian^2 + Noise + Start) | 3297.60 | 1.17 |
| Psi(.)p(Broadcast + Julian + Noise + Start) | 3299.02 | 2.60 |
| Psi(.)p(Broadcast + Julian + Noise) | 3300.01 | 3.59 |
| Psi(.)p(Broadcast + Julian + Noise + Start + Start^2) | 3300.60 | 4.17 |
| Psi(.)p(Broadcast + Julian + Julian^2 + Noise) | 3301.06 | 4.63 |
| Psi(.)p(Broadcast + Julian + Start + Start^2) | 3337.94 | 41.51 |
| Psi(.)p(Broadcast + Julian + Julian^2) | 3338.37 | 41.94 |
| Psi(.)p(Broadcast + Julian + Start) | 3348.87 | 52.44 |
| Psi(.)p(Broadcast + Julian + Julian^2 + Start) | 3349.76 | 53.33 |
| Psi(.)p(Broadcast + Julian) | 3354.66 | 58.24 |
| Psi(.)p(Broadcast + Julian + Julian^2 + Start + Start^2) | 3357.02 | 60.59 |
Phase Two: Choosing variables influencing the occupancy of a point by Flammulated Owls.
| Model | WAIC | Delta |
|---|---|---|
| Psi(Elevation + Elevation^2 + Canopy)p(Broadcast + Julian + Julian^2 + Noise + Start + Start^2) | 3281.80 | 0.00 |
| Psi(Elevation + Elevation^2)p(Broadcast + Julian + Julian^2 + Noise + Start + Start^2) | 3293.16 | 11.37 |
| Psi(Elevation + Canopy + Canopy^2)p(Broadcast + Julian + Julian^2 + Noise + Start + Start^2) | 3293.73 | 11.94 |
| Psi(Elevation + Elevation^2 + Canopy + Canopy^2)p(Broadcast + Julian + Julian^2 + Noise + Start + Start^2) | 3296.07 | 14.27 |
| Psi(Canopy + Canopy^2)p(Broadcast + Julian + Julian^2 + Noise + Start + Start^2) | 3297.81 | 16.02 |
| Psi(Elevation)p(Broadcast + Julian + Julian^2 + Noise + Start + Start^2) | 3300.35 | 18.55 |
| Psi(Elevation + Canopy)p(Broadcast + Julian + Julian^2 + Noise + Start + Start^2) | 3300.38 | 18.58 |
| Psi(Canopy)p(Broadcast + Julian + Julian^2 + Noise + Start + Start^2) | 3303.40 | 21.61 |
Here is a summary of the best supported model from the analysis:
##
## Call:
## spPGOcc(occ.formula = ~(1 | stratum) + elev + I(elev^2) + canopy +
## I(canopy^2), det.formula = topDetForm, data = flowList, cov.model = "exponential",
## NNGP = NNGP, n.neighbors = n.neighbors, n.batch = n.batch,
## batch.length = batch.length, accept.rate = 0.43, n.omp.threads = n.omp.threads,
## verbose = T, n.report = n.report, n.burn = n.burn, n.thin = n.thin,
## n.chains = n.chains, k.fold.threads = k.fold.threads)
##
## Samples per Chain: 10000
## Burn-in: 2000
## Thinning Rate: 4
## Number of Chains: 5
## Total Posterior Samples: 10000
## Run Time (min): 26.5895
##
## Occurrence (logit scale):
## Mean SD 2.5% 50% 97.5% Rhat ESS
## (Intercept) -0.5475 0.4461 -1.4135 -0.5594 0.3408 1.0603 213
## elev 0.1347 0.1950 -0.2275 0.1295 0.5387 1.0175 1243
## I(elev^2) -0.3520 0.1153 -0.5927 -0.3461 -0.1421 1.0834 1010
## canopy 0.0166 0.1135 -0.2164 0.0190 0.2363 1.0155 1467
## I(canopy^2) 0.0101 0.1084 -0.2041 0.0098 0.2232 1.0037 3245
##
## Occurrence Random Effect Variances (logit scale):
## Mean SD 2.5% 50% 97.5% Rhat ESS
## stratum 2.8621 2.0394 0.8482 2.3318 8.0345 1.2447 264
##
## Detection (logit scale):
## Mean SD 2.5% 50% 97.5% Rhat ESS
## (Intercept) -2.0053 0.1724 -2.3459 -2.0052 -1.6738 1.0135 1154
## broad -0.1532 0.1264 -0.3941 -0.1551 0.0952 1.0011 6929
## julian 0.0388 0.0837 -0.1257 0.0397 0.2016 1.0052 2536
## I(julian^2) 0.0028 0.0072 -0.0112 0.0029 0.0168 1.0041 2746
## noise -0.5544 0.0785 -0.7123 -0.5532 -0.4013 1.0002 4612
## start 0.0302 0.0728 -0.1111 0.0315 0.1726 1.0010 6178
## I(start^2) 0.0459 0.0469 -0.0453 0.0456 0.1380 1.0012 5701
##
## Spatial Covariance:
## Mean SD 2.5% 50% 97.5% Rhat ESS
## sigma.sq 0.9945 1.178 0.1984 0.55 4.2164 3.3773 29
## phi 10760.0890 6308.111 439.5363 10762.44 21110.8019 1.0027 1639
## Mean SD 2.5% 50% 97.5% Rhat ESS
## sigma.sq 0.9945 1.178 0.1984 0.55 4.2164 3.3773 29
## phi 10760.0890 6308.111 439.5363 10762.44 21110.8019 1.0027 1639##
##
## Result of the Freeman-Tukey Goodness-of-Fit test (p values >0.05 pass the fit test.).##
## Call:
## ppcOcc(object = top, fit.stat = "freeman-tukey", group = 1)
##
## Samples per Chain: 10000
## Burn-in: 2000
## Thinning Rate: 4
## Number of Chains: 5
## Total Posterior Samples: 10000
##
## Bayesian p-value: 0.0033
## Fit statistic: freeman-tukey
## NULL
Figure 3. Predicted effect of elevation on the probability of point occupancy by Flammulated Owls.
Figure 4. Predicted effect of canopy cover on the probability of point occupancy by Flammulated Owls.
Graphs for probability of detecting an owl, given that one was present:
Figure 5. Predicted effect of interval broadcast type on the probability of detecting at least one Flammulated Owl given that one was present.
Figure 6. Predicted effect of survey date on the probability of detecting at least one Flammulated Owl given that one was present.
Figure 7. Predicted effect of noise on the probability of detecting at least one Flammulated Owl given that one was present.
Figure 8. Predicted effect of survey start time on the probability of detecting at least one Flammulated Owl given that one was present.
West-wide prediction
Random Forest analysis.
We have used a combination of geographic data (e.g., elevation), climate data, and forest structural information (e.g., canopy cover) to build our machine learning models (Table 4).
| Name | Units | Range | Mean | SD | Description |
|---|---|---|---|---|---|
| Elev | m | Elevation above mean sea level | |||
| Canopy | % | Canopy closure | |||
| ppt | mm | Average annual precipitation (rain+melted snow) | |||
| tmin | °C | Average annual minimum temperature | |||
| tmax | °C | Average annual maximum temperature | |||
| tmean | °C | Average annual mean temperature | |||
| tdmean | °C | Average annual mean dewpoint temperature | |||
| vpdmin | hPa | Average annual minimum vapor pressure deficit | |||
| vpdmax | hPa | Average annual maximum vapor pressure deficit | |||
| soltotal | MJ m-2 day-1 | Average annual global shortwave solar radiation received on a horizontal surface | |||
| solslope | MJ m-2 day-1 | Average annual global shortwave solar radiation received on a sloped surface | |||
| solclear | MJ m-2 day-1 | Average annual global shortwave solar radiation received on a horizontal surface under clear sky conditions | |||
| soltrans | dimensionless, 0-1 | Atmospheric transmittance (cloudiness) |
The final dataset used to create the Random Forest model contained 309 presence points. Each run of the Random Forest model was augemented with 618 randomly selected psuedo-absence points from the set of points with no detections, to meet the ‘Rule of Thirds’.
The following Confusion Matrix presents the error rates of the final model as evaluated by data that was set aside for testing (Table 15). This includes a summation of all ten created forests (hence why the point counts are so high). The first row presents the points where no owls were detected, and the error rate indicating how many were predicted to be suitable for Flammulated Owls. Row two represents the points where owls were detected and how many were predicted to be non-suitable for occupancy.| Predicted Absence | Predicted Presence | Class Error | |
|---|---|---|---|
| Psuedo-absence | 5233 | 947 | 0.15 |
| Presence | 1274 | 1816 | 0.41 |
The Receiver Operating Characteristic (ROC) curve maps the true positives (y-axis) against the false positives (x-axis). The further the curve is up to the left, the better (see example in the methods section). The Area Under the Curve (AUC) measures the model performance. The closer to 1.0 the better.
The variable importance graph shows the sorted order of predictor variables by how big a role they played in the final model. The measure is how much the model degraded when that particular variable was removed (relative effect). Variables within these importance charts are usually more separated in their importance. The similararity of importance across these variables are suggestive that different variables are used in different geographic areas, thus balancing each other out overall.
The following prediction maps show the estimate effect of each covariate upon the predicted suitability of the habitat to support Flammulated Owls.These results should be used carefully as machine learning techniques are best at predicting the overall suitability, but are less effective than occupancy models at measuring individual variables influences.
We separated the habitat suitability predictions across the nine-states study area, into five somewhat arbitrary equal interval categories (Table 16). As expected, most of the landscape is not suitable for Flammualted Owls during the breeding season. Should we restrict this to some other landscape designation (e.g., canopy cover above 15%??)? Thoughts?
| Score | Hectares | Acres | Proportion |
|---|---|---|---|
| Very Good (0.8 < X) | 24,576 | 60,729 | 0.00 |
| Good (0.6 < X <= 0.8) | 1,055,040 | 2,607,057 | 0.00 |
| Moderate (0.4 < X <= 0.6) | 48,333,888 | 119,435,454 | 0.20 |
| Poor (0.2 < X <= 0.4) | 167,122,880 | 412,968,993 | 0.69 |
| Very Poor (X <= 0.2) | 25,212,992 | 62,302,564 | 0.10 |
West-wide habitat quality in graph form (same data as Table 6).
West-wide Flammulated Owl habitat suitability prediction plot. Zoom and drag enabled.
Discussion
Stay tuned after all analysis is worked out.
Acknowledgements
We thank the entire Partners in Flight Western Working Group for providing the impetus and coordination for this west-wide survey of Flammulated Owls. In particular, Mike Green, Dave Krueper, Kevin Kritz, and Eric Kirscher were instrumental in raising the core funding from the U.S. Fish and Wildlife Service Division of Migratory Birds and Habitat Programs.
Additional funding for this effort came from many sources, including
- Bureau of Land Management - Shoshone Field Office
- Idaho Department of Fish and Game Wildlife Diversity Program
- Nez Perce-Clearwater National Forest
- Teton Raptor Center
- Wyoming Department of Game and Fish
- …
I need help here…
Most importantly, we thank all the field technicians for persevering through many grueling nights of surveys!
Literature Cited
Cilimburg, A. (2007). Northern region landbird monitoring program: Flammulated Owl protocol. Avian Science Center, University of Montana, Missoula, Montana, USA.
Datta, A., S. Banerjee, A.O. Finley, and A.E. Gelfand. (2016) Hierarchical Nearest-Neighbor Gaussian process models for large geostatistical datasets. Journal of the American Statistical Association, doi: 10.1080/01621459.2015.1044091.
Doser, J. W., Finley A. O., Kéry, M., & Zipkin E. F. (2021). spOccupancy: An R package for single species, multispecies, and integrated spatial occupancy models arXiv preprint arXiv:2111.12163
Evans, J. S., and M. A. Murphy (2019). rfUtilities: An R package for Random Forests Model Selection and Performance Evaluation. https://github.com/jeffreyevans/rfUtilities.
Fylling, M. A., J. D. Carlisle, A. B. Cilimburg, J. A. Blakesley, B. D. Linkhart, and D. W. Holt (2010). Flammulated Owl Survey Protocol. Partners in Flight Western Working Group. 10p.
Gahbauer, M. A., R. A. Miller, N. Paprocki, A. Morici, A. C. Smith, and D. A. Wiggins (2021). Status and monitoring of Short-eared Owls (Asio flammeus) in North and South America. Airo 29:115–142.
Liaw, A., and M. Wiener (2018). randomForest: An R package for Breiman and Cutler’s Random Forests for Classification and Regression. https://www.stat.berkeley.edu/~breiman/RandomForests/.
Linkhart, B. D., and D. A. McCallum (2020). Flammulated Owl (Psiloscops flammeolus), version 1.0. Page in A. F. Poole, editor. Birds of the World. Cornell Lab of Ornithology, Ithaca, NY, USA.
Miller, R. A., N. Paprocki, M. J. Stuber, C. E. Moulton, and J. D. Carlisle (2016). Short-eared Owl (Asio flammeus) surveys in the North American Intermountain West: utilizing citizen scientists to conduct monitoring across a broad geographic scale. Avian Conservation and Ecology 11. http://www.ace-eco.org/vol11/iss1/art3/.
NatureServe 2020
Nichols, J. D., L. L. Bailey, A. F. O’Connell Jr., N. W. Talancy, E. H. C. Grant, A. T. Gilbert, E. M. Annand, T. P. Husband, and J. E. Hines (2008). Multi-scale occupancy estimation and modelling using multiple detection methods. Journal of Applied Ecology 45:1321–1329.
Paluszynska, A., B. Przemyslaw, and J. Yue (2020). randomForestExplainer: An R package for Explaining and Visualizing Random Forests in Terms of Variable Importance. https://github.com/ModelOriented/randomForestExplainer.
Pavlacky Jr., D. C., J. A. Blakesley, G. C. White, D. J. Hanni, and P. M. Lukacs (2012). Hierarchical multi-scale occupancy estimation for monitoring wildlife populations. Journal of Wildlife Management 76:154–162.
Smucker, K., A. Cilimburg, and M. Fylling. (2008). Northern region landbird monitoring program 2008 Flammulated Owl Surveys Final Report. Avian Science Center, Univ. of Montana, Missoula.
Stevens, D. L., and A. R. Olsen (2004). Spatially balanced sampling of natural resources. Journal of the American Statistical Association 99:262–278.
Watanabe, Sumio. 2010. “Asymptotic Equivalence of Bayes Cross Validation and Widely Applicable Information Criterion in Singular Learning Theory.” Journal of Machine Learning Research 11 (12).