Partitioning global change: Assessing the relative importance of changes in climate and land cover for changes in avian distribution

Abstract Understanding the relative impact of climate change and land cover change on changes in avian distribution has implications for the future course of avian distributions and appropriate management strategies. Due to the dynamic nature of climate change, our goal was to investigate the processes that shape species distributions, rather than the current distributional patterns. To this end, we analyzed changes in the distribution of Eastern Wood Pewees (Contopus virens) and Red‐eyed Vireos (Vireo olivaceus) from 1997 to 2012 using Breeding Bird Survey data and dynamic correlated‐detection occupancy models. We estimated the local colonization and extinction rates of these species in relation to changes in climate (hours of extreme temperature) and changes in land cover (amount of nesting habitat). We fit six nested models to partition the deviance explained by spatial and temporal components of land cover and climate. We isolated the temporal components of environmental variables because this is the essence of global change. For both species, model fit was significantly improved when we modeled vital rates as a function of spatial variation in climate and land cover. Model fit improved only marginally when we added temporal variation in climate and land cover to the model. Temporal variation in climate explained more deviance than temporal variation in land cover, although both combined only explained 20% (Eastern Wood Pewee) and 6% (Red‐eyed Vireo) of temporal variation in vital rates. Our results showing a significant correlation between initial occupancy and environmental covariates are consistent with biological expectation and previous studies. The weak correlation between vital rates and temporal changes in covariates indicated that we have yet to identify the most relevant components of global change influencing the distributions of these species and, more importantly, that spatially significant covariates are not necessarily driving temporal shifts in avian distributions.


| INTRODUC TI ON
In recent decades, the distributions of many bird species have changed in Britain (Fuller et al., 1995;Thomas & Lennon, 1999), continental Europe (Böhning-Gaese & Bauer, 1996;Brommer, Lehikoinen, & Valkama, 2012), South Africa (Hockey, Sirami, Ridley, Midgley, & Babiker, 2011), and North America (Sauer, Link, Fallon, Pardieck, & Ziolkowski, 2013). These shifts are of both ecological and conservation interest, and they have stimulated a great deal of research. While various factors may be contributing to distributional shifts, including invasive species (Crowl, Crist, Parmenter, Belovsky, & Lugo, 2008), pollution (Trathan et al., 2015), exploitation (Laursen & Frikke, 2008) and other causes, climate change and land cover changes are likely to be two of the strongest drivers for many species (Barbet-Massin, Thuiller, & Jiguet, 2012;Lemoine, Bauer, Peintinger, & Böhning-Gaese, 2007;Travis, 2003). Understanding the relative impact of these drivers of distributional change has implications for ecological understanding, predictions of future changes, appropriate management strategies, and the allocation of conservation resources. However, species distribution models focused on the effects of climate are far more common than analyses that compare the contributions of climate and land cover (Sirami et al., 2017).
Ecologists have long debated the factors that affect species ranges (Grinnell, 1917;MacArthur, 1972) because of the relevance to ecology, biogeography, and community ecology. In recent decades, global change has given more urgency to the topic (Urban, 2015), and the expectation that both the global environment and species distributions will continue to change has stimulated much research into predicting the future course of change (Huntley et al., 2006). Questions about the relative impact of climate change and land cover change are important because of the implications for these predictions (Sirami et al., 2017). For example, climate change is expected to occur along a latitudinal gradient, while changes in land cover may exhibit less directionality if factors such as urbanization or conversion to or from agricultural use are significant drivers of land cover change (Lemoine et al., 2007;Travis, 2003). As a result, species distributions that respond primarily to climate are more likely to also exhibit directionality (Brommer et al., 2012;Hockey et al., 2011;Parmesan & Yohe, 2003). Whether species distributions respond to land cover or climate change has implications for appropriate responses, such as the design and location of reserves (Araújo, Cabeza, Thuiller, Hannah, & Williams, 2004), habitat restoration, and assisted migration (Hoegh-Guldberg et al., 2008). Accordingly, politically and economically significant decisions may be influenced by these projections. Therefore, there is a need for research into the relative effects of climate change and land cover change on changes in the distribution of species.
When there is interest in projecting changes in species distributions, a crucial insight is that environmental factors that correlate with the current species distribution may not correlate with future changes in distribution (Barbet-Massin et al., 2012). One reason is that projections developed from current distribution patterns assume that species are in equilibrium with their environment (Elith, Kearney, & Phillips, 2010), which may not be true (Zhu, Woodall, & Clark, 2012). A second reason is that dispersal limitations may prevent species from occupying their preferred habitat in the future (Devictor, Julliard, Couvet, & Jiguet, 2008). Therefore, investigating the processes that shape species distributions, rather than current species distribution patterns, is likely to generate greater ecological understanding and better projected distributions (Yackulic, Nichols, Reid, & Der, 2015). For example, the recent distribution of breeding Louisiana Waterthrush (Parkesia motacilla) in North America is best described by mean temperature, but changes in that distribution correlate with mean precipitation as well as mean temperature, indicating that the process cannot be fully described by the pattern observed at one point in time (Clement, Hines, Nichols, Pardieck, & Ziolkowski, 2016). We further note that data on environmental variables typically include both spatial and temporal components. For example, temperature may vary along a north-south gradient, as well as through time.
When the primary research interest is in temporal changes in distributions, it makes sense to isolate the temporal components of environmental variables (which are the essence of global change) when investigating species response to global change.
Here, we assess the relative importance of changes in climate and land cover to changes in distributions for two bird species, the Redeyed Vireo (Vireo olivaceus) and the Eastern Wood Pewee (Contopus virens). We selected just two species with similar biology, but differing population trends, so that we could develop species-specific hypotheses about distributional shifts. We partitioned the spatial and temporal components of our environmental predictors to improve our inferences about the processes affecting species distributions and about the likely consequences of future global change. We used dynamic occupancy modeling to explicitly estimate the vital rates that govern changes in distributions so as to avoid the equilibrium assumption implicit in static species distribution models. Finally, we used an analysis of deviance approach to assess the relative importance of climate and land cover to the changes in distributions for these two species.

