Interactions Among Spatial Scales Constrain Species Distributions in Fragmented Urban Landscapes

Understanding species’ responses to habitat loss is a major challenge for ecologists and conservation biologists, who need quantitative, yet practical, frameworks to design landscapes better able to sustain native species. I here develop one such framework by synthesizing two ecological paradigms: scale-dependence and constraint-like interactions in biological phenomena. I develop a model and apply it to birds around Tucson, USA, investigating the manner in which spatial scales interact to constrain species distributions in fragmented urban landscapes. Species’ responses vary in interesting ways. Surprisingly, most show situations in which habitat at one spatial scale constrains the influence of habitat at another scale. I discuss the implications of this work for conservation in human-dominated landscapes, and the need to recognize constraint-like interactions among processes and spatial scales in ecology.


INTRODUCTION
The distribution of species across ecological landscapes, long a central focus of ecology, now serves as the basis for much research in conservation biology (Scott et al. 2002).Human-caused habitat loss is a leading cause of species endangerment worldwide (Wilcove et al. 1986).Consequently, the future of efforts to conserve biodiversity will be shaped by our ability to understand the varying responses of species to fragmented landscapes, and to translate that understanding into practical conservation solutions.
Researchers have used two general approaches to investigate the distribution of species in fragmented landscapes.A first approach emphasizes the ecological processes underlying the relationship between species' distributions and fragmentation.This approach potentially allows direct investigation of how processes such as dispersal (Hanski et al. 1996), predation (Tomialojc 1982, Haskell et al. 2001), competition (Shochat et al. 2004), survival, andreproduction (Boal andMannan 1999) mediate the response of species to fragmentation.But the application of process-based models is often made impractical by difficulty in generalizing them over heterogeneous landscapes, and by the great cost of data necessary for parameter estimation, e.g., demographic or dispersal data, particularly when dealing with new locations or species.
A second approach uses statistical, "empirical," or "phenomenological" models.This approach involves regression or other statistical modeling techniques applied to species data and one or more environmental variables.These methods generally can be applied to less costly data than process-based models.However, they emphasize patterns that emerge from data rather than a priori biological knowledge, and may suffer from inadequate representation of biological processes (e.g., O'Connor 2002, Van Horne 2002).
To date, much work has been done to document negative impacts of habitat loss and fragmentation in urban areas, and to understand these impacts in terms of underlying biological processes.But so far, recommendations for management remain more heuristic than pragmatic.Abstract and unidirectional recommendations such as "larger habitat areas are better," "animal species X needs more of plant species Y," or "the landscape scale is also important" are of limited practical use.Given many competing uses for urban land, quantitative predictions are essential for efficiently allocating http://www.ecologyandsociety.org/vol11/iss2/art6/conservation resources.For example, quantitative predictions may indicate not only what conservation action is needed at a site, but also how much of it is needed, preventing allocation of scarce resources, e.g., money, beyond that necessary.Within a larger planning context, quantitative biological predictions might additionally indicate sites for which too much effort would be required for restoration, allowing the more effective allocation of resources to other sites.
Managers and policymakers need quantitative, yet practical, frameworks that can guide them in mitigating and reversing the loss of wildlife as urban development continues.In this paper, I develop one such framework using as a study system the birds of Tucson, Arizona, USA.I sought to balance the best features of existing approaches by adapting a statistical model that, though it does not explicitly model biological processes, is nonetheless tailored to account for them implicitly.The resulting model incorporates spatial scale, incomplete occupancy of suitable habitats, and the fact that biological processes may constrain one another as well as augment one another.I use the model to address the questions of whether, and in what manner, species' spatial distributions are shaped by constraint-like interactions among spatial scales.Of importance to practical conservation efforts, the model generates quantitative predictions of bird occurrence from habitat variables that are straightforward both to quantify and to modify within a larger planning context.

Study area
The Tucson metropolitan area encompasses ~1300 km 2 in southern Arizona, USA at around 780 m elevation.Original habitats immediately surrounding Tucson are predominately upland Sonoran desertscrub (Brown and Lowe 1980), comprising various shrubs, trees exhibiting shrub-like growth, and cacti (Fig. 1a).Several natural reserves to the west, north, and east of Tucson retain natural vegetation cover.Developed areas vary, ranging from mostly natural vegetation (Fig. 1b) to entirely impervious cover and non-native vegetation (Fig. 1c).Tucson's population grew 20.1% during the 1990s (United States Census Bureau 2000), while urbanized land area grew even faster.Stressors to the surrounding desert ecosystem include urbanization, aquifer depletion, and surface water diversion (Nabhan and Holdsworth 1998).

Bird data
To collect bird data for this study, I organized the Tucson Bird Count (TBC), a monitoring program in which skilled observers count birds at sites in and around the Tucson metropolitan area.All observers met minimum qualifications and were professional biologists, professional field guides, or birdwatchers with years of experience.I placed one TBC count site randomly within each 1 km 2 cell of a regular grid.During the breeding season each spring, observers counted birds of all species using one 5minute point count at each site.To reduce spatial autocorrelation in observer identity among sites, observers surveyed "routes" of sites arranged as linearly as possible rather than surveying clusters of nearby sites.See Turner (2002Turner ( , 2003) ) for full details of TBC methods.Current distribution maps of all species are available at www.tucsonbirds.org .
Most, though not all, TBC sites were visited every year from 2001-2003.To maximize spatial coverage while reducing the effects of variation among sites in number of years visited, I used data from two randomly selected years for each site, or one when only one existed).Figure 2 maps the 725 sites surveyed for this study.The current study focuses on eight species associated primarily with desertscrub habitats.I extracted presence/absence data from the total data set with each site counting once, i.e., a species is scored as present at a site if counted on either or both years.This produced a data set of 2794 occurrences of eight species (Table 1).

