Owls provide ecosystem services and play crucial roles in the environment making them important to monitor and study. However, standardized methods for most species do not exist, and we lack understanding of the effects of many environmental variables and call-broadcast on detection probability of owls. We performed a multispecies occupancy analysis of owl monitoring data collected from 2004 – 2013 across the state of Maine to examine the effects of environmental variables, conspecific and heterospecific call-broadcast, and general survey protocols on detection of three forest owls: Northern Saw-whet Owl (Aegolius acadicus), Barred Owl (Strix varia), and Great Horned Owl (Bubo virginianus). We found that environmental variables such as cloud cover, precipitation, temperature, time of night, and wind had species-specific effects on detection probability, and ambient noise decreased detection probability for all species. We did not find support for effects of snow cover on detection of any species. We also found that conspecific call-broadcast increased detection of each species, while heterospecific call-broadcast had variable affects. Specifically, we found that Long-eared and Barred Owl broadcast increased the detection of Northern Saw-whet Owl, and our results suggest additional heterospecific effects may exist. Our study showed that compared to the Maine Owl Monitoring Program, surveys examining all three of our focal species can increase efficiency and lower disturbance by only broadcasting Long-eared and Barred Owl calls during a 10-minute survey. We recommend that future owl surveys take into account species-specific effects of conspecific and heterospecific call-broadcast, and use our results when designing survey protocols that include one or more of our focal species. 

Maine Owl Monitoring Program

MOMP consisted of a pilot study during 2002 and 2003 to evaluate methods and study design for a primary study conducted from 2004 to 2013. The main goals of the program were to (1) examine owl abundance and distribution in Maine, (2) evaluate survey methods and monitoring effectiveness, and (3) build a long-term citizen-science project through development of a network of skilled volunteers (Hodgman & Gallo 2004). 
    Survey routes were established across the state (Figure 1) using the Maine DeLorme Atlas (Olathe, KA) grid system to delineate survey blocks that were further divided into four equal quadrants. Within each quadrant, two survey routes were established, and participants were asked to sample their assigned route once during a 6.5-week survey period that started on the first Friday of March each year. All surveys were conducted between 2400 – 0500 hours, and not all routes were sampled equally or continuously across the study years. Because we focus exclusively on detection probability, which is survey-specific (described further below), we assume that our inferences related to survey-level detection are independent of variable route sampling among years.
    Along each route, participants established ten survey locations that were separated by a minimum distance of 1 mile, which we assumed was sufficient for independence among survey points. We believe this assumption is valid, as a distance of 1-mile between points exceeds mean home range size for each of our focal species (Mazur and James 2000, Churchill et al. 2002, Bennett and Bloom 2005), making it unlikely that the same owl was detected at multiple points in a single night. During each survey, participants recorded start time, and quantified temperature, cloud cover, snow cover, wind, noise, and precipitation. Temperature was estimated using handheld thermometers or from regional weather records. Cloud cover represented the percentage of the sky covered by clouds, estimated to the nearest 10%. Snow cover was recorded along a categorical gradient where 1 = no snow cover, 2 = patchy or partial snow cover, and 3 = complete ground cover. Wind speed was approximated following the Beaufort scale from 0 – 5 where 0 = < 1 mph, 1 = 1 – 3 mph, 2 = 4 – 7 mph, 3 = 8 – 12 mph, 4 = 13 – 18 mph, and 5 = 19 – 24 mph winds. Noise was quantified on a scale from 1 – 4, where 1 = relatively quiet, 2 = moderate noise but not affecting ability to hear owls, 3 = loud noise that may affect detection of owls, and 4 = excessive loud noise that probably affects detection of owls. Precipitation was scored categorically as either 0 – 1, where 0 = none or 1 = drizzle. Surveys were not conducted in high winds (> 3 on the Beaufort scale) or during precipitation > 1 on the scale. 
    Surveys began with 3 minutes of passive listening, followed by a call-broadcast 20 seconds in duration and a 1-minute 40-second listening period. At the start of the 5th minute of the survey, a second call of a different species was used, followed by a 5-minute 40-second listening period  . This longer listening period was designed to allow more time for Barred and Northern Saw-whet Owls to respond before playing Great Horned Owl (a potential predator of smaller owls) in order to mediate any possible negative effects of Great Horned Owl call-broadcast on Barred and Northern Saw-whet Owl detection. At minute 11 a final call of a third species was broadcast, followed by another 1-minute 40-second listening period, the end of which marked the completion of the survey. All calls were broadcast for 20-seconds. During 2004 – 2006, Northern Saw-whet, Barred and Great Horned Owl calls were played, in that order. From 2007 – 2013, the Northern Saw-whet Owl call was replaced with a Long-eared Owl call. During each survey, participants recorded all owls heard or seen in each of the 16 sampling intervals, which consisted of 60-second listening intervals (n=10), 20-second call-broadcast intervals (n=3), and 40-second listening intervals (n=3; the remainder of the minute post-broadcast). During initial years, participants used cassette tapes with recordings of the calls of each focal species, while during later years, the option to use a 13-minute recording that was structured to play different owl calls at the designed times was offered. Participants chose to either use this new method, or continue using the cassette tapes. 