| Hypotheses
Our goal was to assess the relative importance of climate and land cover to changes in the distributions of two bird species, the Redeyed Vireo and the Eastern Wood Pewee. We selected these birds because they are migratory, insectivorous, and inhabit similar plant communities (Hamel, 1992), but exhibit divergent trends in relative abundance reported by the BBS (Sauer et al., 2015). The Red-eyed Vireo breeds in the eastern and northern United States and much of southern Canada and winters in South America. It is a common and widespread species that uses a wide variety of forest habitats to nest and glean insects from foliage. Counts of Red-eyed Vireos have been increasing by 0.8% annually since 1966 and by 1.0% annually since 2003, with some localized declines in Florida, Texas, and the Pacific Northwest (Sauer et al., 2015). The Eastern Wood Pewee breeds in the eastern United States and parts of southeastern Canada. In contrast to Red-eyed Vireos, counts of the Eastern Wood Pewee have been declining by 1.5% since 1966 and by 1.2% since 2003 (Sauer et al., 2015). Despite the generally negative trend, there have been localized increases in counts in the Midwest.
To estimate the relative importance of climate and land cover to changes in bird distributions, we generated testable hypotheses expressing relationships between birds and environmental covariates.
A common approach is to generate an extensive set of hypotheses from a suite of environmental metrics. We can then fit the specified models and observe which metrics yield significant results. However, testing many statistical hypotheses has a propensity to identify spurious relationships and such analyses should be considered exploratory (Anderson, Burnham, Gould, & Cherry, 2001;Ioannidis, 2005).
In this study, we developed two a priori hypotheses from ecological principles: 1. Changes in the total annual periods of extreme temperature drive changes in bird distributions. We focus on extreme temperature because we expect these are the periods of greatest stress due to cold or heat stress, increased metabolic costs, reduced resource availability, or other mechanisms (Dawson & Whittow, 2000;Root, 1988). We defined "extreme" relative to the thermoneutral zone (TNZ), which is the range of ambient temperatures under which endotherms can maintain their body temperature without deviating from their basal metabolic rate (Calder & King, 1974). The TNZ varies among species, but has been estimated as 18-38°C for a "generic" bird (Calder & King, 1974). Several small, temperate-zone passerines (similar to our study species), including the Northern Cardinal, Verdin, House Sparrow, and Red-breasted Nuthatch, have been documented to have TNZs similar to that of this "generic" bird (Dawson, 1958;Goldstein, 1974;Hinds & Calder, 1973;Khaliq, Hof, Prinzinger, Böhning-Gaese, & Pfenninger, 2014). Therefore, we expected temperatures below 18°C and above 38°C to be associated with changes in bird distributions.
2. Changes in the amount of appropriate habitat locally available during the breeding season drive changes in bird distributions because inappropriate habitat will not provide sufficient food and shelter for breeding birds (Friggens & Finch, 2015;Iglecia, Collazo, & McKerrow, 2012). For the Eastern Wood Pewee, we defined appropriate habitats as deciduous forest, evergreen forest, and mixed forest (Hamel, 1992). For the Red-eyed Vireo, we defined appropriate habitats as the same forest types and shrub-scrub (Hamel, 1992).
Of course, these hypotheses are not exhaustive, but we selected them as logical, general mechanisms that might underlie avian responses to climate and land cover change. We expected the paucity and specificity of our hypotheses would reduce our chances of obtaining statistically significant results by chance. For example, a null hypothesis that bird distributions will not expand at a specific temperature is harder to reject than the null hypothesis that bird distributions will not expand at any temperature. We view a more specific hypothesis as a more powerful discriminatory tool (Platt, 1964;Popper, 1959), more in keeping with the hypothetico-deductive method (Chamberlin, 1890;Romesburg, 1981) and with a Type 1 error rate closer to the nominal level (Anderson et al., 2001;Ioannidis, 2005).
We stated our hypothesized effects in terms of changes in bird distribution (i.e., birds will increase where habitat is available) rather than static patterns of bird distribution (i.e., birds will be located where habitat is available). We prefer this formulation because it focuses on the ecological processes that underlie species distributions and their dynamics. In metapopulation theory, the vital rates describing these processes are often the local colonization and extinction rates. Models that account for ecological vital rates (i.e., dynamic models) are useful because they are relatively mechanistic, thereby strengthening the generality of findings and contributing more to ecological understanding. In contrast, in order to be useful for forecasting and projection, relationships between environmental covariates and occurrence require an assumption that the species in question is in equilibrium with the environmental factors under study (Elith et al., 2010). The more phenomenological orientation of models that focus on occurrence can be limiting when projecting model results into novel situations, such as future climate and land cover conditions (Cuddington et al., 2013;Gustafson, 2013). In contrast, dynamic models avoid the equilibrium assumption, enabling more robust predictions of occurrence under future conditions (Clement et al., 2016;Yackulic et al., 2015).
Both climate and land cover vary spatially and temporally, and each process could impact bird distributions. It is possible for the level of correlation with bird vital rates to differ across those two dimensions. For example, consider an environmental variable that is completely static through a given time period and a species distribution that shifts during the same time period. In this case, it is not possible for change in the static variable to have caused the distributional shift, because no such change occurred. However, it is still possible for colonization and extinction to correlate with the environmental variable because of a spatial correlation between them.
Here, our goal was to incorporate temporal variation in covariate values into our analysis because much research into global change is motivated by interest in temporal changes in species distributions. Therefore, our hypotheses and models are structured to isolate and evaluate the ability of environmental covariates to account for spatial and temporal variation in local colonization and extinction probabilities.