Environmental data
Several things help to identify factors governing the occurrence of desertscrub birds.First, their natural history is fairly well known (see, e.g., Poole andGill 1992-2002 andreferences therein).For example, each of these species needs desertscrub, and native vegetation can be expected to play a large role in determining where they occur.Second, previous research showed the primacy of native vegetation in determining native bird diversity in humandominated landscapes (e.g., Thomas et al. 1977, Green 1984, see also Germaine et al. 1998).In   Tucson, the total volume of native vegetation retains its status as the strongest correlate of native bird abundance and diversity, despite variation in housing density, exotic vegetation, and other factors (Mills et al. 1989).The present paper seeks more detailed understanding of the relationship between birds and native vegetation cover.
I obtained desertscrub vegetation cover data from a remotely sensed land cover map of the Tucson area produced at the University of Arizona.An initial set of land cover classes was generated through unsupervised classification of 1998 multispectral, i. e., color and near-infrared, high-resolution, i.e., 1 m 2 /pixel, aerial photography of the Tucson area.This was followed by grouping and splitting of classes and additional supervised classification to refine land cover classes.The resulting land cover map comprised 13 land cover classes, e.g., structures, roads, turf grass, thick tree canopy, desertscrub, and was ground-truthed with site visits and by comparison to finer-resolution, i.e., ~0.1 m 2 / pixel, panchromatic imagery.All subsequent analyses use the 'desertscrub' land cover type.
Desertscrub habitat in and around urban areas does not lend itself to the delineation of patches, as often used in studies based on island biogeography theory or patch-based metrics of landscape properties.First, natural desertscrub is not a closed-canopy vegetation type.Subshrubs, grass, and bare ground may intersperse with the densest of desertscrub.Thus, even in protected areas around Tucson, the proportion of an area comprising desertscrub vegetation might not exceed 70%.Second, many realistic approaches to sustaining nature in urban areas will include some degree of non-natural habitat intermingled with natural habitat.The feasibility of such approaches will certainly vary among species, but it makes little sense to limit possible solutions beforehand to contiguous, uniform patches.I thus modeled the response of bird  species to the desertscrub composition of the landscape surrounding each TBC point, calculated as the proportion of pixels within a given radius that are desertscrub.

