MACROECOLOGY OF BREEDING BIRDS OF NEW YORK STATE. INFLUENCES OF
CLIMATE CHANGE, LANDSCAPE MATRIX, AND SPATIAL SCALE.
By
Marta Anna Jarzyna

A DISSERTATION
Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degrees of
Fisheries and Wildlife – Doctor of Philosophy
Ecology, Evolutionary Biology and Behavior – Dual Major
2014

ABSTRACT
MACROECOLOGY OF BREEDING BIRDS OF NEW YORK STATE. INFLUENCES OF
CLIMATE CHANGE, LANDSCAPE MATRIX, AND SPATIAL SCALE.
By
Marta Anna Jarzyna
Climate change has the potential to greatly affect biodiversity and there is convincing evidence
that most taxonomic groups are already responding to the recent climate warming. Even though
the impacts of climate change on biodiversity are felt across a wide range of ecosystems, the
severity and sometimes the direction of these impacts are likely to differ across space because
factors other than climate will influence how biodiversity responds to climate change. For
example, composition and configuration of land cover may play significant roles in the relative
vulnerabilities of biodiversity to changing climate. However, despite some initial investigations,
empirical evidence for the influence of land cover is lacking because studies that integrate both
land cover and climate are rare. The overarching goal of my dissertation was to assess the
relative roles of climate change and land cover in the observed changes in avian biodiversity.
Specifically, I was interested in the influences of land-cover composition and configuration on
relative vulnerabilities of different avian groups to climate change and spatial scales at which
such influences are relevant. Chapter 1 focuses on testing the utility of the spatially-varying
coefficients (SVC) model to quantify the influence of non-stationarity on relationships between
temporal community change in avian assemblages and environmental covariates. Chapter 2
investigates bioclimatic relationships of grassland and forest breeding birds across varying
gradients of grassland and forest habitat. Chapter 3 investigates the relationship between
temporal changes in avian assemblages and the interaction of climate change and land-cover
fragmentation. Chapter 4 explores scale-dependence of temporal changes in avian communities

and investigates relevant environmental drivers of the community change at different spatial
scales.
The collective works in these chapters contribute four primary conclusions for better
understanding of the implications of the interaction of land use and climate change to avian
diversity. First, relationships between changes in biodiversity and environmental factors are
often spatially non-stationary, i.e., the relationship between a response variable and the predictor
covariates varies across the spatial extent of the study. Second, the amount of suitable land cover
can significantly alter species responses to climate change, but this effect of land cover depends
on the type of habitat. Third, impacts of climate change on ecological communities are more
pronounced in regions of unfragmented landscapes than in locations with fragmented habitats.
Fourth, temporal changes in community composition are scale-dependent and so are the
mechanisms driving these changes. Specifically, climate change operates at smaller spatial scales
than landscape fragmentation.
The implications of these findings suggest that current conservation strategies may be
insufficient to protect biodiversity in the face of climate change. Existing biodiversity
conservation generally focuses on conserving large regions with undisturbed and contiguous
habitats that generally support high species diversity. My results suggest that ecological
communities of such habitats will undergo the most drastic compositional changes as a result of
climate change. Thus, it is critical that the existing conservation strategies are placed in the
context of climate change. My work suggests that successful biodiversity conservation needs to
consider both the individual species’ vulnerabilities to climate change based on their habitat
association or life-history strategies and relative vulnerabilities of entire communities based on
the landscape in which these communities persist.

Copyright by
MARTA ANNA JARZYNA
2014

ACKNOWLEDGMENTS

I thank my advisor, Dr. Porter, for his excellent guidance, caring, patience, and providing me
with a superb atmosphere for doing research. I also thank my committee members, Drs. Finley,
Maurer, and Zuckerberg for their support and guidance throughout the entire process. Dr. Finley
patiently taught me hierarchical Bayesian modeling approaches and their applications in a spatial
context. Dr. Maurer introduced me to the fascinating world of macroecology, for which I could
not be more grateful. Dr. Zuckerberg provided excellent suggestions regarding climate change
implications for biodiversity.
I thank the faculty and students in the Quantitative Wildlife Laboratory that provided
insight and support, especially Dr. Snow, who was not only a fair critic of my work but also a
great friend. This research was funded in part by the National Aeronautics and Space
Administration, the Boone and Crockett Club, and the Graduate School at Michigan State
University. I also thank Michigan State University for support.
Last but not least, I am grateful for the love and support from my family and friends.
Special thanks go to my husband, Jim, without whose love and encouragement I would not have
completed this work.

v

TABLE OF CONTENTS

LIST OF TABLES ......................................................................................................................... ix
LIST OF FIGURES ...................................................................................................................... xii
PROLOGUE ................................................................................................................................... 1
ACCOUNTING FOR THE SPACE-VARYING NATURE OF THE RELATIONSHIPS
BETWEEN TEMPORAL COMMUNITY TURNOVER AND THE ENVIRONMENT ............. 6
ABSTRACT ................................................................................................................................ 6
1.1 INTRODUCTION................................................................................................................. 7
1.2 DATA AND METHODS .................................................................................................... 10
1.2.1 Site description ............................................................................................................. 10
1.2.2 Breeding Bird Atlas ...................................................................................................... 11
1.2.3 Temporal turnover ........................................................................................................ 12
1.2.4 Model covariates........................................................................................................... 13
1.2.4.1 Climatic trends ....................................................................................................... 13
1.2.4.2 Habitat fragmentation ............................................................................................ 14
1.2.4.3 Survey effort .......................................................................................................... 16
1.2.5 Statistical analysis......................................................................................................... 16
1.3 RESULTS............................................................................................................................ 20
1.3.1 Temporal turnover ........................................................................................................ 20
1.3.2 Covariates ..................................................................................................................... 22
1.3.3 Statistical analysis......................................................................................................... 22
1.3.3.1 Model fit and predictive ability.............................................................................. 22
1.3.3.2 Coefficient estimates .............................................................................................. 25
1.4 DISCUSSION ..................................................................................................................... 28
1.5 ACKNOWLEDGMENTS................................................................................................... 33
SYNERGISTIC EFFECTS OF CLIMATE AND LAND COVER: ARE GRASSLAND BIRDS
MORE VULNERABLE TO CLIMATE CHANGE? ................................................................... 34
ABSTRACT .............................................................................................................................. 34
2.1 INTRODUCTION............................................................................................................... 35
2.2 METHODS.......................................................................................................................... 38
2.2.1 Breeding Bird Atlas and chosen species....................................................................... 38
2.2.2 Model covariates........................................................................................................... 40
2.2.3 Habitat amount ............................................................................................................. 40
2.2.4 Statistical analysis......................................................................................................... 41
2.3 RESULTS............................................................................................................................ 43
2.3.1 Model comparison ........................................................................................................ 43
2.3.2 Bioclimatic relationships .............................................................................................. 44
2.3.3 Heterogeneity of bioclimatic relationships across habitat amount ............................... 49
2.4 DISCUSSION ..................................................................................................................... 54
2.4.1 Bioclimatic relationships .............................................................................................. 54
vi

2.4.2 Heterogeneity of bioclimatic relationships across habitat amount ............................... 55
2.5 ACKNOWLEDGMENTS................................................................................................... 59
LANDSCAPE FRAGMENTATION AFFECTS RESPONSES OF AVIAN COMMUNITIES TO
CLIMATE CHANGE ................................................................................................................... 61
ABSTRACT .............................................................................................................................. 61
3.1 INTRODUCTION............................................................................................................... 62
3.2 METHODS.......................................................................................................................... 65
3.2.1 Site description ............................................................................................................. 65
3.2.2 Breeding Bird Atlas ...................................................................................................... 66
3.2.3 Temporal turnover ........................................................................................................ 67
3.2.4 Habitat fragmentation ................................................................................................... 68
3.2.5 Climatic trends .............................................................................................................. 69
3.2.6 Survey effort ................................................................................................................. 70
3.2.7 Statistical analysis......................................................................................................... 70
3.3 RESULTS............................................................................................................................ 73
3.3.1 Relationships with climate change and habitat fragmentation ..................................... 73
3.3.1.1 Temporal turnover ................................................................................................. 73
3.3.1.2 Extinction ............................................................................................................... 79
3.3.1.3 Colonization ........................................................................................................... 83
3.3.2 Responses of different migratory groupings ................................................................ 85
3.3.2.1 Temporal turnover ................................................................................................. 85
3.3.2.2 Extinction ............................................................................................................... 87
3.3.2.3 Colonization ........................................................................................................... 89
3.4 DISCUSSION ..................................................................................................................... 91
3.4.1 Relationships with climate change and habitat fragmentation ..................................... 91
3.4.2 Responses of different migratory groupings ................................................................ 93
3.4.3 Conclusions .................................................................................................................. 95
3.5 ACKNOWLEDGMENTS................................................................................................... 95
SCALE-DEPENDENCE OF TEMPORAL CHANGES IN AVIAN COMMUNITIES ............. 97
ABSTRACT .............................................................................................................................. 97
4.1 INTRODUCTION............................................................................................................... 98
4.2 METHODS........................................................................................................................ 100
4.2.1 Breeding Bird Atlas .................................................................................................... 100
4.2.2 Spatial scaling of community change ......................................................................... 101
4.2.3 Factors influencing community change across spatial grains .................................... 102
4.3 RESULTS.......................................................................................................................... 106
4.3.1 Spatial scaling of community change ......................................................................... 106
4.3.2 Factors influencing community change across spatial grains .................................... 113
4.3.2.1 Temporal turnover ............................................................................................... 113
4.3.2.2 Extinction ............................................................................................................. 113
4.3.2.3 Colonization ......................................................................................................... 114
4.4 DISCUSSION ................................................................................................................... 117
4.4.1 Spatial scaling of community change ......................................................................... 117
4.4.2 Factors influencing community change across spatial grains .................................... 119
vii

4.4.3 Conclusions ................................................................................................................ 124
4.5 ACKNOWLEDGMENTS................................................................................................. 124
EPILOGUE ................................................................................................................................. 125
APPENDIX ................................................................................................................................. 130
LITERATURE CITED ............................................................................................................... 162

viii

LIST OF TABLES

Table 1.1 Comparison of the non-spatial, spatially-varying intercept (SVI), and spatially-varying
coefficient (SVC) models for four migratory groupings using Deviance Information Criterion
(DIC), estimated number of parameters pD, and Relative Mean Square Error (RMSE). Predictive
performance of the models was evaluated by calculating RMSEpred of the observed values from
the hold-out set and the predicted values resulting from all three models………………….…....23
Table 1.2 Coefficient estimates resulting from the non-spatial and spatially-varying intercept
(SVI) models for all model covariates: intercept (β0 ), magnitude of the 25-year (1980-2005)
trend in average maximum temperature of the breeding season (βTMAX), magnitude of the 25-year
(1980-2005) trend in average minimum temperature of the breeding season (βTMIN), magnitude of
the 25-year (1980-2005) trend in average total precipitation of the breeding season (βPRECIP),
percent developed land (βDEVEL), edge density (βED), and effort (βEFF). 50% indicates the mean of
the posterior distribution, while 2.5% and 97.5% are the lower and upper quantiles of the
posterior distribution…....……………………………………………………………………..…25
Table 2.1 Coefficient estimates resulting from the spatially-varying intercept (SVI) models for
the intercept (INTERCEPT), 2000-05 average maximum temperature of the breeding season
(TMAX), 2000-05 average minimum temperature of the breeding season (TMIN), 2000-05
average total precipitation of the breeding season (PRECIP), and 2000-05 survey effort (EFF) for
grassland and forest birds. Species occurrence data were retrieved from the 2000-05 New York
State Breeding Bird Atlas………………………………………………………………..………45
Table 2.2 Relationships between amount of habitat (percent of grassland cover, GPLAND, and
percent of forest cover, FPLAND, for grassland and forest birds, respectively) and spatiallyvarying coefficient estimates for 2000-05 average maximum temperature of the breeding season
(TMAX), 2000-05 average minimum temperature of the breeding season (TMIN), and 2000-05
average total precipitation of the breeding season (PRECIP). The “*” symbol indicates that less
than 10% of the spatially-varying coefficient estimates (i.e., TMAX, TMIN, PRECIP coefficient
estimates resulting from the SVC models) had the 95% credible intervals not overlapping zero.
Species occurrence data were retrieved from the 2000-05 New York State Breeding Bird
Atlas………………………………………………………………………………...………..…..51
Table 3.1 Comparison of five competing models of temporal turnover for all four avian
groupings (i.e., all species, long-distance migrants, short-distance migrants, and resident
species). Models included the following covariates: magnitude of the 25-year (1980-2005) trend
in average maximum temperature of the breeding season (TMAX), magnitude of the 25-year
(1980-2005) trend in average minimum temperature of the breeding season (TMIN), magnitude
of the 25-year (1980-2005) trend in average total precipitation of the breeding season (PRECIP),
edge density (ED), percent developed land (DEVEL), interaction between the magnitude of the
25-year (1980-2005) trend in average maximum temperature of the breeding season and edge
density (TMAX*ED), interaction between the magnitude of the 25-year (1980-2005) trend in
average minimum temperature of the breeding season and edge density (TMIN*ED), interaction
between the magnitude of the 25-year (1980-2005) trend in average total precipitation of the
breeding season and edge density (PRECIP*ED), and survey effort (EFF). Model comparison
ix

was done using Deviance Information Criterion (DIC); pD indicates the effective number of
parameters………………………………………………………….…………………………….75
Table 3.2 Comparison of five competing models of extinction for all four avian groupings (i.e.,
all species, long-distance migrants, short-distance migrants, and resident species). Models
included the following covariates: magnitude of the 25-year (1980-2005) trend in average
maximum temperature of the breeding season (TMAX), magnitude of the 25-year (1980-2005)
trend in average minimum temperature of the breeding season (TMIN), magnitude of the 25-year
(1980-2005) trend in average total precipitation of the breeding season (PRECIP), edge density
(ED), percent developed land (DEVEL), interaction between the magnitude of the 25-year
(1980-2005) trend in average maximum temperature of the breeding season and edge density
(TMAX*ED), interaction between the magnitude of the 25-year (1980-2005) trend in average
minimum temperature of the breeding season and edge density (TMIN*ED), interaction between
the magnitude of the 25-year (1980-2005) trend in average total precipitation of the breeding
season and edge density (PRECIP*ED), and survey effort (EFF). Model comparison was done
using Deviance Information Criterion (DIC); pD indicates the effective number of
parameters……………………………………………………………………………………......81
Table 3.3 Comparison of five competing models of colonization for all four avian groupings
(i.e., all species, long-distance migrants, short-distance migrants, and resident species). Models
included the following covariates: magnitude of the 25-year (1980-2005) trend in average
maximum temperature of the breeding season (TMAX), magnitude of the 25-year (1980-2005)
trend in average minimum temperature of the breeding season (TMIN), magnitude of the 25-year
(1980-2005) trend in average total precipitation of the breeding season (PRECIP), edge density
(ED), percent developed land (DEVEL), interaction between the magnitude of the 25-year
(1980-2005) trend in average maximum temperature of the breeding season and edge density
(TMAX*ED), interaction between the magnitude of the 25-year (1980-2005) trend in average
minimum temperature of the breeding season and edge density (TMIN*ED), interaction between
the magnitude of the 25-year (1980-2005) trend in average total precipitation of the breeding
season and edge density (PRECIP*ED), and survey effort (EFF). Model comparison was done
using Deviance Information Criterion (DIC); pD indicates the effective number of
parameters……..............................................................................................................................84
Table A 1.1 List of all bird species (N = 256) from the New York State Breeding Bird Atlas
(BBA)…………………….…...………………….…………………………………………......131
Table A 2.1 Comparison of the spatially-varying intercept (SVI) and spatially-varying
coefficients (SVC) models for grassland and forest birds recorded during the 2000-05 New York
State Breeding Bird Atlas………………………………………………………………………138
Table A 3.1 Coefficient estimates (βs) of each covariate resulting from all tested models of
temporal turnover (starting with the best model as indicated by DIC) for four avian groupings: all
species, long-distance migrants, short-distance migrants, and resident species. Covariates
included intercept (INT), magnitude of the 25-year (1980-2005) trend in average maximum
temperature of the breeding season (TMAX), magnitude of the 25-year (1980-2005) trend in
average minimum temperature of the breeding season (TMIN), magnitude of the 25-year (1980x

2005) trend in average total precipitation of the breeding season (PRECIP), edge density (ED),
percent developed land (DEVEL), interaction between the magnitude of the 25-year (1980-2005)
trend in average maximum temperature of the breeding season and edge density (TMAX*ED),
interaction between the magnitude of the 25-year (1980-2005) trend in average minimum
temperature of the breeding season and edge density (TMIN*ED), interaction between
magnitude of the 25-year (1980-2005) trend in average total precipitation of the breeding season
and edge density (PRECIP*ED), and survey effort (EFF). 50% indicates the mean of the
posterior distribution, while 2.5% and 97.5% are the lower and upper quantiles of the posterior
distribution……………………………………………………………………………………...140
Table A 3.2 Coefficient estimates (βs) of each covariate resulting from all tested models of
extinction (starting with the best model as indicated by DIC) for four avian groupings: all
species, long-distance migrants, short-distance migrants, and resident species. Covariates
included intercept (INT), magnitude of the 25-year (1980-2005) trend in average maximum
temperature of the breeding season (TMAX), magnitude of the 25-year (1980-2005) trend in
average minimum temperature of the breeding season (TMIN), magnitude of the 25-year (19802005) trend in average total precipitation of the breeding season (PRECIP), edge density (ED),
percent developed land (DEVEL), interaction between the magnitude of the 25-year (1980-2005)
trend in average maximum temperature of the breeding season and edge density (TMAX*ED),
interaction between the magnitude of the 25-year (1980-2005) trend in average minimum
temperature of the breeding season and edge density (TMIN*ED), interaction between
magnitude of the 25-year (1980-2005) trend in average total precipitation of the breeding season
and edge density (PRECIP*ED), and survey effort (EFF). 50% indicates the mean of the
posterior distribution, while 2.5% and 97.5% are the lower and upper quantiles of the posterior
distribution……………………………………………………………………………………...145
Table A 3.3 Coefficient estimates (βs) of each covariate resulting from all tested models of
colonization (starting with the best model as indicated by DIC) for four avian groupings: all
species, long-distance migrants, short-distance migrants, and resident species. Covariates
included intercept (INT), magnitude of the 25-year (1980-2005) trend in average maximum
temperature of the breeding season (TMAX), magnitude of the 25-year (1980-2005) trend in
average minimum temperature of the breeding season (TMIN), magnitude of the 25-year (19802005) trend in average total precipitation of the breeding season (PRECIP), edge density (ED),
percent developed land (DEVEL), interaction between the magnitude of the 25-year (1980-2005)
trend in average maximum temperature of the breeding season and edge density (TMAX*ED),
interaction between the magnitude of the 25-year (1980-2005) trend in average minimum
temperature of the breeding season and edge density (TMIN*ED), interaction between
magnitude of the 25-year (1980-2005) trend in average total precipitation of the breeding season
and edge density (PRECIP*ED), and survey effort (EFF). 50% indicates the mean of the
posterior distribution, while 2.5% and 97.5% are the lower and upper quantiles of the posterior
distribution……………………………………………………………………………………...150

xi

LIST OF FIGURES

Figure 1.1 Observed temporal turnover in communities of breeding birds in New York
expressed using Diamond-May index (A) and the fitted temporal turnover values resulting from
the non-spatial model (B), spatially-varying intercept (SVI) model (C), and spatially-varying
coefficients (SVC) model D).…………..………………………………………………………..21
Figure 1.2 Observed temporal turnover of the holdout data set (A) compared with the surfaces of
the predicted values resulting from the non-spatial model (B), spatially-varying intercept (SVI)
model (C), and spatially-varying coefficients (SVC) model (D)……………………..…..……...24
Figure 1.3 Coefficient estimates resulting from the spatially-varying coefficient (SVC) models
of (A) the intercept, (B) the magnitude of the 25-year (1980-2005) trend in average maximum
temperature of the breeding season (βTMAX), (C) the magnitude of the 25-year (1980-2005) trend
in average minimum temperature of the breeding season (βTMIN), (D) the magnitude of the 25year (1980-2005) trend rates in average total precipitation of the breeding season (βPRECIP), (E)
percent developed land (βDEVEL), (F) edge density (βED), and (G) effort (βEFF). Black point
symbols identify grid cells with posterior 90% credible intervals that include zero. The surface
colors identify the direction of the regression coefficients' sign……………………………..….27
Figure 2.1 Mean values of the coefficient estimates of 2000-05 average spring maximum
temperature (TMAX), 2000-05 average spring minimum temperature (TMIN), and 2000-05
average total spring precipitation (PRECIP) for both avian groups: grassland birds (circles) and
forest birds (triangles). 95% confidence intervals indicate whether the mean value of the
coefficient estimates is significantly different from zero. On each plot from left to right:
American kestrel, northern harrier, eastern bluebird, eastern kingbird, horned lark, grasshopper
sparrow, savannah sparrow, vesper sparrow, bobolink, brown-headed cowbird, eastern
meadowlark, blue-gray gnatcatcher, blue-headed vireo, golden-crowned kinglet, great crested
flycatcher, red-breasted nuthatch, tufted titmouse, veery, black-and-white warbler, black-throated
blue warbler, black-throated green warbler, blackburnian warbler, and Canada warbler. Species
occurrence data were retrieved from the 2000-05 New York State Breeding Bird Atlas…….…48
Figure 2.2 Mean values of the associations between the percent habitat (percent grassland cover,
GPLAND, and percent forest cover, FPLAND, for grassland and forest species, respectively) and
coefficient estimates of 2000-05 average spring maximum temperature (TMAX), 2000-05
average spring minimum temperature (TMIN), and 2000-05 average total spring precipitation
(PRECIP) for both avian groups: grassland birds (circles) and forest birds (triangles). 95%
confidence intervals indicate whether the mean value of the associations is significantly different
from zero. On each plot from left to right: American kestrel, northern harrier, eastern bluebird,
eastern kingbird, horned lark, grasshopper sparrow, savannah sparrow, vesper sparrow, bobolink,
brown-headed cowbird, eastern meadowlark, blue-gray gnatcatcher, blue-headed vireo, goldencrowned kinglet, great crested flycatcher, red-breasted nuthatch, tufted titmouse, veery, blackand-white warbler, black-throated blue warbler, black-throated green warbler, blackburnian
warbler, and Canada warbler. Species occurrence data were retrieved from the 2000-05 New
York State Breeding Bird Atlas………………………………………………………………….53
xii

Figure 3.1 Temporal turnover, extinction, and colonization observed between 1980-85 and
2000-05 for communities of (A) all recorded species regardless of their migratory status, (B)
long-distance migrants, (C) short-distance migrants, and (D) resident birds. Temporal turnover
was calculated as a proportion of all species gained or lost between 1980-85 and 2000-05 in a
particular site (i.e., Breeding Bird Atlas block) relative to all species recorded across both time
periods; extinction was calculated as a proportion of species lost between 1980-85 and 2000-05
in a particular site relative to all species present in a block in 1980-85; colonization was
calculated as a proportion of species gained between 1980-85 and 2000-05 in a particular site
relative to all species recorded across both time periods. High values of all three metrics indicate
high temporal turnover, extinction, and colonization…...…………………………………….....74
Figure 3.2 Coefficient estimates (i.e., mean values of the posterior distribution and associated
credible intervals) of the best model for each avian grouping for (A) temporal turnover, (B)
extinction, and (C) colonization. Abbreviations for the covariates are as follows: magnitude of
the 25-year (1980-2005) trend in average maximum temperature of the breeding season
(TMAX), magnitude of the 25-year (1980-2005) trend in average minimum temperature of the
breeding season (TMIN), magnitude of the 25-year (1980-2005) trend in average total
precipitation of the breeding season (PRECIP), edge density (ED), percent developed land
(DEVEL), interaction between the magnitude of the 25-year (1980-2005) trend in average
maximum temperature of the breeding season and edge density (TMAX-ED), interaction
between the magnitude of the 25-year (1980-2005) trend in average minimum temperature of the
breeding season and edge density (TMIN-ED), interaction between magnitude of the 25-year
(1980-2005) trend in average total precipitation of the breeding season and edge density
(PRECIP-ED)…………………………………………………………………...……………......77
Figure 3.3 Effects size plots resulting from the best temporal turnover models for (A) the
interaction of magnitude of the 25-year (1980-2005) trend in average minimum temperature of
the breeding season (TMIN) and edge density (ED) for all species, (B) the interaction of
magnitude of the 25-year (1980-2005) trend in average minimum temperature of the breeding
season (TMIN) and edge density (ED) for long-distance migrants, (C) the interaction of
magnitude of the 25-year (1980-2005) trend in total precipitation of the breeding season
(PRECIP) and edge density (ED) for long-distance migrants, (D) the interaction of magnitude of
the 25-year (1980-2005) trend in average maximum temperature of the breeding season (TMAX)
and edge density (ED) for resident birds, and (E) the interaction of magnitude of the 25-year
(1980-2005) trend in average minimum temperature of the breeding season (TMIN) and edge
density (ED) for resident birds. The interaction terms shown here are the ones whose credible
intervals did not overlap zero. Green points are the observed data points and are plotted to
indicate levels of confidence for the predicted surface……………………………………..……78
Figure 3.4 Effects size plots resulting from the best extinction models for (A) the interaction of
magnitude of the 25-year (1980-2005) trend in average minimum temperature of the breeding
season (TMIN) and edge density (ED) for all species, (B) the interaction of magnitude of the 25year (1980-2005) trend in average minimum temperature of the breeding season (TMIN) and
edge density (ED) for long-distance migrants, (C) the interaction of magnitude of the 25-year
(1980-2005) trend in total precipitation of the breeding season (PRECIP) and edge density (ED)
for long-distance migrants, (D) the interaction of magnitude of the 25-year (1980-2005) trend in
xiii

average minimum temperature of the breeding season (TMIN) and edge density (ED) for shortdistance migrants, and (E) the interaction of magnitude of the 25-year (1980-2005) trend in
average minimum temperature of the breeding season (TMIN) and edge density (ED) for
resident birds. The interaction terms shown here are the ones whose credible intervals did not
overlap zero. Green points are the observed data points and are plotted to indicate levels of
confidence for the predicted surface……………………………………………………..………82
Figure 4.1 Spatial scaling of temporal turnover, extinction, and colonization observed between
1980-85 and 2000-05 for avian communities in New York State across five different spatial
grains (i.e., 5x5 km, 10x10 km, 20x20 km, 40x40 km, and 80x80 km). Temporal turnover was
calculated as a proportion of all species gained or lost between 1980-85 and 2000-05 in a
particular site (i.e., 5x5 km Breeding Bird Atlas block or a scaled-up version of it) relative to all
species recorded across both time periods; extinction was calculated as a proportion of species
lost between 1980-85 and 2000-05 in a particular site relative to all species present in a site in
1980-85; colonization was calculated as a proportion of species gained between 1980-85 and
2000-05 in a particular site relative to all species recorded across both time periods. High values
of all three metrics indicate high temporal turnover, extinction, and colonization.……………107
Figure 4.2 Patterns of temporal turnover observed between 1980-85 and 2000-05 for avian
communities in New York State across five different spatial grains (i.e., 5x5 km, 10x10 km,
20x20 km, 40x40 km, and 80x80 km). Temporal turnover was calculated as a proportion of all
species gained or lost between 1980-85 and 2000-05 in a particular site (i.e., 5x5 km Breeding
Bird Atlas block or a scaled-up version of it) relative to all species recorded across both time
periods. High values indicate high temporal turnover………………………………..……...…108
Figure 4.3 Patterns of extinction observed between 1980-85 and 2000-05 for avian communities
in New York State across five different spatial grains (i.e., 5x5 km, 10x10 km, 20x20 km, 40x40
km, and 80x80 km). Extinction was calculated as a proportion of species lost between 1980-85
and 2000-05 in a particular site (i.e., 5x5 km Breeding Bird Atlas block or a scaled-up version of
it) relative to all species present in a site in 1980-85. High values indicate high proportion of
extinction events……………………………………………………..…………………………111
Figure 4.4 Patterns of colonization observed between 1980-85 and 2000-05 for avian
communities in New York State across five different spatial grains (i.e., 5x5 km, 10x10 km,
20x20 km, 40x40 km, and 80x80 km). Colonization was calculated as a proportion of species
gained between 1980-85 and 2000-05 in a particular site (i.e., 5x5 km Breeding Bird Atlas block
or a scaled-up version of it) relative to all species recorded across both time periods. High values
indicate high proportion of colonization events………………………………………………...112
Figure 4.5 Coefficient estimates (i.e., mean values of the posterior distribution and associated
credible intervals) resulting from the models for (A) temporal turnover, (B) extinction, and (C)
colonization. Abbreviations for the covariates are as follows: magnitude of the 25-year (19802005) trend in average maximum temperature of the breeding season (TMAX), magnitude of the
25-year (1980-2005) trend in average minimum temperature of the breeding season (TMIN),
magnitude of the 25-year (1980-2005) trend in average total precipitation of the breeding season
(PRECIP), edge density (ED), percent developed land (DEVEL), and elevation (ELEV)…….116
xiv

Figure A 1.1 Spatial distribution of (A) magnitude of the 25-year (1980-2005) trend in average
maximum temperatures [ºC], (B) magnitude of the 25-year (1980-2005) trend in average
minimum temperatures [ºC], (C) magnitude of the 25-year (1980-2005) trend in total
precipitation [mm], (D) percent developed land [%], (E) edge density [m/ha], and (F) survey
effort [h]……………………………………………………………...…………………...…….155
Figure A 1.2 Deviations of the spatially-varying coefficient estimates from the global mean for
(A) the intercept, (B) the magnitude of the 25-year (1980-2005) trend in average maximum
temperature of the breeding season (βTMAX), (C) the magnitude of the 25-year (1980-2005) trend
in average minimum temperature of the breeding season (βTMIN), (D) the magnitude of the 25year (1980-2005) trend rates in average total precipitation of the breeding season (βPRECIP), (E)
percent developed land (βDEVEL), (F) edge density (βED), and (G) effort (βEFF)…...……..…….156
Figure A 2.1 The relationships between the amount of habitat (GPLAND and FPLAND for
grassland and forest birds, respectively; left-hand panel of each column) and spatially-varying
coefficient estimates (right-hand panel of each column) of 2000-05 average spring maximum
temperature (TMAX, left column), 2000-05 average spring minimum temperature (TMIN,
middle column), and 2000-05 average total spring precipitation (PRECIP, right column) for
[from top to bottom] American kestrel (AMKE), northern harrier (NOHA), eastern bluebird
(EABL), eastern kingbird (EAKI), horned lark (HOLA), grasshopper sparrow (GRSP), savannah
sparrow (SAVS), vesper sparrow (VESP), bobolink (BOBO), brown-headed cowbird (BHCO),
eastern meadowlark (EAME)], blue-gray gnatcatcher (BGGN), blue-headed vireo (BHVI),
golden-crowned kinglet (GCKI), great crested flycatcher (GCFL), red-breasted nuthatch
(RBNU), tufted titmouse (ETTI), veery (VEER), black-and-white warbler (BAWW), blackthroated blue warbler (BTBW), black-throated green warbler (BTNW), blackburnian warbler
(BLBW), and Canada warbler (CAWA). Dark blue symbols represent coefficient estimates with
95% credible intervals not overlapping zero in locations where species was detected; black
symbols represent coefficient estimates with 95% credible intervals not overlapping zero in
locations where species was not detected; light blue points represent with 95% credible intervals
overlapping zero coefficient estimates in locations where species was not detected; grey points
represent coefficient estimates with 95% credible intervals overlapping zero in locations where
the species was not detected. Black dashed line represents the slope of the relationship between
the coefficient estimates and GPLAND and FPLAND. The right-hand panel shows coefficient
estimates of TMAX, TMIN, and PRECIP resulting from the spatially-varying coefficients
models. Grayed-out locations indicate coefficient estimates with 95% credible intervals
overlapping zero. Species occurrence data were retrieved from the 2000-05 New York State
Breeding Bird Atlas.……………………………………………………...……….……………157
Figure A 3.1 Model validation using plots of predicted vs observed values for temporal turnover,
extinction, and colonization observed between 1980-85 and 2000-05 for communities of (A) all
recorded species regardless of their migratory status, (B) long-distance migrants, (C) shortdistance migrants, and (D) resident birds. Model validation was done for the best model
only……………………………………………………………………………………………..161

xv

PROLOGUE

Climate change has potential to greatly affect biodiversity and there is convincing evidence that
most taxonomic groups are already responding to the recent climate warming (e.g., Jarzyna et al.
2013). Among those ecological responses most often reported are shifting population sizes and
distributions (Jiguet et al. 2010, Zuckerberg et al. 2009), changes in breeding performance
(Ahola et al. 2009, Husby et al. 2009), and changes in the timing of breeding (Root et al. 2003,
Dunn and Winkler 2010), migration (Jones and Cresswell 2010, Tottrup et al. 2010), and
hibernation (Inouye et al. 2000). However, even though the impacts of climate change on
biodiversity are felt across a wide range of ecosystems, the severity and sometimes the direction
of these impacts are bound to differ across space because factors other than climate will
influence how biodiversity responds to climate change. Land cover is among the main factors
with the potential to alter biodiversity’s responses to climate change.
National Aeronautics and Space Administration (NASA) became interested in assessing
the roles of recent climatic and land cover changes in observed landscape-scale shifts in
biodiversity and an evaluation of the relative contributions of each as well as their interaction.
Indeed, research conducted to date suggests that land cover plays a significant role in the relative
vulnerabilities of biodiversity to changing climate. For example, simulation studies have shown
that habitat loss reduces species ability to adjust to changing climatic conditions and thus
exacerbate the detrimental effects of climate change (Travis 2003). Similarly, an increase in the
amount of habitat has been shown to slow down species extinction rates associated with
changing climate (Jeltsch et al. 2011). Habitat fragmentation is also thought to influence species
responses to climate change by affecting the rate of distributional shifts, allowing species in less
fragmented habitats to disperse faster and farther (Opdam and Wascher 2004). However, despite
1

these initial investigations, empirical evidence is lacking because studies that integrate both land
cover and climate are rare.
The overarching goal of my dissertation was to assess the relative roles of climate change
and land cover in the observed changes in avian biodiversity. Specifically, I was interested in the
influences of land-cover composition and configuration on relative vulnerabilities of different
avian groups to climate change and spatial scales at which such influences are relevant. I used
the 1980 and 2000 New York State Breeding Bird Atlas. All chapters for this research were
conducted with approval of the Michigan State University, Institutional Animal Care and Use
Committee for use of digital records (approval date 11-01-2013).
Chapter 1 focuses on testing the utility of the spatially-varying coefficients (SVC) model
to quantify the influence of non-stationarity on relationships between temporal community
change in avian assemblages and changes in climate and land cover. This chapter has been
accepted for the publication in Ecography. Spatial non-stationarity is present when the strength
and nature of the relationship between a response variable and the predictor covariates vary
across the spatial domain (Fortin and Dale 2009, Miller 2012). This spatial variation might be
produced by interactions or feedbacks among covariates or presence of inherently different
mechanisms that impact the outcome in different parts of the study region (Miller 2012). Given
the ubiquity of space-varying relationships in nature (e.g., Bini et al. 2009), it is often unrealistic
to assume that relationships between biodiversity and changes in climate and land cover can be
captured by a set of stationary regression coefficients (i.e., coefficients that are constant across
space). My objective is to compare inferences drawn from the SVC models to that obtained from
space-varying intercept models (i.e., models that account for spatial autocorrelation, but not
spatial non-stationarity) and also models that do not acknowledge any spatial structure beyond

2

what is introduced by the covariates. If the SVC models are superior to the other two tested
models, they will be incorporated in the subsequent chapters to account for spatial nonstationarity of ecological relationships.
Chapter 2 builds on Chapter 1 by incorporating SVC models to investigate bioclimatic
relationships of grassland and forest breeding birds across varying gradients of grassland and
forest habitat. Research suggests that the quantity and quality of available habitat might be
important in shaping bioclimatic relationships of species and thus their responses to climate
change (Opdam and Wascher 2004, Kampichler et al. 2012). For example, habitat loss has been
shown to reduce species ability to adjust to changing climatic conditions and thus exacerbate the
detrimental effects of climate change (Travis 2003, Jeltsch et al. 2011). Thus, we can expect the
bioclimatic relationships of species to be altered by the amount of habitat within a landscape.
Specifically, populations of the same species found in regions with abundant suitable cover
should display weaker responses to climatic variability in comparison with populations found in
regions characterized by less available habitat. Here, I explore whether an increasing amount of
habitat affects avian responses to climatic variability and if the potential effect of more extensive
habitat differs for birds inhabiting grassland and forest lands.
Chapter 3 extends the work done in Chapter 2 by focusing on temporal changes of entire
communities. Specifically, it investigates the relationship between temporal changes in avian
assemblages and the interaction of climate change and land-cover fragmentation. To date, the
majority of research regarding the implications of climate change to biodiversity has focused on
responses of individual species (e.g., La Sorte and Thompson 2007). Variation in the individual
species responses is predicted to lead to disruptions of communities and ecosystems, but the
complex nature of ecological interactions makes it difficult to extrapolate from the scales of

3

individuals and populations to the community or ecosystem level (Walther et al. 2002).
Consequently, in order to fully understand consequences of climate change to ecological
communities it is imperative that a direct link between changes in community composition and
changes in climatic conditions is established. Furthermore, the relationship between community
change and climate change will likely be influenced by the landscape matrix surrounding the
community; thus, landscape context needs to be considered when quantifying responses of
ecological communities to climate change. My objective is to investigate the relationship
between temporal change in avian biodiversity and changes in climatic conditions, and to assess
the role of landscape fragmentation as a factor affecting this relationship.
Chapter 4 builds on Chapter 3 by investigating temporal changes in avian assemblages
across multiple spatial scales. Key biodiversity patterns vary with spatial scale of observation
and mechanisms driving these patterns are inherently scale-dependent (e.g., Storch et al. 2004,
Keil et al. 2011). Thus far, studies investigating scaling of biodiversity and the mechanisms
relevant at each spatial scale have mostly focused on one time step. Measures of temporal
changes in biodiversity have been commonly used as indicators of environmental change (e.g.,
Calvarheiro et al. 2013); thus, it is imperative that we understand how these measures depend on
the spatial scale. Furthermore, in the advent of recent global change, a thorough investigation of
mechanisms relevant to temporal biodiversity dynamics across spatial scales is perhaps more
important than ever. My goal is to evaluate spatial scaling patterns of temporal changes in avian
communities and investigate relevant environmental drivers of the community change at each of
the investigated spatial scales.
The collective works in these four chapters improve our understanding of the impacts of
climate change on avian biodiversity and the role that land-cover composition and configuration