| Data
We obtained data on the detection and nondetection of breeding Red-eyed Vireos and Eastern Wood Pewees from the North American Breeding Bird Survey (Pardieck, Ziolkowski, & Hudson, 2015;www.pwrc.usgs.gov/BBS/RawData/). The Breeding Bird Survey (BBS) is a joint project of the U.S. Geological Survey and the Canadian Wildlife Service, designed to monitor the status and trends of birds breeding in North America (Sauer et al., 2013). Since its inception in 1966, the BBS has expanded to include over 5,000 survey routes across the United States and southern Canada, approximately 3,000 of which are surveyed each year. Each route is surveyed once per year by a skilled volunteer. The date of the survey is timed to occur during peak territorial behavior, typically late May to early July, depending on latitude. Surveyors follow a prescribed route along secondary roads, stopping at approximately 800-m intervals, and performing 3-min roadside point counts, generating 50 counts per 39.4 km route. We selected BBS data because their great geographic and temporal extent is suitable for investigating distributional changes.
We selected a subset of the available BBS data for use in our analysis. Our modeling approach, described below, makes use of stop-specific survey results, but these detailed results are currently available in digital formats only since 1997. Therefore, we used data from 1997 to 2012. Because our study species are found in forested regions of eastern North America, we delineated the study area by aggregating three Omernik Level I classes representing eastern forest ecoregions: Northern Forests, Eastern Temperate Forests, and Tropical Wet Forests (Commission for Environmental Cooperation, 1997;Omernik, 1995Omernik, , 2004Omernik & Griffith, 2014). This study area includes the United States east of the Great Plains (2.9 million km 2 ). This region encompassed the entire breeding range of the Eastern Wood Pewee, except some areas in Canada not covered by our climate and habitat covariates. We selected BBS routes entirely within this study area for analysis. We excluded routes that were not classified as "active" by the BBS and routes that differed from the recommended length (39.4 km) by more than 25%. We only analyzed surveys conducted under acceptable conditions (i.e., acceptable weather, date, time of day, stops; reported by BBS as "runtype = 1").
On the rare occasions that a route had multiple acceptable surveys per year, we used data from the first acceptable survey. This yielded 1,371 routes for our analysis. Finally, we converted counts of birds at each stop along each BBS route to detection/nondetection data suitable for presence-absence modeling. This conversion sacrificed data richness, but allowed us to account for imperfect detection of species. Using the full count data typically requires additional assumptions (Illán et al., 2014) or incorporation of additional data sources (Hooten, Wikle, Dorazio, & Royle, 2007).
Climate data were obtained from the PRISM project (Daly et al., 2008). The PRISM output uses locally weighted regression models to interpolate long-term weather station observations to a 4 km resolution grid. Daily interpolated climate data are available from 1981 to the present. From these data, we calculated two temperature-based indices based on our first hypothesis, described above. The first index was calculated as the total annual hours above 38°C, while the second index was calculated as the total annual hours below 18°C.
Because the temperature-based indices were based on hours of exposure, it was necessary to interpolate the daily PRISM data to estimate hourly values. We considered a simple linear interpolation to impute values between daily maximum and minimum temperatures.
However, for extreme temperatures, this method would likely result in estimated exposure durations that are biased low. Instead, we adopted the method described by Cesaraccio, Spano, Duce, and Snyder (2001; see Figure 3) and Eccel (2010) using the R software package "Interpol.T" (Eccel & Cordano, 2013). In this method, candidate functions are fit to observed hourly temperature data from nearby high-quality ASOS (Automated Station Observing System) weather stations in the study region using a combination of sinusoidal and quadratic functions. For each 4 km PRISM grid cell, the resulting calibration function from the nearest ASOS station is then applied to each day's maximum and minimum temperature to predict the intervening hourly values.
Using the interpolated hourly temperature estimates, we calculated the total annual hours within the zones of heat stress or cold stress for each year from 1981 to 2012. Annual values were calculated over the period June 1-May 31 to represent the breeding cycle. Once the annual residence times were obtained, we applied a 15-year moving average window to the data, resulting in 16 moving average periods starting with the 1981/82-1995/96 period and ending with the 1996/97-2011/12 period. Applying a moving average acts as a high-pass filter that attenuates interannual variability while amplifying any long-term climate change signals. Other period lengths could have been chosen; however, we hypothesized that 15 years represented a reasonable tradeoff between (a) reducing the moving average length so as not to exceed the length of the daily PRISM data series, which begins in 1981, (b) increasing the moving average length so that climatic changes (as opposed to background climatic variation) can be identified, and (c) using a period that is biologically meaningful in that we can reasonably expect any long-term climatic changes in the moving averages would be outside the historical range of variability and thus should affect colonization and extinction rates of sensitive species.
We obtained land cover data from the National Land Cover Database (NLCD) for 2001 (Homer et al., 2007), 2006 (Fry et al., 2011), and 2011 (Homer et al., 2015). We reclassified the given land cover types into habitat or nonhabitat for each bird species, based on our second hypothesis, described above. We calculated the percent of land area that was a suitable habitat within 400 m of each route and used this as a covariate in our analysis. We used 400 m because this is half the distance between stops on BBS routes.