Modeling species occurrence
This modeling effort seeks to encapsulate the behavior of the system, i.e., bird response to desertscrub cover, in a way that aids biological interpretation and generates quantitative predictions given a particular landscape.My approach involved considering several relevant biological features not generally taken into account by standard statistical models, while retaining the relative ease with which one can parameterize such models.To this end, several considerations serve as a guide for improving a base model.proportion desertscrub cover within 1784 m, i.e., landscape scale.Vertical axis: per-site occurrence probability.For visual clarity, vertical axis scaling varies among species according to maximum predicted occurrence probability.Scaling of two horizontal axes remains constant.Four-letter codes identify species; see Table 1 for species names.ASYMP and BILOG represent the asymptotic and bivariate logistic models, respectively.http://www.ecologyandsociety.org/vol11/iss2/art6/observed occurrence.Recording bird occurrence without error would prove quite expensive, perhaps impossible.So, using realistic survey methods, we expect the probability that we observe a species at a site, given that it is present, to lie below unity.Moreover, although we should generally expect a species to occupy suitable habitat, there is little reason to expect it to occupy all suitable habitat all the time.This observation is among the fundamental contributions of metapopulation ecology (Levins 1969, Hanski 1999).With ordinary logistic regression (Eq.1), occurrence must asymptote at unity.Adding an asymptote parameter α results in the following "asymptote model," which allows non-unity occurrence probability in favorable habitat: (2) Consideration 3: Biological processes acting at multiple spatial scales may influence the presence of a species at a given location (Johnson 1980, Wiens 1989, Levin 1992).Ecologists have increasingly incorporated this consideration into modeling of organism-environment relationships for diverse taxa, including plants (Turner et al. 2003), invertebrates (Steffan-Dewenter et al. 2002), mammals (Schaefer and Messier 1995), and birds (Pearson 1993, Lichstein et al. 2002).In birds, some individual-level behavioral and ecological processes such as nest-site selection and foraging may depend on features present within a relatively small area.But other processes, for example pair formation or population-level processes, may affect the distribution of birds over larger areas.I thus computed desertscrub cover at a "landscape" scale, i.e., 1784 m radius, 10 km 2 area, sufficiently large to accommodate a number of territories of any of the focal species (Poole andGill 1992-2002); and a "local," scale, i.e., 56 m, 1 ha area, which lies nearer the area used by single breeding pairs of these species.Other studies have found factors measured at scales similar to these two to have significant, yet differing apparent impacts on bird species occurrence (e.g., Lichstein et al. 2002, Melles et al. 2003).
Consideration 4: Biological processes do not necessarily interact additively.Instead, they may constrain, or enhance, one another.Liebig's "Law of the Minimum" first introduced this notion (Liebig 1840, as cited in Odum 1950).Huston called for increased attention to the limiting or nonlimiting nature of all factors involved in ecological studies.I allowed for the possibility that processes acting at one scale might constrain the influence of processes acting at another scale by using the product of two logistics, described by the following "constraint model": (3) Ecologists have used similar equations to investigate species responses to multiple environmental variables (van Dam et al. 1986, see also Huisman et al. 1993), but not in a multiscale context.In the constraint model, x 1 and x 2 represent desertscrub cover at local and landscape scales, respectively, and subscripts on β and µ parameters indicate the spatial scale to which the parameters correspond.Because terms in the denominator of Eq. 3 interact multiplicatively, unfavorable values for either environmental variable may restrict occurrence regardless of the value of the other environmental variable.The extent to which this constraint behavior occurs, if at all, depends on the particular values of the β and µ parameters.
I also included a fourth model, a bivariate logistic with interaction term, in analyses for comparison with M CONSTRAINT .Multiple logistic regression sees widespread use in the study of species-environment relationships (e.g., many chapters in Scott et al. 2002).A potential advantage of this approach is that, unlike M CONSTRAINT , the bivariate logistic is linear in its parameters, and thus may be fit using existing statistical methods under the generalized linear model (GLM).The bivariate logistic model (M BILOG ) is described as (4) where µ 3 reflects the strength of the interaction between spatial scales.
For each species, I estimated parameters for each model using maximum likelihood methods.I http://www.ecologyandsociety.org/vol11/iss2/art6/assessed the relative fit of each model to data using the Akaike information criterion (AIC; Akaike 1973).For a given species, the model with the lowest AIC value is that which is most likely after penalization for number of parameters.I conducted all analyses in Matlab v7 (Mathworks 2004).
AIC values allow discrimination of relative performance among candidate models, yet in conservation predictive accuracy may be of equal or greater importance.Receiver operating characteristic (ROC) curves, used in other fields since around 1950 but only more recently in ecology (see, e.g., Fielding andBell 1997, Boyce et al. 2002 andreferences therein), are ideally suited for assessing the predictive accuracy of binary models.In particular, integrating an ROC curve to obtain the area under the curve (AUC) provides a statistic with several important properties (Metz 1978, Swets 1988).Unlike many measures of predictive accuracy, the AUC value is independent of both the overall frequency of occurrence and the threshold probability used for predicting presence.Values of AUC vary from 0.5, i.e., no better than random, to 1, i.e., perfect prediction.I followed Swets (1988, see also Boyce et al. 2002) in considering models with AUC ≥ 0.7 as "useful" applications.
Using the same data to estimate parameters and to test predictive accuracy may overestimate model performance.I used k-fold cross-validation (Fielding and Bell 1997) to assess predictive accuracy.In k-fold cross-validation, the data are divided into k roughly equal partitions.One partition is withheld, and the model is fit to the remaining data.Predictive accuracy is then assessed by calculating AUC using only the partition that was not used in fitting.This procedure is repeated for each of the k partitions, and the k values are averaged to obtain that model's AUC.Each AUC result hereafter represents the result of k-fold crossvalidation (k = 10).