Modeling Procedures

We constructed species-specific detection histories, using owl detections during each listening interval (i.e., 20-, 40-, or 60-second listening period) as temporally-repeated observations within each 13-minute survey. We examined the effects of environmental covariates and call-broadcast on species detection probabilities (p) using a multispecies occupancy framework (Rota et al. 2016) in program MARK (White and Burnham 1999, Skelly et al. 2018). Under this formulation, the observation (y) of a particular species i at a specific site j during any given survey k is described as
y_ijk |z_ij=Bernoulli(z_ij p_ijk)
where z reflects a vector of species-specific latent occupancy states that account for the probability of occurrence/co-occurrence for each species (see Rota et al. 2016 for additional details). Estimates of p are then modelled under a binomial distribution which, given the hierarchical nature of the model, are conditional on site-occupancy of the species in question as defined by z. In practice p can be derived under a generalized linear modelling framework as
logit(p_ijk )=β_0+β_1 X_1+β_2 X_2…+β_t X_t
where a logit link distribution is typically used, and species-, survey-, or site-specific heterogeneity in p may be accommodated by the effects (β) of a sequence of categorical or continuous predictor variables. 
    In the context of our survey design and data, we used the discrete listening periods (20-, 40- or 60-second intervals) as each k survey replicate, where p then reflected the probability of detection for each species during any single interval within the 13-minute survey. While occupancy models are often fit using repeated visits to a site, replicated observations within visits are another well-established survey design that may be used in an occupancy context (MacKenzie and Royle 2006). In this case, the occupancy parameter reflects the probability of species presence/absence during the survey, as the assumption of closure among survey replicates is almost certainly met (MacKenzie and Royle 2006, Rota et al. 2009). For our purposes, this survey design implies that occupancy reflected the probability that one or more owls of each species was present at the survey stop during the survey, and p gave the probability that an owl was observed and recorded (primarily through vocalization) during each listening interval, given presence. 
    Our objectives for this paper were strictly related to detection probability so we used a general structure for the occupancy component of the model that allowed independent annual estimates of occupancy for each of the three focal species. In practice, we fit this model structure using a year*species interaction, where year was modeled as a categorical fixed effect, and the interaction allowed occupancy to vary among years independently for each species. We also included an additional term in the model to allow for additive effects of co-occurrence for each species-pair, as well as all Northern Saw-whet, Barred, and Great Horned Owl three species in aggregate. This model accommodated differences in occupancy for each species/year combination, with constant rates of co-occurrence among species for the duration of the study. We elected not to incorporate explicit multi-year dynamics in the model (e.g., Hines et al. 2014) to reduce model complexity given the multi-species design, and because we felt accommodating the temporal dynamics with fixed year effects was appropriate given that our primary objective centered on detection probability rather than occupancy. We held the occupancy structure constant across all models and occupancy parameters reached convergence in all analyses, but for the purpose of this paper we did not interpret the occupancy terms further.
    To assess the fit of our data to the structure of the multi-species occupancy model, we used the occuMulti function in the package unmarked (Fiske and Chandler 2011) in program R to fit a multi-species occupancy model to our data, and then conducted a Chi-squared goodness-of-fit test (MacKenzie and Bailey 2004, Warton et al. 2017) using the parboot function to implement a parametric bootstrap with 500 replications (Fiske and Chandler 2011). We used a model with independent year effects on the single-species occupancy and detection parameters, and an additional intercept-only structure for the species interaction terms. This model form matched the base structure of all models we constructed using Program MARK (i.e., the general model without covariates). Results of these tests suggest the data conformed to the model structure, as bootstrapped replicates did not deviate from the expected distributions given the model (p=0.39), suggesting an adequate fit of the model to the data (Fiske and Chandler 2011).
    We only examined detection probability of Northern Saw-whet, Barred, and Great Horned Owl; all other species recorded during MOMP surveys had a sample size of < 100 detections, and were therefore not considered. We modeled all effects independently (i.e., a unique Beta coefficient) for each species because our goals were to understand species-specific changes in detection in response to environmental and call-broadcast effects. All covariates were z-standardized (mean = 0.0, SD = 1.0). We used Akaike information criterion adjusted for small sample size (AICc) to rank the detection models following a criterion of ∆AICc < 2.0 for models to be competitive. 

Modeling Environmental Covariates