| Analysis
Our goal was to assess the relative importance of climate and land cover to changes in the distribution of Red-eyed Vireos and Eastern Wood Pewees. While our bird data consisted of detection/nondetection data from BBS surveys, our parameters of interest were the probabilities of local colonization and local extinction, as these processes drive changes in distribution. Our predictor variables were our measures of climate and land cover suitability for each species, described above. We used dynamic correlated-detection occupancy models with detection heterogeneity to estimate the relationship between our predictors and bird distribution dynamics (Clement et al., 2016;Hines, Nichols, & Collazo, 2014;Hines et al., 2010). Occupancy models are a class of models that estimate the probability of presence while accounting for imperfect detection of species by analyzing replicated (usually temporally) surveys (MacKenzie et al., 2002(MacKenzie et al., , 2018. Dynamic occupancy models model interannual changes in occupancy by positing a Markov process in which the probability of occurrence at time t is a function of the probability of occurrence at time t−1 (MacKenzie, Nichols, Hines, Knutson, & Franklin, 2003). In contrast to static models, this Markov process recognizes that the distribution of a species is a function of previous environmental conditions and dispersal constraints, as well as current environmental conditions (Dormann et al., 2012;Kéry, Guillera-Arroita, & Lahoz-Monfort, 2013). Dynamic correlated-detection occupancy models are a model extension in which spatially correlated replicated surveys are used to estimate detection probabilities (Hines et al., 2014(Hines et al., , 2010. We selected an occupancy modeling approach because failure to account for imperfect detection can bias estimates of colonization and extinction rates (Ruiz-Gutiérrez & Zipkin, 2011). We selected a dynamic approach because dynamic models offer more ecological realism, more accurate projections, and greater generalizability than static models (Clement et al., 2016;Yackulic et al., 2015). We used correlated-detection models because the BBS generates spatially replicated surveys, and failure to account for spatial correlation of replicates can bias occupancy estimates (Hines et al., 2014(Hines et al., , 2010. Finally, we used a finite mixture model to account for detection heterogeneity because of the variation in habitat and focal species abundance across the study area (Clement et al., 2016). A finite mixture model approximates the heterogeneity of detection probabilities by positing that the population consists of a mixture of routes which have either a relatively high detection probability or a relatively low detection probability.
As detailed above, BBS data consist of numerous routes, each composed of bird observations at 50 stops. A species may be absent from individual stops when it is present on a route, but not vice versa. Detection of a species is taken as proof of presence, but a nondetection may result from a true absence or a false absence. For clarity, we use the term "occupied" to indicate a species is present on a route and the term "available" to indicate a species is locally present at a specific stop at the time of the survey (Nichols, Thomas, & Conn, 2009 π = Pr (probability that a route is a low-detection route); ε t = Pr (route is not occupied in season t + 1 | occupied in season t); and γ t = Pr (route is occupied in season t + 1 | not occupied in season t). Hines et al. (2014Hines et al. ( , 2010 developed a model likelihood that allows these parameters to be modeled as functions of route-and year-specific covariates, while detection parameters can also be influenced by stop-specific covariates. Parameters can then be estimated using maximum-likelihood estimation in program PRESENCE (Hines, 2006).
We developed six specific models of γ and ε for each bird species.
We partitioned spatial and temporal variation in covariates by partitioning covariates into the average value over time, and the annual deviation from that average. In this way, the averaged covariate, x i , included spatial variation only, while the deviation component, Δx iy , included temporal variation only. Using these covariates, our models included (a) a null model in which γ and ε are constant through time and space and (b) a global (most general) model in which γ and ε vary through space according to mean climate and land cover suitability, and through time using a dummy covariate for each year. We also considered intermediate models to assess the explanatory power of our covariates: (c) γ and ε vary spatially with mean climate and land cover suitability, but not with time, (d) γ and ε vary spatially with mean climate and land cover suitability, and temporally with annual changes in climate only, (e) γ and ε vary spatially with mean climate and land cover suitability, and temporally with annual changes in land cover only, and (f) γ and ε vary spatially with mean climate and land cover suitability, and temporally with annual changes in climate and land cover.
Prior to comparing these models, we worked to achieve adequate fit for the other model parameters. Because of the number of model parameters, we considered them sequentially. Initially, we fit the global model for γ and ε, and we modeled θ, θ′, and ψ as functions of climate and land cover covariates. In this initial model, we used a finite mixture model on p to account for heterogeneity, presumably caused by differences among routes in habitat, observers, bird abundance, and other features that were not explicitly modeled (Clement et al., 2016). We also modeled p as a function of climate and land cover covariates and year. Because bird activity often varies with time of day, we also modeled detection as a quadratic function of stop number, which we consider to be a proxy for time of day. From this highly parameterized model, we fit a reduced model by eliminating the covariate with the smallest estimate-to-standard error ratio. We continued to eliminate covariates in this way until AIC increased and then selected the model with the lowest AIC. We used this parameter reduction protocol with ψ, θ, θ′, and p in turn, maintaining the model structure on γ and ε. After identifying parsimonious models for ψ, θ, θ′, and p, we fit the set of six models for γ and ε described above. The sequential approach may not yield identical results as a single analysis, but it is a pragmatic approach for dealing with complex models (MacKenzie et al., 2018). By beginning with highly parameterized models and progressing toward reduced models, we reduced the risk of confounding the occupancy and detection processes (MacKenzie et al., 2018). We also checked for potential multi-collinearity issues by calculating Pearson's correlation coefficient for the different covariates used in the models. We treated an R 2 < 0.5 as indicative of a lack of multi-collinearity.
We also evaluated the goodness of fit of the global model for each species using the naïve colonization and extinction rates to calculate the test statistic of a Hosmer-Lemeshow test (Hosmer & Lemeshow, 2000). In this context, the naïve colonization rate is the number of sites that changed from nondetection sites in year t to detection sites in year t + 1, while the naïve extinction rate is the number of sites that changed from detection sites in year t to nondetection sites in year t + 1. We divided the routes into deciles based on the expected colonization and extinction rates, and compared the predicted rates to the observed rates with chi-square tests (α = 0.05). If necessary, we combined deciles to ensure that the predicted number of detections in a decile was >4. Although the Hosmer-Lemeshow test statistic is often calculated from predicted presence, vital rates seemed to be more appropriate measures in this case because dynamic models estimate changes in presence, rather than presence (Clement et al., 2016). We note that this use of naïve rates focuses on products of model parameters (i.e., γ, ε, p, ψ, ϑ), and therefore, we do not test the fit of the fully parameterized model. We focused on naïve rates because the true vital rates (i.e., γ and ε) cannot be directly observed when detection is imperfect, while the naïve rates can be observed.
We used an analysis of deviance approach (ANODEV; e.g., Skalski, 1996) to evaluate the relative importance of climate and land cover to changes in bird distributions. In this procedure, we compared the deviance explained by our parameter of interest to the deviance explained by the most fully parameterized model in the model set: Dev(M null ) is the deviance of Model 3, which accounts for spatial, but not temporal, variation in colonization and extinction probabilities. Dev(M full ) is the deviance of Model 2, the global model, with both spatial and temporal variation in γ and ε. We then evaluated the deviance explained by climate and/or land cover by using models 4, 5, and 6 for Dev(M cov ). Higher values of R 2 Dev indicate greater explanatory power, with 0 indicating no power, and 1 indicating power equal to that of the fully parameterized model, although R 2 Dev could potentially exceed 1 because our covariate models are not strictly nested within Model 2. Therefore, we calculated and report R 2 Dev for both climate and land cover. This is not the only possible approach to assessing relative importance of explanatory variables, but it is sound and has been used and recommended by others (Grosbois et al., 2008).
The general ANODEV approach to partitioning variation led us to focus on hypothesis testing as a means of assessing the need for extra model parameters (sources of variation). The small set of specific a priori hypotheses also sets this work apart from exploratory studies investigating the adequacy of many different models.
For these reasons, we used likelihood ratio tests as a means of comparing different models, although we note that use of model selection criteria (e.g., AICc) yielded similar inferences.

| RE SULTS
Our covariates indicated a small loss of habitat and increased temperatures during the study period, despite the relatively short time scale. For Eastern Wood Pewee, habitat change was >1% on 38% of routes, with the mean share of suitable habitat near routes declining from 41.6% to 40.9%. For Red-eyed Vireo, habitat change was >1% on 23% of routes, with mean suitable habitat declining from 48.0% to 47.6%. Route temperatures were rarely above the heat stress threshold, but the average time above this threshold per route increased from 0.5 hr per year from 1982 to 1997 to 1.7 hr per year from 1997 to 2012. We also observed a decline in periods of cold stress, from 5,450 to 5,291 hr per route per year.
The Pearson's correlation coefficient did not indicate any multicollinearity problems among our covariates (R 2 <0.2 in all cases).

| Eastern Wood Pewee
For Eastern Wood Pewee, the best-supported model included both climate and habitat measures as covariates of ψ, θ, and θ′. We excluded climate and habitat covariates from the detection probability model because some of those models did not converge. Therefore, the detection model for each species included year and stop number as covariates. As a result, the global model included 85 parameters (Table 1). The Hosmer-Lemeshow test indicated that the global model fit the data adequately, with no evidence of discrepancies between estimated and observed naïve colonization rates (χ 2 = 4.41, df = 8, p = 0.82) and naïve extinction rates (χ 2 = 11.17, df = 8, p = 0.19).
Eastern Wood Pewee were widespread in the study area, with an average probability of occupancy on BBS routes in 1997 of 0.90 in the global model (Figure 1). Occupancy declined slightly to 0.88 by 2012, mirroring the decline in counts reported by the BBS. Initial occupancy was positively related to available habitat and hours of cold stress (Table 1). Probability of detection at the first individual stop during 1997 was low, at 0.13, while probability of detection at the route level was much higher at 0.95 because birds were typically available at many stops within a route. Model 3, which included spatial variation in climate and habitat in the colonization and extinction models, significantly improved model fit, relative to Model 1 (likelihood ratio χ 2 = 44.19, df = 6, p < 0.001). Under this model, estimated colonization rates in 1998 ranged among routes from 0.08 to 0.28 (with a mean of 0.18), while extinction rates varied from 0.01 to 0.04 (with a mean of 0.02; Figure 2). Model 2 differed from Model 3 in that it allowed colonization and extinction rates to vary each year. Model 2 improved model fit relative   (Figure 3).
By 2012, this pattern had reversed, yielding a relatively weak occupancy trend in the southwest (Figure 3). None of these models represented a significant improvement over Model 3 (p > 0.26). In some cases, the sign of the estimated effect of temporal variation in climate and habitat was opposite of the expected effect, but in no instances were these estimates significant (Table 3).

| Red-eyed Vireo
The best-supported model for Red-eyed Vireo included the same covariates as for Eastern Wood Pewee, yielding a global model with 85 parameters (Table 1) (Table 3).
However, we suggest that estimated relationships between spatial variation in environmental covariates and colonization and extinction rates are of greater ecological interest than relatively phenomenological species distribution models or climate envelope models (Clement et al., 2016). The estimated vital rates of colonization and extinction bring us closer to understanding the dynamic process underlying the distribution of these species. Given constant vital rates, we can also project the equilibrium level of occupancy at sites ( i i + i , at site i). We can compare these projected equilibria to current occupancy to project trends for the species under current environmental conditions. We can also use projected changes in key environmental variables to develop predicted range changes in the future. We argue that these projections are more robust than those from static climate envelope models, which assume that species are currently in equilibrium (Yackulic et al., 2015). In contrast, static models can overestimate range changes when the equilibrium does not hold (Morin & Thuiller, 2009).
Given the importance of global change, we were motivated to investigate factors affecting temporal, as well as spatial, variation in vital rates. In this case, spatial variation must be incorporated into such analyses in order to accommodate the different dynamics expected in different portions of a species range. We view the relationship between temporal variation in covariate values and vital rates as the most relevant to understanding the effect of global change on the recent and future distributions of our study species.
Furthermore, a correlation between vital rates and spatially varying covariates may not indicate a similar correlation with temporal changes in those same covariates. For example, a species with strong dispersal constraints might have a strong spatial relationship with a covariate, but a weak temporal relationship due to its immobility (Devictor et al., 2008). Similarly, an estimated spatial relationship may be purely phenomenological due to the nature of spatial variables. When spatial variables follow gradients, it is relatively easy to obtain statistically significant correlations, absent any causal relationship (Grosbois et al., 2008). For example, if colonization rates and temperatures both follow north-south gradients, they may be spatially correlated even if colonization is governed by some other, unidentified factor. In this case, we would expect the weak correlation between temporal changes in temperature and colonization F I G U R E 2 Probability of (a) colonization and (b) extinction for Eastern Wood Pewee in the eastern United States, following 1997, under Model 3. Colonization and extinction vary by the amount of suitable habitat, hours of heat stress, hours of cold stress, all averaged over the 1997-2012 study period rates to be more indicative of the true relationship than the strong correlation with spatial variation in temperature.
BBS surveys indicate different abundance trajectories for our study species, with counts increasing for Red-eyed Vireos and decreasing for Eastern Wood Pewees (Sauer et al., 2015). Because abundance and occupancy tend to be linked (Gaston et al., 2000), we expected some divergence in occupancy dynamics. We found that estimated occupancy was high and stable for both species, although there was greater turnover for Red-eyed Vireos. For example, under Model 3 (spatial variation only in covariates), both average colonization rates (0.49) and average extinction rates (0.04) were at least twice as high for Red-eyed Vireos, compared to Eastern Wood Pewee (0.18 and 0.02, respectively). Red-eyed Vireos exhibited greater spatial variation in colonization and extinction rates (Figure 4), although our selected covariates explained little temporal variation for this species. Eastern Wood Pewee colonization and extinction rates exhibited more of a latitudinal gradient, reflecting a spatial correlation with extreme temperature (Figure 2). Our results also hinted at a greater correlation between temporal changes in temperature and occupancy vital rates for Eastern Wood Pewee, which could play a role in the decline of this species, although the results were not significant (Table 3). Alternatively, it has been suggested that Eastern Wood Pewee may be declining due to a loss of forested wintering habitat (Robbins, Sauer, Greenberg, & Droege, 1989) or due to deermediated changes in forest structure on breeding grounds (deCalesta, 1994), illustrating the complexity of ascertaining cause and effect using observational studies in these natural systems.
Although relatively rare (Sirami et al., 2017), other studies have investigated both climate and land cover effects on bird distributions. In agreement with our finding that temporal variation in climate explained much more of the model deviance than temporal variation in habitat, climate-only species distribution models had better cross-validation support than habitat-only models for 409 European bird species (Barbet-Massin et al., 2012). Similarly, models estimating the number of newly occupied areas by Hooded Warblers (Wilsonia citrina) in Ontario using climate data had better information criterion support than models using forest data (Melles, Fortin, Lindsay, & Badzinski, 2011). In contrast, land-use models outperformed climate models for 18 farmland bird species in the United Kingdom (Eglington & Pearce-Higgins, 2012). Similarly, vegetationbased species distribution models had superior performance than climate-based models for 79 Spanish bird species (Seoane et al., 2004) and for three bird species in New Mexico (Friggens & Finch, 2015). Two studies focused on colonization and extinction rates found mixed results. For Garden Warblers (Sylvia borin) in Britain, vegetation explained more variation in colonization, while temperatures explained more variation in local extinction (Mustin, Amar, & Redpath, 2014). For 122 bird species in Ontario, the best-supported covariates (land cover or climate change) varied among species (Yalcin & Leroux, 2018  be particularly affected by changes in agricultural intensity because of their extensive use of farms (Eglington & Pearce-Higgins, 2012).
Alternatively, it has been argued that spatial scale may play a role in the relative importance of global change factors, with climate more important at large scales (Pearson & Dawson, 2003). We note that two studies that found a larger role for climate both relied on lowerresolution bird surveys (0.5° atlas in Barbet-Massin et al., 2012; 10 km 2 atlas in Melles et al., 2011). In the BBS data that we analyzed, surveys were conducted every 800 m, but we calculated occupancy across each 39.4 km transect. In addition, the cited studies used different data sources and diverse analysis methods, complicating interpretation of the differing results.
While we estimated that climate explained more model deviance than habitat, these estimated relationships were not F I G U R E 3 Annual deviation of colonization and extinction rates from average colonization and extinction rates(̂(t) −̄,̂(t) −̄), for Eastern Wood Pewee in the eastern United States, following 1997 and 2011, under Model 6. Deviations in colonization and extinction rates vary with deviations in the amount of suitable habitat, hours of heat stress, hours of cold stress, relative to average levels over the 1997-2012 study period. (a) Colonization deviation following 1997, (b) extinction deviation following 1997, (c) colonization deviation following 2011, and (d) extinction deviation following 2011 TA B L E 3 Estimated parameters for correlated-detection dynamic occupancy Model 6 relating habitat and climate covariates to occupancy dynamics of breeding Eastern Wood Pewee and Red-eyed Vireo, 1997 Parameter estimates indicate change in log odds of parameters in response to z-transformed covariates. Parameters ψ,ϑ,ϑ′,γ,ε,p1,p2,and π  F I G U R E 4 Probability of (a) colonization and (b) extinction for Red-eyed Vireo in the eastern United States, following 1997, under Model 3. Colonization and extinction vary by the amount of suitable habitat, hours of heat stress, hours of cold stress, all averaged over the 1997-2012 study period significant and did not explain a large fraction of variation. There are several potential explanations for why we did not find a substantial response to changes in habitat and climate. First, our data were lacking in dynamic events to inform our models. The covariates changed relatively little over the time period, making it difficult to detect a response. The lack of change was partly because of the relatively short (16-year) time period studied. Also, we used a moving average for climate data, which tends to reduce differences over short time periods. In addition, the study region has experienced less warming over the past century than other parts of the planet (Walsh et al., 2014). The bird distributions also lacked dynamism. Because extinction rates were low, there were few local extinction events to inform the model. Although colonization rates were high, the fact that initial occupancy was high and extinction rates were low meant there were few opportunities for colonization events.
F I G U R E 5 Annual deviation of colonization and extinction rates from average colonization and extinction rates, for Red-eyed Vireo in the eastern United States, following 1997 and 2011, under Model 6. Deviations in colonization and extinction rates vary with deviations in the amount of suitable habitat, hours of heat stress, hours of cold stress, relative to average levels over the 1997-2012 study period. (a) Colonization deviation following 1997, (b) extinction deviation following 1997, (c) colonization deviation following 2011, and (d) extinction deviation following 2011 Second, measurement error may have obscured responses to global change. The PRISM climate data are estimated by interpolation, rather than direct measurement across the study area (Daly et al., 2008). As such, there may be artifacts or anomalies that impeded our estimation and increased uncertainty. Similarly, NLCD data are remotely sensed and may include some measurement error, which can increase uncertainty in our estimates (Homer et al., 2015).  , precipitation;Illán et al., 2014). It is also possible that the birds respond to finer habitat features than the land cover types we used, or at a different spatial scale (i.e., not 400 m; Pearson & Dawson, 2003). Alternatively, there may be important regional differences or interactions among covariates that we did not capture. Or, it is possible that changes in climate and habitat in their breeding areas simply have little effect on vital rates for these species relative to other factors, such as competition with other bird species, or global change occurring on their wintering grounds in South America (Dormann, 2007). Of course, these misspecification issues are inherent to all observational studies, and not unique to this study.
It would be possible for us to explore some of the misspecification issues just mentioned. For example, we could rerun the same models using new covariates representing a variety of temperature thresholds and habitat buffer sizes to search for a stronger correlation. Experience suggests that persistence would eventually be rewarded with a p-value <0.05 or a low model selection value. However, one reason to avoid such exploration is that it tends to lead to Type I errors, overfit models, and biased parameter estimates (Fieberg & Johnson, 2015). Furthermore, we find the contrast between significant correlations for spatial covariates and nonsignificant correlations for temporal covariates illuminating. In typical climate envelope studies, we would estimate the relationship between covariates and the current distribution of a species, much like our estimate of initial occupancy (Dormann, 2007 Center. The contents of this document are the responsibility of the authors and do not necessarily represent the views of the USGS. Any use of trade, product, or firm names is for descriptive purposes only and does not imply endorsement by the US Government.

D ED I C ATI O N
We dedicate this paper to the late Chandler S. Robbins, "Father" of the North American Breeding Bird Survey (BBS). Data from the BBS form the basis for this specific paper and for hundreds of others focused on topics ranging from local trends to continental macroecology. Chan was a unique combination of bird-lover and scientist whose vision and inquisitive enthusiasm were an inspiration to all of us fortunate enough to have known him.

CO N FLI C T O F I NTE R E S T
None declared.

AUTH O R S' CO NTR I B UTI O N
JAC and JDN conceived the study. AJT and SW analyzed habitat and climate data. MJC, JEH, and JDN performed the occupancy analysis. MJC led the writing of the manuscript. All authors gave final approval for publication.