Asymptotic occurrence
Figure 3 shows observed occurrence data for Gambel's Quail as a function of local desertscrub cover, and fitted responses for both M LOGISTIC and M ASYMP .Visual inspection of this and similar plots for the other species indicates that the fixed asymptote of M LOGISTIC model may cause it to fit poorly, whereas the asymptote model better captures the behavior of the response function.Burnham and Anderson (2002:70) suggest that models with an Akaike information criterion (AIC) value within 2 of the minimum AIC value should also be considered as having "substantial" empirical support.Table 2 reports, for each species, AIC values for each model and the lowest AIC among models (AIC best ).We can evaluate the applicability of the asymptote parameter by comparing only M LOGISTIC and M ASYMP .In such a comparison, M ASYMP would best fit the data for seven of eight species and would have substantial support in all eight cases, while M LOGISTIC would best fit the data for one species and would have substantial support for only two species (Table 2).

Multiple spatial scales
The bivariate logistic (M BILOG ) and constraint (M CONSTRAINT ) models included the potentially interactive effects of desertscrub cover at different spatial scales.For the models using habitat data from only one spatial scale (M LOGISTIC and M ASYMP ), using only landscape-scale data did not improve fits over using only local-scale data.Thus, I hereafter show results for local-scale data only with these two models.Of the four models, M CONSTRAINT provided the best fit in nearly all cases, showing lowest AIC values for six of eight species, and having substantial support (AIC-AIC best ≤ 2) for seven of eight (Table 2).Models M ASYMP and M BILOG each fit just one species best and had substantial support for two species (Table 2).
Table 2 additionally gives mean AUC values, across 10 cross-validation partitions, of the best model for each species as a measure of overall model predictive ability.For six of eight species, the best model showed AUC values ≥ 0.7, indicating useful applications of that model.
Visualizing occurrence probabilities predicted by the best models offers insight into how species respond to habitat at different spatial scales.Figure 4 shows predictions of M CONSTRAINT for two species' responses to local habitat for select values of landscape-level habitat.M CONSTRAINT was the best model for these species.Two-dimensional plots such as this allow quantitative examination of http://www.ecologyandsociety.org/vol11/iss2/art6/species' responses.Alternatively, plotting response surfaces in three dimensions readily enables comparisons of qualitative patterns in response among species.Figure 5 shows fitted response surfaces predicted by the best models for each species with respect to both local and landscape desertscrub cover.

Asymptotic occurrence
In principle, the asymptote (M ASYMP ) and constraint (M CONSTRAINT ) models are implicitly tailored to biological phenomena through judicious, a priori choice of additional parameters.Evaluation on Tucson-area desertscrub birds indicates that this effort was justified: based on the Akaike information criterion (AIC) model selection values, M ASYMP outperformed M LOGISTIC , and M CONSTRAINT outperformed all three other models (Table 2).
That M LOGISTIC showed substantial empirical support for no species, and M BILOG for only two species, highlights the inadequacy of these model's assumption that occurrence must reach 100%.The best fit of M LOGISTIC clearly mischaracterizes the response at low and intermediate habitat values as well as asymptotic ones (Fig. 3).If used in such cases as this, standard logistic models may result in erroneous inferences about or predictions of probability of occurrence.Similar errors may potentially arise if logistic regression is used to model occurrence as a function of patch size, a method employed to estimate minimum area requirements (e.g., Robbins et al. 1989, Vickery et al. 1994).