4

play in the responses of birds to climate change. I believe that by establishing a link between
climate change impacts on biodiversity and land cover, I provide a much-needed and major step
in global change science.

5

ACCOUNTING FOR THE SPACE-VARYING NATURE OF THE RELATIONSHIPS
BETWEEN TEMPORAL COMMUNITY TURNOVER AND THE ENVIRONMENT

Published in Ecography:
Jarzyna M.A., Finley A.O., Porter W.F., Maurer B.A., Beier C.M. and Zuckerberg B. 2014.
Accounting for the space-varying nature of the relationships between temporal community
turnover and the environment. Ecography, 37, 001-011.
ABSTRACT
Non-spatial regression models are rarely adequate for exploring ecological phenomena,
especially in settings where the processes operate at large spatial scales and when model
covariates do not explain all variation present in the ecological response. Given the complexity
of ecological processes, it is often unrealistic to assume a set of stationary regression coefficients
can capture space-varying and scale-dependent relationships between covariates and an outcome
variable. Spatially-varying coefficients (SVC) models fit within a Bayesian inferential
framework provide a statistically robust method to explore potential space-varying and scaledependent impacts of covariates. My study objective was to assess the utility of SVC models for
capturing non-stationary relationships between temporal community dynamics in avian
assemblages and variation in environmental factors. I also wanted to compare the inference
drawn from SVC models to that obtained from space-varying intercept models and also models
that do not acknowledge any spatial structure beyond what is introduced by the covariates. My
analysis examines the temporal turnover, expressed as a proportion, of avian communities across
New York State, USA. Given the expected outcome is non-Gaussian, I detail a generalized linear
model specification of the proposed model structures. My results show the SVC model
outperformed the spatially-varying intercept and non-spatial models in terms of model fit and
model predictive inference. Further, by fitting these models within a Bayesian inferential
framework, I was able to make inferences about the spatial impact of covariates and other
6

process parameters, as well as obtain full posterior predictive inference about the rate of turnover
at new, unobserved locations. I conclude that SVC models provide a flexible framework for
exploring and accounting for non-stationary mechanisms driving ecological patterns.
1.1 INTRODUCTION
Ecological relationships are often explored using non-spatial regression models, even though this
method is rarely adequate in settings where ecological phenomena exhibit spatial structure that
cannot be explained by model covariates (Cressie et al. 2009, Hoeting 2009). This is especially
true for data collected and modeled across large-scale spatial domains where locations in close
proximity will have similar outcomes, a phenomenon called spatial autocorrelation. Spatial
autocorrelation is prevalent in ecology. For example, community composition at any given
location is usually influenced by the species assemblage structure at surrounding localities and
local environmental factors (Legendre 1993) resulting in species abundances being more similar
than expected by chance in locations closer to each other (Lichstein et al. 2002). Ignoring such
spatial autocorrelation results in dependence among model residuals, which violates assumptions
of most regression models and can lead to erroneous parameter estimates and, ultimately,
incorrect ecological inferences and predictions (Dormann et al. 2007, Hoeting 2009, Finley
2011).
Commonly, spatial autocorrelation in ecological phenomena is accommodated by adding
a spatial random effect to the model mean. Such random effects are often specified to follow a
multivariate normal distribution with a mean of zero and spatially structured covariance matrix.
Such specifications provide local adjustment, with spatial structure, to the model intercept
(Diggle et al. 1998, Banerjee et al. 2004). When the outcome is Gaussian, such spatially-varying
intercept models can partition the residual variance into two components - spatial and non-

7

spatial. For a non-Gaussian outcome (e.g., binomial, Poisson), the residual error term is omitted
from the model equation; thus, spatially-varying intercept models integrate all of the residual
variance into the spatial component. Partitioning of residual uncertainty (for a Gaussian
outcome) or attributing the residual uncertainty to a spatial process (for a non-Gaussian outcome)
can improve inference, increase accuracy and precision of model predictions, and reveal missing
covariates.
While adding a space-varying intercept to the model accounts for spatial dependency and
often improves inference, it does not explicitly deal with spatial non-stationarity and its
associated influence on model covariates. Spatial non-stationarity is present when the strength
and nature of the relationship between a response variable and the predictor covariates vary
across the spatial domain (Fortin and Dale 2009, Miller 2012). Interactions or feedbacks among
unobserved and observed covariates or presence of inherently different mechanisms that impact
the outcome in different parts of the study region can produce non-stationarity (Miller 2012).
Non-stationary relationships are common in ecology (e.g., Bini et al. 2009). For example, Foody
(2004) showed that the relationships between avian species richness and total annual
precipitation, mean annual temperature and terrestrial land cover in sub-Saharan Africa were
strongly non-stationary. Other large-scale studies support these findings (e.g., Grotan et al. 2009,
Martin-Queller et al. 2011, McNew et al. 2013). Thus, given this complexity of ecological
phenomena across large spatial scales, it is often unrealistic to assume a set of stationary
regression coefficients (i.e., coefficients that are constant across space) can capture spacevarying and scale-dependent relationships between covariates and outcome variables (Finley
2011). In other words, a global ecological relationship is often affected by local processes
yielding a heterogeneous pattern, and would be described more accurately by local model

8

parameter values that differ from the global values (Miller 2012). Ignoring spatial nonstationarity can have similar consequences to those resulting from ignoring spatial
autocorrelation. In a strongly non-stationary system, estimates of model parameters and
inferences and predictions resulting from models that do not adequately account for spatial
heterogeneity are likely to be flawed. Despite the importance of the problem and apparent
ubiquity of the space-varying nature of ecological phenomena, few methods provide an
opportunity to evaluate and account for potential spatially varying relationships between
ecological responses and environmental covariates. Bayesian spatially-varying coefficients
(SVC) models are one of the more flexible and robust approaches for accommodating nonstationarity (Gelfand et al. 2003). SVC models use a valid probability model that affords full
posterior inferences for model parameters and subsequent prediction of the outcome and
covariate processes at any new location. While spatially-varying coefficient models have been
developed and tested for data that follow a Gaussian distribution (e.g., Finley 2011), models for
non-Gaussian outcomes, e.g., binomial or Poisson distributions, have not yet been widely applied
to ecological questions, with a few exceptions that focus on statistical methodology development
(i.e., Finley et al. 2009, Finley et al. 2011). Here, I consider a SVC model for evaluating the
relationships between the temporal turnover (i.e., changes in community composition across
time) in avian biodiversity and a set of environmental covariates.
I chose to focus on temporal turnover in community composition for two reasons. First,
patterns and changes in biological diversity are an indication of the underlying mechanisms that
control biodiversity and are one of the fundamental fields of investigation in ecology.
Consequently, understanding patterns of biodiversity and their environmental determinants have
been at the center of recent ecological research (e.g., Buckley and Jetz 2008, Jost 2009, Baselga

9

2010, Kraft et al. 2011). While traditionally community turnover has been used to describe
biodiversity patterns across space rather than time, it can also be viewed in terms of a temporal
change in species composition occurring over a specified time period. Indeed, describing
temporal changes in communities is of increasing importance because of the global
environmental threats such as climate change and land-cover change, both of which have a
strong temporal component. Understanding temporal relationships between biodiversity patterns
and environmental variability is necessary to accurately forecast the consequences of
environmental change and design sound conservation strategies.
Second, the often observed space-varying nature of temporal community turnover and the
kinds of inference we seek encourages the use of new, more flexible, tools such as generalized
linear SVC models. Further, because biodiversity dynamics are often summarized using indices
such as Jaccard, Sorensen, or Diamond-May, there is a need for development of methods
appropriate for such non-Gaussian outcomes.
My objective was to test the utility of the SVC model to quantify the influence of nonstationarity on relationships between temporal community dynamics in avian assemblages and
changes in environmental factors such as climate variability and patterns of landscape
fragmentation. Specifically, I sought to assess the merits of a Bayesian spatially-varying
coefficients approach in comparison with non-spatial and spatially-varying intercept models.
1.2 DATA AND METHODS
1.2.1 Site description
The study area is the State of New York, USA. New York covers 128,401 km2, including 4,240
km2 of inland water (excluding of the boundary-water areas of Long Island Sound, New York
Harbor and lakes Ontario and Erie). Most of the state lies between latitudes 42 and 45°N and
10

between longitudes 73.5 and 79.75°W. Climate of New York is affected by the state’s broad
elevation gradient as well as by the proximity to lakes Erie and Ontario, and the Atlantic Ocean.
Adirondack and Catskill Mountains in the east are among the highest regions of the state, with
elevation varying between 600 and 1,500 m. South-western New York is located in the northern
portions of the Allegheny Plateau and its elevation ranges from approximately 300 to 900 m.
North-western New York as well as Hudson River valley and New York City are low-elevation
regions with elevation ranging from approximately 0 to 200 m.
The state's land-cover is approximately 51.2% forest (deciduous, evergreen, and mixed
forests), 13.4% pasture and hay, 8.2% cultivated crops, 2.9% scrub and shrub, 1.0% grassland
and herbaceous, 8.0% wetland cover, 8.7% developed land, and 0.2% barren land (Homer et al.
2004). Forested areas dominate the Adirondack, Catskill and Allegany regions, while agriculture
is prevalent on historic glacial lake plains south of Lake Ontario.
New York offers a broad gradient of landscape fragmentation. The regions of the
Adirondack and Catskill Mountains are the least fragmented. The Hudson River valley runs
through the eastern part of the state and is more fragmented in terms of land-cover. The
landscapes of western New York are more fragmented and characterized mostly by agricultureforest mosaic. Because of these heterogeneous patterns in landscape fragmentation, elevation
gradients, and diversity of ecosystems, we might expect space-varying impacts of covariates to
explain variation in ecological outcomes of interest.
1.2.2 Breeding Bird Atlas
I used the New York State Breeding Bird Atlas (BBA) as our model dataset to characterize
changes in avian communities through time. BBA is a statewide survey that documented the
distribution of breeding birds in New York. The BBA has been conducted in two time periods,

11

1980-85 (hereafter, BBA1980; Andrle and Carroll 1988) and 2000-05 (hereafter, BBA2000;
McGowan and Corwin 2008). For both BBAs, a grid system was used to define the basic unit for
reporting data. The entire State of New York was divided into approximately 1,300 squares, each
measuring 10 km x 10 km. The BBA reporting unit (a block) was one quarter of this square and
measured 5 km x 5 km and a total of 5,335 blocks covered the entirety of New York State. This
data set represents one of the largest and finest resolution atlases in the world (Gibbons et al.
2007).
A total of 242 species were recorded for BBA1980 and 248 species were recorded for
BBA2000 (Appendix 1.1). Avian breeding was recorded at three levels of certainty of breeding
occurrence based on the behavior of birds observed: possible, probable, and confirmed (Andrle
and Carroll 1988, McGowan and Corwin 2008). Observations were made by skilled birders who
spent at least 8 hours in each block, visited all cover types in each block, and included at least
one nighttime visit to document nocturnal species. Observer effort was recorded for each BBA
and reported as a number of person hours (McGowan and Zuckerberg 2008). A block was
considered sufficiently surveyed when at least 76 species were documented (with exceptions of
blocks that might be expected to have fewer species). The BBA represents a presence/absence
dataset, although absence indicates that species could not be found given search criteria
(McGowan and Corwin 2008).
1.2.3 Temporal turnover
While many different indices for quantifying community turnover has been developed over the
years (Gaston et al. 2004, Magurran 2004, Tuomisto 2010a, Tuomisto 2010b), we used
Diamond-May (DM) index (Gaston et al. 2004). DM index is calculated as a proportion and
designates turnover as high when the proportion of species shared between two sites (or two time

12

steps) is low (Gaston et al. 2004). From a statistical perspective, DM is a binomial outcome,
where the number of species lost and gained across space or time is the number of successes,
while the total number of species is the number of trials. DM index has traditionally been used to
quantify spatial turnover, but it can be easily adapted to reflect temporal turnover as follows:
𝐷𝑀 =

𝐸+𝐶

𝐸+𝐶+𝑃

,

Eqn. 1.1

where E is the number of species that went extinct in the block between BBA1980 and

BBA2000 (i.e., extinction), C is the number of species that colonized the block between
BBA1980 and BBA2000 (i.e., colonization), and P is number of species in the same block
common to both BBAs (i.e., persistence). The values of DM are bounded by 0 and 1; values
approaching 0 indicate low temporal turnover in a block between BBA1980 and BBA2000,
while values approaching 1 indicate high temporal turnover in a block between BBA1980 and
BBA2000.
1.2.4 Model covariates
Environmental factors that I deemed especially important in shaping changes in community
composition were trends in climatic variables and landscape fragmentation. I also considered
survey effort (number of survey hours) as a variable potentially affecting the observed temporal
community turnover.
1.2.4.1 Climatic trends
The climate data was derived from the PRISM (Parameter-elevation Regressions on Independent
Slopes Model) climate mapping system (Daly and Gibson 2002). PRISM consists of interpolated
monthly maximum and minimum temperatures and precipitation at a 2.5-arcmin resolution from
1891-2010 for the entire contiguous United States.

13

I calculated the magnitude of the 25-year (1980-2005) trend in average monthly
maximum and minimum temperatures and in total monthly precipitation using Ordinary Least
Squares regression. The slope of the OLS regression indicates the magnitude of the trend and
reflects the amount of change in climatic variables that occurred between 1980 and 2005. I then
interpolated the trend magnitudes for each BBA block and averaged the monthly values to reflect
breeding season (May through September) trend magnitude in maximum and minimum
temperatures (TMAX and TMIN, respectively) and in total precipitation (PRECIP). The trend
magnitudes of TMAX and TMIN were expressed in °C/25 years, while the units of the trend
magnitudes of PRECIP were in mm/25 years.
1.2.4.2 Habitat fragmentation
Even though amount of suitable habitat is often cited as one of the most important factors driving
patterns of species occurrence, it is difficult to quantify suitable habitat while dealing with a
group of species with vastly different habitat requirements. Habitat fragmentation has been
shown to also be a significant driver of songbird community dynamics (e.g., Kennedy et al.
2011). Hence, I chose landscape fragmentation as the primary land-cover factor associated with
community turnover.
Habitat fragmentation variables were derived using the 30x30 m National Land Cover
Data (NLCD, Homer et al. 2004) product. Because there is no land-cover data readily available
for the time period of BBA1980, I used a space-for-time substitution to assess whether temporal
turnover is related to landscape fragmentation (Pickett 1989, Zuckerberg and Porter 2010).
Space-for-time substitution assumes that spatial variation is equivalent to a temporal variation,
and relationships between variables (e.g., species occurrence and land-cover) derived for one
time step across a large spatial extent will be equivalent to that derived for two time steps. Thus,

14