We began model selection for the detection probability component of the analysis by running models with additive or interactive effects of species, ordinal date, and year, and tested quadratic effects of ordinal date. From this model set, we selected the best-performing model and used this as a base structure upon which we built all further models. We first constructed a correlation matrix of all covariates in program R (version 3.4.2), and found that no two covariates were highly correlated (Pearson’s r < 0.5). We then tested individual models with additive effects of the following 7 covariates on detection of each species: time of survey expressed as time since midnight, temperature, cloud cover, snow cover, wind, noise, and precipitation. We combined covariate effects that performed better than the base model into a singular comprehensive model, from which all inference was made. We used the Beta coefficients (β) and their 95% confidence intervals to further evaluate species-specific covariate support (95% CI does not overlap zero), and interpreted relative differences in effects among species (i.e., differences in species-specific β for each supported variable). 

Modeling Effects of Call-broadcast

To determine the effects of broadcast on detection of the three species of owl, we used the best model from our base model selection (combined species, date and year effects) and created individual models with effects for each of the four species’ calls, applied independently (i.e., unique β coefficients) to detection probability for each of the three focal species. We examined the effect of each species’ call-broadcast on detection of each species, where we modeled a unique detection probability for the 20-second period of call-broadcast and each subsequent sampling interval until the next species’ call began. For Great Horned Owl, the effect of conspecific? call-broadcast was modeled to include the 20-second call-broadcast and the remaining sampling intervals until the end of the survey. Effects   of Northern Saw-whet Owl and Long-eared Owl call-broadcast were limited to the years in which they were played (NSWO = 2004 – 2006; LEOW = 2007 – 2013). We accomplished this in practice using a dummy variable in the design matrix, where values of ‘1’ were applied only to sampling intervals during the years in which the call-broadcasts for each respective species were played. The resulting β coefficients for this effect reflected the differences in detection probability between call-broadcast intervals associated with the species during the years in which their calls were played, and all other intervals not associated with the call-broadcast for that species. Although detections of Long-eared Owl were too infrequent to analyze, we still included the effects of Long-eared Owl calls on our three focal species because they were played for seven years of the surveys, and it has been shown that heterospecific calls can affect detection of other owl species (Bosakowski and Smith 1998). 
    Based on preliminary results, we combined all four of call-broadcast structures   (the independent call-broadcast structures of Northern Saw-whet, Long-eared, Barred, and Great Horned Owl into one comprehensive model, and made further inferences and interpretations from this combined model. In this model we also included a general effect of the call broadcast interval (i.e., was any call actively being broadcast) on the detection of all species, because we reasoned that the speaker broadcast could muffle the sound of an immediately responding owl. Finally, we ran a post-hoc model to specifically examine the effects of Barred Owl call-broadcast on Northern Saw-whet Owl detection during the years in which Northern Saw-whet owls were not broadcast (2007 – 2013). This allowed us to evaluate whether the effect of Barred Owl call-broadcast on Northern Saw-whet Owl detection was still supported without the bias of an earlier broadcast of Northern Saw-whet Owl calls, given that Barred Owl calls were always broadcast after Northern Saw-whet owls during 2004 – 2006.

Modeling Survey Duration and Structure

To examine how effective surveys were at detecting species presence, and to determine the most efficient and effective survey structure for detecting each species singly as well as in aggregate, we calculated a cumulative detection probability, p*, as 
 〖p^*〗_i=1-(∏_(1:n)▒〖(1-p_k)〗)
which derived the probability that each i species was detected at least once during a survey of n listening intervals, where p is the probability of detection during a single listening interval (k). We used estimates of pk from the most competitive model identified by our call-broadcast model selection, and calculated p* under a variety of scenarios to explore how potential differences in call-broadcast structure and survey length could affect detection of each species during the survey. Given that models included a year effect, we used estimates of pk from years that were closest to the mean detection probability value for each species. We also assumed mean values for all environmental covariates for this assessment. 
    We first calculated a unique p* for each combination of successive survey lengths (e.g., 1 interval, 2 intervals, etc.) under the existing protocol for each of the three focal species, and we repeated this for both when Northern Saw-whet Owl was played first and when Long-eared Owl was played first. We then explored hypothetical alternative survey structures where we assumed only two species’ calls were broadcast (i.e., either Northern Saw-whet or Long-eared plus Barred Owl but excluding Great Horned Owl), and where we assumed only one species’ call was broadcast (Northern Saw-whet or Long-eared only but excluding the two other species; Table 1). We also explored a hypothetical survey structure consisting of only passive listening, to which we compared all alternative call-broadcast structures. In each of these hypothetical scenarios, we calculated p* for all possible survey lengths. We report change in cumulative detection probability as a function of increased survey duration, and evaluated differences in relative ability to detect each species among survey structures based on the number of intervals required for a species’ detection probability to reach values of both p* > 0.50 and p* > 0.95. In cases where two survey structures achieved the same p* during the same survey interval, we considered the survey structure with the fewest number of call-broadcasts to be more ideal, because use of fewer calls could limit disturbance and chance of human error. Collectively this exercise allowed us to evaluate optimal sampling designs for detecting owl presence during single visit surveys.