Constraints among spatial scales
Many researchers in recent years have investigated the influences of multiple spatial scales in determining species' spatial distributions (e.g., Pearson 1993, McGarigal and McComb 1995, Schaefer and Messier 1995, Lichstein et al. 2002, Steffan-Dewenter et al. 2002, Melles et al. 2003, Turner et al. 2003).However, studies to date seem not to have investigated the potential for constraintlike interactions among scales as is done here with M CONSTRAINT .
M CONSTRAINT provided the best fit to data and had good predictive accuracy for most species.In the two cases in which M CONSTRAINT did not have the lowest Akaike information criterion (AIC) value, M CONSTRAINT 's reponse surface closely resembles that of the best-fitting model.Indeed, the reason the best model in each of these two cases had lowest AIC values was not because the model fit better, i. e., had lower negative log-likelihood than M CONSTRAINT .Rather, these models had lower AIC values because they fit nearly as well as M CONSTRAINT but did so with fewer parameters, and thus incurred lower AIC parameter penalties.This suggests model M CONSTRAINT may be a useful tool for interpretation of species' responses even in some cases in which it lacks the lowest AIC value.
The variation in response surfaces for best models for each species (Fig. 5) indicates that species do not fall neatly into groups that are simply "sensitive" or "not sensitive" to urban development.Rather, some species are constrained by habitat at the local scale only, e.g., Gila Woodpecker, Curve-billed Thrasher, some by the landscape scale only, e.g., Gilded Flicker, and others by both scales, e.g., Cactus Wren, Gambel's Quail, Verdin, Ashthroated Flycatcher, Pyrrhuloxia.Within these categories, species vary in their sensitivity, i.e., the abruptness of their response to habitat cover.
Several features of response surfaces (Fig. 5) have straightforward biological interpretations.One feature common to nearly all species is a downturn in occurrence at the lowest local cover values.Only one species, Gilded Flicker, does not show this pattern.These findings are consistent with the biology of the species: most desertscrub species rarely occur in golf fairways or other areas with little local desertscrub cover.But flickers are an exception, routinely using open ground and turfgrass in foraging for invertebrates.Given the right landscape, this makes them more likely than other species to be found at sites with little local habitat.
Several species show abrupt changes in occurrence at low landscape cover.For example, once landscape cover reaches a minimum level of 10-15%, Gambel's Quail occur frequently, at up to 89% of sites (Figs. 4 and 5).However, below this minimum, quail occurrence quickly drops to almost 0, regardless of local cover.Species with steep changes with respect to landscape cover vary in occurrence on either side of the change, and in the position of the change along the landscape axis, e. g., Pyrrhuloxia require slightly more landscape cover; Verdins slightly less (Fig. 5).However, the http://www.ecologyandsociety.org/vol11/iss2/art6/message of constraint is clear: for species with such steep changes, no amount of compensation in desertscrub cover at the local scale can overcome the limitations of a poor landscape.This result has several implications for conservation and planning.Creating large, contiguous areas of natural habitat within existing urban areas may prove difficult.Thus, more realistic opportunities for urban conservation likely will include some degree of development intermingled with natural habitat.But could such landscapes actually sustain native species?Consider species such as Gambel's Quail and Pyrrhuloxia, whose abrupt responses in Fig. 5 indicate they are indeed sensitive to habitat cover.These species require a greater proportion of habitat cover, at both scales, than now exists in much of Tucson, where mean desertscrub cover of the central 200 km 2 is 5.8%.However, the habitat cover required for these species to approach natural occurrence rates are relatively low, i.e., roughly 10-15%.Thus, although these two species are currently rare in central Tucson, their restoration to Tucson's developed areas may be possible where desertscrub cover can be raised to a modest 10-15% at both local and landscape scales.Such levels may be attainable with habitat restoration in private yards, gardens, and public spaces such as parks and rights-of-way.Due to the interactions among scales, this restoration must take place at several spatial scales.In particular, large areas of contiguous desertscrub may harbor these species, but they are unlikely to occur in nearby urban sites if local habitat is insufficient.Similarly, any local restoration attempts will meet little success without ensuring adequate habitat at the landscape scale.
If implemented, restoration of additional native species to the developed parts of urban areas could provide several benefits.In a broad conservation sense, restoring native species to developed areas would complement the creation of protected areas to increase the overall habitat available for living things (Rosenzweig et al. 2003).Moreover, this restoration to developed areas would bring native species to the areas in which most of Tucson's residents live, increasing opportunities for people to benefit from or develop an appreciation of nature (Turner et al. 2004).
Although some species might be fairly easy to restore to or sustain in developed areas, results of the best models suggest others species need larger areas and/or lower-density development to persist.
Gilded Flicker and Ash-throated Flycatcher both require high cover at the landscape scale, i.e., 30-50% or more (Fig. 5), an amount of habitat incompatible with many urban human land uses.Conservation efforts for species such as these could be more effectively allocated elsewhere, e.g., lessdeveloped areas or large preserves.
It is important to recognize limitations of the modeling approach used here.The potential for constraint-like interactions among spatial scales is likely to be broadly applicable to other locations, contexts, e.g., urban locations and other landscapes, and taxa.However, parameter estimates for particular species should not be assumed to apply elsewhere without external validation.Biotic and abiotic properties of urban landscapes, as well as ecology and behavior of the focal species, their competitors, and predators may vary among locations.

