DOI: 10.5061/dryad.c866t1gkt
title: Data from: Divergent responses of native predators to severe wildfire and biological invasion are mediated by life history

description: This dataset supports an analysis of the relative and combined effects of
 severe wildfires and an early-stage barred owl (Strix varia) invasion on
 three native forest owl species in the Sierra Nevada, California, USA.
 Data were collected between 2018 and 2023 using a regional passive
 acoustic monitoring program, complemented by manipulative (lethal barred
 owl removals) and natural (severe wildfire) experiments. The dataset
 includes species-specific detection histories for flammulated owls
 (Psiloscops flammeolus), great horned owls (Bubo virginianus), and
 northern pygmy owls (Glaucidium californicum), as well as site-level
 covariates describing fire severity, barred owl removal treatment,
 topographic conditions, and survey effort. Derived variables include
 wildfire severity classifications and pre- and post-removal monitoring
 periods. The data can be reused to run Bayesian single-species
 "stacked" single-season occupancy models. All sensitive location
 information of forest owls has been omitted.
description: <strong>Study system</strong> Our
 study area encompassed &gt; 6,000 km2 in the northern Sierra Nevada,
 California, comprised primarily of publicly-managed lands (Lassen and
 Plumas National Forests) and some private lands. The area included a
 topographically complex mountain range that varied in elevation and
 vegetation type (though predominantly mixed conifer forest). Three
 significant environmental changes occurred in this region from 2018 to
 2022: 1) the North Complex Fire ( 1,220 km2 burned in late summer 2020,
 599 km2 [49%] at high-severity); 2) the Dixie Fire (3,740 km2 burned in
 late summer 2021, 2,090 km2 [56%] at high-severity); and 3) a barred owl
 lethal removal experiment (removals conducted across the entire study area
 primarily in summer 2019 described below) (Hofstadter et al. 2022). The
 extent and severity of each fire were both a departure from historical
 fire regimes and emblematic of emerging trends in the Sierra Nevada (Cova
 et al. 2023).  The lethal removal experiment was
 initiated after documentation of rapid population growth in the project
 area between 2017 and 2018 (Wood et al. 2020a); most removals were
 conducted in 2019 (54 birds), with 8 owls removed in 2018, 14 owls removed
 in 2020, and two in 2021 and 2022 (Hofstadter et al. 2022). We lethally
 removed barred owls and hybrids using a 12- or 20-gauge shotgun following
 field protocols established by Diller et al. (2014, 2016), and all
 removals were carried out by trained and permitted personnel from Sierra
 Pacific Industries and the University of Wisconsin (Hofstadter et al.
 2022). Of the 80 total removals, 67 were barred owls and 13 were barred
 owl x spotted owl hybrids. We grouped barred owls and hybrids due to their
 genetic similarity and small hybrid sample size, which prevented a
 separate analysis for barred owls and hybrids.
 <strong>Passive acoustic surveys</strong>
 To characterize the effects of severe fire and invasive barred
 owls on forest owl communities, we conducted passive acoustic surveys in
 the northern Sierra Nevada from early May to late July or the onset of
 wildfire in 2018, 2021, 2022, and 2023, bracketing the periods in which
 the North Complex Fire (August to September, 2020) and Dixie Fire (July to
 October, 2021) burned and barred owls were removed (survey data from 2019
 and 2020 was excluded because lethal removals and acoustic surveys
 occurred simultaneously, violating the assumption of closure). Our study
 relied on nocturnal surveys from 2000 to 0800 hours PDT, aligning with the
 activity periods of most owls and capturing the crepuscular behavior of
 pygmy owls (Sater et al. 2006). Although pygmy owls are primarily diurnal
 hunters, their frequent detections during the crepuscular period justified
 their inclusion in our analysis (Sater et al. 2006, Deshler et al. 2019).
 However, because our surveys did not extend into full daylight hours, some
 pygmy owl detections were likely missed. Passive
 acoustic surveys can provide detection/non-detection data required for
 occupancy modeling (Wood et al. 2019). Each year we surveyed the same 265
 hexagonal cells, each 400 ha, similar in size to spotted owl and barred
 owl territories in the region (Wood et al. 2020a). We systematically
 selected approximately 1 of 5 hexagonal cells from the grid in the Lassen
 and Plumas National Forests for surveying, ensuring they were (1)
 non-contiguous to minimize nonindependence among sites (e.g., detecting
 the same owls in multiple hexagonal cells) and (2) readily accessible by
 road. In each hexagonal cell, we deployed two autonomous recording units
 (ARUs; Swift One Recorder, K. Lisa Yang Center for Conservation
 Bioacoustics, Cornell Lab of Ornithology, Ithaca, NY) at locations with
 high topographic relief and a minimum spacing of 500 m. The elevations of
 sites surveyed ranged from 669 to 2262 m with a mean of 1573 m. In 2018,
 we deployed ARUs three times at each site for approximately one week per
 deployment, with an average interval of 26 days (minimum 14 days) between
 deployments. In 2021, 2022, and 2023, we deployed ARUs continuously for
 five weeks at each site. Although total survey effort was lower and the
 deployment schedule differed in 2018, all deployments occurred during the
 core breeding season (May - July) and at the same spatial locations across
 years. We also accounted for sampling effort differences by including
 sampling effort as a fixed effect in our detection models (see below). We
 deployed ARUs in the same locations each year. Equipped with an
 omnidirectional microphone, the ARUs recorded sound data at a sample rate
 of 32 kHz, 16-bit resolution, and gain of +33 dB.
 <strong>Processing acoustic survey
 data</strong> To identify forest owl
 vocalizations, we used the BirdNET algorithm, a deep convolutional neural
 network designed to identify &gt; 6,000 species by sound, most of them
 birds (Kahl et al. 2021). We used a custom version of BirdNET overfit to
 vocalizations of the species of interest in this study. Every 3 seconds,
 BirdNET produces a unitless numeric prediction score from 0 to 1 for each
 species. The prediction score reflects its relative confidence in
 correctly identifying the species, with higher scores representing higher
 confidence (Wood and Kahl 2024). To determine the
 occurrence of a species-specific vocalization, we selected minimum
 prediction score thresholds that would result in false positive rates of
 &lt; 1% for each species at the scale of a weekly secondary sampling
 period. To do so, we first used the program Raven Pro (K. Lisa Yang Center
 for Conservation Bioacoustics, Cornell Lab of Ornithology) to manually
 review a random subset of at least 200 hours of audio data each year for
 each species, with all selected hours containing at least one BirdNET
 prediction with a score &gt; 0.1. When reviewing the acoustic data, we
 assigned a code of ‘1’ to an hour-long sample if it contained at least one
 confirmed true positive or ‘0’ if it did not contain at least one
 confirmed true positive. In each hour-long sample, we then counted the
 number of BirdNET predictions above a range of scores (0.1, 0.2, …, 0.9;
 0.91, ..., 0.99). Using logistic regression, we estimated the probability
 of a false positive in a random hour of acoustic data as a function of the
 number of BirdNET predictions over each score, where we treated the true
 positive/false positive status as a binary response (i.e., false positive
 =1; true positive = 0) and the number of BirdNET predictions in an hour
 above a given score threshold as the predictor (lme4; Bates et al. 2015).
 We estimated false positive rates for each potential threshold at the
 scale of a secondary sampling period using the equation 1-(1-FP)n, where
 FP was the hourly false positive rate, and n was the number of hours in a
 typical secondary sampling period (n = 84 hours; see below). We selected
 thresholds that yielded an expected &lt; 1% false positive rate for
 each species by year combination. As a result of
 transitioning to different microphones between 2018 and 2021 and
 maintaining the same microphones from 2021 onward, we decided to uniquely
 assess audio data from the 2018 and 2021 seasons, with call rate
 thresholds in 2022 and 2023 matching those from 2021. To increase the
 sample size before and after wildfire and lethal barred owl removals, we
 manually validated predictions for pygmy owls and flammulated owls with a
 15% and 42% false positive rate, respectively. We initially included
 western screech owls and saw-whet owls in our analyses, but had to exclude
 them due to low detection occurrences, both at barred owl lethal removal
 and severely burned sites. <strong>Creating
 encounter histories for occupancy modeling</strong>
 For flammulated and pygmy owls, we created encounter histories
 from the four years of passive acoustic survey data at the level of
 individual ARUs, as these species have relatively small home ranges and
 quiet calls, meaning ARUs separated by ≥ 500m could be treated as
 independent sampling units. For great horned owls, a larger species with a
 larger home range and louder calls, we created encounter histories at the
 level of sampling hexagons (i.e., two ARUs), as multiple ARUs within a
 400-ha hexagon may record calls from the same individual. We divided the
 survey season into eleven one-week-long secondary sampling periods,
 starting on May 7th and ending on July 23rd. We assumed population closure
 throughout the year. If an ARU did not record during a specific secondary
 sampling period, we treated that period as null. However, if an ARU was
 active for only part of a secondary sampling period, we included the data
 from that period and treated survey hours as a detection covariate to
 account for unequal sampling (see below).
 <strong>Environmental predictors</strong>
 To calculate predictor variables, we created circular buffers
 around all sites that represent the approximate size of each species’ home
 range (Giese and Forsman 2003, Bennett and Bloom 2005, Yanco and Linkhart
 2018). For flammulated and pygmy owls, we buffered individual ARUs at 50
 ha and 300 ha, respectively. To represent larger home ranges in great
 horned owls, we delineated buffers of 700 ha around the center of the
 sampled 400-ha hexagons. Next, we utilized severely burned area (with at
 least 75% overstory mortality) as estimated by the Monitoring Trends in
 Burn Severity database (www.mtbs.gov) to
 generate two wildfire predictors: 1) a binary response of either burned
 (1) or unburned (0) indicating whether any severe fire occurred within a
 home range buffer and 2) a continuous value ranging from 0 to 1
 representing the proportion of the area that burned at high-severity
 during the Dixie Fire and North Complex Fire, providing a measure of the
 local extent of severe fire. We created a binary
 predictor to represent the lethal removal of barred owls based on the 57
 lethal removal locations of barred owls and barred owl x spotted owl
 hybrids from 2018 through 2022 that overlapped with our acoustic study
 area (Hofstadter et al. 2022). We first delineated a circular buffer of
 2,004 ha around each barred owl removal location, the approximate home
 range size of this species in the region (Wood et al. 2020a). Then, for
 flammulated and pygmy owls, we classified a site (i.e., ARU) as ‘1’ if
 their species-specific home range buffer overlapped with a barred owl home
 range buffer, and as ‘0’ if no such overlap occurred. For great horned
 owls, we designated sites as '1' if a great horned owl's
 home range buffer centered around a hexagon overlapped a barred owl's
 home range buffer. We coded sites as '0' if the great horned
 owl’s home range buffer did not overlap a barred owl's home range
 buffer.   <strong>Modeling
 site occupancy probability</strong> We developed
 Bayesian single-species “stacked” single-season occupancy models to assess
 the effect of barred owls and wildfire on site occupancy rates for the
 three forest owl species and controlled for imperfect detection (MacKenzie
 et al. 2018, Jones et al. 2020). We stacked the four-year encounter
 histories (2018, 2021, 2022, 2023) because we were more interested in
 regional owl occupancy patterns rather than turnover rates. We considered
 a year and site random effect on occupancy (see below) because the stacked
 model assumes independence among sites and therefore may underestimate
 error in model coefficients. We also chose this data format to increase
 the effective sample size, thereby enabling us to fit a larger set of
 occupancy predictors for less prevalent species and avoid difficulties
 associated with poor model fit.  We first modeled the
 potential effects of five predictors as fixed effects on the probability
 of detecting each species. We considered ordinal date as a predictor to
 accommodate potential changes in vocalization patterns over the breeding
 season (Crozier et al. 2003). We also considered the year as a predictor
 to assess interannual variation in detectability and account for different
 microphones in 2018. We included terrain ruggedness as a predictor because
 owl vocalizations may travel a shorter distance in areas of complex
 topography. We calculated terrain ruggedness by averaging the absolute
 differences between the elevation of a 10 m raster grid cell and the
 values of its eight surrounding cells within a 390 m buffer around each
 ARU location, where elevation data were obtained from a digital elevation
 model (Hijmans et al. 2015). We treated the number of hours ARUs recorded
 per secondary sampling period as a predictor, given that units were
 deployed for variable durations. We used the barred owl
 detection/non-detection history obtained from our passive acoustic surveys
 as a predictor to account for potential changes in vocalization behavior
 in response to barred owl presence (Wood et al. 2020b, Rugg et al. 2023).
 We coded hexagonal cells as ‘1’ if a barred owl was detected at least once
 in a given year at either ARU or ‘0’ if no detection occurred during that
 year. In 2018, barred owls were acoustically detected at 31 hexagonal
 cells, 18 of which overlapped with follow-up barred owl lethal removal
 locations. We did not include background noise as a detection covariate to
 maintain a minimal set of covariates and because BirdNET is robust to high
 ambient noise (Kahl et al. 2021). No pairs of predictors were highly
 correlated (all Pearson’s correlation coefficients &lt; 0.6). We
 incorporated all five predictors into the occupancy models for each
 species to minimize the potential number of models (see below).
 Next, we built six individual occupancy models to assess whether
 barred owl removals, severe wildfires, or interactions between the two
 predictors affected owl occupancy rates. We treated unburned and barred
 owl non-removal sites as control sites for each owl species and employed
 an unbalanced Before-After-Control-Impact design.  To
 test for severe fire effects independent of time since burn, we created
 the Fire Time-Constant Model:                     
 logit(Ψi) = β0 + β1burn + β2allpostburn + μ1[1|year] +
 μ2[1|site]  where burn was an indicator predictor for
 burned sites (i.e., burned = 1; unburned = 0), allpostburn was a
 continuous predictor representing the proportion of the area that burned
 at high-severity at any point post-fire, year was a random effect
 indicator for survey year, and site was a random effect for individual
 sites (i.e., ARUs for flammulated and pygmy owls and hexagonal cells for
 great horned owls).  To test for severe fire effects
 while accounting for different years post-burn, we created the Fire
 Time-Dependent Model:                       logit(Ψi) =
 β0 + β1burn + β2burn1 + β3burn2 + β4burn3 + μ1[1|year] +
 μ2[1|site]   where burn was the same
 as in the Fire Time-Constant Model and burn1, burn2, and burn3 were
 continuous predictors representing the proportion of the area that burned
 at high-severity relative to the number of years post-fire (i.e., all
 pre-fire sites and unburned sites post-fire = 0). To
 test for barred owl effects independent of time since removal, we
 developed the Barred Owl Time-Constant Model:         
                 logit(Ψi) = β0 + β1site-type + β2 allpostlethal +
 μ1[1|year] + μ2[1|site] where site-type was an
 indicator predictor for barred owl removal sites (i.e., removals = 1; no
 removals = 0 based on the criteria described above), allpostlethal was an
 indicator predictor for removal sites at any point post-removal, year was
 a random effect indicator for survey year, and site was a random effect
 for individual sites (i.e., ARUs for flammulated and pygmy owls and
 hexagonal cells for great horned owls).  To test for
 barred owl effects while accounting for years since removal, we built the
 Barred Owl Time-Dependent Model:                       
   logit(Ψi) = β0 + β1site-type + β2 lethal1 + β3lethal2 + β4lethal3^+^ +
 μ1[1|year] + μ2[1|site] where site-type was the same as
 in Barred Owl Time-Constant Model and lethal1, lethal2, and lethal3+ were
 indicator predictors for removal sites one, two, and three or more years
 post-removal (i.e., removal sites post-removal = 1; all pre-removal sites
 and non-removal sites post-removal = 0). We evaluated for the effects of
 barred owl lethal removal parameters by comparing native owl occupancy at
 removal sites to occupancy at those same sites in the years following
 removal. Since lethal removals happened in different years, a direct
 before-and-after comparison between control and impact sites was not
 possible. To address this, we used a contrast analysis to measure overlap
 between posterior distributions at removal sites (i.e.,
 site-type<sub>lethal</sub>) and post-removal sites (i.e.,
 lethal1<sub>lethal</sub>,
 lethal2<sub>lethal</sub>,
 lethal3+<sub>lethal</sub>, and
 allpostlethal<sub>lethal</sub>) (Pastore and Calcagnì 2019).
 We considered differences between posterior distributions to be
 biologically meaningful if their overlap was less than 15%.
 To test for interactions between fire and barred owl removals
 independent of time since burn, we created the Interaction Time-Constant
 Model:                       logit(Ψi) = β0 +
 β1site-type + β2burn + β3allpostlethal*burnYrs + μ1[1|year] +
 μ2[1|site] where site-type and burn were indicator
 predictors for barred owl removal and burned sites (i.e., removal or
 burned site = 1; no removals or unburned site = 0 based on the criteria
 described above) and allpostlethal*burnYrs indicated sites where barred
 owl removals preceded wildfire, measured for all years post-fire (i.e.,
 burned removal site post-burn = 1; all other sites = 0).
 Finally, to test for interactions between fire and barred owl
 removals while accounting for time since burn, we developed the
 Interaction Time-Dependent Model:                    
  logit(Ψi) = β0 + β1site-type + β2burn + β3allpostlethal<em>burn1 +
 β4allpostlethal</em>burn2 + β5allpostlethal*burn3
                                     + μ1[1|year] +
 μ2[1|site] where site-type and burn were the same as in
 Interaction Time-Constant Model and allpostlethal<em>burn1,
 allpostlethal</em>burn2, and allpostlethal<em>burn3 were
 indicator predictors for sites with interactions between barred owl
 removals and severe fire depending on the years since the fire occurred
 (i.e., allpostlethal</em>burn3 = 1 indicates a site where a barred
 owl removal was followed by high-severity fire in the third year after the
 fire, with the covariate representing a binary indicator rather than a
 continuous fire severity measure). We performed occupancy analyses using
 the R package “ubms” (Kellner et al. 2022) and used its predict function
 to estimate and visualize occupancy probabilities.