I quantified landscape fragmentation for one time step only and looked for the relationship
between the temporal turnover and the observed spatial variation in landscape fragmentation.
I used the 2001 NLCD, because it coincided well with the time of BBA2000. The 2001
NLCD consists of 16 land cover classes. Prior to landscape analysis, I consolidated open space
developed, low intensity developed, medium intensity developed, high intensity developed into
one class of developed land; cultivated crops and pasture/hay into one class of agriculture; and
deciduous forest, evergreen forest, mixed forest, and forested wetland into one class of forest.
We decided to consolidate these land-cover types to improve accuracy of classification and to
simplify the environment for purposes of our evaluation. Upon consolidation, I examined nine
land cover classes: water/ice, developed, barren land, forest, scrub/shrub, grassland/herbaceous,
agriculture, and wetlands.
Landscape fragmentation was quantified using FRAGSTATS 4.1 (McGarigal et al. 2002)
and the Geospatial Modelling Environment (GME, http://www.spatialecology.com/gme/). I
recognize that there are multiple ways to measure fragmentation and that careful characterization
of landscape fragmentation for a diverse suite of species often requires multiple metrics.
However, my purpose was not comprehensive ecological analyses, but rather evaluation of the
space-varying qualities of fragmentation. Therefore, I chose a landscape-scale variable to capture
broad-scale variation in habitat fragmentation and best represent requirements of a diverse suite
of species with varying habitat requirements. Specifically, I chose Edge Density (ED) as a
measure of landscape fragmentation because an increase in habitat edge is a primary outcome of
habitat fragmentation (Hargis et al. 1998). ED was also reported as an effective tool for
evaluating landscape fragmentation and performed better than other popular landscape
fragmentation indices (Hargis et al. 1998). However, in situations when landscape consists

15

entirely of one cover type, the ED would be 0 regardless of the type of land-cover present.
Therefore, in order to differentiate between landscapes that consisted entirely of natural land
cover (e.g., forest) and those mostly developed (e.g., urban areas), I also included the proportion
of developed land (DEVEL) metric in our models.
1.2.4.3 Survey effort
Increasing survey effort often results in a higher number of recorded species (Tobler et al. 2008).
Different survey effort between BBA1980 and BBA2000 could decrease the value of the DM
indicating high temporal turnover, which could be a result of different number of recorded
species rather than actual changes in species identities throughout time. Despite the fact that
concerted attention was invested to minimize the effects of survey effort in the BBA surveys, I
wanted to recognize the need to control for survey effort bias. To account for potential survey
effort bias, I calculated the absolute difference in the number of person hours between BBA1980
and BBA2000 (EFF=|EFF1980-EFF2000|) for each BBA block and included it as a covariate in
our models.
1.2.5 Statistical analysis
I evaluated three statistical models, each of which included main effects of all the covariates
(TMAX, TMIN, PRECIP, ED, DEVEL, and EFF). I standardized all the covariates to ease the
comparison and interpretation of the coefficient estimates. Since the temporal turnover (DM) is a
proportion, I used a binomial regression model. I started with a non-spatial regression
specification (hereafter, non-spatial model) and then increased the level of model complexity by
adding spatially structured random effects first to the intercept to form a spatially-varying
intercept (SVI) model, then to the regression parameters associated with the covariates to form
the SVC model.
16

More formally, for a given BBA grid cell the regression outcome variable y(s) is the
numerator in Eqn. 1, where s represents the grid cell's spatial coordinates. Given the denominator
in Eqn. 1, denoted N(s), we assume y(s) follows a binomial distribution. For the given grid cell
π(y(s)│η(s)) ~ Binomial(N(s), p(η(s))), where p(η(s)) is the probability of success at location s
and η(s) is model’s regression equation, which for the SVC model is equal to x(s)β + x(s)w(s)

with x(s) representing the vector comprising an intercept and location specific covariate values.
The parameters to be estimated are the global regression coefficients, β = (β0, βTMAX, βTMIN,
βPRECIP, βED, βDEVEL, βEFF)T, and associated spatially-varying adjustments, w(s) = (w0(s), wTMAX(s),
wTMIN(s), wPRECIP(s), wED(s), wDEVEL(s), wEFF(s))T. A logit link function η(s) = log[p(s)/{1-p(s)}]
was used for this model.
I assume each element in the spatially-varying coefficients vector, w(s), arises from a
spatial Gaussian Process (GP) (see, e.g., Banerjee et al. 2004 or Cressie and Wikle 2011 for more
details). Specifically, the j-th spatially-varying coefficient at location s is wj(s) ~ GP(0, Cj(s, s’)),
where s and s’ are any two locations within the study area, the spatial covariance Cj(s, s’) =
σ2jρj(s, s’; ϕj) with variance parameter σ2j, correlation function ρj(·; ϕj), and spatial decay
parameter ϕj. The spatial decay parameter ϕj controls the decay in spatial correlation; lower
values of the parameter indicate shorter range in spatial correlation. The exponential spatial
correlation was assumed for ρj(·).
Prior distributions on the remaining parameters complete the hierarchical model, see,
e.g., Gelman et al. 2004 for detail on Bayesian model specification. Here, I specified
uninformative prior distributions. The regression coefficients, β j’s followed a Normal
distribution N(0, 100), while the variance components σ2j’s were assigned inverse-Gamma IG(2,
1) priors, where 2 and 1 are shape and scale hyper-parameters, respectively. This specification of

17

the IG provides a prior distribution with infinite variance that is centered on 1. The spatial decay
parameters, ϕj 's, followed an informative Uniform prior, with support that ranged from 1 km to
the maximum distance between any two grid cells. Exploratory analysis using different hyperparameters suggested prior specification had little influence on parameter and predictive
inference. This is not surprising given the large size of the data set and choice of prior
distributions. My sensitivity analysis did suggest that given an IG shape hyper-parameter of 2,
very large scale parameter values, e.g., greater than 10, resulted in poor MCMC chain
convergence by forcing an unreasonable amount variability in the GP. For the SVC model, I ran
three MCMC chains for 60,000 iterations each. I ran 10,000 iterations for the non-spatial and
SVI models, as those required less time to converge. Convergence was diagnosed using the
Gelman-Rubin diagnostic (Gelman and Rubin, 1992).
Given the Bayesian framework, candidate models fit to the observed data were assessed
using the Deviance Information Criterion (DIC; Spiegelhalter et al. 2002). I computed the
expected posterior deviance as D(Ω) = EΩ|Y{-2logL(Data| Ω)}, where Ω is the set of parameters
estimated for each model and L(Data| Ω) is the first stage Gaussian likelihood from the
������� – D(𝛀
� ), where
respective models. I computed the effective number of parameters as pD = 𝐷(𝛀)
� is the posterior mean of the model parameters. The DIC was then computed as D(Ω)+ pD.
𝛀

Lower values of DIC indicate improved fit. Model accuracy was compared using Root Mean
Square Error (RMSE) that was calculated as the square root of the mean squared deviations
between fitted and observed outcomes. Lower values of RMSE indicate improved accuracy.
Candidate model performance for out-of-sample prediction was assessed using a
randomly selected 10% holdout dataset. While there is no convention on how much of the
original data should be held out for model validation, we believe that 10% of the data is

18

sufficient as it amounts to approximately 500 observations. Also, while several other approaches
to verify model’s predictive ability are available, holding-out a set of randomly selected
observations is a common model validation technique (e.g., Pearson et al. 2002, Finley 2011).
An alternative approach to verifying model’s predictive ability would be bootstrapping, where
the original dataset is sampled randomly with replacement (e.g., Austin and Tu 2004).
Bootstrapping requires building and running multiple models and in each case predictive
performance is assessed against the corresponding test data. While bootstrapping is a robust
approach, it is computationally costly, particularly in situations when models are complex. Given
that fitting a single SVC model is very computationally intensive in itself, repeating this
procedure multiple times for bootstrapping would be prohibitive and thus was not feasible in our
case. Hold out set predictions were compared to the observed values using Root Mean Square
Error (RMSEpred) where lower values indicate improved performance.
Prior to statistical analysis, I removed all blocks that did not have continuous land cover
coverage; those were mainly blocks at the periphery of the study area, which extended farther
than the NLCD layer. Also, blocks with more than 50% open water coverage and those that did
not have survey effort data were removed. Ultimately, I used a total of n = 4,271 blocks to fit the
models and nho = 473 holdout blocks to test the models’ predictive ability. The non-spatial and
SVI models were run in R 2.15.1 statistical package (R Development Core Team 2013,
http://www.r-project.org/) using package spBayes (Finley and Banerjee 2013). In close
collaboration with Dr. Andrew Finley, I wrote and compiled the SVC model in C++
programming language and R 2.15.1 (R Development Core Team 2013, http://www.rproject.org/). Summaries of parameter estimates were generated using the R CODA package.

19

1.3 RESULTS
1.3.1 Temporal turnover
The mean grid cell value of the DM index was 0.39, which indicates that on average
approximately 39% of species had either colonized or become extinct and 61% of the species
were common to both BBAs (Fig. 1.1A). Northern regions of the state (i.e., Adirondack
Mountains) experienced the highest temporal community turnover, with values reaching 1.0 (i.e.,
complete turnover in species assemblage) in several locations (Fig. 1.1A). Western (agricultural)
and Eastern (urbanized) parts of the state also showed slightly higher than average temporal
turnover, while the central regions of the state generally underwent the lowest temporal turnover
(Fig. 1.1A). The lowest values of temporal turnover was approximately 0.11.

20

Figure 1.1 Observed temporal turnover in communities of breeding birds in New York
expressed using Diamond-May index (A) and the fitted temporal turnover values resulting from
the non-spatial model (B), spatially-varying intercept (SVI) model (C), and spatially-varying
coefficients (SVC) model (D).

21

1.3.2 Covariates
The north-eastern and south-western parts of the state experienced the largest warming trend,
while the southern and central New York generally underwent decreases in average maximum
temperatures (Appendix 1.2). Minimum temperatures increased across most of the state, though
some small regions have experienced a cooling trend (Appendix 1.2). Northern (i.e., region of
the Adirondack Mountains), western and southern New York experienced drier conditions, while
the rest of the state experienced wetter conditions (Appendix 1.2).
Regions of Adirondack and Catskill Mountains as well as some locations in the south and
south-western parts of the state had the lowest edge density (i.e., lowest habitat fragmentation;
Appendix 1.2). The percentage of developed land was low across most of New York, with
exception of large urban centers such as New York City, Albany, Syracuse, Rochester, and
Buffalo. Survey effort was relatively constant across the study region (Appendix 1.2).
1.3.3 Statistical analysis
1.3.3.1 Model fit and predictive ability
As suggested by the DIC, RMSE, and RMSEpred values provided in Table 1.1, adding spatial
random effects to the posited model improves fit to the observed data, model accuracy, and
predictive performance for new, unobserved, locations. The improvements to fit and model
accuracy suggest there was substantial spatial dependence among the non-spatial models
residuals - a generalized linear model assumption violation - and clear indication that a more
complex spatial model was warranted. The pD in Table 1.1 is the effective number of parameters
and is used in the DIC calculation to penalize more complex models. For example, the nonspatial model has only seven regression coefficients, hence, pD is approximately seven. The
increasing number of spatial random effects added to SVI and SVC, respectively, is reflected by
22

larger penalties. Despite the larger number of effective parameters, the SVC model has lower
DIC than the non-spatial and SVI models.
Table 1.1 Comparison of the non-spatial, spatially-varying intercept (SVI), and spatially-varying
coefficient (SVC) models for four migratory groupings using Deviance Information Criterion
(DIC), estimated number of parameters pD, and Relative Mean Square Error (RMSE). Predictive
performance of the models was evaluated by calculating RMSEpred of the observed values from
the hold-out set and the predicted values resulting from all three models.
Model
SVC
SVI
Non-spatial

ΔDIC
0.00
1,063.84
2,871.31

pD
458.90
106.30
7.20

RMSE
0.076
0.082
0.090

RMSEpred
0.072
0.079
0.088

All candidate models provide a close visual approximation to the observed data (see Fig.
1.1 for the candidate model fitted value surfaces). For the SVI model (Fig. 1.1C), the addition of
the spatial random effect to the model intercept provides a smooth process that captures spatial
discrepancies between observed values and those estimated with the additive linear trend
involving the covariates. The SVC model fitted surface (Fig. 1.1D) provides the closest
approximation to the observed data because the regression coefficients processes are able to
estimate the covariates’ local and regional impact on temporal turnover. For example, the
variation in the temporal turnover in the north-east New York is better represented on the SVC
model fitted surface (Fig. 1D) than on surfaces resulting from the other two models (Fig. 1.1B
and 1.1C).
In terms of model predictive performance, all models’ predicted surfaces approximated
the holdout data relatively well (Fig. 1.2). However, by not relying only on the smooth process
on the intercept, the SVC model is able to accommodate more local-scale variation than the SVI
model (Fig. 1.2D). For example, several locations of high temporal turnover in the western part
of the study area were picked up by the SVC predictive model, but not by either the non-spatial
or the SVI predictive models (Fig. 1.2D).
23

Figure 1.2 Observed temporal turnover of the holdout data set (A) compared with the surfaces of
the predicted values resulting from the non-spatial model (B), spatially-varying intercept (SVI)
model (C), and spatially-varying coefficients (SVC) model (D).

24

1.3.3.2 Coefficient estimates
All regression coefficients in the non-spatial model were statistically significant (Table 1.2), i.e.,
the 95% credible intervals did not include zero. For spatially-varying intercept models, all of the
regression coefficients remained significant (Table 1.2), but their magnitude and credible
intervals changed in comparison with the non-spatial model (the magnitude of the coefficient
estimates was assessed by simply comparing the direction and absolute value of the mean
estimates resulting from both models). The magnitude of β0, βTMAX, βTMIN, βPRECIP, and βED
decreased after the residual spatial dependence was accommodated, while βDEVEL and βEFF
increased in magnitude.
Table 1.2 Coefficient estimates resulting from the non-spatial and spatially-varying intercept
(SVI) models for all model covariates: intercept (β0 ), magnitude of the 25-year (1980-2005)
trend in average maximum temperature of the breeding season (βTMAX), magnitude of the 25-year
(1980-2005) trend in average minimum temperature of the breeding season (βTMIN), magnitude of
the 25-year (1980-2005) trend in average total precipitation of the breeding season (βPRECIP),
percent developed land (βDEVEL), edge density (βED), and effort (βEFF). 50% indicates the mean of
the posterior distribution, while 2.5% and 97.5% are the lower and upper quantiles of the
posterior distribution.
Model

Non-Spatial

SVI

Coefficient
β0
βTMAX
βTMIN
βPRECIP
βDEVEL
βED
βEFF
β0
βTMAX
βTMIN
βPRECIP
βDEVEL
βED
βEFF

Coefficient Estimates
50%
2.5%
-0.492
-0.499
0.063
0.055
0.041
0.034
-0.052
-0.059
0.017
0.010
-0.092
-0.099
-0.010
-0.016
-0.375
-0.435
0.018
0.007
0.026
0.014
-0.015
-0.031
0.032
0.022
-0.047
-0.057
-0.014
-0.020

25

97.5%
-0.486
0.070
0.048
-0.045
0.024
-0.086
-0.004
-0.275
0.029
0.038
-0.001
0.042
-0.036
-0.008

For a SVC model, it is typically not instructive to focus on a covariate's regression
coefficient mean, i.e., β, but rather assess the variability of the coefficient over the study area.
That is, we draw inference by looking to maps of β + w(s) (see Fig. 1.3 for β + w(s) associated
with each covariate). All SVC model coefficient estimates were spatially heterogeneous both in
terms of the magnitude of the coefficient estimate as well as the direction of the relationship
between the covariate and the outcome variable (Fig. 1.3). Here too, grid cells where the 90%
credible intervals of the given β + w(s) include zero are identified with a black point. Or
conversely, those regions with no black points differ significantly from zero in the direction
indicated by the grid cell color. Some caution needs to be taken when considering the hypothesis
tests implied in Fig. 1.3. In a testing situation, where we encounter a very large number of
hypotheses, controlling Type-1 errors would be required, hence raising the question of multiple
comparisons. However, we do not pursue such a formulation here and simply look to Fig. 1.3 to
identify regions of substantial deviations of the spatial coefficients from the global mean (see
also Fig. S2 in the Supplementary Material for maps of deviations of the spatially-varying
coefficient estimates from the global mean).

26

Figure 1.3 Coefficient estimates resulting from the spatially-varying coefficient (SVC) models of (A) the intercept, (B) the magnitude
of the 25-year (1980-2005) trend in average maximum temperature of the breeding season (βTMAX), (C) the magnitude of the 25-year
(1980-2005) trend in average minimum temperature of the breeding season (βTMIN), (D) the magnitude of the 25-year (1980-2005)
trend rates in average total precipitation of the breeding season (βPRECIP), (E) percent developed land (βDEVEL), (F) edge density (βED),
and (G) effort (βEFF). Black point symbols identify grid cells with posterior 90% credible intervals that include zero. The surface
colors identify the direction of the regression coefficients' sign.

27

1.4 DISCUSSION
We found that the SVC model outperformed both the non-spatial and SVI models in terms of
model fit, despite the large number of parameters it had to estimate, and predictive performance.
We feel confident that the increase in model fit is not a result of over parameterization, although,
like any highly flexible random effect model, the SVC models can occasionally over fit the
observed data. However, compared to models with unstructured random effects, the flexible
process imposed by the spatial Gaussian Process reduces the possibility of over fitting.
Additionally, the penalizing qualities of the DIC metric account for the potential favoring the
parameter rich SVC model. Furthermore, over fitting generally results in reduced predictive
performance compared to simpler models; however, we did not find lower predictive
performance of the SVC model. Anecdotally, while working with other data sets, we have
noticed that the propensity to over fit the observed data is greater when the number of observed
locations is small (e.g., n < 100), sparsely sampled, and when the range of spatial dependence is
short. In these settings, the interpolating qualities of the Gaussian Processes reduce the predictive
power of the SVC model. When the data set is large, however, observations densely cover the
domain, and the parameters of the underlying spatial processes associated with the regression
coefficients are well estimated, then the risk of over fitting is reduced and prediction does not
suffer.
The improvement in model fit suggests the relationship between temporal community
turnover and environmental covariates is more spatially-varying in nature than it is stationary.
Indeed, by allowing the regression coefficients to vary spatially over the domain and
accommodate local impact of the covariates, the SVC model was able to explain more localscale variation in temporal community turnover than the other two models. The assumption that

28

regression coefficients are stationary results in poor fit and misleading inference about the
impact of the covariates on temporal turnover. The assumption of scalar regression coefficients is
pervasive in the literature. In the SVC model, any component of model spatial residual pattern
due to non-stationarity of the impact of covariates is apportioned to the coefficients. Thus, the
SVC model, despite being an intrinsically correlative approach, provides the richer opportunity
for ecological interpretation.
We commonly see that once residual spatial dependence is accommodated, e.g., via a
spatial random effect, regression coefficients differ from those resulting from a non-spatial
model. In the non-spatial model, all covariates were significant predictors of temporal
community turnover. Once spatial dependence was taken into account in the SVI model, the
covariates remained significant but the magnitude of their coefficient estimates changed. For
example, the association of the TMAX, TMIN, PRECIP, and DEVEL with the temporal
community turnover was weaker for the SVI model than for the non-spatial model, while that of
ED and EFF became stronger. Such results suggest the non-spatial model violated the
assumption of independent and identically distributed residuals and imply that ecological
inferences drawn from these models will be different. Other research corroborates our findings.
For example, Record et al. (2013) found that climatic variables were significant predictors of
presence-absence of two tree species in non-spatial models, but in the spatial models the
magnitude and sign of some parameter estimates changed. Foody (2004) found the magnitude of
relationships between species richness and the explanatory environmental variables changed
when spatial dependency and non-stationarity were accounted for, even though statistically
significant relationships were established using a conventional global regression analysis.

29

When the space-varying nature of the regression coefficients is accommodated, different
regions within the study area might exhibit conflicting relationships between the outcome and
covariate, causing the study-area-wide mean of the associated regression coefficient to be similar
to those seen in the SVI models. In our case, spatial heterogeneity of the relationships between
covariates and temporal community turnover became apparent after accounting for nonstationarity. Some regions were characterized by negative relationships, other regions displayed
positive relationships, yet another showed no significant relationships at all. For example, a
positive relationship between temporal community turnover and the magnitude of the trend in
minimum temperatures was found in the eastern New York, while this relationship was negative
in the central part of the state. Such spatial pattern might indicate that communities in eastern
New York are responding stronger or faster to increasing temperatures resulting from climate
change than those found in central New York. Similarly, we found that the strongest negative
influences of landscape fragmentation were in eastern part of the state, i.e., in regions dominated
by the Adirondack Mountains and the Hudson River valley. Perhaps in those historically forested
and contiguous habitats, small levels of landscape fragmentation have a much more detrimental
effect on community composition than elsewhere within the state. Furthermore, given that the
strongest influences of these two covariates were detected in similar geographic regions, it is also
plausible that these two ecological factors interact in driving temporal community turnover. The
remaining covariates showed similar heterogeneous patterns. These heterogeneous surfaces
illustrate the complexity of impact and interplay among the covariates in explaining variability in
temporal community turnover.
Other research provides evidence for non-stationarity of ecological relationships (e.g.,
Bini et al. 2009, Martin-Queller et al. 2011, McNew et al. 2013) and thus corroborates findings

30

from our study. For example, Grotan et al. (2009) found the proportion of variance in local
recruitment of great tit explained by spring temperature differed nearly 10-fold among four
Dutch populations. Sæther et al. (2008) showed that population dynamics of eight species of
prairie ducks showed pronounced latitudinal gradient; depending on the location, influence of
increased spring temperature or winter precipitation on population dynamics was either positive
or negative. Anders and Post (2006) found that NAO affected populations of the migrant yellowbilled cuckoo (Coccyzus americanus) only in southern and eastern parts of the US, while El Nino
Southern Oscillation (ENSO) affected populations in northern and central US. Martin-Queller et
al. (2011) found relationships between species richness of woody plants and environmental and
biotic factors to also be spatially heterogeneous. Given this apparent ubiquity of space-varying
relationships in nature, we suggest that researchers studying large-scale ecological phenomena
should be especially cognizant of the need to evaluate and account for potential non-stationarity.
The three models we tested yielded different results despite the fact that the same
covariates were included in all of them. By ignoring spatial dependency in model residuals, the
non-spatial model produced relationships between the covariates and temporal community
turnover. In the SVI model, the relationships between covariate and outcome variable generally
weakened as a result of accommodating the spatial autocorrelation in model covariates via
spatial random effect. By allowing the coefficient estimates to vary spatially, the SVC model
yielded heterogeneous patterns and often conflicting relationships between the outcome and
covariate. The heterogeneity of the model covariates themselves (e.g., the magnitude of the
climatic trends) might be partly responsible for driving such space-varying relationships.
Applying the SVC model enabled us to detect and account for these spatial differences.

31

It must be noted that other approaches exist to evaluate and account for potential nonstationarity of ecological phenomena. One of the more common options is a frequentist method
called geographically-weighted regression (GWR, Fotheringham et al. 2002), which uses spatial
weights to estimate spatially adaptive coefficients. GWR has been the leading statistical method
to account for spatial non-stationarity (e.g., Ma et al. 2012a, Ma et al. 2012b, Miller 2012), in no
small part due to its availability in several popular Geographic Information Systems, e.g., ESRI
products, and statistical computing environments, e.g., spgwr package in R statistical program
(Bivand and Yu 2013). However, recently it has been shown that GWR is not robust to
collinearity among the covariates and the presence of complex spatial correlation structures
(Wheeler and Waller 2009, Finley 2011) and it tends to provide less accurate and more
correlated regression coefficient estimates than those resulting from SVC models (Wheeler and
Cadler 2007). Furthermore, Ma et al. (2012a) used GWR to model species richness derived from
the same dataset as ours and concluded that GWR performed poorly the larger the spatial scale of
investigation. Another central shortcoming of GWR, from an inferential standpoint, is the lack of
a legitimate probability model in the sense that the joint distribution linking the parameters and
the data is not a valid probability distribution. This is problematic for inference because the
standard errors computed from such models may not be justifiable. Asymptotic arguments may
be supplied but are complicated in spatial contexts because of the divergent paradigms of infill
and expanding domain asymptotics. Thus, GWR can be a useful tool for exploratory data
analysis, but generally should not be used in settings where one seeks inference about the
importance of model parameters or in prediction.
Given the growing wealth of space and time indexed ecological data (Kelling et al. 2009,
Reichman et al. 2011, Schimel 2011) deeper insight into ecological processes is contingent upon

32

our ability to specify valid and tractable models able to accommodate often complex spatiotemporal relationships between outcomes and posited environmental drivers. Because the
propensity for non-stationary relationships increases over larger spatial domains and time
periods, need for more sophisticated models also increase as we begin looking at macro-system
questions. This study represents an important case assessing the utility of SVC models for
capturing non-stationary relationships between non-Gaussian biodiversity dynamics and
environmental changes. We suggest that SVC models are especially appropriate for studying
macro-scale (i.e., regional or global) systems where we can expect the ecological processes and
environmental drivers to be strongly heterogeneous across the spatial domain. As such, SVC
models provide a unique opportunity to explore pressing environmental questions about impact
of global environmental change on biodiversity.
1.5 ACKNOWLEDGMENTS
I would like to thank the thousands of volunteers who participated in both New York State
Breeding Bird Atlases. I also thank Kimberley Corwin and Kevin McGowan for supplying atlas
databases and other information. The manuscript benefited from discussions with the members
of the Quantitative Wildlife Laboratory at Michigan State University. I thank James H. Stagge
for technical assistance with hyper-performance computing. I also thank two anonymous
reviewers for their comments on the earlier versions of the chapter. This study received financial
support from NASA and Boone and Crockett Club.

33

SYNERGISTIC EFFECTS OF CLIMATE AND LAND COVER: ARE GRASSLAND
BIRDS MORE VULNERABLE TO CLIMATE CHANGE?

Submitted for publication in Diversity and Distributions:
Jarzyna M.A, Porter W.F, Finley A.O. and Zuckerberg B. In review. Synergistic effects of
climate and land cover: Are grassland birds more vulnerable to climate change? Diversity
and Distributions.
ABSTRACT
Quantity and quality of available habitat will affect the way species respond to climatic
variability. My goal was to investigate bioclimatic relationships of grassland and forest breeding
birds along varying gradients of grassland and forest habitat. Specifically, I sought to investigate
whether: (i) bioclimatic relationships of birds breeding in grasslands differ from those of forest
species, (ii) an increasing amount of habitat affects avian responses to climatic variability, and
(iii) the potential effect of more extensive habitat differs for birds inhabiting grassland and forest
lands. I used the New York State Breeding Bird Atlas as a model dataset and I chose 11
grassland and 12 forest birds. To evaluate bioclimatic relationships along varying amounts of
habitat, I created spatially-varying coefficients models (SVC) by adding spatially-structured
random effects to the regression parameters associated with each of the model covariates. The
resulting spatially-varying coefficient estimates of all climatic covariates were then regressed
against the habitat amount for each individual species. I found that grassland and forest birds
exhibited divergent bioclimatic relationships. On average, grassland birds were more likely to
occur in regions characterized by higher temperatures and lower precipitation when compared to
forest breeding birds. I also found that both groups exhibited heterogeneous bioclimatic
relationships across a range of habitat amount, though the influence of land cover differed for
grassland and forest birds. I conclude that increasing amount of forest cover likely will buffer the
negative consequences of climate change, whereas extensive open grasslands will likely
34

exacerbate the impacts of climatic variability on bird populations. Thus, grassland birds might be
particularly at risk of a changing climate due to the diminished buffering capacity of grassland
systems.
2.1 INTRODUCTION
Global temperature increase resulting from anthropogenic climate change will have far-reaching
consequences for biodiversity, and there is growing evidence that most taxonomic groups are
already responding to recent climate warming (Walther et al. 2002, Parmesan and Yohe 2003,
Jarzyna et al. 2013). However, despite the evidence that climate change is impacting
biodiversity, there is still a need to better understand which species and ecosystems are most
vulnerable, and under what conditions (Lawler 2009).
Deriving bioclimatic relationships is a critical component of determining the potential
impacts of climate change on species and communities (e.g., Guisan and Thuiller 2005, Araújo
and Peterson 2012) and remains the most widely used tool for forecasting changes in species
distributions (e.g., Miller 2010). However, using bioclimatic relationships as a sole predictor of
species’ distributional changes is limited (Pearson et al. 2004, Huntley et al. 2010) because
factors other than climate will affect the way species respond to climate change (Thomas et al.
2004). For example, research suggests the quantity and quality of available habitat might be
important in shaping bioclimatic relationships of species and thus their responses to climate
change (Opdam and Wascher 2004, Kampichler et al. 2012). In simulation studies, habitat loss
has been shown to reduce species ability to adjust to changing climatic conditions and thus
exacerbate the detrimental effects of climate change (Travis 2003), while an increase in the
amount of habitat has been shown to slow down species extinction rates (Jeltsch et al. 2011).
Empirical evidence, however, is lacking because studies that integrate both kinds of impact (land

35

cover and climate) are rare (but see Kampichler et al. 2012). Moreover, no studies that explicitly
evaluate species bioclimatic relationships across a broad range of habitat amount, to the best of
our knowledge, have yet been conducted. Gaining a better understanding of the synergistic
effects of climate and habitat will provide a deeper understanding of which species are more or
less vulnerable to climate change.
Given that decreasing habitat is likely to alter responses of species to climatic variability
and thus potentially exacerbate the negative consequences of climate change (Travis 2003,
Jeltsch et al. 2011), we expect the bioclimatic relationships of species to be altered by the amount
of habitat within a landscape. Populations of the same species found in regions with abundant
suitable cover should display weaker responses to climatic variability in comparison with
populations found in regions characterized by less available habitat. Vulnerability to climate
change can then be assessed based on the observed strength of bioclimatic associations. Species
exhibiting weaker responses to climatic variability and those whose habitats provide modulating
capacity against changing climatic conditions (i.e., where increasing amount of habitat weakens
species’ bioclimatic relationships) would be considered less vulnerable to changing climatic
conditions.
The modifying characteristics of habitat amount on bioclimatic relationships are likely to
differ depending on the type of land cover. Recent research suggests that forest canopies
moderate the impact of climate warming on understory plants by creating cooler microclimatic
conditions during the summer months (Bonan 2008, De Frenne et al. 2013). On the other hand,
open habitats, such as old fields, agricultural, and pasture land (hereafter called grasslands),
generally create warmer, drier, and more variable microclimates than forests (Geiger et al. 2009)
and are therefore less likely to provide a similar moderating effect (Villegas et al. 2010). Given

36

these differences in climatic regimes of both types of cover, and the fact that thermal and habitat
niches of species are often correlated (Barnagaud et al. 2012), we might expect that species
inhabiting grasslands and forests will display divergent bioclimatic relationships. We can also
expect that the presence, absence, and magnitude of the potential mitigating effect of land cover
will be different for these two groups of species.
The comparison of grassland and forest species is timely. The velocity of climate change
has been shown to be greater along latitudinal than elevational gradients (Dobrowski et al. 2013).
Forests are generally found in regions of higher altitudinal gradient than grasslands, which in
turn are generally found in low laying and flat areas. Thus, given a similar latitudinal extent,
grasslands are likely to experience more rapid changes under future climatic conditions than
forested regions. However, how exactly grassland species will respond to this higher rate of
climatic variability is largely unknown. Quantifying the range of bioclimatic relationships across
a gradient of available habitat would be the first step in determining potential impacts of climate
change on these two groups of species.
Here, I used data from the New York State Breeding Bird Atlas to investigate bioclimatic
relationships of grassland and forest birds across a range of habitat. I focus on relationships with
temperature and precipitation, as the main climatic factors driving distributions and breeding
success in bird populations (e.g., Moss et al. 2001). Specifically, I seek to evaluate three
interrelated questions: (i) do bioclimatic relationships of birds breeding in grasslands differ from
those of forest species? (ii) do increasing amounts of habitat modify responses of grassland- and
forest-breeding birds to climatic variability? and (iii) is the effect of increasing habitat amount
different for grassland and forest birds? First, I predict that grassland and forest birds will exhibit
divergent bioclimatic relationships. Specifically, forest birds will have stronger negative

37

relationships with temperatures due to their preference for cooler habitats than grassland birds,
while grassland birds will be found in warmer and drier regions. Second, I predict that an
increasing amount of forest cover will abate the negative responses of forest birds to climatic
variability (Cox et al. 2013). Lastly, I predict no change in the magnitude of bioclimatic
relationships of grassland birds as the amount of their habitat increases. If my predictions
regarding the mitigating qualities of each land cover hold true, we can expect grassland birds to
be particularly vulnerable to changing climate due to the diminished modulating capacity
provided by their habitat.
2.2 METHODS
2.2.1 Breeding Bird Atlas and chosen species
The New York State Breeding Bird Atlas (NYSBBA) is a statewide survey of the distribution of
breeding birds in New York that has been conducted to date in two time periods, 1980-85
(Andrle and Carroll 1988) and 2000-05 (McGowan and Corwin 2008). For both atlases, a grid
system was used to define the basic unit for reporting data; each unit (a block) measured 5 km x
5 km and a total of 5,335 blocks covered the entirety of New York State. The NYSBBA
represents one of the largest and finest resolution atlases in the world (Gibbons et al. 2007) and,
in comparison with other datasets of avian distributions, provides a full state coverage in a form
of a spatial grid that allows for thorough and systematic evaluation of species’ distributional
patterns. The NYSBBA is a detection/non-detection dataset, where non-detection indicates that
species could not be found given search criteria (McGowan and Corwin 2008).
I used the 2000-05 BBA (hereafter, BBA) to quantify species bioclimatic relationships
and to investigate whether birds show spatially-varying (i.e., dependent on the amount of habitat)
responses to climatic variability. I chose 11 grassland birds: American kestrel (Falco sparverius),
38

northern harrier (Circus cyaneus), eastern bluebird (Sialia sialis), eastern kingbird (Tyrannus
tyrannus), horned lark (Eremophila alpestris), grasshopper sparrow (Ammodramus
savannarum), savannah sparrow (Passerculus sandwichensis), vesper sparrow (Pooecetes
gramineus), bobolink (Dolichonyx oryzivorus), brown-headed cowbird (Molothrus ater), and
eastern meadowlark (Sturnella magna). Most of the selected species are grassland generalists
and breed in a diversity of different open habitats, including old fields, pastures, hayfields, and
cultivated fields (DeGraaf and Yamasaki 2001). In 2000-05, the 11 species occupied the
following portion of all BBA blocks: 56% (kestrels), 17% (harriers), 73% (bluebirds), 83%
(kingbirds), 13% (larks), 9% (grasshopper sparrows), 58% (savannah sparrows), 11% (vesper
sparrows), 60% (bobolinks), 75% (cowbirds), and 49% (meadowlarks).
I chose 12 forest species: blue-gray gnatcatcher (Polioptila caerulea), blue-headed vireo
(Vireo solitaries), golden-crowned kinglet (Regulus satrapa), great crested flycatcher
(Myiarchus crinitus), red-breasted nuthatch (Sitta canadensis), tufted titmouse (Baeolophus
bicolor), veery (Catharus fuscescens), black-and-white warbler (Mniotilta varia), blackthroated blue warbler (Setophaga caerulescens), black-throated green warbler (Setophaga
virens), blackburnian warbler (Setophaga fusca), and Canada warbler (Cardellina canadensis).
Most of these forest birds breed in a variety of forest habitats (deciduous, coniferous, and mixed
deciduous-coniferous forests) rather than being restricted to one forest type (DeGraaf and
Yamasaki 2001). In 2000-05, the selected forest birds occupied 20% (gnatcatchers), 52%
(vireos), 22% (kinglets), 79% (flycatchers), 45% (nuthatches), 59% (titmice), and 80% (veeries),
53% (black-and-white warblers), 38% (black-throated blue warblers), 60% (black-throated green
warblers), 40% (blackburnian warblers), and 26% (Canada warblers) of all BBA blocks.

39

2.2.2 Model covariates
The climatic covariates included in all the models were 5-year (2000-05) average values of
maximum and minimum temperatures (TMAX and TMIN, respectively) and average total
precipitation (PRECIP) of the breeding season (April-August). The climate data were derived
from the PRISM (Parameter-elevation Regressions on Independent Slopes Model) climate
mapping system (Daly and Gibson 2002). Though values of TMAX and TMIN are expected to
be correlated with each other (here, Pearson product moment correlation coefficient r = 0.65),
both metrics are biologically meaningful and thus often included together in one model (e.g.,
Stralberg et al. 2009). Correlation between TMAX and PRECIP and TMIN and PRECIP was
much lower (r = -0.18 and r = -0.22, respectively). To account for potential survey effort bias, I
included survey effort (EFF, number of person hours; McGowan and Zuckerberg 2008) as a
variable in all the models.
2.2.3 Habitat amount
Grasslands and forests were identified using the 2001 National Land Cover Data (NLCD)
derived from the Landsat Thematic Mapper satellite data (Homer et al. 2004). I used the 2001
NLCD because it coincided well with the time of the BBA. To derive the amount of the
grasslands, for nine out of 11 species, I consolidated grassland/herbaceous, cultivated crops, and
pasture/hay into one class of grasslands. Due to the more specialized habitat requirements of
some grassland species, I implemented further classification rules to the land cover data. For
bobolinks, I consolidated grassland/herbaceous and pasture/hay into one class (Martin and Gavin
1995, DeGraaf and Yamasaki 2001), while for harriers, I consolidated grassland/herbaceous,
cultivated crops, pasture/hay, and emergent herbaceous wetlands (DeGraaf and Yamasaki 2001,
Smith et al. 2011).

40

To derive the amount of forest within an atlas block, I consolidated deciduous forest,
coniferous forest, and mixed forest into one class of forest. I implemented further land-cover
classification rules for gnatcatchers and titmice, for which I consolidated deciduous and mixed
forest classes (Grubb and Pravasudov 1994, DeGraaf and Yamasaki 2001, Kershner and Ellison
2012), and for kinglets, for which I consolidated coniferous and mixed forest classes (DeGraaf
and Yamasaki 2001, Swanson et al. 2012).
The land-cover was quantified using FRAGSTATS 4.2 (McGarigal et al. 2012) and the
Geospatial Modelling Environment (GME, http://www.spatialecology.com/gme/). We used
Percent Land Cover (PLAND) metric to quantify the amount of grassland (GPLAND) and
forested habitats (FPLAND). GPLAND and FPLAND quantify the proportional abundance of
each patch type in each BBA block (McGarigal et al. 2012).
2.2.4 Statistical analysis
I evaluated two statistical models for each species, both of which included all climate covariates
(TMAX, TMIN, PRECIP) and survey effort (EFF). My first model was a stationary spatial
regression, where I added spatially-structured random effects to the intercept to form a spatiallyvarying intercept (SVI) model. The SVI models accounted for spatial autocorrelation in the
residuals, but did not allow slope coefficients to vary across the spatial domain. Second, I created
a spatially-varying coefficients model (SVC) by adding spatially-structured random effects to the
intercept and slope coefficients associated with each of the climatic covariates and survey effort.
By relaxing the assumption that covariates have the same impact across the domain, the SVC
models allowed for evaluation of potential changes in bioclimatic relationships across space.
For a given BBA grid cell the regression outcome variable y(s) is the detection or nondetection of a particular species (i.e., 1 or 0), where s represents the grid cell's spatial

41

coordinates. We assumed y(s) follows a Bernoulli distribution. For the more general SVC model,
π(y(s)│η(s)) ~ Bernoulli(p(η(s))), where p(η(s)) is the probability of occurrence at s and η(s) is

x(s)β + x(s)w(s) with x(s) representing the vector comprising an intercept and location specific

covariate values. The parameters to be estimated were the global regression coefficients, β = (β0,
βTMAX, βTMIN, βPRECIP, βEFF)T, and associated spatially-varying adjustments, w(s) = (w0(s),
wTMAX(s), wTMIN(s), wPRECIP(s), wEFF(s))T. We used a logit link function η(s) = log[p(s)/{1-p(s)}]
for this model.
I assumed each element in the spatially-varying coefficients vector, w(s), arises from a
spatial Gaussian Process (GP) (Banerjee et al. 2004, Cressie and Wikle 2011). Specifically the jth spatially-varying coefficient at location s is wj(s) ~ GP(0, Cj(s, s’)), where s and s’ are any two

locations within the study area, the spatial covariance Cj(s, s’) = σ2jρj(s, s’; ϕj) with variance
parameter σ2j, correlation function ρj(·; ϕj), and spatial decay parameter ϕj. The exponential
spatial correlation was assumed for ρj(·). The SVI sub-model comprises a single GP on the
intercept, i.e., η(s) is x(s)β + w0(s).

Prior distributions on the remaining parameters completed the Bayesian formulation of
the hierarchical model (e.g., Gelman et al. 2004). The regression coefficients, β j’s followed a
normal distribution N(0, 100), while the variance components σ2j’s were assigned inverse-gamma
IG(2, 1) priors. The spatial decay parameters, ϕj's, followed an informative Uniform prior, with
support that ranged from 1 km to the maximum distance between any two grid cells (i.e.,
approximately 325 km). Parameter estimates were based on post burn-in samples from three
MCMC chains of 10,000 iterations each. Burn-in and convergence were diagnosed using the
Gelman-Rubin diagnostic (Gelman and Rubin 1992). Candidate model fit to the observed data

42

was assessed using the Deviance Information Criterion (DIC; Spiegelhalter et al. 2002). Lower
values of DIC indicate improved model fit.
To investigate whether species bioclimatic relationships are related to the amount of
habitat, I regressed the resulting spatially-varying coefficient estimates of TMAX, TMIN, and
PRECIP against the values of GPLAND and FPLAND (hereafter, TMAX-GPLAND, TMINGPLAND, and PRECIP-GPLAND for grassland birds, and TMAX-FPLAND, TMIN-FPLAND,
and PRECIP-FPLAND for forest birds).
Prior to statistical analysis, I removed all blocks with more than 50% open water
coverage and those that did not have continuous land cover coverage or survey effort data.
Ultimately, I used a total of n = 4,490 blocks for birds of open lands and n = 4,492 blocks for
forest birds to fit the models. I standardized all the covariates to ease the comparison and
interpretation of the coefficient estimates. The SVI models were run in R 2.15.1 statistical
program (R Development Core Team, http://www.r-project.org/) using package spBayes (Finley
and Banerjee 2013). In close collaboration with Dr. Andrew Finley, I wrote the SVC models in
C++ programming language. Summaries of parameter estimates were generated using the R
CODA package (Plummer et al. 2012).
2.3 RESULTS
2.3.1 Model comparison
For all species, the SVC model provided improved fit to the data (Appendix 2.1), suggesting that
birds are responding to climatic variation in a spatially-dependent manner. This spatial
heterogeneity allowed for evaluation of bioclimatic relationships across the changing amount of
habitat.

43

2.3.2 Bioclimatic relationships
For grassland birds, TMAX coefficient estimates resulting from the SVI models were positive
for six species (Table 2.1). For five forest birds, TMAX coefficient estimates were negative for
black-throated blue warblers, black-throated green warblers, blackburnian warblers, vireos,
Canada warblers, and nuthatches (Table 2.1). TMAX estimates were positive for gnatcatchers,
kinglets, flycatchers, titmice, and veeries, suggesting that these birds preferred warmer locations.
As a group, grassland birds occurred in landscapes of higher maximum temperatures than forest
breeding birds (Fig. 2.1).
For most grassland birds, 95% credible intervals of the TMIN coefficient estimates
overlapped zero, with exception of kingbirds for which the estimate was positive (Table 2.1). For
forest birds, TMIN coefficient estimates were negative for seven species, with exception of
kinglets, flycatchers and titmice, for which TMIN was positive. As a group, forest birds in
general occurred in landscapes of lower minimum temperature than grassland birds (Fig. 2.1).
PRECIP coefficient estimates were negative for almost all grassland birds (Table 2.1,
Fig. 2.1). PRECIP coefficient estimates were positive for seven forest birds and negative for one
forest species (Table 1). In general, grassland birds occurred in locations of lower precipitation
than forest birds (Fig. 2.1).

44

Table 2.1 Coefficient estimates resulting from the spatially-varying intercept (SVI) models for
the intercept (INTERCEPT), 2000-05 average maximum temperature of the breeding season
(TMAX), 2000-05 average minimum temperature of the breeding season (TMIN), 2000-05
average total precipitation of the breeding season (PRECIP), and 2000-05 survey effort (EFF) for
grassland and forest birds. Species occurrence data were retrieved from the 2000-05 New York
State Breeding Bird Atlas.
Habitat
Guild

Species

American Kestrel
(Falco sparverius)

Northern Harrier
(Circus cyaneus)

Eastern Bluebird
(Sialia sialis)

Grassland
Birds

Eastern Kingbird
(Tyrannus tyrannus)

Horned Lark
(Eremophila alpestris)

Grasshopper Sparrow
(Ammodramus
savannarum)

Savannah Sparrow
(Passerculus
sandwichensis)

Coefficient
INTERCEPT
TMAX
TMIN
PRECIP
EFF
INTERCEPT
TMAX
TMIN
PRECIP
EFF
INTERCEPT
TMAX
TMIN
PRECIP
EFF
INTERCEPT
TMAX
TMIN
PRECIP
EFF
INTERCEPT
TMAX
TMIN
PRECIP
EFF
INTERCEPT
TMAX
TMIN
PRECIP
EFF
INTERCEPT
TMAX
TMIN
PRECIP
EFF

45

Coefficient Estimates
50%
2.5%
97.5%
-0.277
-0.816
0.353
0.198
-0.044
0.458
0.175
-0.098
0.413
-0.584
-0.827
-0.386
0.281
0.184
0.379
-2.142
-2.704
-1.616
0.086
-0.205
0.364
0.001
-0.268
0.313
-0.438
-0.744
-0.163
0.487
0.387
0.596
0.205
-0.691
1.022
0.531
0.285
0.819
-0.045
-0.343
0.210
-0.399
-0.642
-0.126
0.242
0.152
0.345
2.156
1.623
2.566
0.451
0.124
0.801
0.413
0.094
0.734
-0.770
-1.049
-0.452
0.683
0.489
0.900
-2.588
-3.341
-1.507
-0.008
-0.360
0.322
0.354
-0.016
0.764
-0.617
-0.978
-0.222
0.029
-0.103
0.134
-3.545
-4.075
-2.857
0.235
-0.133
0.613
0.269
-0.087
0.627
-0.291
-0.602
-0.007
0.206
0.112
0.304
0.142
-0.724
0.647
0.288
0.009
0.589
-0.181
-0.491
0.102
-0.577
-0.879
-0.308
0.178
0.096
0.270
Continued on next page

Table 2.1 (cont’d)
Habitat
Guild

Species

Vesper Sparrow
(Pooecetes gramineus)

Bobolink
(Dolichonyx oryzivorus)
Grassland
Birds
Brown-headed Cowbird
(Molothrus ater)

Eastern Meadowlark
(Sturnella magna)

Blue-gray Gnatcatcher
(Polioptila caerulea)

Blue-headed Vireo
(Vireo solitaries)
Forest Birds
Golden-crowned Kinglet
(Regulus satrapa)

Great Crested Flycatcher
(Myiarchus crinitus)

Coefficient
INTERCEPT
TMAX
TMIN
PRECIP
EFF
INTERCEPT
TMAX
TMIN
PRECIP
EFF
INTERCEPT
TMAX
TMIN
PRECIP
EFF
INTERCEPT
TMAX
TMIN
PRECIP
EFF
INTERCEPT
TMAX
TMIN
PRECIP
EFF
INTERCEPT
TMAX
TMIN
PRECIP
EFF
INTERCEPT
TMAX
TMIN
PRECIP
EFF
INTERCEPT
TMAX
TMIN
PRECIP
EFF

46

Coefficient Estimates
50%
2.5%
97.5%
-3.026
-3.768
-2.322
0.025
-0.336
0.431
0.353
-0.038
0.704
-0.277
-0.612
0.035
0.067
-0.060
0.165
-0.680
-1.431
-0.202
0.475
0.219
0.833
-0.334
-0.598
0.056
-0.614
-0.930
-0.311
0.166
0.065
0.267
1.399
0.264
2.059
0.493
0.229
0.808
0.219
-0.121
0.538
-0.544
-0.854
-0.200
0.348
0.197
0.507
0.175
-1.018
1.227
0.326
0.057
0.611
0.020
-0.300
0.295
-0.607
-0.898
-0.289
0.155
0.060
0.236
-2.402
-2.955
-1.620
0.670
0.406
0.952
0.235
-0.042
0.526
-0.457
-0.676
-0.198
0.248
0.156
0.343
-0.684
-1.927
0.182
-0.367
-0.624
-0.094
-0.746
-1.011
-0.459
0.602
0.374
0.830
0.210
0.120
0.309
-2.254
-2.681
-1.812
-0.916
-1.174
-0.678
-0.573
-0.871
-0.249
0.345
0.091
0.589
0.203
0.112
0.296
0.896
-0.003
2.162
0.422
0.194
0.652
0.257
0.020
0.482
-0.231
-0.480
0.040
0.185
0.079
0.303
Continued on next page

Table 2.1 (cont’d)
Habitat
Guild

Species

Red-breasted Nuthatch
(Sitta Canadensis)

Tufted Titmouse
(Baeolophus bicolor)

Veery
(Catharus fuscescens)

Black-and-white Warbler
(Mniotilta varia)
Forest Birds
Black-throated Blue
Warbler
(Setophaga caerulescens)

Black-throated Green
Warbler
(Dendroica virens)

Blackburnian Warbler
(Dendroica fusca)

Canada Warbler
(Cardellina Canadensis)

Coefficient
INTERCEPT
TMAX
TMIN
PRECIP
EFF
INTERCEPT
TMAX
TMIN
PRECIP
EFF
INTERCEPT
TMAX
TMIN
PRECIP
EFF
INTERCEPT
TMAX
TMIN
PRECIP
EFF
INTERCEPT
TMAX
TMIN
PRECIP
EFF
INTERCEPT
TMAX
TMIN
PRECIP
EFF
INTERCEPT
TMAX
TMIN
PRECIP
EFF
INTERCEPT
TMAX
TMIN
PRECIP
EFF

47

Coefficient Estimates
50%
2.5%
97.5%
-1.089
-1.816
-0.437
-0.386
-0.634
-0.150
-0.587
-0.842
-0.333
0.328
0.108
0.548
0.316
0.225
0.404
-0.008
-0.815
0.713
1.119
0.743
1.451
0.247
-0.092
0.557
0.014
-0.311
0.288
0.177
0.064
0.299
1.515
0.799
2.247
0.284
0.024
0.541
-1.000
-1.396
-0.736
-0.122
-0.403
0.222
0.275
0.166
0.389
-0.654
-1.939
0.273
-0.157
-0.337
0.060
-0.236
-0.544
0.002
0.128
-0.101
0.339
0.301
0.200
0.403
-1.402
-2.097
-0.670
-0.540
-1.010
-0.292
-0.799
-1.065
-0.387
0.919
0.657
1.215
0.232
0.140
0.340
0.192
-0.485
0.806
-0.874
-1.173
-0.587
-0.884
-1.182
-0.562
0.752
0.476
0.993
0.188
0.080
0.281
-0.827
-1.308
-0.116
-0.565
-0.884
-0.250
-0.734
-1.067
-0.397
0.845
0.581
1.091
0.225
0.126
0.325
-1.757
-2.247
-1.395
-0.274
-0.553
-0.033
-0.550
-0.849
-0.241
0.414
0.194
0.607
0.143
0.060
0.233

Figure 2.1 Mean values of the coefficient estimates of 2000-05 average spring maximum
temperature (TMAX), 2000-05 average spring minimum temperature (TMIN), and 2000-05
average total spring precipitation (PRECIP) for both avian groups: grassland birds (circles) and
forest birds (triangles). 95% confidence intervals indicate whether the mean value of the
coefficient estimates is significantly different from zero. On each plot from left to right:
American kestrel, northern harrier, eastern bluebird, eastern kingbird, horned lark, grasshopper
sparrow, savannah sparrow, vesper sparrow, bobolink, brown-headed cowbird, eastern
meadowlark, blue-gray gnatcatcher, blue-headed vireo, golden-crowned kinglet, great crested
flycatcher, red-breasted nuthatch, tufted titmouse, veery, black-and-white warbler, black-throated
blue warbler, black-throated green warbler, blackburnian warbler, and Canada warbler. Species
occurrence data were retrieved from the 2000-05 New York State Breeding Bird Atlas.

48

2.3.3 Heterogeneity of bioclimatic relationships across habitat amount
After incorporating covariate-specific spatial random effects in the SVC models, coefficient
estimates of all variables became strongly spatially varying (Appendix 2.2). There was a clear
variation in the values of the TMAX coefficient estimates across the habitat gradient (Table 2.2),
but there was no visible pattern that would differentiate responses of grassland birds from those
of forest species (Fig. 2.2). The strongest relationship was recorded for black-and-white warblers
(posterior mean 0.0089, Table 2.2), followed by larks (posterior mean -0.0082, Table 2.2),
Canada warblers (posterior mean 0.0065, Table 2.2), and harriers (posterior mean 0.0052, Table
2.2). The TMAX-GPLAND and TMAX-FPLAND associations were positive for 13 and
negative for 9 out of 23 species (Table 2.2).
The strongest TMIN-GPLAND and TMIN-FPLAND associations were recorded for
bobolinks (posterior mean -0.0388, Table 2.2), followed by veeries (posterior mean 0.0247,
Table 2.2), flycatchers (posterior mean 0.0231, Table 2.2), and savannah sparrows (posterior
mean -0.0214, Table 2.2). Most of the associations for grassland birds were negative, with
exception of grasshopper sparrows, larks, and vesper sparrows (Table 2.2; Fig. 2.2). The
negative TMIN-GPLAND associations suggested that positive values of TMIN coefficient
estimates were generally found in localities where grasslands were scarce, and declined steadily
towards negative values in regions with abundant grasslands (Appendix 2.2). Thus, only in
regions where grassland habitat was limited, birds were more likely to be found when average
temperatures were higher, whereas we generally saw the reverse relationship in regions with
extensive grasslands.
For forest birds, the TMIN-FPLAND associations were generally positive, with exception
of gnatcatchers, kinglets, nuthatches, and black-throated blue warblers (Table 2.2, Fig. 2.2). The

49

two strongest TMIN-FPLAND associations among forest birds were found for flycatchers and
veeries (Table 2.2), both of which were positive. For these two species, the estimates of TMIN
coefficients were negative in localities with scarce forests and positive in highly forested regions.
Thus, only in regions where forests are scant, flycatchers and veeries are less likely to be found
when average temperatures are higher, whereas we will generally see the reverse relationship in
extensively forested regions. For the remaining species with the positive TMIN-FPLAND
associations, TMIN coefficients were almost exclusively negative (Appendix 2.2). These species
are thus least likely to be found when average temperatures are high, though this relationship is
the strongest in scarcely forested habitats and diminishes with increasing amount of forest.
The eight strongest associations between PRECIP coefficient estimates and land cover
were recorded for grassland birds (in the order of strength: meadowlarks, harriers, bobolinks,
kestrels, kingbirds, vesper sparrows, savannah sparrows, larks; Table 2.2, Fig. 2.2), all of them
positive. Positive associations indicate that PRECIP coefficient estimates were negative in
locations with low GPLAND, suggesting in those places the likelihood of species occurrence
will be low if precipitation is high. With increasing GPLAND, we recorded a weakening of this
relationship (Appendix 2.2). The PRECIP-FPLAND associations for forest species were
generally weaker than those of grassland birds (Fig. 2.2).

50

Table 2.2 Relationships between amount of habitat (percent of grassland cover, GPLAND, and
percent of forest cover, FPLAND, for grassland and forest birds, respectively) and spatiallyvarying coefficient estimates for 2000-05 average maximum temperature of the breeding season
(TMAX), 2000-05 average minimum temperature of the breeding season (TMIN), and 2000-05
average total precipitation of the breeding season (PRECIP). The “*” symbol indicates that less
than 10% of the spatially-varying coefficient estimates (i.e., TMAX, TMIN, PRECIP coefficient
estimates resulting from the SVC models) had the 95% credible intervals not overlapping zero.
Species occurrence data were retrieved from the 2000-05 New York State Breeding Bird Atlas.
Habitat
Guild

Species

Coefficient

TMAX
TMIN
PRECIP
TMAX
Northern Harrier (Circus
TMIN
cyaneus)
PRECIP
TMAX
Eastern Bluebird (Sialia
TMIN
sialis)
PRECIP
TMAX
Eastern Kingbird (Tyrannus
TMIN
tyrannus)
PRECIP
TMAX
Horned Lark (Eremophila
TMIN
alpestris)
PRECIP
TMAX
Grasshopper Sparrow
TMIN*
(Ammodramus savannarum)
PRECIP*
TMAX*
Savannah Sparrow
TMIN
(Passerculus sandwichensis)
PRECIP
TMAX*
Vesper Sparrow (Pooecetes
TMIN
gramineus)
PRECIP
TMAX
Bobolink (Dolichonyx
TMIN
oryzivorus)
PRECIP
TMAX
Brown-headed Cowbird
TMIN
(Molothrus ater)
PRECIP
TMAX
Eastern Meadowlark
TMIN
(Sturnella magna)
PRECIP
American Kestrel (Falco
sparverius)

Grassland
Birds

51

Coefficient Estimates
50%
2.5%
97.5%
-0.0041
-0.0043
-0.0041
-0.0155
-0.0165
-0.0144
0.0158
0.0147
0.0189
0.0053
0.0048
0.0057
-0.0188
-0.0198
-0.0178
0.0281
0.0263
0.0300
-0.0013
-0.0014
-0.0011
-0.0111
-0.0123
-0.0100
0.0033
0.0029
0.0037
0.0015
0.0014
0.0016
-0.0133
-0.0144
-0.0123
0.0141
0.0133
0.0149
-0.0082
-0.0085
-0.0079
0.0051
0.0050
0.0053
0.0055
0.0053
0.0057
-0.0020
-0.0021
-0.0019
0.0025
0.0024
0.0026
-0.00024 -0.00025 -0.00022
-0.0003
-0.0004
-0.0002
-0.0214
-0.0223
-0.0205
0.0105
0.0093
0.0177
0.0012
0.0009
0.0015
0.0023
0.0021
0.0025
0.0106
0.0101
0.0111
0.0026
0.0024
0.0029
-0.0576
-0.0616
-0.0534
0.0256
0.0234
0.0278
-0.0041
-0.0044
-0.0037
-0.0177
-0.0186
-0.0168
-0.0021
-0.0026
-0.0015
0.0048
0.0046
0.0050
-0.0099
-0.0106
-0.0093
0.0316
0.0302
0.0328
Continued on next page

Table 2.2 (cont’d)
Habitat
Guild

Species
Blue-gray Gnatcatcher
(Polioptila caerulea)
Blue-headed Vireo (Vireo
solitaries)
Golden-crowned Kinglet
(Regulus satrapa)
Great Crested Flycatcher
(Myiarchus crinitus)
Red-breasted Nuthatch
(Sitta Canadensis)
Tufted Titmouse
(Baeolophus bicolor)

Forest Birds
Veery (Catharus
fuscescens)
Black-and-white Warbler
(Mniotilta varia)
Black-throated Blue
Warbler (Setophaga
caerulescens)
Black-throated Green
Warbler (Dendroica virens)
Blackburnian Warbler
(Dendroica fusca)
Canada Warbler (Cardellina
Canadensis)

Coefficient
TMAX
TMIN
PRECIP
TMAX
TMIN
PRECIP
TMAX
TMIN
PRECIP
TMAX
TMIN
PRECIP
TMAX
TMIN
PRECIP
TMAX
TMIN
PRECIP
TMAX
TMIN
PRECIP
TMAX
TMIN
PRECIP
TMAX
TMIN
PRECIP
TMAX
TMIN
PRECIP
TMAX
TMIN
PRECIP
TMAX
TMIN
PRECIP

52

Coefficient Estimates
50%
2.5%
97.5%
0.00007
0.00001
0.00014
0.0001
-0.0002
0.0004
-0.00005 -0.00004
0.0003
0.0025
0.0022
0.0029
0.0049
0.0046
0.0051
-0.0032
-0.0035
-0.0029
0.0017
0.0013
0.0020
-0.0013
-0.0019
-0.0008
-0.0017
-0.0033 -0.00009
-0.0017
-0.0021
-0.0014
0.0231
0.0221
0.0241
0.0048
0.0046
0.0050
0.0010
0.0008
0.0012
-0.0013
-0.0015
-0.0011
0.00075
0.0003
0.0012
0.00045
0.00009
0.00089
0.0094
0.0079
0.0108
0.00006 -0.00006
0.00018
0.0042
0.0039
0.0045
0.0247
0.0234
0.0261
-0.0011
-0.0013
-0.0008
0.0089
0.0084
0.0094
0.0013
0.0011
0.0014
-0.0008
-0.0011
-0.0005
-0.0020
-0.0021
-0.0019
-0.00011 -0.00014 -0.00009
0.00164
0.00155
0.00172
-0.0004
-0.0008 -0.00001
0.0037
0.0033
0.0041
-0.0022
-0.0024
-0.0020
0.0037
0.0034
0.0039
0.0119
0.0113
0.0125
-0.0025
-0.0038
-0.0013
0.0065
0.0060
0.0070
0.00603
0.00558
0.00648
-0.0024
-0.0028
-0.0021

Figure 2.2 Mean values of the associations between the percent habitat (percent grassland cover,
GPLAND, and percent forest cover, FPLAND, for grassland and forest species, respectively) and
coefficient estimates of 2000-05 average spring maximum temperature (TMAX), 2000-05
average spring minimum temperature (TMIN), and 2000-05 average total spring precipitation
(PRECIP) for both avian groups: grassland birds (circles) and forest birds (triangles). 95%
confidence intervals indicate whether the mean value of the associations is significantly different
from zero. On each plot from left to right: American kestrel, northern harrier, eastern bluebird,
eastern kingbird, horned lark, grasshopper sparrow, savannah sparrow, vesper sparrow, bobolink,
brown-headed cowbird, eastern meadowlark, blue-gray gnatcatcher, blue-headed vireo, goldencrowned kinglet, great crested flycatcher, red-breasted nuthatch, tufted titmouse, veery, blackand-white warbler, black-throated blue warbler, black-throated green warbler, blackburnian
warbler, and Canada warbler. Species occurrence data were retrieved from the 2000-05 New
York State Breeding Bird Atlas.

53

2.4 DISCUSSION
My findings supported the notion that grassland and forest birds exhibit divergent bioclimatic
relationships. Forest breeding birds were more likely to occur in cooler and wetter regions, while
grassland birds preferred warmer and drier conditions. I also found support for our second
prediction that increasing amounts of forests will mitigate the negative responses of forest
breeding birds to climatic variability. Contrary to my prediction that the magnitude of
bioclimatic relationships would remain constant for grassland birds, I found that the bioclimatic
relationships of grassland birds are highly heterogeneous across the habitat gradient. Notably, the
evidence suggests that the negative consequences of climate change are likely to be intensified in
regions with extensive grasslands. Taken together, my results support our initial hypothesis that
grassland birds will be more vulnerable to a warming climate than birds of forested habitats.
2.4.1 Bioclimatic relationships
Space-varying intercept models allowed comparison of the spatially-homogeneous bioclimatic
relationships of grassland and forest birds. In general, there were differences in the way
grassland and forest birds responded to each of the climatic covariate. Six grassland birds
responded positively to maximum temperatures, indicating that in general these birds prefer
warmer regions. The likelihood of occurrence of seven forest birds was negatively associated
with maximum temperatures, suggesting the higher likelihood of forest birds in cooler locations.
Most forest species also had a strong negative response to minimum temperatures, while almost
none of the grassland birds responded to this covariate. The spatially-homogeneous bioclimatic
relationships suggest that only a portion of grassland birds might benefit from increasing
temperatures resulting from climate change. On the other hand, evidence from stationary models
suggests that most forest species will be negatively affected as temperatures continue to rise. My

54

results therefore corroborate other findings, which suggest that forest birds are generally colddwelling species (Clavero and Brotons 2010) and thus, all else being equal, potentially more
vulnerable to rising temperatures than grassland birds. However, caution is advised here as the
results of the space-varying intercept models do not take into account the potential influences of
habitat and thus lead to conclusions that differ from those arising from the models accounting for
the spatial heterogeneity of bioclimatic relationships.
I found that grassland birds had a strong negative response to precipitation. This finding
fits well with the notion that open habitats generally exhibit drier microclimatic conditions
(Geiger et al. 2009), thus supporting fauna more sensitive to high precipitation regimes (Skagen
and Adams 2012). On the contrary, the likelihood of occurrence of forest birds was positively
affected by precipitation. Thus, my results suggest that precipitation will be an important
climatic variable, albeit less predictable than temperature, in shaping responses of breeding birds
to climate change. Furthermore, given the disparate relationships of grassland and forest birds
with precipitation, we can expect that changing precipitation patterns will likely affect these two
avian groups in a contrasting way. In situations when climate change causes precipitation to
increase, grassland birds will likely be adversely affected, while forest species will generally
benefit from increasing rainfall.
2.4.2 Heterogeneity of bioclimatic relationships across habitat amount
Evaluation of spatial variation in bioclimatic relationships using space-varying coefficients
models allowed a richer ecological interpretation and painted a different picture of potential
vulnerabilities of grassland and forest birds to climate change. My results supported the
expectation that bioclimatic relationships of forest birds exhibit strong spatial variation, which
was associated with the amount of available forest cover. As I predicted, the relationships

55

between forest birds and maximum and minimum temperatures showed mostly positive
associations with the amount of forest cover. Regions with limited forest cover were
characterized by the strongest negative relationships with temperature, and these relationships
became weaker or positive as the amount of forest increased. Thus, more extensively forested
regions might successfully mitigate the increasing ambient temperatures associated with climate
change. Such modulating effect is likely due to the cooler microclimatic conditions created by
the forest canopy closure, which decreases summer ground-level temperatures via increased
shading (De Frenne et al. 2013). These counteracting effects of forest canopies on climate
warming are likely to persist into the future as long as the structure of forest canopies remains
undisturbed. My results are also in agreement with studies suggesting that expansion of forests
will potentially compensate or reverse responses of communities to climate change (Clavero et
al. 2011) by leading to the dominance of cold-dwelling assemblages.
The amount of forest cover in general did not affect the occurrence-precipitation
relationships of forest birds. This finding corroborates the results from the spatially-stationary
framework, where positive relationships of forest birds with precipitation were recorded. The
lack of spatial heterogeneity in the occurrence-forest cover relationship affirms that forest birds
will be more likely to occupy wetter regions regardless of the amount of habitat they have
available. Thus, the effects of precipitation changes resulting from climate change on forest
species will likely be homogenous throughout the species’ geographic extent.
I found that bioclimatic relationships of grassland birds were spatially heterogeneous, but
the associations between occurrence-climate relationships and habitat were different than the
ones found for forest birds. Relationships of grassland species with maximum and minimum
temperatures generally showed mostly negative associations with the amount of grassland cover.

56

Interestingly, responses of grassland birds to maximum temperatures across the increasing
amount of their habitat were not consistent across all species within the group, with some species
showing negative relationships (kestrels, cowbirds, bluebirds, grasshopper sparrows, and larks)
and others having positive associations (bobolinks, kingbirds, meadowlarks, and harriers). For
the species exhibiting positive associations, we see an increase in the strength of the relationship
between likelihood of species occurrence and maximum temperatures as the amount of
grasslands increases. This indicates that birds in localities with abundant grasslands might
benefit from climate change-induced temperature increases more than those found in regions
with scarce suitable cover. Species exhibiting negative associations will likely benefit from
temperature increases in regions with scarce grasslands.
The associations between species responses to minimum temperatures and the amount of
grassland cover were more consistent across all grassland species. I found that likelihood of
occurrence of grassland birds is greater when minimum temperatures are higher, but only in
locations where their habitat is relatively limited. On the other hand, locations characterized by
extensive grasslands show decreasing probability of occurrence associated with higher minimum
temperatures. Thus, I suggest that extensive grasslands might exacerbate, rather than mitigate,
the negative consequences climate change might have on birds of those habitats. One possible
explanation may be that smaller patches of grassland habitat are likely to be surrounded by cover
types with cooler microclimates, thus creating more benign environmental conditions in the face
of increasing temperatures. On the other hand, large and unfragmented grassland patches will be
deprived of the potential counteracting effects of surrounding patches on climate warming and
thus might be more susceptible to rapid temperature increases. Hence, species found in those
habitats will likely experience higher rates of extreme weather events (e.g., heat waves,

57

flooding), potentially making them more vulnerable to the impacts of climate change. This
potential influence of landscape matrix on species bioclimatic relationships is not surprising.
Indeed, studies have shown that the quality of the surrounding matrix influences species
persistence (e.g., Prugh et al. 2008, Willis and Bhagwat 2009) and has a strong influence on
species responses to environmental conditions (Pervedello and Vieira 2010). My results thus
reinforce the need to acknowledge the importance of the entire landscape when evaluating
bioclimatic relationships and their interactions with land cover, rather than focusing only on the
perceived suitable habitat.
Grassland birds also showed strong variation in precipitation relationships across their
habitat gradient. For the majority of grassland species, their occurrence was negatively affected
by high precipitation regimes in regions with scarce grassland habitat. This relationship
weakened with increasing amount of grasslands. Thus, increasing rainfall would be most
detrimental to grassland birds in regions with low amounts of their habitat, whereas populations
in regions with abundant grasslands should be less affected by increasing precipitation. Given the
uncertainty of future precipitation trends (IPCC 2013), forecasting how grassland birds will be
affected by changing precipitation patterns is difficult. However, regardless of whether rainfall
increases or decreases, my study indicates that the direction and magnitude of responses of these
birds will be heavily dependent on the quantity of available habitat.
In summary, I found evidence that grassland and forest birds exhibit divergent
bioclimatic relationships. I also found that the responses of both avian groups to climatic
variability are affected by the amount of habitat, though the influence of land cover differs for
grassland and forest birds. An increasing amount of forest cover might mitigate negative effects
of climate change, especially rising temperatures, in populations of forest breeding birds.

58

However, contrary to my expectations, I found that grasslands will likely exacerbate, rather than
mitigate, the negative consequences of increasing temperatures in locations where this habitat is
prevalent. Our work thus shows that the interaction of climate and land cover can produce
unexpected relationships and needs to be carefully evaluated to make sound forecasts of
consequences of climate change to natural communities.
Grassland birds already show almost universal pattern of decline across both North
American and European continents (Reif 2013). The decline in area of agricultural lands
associated with abandonment of agricultural fields and their conversion to woody successional
stages is one of the most important factors responsible for decline of these birds (Fuller et al.
1995, Brennan and Kuvlesky 2005). Recent studies suggest that climate change will likely add to
declines of grassland birds (e.g., Kleijn et al. 2010). My work corroborates the results from these
studies but also stresses the importance of considering both climate and land cover as synergistic
factors shaping avian responses. In terms of grassland birds, the impacts of both climate and land
cover are difficult to reconcile. While regions of extensive grasslands generally benefit
populations of grassland birds by providing suitable breeding habitat, they are also likely to
increase the vulnerability of birds to changing climatic conditions. Thus, climate change will
threaten grassland birds mostly in habitats, in which they would otherwise be able to persist. My
findings therefore add another dimension to the knowledge of potential vulnerabilities of
grassland birds to the environmental change and stress the importance of considering
implications of climate change while designing future conservation strategies.
2.5 ACKNOWLEDGMENTS
I would like to thank the volunteers who participated in both New York State Breeding Bird
Atlases. I also thank Kimberley Corwin and Kevin McGowan for supplying atlas databases. The

59

manuscript benefited from discussions with the members of the Quantitative Wildlife Laboratory
at Michigan State University. I thank Dr. Colin M. Beier, Daniel Bishop, and John Wiley for
supplying climate data. This study received financial support from NASA and the Boone and
Crockett Club.

60

LANDSCAPE FRAGMENTATION AFFECTS RESPONSES OF AVIAN
COMMUNITIES TO CLIMATE CHANGE

In preparation for publication in Global Change Biology:
Jarzyna M.A. (following order of authors TBD, currently alphabetical) Finley A.O., Maurer B.A,
Porter W.F. and Zuckerberg B. Landscape fragmentation affects responses of avian
communities to climate change. To be submitted to Global Change Biology.
ABSTRACT
Accurate forecasting of the consequences of climate change is contingent upon our
understanding of the relationship between biodiversity patterns and climatic variability. While a
range of implications of climate change to individual species have been well documented,
understanding of the community-level consequences of changing climate is still limited.
Additionally, factors other than climate will affect the way communities respond to climate
change. My objectives were to investigate the relationship between temporal turnover in avian
biodiversity and changes in climatic conditions, and to assess the roles of landscape
fragmentation and migratory behavior as factors affecting this relationship. Using data from the
1980 and 2000 New York State Breeding Bird Atlas, I quantified temporal turnover in the entire
avian biodiversity and in diversity of different migratory groupings. I applied Bayesian spatiallyvarying intercept models to evaluate the relationship between temporal turnover and covariates
of temporal trends in climatic conditions as well as landscape fragmentation. The DIC statistic
designated the models including interaction terms between climate change and landscape
fragmentation to be superior to the models without the interaction terms, suggesting that the
relationship between temporal turnover and changes in climatic conditions was affected by the
level of landscape fragmentation. Specifically, I found weaker associations between temporal
turnover and climatic change in regions with prevalent habitat fragmentation. I suggest that avian
communities in fragmented landscapes are more robust to climate change than communities
61

found in contiguous habitats because they comprise of species with wider thermal niches and are
thus less susceptible to temperature increases. I conclude that highly fragmented regions are like
to undergo less pronounced changes in composition and structure of faunal communities as a
result of climate change, whereas those changes are likely to be greater in contiguous and
unfragmented habitats.
3.1 INTRODUCTION
Accurate forecasting of the consequences of climate change and sound conservation strategies
resulting from these forecasts are contingent upon our understanding of the relationship between
biodiversity patterns and climatic variability. To date, the majority of research regarding the
implications of climate change to biodiversity has evaluated responses of individual species (e.g.,
Parmesan and Yohe 2003, La Sorte and Thompson 2007, Marini et al. 2009, Zuckerberg et al.
2009). Variation in the individual species responses is predicted to lead to disruptions of
communities and ecosystems (Brotons and Jiguet 2010), but the complex nature of ecological
interactions makes it difficult to extrapolate from the scales of individuals and populations to the
community or ecosystem level (Walther et al. 2002). However, in order to fully understand
consequences of climate change to ecological communities it is imperative that a direct link
between changes in community composition and changes in climatic conditions is established.
Measuring species turnover has potential to offer a useful and simple way of evaluating
the implications of climate change to ecological communities (e.g., Tuomisto 2010a, Tuomisto
2010b, Barton et al. 2013). Species turnover can be quantified in multiple ways but generally is
an aggregate of gains and losses of species comprising a community (e.g., Tuomisto 2010b).
Traditionally, turnover has been used to describe biodiversity patterns across space rather than
time (i.e., spatial turnover), but it can also be viewed in terms of a temporal change in

62

community composition that has occurred at the same site across a particular time period (i.e.,
temporal turnover). Given that spatial turnover in species composition is thought to be at least
partially related to spatial variation in environmental variables (Gaston et al. 2007, Buckley and
Jetz 2008, Gutiérrez-Cánovas et al. 2013, Siefert et al. 2013), it is expected that temporal
turnover is related to temporal changes such as climate and thus could be a useful indicator of
impacts of climate change on biodiversity.
Even though impacts of climate change on biodiversity will be felt across a wide range of
ecosystems, different communities are likely to respond to climate change in a disparate way
because a wide range of environmental factors will influence species’ vulnerabilities to climate
change. For example, past research suggests that the interaction between climate and land cover
will affect the way in which biodiversity responds to climate change (Opdam and Wascher 2004;
Jeltsch et al. 2011). Specifically, habitat fragmentation is thought to affect the rate of
distributional shifts associated with climate change, allowing species in less fragmented habitats
to disperse faster and farther (Opdam and Wascher 2004). If species in contiguous (i.e.,
unfragmented) habitats respond quicker to climate change than species in fragmented habitats,
then changes in communities should also follow the same pattern. Thus, higher rates of temporal
turnover in contiguous habitats are expected compared to temporal turnover in fragmented
habitats. Stronger association of the temporal turnover with changing climatic conditions in
contiguous habitats is also expected.
Beyond influences of land cover, some behavioral traits are likely to affect species’
responses to climate change (e.g., Jones and Cresswell 2010). Migratory strategy is one of such
traits (Lehikoinen and Sparks 2010). For example, some resident species of birds are expected to
respond relatively fast to changing climate on their breeding and wintering grounds because they

63

are exposed to year-round climatic conditions of the same site. Similarly, short-distance migrants
are likely to adjust the timing of migration to changing temperatures because they are thought to
rely more on weather conditions. On the other hand, long-distance migratory birds are cued by
photoperiod (Forchhammer et al. 2002, Saino et al. 2009, DeLeon et al. 2011), thus are unlikely
to be able to track changing climatic conditions as closely as short-distance migrants or resident
species (e.g., Lemoine et al. 2007, Doxa et al. 2012). Indeed, resident birds have undergone
increases in occupancy and abundance in recent decades, while migratory birds (long-distance
migrants in particular) have experienced population losses (e.g., Lemoine and Böhning-Gaese
2003). However, the consensus what proportion of the decline in migratory bird populations, if
any, can be attributed to climate change has not been reached. If migratory traits indeed affect
birds’ responses to climate change, we can expect different migratory groups to show different
associations with changing climate. Specifically, temporal turnover in long-distance migrant
communities is expected to show weaker association with changes in climatic conditions than
temporal turnover in short-distance migrants or resident birds communities.
Despite some initial investigations into recent temporal changes in community
composition (e.g., Kampichler et al. 2012), the relationship between temporal turnover in
biodiversity and the interaction of climatic changes, land-cover, and migratory behavior, to the
best of my knowledge, has not yet been evaluated. The paucity of research in this arena has
likely been caused by the scarcity of long-term, large-scale ecological data. I intend to fill this
research gap by investigating the drivers of temporal turnover in avian biodiversity across New
York using Breeding Bird Atlas data (Andrle and Carroll 1988, McGowan and Corwin 2008).
Specifically, I sought to evaluate two hypotheses: (1) temporal turnover in avian assemblages is
related to climate change, but the strength of the relationship between temporal turnover and

64

changes in climatic conditions varies with the degree of habitat fragmentation, and (2) temporal
turnover in communities of long-distance migrants show weaker association with changing
climatic conditions than temporal turnover in short-distance migrants and resident birds
communities.
3.2 METHODS
3.2.1 Site description
The study area is the State of New York, USA. New York covers 128,401 km2, including 4,240
km2 of inland water. Climate of New York is affected by the state’s broad elevation gradient as
well as by the proximity to lakes Erie and Ontario, and the Atlantic Ocean. Adirondack and
Catskill Mountains in the east are among the highest regions of the state, with elevation varying
between 600 and 1,500 m. The elevation of the south-western New York ranges from
approximately 300 to 900 m, while north-western New York as well as Hudson River valley and
New York City are low-lying with elevation ranging from approximately 0 to 200 m.
In general, New York state is covered by approximately 51.2% forest (deciduous,
evergreen, and mixed forests), 13.4% pasture and hay, 8.2% cultivated crops, 2.9% scrub and
shrub, 1.0% grassland and herbaceous, 8.0% wetland cover, 8.7% developed land, and 0.2%
barren land (Homer et al. 2004). New York offers a broad gradient of landscape fragmentation.
The regions of the Adirondack and Catskill Mountains are the least fragmented, with expansive
regions of forest land. The land cover of the Hudson River valley is fragmented by residential
and commercial development, whereas the landscapes of western New York are characterized
mostly by agriculture-forest mosaic.

65

3.2.2 Breeding Bird Atlas
I used the New York State Breeding Bird Atlas (BBA) to characterize changes in avian
communities through time. BBA is a statewide survey that documented the distribution of
breeding birds in New York. To date, BBA has been conducted in two time periods, 1980-85
(hereafter, BBA1980; Andrle and Carroll 1988) and 2000-05 (hereafter, BBA2000; McGowan
and Corwin 2008). For both BBAs, a grid system was used to define the basic unit for reporting
data. The entire State of New York was divided into approximately 1,300 squares, each
measuring 10 km x 10 km. The BBA reporting unit (a block) was one quarter of this square and
measures 5 km x 5 km and a total of 5,335 blocks covered the entirety of New York State. This
data set represents one of the largest and finest resolution atlases in the world (Gibbons et al.
2007).
A total of 242 species were recorded for BBA1980 and 248 species were recorded for
BBA2000 (Andrle and Carroll 1988, McGowan and Corwin 2008). Observations were made by
skilled birders who spent at least 8 hours in each block, visited all cover types in each block, and
included at least one nighttime visit to document nocturnal species. Observer effort was recorded
for each BBA and reported as a number of person hours (measured as the sum of the number of
hours spent in each block x the number of people surveying each block; McGowan and
Zuckerberg 2008). The BBA represents a detection/non-detection dataset, although nondetection indicates that species could not be found given search criteria (McGowan and Corwin
2008).

66

3.2.3 Temporal turnover
I quantified temporal turnover (TURN) as an aggregate of species gains (hereafter called
colonization) and species losses (hereafter called extinction) within a particular site. Specifically,
temporal turnover was calculated as follows:
𝑇𝑈𝑅𝑁 =

𝐸+𝐶

𝐸+𝐶+𝑃

,

Eqn. 3.1

where E is the number of species that went extinct in the block between BBA1980 and
BBA2000, C is the number of species that colonized the block between BBA1980 and
BBA2000, and P is number of species in the same block common to both BBAs (i.e.,
persistence).
In addition to temporal turnover, I quantified its components: proportion of colonization
events (i.e., proportion of species gained a particular site between BBA1980 and BBA2000;
COL) and proportion of extinction events (i.e., proportion of species lost a particular site
between BBA1980 and BBA2000; EXT) in each block between BBA1980 and BBA2000. COL
and EXT were calculated as follows:
𝐶𝑂𝐿 =

𝐸𝑋𝑇 =

𝐶

𝐸+𝐶+𝑃
𝐸

𝐸+𝑃

,

Eqn. 3.2

,

Eqn. 3.3

For TURN and COL, I used the sum of E, C, and P as the denominator of the Eqns. 3.1

and 3.2. I recognize that species already present in a block in BBA1980 could not have colonized
that block and that species other than those included in E, C, and P metrics could have
potentially colonized a site. However, I was simply interested in the proportion of species that
colonized the site relative to all the species found within the block; thus I deem the sum of E, C,
and P as an appropriate choice for the denominator of both equations. For EXT, I used the sum
of E and P as the denominator, because only species present in a block in BBA1980 could have
67

gone extinct. The values of TURN, COL, and EXT are bounded by 0 and 1; values approaching
0 indicate low temporal turnover, colonization, or extinction in a block between BBA1980 and
BBA2000, while values approaching 1 indicate high temporal turnover, colonization, or
extinction in a block between BBA1980 and BBA2000.
I calculated TURN, EXT, and COL for four groups of birds: all species regardless of
their migratory strategy (hereafter, all species), long-distance migrants, short-distance migrants,
and resident species. I used The Birds of North America Online to classify birds to migratory
groupings (http://bna.birds.cornell.edu/bna/, Appendix 1.1). I considered a species a resident if at
least a part of its wintering range is located in New York. I considered a species a long-distance
migrant if it crosses a body of water (i.e., either the Gulf of Mexico or both the Gulf of Mexico
and the Caribbean Sea) to reach its wintering grounds.
3.2.4 Habitat fragmentation
Habitat fragmentation was identified using the National Land Cover Data (NLCD) derived from
the Landsat Thematic Mapper satellite data. Even though the NLCD is available for three time
periods (i.e., 1992, 2001, and 2006), there is no land-cover data readily available for the time
period of BBA1980. Thus, I used a space-for-time substitution to assess whether temporal
turnover is related to landscape fragmentation (e.g., Pickett 1989). Space-for-time substitution
assumes that spatial variation is equivalent to a temporal variation, and relationships between
variables (e.g., species turnover and land-cover) derived for one time step across a large spatial
extent will be equivalent to that derived for two time steps. I used the 2001 NLCD, because it
coincided well with the time of BBA2000.
Prior to landscape analysis, I consolidated the following land cover classes: open water
and perennial ice and snow classes into one class of water/ice; open space developed, low

68

intensity developed, medium intensity developed, high intensity developed into one class of
developed land; cultivated crops and pasture/hay into one class of agriculture; and deciduous
forest, evergreen forest, mixed forest, and forested wetland into one class of forest.
Consolidation of these land-cover types improved accuracy of classification and simplified the
environment for purposes of our evaluation. Upon consolidation, I examined the following land
cover classes: water/ice, developed, barren land, forest, scrub/shrub, grassland/herbaceous,
agriculture, and wetlands.
I used FRAGSTATS 4.1 (McGarigal et al. 2012) and the Geospatial Modelling
Environment (GME, http://www.spatialecology.com/gme/) to quantify landscape fragmentation.
Because I focused my analysis on a diverse suite of species with varying habitat requirements, I
chose a landscape-scale variable to capture broad-scale variation in habitat fragmentation. I
chose Edge Density (ED) as a measure of landscape fragmentation because an increase in habitat
edge is a primary outcome of habitat fragmentation (Hargis et al. 1998). ED was also reported as
an effective tool for evaluating landscape fragmentation and performed better than other popular
landscape fragmentation indices (Hargis et al. 1998). However, in situations when landscape
consists entirely of one cover type, the ED would be 0 regardless of the type of land-cover
present. Therefore, in order to differentiate between landscapes that consisted entirely of natural
land cover (e.g., forest) and those mostly developed (e.g., cities), I also calculated the proportion
of developed land (DEVEL) metric and included it as a covariate in our models.
3.2.5 Climatic trends
The climate data was derived from the PRISM (Parameter-elevation Regressions on Independent
Slopes Model) climate mapping system (Daly and Gibson 2002). PRISM consists of interpolated

69

monthly maximum and minimum temperatures and precipitation at a 2.5-arcmin resolution from
1891-2010 for the entire contiguous United States.
I calculated the magnitude of the 25-year (1980-2005) trend in average monthly
maximum and minimum temperatures and in total monthly precipitation using Ordinary Least
Squares regression. The slope of the OLS regression indicates the magnitude of the trend and
reflects the amount of change in climatic variables that occurred between 1980 and 2005. I then
interpolated the trend magnitudes for each BBA block and averaged the monthly values to reflect
breeding season (May through September) trend magnitude in maximum and minimum
temperatures (TMAX and TMIN, respectively) and in total precipitation (PRECIP). The trend
magnitudes of TMAX and TMIN were expressed in °C/25 years, while the units of the trend
magnitudes of PRECIP were in mm/25 years.
3.2.6 Survey effort
Increasing survey effort often results in a higher number of recorded species (Tobler et al. 2008).
Drastically different survey effort between BBA1980 and BBA2000 could increase the values of
TURN, EXT, and COL indices indicating high temporal turnover, which could be a result of
different number of recorded species rather than actual changes in species identities throughout
time. To account for potential survey effort bias, I calculated the absolute difference in the
number of person hours between BBA1980 and BBA2000 (EFF = |EFF1980 − EFF2000|) for

each BBA block and included it as a covariate in the models.
3.2.7 Statistical analysis

For each of the four groupings (all species, long-distance migrants, short-distance migrants, and
resident species), I evaluated five competing statistical models. All five models included all
climatic variables (TMAX, TMIN, and PRECIP) as well as survey effort (EFF), one model
70

additionally included ED, one included interaction terms between all climatic covariates and ED,
one included both ED and DEVEL, and one included interaction terms between all climatic
covariates and ED as well as the main effect of DEVEL. I standardized all the covariates to ease
the comparison and interpretation of the coefficient estimates.
We used a Bayesian model estimation procedure. Each response variable (temporal
turnover TURN, extinction EXT, or colonization COL) at each location (s) was modeled as a
binomial random variable. Let y(s) be a response variable, then y(s) ~ B(ρ(s),N(s)) where N(s) is
the total number of species found at least once in a block from the first to the second atlas count
(for TURN and COL) or the total number of species found in first atlas count (for EXT). For
each BBA block, let p(s) = y(s)/N(s). Also, let x(s) be a vector containing the explanatory
variables for each model. All models had spatially varying intercepts, which were created by
adding spatially structured random effects to the model intercept. A generalized linear model
approach was used to relate responses to explanatory variables. Using the logit link η(s) =
log[p(s)/{1-p(s)}] the SVI model was
𝜂(𝑠) = 𝑤(𝐬) + 𝜷𝐱(𝑠)

Eqn. 3.4

where β is a vector of regression coefficients and w(s) is a random effect representing

spatially varying adjustments to the intercept. We assume each element in w(s) arises from a
spatial Gaussian Process (GP) (see, e.g., Banerjee et al. 2004 or Cressie and Wikle 2011 for more
details). Specifically, w(s) ~ GP(0, C(s, s’)), where s and s’ are any two locations within the

study area, the spatial covariance C(s, s’) = σρ(s, s’; ϕ) with variance parameter σ2, correlation
function ρ(·; ϕ), and spatial decay parameter ϕ. The exponential spatial correlation was assumed
for ρ(·).

71

Prior distributions on the parameters complete the hierarchical model specification (e.g.,
Gelman et al. 2004). The global regression coefficients, β’s followed a Normal distribution N(0,
100), while the variance component σ2 was assigned inverse-Gamma IG(2, 1) priors. The spatial
decay parameter, ϕ, followed an informative Uniform prior, with support that ranged from 1 km
to the maximum distance between any two grid cells (i.e., approximately 350 km). I ran three
MCMC chains for 20,000 iterations each. Convergence was diagnosed using the Gelman-Rubin
diagnostic (Gelman and Rubin 1992). As appropriate for a Bayesian framework, candidate model
fit to the observed data was assessed using the Deviance Information Criterion (DIC;
Spiegelhalter et al. 2002). Lower values of DIC indicate improved fit.
Prior to statistical analysis, I removed all blocks that did not have continuous land cover
coverage, blocks with more than 50% open water coverage, and those that did not have survey
effort data. I withheld 10% of the data for model validation. Ultimately, I used a total of N =
4,271 blocks to fit the models and 473 blocks to validate the models. The models were run in R
2.15.1 statistical package (R Development Core Team 2013, http://www.r-project.org/) using
package spBayes (Finley and Banerjee 2013). In a Bayesian framework, parameter estimates
have valid posterior distributions and inferences are made using mean or median as well as the
credible intervals of these posterior distributions. Credible intervals overlapping zero imply that
the parameter in question is not different from zero. Summaries of parameter estimates were
generated using the R CODA package (Plummer et al. 2012).

72

3.3 RESULTS
3.3.1 Relationships with climate change and habitat fragmentation
3.3.1.1 Temporal turnover
For all species, the mean value of TURN was 0.384 (95% confidence intervals resulting from
10,000 non-parametric bootstrap samples: LCI 0.381, UCI 0.387). This indicates that on average
approximately 62% of the species were common to both BBAs and 38% of species had either
colonized or become extinct within the average block (Fig. 3.1). The north-eastern part of the
state, which is the location of the Adirondack Mountains, experienced the highest rates of
temporal turnover across BBA1980 and BBA2000 (Fig. 3.1). Central and north-western New
York underwent lowest temporal turnover (Fig. 3.1).
Model selection indicated that the variation in temporal turnover of all species was best
accounted for by the most complex model (Table 3.1). The selection of the most complex model
points to the importance of the interactions between landscape fragmentation and change in
climatic conditions in explaining temporal turnover in species composition between BBA1980
and BBA2000. Model validation indicated relatively good abilities of the best model to predict
temporal turnover in all species, though low values of TURN tended to be over-predicted, while
high value tended to be under-predicted (Appendix 3.1).

73

Figure 3.1 Temporal turnover, extinction, and colonization observed between 1980-85 and 2000-05 for communities of (A) all
recorded species regardless of their migratory status, (B) long-distance migrants, (C) short-distance migrants, and (D) resident birds.
Temporal turnover was calculated as a proportion of all species gained or lost between 1980-85 and 2000-05 in a particular site (i.e.,
Breeding Bird Atlas block) relative to all species recorded across both time periods; extinction was calculated as a proportion of
species lost between 1980-85 and 2000-05 in a particular site relative to all species present in a block in 1980-85; colonization was
calculated as a proportion of species gained between 1980-85 and 2000-05 in a particular site relative to all species recorded across
both time periods. High values of all three metrics indicate high temporal turnover, extinction, and colonization.

74

Table 3.1 Comparison of five competing models of temporal turnover for all four avian
groupings (i.e., all species, long-distance migrants, short-distance migrants, and resident
species). Models included the following covariates: magnitude of the 25-year (1980-2005) trend
in average maximum temperature of the breeding season (TMAX), magnitude of the 25-year
(1980-2005) trend in average minimum temperature of the breeding season (TMIN), magnitude
of the 25-year (1980-2005) trend in average total precipitation of the breeding season (PRECIP),
edge density (ED), percent developed land (DEVEL), interaction between the magnitude of the
25-year (1980-2005) trend in average maximum temperature of the breeding season and edge
density (TMAX*ED), interaction between the magnitude of the 25-year (1980-2005) trend in
average minimum temperature of the breeding season and edge density (TMIN*ED), interaction
between the magnitude of the 25-year (1980-2005) trend in average total precipitation of the
breeding season and edge density (PRECIP*ED), and survey effort (EFF). Model comparison
was done using Deviance Information Criterion (DIC); pD indicates the effective number of
parameters.
Grouping

All Species

Long-distance
Migrants

Short-distance
Migrants

Resident
Species

Model

ΔDIC

pD

TMAX*ED + TMIN*ED + PRECIP*ED +
DEVEL + EFF
TMAX + TMIN + PRECIP + ED + DEVEL +
EFF
TMAX*ED + TMIN*ED + PRECIP*ED + EFF
TMAX + TMIN + PRECIP + ED + EFF
TMAX + TMIN + PRECIP + EFF
TMAX*ED + TMIN*ED + PRECIP*ED +
DEVEL + EFF
TMAX + TMIN + PRECIP + ED + DEVEL +
EFF
TMAX*ED + TMIN*ED + PRECIP*ED + EFF
TMAX + TMIN + PRECIP + ED + EFF
TMAX + TMIN + PRECIP + EFF
TMAX + TMIN + PRECIP + ED + DEVEL +
EFF
TMAX*ED + TMIN*ED + PRECIP*ED +
DEVEL + EFF
TMAX + TMIN + PRECIP + ED + EFF
TMAX*ED + TMIN*ED + PRECIP*ED + EFF
TMAX + TMIN + PRECIP + EFF
TMAX*ED + TMIN*ED + PRECIP*ED +
DEVEL + EFF
TMAX*ED + TMIN*ED + PRECIP*ED + EFF
TMAX + TMIN + PRECIP + ED + DEVEL +
EFF
TMAX + TMIN + PRECIP + ED + EFF
TMAX + TMIN + PRECIP + EFF

0.0

106.7

20.0

104.3

33.9
55.9
109.9
0.0

106.2
103.6
102.6
96.1

10.8

93.6

53.9
69.1
72.3
0.0

95.1
93.3
92.3
88.9

2.3

92.5

22.8
25.9
57.1
0.0

89.3
92.4
83.5
95.0

1.6
7.5

94.5
91.4

12.6
103.6

91.4
93.6

75

Parameter estimates of the top model indicated that temporal turnover was positively
related to TMAX, suggesting that higher temporal turnover was observed in regions where
maximum temperatures increased the most (Fig. 3.2; also see Appendix 3.2 for parameter
estimates of all models of temporal turnover). Temporal turnover was also positively associated
with TMIN, but the effect of minimum temperatures on TURN was affected by the value of ED
(Fig. 3.3). In regions of low landscape fragmentation, large increases in minimum temperatures
were associated with high temporal turnover, but the effect of TMIN diminished or reversed in
locations with high landscape fragmentation. There was a negative association of temporal
turnover and ED, but the presence of multiple ED interaction terms in the model makes it
difficult to interpret the main effects of ED on temporal turnover. There was a positive
relationship between temporal turnover and DEVEL (Fig. 3.2), indicating that highly urbanized
regions generally experienced increased temporal turnover.

76

Figure 3.2 Coefficient estimates (i.e., mean values of the posterior distribution and associated credible intervals) of the best model for
each avian grouping for (A) temporal turnover, (B) extinction, and (C) colonization. Abbreviations for the covariates are as follows:
magnitude of the 25-year (1980-2005) trend in average maximum temperature of the breeding season (TMAX), magnitude of the 25year (1980-2005) trend in average minimum temperature of the breeding season (TMIN), magnitude of the 25-year (1980-2005) trend
in average total precipitation of the breeding season (PRECIP), edge density (ED), percent developed land (DEVEL), interaction
between the magnitude of the 25-year (1980-2005) trend in average maximum temperature of the breeding season and edge density
(TMAX-ED), interaction between the magnitude of the 25-year (1980-2005) trend in average minimum temperature of the breeding
season and edge density (TMIN-ED), interaction between magnitude of the 25-year (1980-2005) trend in average total precipitation of
the breeding season and edge density (PRECIP-ED).

77

Figure 3.3 Effects size plots resulting from the best temporal turnover models for (A) the interaction of magnitude of the 25-year
(1980-2005) trend in average minimum temperature of the breeding season (TMIN) and edge density (ED) for all species, (B) the
interaction of magnitude of the 25-year (1980-2005) trend in average minimum temperature of the breeding season (TMIN) and edge
density (ED) for long-distance migrants, (C) the interaction of magnitude of the 25-year (1980-2005) trend in total precipitation of the
breeding season (PRECIP) and edge density (ED) for long-distance migrants, (D) the interaction of magnitude of the 25-year (19802005) trend in average maximum temperature of the breeding season (TMAX) and edge density (ED) for resident birds, and (E) the
interaction of magnitude of the 25-year (1980-2005) trend in average minimum temperature of the breeding season (TMIN) and edge
density (ED) for resident birds. The interaction terms shown here are the ones whose credible intervals did not overlap zero. Green
points are the observed data points and are plotted to indicate levels of confidence for the predicted surface.

78

3.3.1.2 Extinction
The mean value of EXT for all species was 0.210 (LCI 0.206, UCI 0.213), which indicates that
on average approximately 21% of species of the original assemblage had become extinct within
the average block (Fig. 3.1). North-eastern part of the state (i.e., the Adirondack Mountains) and
regions surrounding the Finger Lakes displayed higher than average rates extinction (Fig. 3.1).
Extinction was lower in remaining parts of New York.
Model selection indicated that the most complex model provided the best fit to the
variation in the proportion of extinction events (Table 3.2). Selection of the most complex model
indicates that interactions between land cover fragmentation and change in climatic conditions
contribute to explaining extinction rates observed in avian communities. Model validation
indicated the best model in general under-predicted the low values of extinction and overpredicted high values of extinction (Appendix 3.1).
Parameter estimates of the top model indicated that extinction in communities of all
species was positively related to TMIN, though the effect of minimum temperatures on
extinction rates was affected by ED (Fig. 3.4; also see Appendix 3.3 for parameter estimates
resulting from all models of extinction). In regions of low landscape fragmentation, large
increases in minimum temperatures were associated with high proportion of extinction events,
but the effect of TMIN diminished or reversed in locations with high landscape fragmentation
(Fig. 3.4). Extinction was positively associated with PRECIP (Fig. 3.2), indicating that regions
with increasing precipitation experienced higher rates of extinction. This positive relationship
between extinction and PRECIP was not influenced by ED. There was a negative association of
extinction and ED, but the presence of multiple ED interaction terms in the model makes it
difficult to interpret the main effects of ED on extinction in long-distance migrant communities.

79

There was a positive relationship between extinction and DEVEL (Fig. 3.2), indicating increased
proportion of extinction events in highly urbanized regions.

80

Table 3.2 Comparison of five competing models of extinction for all four avian groupings (i.e.,
all species, long-distance migrants, short-distance migrants, and resident species). Models
included the following covariates: magnitude of the 25-year (1980-2005) trend in average
maximum temperature of the breeding season (TMAX), magnitude of the 25-year (1980-2005)
trend in average minimum temperature of the breeding season (TMIN), magnitude of the 25-year
(1980-2005) trend in average total precipitation of the breeding season (PRECIP), edge density
(ED), percent developed land (DEVEL), interaction between the magnitude of the 25-year
(1980-2005) trend in average maximum temperature of the breeding season and edge density
(TMAX*ED), interaction between the magnitude of the 25-year (1980-2005) trend in average
minimum temperature of the breeding season and edge density (TMIN*ED), interaction between
the magnitude of the 25-year (1980-2005) trend in average total precipitation of the breeding
season and edge density (PRECIP*ED), and survey effort (EFF). Model comparison was done
using Deviance Information Criterion (DIC); pD indicates the effective number of parameters.
Grouping

All Species

Long-distance
Migrants

Short-distance
Migrants

Resident
Species

Model

ΔDIC

pD

TMAX*ED + TMIN*ED + PRECIP*ED +
DEVEL + EFF
TMAX*ED + TMIN*ED + PRECIP*ED + EFF
TMAX + TMIN + PRECIP + ED + DEVEL +
EFF
TMAX + TMIN + PRECIP + ED + EFF
TMAX + TMIN + PRECIP + EFF
TMAX*ED + TMIN*ED + PRECIP*ED +
DEVEL + EFF
TMAX + TMIN + PRECIP + ED + DEVEL +
EFF
TMAX*ED + TMIN*ED + PRECIP*ED + EFF
TMAX + TMIN + PRECIP + EFF
TMAX + TMIN + PRECIP + ED + EFF
TMAX*ED + TMIN*ED + PRECIP*ED +
DEVEL + EFF
TMAX + TMIN + PRECIP + ED + DEVEL +
EFF
TMAX + TMIN + PRECIP + ED + EFF
TMAX*ED + TMIN*ED + PRECIP*ED + EFF
TMAX + TMIN + PRECIP + EFF
TMAX*ED + TMIN*ED + PRECIP*ED +
DEVEL + EFF
TMAX*ED + TMIN*ED + PRECIP*ED + EFF
TMAX + TMIN + PRECIP + ED + EFF
TMAX + TMIN + PRECIP + ED + DEVEL
+EFF
TMAX + TMIN + PRECIP + EFF

0.0

118.2

33.4
45.1

117.9
115.3

75.7
165.8
0.0

114.6
115.9
110.6

22.8

107.3

49.3
78.6
78.7
0.0

108.9
106.7
106.6
97.1

4.7

94.4

15.9
17.4
77.15
0.0

93.4
97.3
94.8
104.5

2.3
20.5
20.9

102.5
99.5
101.5

93.4

100.7

81

Figure 3.4 Effects size plots resulting from the best extinction models for (A) the interaction of magnitude of the 25-year (1980-2005)
trend in average minimum temperature of the breeding season (TMIN) and edge density (ED) for all species, (B) the interaction of
magnitude of the 25-year (1980-2005) trend in average minimum temperature of the breeding season (TMIN) and edge density (ED)
for long-distance migrants, (C) the interaction of magnitude of the 25-year (1980-2005) trend in total precipitation of the breeding
season (PRECIP) and edge density (ED) for long-distance migrants, (D) the interaction of magnitude of the 25-year (1980-2005) trend
in average minimum temperature of the breeding season (TMIN) and edge density (ED) for short-distance migrants, and (E) the
interaction of magnitude of the 25-year (1980-2005) trend in average minimum temperature of the breeding season (TMIN) and edge
density (ED) for resident birds. The interaction terms shown here are the ones whose credible intervals did not overlap zero. Green
points are the observed data points and are plotted to indicate levels of confidence for the predicted surface.

82

3.3.1.3 Colonization
For all species, the mean value of COL was 0.214 (LCI 0.211, UCI 0.218), indicating that on
average 21% of all the species present in a block across both time periods are the ones that had
colonized the site. Colonization was in general higher in eastern and western parts of the state,
while central and southern New York underwent lower colonization rates (Fig. 3.1). However,
the spatial pattern of colonization rates was less clumped than the one displayed by temporal
turnover or extinction.
Model selection procedure selected more than one model as the top model accounting for
the variation in the proportion of colonization events. The top two models included the most
complex model (i.e., containing all the interaction terms) and the model with main effects of all
covariates but without the interaction terms (Table 3.3). However, despite the fact that the most
complex model was selected among the top models, none of the interaction terms was associated
with COL (see Appendix 3.4 for parameter estimates resulting from all models of colonization).
Model validation indicated the top models under-predicted the low values of colonization and
over-predicted high values of colonization (Appendix 3.1). Colonization was positively related to
TMAX and negatively related to PRECIP, indicating that highest proportions of colonization
events were found in regions where maximum temperatures increased the most and precipitation
decreased the most (Fig. 3.2). There was a positive relationship between colonization and
DEVEL (Fig. 3.2), which indicates that highly urbanized regions experienced highest
colonization rates.

83

Table 3.3 Comparison of five competing models of colonization for all four avian groupings
(i.e., all species, long-distance migrants, short-distance migrants, and resident species). Models
included the following covariates: magnitude of the 25-year (1980-2005) trend in average
maximum temperature of the breeding season (TMAX), magnitude of the 25-year (1980-2005)
trend in average minimum temperature of the breeding season (TMIN), magnitude of the 25-year
(1980-2005) trend in average total precipitation of the breeding season (PRECIP), edge density
(ED), percent developed land (DEVEL), interaction between the magnitude of the 25-year
(1980-2005) trend in average maximum temperature of the breeding season and edge density
(TMAX*ED), interaction between the magnitude of the 25-year (1980-2005) trend in average
minimum temperature of the breeding season and edge density (TMIN*ED), interaction between
the magnitude of the 25-year (1980-2005) trend in average total precipitation of the breeding
season and edge density (PRECIP*ED), and survey effort (EFF). Model comparison was done
using Deviance Information Criterion (DIC); pD indicates the effective number of parameters.
Grouping

All Species

Long-distance
Migrants

Short-distance
Migrants

Resident
Species

Model

ΔDIC

pD

TMAX*ED + TMIN*ED + PRECIP*ED +
DEVEL + EFF
TMAX + TMIN + PRECIP + ED + DEVEL +
EFF
TMAX + TMIN + PRECIP + EFF
TMAX*ED + TMIN*ED + PRECIP*ED + EFF
TMAX + TMIN + PRECIP + ED + EFF
TMAX*ED + TMIN*ED + PRECIP*ED +
DEVEL + EFF
TMAX + TMIN + PRECIP + ED + DEVEL +
EFF
TMAX*ED + TMIN*ED + PRECIP*ED + EFF
TMAX + TMIN + PRECIP + ED + EFF
TMAX + TMIN + PRECIP + EFF
TMAX + TMIN + PRECIP + EFF
TMAX + TMIN + PRECIP + ED + DEVEL +
EFF
TMAX*ED + TMIN*ED + PRECIP*ED +
DEVEL + EFF
TMAX + TMIN + PRECIP + ED + EFF
TMAX*ED + TMIN*ED + PRECIP*ED + EFF
TMAX + TMIN + PRECIP + ED + DEVEL
+EFF
TMAX + TMIN + PRECIP + ED + EFF
TMAX*ED + TMIN*ED + PRECIP*ED + EFF
TMAX*ED + TMIN*ED + PRECIP*ED +
DEVEL + EFF
TMAX + TMIN + PRECIP + EFF

0.0

117.0

1.25

114.1

7.5
8.5
11.5
0.0

111.6
116.7
114.4
112.0

0.19

109.0

10.5
14.4
22.1
0.0
1.3

110.8
108.4
106.9
92.3
95.9

4.2

97.8

5.3
8.1
0.0

95.8
98.2
104.8

0.5
1.3
5.2

104.9
108.0
109.8

26.6

103.7

84

3.3.2 Responses of different migratory groupings
3.3.2.1 Temporal turnover
The mean temporal turnover in long-distance migratory communities was 0.381 (LCI 0.377, UCI
0.384). Short-distance migrants displayed lower average temporal turnover than long-distance
migrants, with the mean TURN value of 0.349 (LCI 0.345, UCI 0.353). Resident birds
underwent higher temporal turnover than other migratory groupings, with the mean value of
TURN of 0.423 (LCI 0.419, UCI 0.427). For all three migratory groupings, the Adirondack
Mountains underwent the highest temporal turnover across BBA1980 and BBA2000 (Fig. 3.1).
Additionally, resident birds also displayed high temporal turnover in the Catskills Mountains, as
well as in the western parts of New York. This high-turnover pattern, however, was not found for
long- or short-distance migratory birds in these regions (Fig. 3.1). Rather, migratory birds
showed higher temporal turnover also in more urbanized areas surrounding New York City;
temporal turnover in resident species was low in this region. The remaining regions of the state
generally underwent the lower temporal turnover (Fig. 3.1).
The most complex model was selected as the top model for long-distance migrants and
resident species (Table 3.1). For resident birds, the second best model (i.e., the model with
interaction terms but without main effects of DEVEL) can also be considered viable given the
values of ΔDIC (Table 3.1). The selection of the most complex model indicates that interactions
between landscape fragmentation and some aspects of climatic change are important to
explaining the temporal turnover of in communities of long-distance migrant and resident birds.
For short-distance migrants, the top selected model was the one with main effects of all
covariates, but lacking the interaction terms (Table 3.1), suggesting that interactions between
land cover and climate do not contribute importantly to explaining the changes in the

85

communities of short-distance migrants. Model validation indicated that temporal turnover was
best predicted for long-distance migrants and resident birds (Appendix 3.1 for model validation
plots), while the top model for short-distance migratory birds performed relatively poorly.
Based on the results from the top model for long-distance migrants, temporal turnover in
long-distance migrant communities was positively associated with TMIN (Fig. 3.2), but the
effect of minimum temperatures on TURN was affected by the value of ED (Fig. 3.3). Similarly
to results for all species, we found that in regions of low landscape fragmentation, large increases
in minimum temperatures were associated with high temporal turnover, but the effect of TMIN
diminished or reversed in locations with high landscape fragmentation. Temporal turnover was
also related to the interaction between PRECIP and ED (Fig. 3.2). In regions of low landscape
fragmentation, increasing precipitation is associated with decreasing temporal turnover; on the
other hand, in regions of high landscape fragmentation, temporal turnover increases with
increasing precipitation. The effect size of the PRECIP-ED interaction on temporal turnover,
however, is lower than that of TMIN-ED interaction (Fig. 3.3). There was a positive relationship
between temporal turnover and DEVEL (Fig. 3.2), indicating that high temporal turnover is
associated with highly urbanized regions.
Parameter estimates of the top model for short-distance migrants indicated that temporal
turnover in short-distance migrant communities was positively related to TMAX. Thus, higher
temporal turnover was observed in regions where maximum temperatures increased the most
(Fig. 3.2). The relationship between temporal turnover in short-distance migratory communities
and ED was negative, while the relationship between temporal turnover and DEVEL was
positive. These parameter estimates suggest that low temporal turnover was generally found in

86

regions with high landscape fragmentation, though highly urbanized regions generally
experienced increased temporal turnover.
Parameter estimates of the top model for resident species indicated that the effect of both
maximum and minimum temperatures on TURN in resident birds was affected by the value of
ED (Fig. 3.2, Fig. 3.3). In regions of low landscape fragmentation, large increases in maximum
and minimum temperatures were associated with high temporal turnover, but the effect of TMIN
diminished or reversed in locations with high landscape fragmentation. Temporal turnover was
also negatively associated with PRECIP (Fig. 3.2), indicating that higher temporal turnover was
observed in regions where precipitation decreased the most. This relationship between temporal
turnover and PRECIP was not affected by the level of landscape fragmentation. There was a
negative association of temporal turnover and ED, but the presence of multiple ED interaction
terms in the model makes it difficult to interpret the main effects of ED on temporal turnover.
TURN in resident birds communities was negatively related to DEVEL, indicating that low
temporal turnover in resident birds was generally found in highly developed regions.
3.3.2.2 Extinction
Average extinction in long-distance migrants equaled 0.226 (LCI 0.221, UCI 0.230). Shortdistance migrants displayed slightly lower average extinction than long-distance migrants, with
the mean values of EXT of 0.187 (LCI 0.183, UCI 0.190). The mean EXT for resident birds was
0.203 (LCI 0.198, UCI 0.207). For all three migratory groups, higher than average extinction was
displayed in the north-eastern part of the state, approximately in the location of the Adirondack
Mountains (Fig. 3.1). Long- and short-distance migrants also underwent higher extinction in the
Finger Lakes region (Fig. 3.1).

87

The most complex model was selected as the top model for all three migratory groupings
(Table 3.2). Selection of the most complex model suggests that interactions between landscape
fragmentation and trends in climatic factors contribute to explaining extinction rates observed in
disparate migratory groupings. Similarly to models of temporal turnover, model validation
indicated that extinction was best predicted for long-distance migrants and resident birds (see
Appendix 3.1 for model validation plots), while the top model for short-distance migratory birds
performed relatively poorly.
Based on the results from the top model for long-distance migrants, extinction was
negatively associated with TMAX (Fig. 3.2), indicating that higher extinction in long-distance
migrant communities was observed in regions where maximum temperatures decreased the most.
The relationship between extinction and TMIN was affected by ED (Fig. 3.4). In regions of low
landscape fragmentation, large increases in minimum temperatures were associated with high
proportion of extinction events, but the effect of TMIN diminished in regions with pervasive
landscape fragmentation (Fig. 3.4). PRECIP was positively associated with extinction (Fig. 3.2),
though this effect was also affected by the level of landscape fragmentation (Fig. 3.4). In regions
of low landscape fragmentation, increasing precipitation was associated with slowly decreasing
extinction rates; while in highly-fragmented regions, extinction increased sharply with increasing
precipitation (Fig. 3.4). There was a positive association between extinction and DEVEL (Fig.
3.2), suggesting increased extinction in long-distance migrants in highly urbanized regions.
Parameter estimates of the top model for short-distance migrants suggested the
importance of the interaction of TMIN and ED (Fig. 3.2). In regions of low landscape
fragmentation, large increases in minimum temperatures were associated with high proportion of
extinction events in short-distance migratory birds, but the effect of TMIN reversed and

88

diminished in highly fragmented regions (Fig. 3.4). Extinction was also positively related to
DEVEL (Fig. 3.2), suggesting that high rates of extinction were expected in highly urbanized
regions.
Based on the results from the top model for resident birds, extinction was positively
associated with TMAX (Fig. 3.2), indicating that extinction increased with increasing maximum
temperatures. Extinction was related to the interaction of TMIN and ED (Fig. 3.2). Specifically,
in regions of low landscape fragmentation, large increases in minimum temperatures were
associated with high proportion of extinction events, but the effect of TMIN reversed and
diminished in highly fragmented locations (Fig. 3.4). Extinction was also negatively related to
DEVEL, indicating that highly developed regions experienced lower rates of extinction in
resident birds.
3.3.2.3 Colonization
Mean colonization for long-distance migratory birds was 0.192 (LCI 0.188, UCI 0.195), while
short-distance migratory birds showed slightly lower rates of colonization, with a mean COL
value of 0.155 (LCI 0.152, UCI 0.159). In comparison, resident birds displayed the highest
colonization rates, with the mean value of COL of 0.268 (LCI 0.264, UCI 0.272). Resident
species underwent high colonization throughout most of the state, with exception of central and
southern parts of New York (Fig. 3.1). Long- and short-distance migrants showed relatively low
colonization rates throughout New York, with exception of several isolated locations in the
Adirondack Mountains (Fig. 3.1).
More than one model was selected as the top model for all three avian groupings (Table
3.3). For long-distance migrants, the most complex model along with the model with main
effects of all covariates but without the interaction terms were selected as top models (Table 3.3).

89

For short-distance migrants, the simplest model along with the model with main effects of all
covariates but without the interaction terms were selected as the top models (Table 3.3). For
resident birds, the top three models included the model with main effects of all covariates but
without the interaction terms (the best model), the model with main effects excluding the
DEVEL covariate (second best model), and the model with interaction terms excluding the main
effects of DEVEL (third best model). Model validation indicated that colonization was best
predicted for resident birds (see Appendix 3.1 for model validation plots), while the top models
for long- and short-distance migratory birds performed relatively poorly.
Given that none of the interactions terms from the top model was associated with
colonization in long-distance migrants, we make inferences based on the model without the
interaction terms. Proportion of colonization events in long-distance migrant communities was
positively associated with TMAX and negatively associated with PRECIP (Fig. 3.2), indicating
that higher colonization was observed in regions where maximum temperatures increased the
most and where precipitation decreased the most. There was a positive relationship between
colonization and both ED and DEVEL (Fig. 3.2), suggesting higher colonization rates in regions
with high landscape fragmentation and those highly urbanized.
For short-distance migrants, we did not find any relationships between proportion of
colonization events and any of the tested covariates in the best model (Fig. 3.2). Based on the
results from the top model for resident birds, colonization was negatively associated with
PRECIP (Fig. 3.2), indicating higher colonization in regions where precipitation decreased the
most. Colonization in resident bird communities was also negatively related to ED (Fig. 3.2),
indicating that low colonization in resident birds was generally found in regions with pervasive
landscape fragmentation.

90

3.4 DISCUSSION
My findings supported my first hypothesis that changes in avian community composition are
related to changes in climate, but the strength of the relationship between temporal turnover and
climate change is affected by the level of landscape fragmentation. Associations were weaker
between community change and climatic change in regions with prevalent habitat fragmentation.
However, my findings did not support my second hypothesis regarding the differences in
responses of long-distance migrants and other migratory groupings to climate change.
Associations of temporal turnover and climate change were generally similar among all
migratory groupings.
3.4.1 Relationships with climate change and habitat fragmentation
There are at least two alternative interpretations for the pattern of associations between
community dynamics, climatic change, and landscape fragmentation. One, communities in
contiguous landscapes might be able to respond more quickly to changes in temperatures and
precipitation, perhaps due to potential for unimpeded movement. Such interpretation would fit
well with research suggesting that unfragmented habitats allow higher rates of species
distributional shifts associated with climate change (e.g., Opdam and Washer 2004). If that were
the case, communities in fragmented habitats would be likely to lag behind changing climatic
conditions and, therefore, be less capable of adapting in the face of climate change. This lag
would initially be reflected in communities being compositionally similar across time and thus
showing no temporal turnover. However, through time, it is likely that species will eventually
disappear from locations in which they can no longer find suitable climatic conditions. Such
localized species extinctions will ultimately result in higher community temporal turnover and

91

higher extinction rates, perhaps on par with those currently observed in unfragmented
landscapes.
How plausible is the unimpeded movement explanation in New York State? Using the
same dataset, Zuckerberg et al. (2009) documented shifts in breeding ranges of 129 species of
birds. They found that the birds exhibiting shifts in their range boundaries all had different
breeding habitat associations, suggesting that land cover was not a factor in the observed range
expansion. Specifically, they found no evidence for more pronounced range shifts of forest
breeding birds in unfragmented habitats such as those of the Adirondack and Catskill Mountains.
Given their findings, it seems unlikely that contiguous habitats facilitating the climate changedriven range shifts may be at the root of the higher rates of temporal turnover, extinction, and
colonization I am observing.
An alternative interpretation is that avian communities in fragmented landscapes are
more robust to changes in temperatures and precipitation resulting from climate change than
communities found in contiguous habitats. If that were the case, we would expect little or no
temporal change in community composition and weaker associations with changing climatic
conditions in fragmented landscapes. Indeed, heterogeneous and fragmented habitats are often
composed of a large proportion of habitat generalists (e.g., Tscharntke et al. 2012, Estavillo et al.
2013). Habitat generalists tend to have wider thermal breadths than habitat specialists
(Barnagaud et al. 2012), which potentially allows them to tolerate greater changes in climatic
conditions. On the other hand, narrower thermal niches of habitat specialists (Barnagaud et al.
2012) might make them vulnerable to even small increases in temperatures. Thus, in this
scenario, we would expect communities in contiguous or less fragmented habitats (i.e.,
comprised mainly of habitat specialists) to continue to undergo relatively large compositional

92

changes as a result of climate change, while those found in fragmented regions (i.e., comprised
largely of generalists) to show relatively small changes, until perhaps some thermal threshold is
reached.
A post hoc analysis revealed that temporal turnover in generalist birds was lower than
that of birds with any other habitat associations (mean temporal turnover of 0.238 [LCI 0.234,
UCI 0.242] for generalist birds in comparison with 0.384 [LCI 0.381, UCI 0.387] for all breeding
birds, 0.444 [LCI 0.440, UCI 0.448] for forest, 0.317 [LCI 0.131, UCI 0.321] for shrubland,
0.439 [LCI 0.432, UCI 0.446] for grassland, and 0.536 [LCI 0.532, UCI 0.541] for wetland
breeding birds). Rates of extinction and colonization events were also lower for generalists than
for any other habitat grouping (extinction: 0.122 [LCI 0.119, UCI 0.126; generalists], 0.210 [LCI
0.206, UCI 0.213; all breeding birds], 0.233 [LCI 0.228, UCI 0.237; forest], 0.181 [LCI 0.177,
UCI 0.186; shrubland], 0.331 [LCI 0.323, UCI 0.339; grassland], 0.307 [LCI 0.301, UCI 0.313;
wetland]; colonization: 0.132 [LCI 0.128, UCI 0.135; generalists], 0.214 [LCI 0.211, UCI 0.218;
all breeding birds], 0.265 [LCI 0.261, UCI 0.269; forest], 0.159 [LCI 0.156, UCI 0.163;
shrubland], 0.148 [LCI 0.143, UCI 0.152; grassland], 0.309 [LCI 0.304, UCI 0.314; wetland]).
Lower values of temporal turnover, extinction, and colonization in these communities perhaps
suggest that generalist birds are indeed more robust to climatic changes, lending support to the
second interpretation of the pattern of stronger associations between community change and
climatic change in unfragmented regions.
3.4.2 Responses of different migratory groupings
I found differences in the rates of temporal turnover in communities of resident, long-, and shortdistance migratory birds. For example, communities of long-distance migrants showed higher
proportion of extinction events than the other two avian groupings. Resident birds displayed

93

higher temporal turnover and higher proportion of colonization events than long- and shortdistance migratory birds. Despite these differences in community change, I found little evidence
to support our second prediction that temporal turnover in communities of long-distance
migrants would show weaker associations with changing climatic conditions. In fact, I found that
the strength of the associations of temporal turnover with climatic changes was similar across all
avian groupings. Thus, it is unlikely that the disparities in the observed rates of community
change in this ecological system can be explained by the variation in responses of different
migratory groupings to changing climatic conditions.
Similarly, I found few differences in responses of different migratory groupings to
landscape fragmentation, though strong interaction between climatic and habitat fragmentation
variables might have obscured a potential relationship. Interestingly, however, I found that
increasing amount of urban development was associated with high temporal turnover in longand short-distance migratory communities, but resident species underwent lower temporal
turnover in highly urbanized regions. It is likely that highly developed and residential areas
might be beneficial to some resident birds, especially during the wintering months when
supplemental feeding is provided. Thus, it is plausible that the observed higher proportion of
colonization events in resident birds is driven by factors other than climatic changes.
Alternatively, higher proportion of colonization events in resident species can be
associated with less benign climatic conditions during winter months rather than climatic
changes that occurred during the breeding season. Winter temperatures in higher latitudes have
shown larger increases than temperatures during spring or summer months and are expected to
continue to rise more drastically (IPCC 2013). In contrast, temperatures in lower latitudes have
not increased as drastically (IPCC 2013); thus climatic conditions for migratory birds on their

94

wintering grounds remained relatively stable. Milder winters in higher latitudes combined with
increases in supplemental feeding would contribute to higher survival of resident birds, thus
potentially increasing chances of colonization of new sites by these birds. It is likely then that
future avian communities will become more homogenous by increasing in proportion of resident
species and decreasing in proportion of migratory birds.
3.4.3 Conclusions
My study provides new insights into the drivers of temporal turnover in avian communities. I
showed that while changes in avian community composition are closely related to changes in
climatic conditions, this association is strongly influenced by the level of landscape
fragmentation. Specifically, highly fragmented regions are likely to undergo smaller changes in
composition and structure as a result of climate change, whereas those changes are likely to be
more pronounced in contiguous and unfragmented habitats. I suggest that relative contributions
of habitat generalists and specialists in ecological communities might be one of the leading
factors explaining these divergent patterns.
My findings have significant implications for conservation theory and practice. Current
conservation strategies commonly focus on conserving large regions with undisturbed and
contiguous habitats that generally support high species diversity. My research suggests that
faunal communities of such habitats will undergo the most drastic compositional changes as a
result of climate change. Thus, my research stresses the importance of putting the existing
conservation strategies in the context of climate change.
3.5 ACKNOWLEDGMENTS
I would like to thank the thousands of volunteers who participated in both New York State
Breeding Bird Atlases. I also thank Kimberley Corwin and Kevin McGowan for supplying atlas
95

databases and other information. The manuscript benefited from discussions with members of
Quantitative Wildlife Laboratory at Michigan State University. I thank Dr. Colin M. Beier,
Daniel Bishop, and John Wiley for supplying the climate dataset. This study received financial
support from NASA and Boone and Crockett Club.

96

SCALE-DEPENDENCE OF TEMPORAL CHANGES IN AVIAN COMMUNITIES

In prep for publication:
Jarzyna, M.A. et al. (the co-authors to be determined). Temporal changes in avian communities
and processes driving these changes are scale-dependent. Potential outlets: Global Change
Biology, Global Ecology and Biogeography, Nature Climate Change, Journal of
Biogeography.
ABSTRACT
It has long been acknowledged that key biodiversity patterns vary with spatial scale and the
mechanisms driving these patterns are inherently scale-dependent. Studies investigating scaling
of biodiversity and the mechanisms pertinent at each spatial scale have mostly focused on one
time step. My goal was to evaluate spatial scaling patterns of avian assemblages through time
and to identify environmental drivers of the change. I used temporal turnover, proportion of
extinction events, and proportion of colonization events as measures of temporal change in avian
communities at 5 spatial resolutions (5x5, 10x10, 20x20, 40x40, and 80x80 km). I applied
Bayesian spatially-varying intercept models to evaluate the relationships among community
change and change in climatic conditions, landscape fragmentation, and elevation at each of the
spatial scales. I found that temporal turnover, extinction, and colonization declined with
increasing spatial scale, but this decline was steeper for extinction than for colonization.
Specifically, mean temporal turnover decreased from 0.384 at 5x5 km to 0.144 at 80x80 km
spatial scale. Mean extinction decreased from 0.210 at 5x5 km to 0.036 at 80x80km scale. Mean
colonization decreased from 0.214 at 5x5 km to 0.110 at 80x80 km scale. I also found that the
influence of different environmental drivers to changes in avian communities was scaledependent. Specifically, climate change covariates were important at smaller scales (i.e., 5x5 and
10x10 km), while landscape characteristics were relevant across coarser scales. My findings

97

suggest that climate change will initially affect communities at fine spatial scales, though the
magnitude of these changes is likely to be facilitated by landscape fragmentation at larger scales.
4.1 INTRODUCTION
Patterns of biological diversity through space and time are an indication of the underlying
mechanisms that control biodiversity and as such are of key interest to ecology. It has long been
acknowledged that key biodiversity patterns vary with spatial scale of observation (e.g., Storch et
al. 2004, Keil et al. 2011). This scale-dependence is perhaps best demonstrated in patterns of
species richness (e.g., Arrhenius 1921, Rahbek 2005), but other attributes of biodiversity have
also been shown to be contingent upon the spatial scale of investigation. For example, spatial
turnover in species composition (i.e., an aggregate of gains and losses of species comprising a
community) has been demonstrated to have scaling properties (Gaston et al. 2007), as have
patterns of extinction risk (Hartley and Kunin 2003), immigration and emigration (Englund and
Hamback 2007), colonization and biological invasions (Menendez and Thomas 2000, Davies et
al. 2005, Powell et al. 2013). Given such pervasive scale-dependence of biodiversity patterns, it
is to be expected that mechanisms controlling these patterns are also conditional on the scale of
investigation (e.g., Ricklefs 1986, Currie 2004, Storch et al. 2004). For example, many attributes
of biological communities at the continental scales are thought to be driven primarily by
evolutionary factors (Ricklefs 2004, Keil and Jetz 2014), whereas climate variability and
regional land cover dynamics are pertinent across spatial grains of tens or hundreds of kilometers
(Willis and Whittaker 2002, Field et al. 2009). At yet finer grains, specific habitat characteristics
and intra- and inter-specific interactions are thought to become the primary drivers of
biodiversity (Wiens 1989, Tilman 2004, Belmaker and Jetz 2011). However, while research to
date has provided insight into scale-dependence of biodiversity patterns, the majority of studies

98

were limited to one time step and few have looked at spatial scaling of temporal biodiversity
change.
Measures of temporal changes in biodiversity have been commonly used as indicators of
environmental change (e.g., Calvarheiro et al. 2013, Chapter 3 of this volume); thus, it is
important that we understand how these measures depend on the spatial scale. Furthermore,
given the widespread use of biodiversity as a measure of ecosystem response to global change, a
thorough investigation of mechanisms relevant to temporal biodiversity dynamics across spatial
scales is perhaps more important than ever. A few examples of scaling of temporal changes in
biodiversity do exist. For example, Keil et al. (2011) found that temporal change in species
richness decreased with increasing spatial grain, while Calvarheiro et al. (2013) showed spatial
turnover in species composition to be decreasing across time, indicating ongoing homogenization
of biological communities. Despite these initial investigations into the scaling properties of
temporal biodiversity change, limitations to the understanding of this phenomenon persist. First,
studies that explicitly look at scale-dependence of temporal biodiversity change generally focus
on species richness as a measure of biodiversity. However, consequences of biodiversity change
will extend much beyond changes in species richness (Barbet-Massin and Jetz 2014). Thus, it is
important to examine measures that explicitly take species identities into account, such as
temporal turnover and the underlying processes contributing to temporal turnover, i.e., extinction
and colonization. Second, to the best of my knowledge, none of the studies has yet investigated
the relative importance of environmental factors influencing temporal changes in biodiversity at
different spatial scales. Indeed, climate and land-cover change are currently considered the two
biggest threats to biodiversity (e.g., de Chazal and Rounsevell 2009), but the spatial scales at
which these two factors are relevant to biodiversity are not yet fully known.

99

Here, I used data from the New York State Breeding Bird Atlas to evaluate spatial scaling
patterns of temporal changes in avian assemblages. I define temporal change through measures
of temporal turnover, proportion of extinction events, and proportion of colonization events.
Though spatial scaling often refers to either the resolution (i.e., grain) of the investigation or the
geographic extent of the study, I focus on the resolution only and refer to it as spatial scale or
spatial grain interchangeably throughout the paper. I sought to evaluate two interrelated
questions: (1) how do patterns of temporal changes in avian assemblages vary across spatial
scales? and (2) at what spatial grains do different environmental factors influence temporal
turnover, extinction, and colonization?
4.2 METHODS
4.2.1 Breeding Bird Atlas
The New York State Breeding Bird Atlas (BBA) is a statewide survey that documented the
distribution of breeding birds in New York. To date, BBA has been conducted in two time
periods, 1980-85 (hereafter, BBA1980; Andrle and Carroll 1988) and 2000-05 (hereafter,
BBA2000; McGowan and Corwin 2008). For both BBAs, a grid system was used to define the
basic unit for reporting data; the BBA reporting unit (a block) measures 5x5 km and a total of
5,335 blocks covered the entirety of New York State. This data set represents one of the largest
and finest resolution atlases in the world (Gibbons et al. 2007).
A total of 242 species were recorded for BBA1980 and 248 species were recorded for
BBA2000 (Andrle and Carroll 1988, McGowan and Corwin 2008, see Appendix 1.1).
Observations were made by skilled birders who spent at least 8 hours in each block, visited all
cover types in each block, and included at least one nighttime visit to document nocturnal
species. Observer effort was recorded for each BBA and reported as a number of person hours
100

(measured as the sum of the number of hours spent in each block x the number of people
surveying each block; McGowan and Zuckerberg 2008). The BBA represents a detection/nondetection dataset and non-detection indicates that species could not be found given search criteria
(McGowan and Corwin 2008).
4.2.2 Spatial scaling of community change
To quantify temporal change in avian assemblages, I calculated temporal turnover (TURN) as an
aggregate of species gains (hereafter called colonization) and species losses (hereafter called
extinction):
𝑇𝑈𝑅𝑁 =

𝐸+𝐶

𝐸+𝐶+𝑃

,

Eqn. 4.1

where E is the number of species that went extinct in the block between BBA1980 and

BBA2000, C is the number of species that colonized the block between BBA1980 and
BBA2000, and P is number of species in the same block common to both BBAs (i.e.,
persistence).
In addition to temporal turnover, I quantified its components: proportion of extinction
events (i.e., proportion of species lost at a particular site between BBA1980 and BBA2000;
EXT) and proportion of colonization events (i.e., proportion of species gained at a particular site
between BBA1980 and BBA2000; COL) in each block between BBA1980 and BBA2000. Here,
I use terms extinction and colonization to refer to the gain or loss of a species from a defined
point in space (in this case, a spatial grain of 5x5 km, 10x10 km, 20x20 km, 40x40 km, or 80x80
km; Gaston and Blackburn 2002). EXT and COL were calculated as follows:
𝐸𝑋𝑇 =
𝐶𝑂𝐿 =

𝐸

𝐸+𝑃
𝐶

,

𝐸+𝐶+𝑃

Eqn. 4.2
,

Eqn. 4.3

101

For TURN and COL, we used the sum of E, C, and P as the denominator of the Eqns. 1
and 2. We recognize that species already present in a block in BBA1980 could not have
colonized that block and that species other than those included in E, C, and P metrics could have
potentially colonized the site. However, we were simply interested in the proportion of species
that colonized the site relative to all the species found within the block; thus we deem the sum of
E, C, and P as an appropriate choice for the denominator of both equations. For EXT, we used
the sum of E and P as the denominator, because only species present in a block in BBA1980
could have potentially gone extinct. The values of TURN, EXT, and COL are bounded by 0 and
1; values approaching 0 indicate low temporal turnover, extinction, or colonization while values
approaching 1 indicate high temporal turnover, extinction, or colonization in a block between
BBA1980 and BBA2000. To account for unequal sampling effort, we followed Chao et al.
(2005) and used individual-based probabilistic approach to calculate E, C, and P.
We calculated TURN, EXT, and COL at the following spatial grains: 5x5 km, 10x10 km,
20x20 km, 40x40 km, and 80x 80 km. To estimate the mean of each metric and construct the
95% confidence intervals, we applied a non-parametric bootstrap procedure (10,000
resamplings) at each of the spatial grains.
4.2.3 Factors influencing community change across spatial grains
We included temporal trends in climatic conditions across the 1980-2005 time period and
landscape fragmentation as potential correlates of temporal turnover, extinction, and
colonization. We included elevation in all the models to account for topographical variation of
New York State. To quantify trends in climatic conditions, we used the PRISM (Parameterelevation Regressions on Independent Slopes Model) climate mapping system (Daly and Gibson
2002). I calculated the magnitude of the 25-year (1980-2005) trend in average maximum and

102

minimum temperatures of the breeding season (TMAX and TMIN, respectively; expressed in
°C/25 years) and in total monthly precipitation of the breeding season (PRECIP; expressed in
mm/25 years) using Ordinary Least Squares regression.
I used the National Land Cover Data (NLCD) derived from the Landsat Thematic
Mapper satellite data to quantify landscape fragmentation. There is no NLCD available for the
time period of BBA1980; thus, I used a NLCD2001 and space-for-time substitution to assess
whether temporal turnover, extinction, and colonization are related to landscape fragmentation
(for more information on space-for-time substitution see Pickett 1989). Because I focused our
analysis on a diverse suite of species with varying habitat requirements, I chose a landscapescale variable to capture broad-scale variation in habitat fragmentation. Edge Density (ED) is a
good measure of landscape fragmentation because an increase in habitat edge is a primary
outcome of habitat fragmentation (Hargis et al. 1998). However, in situations when landscape
consists entirely of one cover type, the ED would be 0 regardless of the type of land-cover
present. Therefore, to differentiate between landscapes that consisted entirely of natural land
cover (e.g., forest) and those mostly developed (e.g., cities), I also calculated the proportion of
developed land (DEVEL). Landscape analysis was conducted using FRAGSTATS 4.1
(McGarigal et al. 2012) and the Geospatial Modelling Environment (GME,
http://www.spatialecology.com/gme/).
I used Digital Elevation Models for New York State to quantify mean elevation in each
block and for each spatial grain (ELEV). To account for potential survey effort bias, I calculated
the absolute difference in the number of person hours between BBA1980 and BBA2000
(EFF = |EFF1980 − EFF2000|) for each BBA block and included it as a covariate in my

103

models. I standardized and centered all the covariates to ease the interpretation of the parameter
estimates.
To estimate the importance of environmental factors in driving temporal turnover,
extinction, and colonization at each of the spatial grains (i.e., 5x5, 10x10, 20x20, and 40x40 km),
I built a statistical model that included the main effects of all the covariates (TMAX, TMIN,
PRECIP, ED, DEVEL, ELEV, and EFF). I did not run a model for the largest spatial grain (i.e.,
80x80 km) because of the insufficient sample size. We used a Bayesian model estimation
procedure. Each response variable (temporal turnover TURN, extinction EXT, or colonization
COL) at each location (s) was modeled as a binomial random variable. Let y(s) be a response
variable, then y(s) ~ B(ρ(s),N(s)) where N(s) is the total number of species found at least once in
a block from the first to the second atlas count (for TURN and COL) or the total number of
species found in first atlas count (for EXT). For each BBA block, let p(s) = y(s)/N(s). Also, let
x(s) be a vector containing the explanatory variables for each model. All models had spatially
varying intercepts, which were created by adding spatially structured random effects to the
model intercept. A generalized linear model approach was used to relate responses to
explanatory variables. Using the logit link η(s) = log[p(s)/{1-p(s)}] the SVI model was
𝜂(𝑠) = 𝑤(𝐬) + 𝜷𝐱(𝑠)

Eqn. 4.4

where β is a vector of regression coefficients and w(s) is a random effect representing

spatially varying adjustments to the intercept. We assume each element in w(s) arises from a
spatial Gaussian Process (GP) (see, e.g., Banerjee et al. 2004 or Cressie and Wikle 2011 for more
details). Specifically, w(s) ~ GP(0, C(s, s’)), where s and s’ are any two locations within the

study area, the spatial covariance C(s, s’) = σρ(s, s’; ϕ) with variance parameter σ2, correlation

104

function ρ(·; ϕ), and spatial decay parameter ϕ. The exponential spatial correlation was assumed
for ρ(·).
Prior distributions on the parameters complete the hierarchical model specification (e.g.,
Gelman et al. 2004). The global regression coefficients, β’s followed a Normal distribution N(0,
100), while the variance component σ2 was assigned inverse-Gamma IG(2, 1) prior. The spatial
decay parameter, ϕ, followed an informative Uniform prior, with support that ranged from 1 km
to the maximum distance between any two grid cells. I ran three MCMC chains for 20,000
iterations each. Convergence was diagnosed using the Gelman-Rubin diagnostic (Gelman and
Rubin 1992).
Prior to statistical analysis, I removed all blocks that did not have continuous land cover
coverage, blocks with more than 50% open water coverage, and those that did not have survey
effort data. Ultimately, I used a total of N5x5 = 4,742 blocks to fit the models for the 5x5 km
grain; N10x10 = 1,186 blocks to fit the models for the 10x10 km grain; N20x20 = 226 blocks to fit
the models for the 20x20 km grain; and N40x40 = 54 blocks to fit the models for the 40x40 km
grain. The 80x80 km spatial grain included only N80x80 = 7 blocks. All models were run in R
2.15.1 statistical program (R Development Core Team 2013, http://www.r-project.org/) using
package spBayes (Finley and Banerjee 2013). In a Bayesian framework, parameter estimates have
valid posterior distributions and inferences are made using mean or median as well as the credible
intervals of these posterior distributions. Credible intervals overlapping zero imply that the parameter in
question is not different from zero. Summaries of parameter estimates were generated using the R

CODA package (Plummer et al. 2012).

105

4.3 RESULTS
4.3.1 Spatial scaling of community change
Mean TURN decreased with increasing spatial grain (Fig. 4.1A). At 5x5 km spatial grain, mean
value of TURN was 0.384, indicating that on average 38% of species had either colonized or
become extinct within the average 5x5 km block and approximately 62% of the species persisted
across both time periods. At 10x10 km, mean value of TURN decreased to 0.266. At 20x20 km
grain, mean TURN further decreased to 0.201; while at 40x40 km it was 0.160. At the largest
spatial grain (i.e., 80x80 km), mean value of TURN was 0.144, which indicates that on average
only 14% of species had either colonized or become extinct within the average 80x80 km block
and approximately 86% of the species persisted over time. In terms of its spatial distribution, the
highest temporal turnover at the 5x5 km and 10x10 km spatial grains was in the north-eastern
part of the state, where the Adirondack Mountains are (Fig. 4.2). At the 20x20 km grain, the
spatial patterns of temporal turnover became more homogenous, with central New York showing
lower temporal turnover than eastern or western regions. At the 40x40 km spatial grain, most
sites showed temporal turnover between 0.1 and 0.2, with a few exceptions in the north-eastern
New York. At the 80x80 km spatial scale, temporal turnover exceeded 0.2 in only one block
(Fig. 4.2).

106

Figure 4.1 Spatial scaling of temporal turnover, extinction, and colonization observed between 1980-85 and 2000-05 for avian
communities in New York State across five different spatial grains (i.e., 5x5 km, 10x10 km, 20x20 km, 40x40 km, and 80x80 km).
Temporal turnover was calculated as a proportion of all species gained or lost between 1980-85 and 2000-05 in a particular site (i.e.,
5x5 km Breeding Bird Atlas block or a scaled-up version of it) relative to all species recorded across both time periods; extinction was
calculated as a proportion of species lost between 1980-85 and 2000-05 in a particular site relative to all species present in a site in
1980-85; colonization was calculated as a proportion of species gained between 1980-85 and 2000-05 in a particular site relative to all
species recorded across both time periods. High values of all three metrics indicate high temporal turnover, extinction, and
colonization.

107

Figure 4.2 Patterns of temporal turnover observed between 1980-85 and 2000-05 for avian communities in New York State across
five different spatial grains (i.e., 5x5 km, 10x10 km, 20x20 km, 40x40 km, and 80x80 km). Temporal turnover was calculated as a
proportion of all species gained or lost between 1980-85 and 2000-05 in a particular site (i.e., 5x5 km Breeding Bird Atlas block or a
scaled-up version of it) relative to all species recorded across both time periods. High values indicate high temporal turnover.

108

I detected a decrease in mean EXT with increasing spatial grain (Fig. 4.1B). At 5x5 km
spatial grain, mean value of EXT was 0.210, indicating that on average 21% of species had
become extinct within the average 5x5 km block between 1980 and 2005. At 10x10 km spatial
scale, mean value of EXT declined to 0.132. At 20x20 km grain, mean EXT was 0.096; while at
40x40 km, it decreased further and equaled 0.071. At the largest spatial grain (i.e., 80x80 km),
mean value of EXT was 0.036, which indicates that on average only 4% of species had become
extinct within the average 80x80 km block. Similarly to TURN, the highest proportion of
extinction events at the 5x5 km and 10x10 km spatial grains was found in the location of the
Adirondack Mountains (Fig. 4.3). At the 20x20 km and 40x40 spatial grains, EXT also was
higher in the north-east. At the 80x80 km grain, extinction became more spatially homogeneous
and never exceeded 0.1 (Fig. 4.3).
In general, mean COL also decreased with increasing spatial grain (Fig. 4.1C). At 5x5
km, mean COL took a value of 0.214, indicating that approximately 21% of all species found
across both time periods in a 5x5 km block were detected only during the second BBA. At the
10x10 km spatial scale, mean COL was 0.152. Starting at the 20x20 km grain, the decline in the
mean COL values tapered off and the values stabilized around 0.1. Mean COL was 0.116 at
20x20 km, 0.094 at 40x40 km, and 0.110 at 80x80 km spatial grain. High proportion of
colonization events for the 5x5 km and 10x10 km spatial grains were found in the north-eastern
and western parts of the state (Fig. 4.4). However, unlike in the case of temporal turnover and
extinction, COL values were less spatially clumped and more dispersed (Fig. 4.4). At 20x20 km,
colonization rates were homogeneously distributed throughout the state; at the 40x40 km, lower
colonization values were found in north-eastern (i.e., Adirondack Mountains) and south-central

109

(i.e., Allegany Plateau) New York. At the 80x80 km grain, colonization rates were generally low,
with the exception of the most northerly block (Fig. 4.4).

110

Figure 4.3 Patterns of extinction observed between 1980-85 and 2000-05 for avian communities in New York State across five
different spatial grains (i.e., 5x5 km, 10x10 km, 20x20 km, 40x40 km, and 80x80 km). Extinction was calculated as a proportion of
species lost between 1980-85 and 2000-05 in a particular site (i.e., 5x5 km Breeding Bird Atlas block or a scaled-up version of it)
relative to all species present in a site in 1980-85. High values indicate high proportion of extinction events.

111

Figure 4.4 Patterns of colonization observed between 1980-85 and 2000-05 for avian communities in New York State across five
different spatial grains (i.e., 5x5 km, 10x10 km, 20x20 km, 40x40 km, and 80x80 km). Colonization was calculated as a proportion of
species gained between 1980-85 and 2000-05 in a particular site (i.e., 5x5 km Breeding Bird Atlas block or a scaled-up version of it)
relative to all species recorded across both time periods. High values indicate high proportion of colonization events.

112

4.3.2 Factors influencing community change across spatial grains
4.3.2.1 Temporal turnover
At the 5x5 km spatial grain, TMAX was positively associated with TURN, indicating that
increasing maximum temperatures were associated with high temporal turnover (Fig. 4.5). There
was a negative association of temporal turnover and ED, suggesting that temporal turnover was
higher in less fragmented regions. There was a positive relationship between temporal turnover
and DEVEL as well as temporal turnover and ELEV (Fig. 4.5), indicating that highly urbanized
regions and locations in higher elevations generally experienced increased temporal turnover.
At the 10x10 km spatial grain, TMAX was positively related to TURN (Fig. 4.5). I also
found a negative relationship between temporal turnover and ED and a positive relationship
between temporal turnover and DEVEL (Fig. 4.5). I detected no relationship between TURN and
ELEV at this spatial grain. At the 20x20 km spatial grain, I found no relationships between
TURN and any of the climate change covariates. I found a negative relationship between TURN
and ED and a positive relationship between TURN and DEVEL (Fig. 4.5). I also detected a
negative relationship between TURN and ELEV, suggesting that temporal turnover was higher
in lower elevations. At 40x40 km grain, there I found no relationships between temporal turnover
and any of the covariates (Fig. 4.5).
4.3.2.2 Extinction
At the 5x5 km spatial grain, I found positive relationships between EXT and TMIN (Fig. 4.5),
indicating that increasing minimum temperatures were associated with high proportion of
extinction events. Extinction was also positively related to PRECIP (Fig. 4.5), suggesting higher
extinction in regions where precipitation increased the most. I found a negative association of

113

extinction and ED and a positive association of extinction and DEVEL, indicating that extinction
was highest in contiguous landscapes and in highly urbanized regions. There was a relationship
between extinction and ELEV (Fig. 4.5).
At the 10x10 km spatial grain, TMIN was positively associated with EXT, but I no longer
found relationships between extinction and TMAX or PRECIP. Extinction was negatively
associated with ED and positively related to DEVEL. I found no relationship between extinction
and ELEV at this spatial grain. At the 20x20 km spatial grain, none of the climate change
variables was associated with EXT. I found a negative relationship between EXT and ED (Fig.
4.5) and a positive relationship between EXT and DEVEL. I found no relationship between
extinction and ELEV at this spatial grain. At the 40x40 km spatial grain, I found no relationships
between extinction and any of the tested covariates (Fig. 4.5).
4.3.2.3 Colonization
At the 5x5 km spatial grain, colonization was positively related to TMAX, indicating that regions
where temperatures increased the most experienced highest proportion of colonization events
(Fig. 4.5). Colonization was also negatively related to PRECIP (Fig. 4.5), indicating that regions
with decreasing precipitation on average experienced higher proportion of colonization events.
COL was positively associated with DEVEL (Fig. 4.5), indicating that colonization was
generally higher in highly urbanized regions.
At the 10x10 km spatial grain, I found a positive relationship between colonization
and DEVEL. However, I did not find associations between COL and climatic change variables,
ED, or ELEV at this spatial scale. At the 20x20 km spatial grain, I found a positive relationship
between colonization and TMAX and a negative relationship between colonization and TMIN
(Fig. 4.5). I also found a positive association of colonization and ED, indicating that highest

114

colonization occurred in regions of highest fragmentation. I found no relationships between
colonization and DEVEL or colonization and ELEV at 20x20 km. At the 40x40 km spatial grain,
I found no relationships between colonization and any of the tested covariates (Fig. 4.5).

115

Figure 4.5 Coefficient estimates (i.e., mean values of the posterior distribution and associated credible intervals) resulting from the
models for (A) temporal turnover, (B) extinction, and (C) colonization. Abbreviations for the covariates are as follows: magnitude of
the 25-year (1980-2005) trend in average maximum temperature of the breeding season (TMAX), magnitude of the 25-year (19802005) trend in average minimum temperature of the breeding season (TMIN), magnitude of the 25-year (1980-2005) trend in average
total precipitation of the breeding season (PRECIP), edge density (ED), percent developed land (DEVEL), and elevation (ELEV).

116

4.4 DISCUSSION
I sought to address the issue of scale-dependence of temporal changes in avian assemblages by
evaluating temporal turnover, extinction, and colonization across different spatial grains and
investigating potential environmental drivers of these changes in assemblages at each of the
scales. I found that temporal turnover, extinction, and colonization declined as the spatial grain
increased, indicating that community change is smaller at coarse scales. This decline with
increasing spatial grain, however, was faster for temporal turnover and extinction than it was for
colonization. I also found that the relevance of different environmental drivers to changes in
avian communities was scale-dependent. Specifically, I found that climate change covariates
were pertinent at smaller spatial grains than landscape characteristics.
4.4.1 Spatial scaling of community change
I found that temporal changes in community composition decreased with increasing spatial scale
of investigation. This result is not surprising given that increasing spatial grains will encompass
larger areas and higher number of species, thus reducing both the likelihood of individual species
extinction or colonization (e.g., Gaston et al. 2004) and the proportion of extinction and
colonization events. However, while this decline was relatively steady for temporal turnover and
extinction, the decline in the proportion of colonization events tapered off at 20x20 km grain and
remained relatively constant across the three coarsest spatial grains (i.e., 20x20 km, 40x40 km,
and 80x80 km). While there is little direct evidence that patterns of extinction and colonization
exhibit different spatial scaling, there are reasons to suspect that it might be so. For example,
Wilson et al. (2004) suggested that processes of range retraction and expansion, and perhaps by
extension extinction and colonization, affect species distributional patterns in a disparate way.
Declining species generally leave sparse distributions because they tend to retract their ranges to
117

either optimal habitats or regions where extinction forces are not operating (Johnson 1998).
Thus, it is expected that contracting ranges would result in substantial species loss at fine spatial
grains and relatively little loss at coarse grains (Wilson et al. 2004, Keil et al. 2011). The results
from my study lend support to this expectation because I showed that there is substantial
extinction at fine grains and relatively little species loss at coarser scales.
On the other hand, species expanding their ranges tend to form more aggregated
distributions because range expansion is often dispersal-limited (Shigesada and Kawasaki 1997).
Therefore, it is expected that expanding range would involve substantial gain at fine scales as
coarse-grain areas are slowly colonized (Wilson et al. 2004). My results showing substantial
colonization at the finest grain fit well with this expectation. However, I suggest that low
colonization at coarse spatial grains proposed by Wilson et al. (2004) is characteristic of species
with relatively low mobility and low dispersal abilities, such as investigated in their study
butterflies. Birds, as highly mobile and far-dispersing organisms, would likely display a different
expansion signature, characterized by relatively high rates of colonization even at coarse grains.
Indeed, Gaston and Blackburn (2002) found that birds which exhibited higher colonization rates
in the UK tended to disperse farther during both their juvenile and adult life stages. Other
examples of a positive relationship between colonization rate and dispersal abilities abound (e.g.,
Brown and Kodric-Brown 1977, Juliano 1983, Mouquet and Loreau 2003). Consequently, I
propose that far-dispersing birds are able to colonize distant coarse-grain areas, resulting in
relatively high and constant proportion of colonization events across coarser grains explored in
this study (i.e., 20x20, 40x40, and 80x80 km).

118

4.4.2 Factors influencing community change across spatial grains
I found that the main drivers of avian biodiversity change are not only often specific to the
spatial scale of investigation, but also differ in their likely impacts among the different metrics of
community change. For example, changes in climatic conditions were relevant to both temporal
turnover and extinction only at finer spatial grains of 5x5 and 10x10 km. I detected no impacts of
climate change on these two measures of biodiversity change across coarser spatial scales. On
the other hand, colonization was associated with changes in climate both at fine grains as well as
at the intermediate spatial scales (i.e., 20x20 km). This association of colonization with climate
change at coarser grains might further explain why I detected stabilizing of colonization values at
coarser spatial grains. Perhaps changing climatic conditions in combination with dispersal
abilities are facilitating slightly higher colonization rates at coarser spatial grains in comparison
with those of temporal turnover and extinction.
Research to date suggests that while patterns of species richness are associated with
climate mostly at coarse spatial grains (> 500 km2; Field et al. 2009, McGill 2010), species
turnover is related to climatic variability at much smaller grains. For example, Gaston et al.
(2004) found that the single best environmental predictor of spatial turnover at 10x10 km spatial
grain was temperature. In their study, the influence of temperature decreased as the spatial scale
increased and became a non-significant predictor of turnover at spatial grains larger than 30x30
km (Gaston et al. 2004). My analysis corroborates their study and extends their results by
demonstrating that influences of temporal climate change diminish at spatial scales larger than
20x20 km.
In general, higher rates of temporal turnover, extinction, and colonization were associated
with increases in temperatures recorded during the same time period at the two smallest spatial

119

scales (i.e., 5x5 and 10x10 km). This association indicates that climate change might have
partially driven the observed changes in avian communities across New York. As temperatures
increased, thermal requirements of some species were no longer met while new thermal niches
became available for other species, resulting in higher extinction and colonization, respectively.
As temperatures continue to rise in the next decades, we can expect further reshuffling of
community composition through processes of extinction and colonization and perhaps even
development of novel species assemblages (i.e., assemblages with no modern analogue;
Stralberg et al. 2009). My analysis suggests that these novel or close-to-novel communities
would initially form at the smallest spatial grains (i.e., up to 100 km2), but given enough time
coarser spatial scales might also undergo pronounced changes in community composition. Other
studies corroborate my results. For example, Stralberg et al. (2009) predicted that up to 57% of
California might be occupied by novel bird assemblages by 2070 as a result of individualistic
shifts in species distributions driven by changing climate. The percentage of novel assemblages,
however, was highly dependent on the spatial scale of investigation and decreased with
increasing spatial grain. At the largest spatial scales, novel species assemblages were predicted to
occupy only up to 3.5% of California (Stralberg et al. 2009). My results suggest that climate
change might induce similar scale-dependent changes in community composition across New
York.
Interestingly, I found that colonization was negatively related to minimum temperatures
at the 20x20 km spatial grain. This negative association with minimum temperatures indicates,
perhaps counterintuitively, that avian communities underwent highest colonization rates in
regions where temperatures changed the least, while colonization was lowest in regions that
experienced substantial increases in minimum temperatures. I also found a positive association

120

of maximum temperatures with the proportion of colonization events at this spatial scale. Given
that minimum temperatures have shown larger increases than maximum temperatures and are
expected to continue to rise more drastically (IPCC 2013), we might expect that changes in
minimum temperatures will have stronger impact on avian communities than increasing
maximum temperatures. Thus, my results suggest that future changes in temperatures may
contribute to lower levels of immigration and colonization from surrounding assemblages,
further indicating that climate change-induced temperature increases might reduce colonization
rates at the 20x20 km spatial grain. This result might provide further support for the findings
from my previous works (Chapters 1 and 2 of this volume). Specifically, I showed that
influences of climate change are often spatially non-stationary and depend on other
environmental factors, such as land-cover. Thus, while the global coefficient estimate might be
indicating a negative association between changing temperatures and colonization, it is likely
that this association varies through space and takes both positive and negative values.
My analysis suggests that land-cover is pertinent to changes in avian communities at all
spatial grains, with exception of the largest scale (i.e., 40x40 km). These findings indicate that
landscape matrix matters over larger spatial scales than changing climatic conditions. Research
to date provides conflicting accounts of the relative importance of climate and land-cover as
drivers of biodiversity change. Some studies showed that landscape matrix indeed influences
species turnover at larger spatial scales (e.g., Krawchuk and Taylor 2003), while other research
suggests that landscape structure is a valid predictor at smaller spatial scales than climatic
variability (e.g., Rahbek 2005). My results point to the importance of landscape matrix as a
potential facilitator of species responses to climate change. While climate change will likely
affect ecological communities at smaller spatial grains, landscape matrix might influence the

121

strength of these responses by facilitating or limiting the dispersal and movement of species
through the landscape. Indeed, I showed that responses of avian communities to climate change
are influenced by the level of landscape fragmentation (Chapter 3 of this volume). My results
add another dimension by showing that these landscape influences are likely to operate at scales
larger than climate change.
In general, highly fragmented regions experienced smaller changes in community
composition (especially temporal turnover and extinction), while these changes were more
pronounced in contiguous landscapes. I suggest that there are at least two alternative
interpretations for this pattern. It is plausible that communities in contiguous landscapes might be
able to respond quicker to changes in the environment due to potential for unimpeded movement.
Indeed, research suggests that unfragmented habitats might allow higher rates of species
distributional shifts associated with climate change (Opdam and Washer 2004, Jeltsch et al.
2011), which would ultimately result in more pronounced changes in community composition.
Alternatively, avian communities in fragmented landscapes may be more robust to
environmental disturbance than communities found in contiguous habitats because they comprise
largely of habitat generalists (e.g., Tscharntke et al. 2012, Estavillo et al. 2013), which tend to
have wider habitat and thermal breadths than habitat specialists (Barnagaud et al. 2012). These
wider niches potentially allow generalist birds to tolerate greater changes in environmental
conditions, thus resulting in lower rates of community change. I suggested the latter as more
plausible explanation for this ecological system (Chapter 3 of this volume) because communities
comprised of generalist species underwent smaller changes than those comprised of habitat
specialists.

122

I found that high levels of development were associated with higher values of temporal
turnover, extinction or colonization and this relationship held across most spatial scales.
Developed areas generally provide unsuitable conditions for a wide variety of different species
due to increasing pressure of human activities on the remaining patches of natural habitat (e.g.,
Huste and Boulinier 2007). This combination results in high extinction rates. This expectation is
supported by my study because I show that development influenced extinction to a higher degree
and over larger spatial grains than it did colonization. We can thus expect that increasing human
development of natural habitats will result in increased rates of extinction in the future.
I found temporal turnover to be associated with elevation, but this relationship was not
found for either extinction or colonization. At the finest spatial grain (i.e., 5x5 km), the
association of elevation with temporal turnover was positive, indicating that community change
was highest in high elevation regions. Gaston et al. (2007) found a similar relationship for world
avifauna. A plausible explanation for this increase in community change with elevation is the
relationship between turnover and species richness. Communities that are less species rich tend
to undergo higher rates of turnover (e.g., Lennon et al. 2001). In New York, high-elevation
regions of the Adirondack or Catskill Mountains generally support lower species richness than
low-laying agricultural plains, thus potentially exhibiting higher temporal community change.
Perhaps surprisingly, the influence of elevation on temporal turnover either disappeared or
became negative with increasing spatial grain. The negative association at coarser grains
indicates that high elevation regions experienced smallest changes in composition of avian
assemblages. Perhaps the topographical variability at the coarser spatial grains is no longer
associated with species richness, thus resulting in the lack of or reverse relationship between
elevation and temporal turnover.

123

4.4.3 Conclusions
My study contributes new insights into the patterns of temporal change in avian communities
across different spatial grains and drivers relevant to these changes. Changes in community
composition are scale-dependent and in general diminish with increasing spatial grain. This
decline is steeper for extinction than for colonization, likely due to two factors: (i) differences in
spatial patterns of range retraction and expansion and (ii) high mobility and dispersal abilities of
birds. The influences of climate change, landscape matrix, and elevation on temporal changes in
avian assemblages are scale-dependent. Impacts of climate change are relevant at small spatial
grains, whereas influences of landscape matrix remain important also at larger spatial grains. I
suggest that, as climate change progresses, avian assemblages will undergo increasing reshuffling at local scales, perhaps resulting in novel species assemblages. The magnitude of these
biodiversity changes resulting from climate change, however, is likely to be facilitated by
landscape fragmentation at larger scales. By establishing a link between changes in biodiversity
and environmental processes controlling these changes across different spatial scales, my study
provides a major contribution to the global change science.
4.5 ACKNOWLEDGMENTS
I would like to thank the thousands of volunteers who participated in both New York State
Breeding Bird Atlases. I also thank Kimberley Corwin and Kevin McGowan for supplying atlas
databases and other information. The manuscript benefited from discussions with members of
Quantitative Wildlife Laboratory at Michigan State University. This study received financial
support from NASA and Boone and Crockett Club.

124

EPILOGUE

In my dissertation, I sought to address the need expressed by NASA to assess the relative roles of
climate change and land cover in the observed changes in avian biodiversity. Specifically, my
goal was to examine influences of land-cover composition and configuration on the way avian
biodiversity responds to climate change and to investigate spatial scales relevant to the
understanding of the impacts of global change. In Chapter 1, I tested the utility of the spatiallyvarying coefficients (SVC) models to quantify the influence of spatial non-stationarity on
relationships between temporal community change and changes in climate and land cover. In
Chapter 2, I evaluated bioclimatic relationships of grassland and forest breeding birds across
varying gradients of their habitat. In Chapter 3, I investigated relationships between temporal
changes in community composition and the interaction of climate change and land cover
fragmentation. In Chapter 4, I investigated temporal community change across different spatial
scales and determined factors relevant to changes in avian assemblages at each of the scales.
The collective works in these chapters contribute four primary and novel conclusions for
better understanding of the implications of the interaction of land use and climate change to
avian diversity. First, relationships between changes in biodiversity and environmental factors
are often spatially non-stationary, i.e., the relationship between a response variable and the
predictor covariates varies across the spatial extent of the study. Second, the amount of suitable
land cover can significantly alter species responses to climate change, but this effect of land
cover depends on the type of habitat. Specifically, extensive forest land has potential to mitigate
the negative consequences of climate change to forest breeding birds, whereas expansive open
lands are likely to exacerbate climate change impacts on populations of grassland birds. Third,
impacts of climate change on ecological communities are more pronounced in regions of
125

unfragmented landscapes than in locations with fragmented habitats. The differences are likely
caused by the fact that communities in fragmented habitats generally tend to be comprised of
species with wider climatic niches, thus being less vulnerable to changing climatic conditions.
Fourth, temporal changes in community composition are scale-dependent and so are the
mechanisms driving these changes. Specifically, climate change operates at smaller spatial scales
than landscape fragmentation. Consequently, as temperatures continue to rise in the next
decades, we can expect reshuffling of community composition and perhaps development of
completely novel species assemblages, which will initially form at the smallest spatial grains.
Chapter 1 addressed the issue of spatial non-stationarity in ecological relationships. It is
widely recognized that the assumption of stationary regression coefficients results in poor fit and
misleading inference about the impact of the covariates on the response variable. However, the
current methods for accounting for non-stationarity tend to provide inaccurate and correlated
regression coefficient estimates because they are not robust to collinearity among the covariates
and the presence of complex spatial correlation structures. I tested the usefulness of spatiallyvarying coefficient (SVC) models to detect and account for spatial non-stationarity between
temporal community changes in avian assemblages and changes in climate and land cover. The
SVC method provides an improvement over methods that ignore spatial non-stationarity in terms
of model fit and model predictive performance. Further, by fitting these models within a
Bayesian inferential framework, I was able to make inferences about the spatial impact of
covariates and other process parameters, as well as obtain full posterior predictive inference
about the rate of turnover at new, i.e., unobserved locations.
Chapter 2 addressed the issue of variation in bioclimatic relationships of forest and
grassland birds across gradients of forest and grassland habitats. Simulation studies suggest that

126

quantity and quality of available habitat affect the way species respond to climatic variability,
but empirical evidence is lacking because studies that integrate both land cover and climate are
rare. I identified the spatial variation in bioclimatic relationships of forest and grassland breeding
birds and related this variation to changing amount of forest and grassland cover. My analysis
indicated that increasing amounts of forests weaken the negative responses of forest breeding
birds to climatic variability, while expansive grasslands strengthen these negative relationships.
My work suggests that extensive forests will mitigate consequences of climate change to avian
populations, but these negative consequences are likely to be intensified in regions with
extensive grasslands. It is thus likely that the ongoing field abandonment and reforestation in
New York will compensate or reverse the negative consequences of climate change by leading to
the dominance of cold-dwelling assemblages (Clavero et al. 2011).
Chapter 3 addressed the lack of understanding in how climate change and landscape
fragmentation interact to influence temporal changes in communities. Landscape fragmentation
is likely to influence the way communities respond to climate change, for example by affecting
the rate of species distributional shifts associated with climate change. My findings suggested
that avian communities in regions with prevalent habitat fragmentation were affected by changes
in climate to a lesser degree than communities found in unfragmented landscapes. I suggested
that avian communities in fragmented landscapes must be more robust to climate change than
communities found in contiguous habitats because they comprise of species with wider thermal
niches and are thus less susceptible to temperature increases.
Chapter 4 addressed the issue of spatial scaling of temporal community change measured
as temporal turnover, extinction, and colonization. It is widely recognized that patterns of
biodiversity vary with spatial scale of investigation, but very few studies have investigated such

127

scaling patterns across more than one time step. My analysis suggested that temporal changes in
community composition diminish with increasing spatial scales of investigation. The decline,
however, is steeper for extinction than for colonization, likely due to two factors: (i) differences
in spatial patterns of range retraction and expansion and (ii) high mobility and dispersal abilities
of birds. I showed that processes driving temporal changes in community composition also vary
with spatial scale. Impacts of climate change are relevant at small spatial grains, whereas
influences of landscape matrix remain important also at larger spatial grains. I suggested that, as
climate change progresses, avian assemblages will undergo increasing re-shuffling at local
scales, perhaps resulting in novel or close-to-novel species assemblages. The magnitude of these
changes associated with climate change, however, is likely to be facilitated by landscape
fragmentation at larger scales.
The collective works offer several novel contributions to the understanding of the impacts
of the interaction of climate change and land cover on biodiversity. First, I showed that changes
in ecological communities will be more severe in open, contiguous and unfragmented habitats
than in regions with expansive landscape fragmentation. Second, contrary to the popular belief, I
showed that climate change will operate at relatively small spatial scales. The implications of
these findings suggest that current conservation strategies may be insufficient to protect
biodiversity in the face of climate change. Existing biodiversity conservation generally focuses
on conserving large regions with undisturbed and contiguous habitats that generally support high
species diversity. My results suggest that ecological communities of such habitats will undergo
the most drastic compositional changes as a result of climate change. Thus, it is critical that the
existing conservation strategies are placed in the context of climate change. My work suggests
that in order to develop successful biodiversity conservation, one needs to consider both the

128

individual species’ vulnerabilities to climate change based on their habitat association or lifehistory strategies and relative vulnerabilities of entire communities based on the landscape in
which these communities persist. I believe that my works provide the first step into developing
such sound conservation strategies.

129

APPENDIX

130

Table A 1.1 List of all bird species (N = 256) from the New York State Breeding Bird Atlas
(BBA) used in the analysis.
BBA Code
ABDU
ACFL
AGWT
ALFL
AMBI
AMCO
AMCR
AMGO
AMKE
AMOY
AMRE
AMOR
AMWI
AMWO
AWPE
BAEA
BANS
BAOR
BARS
BAWW
BBCU
BBWA
BBWO
BCCH
BCNH
BDOW
BEKI
BGGN
BHCO
BHPA
BHVI
BITH
BLBW
BLGR
BLJA
BLPW
BLRA
BLSK
BLTE

Common Name
American Black Duck
Acadian Flycatcher
Green-winged Teal
Alder Flycatcher
American Bittern
American Coot
American Crow
American Goldfinch
American Kestrel
American Oystercatcher
American Redstart
American Robin
American Wigeon
American Woodcock
American White Pelican
Bald Eagle
Bank Swallow
Baltimore Oriole
Barn Swallow
Black-and-white Warbler
Black-billed Cuckoo
Bay-breasted Warbler
Black-backed Woodpecker
Black-capped Chickadee
Black-crowned Night-Heron
Barred Owl
Belted Kingfisher
Blue-gray Gnatcatcher
Brown-headed Cowbird
Black-hooded Parakeet
Blue-headed Vireo
Bicknell’s Thrush
Blackburnian Warbler
Blue Grosbeak
Blue Jay
Blackpoll Warbler
Black Rail
Black Skimmer
Black Tern

Scientific Name
Anas rubripes
Empidonax virescens
Anas crecca
Empidonax alnorum
Botaurus lentiginosus
Fulica americana
Corvus brachyrhynchos
Spinus tristis
Falco sparverius
Haematopus palliatus
Setophaga ruticilla
Turdus migratorius
Anas americana
Scolopax minor
Pelecanus erythrorhynchos
Haliaeetus leucocephalus
Riparia riparia
Icterus galbula
Hirundo rustica
Mniotilta varia
Coccyzus erythropthalmus
Setophaga castanea
Picoides arcticus
Poecile atricapillus
Nycticorax nycticorax
Strix varia
Megaceryle alcyon
Polioptila caerulea
Molothrus ater
Nandayus nenday
Vireo solitarius
Catharus bicknelli
Setophaga fusca
Passerina caerulea
Cyanocitta cristata
Setophaga striata
Laterallus jamaicensis
Rynchops niger
Chlidonias niger
Continued on next page

131

Table A 1.1 (cont’d)
BBA Code
BLVU
BNOW
BOCH
BRBL
BRCR
BRTH
BRWA

Common Name
Black vulture
Barn Owl
Boreal Chickadee
Brewer’s Blackbird
Brown Creeper
Brown Thrasher
Brewster’s Warbler

BTBW
BTGR
BTNW
BUFF
BHWA
BWTE
BWWA
CAEG
CAGO
CANV
CARW
CATE
CAWA
CCSP
CEDW
CERW
CHSP
CHSW
CLRA
CLSW
CMWA
COEI
COGO
COGR
COHA
COLO
COME
COMO
CONI
CORA
COSN
COTE
COYE
CSWA

Black-throated Blue Warbler
Boat-tailed Grackle
Black-throated Green Warbler
Bufflehead
Broad-winged Hawk
Blue-winged Teal
Blue-winged Warbler
Cattle Egret
Canada Goose
Canvasback
Carolina Wren
Caspian Tern
Canada Warbler
Clay-colored Sparrow
Cedar Waxwing
Cerulean Warbler
Chipping Sparrow
Chimney Swift
Clapper Rail
Cliff Swallow
Cape May Warbler
Common Eider
Common Goldeneye
Common Grackle
Cooper’s Hawk
Common Loon
Common Merganser
Common Moorhen
Common Nighthawk
Common Raven
Wilson’s Snipe
Common Tern
Common Yellowthroat
Chestnut-sided Warbler

Scientific Name
Coragyps atratus
Tyto alba
Poecile hudsonicus
Euphagus cyanocephalus
Certhia Americana
Toxostoma rufum
Vermivora chrysoptera x V.
cyanoptera
Setophaga caerulescens
Quiscalus major
Setophaga virens
Bucephala albeola
Buteo platypterus
Anas discors
Vermivora cyanoptera
Bubulcus ibis
Branta Canadensis
Aythya valisineria
Thryothorus ludovicianus
Hydroprogne caspia
Cardellina Canadensis
Spizella pallid
Bombycilla cedrorum
Setophaga cerulean
Spizella passerine
Chaetura pelagic
Rallus longirostris
Petrochelidon pyrrhonota
Setophaga tigrina
Somateria mollissima
Bucephala clangula
Quiscalus quiscula
Accipiter cooperii
Gavia immer
Mergus merganser
Gallinula galeata
Chordeiles minor
Corvus corax
Gallinago delicate
Sterna hirundo
Geothlypis trichas
Setophaga pensylvanica
Continued on next page
132

Table A 1.1 (cont’d)
BBA Code
CWWI
DCCO
DICK
DOWO
EABL
EAKI
EAME
EAPH
EASO
EATO
EAWP
ECDO
ETTI
EUST
EVGR
FICR
FISP
FOTE
GADW
GBBG
GBHE
GBTE
GCFL
GCKI
GHOW
GLIB
GOEA
GRAJ
GRCA
GREG
GRHE
GRPA
GRSC
GRSP
GWWA
HAWO
HERG
HESP
HETH
HOFI
HOLA
HOME

Common Name
Chuck-will’s-widow
Double-crested Cormorant
Dickcissel
Downy Woodpecker
Eastern Bluebird
Eastern Kingbird
Eastern Meadowlark
Eastern Phoebe
Eastern Screech-Owl
Eastern Towhee
Eastern Wood-Pewee
Eurasian Collared-Dove
Tufted Titmouse
European Starling
Evening Grosbeak
Fish Crow
Field Sparrow
Forster’s Tern
Gadwall
Great Black-backed Gull
Great Blue Heron
Gull-billed Tern
Great Crested Flycatcher
Golden-crowned Kinglet
Great Horned Owl
Glossy Ibis
Golden Eagle
Gray Jay
Gray Catbird
Great Egret
Green Heron
Gray Partridge
Greater Scaup
Grasshopper Sparrow
Golden-winged Warbler
Hairy Woodpecker
Herring Gull
Henslow’s Sparrow
Hermit Thrush
House Finch
Horned Lark
Hooded Merganser

Scientific Name
Antrostomus carolinensis
Phalacrocorax auritus
Spiza Americana
Picoides pubescens
Sialia sialis
Tyrannus tyrannus
Sturnella magna
Sayornis phoebe
Megascops asio
Pipilo erythrophthalmus
Contopus virens
Streptopelia decaocto
Baeolophus bicolor
Sturnus vulgaris
Coccothraustes vespertinus
Corvus ossifragus
Spizella pusilla
Sterna forsteri
Anas strepera
Larus marinus
Ardea Herodias
Gelochelidon nilotica
Myiarchus crinitus
Regulus satrapa
Bubo virginianus
Plegadis falcinellus
Aquila chrysaetos
Perisoreus Canadensis
Dumetella carolinensis
Ardea alba
Butorides virescens
Perdix perdix
Aythya marila
Ammodramus savannarum
Vermivora chrysoptera
Picoides villosus
Larus argentatus
Ammodramus henslowii
Catharus guttatus
Haemorhous mexicanus
Eremophila alpestris
Lophodytes cucullatus
Continued on next page
133

Table A 1.1 (cont’d)
BBA Code
HOSP
HOWA
HOWR
INBU
KEWA
KILL
KIRA
LAGU
LAWA

Common Name
House Sparrow
Hooded Warbler
House Wren
Indigo Bunting
Kentucky Warbler
Killdeer
King Rail
Laughing Gull
Lawrence’s Warbler

LBHE
LEBI
LEFL
LEOW
LESC
LETE
LISP
LOSH
LOWA
MALL
MAWA
MAWR
MBDH

Little Blue Heron
Least Bittern
Least Flycatcher
Long-eared Owl
Lesser Scaup
Least Tern
Lincoln’s Sparrow
Loggerhead Shrike
Louisiana Waterthrush
Mallard
Magnolia Warbler
Marsh Wren
Mallard x American Black Duck
Hybrid
Merlin
Mourning Dove
Monk Parakeet
Mourning Warbler
Mute Swan
Yellow-rumped Warbler
Nashville Warbler
Northern Bobwhite
Northern Cardinal
Northern Goshawk
Northern Harrier
Northern Mockingbird
Northern Parula
Northern Pintail
Northern Waterthrush
Northern Rough-winged Swallow
Northern Shoveler
Nelson’s Sharp-tailed Sparrow

MERL
MODO
MOPA
MOWA
MUSW
MYWA
NAWA
NOBO
NOCA
NOGO
NOHA
NOMO
NOPA
NOPI
NOWA
NRWS
NSHO
NSTS

134

Scientific Name
Passer domesticus
Setophaga citrine
Troglodytes aedon
Passerina cyanea
Geothlypis Formosa
Charadrius vociferous
Rallus elegans
Leucophaeus atricilla
Vermivora chrysoptera x V.
cyanoptera
Egretta caerulea
Ixobrychus exilis
Empidonax minimus
Asio otus
Aythya affinis
Sternula antillarum
Melospiza lincolnii
Lanius ludovicianus
Parkesia motacilla
Anas platyrhynchos
Setophaga magnolia
Cistothorus palustris
Anas platyrhynchos x A.rubripes
Falco columbarius
Zenaida macroura
Myiopsitta monachus
Geothlypis Philadelphia
Cygnus olor
Setophaga coronate
Oreothlypis ruficapilla
Colinus virginianus
Cardinalis cardinalis
Accipiter gentilis
Circus cyaneus
Mimus polyglottos
Setophaga Americana
Anas acuta
Parkesia noveboracensis
Stelgidopteryx serripennis
Anas clypeata
Ammodramus nelson
Continued on next page

Table A 1.1 (cont’d)
BBA Code
NSWO
OROR
OSFL
OSPR
OVEN
PBGR
PEFA
PHVI
PIPL
PISI
PIWA
PIWO
PRAW
PROW
PUFI
PUMA
RBGR
RBGU
RBME
RBNU
RBWO
RCKI
RECR
REDH
REVI
RHWO
RNDU
RODO
ROST
RPHE
RSHA
RTHA
RTHU
RUBL
RUDU
RUGR
RWBL
SACR
SAVS
SCJU
SCTA
SEOW

Common Name
Northern Saw-whet Owl
Orchard Oriole
Olive-sided Flycatcher
Osprey
Ovenbird
Pied-billed Grebe
Peregrine Falcon
Philadelphia Vireo
Piping Plover
Pine Siskin
Pine Warbler
Pileated Woodpecker
Prairie Warbler
Prothonotary Warbler
Purple Finch
Purple Martin
Rose-breasted Grosbeak
Ring-billed Gull
Red-breasted Merganser
Red-breasted Nuthatch
Red-bellied Woodpecker
Ruby-crowned Kinglet
Red Crossbill
Redhead
Red-eyed Vireo
Red-headed Woodpecker
Ring-necked duck
Rock Pigeon
Roseate Tern
Ring-necked Pheasant
Red-shouldered Hawk
Red-tailed Hawk
Ruby-throated Hummingbird
Rusty Blackbird
Ruddy Duck
Ruffed Grouse
Red-winged Blackbird
Sandhill Crane
Savannah Sparrow
Dark-eyed Junco
Scarlet Tanager
Short-eared Owl

Scientific Name
Aegolius acadicus
Icterus spurious
Contopus cooperi
Pandion haliaetus
Seiurus aurocapilla
Podilymbus podiceps
Falco peregrinus
Vireo philadelphicus
Charadrius melodus
Spinus pinus
Setophaga pinus
Dryocopus pileatus
Setophaga discolor
Protonotaria citrea
Haemorhous purpureus
Progne subis
Pheucticus ludovicianus
Larus delawarensis
Mergus serrator
Sitta Canadensis
Melanerpes carolinus
Regulus calendula
Loxia curvirostra
Aythya Americana
Vireo olivaceus
Melanerpes erythrocephalus
Aythya collaris
Columba livia
Sterna dougallii
Phasianus colchicus
Buteo lineatus
Buteo jamaicensis
Archilochus colubris
Euphagus carolinus
Oxyura jamaicensis
Bonasa umbellus
Agelaius phoeniceus
Grus Canadensis
Passerculus sandwichensis
Junco hyemalis
Piranga olivacea
Asio flammeus
Continued on next page
135

Table A 1.1 (cont’d)
BBA Code
SESP
SEWR
SNEG
SORA
SOSP
SPGR
SPSA
SSHA
SSTS
SUTA
SWSP
SWTH
TEWA
TRES
TRHE
TRUS
TTWO
TUVU
UPSA
VEER
VESP
VIRA
WAVI
WBNU
WEKI
WEME
WEVI
WEWA
WFIB
WIFL
WILL
WIPH
WITU
WIWA
WIWR
WODU
WOTH
WPWI
WTSP
WWCR
WWTE
YBCH

Common Name
Seaside Sparrow
Sedge Wren
Snowy Egret
Sora
Song Sparrow
Spruce Grouse
Spotted Sandpiper
Sharp-shinned Hawk
Saltmarsh Sharp-tailed Sparrow
Summer Tanager
Swamp Sparrow
Swaison’s Thrush
Tennessee Warbler
Tree Swallow
Tricolored Heron
Trumpeter Swan
American Three-toed Woodpecker
Turkey Vulture
Upland Sandpiper
Veery
Vesper Sparrow
Virginia Rail
Warbling Vireo
White-breasted Nuthatch
Western Kingbird
Western Meadowlark
White-eyed Vireo
Worm-eating Warbler
White-faced Ibis
Willow Flycatcher
Willet
Wilson’s Phalarope
Wild Turkey
Wilson’s Warbler
Winter Wren
Wood Duck
Wood Thrush
Whip-poor-will
White-throated Sparrow
White-winged Crossbill
White-winged Tern
Yellow-breasted Chat

136

Scientific Name
Ammodramus maritimus
Cistothorus platensis
Egretta thula
Porzana Carolina
Melospiza melodia
Falcipennis Canadensis
Actitis macularius
Accipiter striatus
Ammodramus caudacutus
Piranga rubra
Melospiza Georgiana
Catharus ustulatus
Oreothlypis peregrine
Tachycineta bicolor
Egretta tricolor
Cygnus buccinators
Picoides dorsalis
Cathartes aura
Bartramia longicauda
Catharus fuscescens
Pooecetes gramineus
Rallus limicola
Vireo gilvus
Sitta carolinensis
Tyrannus verticalis
Sturnella neglecta
Vireo griseus
Helmitheros vermivorum
Plegadis chihi
Empidonax traillii
Tringa semipalmata
Phalaropus tricolor
Meleagris gallopavo
Cardellina pusilla
Troglodytes hiemalis
Aix sponsa
Hylocichla mustelina
Antrostomus vociferous
Zonotrichia albicollis
Loxia leucoptera
Chlidonias leucopterus
Icteria virens
Continued on next page

Table A 1.1 (cont’d)
BBA Code
YBCH
YBCU
YBFL
YBSA
YCNH
YPWA
YSFL
YTVI
YTWA
YWAR

Common Name
Yellow-breasted Chat
Yellow-billed Cuckoo
Yellow-bellied Flycatcher
Yellow-bellied Sapsucker
Yellow-crowned Night Heron
Palm Warbler
Northern Flicker
Yellow-throated Vireo
Yellow-throated Warbler
Yellow Warbler

Scientific Name
Icteria virens
Coccyzus americanus
Empidonax flaviventris
Sphyrapicus varius
Nyctanassa violacea
Setophaga palmarum
Colaptes auratus
Vireo flavifrons
Setophaga dominica
Setophaga petechia

137

Table A 2.1 Comparison of the spatially-varying intercept (SVI) and spatially-varying
coefficients (SVC) models for grassland and forest birds recorded during the 2000-05 New York
State Breeding Bird Atlas.
Habitat
Guild

Species
Common Name
American
Kestrel
Northern
Harrier
Eastern
Bluebird
Eastern
Kingbird
Horned Lark

Grassland Grasshopper
Birds
Sparrow
Savannah
Sparrow
Vesper Sparrow
Bobolink

Forest
Birds

Brown-headed
Cowbird
Eastern
Meadowlark
Blue-gray
Gnatcatcher
Blue-headed
Vireo
Goldencrowned
Kinglet
Great Crested
Flycatcher
Red-breasted
Nuthatch

Model

Scientific
Name
Falco
sparverius

Spatially-Varying Coefficients
Spatially-Varying Intercept
Spatially-Varying Coefficients
Circus cyaneus
Spatially-Varying Intercept
Spatially-Varying Coefficients
Sialia sialis
Spatially-Varying Intercept
Tyrannus
Spatially-Varying Coefficients
tyrannus
Spatially-Varying Intercept
Spatially-Varying Coefficients
Eremophila
alpestris
Spatially-Varying Intercept
Spatially-Varying Coefficients
Ammodramus
savannarum
Spatially-Varying Intercept
Passerculus
Spatially-Varying Coefficients
sandwichensis Spatially-Varying Intercept
Pooecetes
Spatially-Varying Coefficients
gramineus
Spatially-Varying Intercept
Spatially-Varying Coefficients
Dolichonyx
oryzivorus
Spatially-Varying Intercept
Spatially-Varying Coefficients
Molothrus ater
Spatially-Varying Intercept
Sturnella
Spatially-Varying Coefficients
magna
Spatially-Varying Intercept
Polioptila
Spatially-Varying Coefficients
caerulea
Spatially-Varying Intercept
Spatially-Varying Coefficients
Vireo solitaries
Spatially-Varying Intercept
Regulus
Spatially-Varying Coefficients
satrapa
Spatially-Varying Intercept
Myiarchus
crinitus
Sitta
canadensis

ΔDIC

0
175.575
0
156.97
0
132.33
0
79.087
0
42.339
0
30.121
0
169.515
0
52.299
0
240.449
0
87.519
0
162.862
0
73.544
0
112.179
0
18.148

Spatially-Varying Coefficients
0
Spatially-Varying Intercept
99.505
Spatially-Varying Coefficients
0
Spatially-Varying Intercept
114.246
Continued on next page

138

Table A 2.1 (cont’d)
Habitat
Guild

Species
Common Name
Tufted
Titmouse
Veery

Forest
Birds

Black-andwhite Warbler
Black-throated
Blue Warbler
Black-throated
Green Warbler
Blackburnian
Warbler

Scientific
Name
Baeolophus
bicolor
Catharus
fuscescens
Mniotilta varia

Setophaga
caerulescens
Dendroica
virens
Dendroica
fusca
Cardellina
Canada Warbler
canadensis

Model

Spatially-Varying Coefficients
Spatially-Varying Intercept
Spatially-Varying Coefficients
Spatially-Varying Intercept
Spatially-Varying Coefficients
Spatially-Varying Intercept
Spatially-Varying Coefficients
Spatially-Varying Intercept
Spatially-Varying Coefficients
Spatially-Varying Intercept
Spatially-Varying Coefficients
Spatially-Varying Intercept
Spatially-Varying Coefficients
Spatially-Varying Intercept

139

ΔDIC

0
138.906
0
133.793
0
115.553
0
50.915
0
82.678
0
154.844
0
96.563

Table A 3.1 Coefficient estimates (βs) of each covariate resulting from all tested models of
temporal turnover (starting with the best model as indicated by DIC) for four avian groupings: all
species, long-distance migrants, short-distance migrants, and resident species. Covariates
included intercept (INT), magnitude of the 25-year (1980-2005) trend in average maximum
temperature of the breeding season (TMAX), magnitude of the 25-year (1980-2005) trend in
average minimum temperature of the breeding season (TMIN), magnitude of the 25-year (19802005) trend in average total precipitation of the breeding season (PRECIP), edge density (ED),
percent developed land (DEVEL), interaction between the magnitude of the 25-year (1980-2005)
trend in average maximum temperature of the breeding season and edge density (TMAX*ED),
interaction between the magnitude of the 25-year (1980-2005) trend in average minimum
temperature of the breeding season and edge density (TMIN*ED), interaction between
magnitude of the 25-year (1980-2005) trend in average total precipitation of the breeding season
and edge density (PRECIP*ED), and survey effort (EFF). 50% indicates the mean of the
posterior distribution, while 2.5% and 97.5% are the lower and upper quantiles of the posterior
distribution.
Avian
Grouping

Model

TMAX*ED +
TMIN*ED +
PRECIP*ED +
DEVEL + EFF

All Species

TMAX + TMIN +
PRECIP + ED +
DEVEL + EFF

TMAX*ED +
TMIN*ED +
PRECIP*ED + EFF

Coefficient
βINT
βTMAX
βTMIN
βPRECIP
βED
βDEVEL
βTMAX-ED
βTMIN-ED
βPRECIP-ED
βEFF
β0
βTMAX
βTMIN
βPRECIP
βED
βDEVEL
βEFF
βINT
βTMAX
βTMIN
βPRECIP
βED
βTMAX-ED
βTMIN-ED
βPRECIP-ED
βEFF

140

Coefficient Estimates
50%
2.5%
97.5%
-0.368
-0.442
-0.294
0.018
0.005
0.030
0.020
0.008
0.032
-0.014
-0.034
0.005
-0.043
-0.053
-0.032
0.034
0.023
0.045
-0.003
-0.012
0.006
-0.017
-0.025
-0.008
0.070
-0.002
0.018
-0.016
-0.022
-0.009
-0.397
-0.466
-0.318
0.019
0.008
0.031
0.024
0.012
0.037
-0.019
-0.038
-0.002
-0.041
-0.054
-0.030
0.035
0.024
0.046
-0.016
-0.022
-0.009
-0.351
-0.433
-0.267
0.014
0.002
0.025
0.020
0.007
0.032
-0.012
-0.028
0.004
-0.046
-0.056
-0.035
-0.001
-0.011
0.008
-0.015
-0.025
-0.006
0.012
0.002
0.021
-0.015
-0.022
-0.008
Continued on next page

Table A 3.1 (cont’d)
Avian
Grouping

Model

TMAX + TMIN +
PRECIP + ED + EFF
All Species
TMAX + TMIN +
PRECIP + EFF

TMAX*ED +
TMIN*ED +
PRECIP*ED +
DEVEL + EFF

Long-distance
Migrants

TMAX + TMIN +
PRECIP + ED +
DEVEL + EFF

TMAX*ED +
TMIN*ED +
PRECIP*ED + EFF

Coefficient
β0
βTMAX
βTMIN
βPRECIP
βED
βEFF
β0
βTMAX
βTMIN
βPRECIP
βEFF
βINT
βTMAX
βTMIN
βPRECIP
βED
βDEVEL
βTMAX-ED
βTMIN-ED
βPRECIP-ED
βEFF
β0
βTMAX
βTMIN
βPRECIP
βED
βDEVEL
βEFF
βINT
βTMAX
βTMIN
βPRECIP
βED
βTMAX-ED
βTMIN-ED
βPRECIP-ED
βEFF

141

Coefficient Estimates
50%
2.5%
97.5%
-0.366
-0.431
-0.299
0.014
0.002
0.028
0.022
0.010
0.034
-0.017
-0.033
-0.002
-0.043
-0.054
-0.031
-0.015
-0.022
-0.009
-0.379
-0.451
-0.327
0.019
0.007
0.031
0.023
0.010
0.034
-0.027
-0.041
-0.009
-0.015
-0.021
-0.009
-0.417
-0.507
-0.188
0.016
-0.003
0.034
0.033
0.015
0.050
-0.001
-0.028
0.026
0.008
-0.008
0.024
0.074
0.058
0.092
0.002
-0.011
0.016
-0.016
-0.030
-0.003
0.019
0.003
0.030
-0.016
-0.026
-0.006
-0.433
-0.523
-0.303
0.015
-0.003
0.032
0.035
0.018
0.049
-0.012
-0.038
0.011
0.011
-0.006
0.027
0.074
0.057
0.092
-0.016
-0.026
-0.006
-0.370
-0.447
-0.285
0.009
-0.009
0.026
0.027
0.012
0.044
0.006
-0.018
0.030
0.005
-0.010
0.022
0.006
-0.008
0.019
-0.015
-0.028
-0.003
0.024
0.009
0.038
-0.014
-0.024
-0.003
Continued on next page

Table A 3.1 (cont’d)
Avian
Grouping

Model

TMAX + TMIN +
PRECIP + ED + EFF
Long-distance
Migrants
TMAX + TMIN +
PRECIP + EFF

TMAX + TMIN +
PRECIP + ED +
DEVEL + EFF

Short-distance
Migrants

TMAX*ED +
TMIN*ED +
PRECIP*ED +
DEVEL + EFF

TMAX + TMIN +
PRECIP + ED + EFF

Coefficient
β0
βTMAX
βTMIN
βPRECIP
βED
βEFF
β0
βTMAX
βTMIN
βPRECIP
βEFF
β0
βTMAX
βTMIN
βPRECIP
βED
βDEVEL
βEFF
βINT
βTMAX
βTMIN
βPRECIP
βED
βDEVEL
βTMAX-ED
βTMIN-ED
βPRECIP-ED
βEFF
β0
βTMAX
βTMIN
βPRECIP
βED
βEFF

142

Coefficient Estimates
50%
2.5%
97.5%
-0.277
-0.407
-0.125
0.009
-0.008
0.026
0.032
0.014
0.050
-0.010
-0.036
0.017
0.010
-0.006
0.027
-0.014
-0.025
-0.005
-0.405
-0.493
-0.263
0.007
-0.013
0.023
0.031
0.014
0.050
-0.008
-0.033
0.014
-0.014
-0.024
-0.003
-0.598
-0.656
-0.513
0.029
0.005
0.051
0.015
-0.004
0.037
0.001
-0.028
0.027
-0.068
-0.089
-0.046
0.058
0.039
0.079
-0.008
-0.022
0.004
-0.584
-0.640
-0.529
0.027
0.003
0.050
0.016
-0.008
0.035
0.000
-0.033
0.032
-0.070
-0.093
-0.046
0.058
0.037
0.079
0.004
-0.014
0.022
-0.014
-0.032
0.004
-0.004
-0.022
0.017
-0.009
-0.022
0.003
-0.579
-0.656
-0.518
0.019
-0.004
0.042
0.009
-0.015
0.033
0.002
-0.030
0.033
-0.072
-0.094
-0.046
-0.007
-0.020
0.006
Continued on next page

Table A 3.1 (cont’d)
Avian
Grouping

Model

TMAX*ED +
TMIN*ED +
PRECIP*ED + EFF
Short-distance
Migrants

TMAX + TMIN +
PRECIP + EFF

TMAX*ED +
TMIN*ED +
PRECIP*ED +
DEVEL + EFF

Resident
Species

TMAX*ED +
TMIN*ED +
PRECIP*ED + EFF

TMAX + TMIN +
PRECIP + ED +
DEVEL + EFF

Coefficient
βINT
βTMAX
βTMIN
βPRECIP
βED
βTMAX-ED
βTMIN-ED
βPRECIP-ED
βEFF
β0
βTMAX
βTMIN
βPRECIP
βEFF
βINT
βTMAX
βTMIN
βPRECIP
βED
βDEVEL
βTMAX-ED
βTMIN-ED
βPRECIP-ED
βEFF
βINT
βTMAX
βTMIN
βPRECIP
βED
βTMAX-ED
βTMIN-ED
βPRECIP-ED
βEFF
β0
βTMAX
βTMIN
βPRECIP
βED
βDEVEL
βEFF

Continued on next page
143

Coefficient Estimates
50%
2.5%
-0.562
-0.632
0.021
-0.002
0.011
-0.011
0.003
-0.029
-0.072
-0.094
0.008
-0.011
-0.012
-0.031
0.002
-0.019
-0.007
-0.020
-0.501
-0.576
0.032
0.009
0.010
-0.013
-0.008
-0.048
-0.009
-0.024
-0.232
-0.331
0.026
0.005
0.009
-0.013
-0.041
-0.071
-0.105
-0.128
-0.027
-0.047
-0.016
-0.032
-0.017
-0.034
0.002
-0.015
-0.024
-0.036
-0.169
-0.349
0.024
0.004
0.009
-0.012
-0.044
-0.074
-0.104
-0.123
-0.018
-0.035
-0.018
-0.035
-0.001
-0.018
-0.025
-0.037
-0.244
-0.350
0.028
0.006
0.014
-0.005
-0.045
-0.071
-0.103
-0.122
-0.027
-0.046
-0.024
-0.037

97.5%
-0.456
0.042
0.033
0.039
-0.050
0.026
0.006
0.021
0.005
-0.385
0.053
0.031
0.024
0.004
-0.101
0.046
0.028
-0.016
-0.086
-0.010
-0.002
-0.002
0.021
-0.012
-0.072
0.047
0.033
-0.014
-0.085
-0.003
-0.002
0.017
-0.014
-0.097
0.051
0.035
-0.019
-0.083
-0.009
-0.012

Table A 3.1 (cont’d)
Avian
Grouping

Model

TMAX + TMIN +
PRECIP + ED + EFF
Resident
Species
TMAX + TMIN +
PRECIP + EFF

Coefficient
β0
βTMAX
βTMIN
βPRECIP
βED
βEFF
β0
βTMAX
βTMIN
βPRECIP
βEFF

144

Coefficient Estimates
50%
2.5%
-0.202
-0.337
0.032
0.011
0.015
-0.005
-0.041
-0.070
-0.102
-0.121
-0.025
-0.037
-0.274
-0.429
0.042
0.021
0.012
-0.009
-0.064
-0.093
-0.025
-0.037

97.5%
-0.112
0.052
0.034
-0.008
0.082
-0.013
-0.201
0.063
0.036
-0.035
-0.014

Table A 3.2 Coefficient estimates (βs) of each covariate resulting from all tested models of
extinction (starting with the best model as indicated by DIC) for four avian groupings: all
species, long-distance migrants, short-distance migrants, and resident species. Covariates
included intercept (INT), magnitude of the 25-year (1980-2005) trend in average maximum
temperature of the breeding season (TMAX), magnitude of the 25-year (1980-2005) trend in
average minimum temperature of the breeding season (TMIN), magnitude of the 25-year (19802005) trend in average total precipitation of the breeding season (PRECIP), edge density (ED),
percent developed land (DEVEL), interaction between the magnitude of the 25-year (1980-2005)
trend in average maximum temperature of the breeding season and edge density (TMAX*ED),
interaction between the magnitude of the 25-year (1980-2005) trend in average minimum
temperature of the breeding season and edge density (TMIN*ED), interaction between
magnitude of the 25-year (1980-2005) trend in average total precipitation of the breeding season
and edge density (PRECIP*ED), and survey effort (EFF). 50% indicates the mean of the
posterior distribution, while 2.5% and 97.5% are the lower and upper quantiles of the posterior
distribution.
Avian
Grouping

Model

TMAX*ED +
TMIN*ED +
PRECIP*ED +
DEVEL + EFF

All Species

TMAX*ED +
TMIN*ED +
PRECIP*ED + EFF

TMAX + TMIN +
PRECIP + ED +
DEVEL + EFF

Coefficient
βINT
βTMAX
βTMIN
βPRECIP
βED
βDEVEL
βTMAX-ED
βTMIN-ED
βPRECIP-ED
βEFF
βINT
βTMAX
βTMIN
βPRECIP
βED
βTMAX-ED
βTMIN-ED
βPRECIP-ED
βEFF
β0
βTMAX
βTMIN
βPRECIP
βED
βDEVEL
βEFF

145

Coefficient Estimates
50%
2.5%
97.5%
-1.016
-1.088
-0.923
0.007
-0.008
0.022
0.044
0.027
0.061
0.036
0.014
0.061
-0.079
-0.093
-0.063
0.043
0.029
0.058
-0.002
-0.014
0.011
-0.040
-0.052
-0.025
0.010
-0.003
0.026
-0.031
-0.040
-0.022
-1.039
-1.105
-0.963
0.000
-0.017
0.019
0.040
0.025
0.055
0.039
0.016
0.063
-0.082
-0.098
-0.068
0.000
-0.013
0.013
-0.038
-0.052
-0.026
0.014
0.001
0.027
-0.030
-0.039
-0.020
-0.948
-1.086
-0.832
0.010
-0.007
0.026
0.053
0.039
0.070
0.029
0.008
0.051
-0.074
-0.088
-0.059
0.043
0.029
0.056
-0.030
-0.039
-0.021
Continued on next page

Table A 3.2 (cont’d)
Avian
Grouping

Model

TMAX + TMIN +
PRECIP + ED + EFF
All Species
TMAX + TMIN +
PRECIP + EFF

TMAX*ED +
TMIN*ED +
PRECIP*ED +
DEVEL + EFF

Long-distance
Migrants

TMAX + TMIN +
PRECIP + ED +
DEVEL + EFF

TMAX*ED +
TMIN*ED +
PRECIP*ED + EFF

Coefficient
β0
βTMAX
βTMIN
βPRECIP
βED
βEFF
β0
βTMAX
βTMIN
βPRECIP
βEFF
βINT
βTMAX
βTMIN
βPRECIP
βED
βDEVEL
βTMAX-ED
βTMIN-ED
βPRECIP-ED
βEFF
β0
βTMAX
βTMIN
βPRECIP
βED
βDEVEL
βEFF
βINT
βTMAX
βTMIN
βPRECIP
βED
βTMAX-ED
βTMIN-ED
βPRECIP-ED
βEFF

146

Coefficient Estimates
50%
2.5%
97.5%
-0.984
-1.079
-0.887
0.004
-0.012
0.021
0.052
0.037
0.066
0.029
0.007
0.051
-0.075
-0.091
-0.061
-0.030
-0.038
-0.021
-1.014
-1.114
-0.843
0.013
-0.006
0.030
0.053
0.038
0.070
0.012
-0.012
0.033
-0.031
-0.040
-0.022
-1.063
-1.176
-0.977
-0.024
-0.049
-0.001
0.056
0.033
0.083
0.044
0.010
0.080
-0.029
-0.050
-0.006
0.087
0.065
0.109
-0.008
-0.026
0.011
-0.032
-0.050
-0.013
0.026
0.007
0.046
-0.024
-0.037
-0.010
-1.101
-1.171
-1.044
-0.021
-0.047
0.004
0.069
0.047
0.093
0.030
0.000
0.060
-0.023
-0.044
0.001
0.087
0.066
0.110
-0.024
-0.037
-0.011
-0.985
-1.065
-0.896
-0.031
-0.056
-0.009
0.054
0.030
0.084
0.051
0.012
0.082
-0.034
-0.056
-0.012
-0.005
-0.023
0.013
-0.028
-0.046
-0.010
0.035
0.017
0.054
-0.023
-0.035
-0.010
Continued on next page

Table A 3.2 (cont’d)
Avian
Grouping

Model

TMAX + TMIN +
PRECIP + EFF
Long-distance
Migrants
TMAX + TMIN +
PRECIP + ED + EFF

TMAX*ED +
TMIN*ED +
PRECIP*ED +
DEVEL + EFF

Short-distance
Migrants

TMAX + TMIN +
PRECIP + ED +
DEVEL + EFF

TMAX + TMIN +
PRECIP + ED + EFF

Coefficient
β0
βTMAX
βTMIN
βPRECIP
βEFF
β0
βTMAX
βTMIN
βPRECIP
βED
βEFF
βINT
βTMAX
βTMIN
βPRECIP
βED
βDEVEL
βTMAX-ED
βTMIN-ED
βPRECIP-ED
βEFF
β0
βTMAX
βTMIN
βPRECIP
βED
βDEVEL
βEFF
β0
βTMAX
βTMIN
βPRECIP
βED
βEFF

147

Coefficient Estimates
50%
2.5%
97.5%
-1.009
-1.134
-0.850
-0.023
-0.046
-0.004
0.066
0.042
0.091
0.029
-0.003
0.061
-0.022
-0.035
-0.010
-1.125
-1.398
-0.962
-0.027
-0.051
-0.005
0.065
0.042
0.090
0.031
-0.004
0.061
-0.024
-0.044
-0.002
-0.023
-0.037
-0.010
-1.302
-1.441
-1.215
0.031
-0.001
0.060
0.024
-0.007
0.059
0.040
-0.005
0.078
-0.126
-0.153
-0.097
0.067
0.042
0.096
0.014
-0.010
0.040
-0.037
-0.062
-0.010
-0.013
-0.038
0.013
-0.032
-0.050
-0.013
-1.307
-1.423
-1.181
0.031
0.001
0.062
0.027
-0.001
0.059
0.039
0.000
0.087
-0.125
-0.156
-0.094
0.068
0.038
0.096
-0.032
-0.050
-0.015
-1.179
-1.285
-1.108
0.023
-0.010
0.053
0.027
-0.003
0.057
0.039
-0.002
0.077
-0.127
-0.157
-0.099
-0.031
-0.049
-0.012
Continued on next page

Table A 3.2 (cont’d)
Avian
Grouping

Model

TMAX*ED +
TMIN*ED +
PRECIP*ED + EFF
Short-distance
Migrants

TMAX + TMIN +
PRECIP + EFF

TMAX*ED +
TMIN*ED +
PRECIP*ED +
DEVEL + EFF

Resident
Species

TMAX*ED +
TMIN*ED +
PRECIP*ED + EFF

TMAX + TMIN +
PRECIP + ED + EFF

Coefficient
βINT
βTMAX
βTMIN
βPRECIP
βED
βTMAX-ED
βTMIN-ED
βPRECIP-ED
βEFF
β0
βTMAX
βTMIN
βPRECIP
βEFF
βINT
βTMAX
βTMIN
βPRECIP
βED
βDEVEL
βTMAX-ED
βTMIN-ED
βPRECIP-ED
βEFF
βINT
βTMAX
βTMIN
βPRECIP
βED
βTMAX-ED
βTMIN-ED
βPRECIP-ED
βEFF
β0
βTMAX
βTMIN
βPRECIP
βED
βEFF

148

Coefficient Estimates
50%
2.5%
97.5%
-1.300
-1.389
-0.981
0.022
-0.009
0.059
0.021
-0.015
0.054
0.035
-0.007
0.073
-0.129
-0.159
-0.101
0.016
-0.010
0.042
-0.034
-0.060
-0.008
-0.007
-0.034
0.019
-0.030
-0.050
-0.011
-1.223
-1.397
-1.014
0.035
0.005
0.065
0.029
-0.006
0.062
0.008
-0.031
0.056
-0.032
-0.050
-0.014
-1.209
-1.351
-1.015
0.035
0.003
0.068
0.019
-0.011
0.052
0.035
-0.009
0.071
-0.144
-0.174
-0.133
-0.032
-0.058
-0.005
-0.010
-0.033
0.014
-0.047
-0.073
-0.019
0.010
-0.019
0.034
-0.041
-0.058
-0.024
-1.210
-1.347
-1.055
0.039
0.011
0.068
0.023
-0.010
0.054
0.036
-0.010
0.082
-0.143
-0.172
-0.115
-0.010
-0.033
0.011
-0.048
-0.071
-0.024
0.008
-0.018
0.034
-0.042
-0.061
-0.025
-1.206
-1.300
-1.034
0.046
0.012
0.079
0.034
0.005
0.061
0.021
-0.016
0.064
-0.131
-0.160
-0.104
-0.041
-0.059
-0.024
Continued on next page

Table A 3.2 (cont’d)
Avian
Grouping

Model

TMAX + TMIN +
PRECIP + ED +
DEVEL + EFF
Resident
Species
TMAX + TMIN +
PRECIP + EFF

Coefficient
β0
βTMAX
βTMIN
βPRECIP
βED
βDEVEL
βEFF
β0
βTMAX
βTMIN
βPRECIP
βEFF

149

Coefficient Estimates
50%
2.5%
-1.194
-1.308
0.040
0.007
0.033
0.000
0.026
-0.017
-0.136
-0.165
-0.034
-0.063
-0.041
-0.059
-1.253
-1.394
0.059
0.029
0.035
0.005
-0.009
-0.045
-0.044
-0.061

97.5%
-1.044
-0.074
0.065
0.069
-0.110
-0.007
-0.025
-1.135
0.092
0.064
0.033
-0.027

Table A 3.3 Coefficient estimates (βs) of each covariate resulting from all tested models of
colonization (starting with the best model as indicated by DIC) for four avian groupings: all
species, long-distance migrants, short-distance migrants, and resident species. Covariates
included intercept (INT), magnitude of the 25-year (1980-2005) trend in average maximum
temperature of the breeding season (TMAX), magnitude of the 25-year (1980-2005) trend in
average minimum temperature of the breeding season (TMIN), magnitude of the 25-year (19802005) trend in average total precipitation of the breeding season (PRECIP), edge density (ED),
percent developed land (DEVEL), interaction between the magnitude of the 25-year (1980-2005)
trend in average maximum temperature of the breeding season and edge density (TMAX*ED),
interaction between the magnitude of the 25-year (1980-2005) trend in average minimum
temperature of the breeding season and edge density (TMIN*ED), interaction between
magnitude of the 25-year (1980-2005) trend in average total precipitation of the breeding season
and edge density (PRECIP*ED), and survey effort (EFF). 50% indicates the mean of the
posterior distribution, while 2.5% and 97.5% are the lower and upper quantiles of the posterior
distribution.
Avian
Grouping

Model

TMAX*ED +
TMIN*ED +
PRECIP*ED +
DEVEL + EFF

All Species
TMAX + TMIN +
PRECIP + ED +
DEVEL + EFF

TMAX + TMIN +
PRECIP + EFF

Coefficient
βINT
βTMAX
βTMIN
βPRECIP
βED
βDEVEL
βTMAX-ED
βTMIN-ED
βPRECIP-ED
βEFF
β0
βTMAX
βTMIN
βPRECIP
βED
βDEVEL
βEFF
β0
βTMAX
βTMIN
βPRECIP
βEFF

150

Coefficient Estimates
50%
2.5%
97.5%
-1.337
-1.408
-1.276
0.020
0.006
0.035
-0.011
-0.025
0.003
-0.058
-0.076
-0.038
-0.005
-0.017
0.009
0.021
0.007
0.033
0.000
-0.011
0.011
0.011
0.000
0.022
0.011
-0.001
0.023
0.002
-0.006
0.010
-1.323
-1.422
-1.269
0.017
0.003
0.032
-0.011
-0.025
0.002
-0.059
-0.077
-0.040
-0.004
-0.019
0.012
0.022
0.009
0.035
0.002
-0.006
0.010
-1.320
-1.401
-1.260
0.017
0.000
0.031
-0.013
-0.026
0.001
-0.061
-0.077
-0.042
0.003
-0.005
0.010
Continued on next page

Table A 3.3 (cont’d)
Avian
Grouping

Model

TMAX*ED +
TMIN*ED +
PRECIP*ED + EFF
All Species

TMAX + TMIN +
PRECIP + ED + EFF

TMAX*ED +
TMIN*ED +
PRECIP*ED +
DEVEL + EFF
Long-distance
Migrants

TMAX + TMIN +
PRECIP + ED +
DEVEL + EFF

Coefficient
βINT
βTMAX
βTMIN
βPRECIP
βED
βTMAX-ED
βTMIN-ED
βPRECIP-ED
βEFF
β0
βTMAX
βTMIN
βPRECIP
βED
βEFF
βINT
βTMAX
βTMIN
βPRECIP
βED
βDEVEL
βTMAX-ED
βTMIN-ED
βPRECIP-ED
βEFF
β0
βTMAX
βTMIN
βPRECIP
βED
βDEVEL
βEFF

151

Coefficient Estimates
50%
2.5%
97.5%
-1.355
-1.412
-1.247
0.016
0.000
0.030
-0.013
-0.028
0.003
-0.056
-0.075
-0.037
-0.005
-0.019
0.010
0.001
-0.011
0.012
0.011
0.000
0.022
0.013
0.002
0.026
0.003
-0.006
0.011
-1.317
-1.364
-1.254
0.015
0.001
0.031
-0.012
-0.025
0.002
-0.060
-0.079
-0.041
-0.006
-0.019
0.007
0.002
-0.006
0.010
-1.410
-1.471
-1.296
0.046
0.026
0.068
-0.006
-0.027
0.016
-0.045
-0.078
-0.013
0.043
0.018
0.063
0.041
0.020
0.062
0.015
-0.002
0.034
0.008
-0.009
0.026
0.015
-0.004
0.034
-0.002
-0.016
0.011
-1.423
-1.503
-1.358
0.043
0.017
0.066
-0.007
-0.029
0.014
-0.056
-0.083
-0.023
0.042
0.022
0.065
0.042
0.020
0.064
-0.003
-0.016
0.011
Continued on next page

Table A 3.3 (cont’d)
Avian
Grouping

Model

TMAX*ED +
TMIN*ED +
PRECIP*ED + EFF

Long-distance
Migrants
TMAX + TMIN +
PRECIP + ED + EFF

TMAX + TMIN +
PRECIP + EFF

TMAX + TMIN +
PRECIP + EFF
Short-distance
Migrants
TMAX + TMIN +
PRECIP + ED +
DEVEL + EFF

Coefficient
βINT
βTMAX
βTMIN
βPRECIP
βED
βTMAX-ED
βTMIN-ED
βPRECIP-ED
βEFF
β0
βTMAX
βTMIN
βPRECIP
βED
βEFF
β0
βTMAX
βTMIN
βPRECIP
βEFF
β0
βTMAX
βTMIN
βPRECIP
βEFF
β0
βTMAX
βTMIN
βPRECIP
βED
βDEVEL
βEFF

152

Coefficient Estimates
50%
2.5%
97.5%
-1.432
-1.518
-1.369
0.042
0.019
0.067
-0.008
-0.028
0.012
-0.047
-0.078
-0.016
0.041
0.022
0.060
0.019
0.000
0.036
0.008
-0.008
0.026
0.018
-0.003
0.036
0.000
-0.013
0.011
-1.421
-1.491
-1.361
0.037
0.012
0.059
-0.010
-0.031
0.013
-0.057
-0.083
-0.028
0.040
0.018
0.064
-0.001
-0.014
0.011
-1.418
-1.484
-1.355
0.033
0.012
0.056
-0.010
-0.030
0.012
-0.047
-0.077
-0.019
-0.001
-0.015
0.012
-1.436
-1.507
-1.362
0.009
-0.019
0.040
-0.006
-0.034
0.019
-0.029
-0.066
0.009
0.014
-0.003
0.029
-1.429
-1.530
-1.361
0.014
-0.014
0.040
-0.006
-0.036
0.021
-0.028
-0.063
0.015
-0.003
-0.028
0.023
0.027
0.000
0.052
0.013
-0.003
0.028
Continued on next page

Table A 3.3 (cont’d)
Avian
Grouping

Model

TMAX*ED +
TMIN*ED +
PRECIP*ED +
DEVEL + EFF

Short-distance
Migrants

TMAX + TMIN +
PRECIP + ED + EFF

TMAX*ED +
TMIN*ED +
PRECIP*ED + EFF

TMAX + TMIN +
PRECIP + ED +
DEVEL + EFF
Resident
Species
TMAX + TMIN +
PRECIP + ED + EFF

Coefficient
βINT
βTMAX
βTMIN
βPRECIP
βED
βDEVEL
βTMAX-ED
βTMIN-ED
βPRECIP-ED
βEFF
β0
βTMAX
βTMIN
βPRECIP
βED
βEFF
βINT
βTMAX
βTMIN
βPRECIP
βED
βTMAX-ED
βTMIN-ED
βPRECIP-ED
βEFF
β0
βTMAX
βTMIN
βPRECIP
βED
βDEVEL
βEFF
β0
βTMAX
βTMIN
βPRECIP
βED
βEFF

153

Coefficient Estimates
50%
2.5%
97.5%
-1.457
-1.564
-1.383
0.016
-0.016
0.043
-0.006
-0.033
0.022
-0.027
-0.066
0.009
-0.003
-0.029
0.026
0.026
-0.004
0.052
0.000
-0.021
0.020
0.009
-0.012
0.031
0.007
-0.016
0.030
0.013
-0.002
0.028
-1.448
-1.583
-1.371
0.011
-0.018
0.045
-0.011
-0.036
0.015
-0.026
-0.064
0.009
-0.006
-0.031
0.019
0.014
-0.002
0.028
-1.475
-1.552
-1.396
0.013
-0.018
0.044
-0.010
-0.036
0.021
-0.028
-0.071
0.015
-0.002
-0.028
0.022
0.001
-0.021
0.023
0.013
-0.010
0.035
0.012
-0.011
0.034
0.014
-0.001
0.030
-0.989
-1.137
-0.829
0.000
-0.025
0.026
-0.014
-0.038
0.011
-0.090
-0.126
-0.050
-0.062
-0.085
-0.040
-0.016
-0.039
0.006
-0.006
-0.019
0.006
-0.976
-1.113
-0.864
0.000
-0.024
0.025
-0.012
-0.037
0.012
-0.091
-0.122
-0.057
-0.061
-0.084
-0.037
-0.007
-0.019
0.007
Continued on next page

Table A 3.3 (cont’d)
Avian
Grouping

Model

TMAX*ED +
TMIN*ED +
PRECIP*ED +
DEVEL + EFF

Resident
Species
TMAX*ED +
TMIN*ED +
PRECIP*ED + EFF

TMAX + TMIN +
PRECIP + EFF

Coefficient
βINT
βTMAX
βTMIN
βPRECIP
βED
βDEVEL
βTMAX-ED
βTMIN-ED
βPRECIP-ED
βEFF
βINT
βTMAX
βTMIN
βPRECIP
βED
βTMAX-ED
βTMIN-ED
βPRECIP-ED
βEFF
β0
βTMAX
βTMIN
βPRECIP
βEFF

154

Coefficient Estimates
50%
2.5%
-1.072
-1.200
-0.002
-0.025
-0.012
-0.035
-0.090
-0.127
-0.060
-0.086
-0.015
-0.038
-0.014
-0.033
0.013
-0.007
0.004
-0.014
-0.006
-0.020
-0.974
-1.079
-0.001
-0.024
-0.011
-0.033
-0.090
-0.121
-0.058
-0.080
-0.015
-0.034
0.013
-0.007
0.003
-0.017
-0.006
-0.021
-0.951
-1.170
0.006
-0.016
-0.011
-0.034
-0.099
-0.128
-0.007
-0.022

97.5%
-0.950
0.025
0.012
-0.060
-0.036
0.007
0.006
0.031
0.025
0.008
-0.819
0.023
0.012
-0.055
-0.035
0.005
0.031
0.022
0.007
-0.681
0.031
0.013
-0.070
0.007

Figure A 1.1 Spatial distribution of (A) magnitude of the 25-year (1980-2005) trend in average maximum temperatures [ºC], (B)
magnitude of the 25-year (1980-2005) trend in average minimum temperatures [ºC], (C) magnitude of the 25-year (1980-2005) trend
in total precipitation [mm], (D) percent developed land [%], (E) edge density [m/ha], and (F) survey effort [h].

155

Figure A 1.2 Deviations of the spatially-varying coefficient estimates from the global mean for (A) the intercept, (B) the magnitude of
the 25-year (1980-2005) trend in average maximum temperature of the breeding season (βTMAX), (C) the magnitude of the 25-year
(1980-2005) trend in average minimum temperature of the breeding season (βTMIN), (D) the magnitude of the 25-year (1980-2005)
trend rates in average total precipitation of the breeding season (βPRECIP), (E) percent developed land (βDEVEL), (F) edge density (βED),
and (G) effort (βEFF).

156

Figure A 2.1 The relationships between the amount of habitat (GPLAND and FPLAND for
grassland and forest birds, respectively; left-hand panel of each column) and spatially-varying
coefficient estimates (right-hand panel of each column) of 2000-05 average spring maximum
temperature (TMAX, left column), 2000-05 average spring minimum temperature (TMIN,
middle column), and 2000-05 average total spring precipitation (PRECIP, right column) for
[from top to bottom] American kestrel (AMKE), northern harrier (NOHA), eastern bluebird
(EABL), eastern kingbird (EAKI), horned lark (HOLA), grasshopper sparrow (GRSP), savannah
sparrow (SAVS), vesper sparrow (VESP), bobolink (BOBO), brown-headed cowbird (BHCO),
eastern meadowlark (EAME)], blue-gray gnatcatcher (BGGN), blue-headed vireo (BHVI),
golden-crowned kinglet (GCKI), great crested flycatcher (GCFL), red-breasted nuthatch
(RBNU), tufted titmouse (ETTI), veery (VEER), black-and-white warbler (BAWW), blackthroated blue warbler (BTBW), black-throated green warbler (BTNW), blackburnian warbler
(BLBW), and Canada warbler (CAWA). Dark blue symbols represent coefficient estimates with
95% credible intervals not overlapping zero in locations where species was detected; black
symbols represent coefficient estimates with 95% credible intervals not overlapping zero in
locations where species was not detected; light blue points represent with 95% credible intervals
overlapping zero coefficient estimates in locations where species was not detected; grey points
represent coefficient estimates with 95% credible intervals overlapping zero in locations where
the species was not detected. Black dashed line represents the slope of the relationship between
the coefficient estimates and GPLAND and FPLAND. The right-hand panel shows coefficient
estimates of TMAX, TMIN, and PRECIP resulting from the spatially-varying coefficients
models. Grayed-out locations indicate coefficient estimates with 95% credible intervals
overlapping zero. Species occurrence data were retrieved from the 2000-05 New York State
Breeding Bird Atlas.

Continued on next page
157

Figure A 2.1 (cont’d)

Continued on next page

158

Figure A 2.1 (cont’d)

Continued on next page

159

Figure A 2.1 (cont’d)

160

Figure A 3.1 Model validation using plots of predicted vs observed values for temporal turnover,
extinction, and colonization observed between 1980-85 and 2000-05 for communities of (A) all
recorded species regardless of their migratory status, (B) long-distance migrants, (C) shortdistance migrants, and (D) resident birds. Model validation was done for the best model only.

161

LITERATURE CITED

162

LITERATURE CITED

Anders A.D. and Post E. 2006. Distribution-wide effects of climate on population densities of a
declining migratory landbird. Journal of Animal Ecology, 75, 221-227.
Andrle R.F. and Carroll J.R. 1988. The Atlas of Breeding Birds in New York State. Cornell
University Press, Ithaka, NY.
Araújo M.B. and Peterson T.A. 2012. Uses and misuses of bioclimatic envelope modeling.
Ecology, 93(7), 1527-1539.
Arrhenius O. 1921. Species and area. Journal of Ecology, 9, 95-99.
Austin P.C. and Tu J.V. 2004. Bootstrap methids for developing predictive models. The
American Statistician, 58(2), 131-137.
Banerjee S., Carlin B.P. and Gelfand A.E. 2004. Hierarchical modeling and analysis for spatial
data. Chapman & Hall/CRC, Boca Raton, FL.
Barbet-Massin M. and Jetz W. In review. The effects of range changes on the functional
turnover, structure and diversity of bird assemblages under future climate scenarios.
Barnagaud J.-Y., Devictor V., Jiguet F., Barbet-Massin M., Le Viol I. and Archaux F. 2012.
Relating habitat and climatic niches in birds. PLoS ONE, 7(3), e32819.
Barton P.S., Cunningham S.A., Manning A.D., Gibb H., Lindenmayer D.B. and Didham R.K.
2013. The spatial scaling of beta diversity. Global Ecology and Biogeography, 22, 639-647.
Baselga A. 2010. Partitioning the turnover and nestedness components of beta diversity. Global
Ecology and Biogeography, 19, 134-143.
Belmaker J. and Jetz W. 2011. Cross-scale variation in species richness-environment
associations. Global Ecology and Biogeography, 20, 464-474.
Bini L.M., Diniz-Filho J.A.F., Rangel T.F.L.V.N., Akre T.S.B., Albaladejo R.G., Albuquerque
F.S., Aparicio A., Araújo M.B., Baselga A., Beck J., Bellocq M.I., Böhning-Gaese K., Borges
P.A.V., Castro-Parga I., Chey V.K., Chown S.L., De Marco P.Jr, Dobkin D.S., Ferrer-Castan D.,
Field, R., Filloy J., Fleishman E., Gomez J.F., Hortal J., Iverson J.B., Kerr J.T., Kissling W.D.,
Kitching I.J., Leon-Cortes J.L., Lobo J.M., Montoya D., Morales-Castilla I., Moreno J.C.,
Oberdoff T., Olalla-Tarraga M.A., Pausas J.G., Qian H., Rahbek C., Rodriguez M.A., Rueda M.,
Ruggiero A., Sackmann P., Sanders N.J., Terribile L.C., Vetaas O.R. and Hawkins B.A. 2009.
Coefficient shifts in geographical ecology: an empirical evaluation of spatial and non-spatial
regression. Ecography, 32, 193-204.

163

Bivand R. and Yu D. 2013. spgwr: Geographically weighted regression. R package version 0.624.
Bonan G.B. 2008. Forests and climate change: forcings, feedbacks, and the climate benefits of
forests. Science, 320(5882), 1444-1449.
Brennan L.A. and Kuvlesky W.P. 2005. North American grassland birds: an unfolding
conservation crisis? Journal of Wildlife Management, 69, 1-13.
Brotons L. and Jiguet F. 2010. Bird communities and climate change. In: Effects of climate
change on birds (eds: Moller A.P., Fielder W. and Berthold P.). Oxford University Press,
Oxford.
Buckley L.B. and Jetz W. 2008. Linking global turnover of species and environments.
Proceedings of the National Academy of Sciences of the United States of America, 105, 1783617841.
Carvalheiro L.G., Kunin W.E., Keil P., Aguirre-Gutiérrez J., Ellis W.N., Fox R., Groom Q.,
Hennekens S., van Landuyt W., Maes D., van de Meutter F., Michez D., Rasmont P., Ode B.,
Potts S.G., Reemer M., Roberts S.P.M., Schaminee J., Wallis de Vries M.F. and Biesmeijer J.C.
2013. Species richness declines and biotic homogenization have slowed down for NW-European
pollinators and plants. Ecology Letters, 16, 870-878.
Chao A., Cahzdon R.L., Colwell R.K. and Shen T. 2005. A new statistical approach for assessing
similarity of species composition with incidence and abundance data. Ecology Letters, 8, 148159.
Clavero M. and Brotons L. 2010. Functional homogenization of bird communities along habitat
gradients: accounting for niche multidimensionality. Global Ecology and Biogeography, 19,
684–696.
Clavero M., Villero D. and Brotons L. 2011. Climate change or land use dynamics: Do we know
what climate change indicators indicate? PLoS ONE, 6(4), e18581.
Cox W.A., Thompson III F.R., Reidy J. and Faaborg J. 2013. Temperature can interact with
landscape factors to affect songbird productivity. Global Change Biology, 19, 1064-1074.
Cressie N., Calder C.A, Clark J.S., ver Hoef J.M. and Wikle, C.K. 2009. Accounting for
uncertainty in ecological analysis: the strengths and limitations of hierarchical statistical
modeling. Ecological Applications, 19, 553-570.
Cressie N. and Wikle C.K. 2011. Statistics for Spatio-Temporal Data. John Wiley and Sons,
Hoboken, NJ.

164

Currie D.J. 2004. Regional-to-global patterns of biodiversity, and what they have to say about
mechanisms. In: Scaling Biodiversity (eds: Storch D., Marquet P.A. and Brown J.H.). Cambridge
University Press, Cambridge, UK.
Daly C. and Gibson W. 2002. Parameter-estimation on Independent Slopes Model (PRISM). The
PRISM Climate Group, Oregon.
Davies K.F. Chesson P., Harrison S., Inouye B.D., Melbourne B.A. and Rice K.J. 2005. Spatial
heterogeneity explains the scale dependence of the native-exotic diversity relationship. Ecology,
86, 1602-1610.
de Chazal J. and Rounsevell M.D.A. 2009. Land-use and climate change within assessments of
biodiversity change: A review. Global Environmental Change, 19, 306-315.
De Frenne P., Rodríguez-Sánchez F., Coomes D.A., Baeten L., Verstraeten G., Vellend M.,
Bernhardt-Romermann M., Brown C.D., Brunet J., Cornelus J., Decocq G.M., Dierschke H.,
Eriksson O., Gilliam F.S., Hedl R., Heinken T., Hermy M., Hommel P., Jenkins M.A., Kelly
D.L., Kirby K.J., Mitchell F.J.G., Naaf T., Newman M., Peterken G., Petrik P., Schultz J.,
Sonnier G., Van Calster H., Waller D.M., Walther G.-R., White P.S., Woods K.D., Wulf M.,
Graae B.J. and Verheyen K. 2013 Microclimate moderates plant responses to macroclimate
warming. Proceedings of the Royal Society B-Biological Sciences, 110(46), 18561–18565.
DeGraaf R.M. and Yamasak, M. 2001. New England Wildlife; Habitat, Natural History, and
distribution. University Press of New England, Hanover, NH.
DeLeon R.L., DeLeon E.E. and Rising G.R. 2011. Influence of climate change on avian
migrants' first arrival dates. Condor, 113, 915-923.
Diggle P.J., Tawn J.A. and Moyeed R.A. 1998. Model-based geostatistics. Journal of the Royal
Statistical Society Series C Applied Statistics, 47, 299-326.
Dobrowski S.Z., Abatzoglou J., Swanson A.K., Greenberg J.A., Mynsberge A.R., Holden Z.A.
and Schwartz M.K. 2013. The climate velocity of the contiguous United States during the 20th
century. Global Change Biology, 19, 241-251.
Dormann C.F., McPherson J.M., Araújo M.B., Bivand R., Bolliger J., Carl G., Davies R.G.,
Hirzel A., Jetz W., Kissling D.W., Kuhn I., Ohlemuller R., Peres-Neto P.R., Reineking B.,
Schrider B., Schurr F.M. and Wilson R. 2007. Methods to account for spatial autocorrelation in
the analysis of species distributional data: a review. Ecography, 30, 609-628.
Doxa A., Robert A., Crivelli A., Catsadorakis G., Naziridis T., Nikolaou H., Jiguet F. and
Theodorou K. 2012 Shifts in breeding phenology as a response to population size and climatic
change: A comparison between short- and long-distance migrant species. Auk, 129, 753-762.
Englund G. and Hamback P.A. 2007. Scale dependence of immigration rates: models, metrics
and data. Journal of Animal Ecology, 76, 30-35.
165

Estavillo C., Pardini R. and da Rocha P.L.B. 2013. Forest loss ad the biodiversity threshold: Am
evaluation considering species habitat requiremens and the use of matrix habitat. PLoS ONE,
8(12), e823369.
Field R., Hawkins B.A., Cornell H.V., Currie D.J., Diniz‐Filho J.A.F., Guégan J., Kaufman
D.M., Kerr J.T., Mittelbach G.G., Oberdorff T., O’Brien E.M. and Turner J.R.G. 2009. Spatial
species-richness gradients across scales: a meta-analysis. Journal of Biogeography, 36, 132–147.
Finley A.O., Banerjee S. and McRoberts R.E. 2009. Hierarchical spatial models for predicting
tree species assemblages across large domains. Annals of Applied Statistics, 3, 1052-1079.
Finley A.O. 2011. Comparing spatially-varying coefficients models for analysis of ecological
data with non-stationary and anisotropic residual dependence. Methods in Ecology and Evoluion,
2, 143-154.
Finley A.O., Banerjee S. and MacFarlane D.W. 2011. A hierarchical model for quantifying forest
variables over large heterogeneous landscapes with uncertain forest areas. Journal of the
American Statistical Association, 106, 31-48.
Finley A.O. and Banerjee S. 2013. spBayes: Univariate and multivariate spatial-temporal
modeling. R package version 0.3-7.
Foody G.M. 2004. Spatial nonstationarity and scale-dependency in the relationship between
species richness and environmental determinants for the sub-Saharan endemic avifauna. Global
Ecology and Biogeography, 13, 315-320.
Fortin M.J. and Dale M.R.T. 2009. Spatial autocorrelation in ecological studies: A legacy of
solutions and myths. Geographical Analysis, 41, 392-397.
Fotheringham A.S., Brunsdon C. and Charlton M. 2002. Geographically Weighted Regression:
The analysis of spatially varying relationships. John Wiley and Sons, Hoboken, NJ.
Forchhammer M.C., Post E. and Stenseth N.C. 2002. North Atlantic Oscillation timing of longand short-distance migration. Journal of Animal Ecology, 71, 1002-1014.
Fuller R.J., Gregory R.D., Gibbons D.W., Marchant J.H., Wilson J.D., Baillie S.R. and Carter N.
1995. Population declines and range contractions among lowland farmland birds in Britain.
Conservation Biology, 9(6), 1425-1441.
Gaston K.J. and Blackburn T.M. 2002. Large-scale dynamics in colonization and extinction for
breeding birds in Britain. Journal of Animal Ecology, 71, 390-399.
Gaston K.J., Evans K.L. and Lennon J.J. 2004. The scaling of spatial turnover: pruning the
thicket. In: Scaling Biodiversity (eds: Storch D., Marquet P.A. and Brown J.H.). Cambridge
University Press, Cambridge, UK.
166

Gaston K.J., Davies R.G., Orme C.D.L., Olson V.A., Thomas G.H., Ding T.-S., Rasmussen P.C.,
Lennon J.J., Bennett P.M., Owens I.P.F. and Blackburn T.M. 2007. Spatial turnover in the global
avifauna. Proceedings of the Royal Society B-Biological Sciences, 274, 1567-1574.
Geiger R., Aron R.H. and Todhunter P. 2009. The climate near the ground. Rowman &
Littlefield Publishing Group, Lanham, MD.
Gelfand A.E., Hyon-Jung K., Sirmans C.F. and Banerjee S. 2003. Spatial modeling with spatially
varying coefficient processes. Journal of the American Statistical Association, 98, 387-396.
Gelman A. and Rubin D.B. 1992. Inference from iterative simulation using multiple
sequences. Statistical Science, 7, 457–511.
Gelman A., Carlin J.B., Stern H.S. and Rubin D.B. 2004. Bayesian Data Analysis, 2nd edn.
Chapman & Hall/CRC, Boca Raton, FL.
Gibbons D.W., Donald P.F., Bauer H.-G., Fornasari L. and Dawson I.K. 2007. Mapping avian
distributions: the evolution of bird atlases. Bird Study, 54, 324-334.
Grotan V., Sæther B.E., Engen S., van Balen J.H., Perdeck A.C. and Visser, M.E. 2009. Spatial
and temporal variation in the relative contribution of density dependence, climate variation and
migration to fluctuations in the size of great tit populations. Journal of Animal Ecology, 78, 447459.
Grubb Jr T.C. and Pravasudov V.V. 1994. Tufted Titmouse (Baeolophus bicolor). In: The Birds
of North America Online (ed: Poole A.). Cornell Lab of Ornithology, Ithaca, NY. Retrieved from
the Birds of North America Online: http://bna.birds.cornell.edu/bna/species/086
Guisan A. and Thuiller W. 2005. Predicting species distribution: offering more than simple
habitat models. Ecology Letters, 8, 993-1009.
Gutiérrez-Cánovas C., Millan A., Velasco J., Vaughan I.P. and Ormerod S.J. 2013. Contrasting
effects of natural and anthropogenic stressors on beta diversity in river organisms. Global
Ecology and Biogeography, 22, 796-805.
Hargis C.D., Bissonette J.A. and David J.L. 1998. The behavior of landscape metrics commonly
used in the study of habitat fragmentation. Landscape Ecology, 13, 167-186.
Hartley S. and Kunin W.E. 2003. Scale dependency of rarity, extinction risk, and conservation
priority. Conservation Biology. 17, 1559-1570.
Hoeting J.A. 2009. The importance of accounting for spatial and temporal correlation in analyses
of ecological data. Ecological Applications, 19, 574-577.

167

Homer C., Huang C., Yang L., Wylie B. and Coan M. 2004. Development of a 2001 National
Landcover Database for the United States. Photogrammetric Engineering and Remote Sensing,
70(7), 829-840.
Huntley B., Barnard P., Altwegg R., Chambers L., Coetzee B.W.T., Gibson L., Hockey P.A.R.,
Hole D.G., Midgley G.F., Underhill L.G. and Willis S.G. 2010. Beyond bioclimatic envelopes:
dynamic species’ range and abundance modelling in the context of climatic change. Ecography,
33, 621-626.
Huste A. and Bouliner T. 2007. Determinants of local extinction and turnover rates in urban bird
communities. Ecological Applications, 17(1), 168-180.
Intergovernmental Panel on Climate Change. 2013. Working Group I Contribution to the IPCC
Fifth Assessment Report (AR5). Climate Change 2013: The Physical Science Basis. Cambridge
University Press, Cambridge, UK.
Jarzyna M.A., Zuckerberg B. and Porter W.F. 2013. Climate Change and Wildlife. In: Wildlife
Management and Conservation: Contemporary Principles and Practices (eds: Krausman P.R.
and Cain III J.W.). Johns Hopkins University Press, Baltimore, MD.
Jeltsch F., Moloney K.A., Schwager M., Körner K. and Blaum N. 2011. Consequences of
correlations between habitat modifications and negative impact of climate change for regional
species survival. Agriculture, Ecosystems and Environment, 145, 49-58.
Johnson C.N. 1998. Species extinction and the relationship between distribution and abundance.
Nature, 394, 272-274.
Jones T. and Cresswell W. 2010. The phenology mismatch hypothesis: are declines of migrant
birds linked to uneven global climate change? Journal of Animal Ecology, 79, 98-108.
Jost L. 2007. Partitioning diversity into independent alpha and beta components. Ecology, 88,
2427-2439.
Juliano S.A.1983. Body size, dispersal ability, and range size in North American species of
Brachinus (Coleoptera: Carabidae). The Coleopterists Bulletin, 37(3), 232-238.
Kampichler C., van Turnhout C.A.M., Devictor V. and van der Jeugd H.P. 2012. Large-scale
changes in community composition: determining land use and climate change signals. PLoS
ONE, 7(4), e35272.
Keil P., Biesmeijer J.C., Barendregt A., Reemer M. and Kunin W.E. 2010. Biodiversity change
is scale-dependent: an example from Dutch and UK hoverflies (Diptera, Syrphidae). Ecography,
34, 392-401.
Keil P. and Jetz W. 2014. Downscaling the environmental associations and spatial patterns of
species richness. Ecological Applications, 4, 82-94.
168

Kelling S., Hochachka W.M., Fink D., Riedewald M., Caruana R. and Hooker G. 2009. DataIntensive Science: A new paradigm for biodiversity studies. BioScience, 59, 613-620.
Kennedy C.M., Campbell Grant E.H., Neel M.C., Fagan W.F. and Marra P.P. 2011. Landscape
matrix mediates occupancy dynamics of Neotropical avian insectivores. Ecological Applications,
21, 1837-1850.
Kershner E.L. and Ellison W.G. 2012. Blue-gray Gnatcatcher (Polioptila caerulea). In: The
Birds of North America Online (ed: Poole A.). Cornell Lab of Ornithology, Ithaca, NY.
Retrieved from the Birds of North America Online: http://bna.birds.cornell.edu/bna/species/023
Kleijn D., Schekkerman H., Dimmers W.J., Van Kats R.J.M., Melman D. and Teunissen W.A.
2010. Adverse effects of agricultural intensification and climate change on breeding habitat
quality of Black-tailed Godwits Limosa l. limosa in the Netherlands. Ibis, 152(3), 475-486.
Kraft N.J.B., Comita L.S., Chase J.M., Sanders N.J., Swenson N.G., Crist T.O., Stegen J.C.,
Vellend M., Boyle B., Anderson M.J., Cornell H.V., Davies K.F., Freestone A.L., Inouye B.D.,
Harrison S.P. and Myers J.A. 2011. Disentangling the drivers of beta diversity along latitudinal
and elevational gradients. Science, 333, 1755-1758.
Krawchuk M.A. and Taylor P.D. 2003. Changing importance of habitat structure across multiple,
spatial scales for three species of insects. Oikos, 103, 153-161.
La Sorte F.A. and Thompson F.R. 2007. Poleward shifts in winter ranges of North American
birds. Ecology, 88, 1803-1812.
Lawler J.J. 2009. Climate change adaptation strategies for resource management and
conservation planning. Annals of the New York Academy of Sciences, 1162, 79-98.
Legendre P. 1993. Spatial autocorrelation – trouble or new paradigm. Ecology, 74, 1659-1673.
Lehikoinen E. and Sparks T.H. 2010. Changes in migration. In: Effects of climate change on
birds (eds: Moller A.P., Fielder W. and Berthold P.). Oxford University Press, Oxford.
Lemoine N. and Böhning-Gaese K. 2003. Potential impacts of global climate change on species
richness of long-distance migrants. Conservation Biology, 17(2), 577-586.
Lemoine N., Schaefer H.-C. and Böhning-Gaese K. 2007. Species richness of migratory birds is
influenced by global climate change. Global Ecology and Biogeography, 16, 55-64.
Lennon J.J., Koleff P., Greenwood J.J.D. and Gaston K.J. 2001. The geographical structure of
British bird distributions: diversity, spatial turnover and scale. Journal of Animal Ecology, 70,
966–979.

169

Lichstein J.W., Simons T.R., Shriner S.A., Franzreb K.E. 2002. Spatial autocorrelation and
autoregressive models in ecology. Ecological Monographs, 72, 445-463.
Ma Z., Zuckerberg B., Porter W.F. and Zhang L. 2012a. Spatial Poisson models for examining
the influence of climate and land cover pattern on bird species richness. Forest Science, 58, 6174.
Ma Z., Zuckerberg B., Porter W.F. and Zhang L. 2012b. Use of localized descriptive statistics
for exploring the spatial pattern changes of bird species richness at multiple spatial scales.
Applied Geography, 32, 185-194.
Magurran A.E. 2004. Measuring biological diversity. Blackwell Publishing, Oxford, UK.
Marini M.A., Barbet-Massin M., Lopes L.E. and Jiguet F. 2009. Predicted climate-driven bird
distribution changes and forecasted conservation conflicts in a neotropical savanna.
Conservation Biology, 23, 1558-1567.
Martin S.G. and Gavin T.A. 1995. Bobolink (Dolichonyx oryzivorus). In: The Birds of North
America Online (ed: Poole A.). Cornell Lab of Ornithology, Ithaca, NY. Retrieved from the
Birds of North America Online: http://bna.birds.cornell.edu/bna/species/176
Martin-Queller E., Gil-Tena A. and Saura S. 2011. Species richness of woody plants in the
landscapes of Central Spain: the role of management disturbances, environment and nonstationarity. Journal of Vegetation Science, 22, 238-250.
McGarigal K., Cushman S.A. and Ene, E. 2012. FRAGSTATS v4: Spatial Pattern Analysis
Program for Categorical and Continuous Maps. Computer software program produced by the
authors at the University of Massachusetts, Amherst.
McGill B.J. 2010. Matters of scale. Science, 328, 575-576.
McGowan K.J. and Zuckerberg B. 2008. Summary of results. In: The Second Atlas of Breeding
Birds in New York State (eds: McGowan K.J. and Corwin K.). Cornell University Press, Ithaka,
NY.
McGowan K.J. and Corwin K. 2008. The Second Atlas of Breeding Birds in New York State.
Cornell University Press, Ithaka, NY.
McNew L.B., Gregory A.J. and Sandercock B.K. 2013. Spatial heterogeneity in habitat selection:
Nest site selection by greater prairie-chickens. Journal of Wildlife Management, 77, 791-801.
Menendez R. and Thomas C.D. 2000. Metapopulation structure depends on spatial scale in the
host-specific moth Wheeleria spilodactylus (Lepidoptera: Pterophoridae). Journal of Animal
Ecology, 69, 935-951.
Miller J. 2010. Species distribution modeling. Geography Compass, 4, 490–509.
170

Miller J.A. 2012. Species distribution models: Spatial autocorrelation and non-stationarity.
Progress in Physical Geography, 36, 681-692.
Moss R., Oswald J. and Baines D. 2001. Climate change and breeding success: decline of the
capercaillie in Scotland. Journal of Animal Ecology, 70, 47-61.
Mouquet N. and Loreau M. 2003. Community patterns insource-sink metacommunities. The
American Naturalist, 162(5), 544-557.
Opdam P. and Wascher D. 2004. Climate change meets habitat fragmentation: linking landscape
and biogeographical scale levels in research and conservation. Biological Conservation, 117,
285-297.
Parmesan C. and Yohe G. 2003. A globally coherent fingerprint of climate change impacts
across natural systems. Nature, 421, 37-42.
Pearson R.G., Dawson T.P., Berry P.M. and Harrison P.A. 2002. SPECIES: A spatial evaluation
of climate impact on the envelope of species. Ecological Modelling, 154, 289-300.
Pearson R.G., Dawson T.P. and Liu C. 2004. Modelling species distributions in Britain: a
hierarchical integration of climate and land-cover data. Ecography, 27, 285-298.
Pickett S.T.A. 1989. Space-for-time substitution as an alternative to long-term studies. In: Longterm studies in ecology. Approaches and alternatives (ed: Likens G.E.). Springler-Verlag New
York, Inc., New York, NY.
Plummer M., Best N., Cowles K., Vines K., Sarkar D. and Almond R. 2012. coda: Output
analysis and diagnostics for MCMC. R package version 0.16-1.
Powell K.I., Chase J.M. and Knight T.M. 2013. Invasive plants have scale-dependent effects on
diversity by altering species-area relationships. Science, 339, 316-318.
Prevedello J.A. and Vieira M.V. 2010. Does the type of matrix matter? A quantitative review of
the evidence. Biodiversity Conservation, 19, 1205–1223.
Prugh L.R., Hodges K.E., Sinclair A.R.E. and Brashares J.S. 2008. Effect of habitat area and
isolation on fragmented animal populations. Proceeding of the National Academy of Sciences of
the United States of America, 105(52), 20770-20775.
Rahbek C. 2005. The role of spatial scale and the perception of large-scale species-richness
patterns. Ecology Letters, 8, 224-239.
Record S., Fitzpatrick M.C., Finley A.O., Veloz S. and Ellison A.M. 2013. Should species
distribution models account for spatial autocorrelation? A test of model projections across eight
millennia of climate change. Global Ecology and Biogeography, 22, 760-771.
171

Reichman O.J., Jones M.B. and Schildhauer M.P. 2011. Challenges and opportunities of open
data in ecology. Science, 331, 703-705.
Reif J. 2013. Long-term trends in bird populations: a review of patterns and potential drivers in
North America and Europe. Acta Ornithologica, 48(1), 1-16.
Ricklefs R.E. 1986. Community diversity – relative roles of local and regional processes.
Science, 235, 167-171.
Ricklefs R.E. 2004. A comprehensive framework for global patterns in biodiversity. Ecology
Letters, 7, 1-15.
Saino N., Rubolini D., Lehikoinen E., Sokolov L.V., Bonisolo-Alquati A., Ambrosini R.,
Boncoraglio G. and Moller A.P. 2009. Climate change effects on migration phenology may
mismatch brood parasitic cuckoos and their hosts. Biology Letters, 5, 539-541.
Sæther B.E., Lillegard M., Grotan V., Drever M.C., Engen S., Nudds T.D. and Podruzny K.M.
2008. Geographical gradients in the population dynamics of North American prairie ducks.
Journal of Animal Ecology, 77, 869-882.
Schimel D. 2011. The era of continental-scale ecology. Frontiers in Ecology and the
Environment, 9, 311.
Shigesada N. and Kawasaki K. 1997. Biological invasions: Theory and practice, Oxford
University Press, Oxford, UK.
Siefert A., Ravenscroft C., Weiser M.D. and Swenson N.G. 2013. Functional beta-diversity
patterns reveal deterministic community assembly processes in eastern North American trees.
Global Ecology and Biogeography, 22, 682-691.
Skagen S.K. and Adams A.A.Y. 2012. Weather effects on avian breeding performance and
implications of climate change. Ecological Application, 22(4), 1131-1145.
Smith K.G., Wittenberg S.R., Macwhirter R.B. and Bildstein K.L. 2011. Northern Harrier
(Circus cyaneus). In: The Birds of North America Online (ed: Poole A.). Cornell Lab of
Ornithology, Ithaca, NY. Retrieved from the Birds of North America Online:
http://bna.birds.cornell.edu/bna/species/602
Spiegelhalter D.J., Best N.G., Carlin B.P. and van der Linde A. 2002. Bayesian measures of
model complexity and fit (with discussion). Journal of the Royal Statistical Society, Series B, 64,
583–639.
Storch D., Marquet P.A. and Brown J.H. 2004. Scaling Biodiversity. Cambridge University
Press, Cambridge, UK.

172

Stralberg D., Jongsomjit D., Howell C.A., Snyder M.A., Alexander J.D., Wiens, J.A. and Root T.
2009. Re-shuffling of species with climate disruption: A no-analog future for California birds?
PLoS ONE, 4(9), e6825.
Swanson D.L., Ingold J.L. and Galati R. 2012. Golden-crowned Kinglet (Regulus satrapa). In:
The Birds of North America Online (ed: Poole A.). Cornell Lab of Ornithology, Ithaca, NY.
Retrieved from the Birds of North America Online: http://bna.birds.cornell.edu/bna/species/301
Thomas C.D., Bodsworth E.J., Wilson R.J., Simmons A.D., Davies Z.G., Musche M. and
Conradt L. 2001. Ecological and evolutionary processes at expanding range margins. Nature,
411(6837), 577–581.
Thomas C.D., Cameron A., Green R.E., Bakkenes M., Beaumont L.J., Collingham Y.C., Erasus
B.F.N., de Siqueira M.F., Grainger A., Hannah L., Hughes L., Huntley B., van Jaarsveld A.S.,
Midgley G.F., Miles L., Ortega-Huerta M.A., Peterson A.T., Phillips O.L. and Williams S.E.
2004. Extinction risk from climate change. Nature, 427, 145-148.
Tilman D. 2004. Niche tradeoffs, neutrality, and community structure: A stochastic theory of
resource competition, invasion, and community assembly. Proceedings of the National Academy
of Sciences of the USA, 101, 10854-10861.
Tobler M.W., Carrillo-Percastegui S.E., Leite Pitman R., Mares R. and Powell G. 2008. An
evaluation of camera traps for inventorying large- and medium-sized terrestrial rainforest
mammals. Animal Conservation, 11, 169-178.
Travis J.M.J. 2003. Climate change and habitat destruction: a deadly anthropogenic cocktail.
Proceedings of the Royal Society B-Biological Sciences, 270, 467–473.
Tscharntke T., Tylianakis J.M., Rand T.A., Didham R.K., Fahrig L., Batary P., Bengtsson J.,
Clough Y., Crist T.O., Dormann C.F., Ewers R.M., Frund J., Holt R.D., Holzschuh A., Klein
A.M., Kleijn D., Kremen C., Landis D.A., Laurance W., Lindenmayer D., Scherber C., Sodhi N.,
Steffan-Dewenter I., Thies C., van der Putten W.H. and Westphal C. 2012. Landscape
moderation of biodiversity patterns and processes – eight hypotheses. Biological Reviews, 87(3),
661-685.
Tuomisto H. 2010a. A diversity of beta diversities: straightening up a concept gone awry. Part 1.
Defining beta diversity as a function of alpha and gamma diversity. Ecography, 33, 2-22.
Tuomisto H. 2010b. A diversity of beta diversities: straightening up a concept gone awry. Part 2.
Quantifying beta diversity and related phenomena. Ecography, 33, 23-45.
Villegas J.C., Breshears D.D., Zou C.B. and Royer P.D. 2010. Seasonally pulsed heterogeneity
in microclimate: Phenology and cover effects along deciduous grassland–forest continuum.
Vadose Zone Journal, 9(3), 537-547.

173

Walther G.R., Post E., Convey P., Menzel A., Parmesan C., Beebee T.J.C., Fromentin J.-M.,
Hoegh-Guldberg O. and Bairlein F. 2002. Ecological responses to recent climate change. Nature,
416, 389-395.
Wheeler D.C. and Cadler C.A. 2007. An assessment of coefficient accuracy in linear regression
models with spatially varying coefficients. Journal of Geographical Systems, 9,145-166.
Wheeler D.C. and Waller L.A. 2009. Comparing spatially varying coefficient models: a case
study examining violent crime rates and their relationships to alcohol outlets and illegal drug
arrests. Journal of Geographical Systems, 11, 1-22.
Wiens J.A. 1989. Spatial scaling in ecology. Functional Ecology, 3, 385-397.
Willis K.J. and Bhagwat S.A. 2009. Biodiversity and climate change. Science, 326, 806-807.
Wilson R.J., Thomas C.D., Fox R., Roy D.B. and Kunin W.E. 2004. Spatial patterns in species
distributions reveal biodiversity change. Nature, 432, 393-396.
Zuckerberg B., Woods A. and Porter W.F. 2009. Poleward shifts in breeding bird distributions in
New York State. Global Change Biology, 15, 1866-1883.
Zuckerberg B. and Porter W.F. 2010. Thresholds in the long-term responses of breeding birds to
forest cover and fragmentation. Biological Conservation, 143, 952-962.

174