CONCLUSION
M ASYMP and M CONSTRAINT model relevant factors implicitly rather than explicitly, operating under the assumption that a judiciously chosen, reduced set of variables can act as effective surrogates for biological processes.The simplifying nature of this assumption eases parameterization in comparison with more process-based models, and relies on data, e.g., occurrence, obtainable with relatively modest effort.M CONSTRAINT , in particular, by allowing for the possibility that factors at different scales constrain one another, provides insights into biological phenomena not addressed by conventional statistical approaches.These insights are useful for understanding the distribution of species in complex urban environments, and offer improved predictive ability for evaluating alternative future scenarios for sustaining biodiversity in urban landscapes.Given the observed differences among species, and the intense competing human uses of urban lands, the best models for each species should be used to predict distributions under alternative future development or restoration scenarios.Mapped predictions of these models in planned scenarios would reveal whether particular alternatives include sufficient habitat at relevant spatial scales for each species.

Fig. 2 .
Fig. 2. Map of Tucson study area showing roads a gray lines, interstate highways as divided lines, nature reserves in gray shading, and 725 survey points as dots.Inset shows position of Tucson in western United States.

Fig. 3 .
Fig. 3. Maximum likelihood fits of models M LOGISTIC shown as a dotted line and M ASYMP , a solid line to Tucson Bird Count data for Gambel's Quail (GAQU) response to local desertscrub cover, i.e., proportion of desertscrub cover within 56 m.Open circles: average observed GAQU occurrence in bins of 30 sites having adjacent desertscrub cover values.

Fig. 4 .Fig. 5 .
Fig. 4. Maximum likelihood fits of model M CONSTRAINT reveal differences among species in responses to habitat at different spatial scales.Shown are (a) Gambel's Quail (GAQU) and (b) Ash-throated Flycatcher (ATFL), two species for which M CONSTRAINT was the best model.Each line shows the fitted response to local habitat, i.e., x 1 , the proportion desertscrub cover within 56 m, for a selected value of landscape habitat and x 2 , the proportion desertscrub cover within 1784 m.

Table 1 .
Desertscrub study species and overall frequency of occurrence in Tucson Bird Count data.† Number of sites at which species recorded.‡ Percent of 725 total sites at which species recorded.

Table 2 .
Maximum likelihood parameter estimates and model selection criteria for eight Tucson-area desertscrub species.See Table1for species names.† Akaike Information Criterion.Boldface indicates models considered to have "substantial empirical support," i.e., AIC-AIC best ≤ 2. ‡ Lowest AIC value among models, indicating highest likelihood given the data after penalization for number of parameters.§ Mean area under the receiver operating curve (AUC) based on 10-fold cross-validation of the best model.M LOGISTIC, M ASYMP, M BILOG, and M CONSTRAINT represent the logistic, asymptote, bivariate logistic, and constraint models, respectively.