description: # Data from: Divergent responses of native predators to severe wildfire
 and biological invasion are mediated by life history Dataset DOI:
 [10.5061/dryad.c866t1gkt](10.5061/dryad.c866t1gkt) ## Description of the
 data and file structure This dataset contains detection histories and
 associated covariates for single-species "stacked" single-season
 occupancy models of forest owls in the northern Sierra Nevada, California.
 Data were collected during passive acoustic monitoring surveys conducted
 in 2018, 2021, 2022, and 2023 at multiple fixed survey sites. The surveys
 targeted three focal owl species: great horned owl (*Bubo virginianus*,
 GHO), flammulated owl (*Psiloscops flammeolus*, FLM), and northern pygmy
 owl (*Glaucidium gnoma*, NPO). Each dataset (data_stacked) combines weekly
 detection/non-detection records from autonomous recording units (ARUs)
 with site-level and survey-level covariates describing topography, barred
 owl presence/removal treatments, wildfire history, and survey effort. The
 structure supports Bayesian occupancy modeling in R. ### Files and
 variables #### File: barry_etal_owl_occupancy.R **Description:** R script
 to fit a Bayesian single-species "stacked" single-season
 occupancy model using the data. #### Files: FLM_data_stacked.csv,
 NPO_data_stacked.csv, GHO_data_stacked.csv **Description:** Encounter
 histories and covariates for flammulated owls (FLM), northern pygmy owls
 (NPO), and great horned owls (GHO) ##### Variables * v1 -v11: Weekly
 detection histories for each site during the breeding season (0 = not
 detected, 1 = detected, NA = not surveyed). * sp: Species code (FLM, NPO,
 GHO). * Site: Site identifier. * topo_tr: Terrain ruggedness index. *
 BOpresence: Barred owl detected during the season (1 = yes, 0 = no). *
 RemovalTreat: Site within a barred owl lethal removal treatment area (1 =
 yes, 0 = no). * FireTreat: Site within a high-severity wildfire footprint
 (1 = yes, 0 = no). * pSevYrs: Proportion of site area burned at high
 severity at any point post-fire. * Year: Survey year code (0 = 2018; 1 =
 2021; 2 = 2022; 3 = 2023). * SeasonBA: Pre/post-barred owl removal period
 (0 = pre-removal, 1 = post-removal). * PostRemovalImpact: Binary indicator
 for whether the site was in the post-removal period for barred owls (1 =
 post-removal, 0 = pre-removal). * B1-B3: Binary indicators for whether the
 site was in year 1, year 2, or year 3 following barred owl removal. *
 pSev1-pSev3: Proportion of the site area burned at high severity relative
 to the number of years post-fire (i.e., all pre-fire sites and unburned
 sites post-fire = 0). * PostRemovalpSev1-PostRemovalpSev3: Indicator
 predictors for sites with interactions between barred owl removals and
 severe fire depending on the years since the fire occurred (i.e.,
 allpostlethal*burn3 = 1 indicates a site where a barred owl removal was
 followed by high-severity fire in the third year after the fire). *
 PostRemovalpSevYrs: Sites where barred owl removals preceded wildfire,
 measured for all years post-fire (i.e., burned removal site post-burn = 1;
 all other sites = 0). * date.1.1-date.1.11: Scaled ordinal date covariates
 for each week of the season. * ARUhours1-ARUhours11: Total recording hours
 for each secondary sampling period, used to account for unequal sampling
 effort. ## Code/software Analyses were run in **R version 4.2.2
 (2022-10-31)**. Core packages: * unmarked – data structures
 (unmarkedFrameOccu) for occupancy models - ubms – Bayesian occupancy
 models via Stan (stan_occu) * rstan – Stan backend for MCMC sampling -
 dplyr – data wrangling
 (filter, rename, group_by, summarise, mutate, arrange) * ggplot2 –
 plotting predictions and effects - overlapping – posterior overlap metrics
 (e.g., overlap between groups) ## Access information Other publicly
 accessible locations of the data: * This dataset is available exclusively
 through Dryad and is not archived in other public repositories. Data was
 derived from the following sources: * Wildfire perimeter and severity data
 were obtained from the Monitoring Trends in Burn Severity (MTBS) program
 ([https://www.mtbs.gov](https://www.mtbs.gov)), which is publicly
 available under U.S. federal open data policy. - Barred owl removal
 treatment locations and survey data were collected as part of the barred
 owl removal experiment described in Hofstadter et al. (2022). * Terrain
 ruggedness data were derived from U.S. Geological Survey (USGS) digital
 elevation models (public domain).