Temporal trends

We modeled temporal variation in Spotted Owl vital rates in a variety of ways, including annual time effects (t), linear effects (T), log-linear effects (lnT), and spline models (SPLINE). We used spline models rather than quadratic time effect models because spline models could capture the quadratic effect were it to be present, but also allowed more overall flexibility for complex temporal patterns over long time periods (Bonner et al. 2009). In particular, spline models could accommodate a slowly varying external influence (e.g., regional climatic pattern) or a single “blip” in the middle of the time series. If we observed strong support for a spline model across multiple study areas, this might indicate that we were omitting some fundamental covariate. We used a penalized cubic spline model with knots every 5 yr going backward from the year with the last estimable parameter and with the first interval ≥5 years. For the fecundity analyses only, we also investigated a temporal even–odd year effect (EO), in which years of high reproductive output alternated with years of low reproductive output.

Barred Owl covariate

We did not specifically survey Barred Owls along with Spotted Owls. However, Barred Owls frequently responded when we used playback calls or vocal imitations to locate Spotted Owls, and we recorded all such detections. Based on a calling experiment conducted by Wiens et al. (2011), we estimated that cumulative annual detection probabilities of Barred Owls were >85% at territories in which we conducted ≥3 nocturnal surveys for Spotted Owls. Based on these detection data, we created a Barred Owl covariate that was specific to year and study area and reflected the proportion of Northern Spotted Owl territories (i.e. Thiessen polygons) in which Barred Owls were detected ≥1 time per year in each study area (Anthony et al. 2006, Forsman et al. 2011; Appendix C Figure 13). We used this covariate to model the Barred Owl effect on fecundity and apparent survival in individual study areas and in the meta-analysis. For the occupancy analysis, we used detections of Barred Owls from the multiple Spotted Owl surveys conducted each year to estimate the effects of each species on the other, including probability of detection and rates of colonization and extinction (Bailey et al. 2009, Yackulic et al. 2014).

Habitat covariates

Habitat covariates were year-specific covariates applied at the study area scale for survival and fecundity analyses (Figure 2B). For occupancy analyses in individual study areas, habitat covariates were developed as year- and territory-specific covariates for individual owl territories (i.e. Thiessen polygons; Figure 2A). Covariates developed to represent the amount and spatial distribution of Northern Spotted Owl habitat within study areas included: (1) the proportion of Northern Spotted Owl nesting and roosting habitat in each study area each year (HAB), (2) habitat disturbance, or the mean percentage of nesting and roosting habitat lost during 3-yr intervals prior to and including each survey year (HD), (3) a neighborhood focal statistic defined as the proportion of 30 × 30 m pixels in each study area with ≥50% nesting and roosting habitat within 800 m of each pixel (CORE), and (4) the proportion of edge habitat within each study area (EDGE), with edge as the amount of nesting and roosting habitat within 100 m of all other cover types (Appendix D Figure 14). We hypothesized that higher proportions of HAB and CORE would have positive effects on Spotted Owl vital rates, and that EDGE would have a variable effect, depending on the amount of nesting and roosting habitat remaining (EDGE and HAB interaction) and differences in the dominant types of prey available in each study area (Franklin et al. 2000). We predicted that study areas with more habitat disturbance (i.e. more habitat loss) would have lower survival rates than areas experiencing less habitat disturbance. These covariates were developed for all analyses as follows:

Weather and climate

We used a variety of covariates to investigate any possible effects of weather and climate on vital rates of Northern Spotted Owls (Table 1). All weather and climate covariates were time-specific, linear effects applied at the scale of individual study areas. These variables included measures of seasonal and annual weather, as well as long-term climatic conditions. Specific covariates included total precipitation and mean minimum temperature during various life-history stages, the Palmer Drought Severity Index (PDSI), the Southern Oscillation Index (SOI), and the Pacific Decadal Oscillation (PDO; Franklin et al. 2000, Seamans et al. 2002, LaHaye et al. 2004, Olson et al. 2004, Dugger et al. 2005, Glenn 2010, 2011a, 2011b, Forsman et al. 2011). Mean temperature and precipitation data were obtained from PRISM (Parameter Elevated Regression on Independent Slope Models) maps of each study area (PRISM Climate Group,  http://prism.oregonstate.edu). PRISM maps were raster-based digital maps with 4-km² resolution for mean monthly temperature (minimum and maximum; °C) and precipitation (cm), developed from weather station data and a digital elevation model (Daly 2006). From the mean monthly PRISM maps we calculated total precipitation and mean minimum monthly temperature for seasons that corresponded with important life-history stages of Spotted Owls as follows: winter (W; November 1–February 28), early nesting season (EN; March 1–April 30), late nesting season (LN; May 1–June 30), and the entire annual cycle (A; July 1–June 30). Temperature and precipitation values for each study area and time period were obtained by computing the average values of PRISM raster cells that fell within the study area boundaries.

The PDSI is a measurement of moisture conditions standardized for comparison across regions ( http://www.ncdc.noaa.gov/temp-and-precip/climatological-rankings/). The PDSI is calculated using precipitation, temperature, and soil moisture data, which allows the derivation of the basic components of the water balance, including evapotranspiration, soil recharge, runoff, and moisture loss from the surface layer (Alley 1984). We considered the PDSI an index of primary productivity that had the potential to influence the abundance of Spotted Owl prey. Values ranged from −6 (extremely dry) to +6 (extremely wet), with 0 representing normal water balance conditions. The PDSI was calculated separately for each climatic region in Washington, Oregon, and California. Most study areas fell within a single climatic region. When study areas included multiple climatic regions, we used a weighted average PDSI based on the proportion of the study area that fell within each region. We also averaged monthly values of the SOI ( http://www.cpc.ncep.noaa.gov/data/indices/soi), and the PDO ( http://jisao.washington.edu/pdo/PDO.latest) to generate annual measures (July 1 to June 30) that reflected region-wide climatic patterns that affected all study areas.

Land ownership, region, latitude, and prey species richness

In our meta-analyses, we evaluated whether Spotted Owl vital rates or rates of population change varied in relation to land ownership, region, latitude, and prey species richness. Land ownership (OWN) was a categorical variable that divided the 11 study areas into 3 ownership categories depending on whether primary ownership was federal, private, or a relatively equal mix of federal and private (Table 3). The region (REG) covariate classified each study area into 1 of 6 geographic regions based on state boundaries and major vegetation types (Table 3). Latitude (LAT) was a continuous variable measured at the center of each study area. The prey diversity index (PREY) was a discrete variable that characterized the maximum number of potential mammalian prey species (range: 6–17) that were available to Spotted Owls in each study area. We estimated the PREY covariate by using extensive data on the diets of Spotted Owls (Cutler and Hays 1991, Carey et al. 1992, Zabel et al. 1995, Ward et al. 1998, Forsman et al. 2001, 2004, Rosenberg et al. 2003) and species distribution maps in NatureServe ( http://www.natureserve.org/conservation-tools/data-maps-tools/digital-distribution-maps-mammals-western-hemisphere), summarized across the range of the Northern Spotted Owl in 50-km hexagons. Long-term data on prey abundance were not available for any of the Spotted Owl study areas, so the PREY covariate was a simple attempt to address variation in prey species richness among study areas.

Reproduction covariate

We used a covariate that was specific to year and study area, mean number of young fledged per female (NYF) in the current year t, to test whether reproductive effort affected adult survival in the interval between year t and year t + 1 (Franklin et al. 1996, Anthony et al. 2006, Forsman et al. 2011; Appendix E Figure 15). We also investigated the effect of reproduction on detection probability in year t, because breeding birds are generally easier to detect than nonbreeders (Anthony et al. 2006, Forsman et al. 2011, Stoelting et al. 2015).

Barred Owl removal study

Beginning in 2009, a paired before–after control–impact (BACI) study design was implemented in the GDR study area, where lethal removal of Barred Owls was the treatment effect on Northern Spotted Owl vital rates (Diller et al. 2014). The GDR demographic study area was partitioned into treatment (Barred Owls lethally removed) and control (Barred Owls undisturbed) areas to estimate the response of Spotted Owl fecundity, survival, and rate of population change to the removal activities. To account for geographical differences in the history of timber harvesting, physiographical patterns, and density of Barred Owl and Spotted Owl territories, the GDR study area was divided into 3 paired treatment and control areas totaling 84,205 and 72,711 ha, respectively. Within these treatment areas, investigators attempted to remove all Barred Owls detected (Diller et al. 2014). For analyses involving individual study areas, a BACI (e.g., Stewart-Oaten et al. 1986) design was incorporated for the GDR study area, with parameters estimated separately for treatment and control areas both before and after removals began, unless otherwise noted. For the meta-analyses conducted with data from all study areas combined, only data from control and treatment areas prior to the Barred Owl removals (up to 2008) and control areas after removals began (2009–2013) were included, so that the GDR data were comparable with data from the other study areas.

Annual rate of population change (λ)

We estimated annual rates of population change (λ) in individual study areas using the λ- and recruitment- (f) parameterizations of the temporal symmetry models (Pradel 1996) implemented in program MARK. For this analysis we used all banded, territorial birds (S1, S2, adults) combined into a single age class. We ran 5 random effects models on λ_t, including the intercept-only (no effect), general time (t), linear time trend (T), log-linear time trend (lnT), and spline models (SPLINE; with knots every 5 yr backward from 2013, such that the first interval was ≥5 yr; Bonner et al. 2009). We dealt with expansions or contractions in areas surveyed (Table 2) using the design matrix in program MARK, such that all estimates of λ reflected changes in owl numbers and were not confounded with sampling changes (Anthony et al. 2006, Forsman et al. 2011; see Appendix H for further details). Start years varied by study area, but estimates were generated from the start date through 2012 in all cases (Table 2). However, for general models with time-specific capture and survival probabilities, the first and last estimates (λ_1, λ_k−1) were confounded with other parameters, and the second estimate was frequently biased (λ₂; Hines and Nichols 2002). Thus, we only present estimates from 2 yr after the start date through 2011, and used these estimates in random effects models.

Estimates of realized population change

We estimated realized population change (ƊĮ_t), which portrays the population trajectory (Δ_t = N_t/N_x) in each year of the study (N_t) relative to population size in the first year (N_x) that λ_t was estimated (Franklin et al. 2004).

Annual estimates of λ̂_t were based on the full fixed effects model [φ(t) p(t) f(t); i.e. time-dependent (t) survival (φ), capture probability (p), and recruitment (f)], and annual estimates of realized population change (ƊĮ_t) were computed as:

We estimated 95% confidence limits (CL) for ƊĮ_t using a parametric bootstrap algorithm (Franklin et al. 2004). Estimates of annual survival (Φ̂_t), recruitment (fĮ_t), and recapture probabilities (pĮ_t) from the full fixed effects model [φ(t) p(t) f(t)], and an estimate of initial abundance (NĮ_x), were used to stochastically generate 1,000 sets of individual capture histories. These simulated capture history datasets were analyzed to obtain 1,000 estimates of λ_t and Δ_t, and these estimates were used to generate empirical confidence intervals based on the i^th and j^th values of Δ_t arranged in ascending order, where i = 25 (0.025*1,000) and j = 975 (0.975*1,000).

Meta-analysis of annual rate of population change

We conducted the meta-analysis of the annual finite rate of population change using the same data that we used to estimate λ_t for individual study areas. However, we only used data for 1992–2013 so that we could make inferences based on the same years for all study areas. As in the analysis of individual areas, we used ĉ = 1 for modeling all rates of population change (Anthony et al. 2006, Forsman et al. 2011; see Appendix H). We used the global model [φ(g*t) p(g*t) f(g*t)] as the basis for random effects modeling of covariate effects on recruitment, where g indicated individual study areas. We only included climate and habitat covariates that we predicted would have effects on recruitment. In some cases we modeled 1- or 2-yr lag effects of climate covariates, because we hypothesized that this was likely the most appropriate relationship if climate was associated with annual reproductive output (NYF) and it took fledged young at least 1 or 2 yr to be recruited into the territorial population.

Territory occupancy modeling

We investigated the co-occurrence dynamics of Northern Spotted Owls and Barred Owls based on 19 yr of detection data for both species (1995–2013) in 10 study areas, and 15 yr of detection data (1999–2013) in the GDR study area. We excluded data from the Barred Owl removal treatment areas in the GDR study area after 2008. We created detection histories that consisted of a sequence of detections (1) and nondetections (0) for both species within and among years on all study areas. We applied these data to the multiseason (robust design) extension of the conditional, 2-species occupancy model (MacKenzie et al. 2004, 2006) following Miller et al. (2012) and Yackulic et al. (2014), and used program MARK to estimate occupancy parameters and model selection results. Model parameters included initial occupancy (Ψ₁), colonization (γ_i), extinction (ϵ_i), and detection probabilities (p_ij) for both species as potential functions of presence of the other species. For initial occupancy, we used the parameterization of Richmond et al. (2010) because it is more stable than the parameterization of the original 2-species models developed by MacKenzie et al. (2004, 2006), which can fail to converge when covariates are included. Based on prior research (Van Lanen et al. 2011, Wiens et al. 2014), we assumed that the Barred Owl was the dominant species (coded as “A”) and that the Northern Spotted Owl was the subordinate species (coded as “B”). Although occupancy dynamics parameters for both Spotted Owls and Barred Owls were generated in this analysis, here we focus on the patterns of occupancy and occupancy dynamics (extinction and colonization rates) for Spotted Owls only, in relation to the presence or absence of Barred Owls. The specific parameters of interest were: (1) initial probability of occupancy of Northern Spotted Owls when Barred Owls were absent () and when Barred Owls were present (), (2) the probability that a territory unoccupied by a Spotted Owl in year i was occupied by a Spotted Owl in year i + 1 (i.e. colonization) when Barred Owls were present () and when Barred Owls were absent (), (3) the probability that a territory occupied by a Spotted Owl in year i was unoccupied in year i + 1 (i.e. local extinction) when Barred Owls were present () and when Barred Owls were absent (), and (4) annual probability of territory occupancy by Northern Spotted Owls when Barred Owls were present () and when Barred Owls were absent (), which was derived using the best model structure for detection, extinction, and colonization rates.

We analyzed each study area separately using fixed effects models and an iterative model selection process identified a priori (see Appendix H for details). We modeled colonization and extinction rates for both species with linear time trends (T), year-specific effects (t), no temporal effects (intercept-only), the presence vs. absence of the other species, and the effects of habitat covariates. Finally, the effects of 2- and 3-yr lags in Spotted Owl annual reproductive output were also modeled for Spotted Owl colonization probabilities (see Yackulic et al. 2014 for details on this general approach).

Fecundity

Analyses for individual study areas were conducted on the number of young produced per territorial female (NYF), but our results are presented as fecundity, defined as the number of female young produced per territorial female per year. This was calculated as NYF/2, because the sex ratio of juvenile owls at hatching is approximately 1:1 (Fleming et al. 1996). Spotted Owls are strongly territorial, have high site fidelity, and are detectable even when they are not breeding (Franklin et al. 1996, Reid et al. 1999). Thus, we assumed that sampling over the course of an entire breeding season was not biased toward birds that reproduced, and that the sample of owls used in our analyses was representative of the territorial population. Owls that were recruited into the banded population were assigned to 1 of 3 discrete age classes based on their age at first capture as a territorial bird (S1 = 1 yr old, S2 = 2 yr old, Adult = ≥3 yr old; Table 2). We determined age classes based on known age of birds first banded as juveniles, or plumage attributes of birds first banded as nonjuveniles (Forsman 1981, Moen et al. 1991, Franklin et al. 1996).

Mean annual NYF was computed by age class and then averaged across years for estimates of age-specific reproductive output. Standard errors were calculated as the standard errors of the averages among years, which gave equal weight to all years regardless of the number of owls sampled (Anthony et al. 2006, Forsman et al. 2011). This approach essentially treated year as a random effect, with year effects being large relative to within-year sampling variation.

We developed an a priori model set and used a linear mixed model approach implemented with PROC MIXED in SAS (SAS Institute 2008) to investigate patterns of variation and hypothesized relationships between covariates and NYF (see Appendix H for details). Models included the effects of age (S1, S2, Adult), annual time variation (t), linear or quadratic time trends (T, TT), an autoregressive time effect (AR1), the Barred Owl (BO) covariate, a temporal even–odd year effect (EO) in which years of high reproductive output alternated with years of low reproductive output, and the weather, climate, and habitat covariates described previously.

Meta-analysis of fecundity

We restricted the meta-analysis of fecundity to adult females because the sample size of younger age classes was small (<10%), particularly in the most recent years of study. We used the same covariates as in the individual study area analyses to generate an a priori model set, with the addition of models investigating the effects of latitude (LAT), region (REG), land ownership (OWN), and prey species richness (PREY) as fixed random variables. We used mixed models to analyze mean NYF per year, and treated sampling units (study areas within years) as random effects (Anthony et al. 2006, Forsman et al. 2011).

Apparent survival

We used capture–recapture (resighting) data and Cormack-Jolly-Seber open population models (Lebreton et al. 1992) to estimate recapture probabilities (p) and annual apparent survival probabilities (φ) of nonjuvenile, territorial owls (Table 2). Annual estimates of survival corresponded roughly to the interval from June 15 in year t to June 14 in year t + 1, which reflected the approximate midpoint of the annual field season during which demographic (mark–resighting) data were collected (March–August). Estimates and model selection results to investigate the effects of Barred Owls, reproduction, habitat, weather, climate, and time effects on apparent survival of Spotted Owls were generated using the Method of Moments random effects module in program MARK (White and Burnham 1999; Appendix F, H).

Meta-analysis of apparent survival

We used the same general protocol for the meta-analysis of apparent survival as for the analysis of apparent survival in individual study areas (for details see Appendix H). We ran random effects models in program MARK (White et al. 2001) to investigate the effect of covariates (i.e. time, Barred Owls, cost of reproduction, weather, climate, habitat, latitude, region, and prey species richness), always excluding the last confounded estimate of survival (φ_K; Burnham and White 2002, Burnham 2013; Appendix F).

Individual study areas

We estimated annual rates of population change (λ) using capture histories for 5,992 territorial owls from all age classes (S1, S2, Adult; Table 2). We used a base model for random effects modeling that included general time effects on survival [φ(t)] and lambda (λ(t); Table 4) in all study areas except GDR, where the recruitment [f(t)] and lambda [λ(t)] parameterization was used to facilitate convergence. The best fixed effects structure for capture rates included annual time effects [p(t)] in 4 study areas (RAI, COA, CAS, and HUP), additive effects of sex and annual time [p(sex + t)] in 6 areas (CLE, HJA, TYE, KLA, NWC, and GDR; capture rates higher for males), and an interaction between sex and time [p(sex*t)] in 1 area (OLY).

TABLE 4.

Estimates, standard errors (SE), and lower (LCL) and upper (UCL) 95% confidence limits of mean annual rate of population change based on a reverse Jolly-Seber model (λ̂_RJS) and temporal process variance (Σ̂_temporal) for Northern Spotted Owls in 11 study areas in Washington, Oregon, and California, USA, 1985–2013. Estimates of λ̂_RJS were generated using the intercept-only random effects model [RE(.)]. Estimates of temporal process variance were based on the best random effects models using time-specific estimates of survival (φ), capture probability (p), and rate of population change (λ) or recruitment (f).

The best random effects model for 7 of the 11 study areas included a negative linear time trend on λ (RE(T); Table 4), with 95% CIs of covariate coefficients widely overlapping 0 only for CLE and OLY (Figure 3), suggesting that annual rates of decline were increasing over time in many areas. The spline model [RE(SPLINE)] performed best for HJA and NWC, although the linear time trend model was also competitive for NWC, and the coefficient was negative with a 95% CI that did not overlap zero. The intercept-only model [RE(.)] received the most support for RAI and CAS.

FIGURE 3.

Estimates of covariate coefficients and 95% confidence intervals for the linear time trend (T) on λ from the best random effects model containing a linear time trend, for Northern Spotted Owls in 11 study areas in Washington, Oregon, and California, USA. The estimate for the GDR study area represents the parallel (additive) trend on control and treatment areas before Barred Owl removals began (1990–2008). See Table 2 for study area abbreviations.

Mean estimates of λ from the RE(.) models suggested declining population trends (i.e. λ̂ < 1.0) in almost all study areas, with strong evidence of declines in CLE, RAI, OLY, COA, HJA, NWC, HUP, GDR-CB, GDR-TB, and GDR-CA, and less evidence of declines in TYE, KLA, and CAS (Figure 4). The only estimate of λ that suggested an increasing population was observed in GDR treatment areas after Barred Owl removals began in 2009 (GDR-TA; λ̂ = 1.03), although the 95% CI widely overlapped 1.0. Estimated annual rates of decline were variable (Table 4), but were lowest in the GDR control areas before Barred Owl removals began in treatment areas in 2009 (1.2% annual decline), and highest in the CLE study area in Washington (8.4% annual decline) and in GDR control areas after 2009 (12.0% annual decline). The weighted mean estimate of λ for all study areas (excluding GDR-TB and GDR-TA) was 0.962 ± 0.019 (95% CI: 0.925 to 0.999), indicating an estimated decline of 3.8% per year across the range of the Northern Spotted Owl.

FIGURE 4.

Estimated mean rates of population change (λ̂) and 95% confidence intervals for Northern Spotted Owls from the random effects (RE) intercept-only model [RE λ(.)] in each of 11 study areas in Washington, Oregon, and California, USA, 1985–2013. Estimates for the GDR study area are presented separately for control and treatment areas before (1990–2008) and after (2009–2013) Barred Owls were removed (GDR-CB = control before removal, GDR-TB = treatment before removal, GDR-CA = control after removal, GDR-TA = treatment after removal). See Table 2 for study area abbreviations.

Realized population change

Our estimates of realized population change indicated that populations in Washington declined by 55–77% (Figure 5A). Declines in Oregon were more variable, ranging from 31% in TYE to 68% in COA (Figure 5B), and in 2 cases (KLA and TYE) the 95% CIs for realized population change widely overlapped 1.0 for most or all of the last several years, indicating uncertainty about annual rates of population change for these areas. In California, declines ranged from 32% to 55%, except in the treatment areas of GDR (GDR-T), where the estimated overall population decline was only 9% (Figure 5C). Realized population change estimates for HUP and GDR-T included confidence limits that overlapped 1.0 in many years, indicating uncertainty about annual rates of population change in these areas.

FIGURE 5.

Annual estimates of realized population change (Δ_t) with 95% confidence intervals for Northern Spotted Owls at (A) 3 study areas in Washington, (B) 5 study areas in Oregon, and (C) 3 study areas in California, USA. Estimates for the GDR study area are presented separately for control and treatment areas in relation to Barred Owl removals beginning in 2009. See Table 2 for study area abbreviations.

FIGURE 5.

Continued.

Meta-analysis of annual rate of population change: recruitment

As described above, we focused the meta-analysis of population change on recruitment parameterization, in which λ was written as the sum of apparent survival and recruitment rate. The fixed effects model on which all random effects models for recruitment were based included area by year interactions on survival, capture probabilities, and recruitment [φ(Area*t) p(Area*t) f(Area*t)]. The best random effects model had nearly all of the model weight and included interactions between total winter precipitation (WP) and mean winter minimum temperature (WMT; Table 5). As predicted, recruitment was negatively affected by WP as a main effect (β̂ = −0.0003 ± 0.0001 SE, 95% CI: −0.0005 to −0.0001). However, contrary to our prediction, the main effect of WMT was also negatively related to recruitment rates (β̂ = −0.0066 ± 0.0037 SE, 95% CI: −0.0140 to 0.0006), although the interaction between precipitation and temperature was positive (β̂ = 0.0001 ± 0.0000 SE, 95% CI: −0.0000 to 0.0001). Thus, recruitment was highest when both total precipitation (29 cm) and mean minimum winter temperature (−9.5°C) were lowest, and higher WP resulted in lower recruitment rates when WMT was low (<4°C; Figure 6). When the total precipitation level was near its average for all years and study areas (112 cm), recruitment rates were nearly constant across the range of mean minimum temperatures.

TABLE 5.

Model selection results from the meta-analysis of the finite rate of population change (λ) of adult Northern Spotted Owls in 11 demographic study areas in Washington, Oregon, and California, USA, 1994–2013. Random effects models (RE) of recruitment were run using a general fixed effects base model [φ(Area*t) p(Area*t) f(Area*t)] that included the interaction between study area (Area) and general time effects (t) on apparent survival (φ), capture probability (p), and recruitment (f). Models were ranked according to Akaike's Information Criterion adjusted for small sample size (AIC_c). The model deviance (Deviance), number of parameters (K), difference in AIC_c between each model and the model with the lowest AIC_c (ΔAIC_c), and Akaike model weights (w_i) are given for all models.

FIGURE 6.

Predicted estimates of recruitment of Northern Spotted Owls from the best random-effects (RE) model from the meta-analysis of lambda, using the survival (φ), recruitment (f), and capture probability (p) parameterization with study area (Area) and general time (t) fixed effects [φ(Area*t) p(Area*t) f(Area*t) RE f(WP*WMT)]. Estimates of recruitment are plotted across the range of mean minimum winter temperatures (WMT) from the data, for the minimum (29 cm), mean (112 cm), and maximum (297 cm) levels of total winter precipitation (WP) across all study areas and years.

Individual study areas

We estimated fecundity using 12,969 records in which we determined the number of young produced by territorial females of known age, 91% of which were females ≥3 years old (i.e. adults; Table 9). Female age was an important factor affecting fecundity in all study areas, with mean fecundity generally lowest for 1-yr-olds, intermediate for 2-yr-olds, and highest for adults. In most study areas, the mean annual fecundity of adult females was between 0.18 and 0.34. The one notable exception was the CLE study area in Washington, where mean annual fecundity for adult females (x̄ = 0.570 ± 0.045) was nearly twice as high as in any other area.

TABLE 9.

Estimates of means (x̄) with standard errors (SE) of age-specific fecundity (number of female young produced per female per year) of Northern Spotted Owls in 11 study areas in Washington, Oregon, and California, USA, 1985–2013.

Female age was strongly associated with NYF and occurred in top or competitive models for all study areas (Table 10). However, there was considerable model uncertainty for other effects. Of the 11 study areas, 5 had top models or competitive models that included negative linear (T) or quadratic (TT) time trends on fecundity, including 1 area in Washington (CLE), 3 areas in Oregon (COA, HJA, and TYE), and 1 area in California (NWC; Table 10, 11). The 95% CIs for the covariate coefficients (β̂) from these models excluded 0 for TYE and NWC, suggesting strong support for a declining trend in fecundity in these areas (Table 11). There was less support for declining trends in CLE, COA, and HJA, as 95% CIs slightly overlapped 0 for these study areas. Annual variation in fecundity was particularly high in study areas in Washington, which may have made it more difficult to detect trends in that region (Figure 9). For example, it was common for there to be years of no reproduction in the RAI and OLY study areas in Washington, whereas years with no reproduction were rare in study areas in Oregon, and were never observed in any of the California study areas.

TABLE 10.

Model selection results for models with ΔAIC_c < 2 from the analysis of mean age-specific number of young fledged per year per female (NYF) for Northern Spotted Owls in 11 study areas in Washington, Oregon, and California, USA, 1985–2013. Models were ranked according to Akaike's Information Criterion adjusted for small sample size (AIC_c), and the difference in AIC_c between each model and the model with the lowest AIC_c (ΔAIC_c), the number of parameters (K), Akaike weight (w_i), and model deviance (−2logL) are included for each candidate model.

TABLE 11.

Best model containing a linear (T), quadratic (TT), or autoregressive (AR(1)) time effect on the mean annual number of young fledged per adult female Northern Spotted Owl (NYF) in each of 11 study areas in Washington, Oregon, and California, USA, 1985–2013. Model covariate coefficients (β̂), standard errors, and 95% confidence limits (lower: LCL; upper: UCL) are presented for the linear and/or quadratic term in each model.

FIGURE 9.

Annual fluctuations in mean fecundity (number of female young fledged per female) of Northern Spotted Owls in (A) 3 study areas in Washington, (B) 5 study areas in Oregon, and (C) 3 study areas in California, USA. Mean fecundity was graphed separately for the areas within the Green Diamond (GDR) study area where Barred Owls were removed (2009–2013; GDR–Treatment) and where Barred Owls were not removed (1990–2013; GDR–Control). See Table 2 for study area abbreviations.

Barred Owl presence (BO) was included in a top or competing fecundity model in only 2 of the 11 study areas (COA and KLA; Table 10, 12) and the relationship was negative, with 95% CIs for the covariate coefficients not overlapping 0 in both cases. It is important to note that the proportion of each study area occupied by Barred Owls gradually increased over time (Appendix C Figure 13), so the temporal effect and BO effect were highly correlated and not easily separated. This may explain the lack of effect (or counterintuitive effect) of Barred Owls on fecundity in some areas.

TABLE 12.

Best model including the effect of Barred Owls (BO) on the mean annual number of young fledged per adult female Northern Spotted Owl (NYF) for each of 11 study areas in Washington, Oregon, and California, USA, 1985–2013. Model covariate coefficients (β̂), standard errors, and 95% confidence limits (lower: LCL; upper: UCL) are included for the Barred Owl effect in each model.

Habitat covariates (HAB and CORE) were included in the top or competitive models for fecundity in 7 study areas (CLE, COA, HJA, TYE, CAS, NWC, GDR, GDR-C, and GDR-T; Table 10, 13). More nesting and roosting habitat (HAB or CORE) was associated with higher NYF in these areas; however, only in TYE, NWC, GDR, and GDR‐C did the 95% CIs for covariate coefficients not overlap 0, suggesting that there was little support for strong associations between habitat and NYF in most study areas (Table 13).

TABLE 13.

Best model containing the effect of habitat on the mean annual number of young fledged per adult female Northern Spotted Owl (NYF) for each of 11 study areas in Washington, Oregon, and California, USA, 1985–2013. Model covariate coefficients (β̂), standard errors (SE), and 95% confidence limits (lower: LCL; upper: UCL) are included for the habitat effect in each model.

Climate covariates occurred in competitive models for 8 of 11 study areas (Table 10, 14), but the best covariate and the direction of the effect varied among areas (Table 14). For example, the effect of mean monthly minimum temperature during the early nesting season (ENT) occurred in the best model or a competitive model for 5 areas (CLE, RAI, HJA, CAS, and HUP), but the effect was positive in 3 areas (CLE, HJA, and CAS) and negative in 2 others (RAI and HUP; Table 14), with 95% CIs that did not overlap 0 for CLE, CAS, and HUP. Precipitation during the early nesting season (ENP) occurred in competitive models for only 2 study areas (RAI and HUP), and in both areas ENP interacted with ENT. The 95% CIs around the ENP covariate coefficients excluded or slightly overlapped 0 in both cases, and the association was negative (as predicted) in 1 area (HUP), but positive in the other (RAI). Precipitation during the late nesting season (LNP) was included in a best or competitive model for only 1 study area (KLA), where it occurred as an interaction with the presence of Barred Owls, and where the 95% CI for the effect of LNP did not overlap 0. Mean minimum winter temperature (WMT) and total winter precipitation (WP) were included in top or competitive models for only 2 study areas (NWC and GDR). In both cases, higher minimum temperatures during winter were associated with higher fecundity, and more winter precipitation was associated with lower fecundity, but all 95% CIs for the WMT and WP covariate coefficients overlapped 0, except for the effects of WMT in the Barred Owl treatment areas in the GDR study area. The Southern Oscillation Index (SOI) occurred in a competitive model for only 1 study area (OLY), where the 95% CI around the covariate coefficient slightly overlapped 0, suggesting weak support for this effect (Table 10, 14).

TABLE 14.

Best model containing the effect(s) of climate or weather on the mean annual number of young fledged per adult female Northern Spotted Owl (NYF) for each of 11 study areas in Washington, Oregon, and California, USA, 1985–2013. Model covariate coefficients (β̂), standard errors (SE), and 95% confidence limits (lower: LCL; upper: UCL) are included for the climate or weather effect(s) in each model.

Estimation of spatial (site-to-site), temporal (year-to-year), and residual variance in the territory-specific data indicated that the proportion of variance in NYF attributable to territories (spatial variance) was generally low (<6%; Table 15). The proportion of variance attributable to fluctuations over time (temporal variance) ranged from 5% to 20%, while the proportion of unexplained variation (residual variance) was generally high (>77%). As a consequence, the variation in NYF that was able to be explained by time and territory was overwhelmed by unexplained residual variation (i.e. other factors not considered).

TABLE 15.

Variance components of the mean annual number of young fledged per adult female Northern Spotted Owl (NYF) from the mixed-model analysis of year- and territory-specific effects in 11 study areas in Washington, Oregon, and California, USA, 1985–2013.

Meta-analysis of fecundity

The meta-analysis of fecundity produced 6 competitive models, all of which included the additive effects of region and time (Table 16). These models all suggested that fecundity varied by time, and was parallel across regions (Figure 10). A linear time trend (T) in fecundity was not supported because of the complex pattern in fecundity over time that resulted due to the dissipation of the even–odd year effect in most study areas after about 1999. Model weights were fairly evenly distributed among the 6 competitive models, but the 2 models with the highest model weights were a model that included the negative effect of core habitat (CORE; β̂ = −0.14, 95% CI: −0.25 to −0.02), and a model that included the additive negative effects of Barred Owls (BO; β̂ = −0.14, 95% CI: −0.30 to 0.01) and the amount of edge habitat (EDGE; β̂ = −0.60, 95% CI: −1.20 to 0.00; Table 16). Three other competitive models included 1 or more of these 3 covariates (CORE, BO, EDGE), and the amount of nesting and roosting habitat (HAB; β̂ = −0.23, 95% CI: −0.46 to 0.01) was also included in a competitive model. However, only for core habitat did the 95% CI for the covariate coefficient not overlap 0, although this relationship suggested that more core habitat was associated with decreased fecundity, contrary to predictions. None of the models that included the effects of land ownership, latitude, prey species richness, habitat disturbance, or climate were competitive, indicating that these covariates did not explain variation in fecundity across the range of the Northern Spotted Owl. The region effect on fecundity (Figure 10) appeared to be largely related to the pattern previously described in the individual study area analysis (Table 9, Figure 9), in which the east slope of the Cascades in Washington (region: Washington mixed conifer; 1 study area: CLE) had much higher fecundity than all other areas (Figure 9, 10).

TABLE 16.

Model selection results for models with ΔAIC_c < 5 from the meta-analysis of the mean annual number of young fledged per adult female Northern Spotted Owl (NYF) in 11 study areas in Washington, Oregon, and California, USA, 1985–2013. Models were ranked according to Akaike's Information Criterion adjusted for small sample size (AIC_c) and the difference in AIC_c between each model and the model with the lowest AIC_c (ΔAIC_c). The number of parameters (K), Akaike weight (w_i), and model deviance (−2logL) are included for each candidate model.

FIGURE 10.

Annual fluctuations in mean fecundity (number of female young fledged per female) of Northern Spotted Owls in 6 geographic regions from the meta-analysis of 11 study areas in Washington, Oregon, and California, USA, 1985–2013.

Individual study areas

We used the encounter histories of 5,090 owls (excluding 1-yr-olds [S1]; see Appendix H) to estimate apparent survival in individual study areas (Table 2, 17). The best fixed effects model that we used as the basis for random effects modeling included time effects on apparent survival and capture rates for 6 areas (φ(t) p(t); RAI, OLY, COA, TYE, CAS, and HUP), and a time effect on survival and an additive time and sex effect (males higher than females) on capture rates for 5 areas (φ(t) p(sex + t); CLE, HJA, KLA, NWC, and GDR). For the GDR study area the fixed effects model structure also included an interactive effect of the Barred Owl treatment effect and time on survival, and additive effects of Barred Owl treatment, sex, and time on recapture rates [φ(Trt*t) p(sex + t + Trt)]. Survival and captures rates in the GDR study area were higher in treatment areas after Barred Owl removals began in 2009, and males had higher recapture rates than females. The best random effects model varied widely by study area, with temporal variation best modeled by a trend or explained by a variety of factors, including Barred Owls, climate, and habitat (Appendix I Figure 16). Mean estimates of apparent survival ranged from a low of 0.804 ± 0.032 in GDR control areas after Barred Owl removals began in 2009, to a high of 0.870 ± 0.009 in the HJA study area and 0.870 ± 0.021 in the GDR treatment areas after Barred Owl removals began.

TABLE 17.

Estimates of mean apparent survival () and temporal process variance (Σ̂_temporal), with associated standard errors and 95% lower (LCL) and upper (UCL) confidence limits, for adult female Northern Spotted Owls in 11 study areas in Washington, Oregon, and California, USA, 1985–2013. Estimates of mean apparent survival were generated using the intercept-only random effects model [RE(.)]. Estimates of temporal process variance were based on the best random effects model (Mean, T, lnT, TT, or SPLINE) for each study area using time-specific estimates of apparent survival (φ) and capture probability (p).

Five of the 11 study areas included either a negative linear (T) or log-linear (lnT) time trend on survival in the best model (GDR) or competitive models (CLE, RAI, CAS, and HUP), but the effect was strong in only 1 area (RAI) as evidenced by a 95% CI for the covariate coefficient that did not overlap 0 (Table 18). However, in 9 of the 10 study areas where it was investigated, the Barred Owl covariate (BO), which exhibited increasing positive trends over time in all areas (Appendix C Figure 13), was included in the random effects structure in the best model or a competitive model (Table 19). The 95% CIs around the covariate coefficients for the Barred Owl effect did not overlap 0 for 4 areas (RAI, COA, HJA, and NWC) and barely overlapped 0 for 3 others (CLE, CAS, and HUP), and in all cases coefficients (or interactions) were negative, suggesting that the presence of Barred Owls was negatively associated with apparent survival in many study areas. Although the Barred Owl covariate was not modeled for GDR (because it was confounded with treatment effects), the best random effects survival model for GDR included higher apparent survival in treatment areas where Barred Owls were removed, although the 95% CI around the covariate coefficient for the treatment effect included 0 (After*Trt*T interaction; Table 19). Based on the best survival models that included either time trends or the negative effect of Barred Owl detections, we concluded that there was strong support for declining apparent survival in at least 8 of 11 study areas (CLE, RAI, HJA, TYE, CAS, NWC, HUP, and GDR).

TABLE 18.

Best model containing a time effect on apparent survival of nonjuvenile Northern Spotted Owls for each of 11 study areas in Washington, Oregon, and California, USA, 1985–2013. Model covariate coefficients (β̂), standard errors (SE), and 95% confidence limits (lower: LCL; upper: UCL) are included for the specified effect in each model.

TABLE 19.

Best model containing a Barred Owl (BO) effect on apparent survival of nonjuvenile Northern Spotted Owls for each of 11 study areas in Washington, Oregon, and California, USA, 1985–2013. Model covariate coefficients (β̂), standard errors (SE), and 95% confidence limits (lower: LCL; upper: UCL), are included for the specified Barred Owl effect in each model.

The effect of population-level reproductive rates (R) occurred in top models of apparent survival for 2 areas (RAI and OLY), and in competitive models for 3 areas (CLE, CAS, and HUP; Table 20). However, 95% CIs around the covariate coefficients for R overlapped 0 in all but 1 study area (RAI), suggesting that there was no support for an effect of R on survival in most study areas.

TABLE 20.

Best model containing a reproductive (R) effect on apparent survival of nonjuvenile Northern Spotted Owls for each of 11 study areas in Washington, Oregon, and California, USA, 1985–2013. Model covariate coefficients (β̂), standard errors (SE), and 95% confidence limits (lower: LCL; upper: UCL) are included for the reproductive effect in each model.

The amount of nesting and roosting habitat (HAB) did not occur in top or competitive models of apparent survival for any of the 10 areas for which it was investigated (Table 21). The amount of core habitat was the best habitat covariate overall for 1 area (KLA) and occurred in competitive models for 2 other areas (RAI and HUP); however, 95% CIs for all CORE covariate coefficients widely overlapped 0, suggesting little support for this effect. Where investigated, the amount of interface between nesting and roosting habitat and other cover types (EDGE) occurred in a top or competitive model for only 2 of 9 areas (CLE and COA). In both cases the 95% CIs around the EDGE covariate coefficients slightly overlapped 0, but, similarly to the fecundity analysis, the coefficients were positive, which was contrary to what we predicted. The annual amount of nesting and roosting habitat disturbance (HD) occurred in a top or competitive model for 2 of 10 areas (TYE and CAS), and was the best habitat covariate overall for the GDR study area. In both cases, higher estimates of the amount of reduction in nesting and roosting habitat were associated with lower survival, and in the TYE area this effect was included as part of a complex interaction with the proportion of territories in which Barred Owls were detected.

TABLE 21.

Best model containing a habitat effect on apparent survival of nonjuvenile Northern Spotted Owls for each of 11 study areas in Washington, Oregon, and California, USA, 1985−2013. Model covariate coefficients (β̂), standard errors (SE), and 95% confidence limits (lower: LCL; upper: UCL) are included for the specified habitat effect(s) in each model.

Climate covariates included in the analysis of survival (PDSI, SOI, PDO, ENP, ENT, LNP, LNT, WP, and WMT) occurred in top or competitive models for all study areas except GDR, but there was little consistency among areas as to which covariate was important (Table 22). In the CLE study area, the Palmer Drought Severity Index (PDSI) was positively related to apparent survival, and 95% CI for the covariate coefficient slightly overlapped 0. The Southern Oscillation Index (SOI) occurred in the top model for COA and as an additive effect with PDSI for CLE. In both cases the covariate coefficients were negative and the 95% CIs either did not overlap 0 (COA) or only slightly overlapped 0 (CLE). The Pacific Decadal Oscillation (PDO) covariate occurred in competitive models for 2 areas and in both areas the relationship with apparent survival was positive, with 95% CI around covariate coefficients that either did not overlap 0 (RAI), or slightly overlapped 0 (OLY). Precipitation during the early nesting period (ENP) occurred as a main effect in a competitive model for OLY, with the 95% CI slightly overlapping 0. In addition, ENP occurred in the best model for NWC as part of an interaction with temperature during the early nesting period (ENT). In both cases, increased precipitation during the early nesting period (ENP) was associated with decreased survival rates, but for NWC there was an ameliorating effect of temperature (ENT) associated with the interaction. Mean minimum temperature during the early nesting season (ENT) was negatively associated with survival for NWC, where it occurred in an interaction with ENP, and although the 95% CI for the ENT coefficient overlapped 0, the 95% CIs for the ENP and the ENT*ENP interaction term coefficients did not. In contrast, the effect of ENT was positive for GDR, with higher temperatures during the early nesting season associated with higher survival. Total precipitation during the late nesting season (LNP) was in a top or competitive model only for the KLA area, where it occurred in an interaction with late nesting season temperature (LNT). The main effects of LNT and LNP were positively related to survival, with 95% CIs around the covariate coefficients not overlapping 0, but the interaction was negative, suggesting that both high precipitation and high temperature in combination had a negative effect on survival. Mean minimum temperature during the late nesting season (LNT) was also in a competitive model for TYE; the effect was positive, with a 95% CI for the covariate coefficient that did not overlap 0. Mean minimum winter temperature (WMT) was in the top or competitive model for 3 study areas (HJA, CAS, and HUP), and the fact that 95% CI for the WMT coefficient did not overlap 0 suggested strong support for a positive relationship with apparent survival in the Oregon Cascades (HJA and CAS), but less support for an association with survival in the HUP area in northern California, where the 95% CI around the WMT covariate coefficient slightly overlapped 0. Total winter precipitation (WP) did not occur in a top or competitive model for any study area.

TABLE 22.

Best model containing a climate effect on apparent survival of nonjuvenile Northern Spotted Owls for each of 11 study areas in Washington, Oregon, and California, USA, 1985−2013. Model covariate coefficients (β̂), standard errors (SE), and 95% confidence limits (lower: LCL; upper: UCL) are included for the specified climate effect(s) in each model.

Meta-analysis of apparent survival

We based the meta-analysis of apparent survival on the same encounter histories that we used for birds in the individual study areas (5,090 nonjuvenile owls [≥2 yr of age]; Table 2). There was no indication that survival rates of owls in the 8 study areas associated with the monitoring program (NWFP in Table 23) differed from survival rates of owls in the 3 nonfederal study areas (β̂ = 0.007 ± 0.092 SE, 95% CI: −0.012 to 0.024; Table 23, 24). Therefore, we proceeded with the meta-analysis of all 11 study areas combined. The best random effects model included a common intercept across all study areas, positive effects of PDO, and negative effects of SOI (Table 23, 24). The 95% CI around PDO and SOI slightly overlapped 0 (Table 24), suggesting some support for higher survival when the PDO was in a warming phase (warmer, drier winters in the Pacific Northwest) and when the SOI was negative, with a strongly negative SOI indicative of El Niño events (warmer, drier winters in the Pacific Northwest; Figure 11). Compared with the random effects model with no covariates [RE(.)], the best model with PDO and SOI explained ~12% of the variation in apparent survival. The best model that included the Barred Owl covariate was ranked 4^th and also included PDO (Table 23). Compared with the model that included only the random effect of area [RE(Area)], the model that included Area, PDO, and BO explained ~13% of the variation in survival (Figure 12). The Barred Owl covariate was negatively associated with apparent survival across all study areas as predicted (β̂ = −0.057 ± 0.030 SE, 95% CI: −0.117 to −0.001; Table 24).

TABLE 23.

Model selection criteria for a priori random effects models (RE) in the meta-analysis of apparent survival of adult Northern Spotted Owls in 11 demographic study areas in Washington, Oregon, and California, USA, 1985−2013. The best fixed effects model, which included the interaction between area and year on survival and interactions among sex, study area, and year on detection rates [φ(Area*t) p(Sex*Area*t)], was used for all random effects modeling. Models were ranked according to Akaike's Information Criterion adjusted for small sample size (AIC_c). The model deviance (Deviance), number of parameters (K), difference in AIC_c between each model and the model with the lowest AIC_c (ΔAIC_c), and Akaike model weights (w_i) are given for all models.

TABLE 24.

Best random effects (RE) model that included the specified covariate in the meta-analysis of apparent survival of nonjuvenile Northern Spotted Owls in each of 11 study areas in Washington, Oregon, and California, USA, 1985–2013. Model covariate coefficients (β̂), standard errors (SE), and 95% confidence limits (lower: LCL; upper: UCL) are included for the effect listed (β̂ ^c label).

FIGURE 11.

Estimates of the effects of the Pacific Decadal Oscillation (PDO) and Southern Oscillation Index (SOI) on apparent survival of Northern Spotted Owls generated from the best random effects (RE) model from the meta-analysis [RE φ(PDO + SOI)] based on the fixed effects model [φ(Area*t) p(sex*Area*t)] with study area (Area) and general time variation (t) on apparent survival (φ), and study area, time, and sex effects on capture probability (p). Random effects model estimates are plotted with shrinkage estimates (S-tilde) for the (A) 3 study areas in Washington, (B) 5 study areas in Oregon, and (C) 3 study areas in California, USA. Only sites where Barred Owls were not removed (control areas) were included for the Green Diamond study area (GDR-C). See Table 2 for study area abbreviations.

FIGURE 12.

Estimates of the effects of Barred Owls (BO), the Pacific Decadal Oscillation (PDO), and study area (Area) on apparent survival (φ) of Northern Spotted Owls generated from the best meta-analysis random effects model (RE; ranked 4th) that included a Barred Owl effect [RE φ(Area + BO + PDO)] and was based on a fixed effect model [φ(Area*t) p(sex*Area*t)] with study area and general time variation (t) on apparent survival, and study area, time, and sex on capture probability (p). Random effects model estimates are plotted with original time-dependent survival estimates (MLE) and shrinkage estimates (S-tilde) for (A) 1 study area in Washington (CLE), (B, C) 2 study areas in Oregon (HJA, CAS), and (D) 1 study area in California (NWC). See Table 2 for study area abbreviations.

The best annual time trend model included the additive effects of study area and a cubic spline with 4 knots (Table 23). Negative trends in survival were evident from all trend models (Table 24). However, the trend models were all ranked relatively low, indicating that other covariates (e.g., PDO, SOI, and BO) better captured the overall temporal variation in survival across study areas (Table 23).