EVALUATION OF STATISTICAL METHODS FOR ESTIMATING LONG-TERM POPULATION CHANGE FROM EXTENSIVE WILDLIFE SURVEYS by Len Thomas B.Sc, The University of Sheffield, U.K., 1990 M.Sc, The University of York, U.K., 1991 A THESIS SUBMITTED IN PARTIAL FULFILLMENT OF THE REQUIREMENTS FOR THE DEGREE OF Doctor of Philosophy in THE FACULTY OF GRADUATE STUDIES Department of Forestry and Centre for Applied Conservation Biology We accept this thesis as conforming to the required standard THE UNIVERSITY OF BRITISH COLUMBIA September 1997 ©Len Thomas, 1997 In presenting this degree at the thesis in University of partial fulfilment of of department this or thesis for by his or requirements British Columbia, I agree that the freely available for reference and study. I further copying the representatives. an advanced Library shall make it agree that permission for extensive scholarly purposes may be her for It is granted by the understood that head of copying my or publication of this thesis for financial gain shall not be allowed without my written permission. . , Department of The University of British Columbia Vancouver, Canada DE-6 (2/88) Abstract Monitoring long-term changes in population abundance is an integral part of conservation-oriented research and management. Many extensive monitoring programs are based on annual counts at a large number of permanent survey sites. However, estimating population change from the resulting data is challenging and controversial. At least eight statistical methods have been proposed and used, and in this thesis I evaluated these methods using data from the North American Breeding Bird Survey (BBS). I began by demonstrating the importance of analysis method for the inferences drawn about species trends. I compared the results when three methods of trend analysis were applied to BBS data for 115 species in British Columbia, and found differences in the direction, magnitude and statistical significance of the estimates. I then explored the issues surrounding the analysis of these data. I showed that even under ideal survey conditions subjective decisions must be made about whether to estimate annual indices or population trends and, for trend methods, which model to use for trend. In reality the analysis is further complicated by measurement error and missing data. I suggested the use of simulation studies to determine the accuracy of the methods, and provided a simple example using three geographic scenarios of population decline. I next made a preliminary assessment of the reliability of eight candidate analysis methods by applying them to continental data for four case study species (Cerulean Warbler, Carolina Wren, Chestnut-collared Longspur, and Belted Kingfisher). A l l methods except one (Mountford moving windows) produced seemingly plausible results in most cases. However each method demonstrated obvious flaws in some cases, and I used these to suggest potential areas of further methodological development. I concluded that loglinear modeling is currently the best method for estimating annual indices. For trend analysis, I could not determine whether route regression or direct estimation (based on generalized modeling) is generally superior. Large, long-term population changes in well-surveyed species are almost certain to be detected by any of these methods. However, their accuracy and robustness remains uncertain, and further research is required both to develop the methods and to evaluate them through comprehensive simulation studies. Table of Contents Abstract ii List of Tables vi List of Figures viii List of Abbreviations x Acknowledgements xi Foreword xii Chapter 1 General Introduction Introduction 1 The Breeding Bird Survey 2 Methods of Estimating Long-term Population Change 3 Chain Index Methods 3 Linear Route Regression 5 Nonlinear Route Regression 6 Linear Modeling 7 Overview of Thesis 8 Figures Chapter 2 The Importance of Analysis Method for Population Trend Estimates Introduction 14 Methods 15 Breeding Bird Survey , 15 Data Selection 15 Trend Estimation 16 Comparison of Trend Estimates 19 Results 20 Discussion 22 Causes of the Differences 22 Implications of the Results 25 Priorities for Further Work 26 Summary 27 Tables 29 Figures - iii - Chapter 3 Issues in the Estimation of Long-term Population Change Introduction 36 Analysis of Survey Data Under Ideal Conditions 36 Selection of the Unit of Analysis 37 Model of Trend 37 Additional Factors the Complicate the Analysis of Survey Data 39 Incomplete Sampling Frame 39 Measurement Error 40 Missing Data 41 Evaluating the Methods Using Monte Carlo Simulations 42 Summary 46 Tables 47 Figures Chapter 4 Four Case Studies Introduction 55 Methods 57 Data Selection 57 Data Analysis 58 Comparison of Estimates 62 Results 63 Cerulean Warbler 63 Carolina Wren 68 Chestnut-collared Longspur 71 Belted Kingfisher 73 Summary of Performance of Methods 75 Discussion 76 Linear Route Regression 76 Nonlinear Route Regression 79 Loglinear Modeling 81 Moving-window Mountford's Method 86 Generality of Findings 88 Summary 88 Tables 90 Figures - IV - Chapter 5 General Conclusions Introduction 119 Evaluation of Analysis Methods 119 Annual Index Methods 119 Trend Analysis Methods 121 Future Research 123 Development of Current Methods 123 Further Evaluation of Methods 124 Alternative Approaches 126 Implications of Research 128 Tables 130 Figure 133 Literature Cited 135 Appendix 1 Computer Programs For Analysis of Population Change 145 Appendix 2 Breeding Bird Survey Data 147 -v- List of Tables 2.1 Summary of the three Breeding Bird Survey trend estimation methods compared in Chapter 2 2.2 Pairwise comparisons among three Breeding Bird Survey trend estimation methods of the proportion of 115 species with trends in different directions or with different levels of statistical significance 2.3 29 31 Estimates of the effect of sample size and species abundance on the proportion of differences in the direction and statistical significance of species trends between three Breeding Bird Survey analysis methods 32 3.1 Methods that have been used to estimate annual indices from bird count data 47 3.2 Methods that have been used to estimate population trends from bird count data 48 3.3 Parameters used in the simulation study of the performance of three Breeding Bird Survey trend-estimation methods in the context of overall population declines and under three different geographic scenarios of population change 49 3.4 Bias and precision of two trend estimators (U.S. National Biological Service and Canadian Wildlife Service linear route-regression) from 100 simulations of Breeding Bird Survey data at each of three geographic scenarios of population decline 50 4.1 List of analysis methods used in Chapter 4, with their associated data selection criteria 90 4.2 Summary of the Breeding Bird Survey data used in this chapter 91 4.3 Estimates of overall population trend for the Cerulean Warbler between 1966 and 1995 (data selected using stringent data selection criteria) 92 Proportion of Breeding Bird Survey routes with estimated declines in abundance for four species using six route-regression trend estimators 93 4.4 4.5 Comparison of loglinear models for the analysis of population trends in Cerulean Warbler using Breeding Bird Survey data from 1966-1995 (data selected using stringent data selection criteria) 94 4.6 Goodness-of-fit tests for Mountford moving windows applied to Breeding Bird Survey data for Cerulean Warbler 1966-1995 4.7 4.8 95 Estimates of overall population trend for the Carolina Wren between 1967 and 1995 (data selected using stringent data selection criteria) 96 Comparison of loglinear models for the analysis of population trends in the Carolina Wren using Breeding Bird Survey data from 1967-1995 97 - vi - 4.9 Estimates of overall population trend for the Chestnut-collared Longspur between 1967 and 1995 (data selected using stringent data selection criteria) 98 4.10 Estimates of overall population trend for the Belted Kingfisher between 1968 and 1995 (data selected using stringent data selection criteria) 99 4.11 Comparison of loglinear models for the analysis of population trends in the Belted Kingfisher using Breeding Bird Survey data from 1968-1995 100 4.12 Summary of problems noted in case studies 101 4.13 Effect of incorporating observer differences on trend estimates for four species from three analyses 102 5.1 Summary evaluation of methods that have been used to estimate annual indices from bird count data 130 5.2 Summary evaluation of methods that have been used to estimate trends from bird count data 131 A l . 1 List of computer programs used in this thesis A2.1 Number of Breeding Bird Survey routes and surveys (1966 - 1995) in physiographic strata excluded from the case studies of Chapter 4 - vii - 146 148 List of Figures 1.1 Map showing the location of all North American Breeding Bird Survey routes that have been surveyed one or more times 11 1.2 Contour map showing the number of annual surveys on Breeding Bird Survey routes between 1966 and 1995 12 1.3 Number of Breeding Bird Survey surveys performed each year between 1966 and 1995. ...13 2.1 Number of species with increasing and decreasing trends from three analysis methods using the same Breeding Bird Survey data for 115 species from British Columbia 34 2.2 Plot of the magnitude of trends estimated from U.S. National Biological Service route regression against the magnitude of trends from Canadian Wildlife Service route regression for 115 species 35 3.1 Four models of trend applied to the same hypothetical count data 52 3.2 Three geographic scenarios of population decline 53 3.3 The direction and statistical significance of trends estimated by three methods of trend analysis from 100 simulations of Breeding Bird Survey data 54 4.1 Counts of Cerulean Warblers on each Breeding Bird Survey route between 1966 and 1995 105 4.2 Relationship between measures of variation in counts on a Breeding Bird Survey route and the mean count on that route for four species 106 4.3 Population indices for Cerulean Warblers between 1966 and 1995, calculated using eight different estimators 107 Population indices for the Cerulean Warbler between 1966 and 1995 calculated using model 3 of program T R I M 108 Mountford moving-window estimates of population indices for Cerulean Warblers between 1966 and 1996 109 Population indices for the Cerulean Warbler between 1966 and 1995, calculated using eight different estimators, with data selected using tolerant data selection criteria 110 Population indices for the Carolina Wren between 1967 and 1995, calculated using eight different estimators, with data selected using stringent data selection criteria Ill 4.4 4.5 4.6 4.7 4.8 Mountford moving-window estimates of population indices for Carolina Wren between 1966 and 1996 for "without observers" models at four window sizes 112 - viii - 4.9 Population indices for the Chestnut-collared Longspur between 1967 and 1995, calculated using eight different estimators, with data selected using stringent data selection criteria 113 4.10 Population indices for the Chestnut-collared Longspur between 1967 and 1995, calculated using seven different estimators (program T R I M failed to produce log-linear model estimates), with data selected using tolerant data selection criteria 114 4.11 Population indices for the Belted Kingfisher between 1968 and 1995, calculated using eight different estimators, with data selected using stringent data selection criteria 115 4.12 Population indices for the Belted Kingfisher between 1968 and 1995, calculated using eight different estimators, with data selected using tolerant data selection criteria 116 4.13 Relationship between mean abundance and estimated route trend for the linear route regression estimators using tolerantly selected data for Belted Kingfisher between 1968 and 1995 117 4.14 Residual analysis from loglinear models 118 5.1 134 A flow-chart guide to choosing among the analysis methods - ix - List of Abbreviations AIC Akaike's Information Criterion BBS Breeding Bird Survey CI Confidence interval CWS Canadian Wildlife Service GAM Generalized additive modeling GLM Generalized linear modeling NNRR Nonlinear nonparametric route regression NSRR Nonlinear semiparametric route regression QAIC Quasi-likelihood Akaike's Information Criterion RR Route regression SE Standard error TRIM Trends and indices for monitoring data (a computer program for loglinear modeling) USNBS U.S. National Biological Service - X - Acknowledgements I thank my research supervisor, Kathy Martin, for her support and encouragement throughout this thesis work, and the many other side projects in which I became entangled over the past six years. I also gratefully acknowledge the help of my supervisory committee, Val LeMay, Jamie Smith, John Smith and Carl Walters, each of whom passed on valuable nuggets of wisdom at various times. Many others helped to make the University of British Columbia (UBC) a lively and stimulating place to study, especially members of the Centre for Applied Conservation Biology and the Centre for Biodiversity Research. This work would not have been possible without the vision of Chandler Robbins, who originated the Breeding Bird Survey (BBS), and the efforts of the thousands of skilled volunteer observers who collect BBS data each year. In addition, various members of the two federal agencies that co-ordinate the BBS provided support and advice: from the Canadian Wildlife Service I thank Peter Blancher, Brian Collins, Ricky Dunn, Connie Downes and Bob Elner, and from the Biological Resources Division, U.S. Geological Survey, I thank Sam Droege, B i l l Link, Bruce Peterjohn and John Sauer. Many of the computer programs I used to analyze B B S data were originally written by others (see Appendix 1), and I am grateful to these people for sharing their code. Financial support for this work came from the Canadian Commonwealth Scholarship and Fellowship Plan, the Natural Sciences and Engineering Research Council of Canada and the Canadian Wildlife Service (Pacific and Yukon Region and National Wildlife Research Centre, Ottawa). I am grateful for the helpful suggestions made by those who reviewed earlier versions of Chapters 2 and 3 (see Foreword): Brian Collins, Connie Downes, Erica Dunn, John Hagan and John Sauer for Chapter 2, and Brian Collins, Peter Kareiva, David Ward and two anonymous referees for Chapter 3. M y ideas for Chapters 4 and 5 were improved by discussions with Steve Buckland, Paul Gustafson, Nancy Heckman and Bent Jorgensen. Tina Schivatcheva, Susan Shirkoff and Pierre Vernier helped me produce Figures 1.1 and 1.2 from GIS data provided by John Sauer. Although much of this work was done in various offices and cubbyholes at U B C , some of the best insights came while I was away in the mountains. I am therefore indebted to the good friends who dragged me out on many backcountry adventures over the years - you know who you are! Three very special friends must, however, be named. Rachel Holt, your contribution is too important to discuss without a glass of single malt in hand. Wes Hochachka, hi, and thanks for all the hoochery. Last, and most important, I dedicate this thesis to Katrine Vogt. - xi - Foreword Parts of this thesis have previously been published. A version of Chapter 2 was published as Thomas, L., and K. Martin. 1996. The importance of analysis method for Breeding Bird Survey population trend estimates. Conservation Biology 10: 479-490. and a version of chapter 3 was published as Thomas, L. 1996. Monitoring long-term population change: why are there so many analysis methods? Ecology 77: 49-58. The first paper was co-authored with my research supervisor, Dr. Kathy Martin, who played an important part in the early development of the work and provided substantial editorial input. I initiated the work, performed the analysis and drafted the manuscript. I agree that the contributions of the candidate to our co-authored paper are as stated. (Dr. K . Martin) - xii - Chapter 1 General Introduction Introduction Population monitoring can play a critical role in conservation biology, by providing the information necessary to identify conservation problems at an early stage and suggest possible solutions (Hellawell 1991; James and McCulloch 1995). However, to be effective, population parameters must be measured with sufficient accuracy to determine whether biologically important changes have occurred. Many regional, national and international population monitoring programs are based on a large number of permanent sampling sites that are surveyed at regular intervals using standardized counting techniques. Unfortunately, estimating long-term population changes from the resulting data is not straightforward. A number of statistical methods have been proposed and used but there is currently no consensus as to their relative merits (Geissler and Noon 1981; Sauer and Droege 1990b; Greenwood et al. 1993; ter Braak et al. 1994; Peterjohn et al. 1995; James et al. 1996). In this thesis I compare the various methods and evaluate them based on theoretical considerations and on their performance when applied to both real and simulated data. I focus on the estimation of change in abundance over time for a single species within a defined sampling area; I do not consider summary analyses that involve groups of species (James et al. 1996; Link and Sauer 1995, 1996) or spatial statistics (Flather and Sauer 1996; Villard and Maurer 1996). The data that I use come from the North American Breeding Bird Survey (BBS), which I describe in the next section. I then briefly introduce the statistical methods that have been used to estimate population change from these data. Lastly, I outline the contents of the rest of this dissertation. The Breeding Bird Survey The BBS is the most important source of information about the status of landbird populations in North America. The survey is operated cooperatively by U.S. National Biological Service (USNBS) and the Canadian Wildlife Service (CWS), and currently 1 comprises about 3700 39.2km long permanent survey routes located along secondary roads within all continental states in the U.S. and provinces and territories in Canada (Fig. l . l ) . 2 Routes were originally selected according to a stratified random sampling design, with a predetermined number of routes located in each degree block of latitude and longitude (approx. 112km by 80km; see Robbins and van Velzen 1970 for details of the survey design). Some routes (approx. 200) have since become unusable and have been replaced by newly selected routes. Additional routes have also been added in some states at various times (B. G. Peterjohn, USNBS Patuxent Environmental Science Center, pers. comm.). The density of routes is relatively uniform within states and the southern parts of the Canadian provinces, but varies between states and provinces from 8 routes per degree block in some eastern states to 0.5 routes per degree block in most western states and provinces (Fig 1.1). Analyses of BBS data therefore often post-stratify by state and province, and many also poststratify by ecologically similar physiographic strata (Fig 1.1; Bystrak 1981). Surveys are carried out between late May and early July, depending upon latitude. Starting half an hour before local sunrise and in reasonable weather, a skilled volunteer observer drives the route once, stopping every 0.8 km and recording all the birds seen within a 400 meter radius or heard at any distance in a three minute period (see Robbins et al. 1986 Appendix B for a detailed description of the survey protocol). For each species, the number of birds recorded is summed over the 50 stops to give the survey count on the route in that year. Approximately 70% of the 3700 active routes are surveyed in each year, although survey 1 The U.S. federal agency that administers the BBS has changed name (and sometimes mandate) several times since the survey began (Bureau of Sport Fisheries and Wildlife 1966-1974; Fish and Wildlife Service 1974-1993; National Biological Survey 1993-1995; National Biological Service 1995-1996; Biological Resources Division, Geological Survey 1996-present). Throughout this dissertation I will use the name by which it is currently most commonly known: the U.S. National Biological Service (USNBS). 2 All figures and tables appear at the end of the chapter in which they are referenced. coverage varies both spatially and temporally. Routes in some eastern and central states have been surveyed in almost every year, while those in sparsely populated parts of the western states and northern provinces and territories have rarely been visited (Fig. 1.2). Participation in the survey has increased since its inception in 1966 (Fig 1.3). While there are a number of limitations inherent in the survey methodology (Chapter 2; Robbins and van Velzen 1970; Bystrak 1981; Robbins et al. 1986; Droege 1990; Peterjohn et al. 1995) , the BBS is the only systematic surveillance program that covers the breeding ranges of most North American landbird species. Summaries containing species trends have been published regularly by the USNBS and CWS (e.g., Peterjohn et al. 1994; Downes and Collins 1996) , as well as more detailed accounts for some species (e.g., Sauer and Droege 1990a, d; Sauer et al. 1994a; Peterjohn and Sauer 1995) and analyses of multispecies patterns (e.g., Robbins et al. 1989; Sauer and Droege 1992). BBS data and trend analyses have been used by third parties for conservation prioritization schemes (e.g., Hunter et al. 1993), for population status reviews (e.g., Terborgh 1989; Askins et al. 1990; Askins 1993), for tests of hypotheses of causation (e.g., Bohning-Gaese et al. 1993; Johnson and Schwartz 1993; Heckert 1997) and as a benchmark in the validation of other surveillance programs (e.g., Holmes & Sherry 1988; Hagan et al. 1992; Hussell et al. 1992; Witham & Hunter 1992). Morton (1992: 585) summed up the general feeling among ornithologists by stating: "Population trend determinations by the BBS have become the accepted measure of North American breeding bird populations." Methods of estimating long-term population change Chain index methods Missing data are a ubiquitous feature of long-term volunteer-based surveillance programs like the BBS. Because of this, simple averages of the counts in each year are at best inefficient, and will be biased if there is an association between the missing data and bird abundance (Chapter 3; Geissler and Noon 1981; James et al. 1990; ter Braak et al. 1994). The earliest solution to this problem was the use of the chain index method (Bailey 1967; Ogilvie 1967), also known as the ratio method or proportional base-year adjustment. For each pair of adjacent years, counts from the routes that were surveyed in both years are summed and divided by one another to yield a ratio that estimates the change in abundance over that 2-year period. Successive ratios are "chained" (multiplied) together to produce indices of abundance relative to some "base-year" (arbitrary starting point). Variations on this method have been widely used (e.g., Erskine 1978; Marchant et al. 1990; O'Connor 1992a) but are no longer recommended for two reasons. Firstly, the chaining of successive indices can lead to an accumulation of errors, which produces very imprecise index estimates (Geissler and Noon 1981; Moss 1985). Secondly, ratio estimates can be positively biased (Barker and Sauer 1992; Link et al. 1994; Underhill and Prys-Jones 1994). A n improvement to the chain index method was suggested by Mountford (1982, 1985). This method uses not only ratios of counts in adjacent years, but also those from non-adjacent years to fit the model count = route effect * year effect + error where the error terms are random variables with means of zero, variances proportional to the expected count, and decreasing covariances between successive error terms on a route with increasing separation in time. The model is based on the assumption that yearly population changes are constant across all routes (no interaction effect), and this can be checked using a goodness-of-fit test (Mountford 1982). Unfortunately the method appears to be very sensitive to this assumption, which means that only relatively short time periods can be used in each analysis (6-12 years Peach and Baillie 1994; ter Braak et al. 1994; Peterjohn et al. 1995), making it of little use in the analysis of long-term surveys such as the BBS (30 years). Peach and Baillie (1994) have proposed a "moving windows Mountford method", in which successive overlapping Mountford models are fit to the data, each of which is of a time period short enough to pass the goodness-of-fit test. These successive models are then in effect chained together to produce indices for the entire time period of interest. This ad hoc solution has yet to be rigorously evaluated. Linear route regression The development of linear route regression was begun in response to the problems with chain indexing (Geissler and Noon 1981). Route regression (RR) methods are two-stage estimators in which trends are first estimated for each route and then the route trends are combined into an overall estimate for the region of interest (this is more generally known as multilevel or heirarchical modeling, Goldstein 1995). Linear R R is motivated by the exponential growth model expected count = route effect * route trend (1) year Most published analyses fit this model using linear least-squares regression on log-transformed data, using the equation log (count + constant) = log (observer effect) + log (route trend) * year + error where the error terms are independent random variables with identical Gaussian distributions and zero means (Geissler and Sauer 1990). The inclusion of separate intercepts for each observer on a route allows differences in detection efficiency between observers to be taken into account; in the chain index methods the only way to do this is to treat each observer as a different route. Regional trends are calculated as the weighted mean of the route trends, using weightings for the area represented by the route, the relative abundance of the species on the route, and the expected precision of the trend estimate. The merit of the weightings, and the appropriate method of calculating them (especially the abundance weighting), is controversial (James et al. 1990; Peterjohn et al. 1994). Variations of this method have been used by the U S N B S and CWS analyze BBS One to data since 1986 (e.g., Peterjohn and Sauer 1993; Downes and Collins 1996). important disadvantage of using least-squares regression on log-transformed counts is that a constant must be added to the counts before transformation, because the log of a zero count is undefined. This biases the trend estimates towards zero (i.e., no change, Geissler and Link 1988; Collins 1990; Link and Sauer 1994). Another problem is that trends are estimated on the log scale and back-transformation to the arithmetic scale is problematic (Geissler and Link 1988). To bypass these problems Link and Sauer (1994) recently proposed a method based on estimating equations (Godambe 1991) that directly fits the exponential growth model (equation (1), above). No transformation is performed, so back-transformation is not a problem, and no constant need be added. In addition the error distribution is not specified, so the method is very robust (assuming that the model is correct). This method of estimating route trends was used by the USNBS in their most recent published report (Peterjohn et al. 1994), but appropriate weightings for combining route trends based on this estimator have yet to be published. Nonlinear route regression A recurring criticism of linear R R methods is that they assume trends stay the same throughout the period analyzed, while in reality we would expect them to change over time. This problem is addressed within the linear R R framework by estimating trends for subsets of the entire time period and/or by calculating "residual" annual indices about the trend line (e.g., Peterjohn et al. 1995). Another solution, first proposed by James et al. (1990; see also Taub 1990) is to fit non-linear smoothers (Hastie and Tibshirani 1990 Chapter 2) to the route counts, and then average the smoothed route indices to produce regional trend estimates. In the James et al. method (nonlinear nonparametric route regression, NNRR), the loess smoother is used to smooth indices at both the route and physiographic stratum levels. A recent extension to the method (nonlinear semiparametric route regression, NSRR, James et al. 1996) allows differences between observers to be taken into account. One disadvantage of this approach is that it requires the specification of the level of smoothing required, and this is set subjectively (James et al. 1990). Oversmoothing may obscure abrupt but biologically important changes in abundance (Peterjohn et al. 1995), while undersmoothing leaves too much noise in the indices (James et al. 1990). The choice of weightings to use when combining route trends is again somewhat controversial, although there is no need for an abundance weighting with this method. A second nonlinear R R method, rank-trend analysis, was first applied to multi-site count data by Titus (1990; see also Titus et al. 1990). Counts on each route are ranked, a rank-trend statistic is calculated, and these route statistics are combined to give an overall statistic that indicates whether the counts show a tendency to increase or decline through time. Differences among observers are not accounted for, and no weighting scheme has been developed. The method makes no assumptions about the distribution of counts (except that the distribution is the same for all counts), but is relatively uninformative because it gives no information about the magnitude of the estimated population change. Linear modeling In contrast to R R methods, which use a two-stage (hierarchical or mixed-model) estimation process, linear modeling can be used to estimate the parameters directly from the count data. Sauer and Geissler (1990) used least-squares regression to fit a complex model containing observer, route-specific trend and year effects to log-transformed counts. As with least-squares linear R R a constant must be added to the counts before log-transformation and this biases the estimates. This problem is avoided by ter Braak et al. (1994), who proposed the use of loglinear modeling (Poisson regression), a type of generalized linear modeling ( G L M , McCullagh and Nelder 1989) in which the counts are assumed to be Poisson random variables and the logarithm of the expected count is modeled as a linear function of the parameters of interest. The models are fit using the method of maximum likelihood, which enables the use of a suite of well-established techniques for data-based model selection, model checking and parameter estimation in G L M s (McCullagh and Nelder 1989; Burnham and Anderson 1992). Fitting linear and loglinear models using standard computer software can be very computer memory intensive, especially for large datasets like the BBS which often require many parameters to be modeled adequately (Sauer and Geissler 1990). Recently, Pannekoek and van Strien (1996) have developed PC software (program TRIM) that enables large models to be fit and has modest memory requirements. This program can also be used to model serial correlation (autocorrelation) between counts on a site and overdispersion in counts relative to the Poisson distribution. The program can fit both annual index models expected count = route effect * year effect and exponential trend models expected count = route effect * trend year and the year effect and trend parameters can depend upon categorical covariables (such as region, habitat type, etc.). Another method that has been used to fit linear and loglinear models is impution of missing counts. Both applications of imputing to BBS data (Moses and Rabinowitz 1990; Bohning-Gaese et al. 1993) used linear models and would be improved by the use of loglinear models (e.g., Underhill and Prys-Jones 1994), so avoiding the need to transform the data and backtransform the estimates. However, imputing is not generally recommended because the estimation algorithm (expectation-maximization) is notoriously slow to converge when there are many missing data and identical results can be obtained using loglinear Poisson regression (ter Braaketal. 1994). Overview of thesis In this chapter I have introduced the North American Breeding Bird Survey and the statistical methods that I will be considering in subsequent chapters. Before going on to examine the methods in detail, however, it is important to determine whether differences in analysis method have an important effect on the inferences drawn about population trends. I do this in Chapter 2, by applying three analysis methods to BBS data from 115 species in British Columbia and quantifying the differences in trend estimates. I find that even relatively small differences in methods of data selection and analysis can produce important discrepancies in the results. I discuss possible causes of the observed differences, and their implications for users of B B S trend information. In Chapter 3 I explore the issues involved in the estimation of long-term population change from extensive survey data. I consider first how the data could be analyzed in an idealized survey, and then discuss the effect of various practical constraints on the analysis. I then suggest the use of Monte Carlo simulations to empirically test the performance of the methods under plausible scenarios of population change, and provide an example of the approach using the same analysis methods as in the previous chapter. Case studies are a useful way of examining the performance of the analysis methods in detail. In Chapter 4 I present studies of four species, using eight of the most promising analysis methods applied to continent-wide BBS data. The four study species are chosen to represent extremes with respect to two important factors that influence the reliability of the analysis methods: abundance (i.e., average count per route) and sample size (number of routes containing sightings). I place the results in the context of previous work, consider what inferences can be made about the reliability of the methods, and discuss possible improvements to the methods. In the final chapter, I bring this work together by giving a summary evaluation of the methods. I then recommend some future directions for research and suggest three promising new approaches to trend analysis. Lastly, I discuss the implications of my findings for the use of extensive wildlife surveys like the BBS in estimating long-term population change. Figure Legends Figure 1.1. Map showing the location of all North American Breeding Bird Survey (BBS) routes that have been surveyed one or more times (n = 3924; approximately 200 of these routes are no longer suitable for use, see text). Thick black lines show state and provincial boundaries and thin green lines show boundaries between the 99 physiographic strata (Bystrak 1981) that are often used in analyses of BBS data. Figure 1.2. Contour map showing the number of annual surveys on BBS routes between 1966 and 1995. On average, routes in central and eastern states have been surveyed more times than those in western and northern states and provinces. Contours were generated using inverse distancing (Environmental Systems Research Institute Inc. 1996). Figure 1.3. Number of BBS surveys performed each year between 1966 and 1995. -10- Figure 1.1 -11 - Number of surveys per route between 1966 and 1995 0-10 10-20 20-30 - 12- Figure 1.3 Chapter 2 The Importance of Analysis Method for Population Trend Estimates Introduction There is widespread concern that many North American landbird populations are undergoing long-term declines (Terborgh 1989, 1992; Askins et al. 1990; Askins 1993; Finch and Stangel 1993). For these species, the primary source of population information at the continental scale is the North American Breeding Bird Survey (BBS). Trend determinations from the BBS are widely used by government agencies, conservation organizations and independent researchers to set priorities for conservation research and management (e.g., Hunter et al. 1993; U.S. Fish and Wildlife Service 1995) and to test possible hypotheses of causation (e.g., Sauer et al. 1996; Heckert 1997). Unfortunately, estimating trends from the BBS is a not straightforward; a number of statistical methods have been proposed and applied but there is currently no consensus as to their relative merits (Chapter 1). Here, I focus on the practical consequences of the application of different analysis methods for those that use BBS trend information. I compare trend estimates produced by different analysis methods applied to the same data, and discuss the implications of any observed differences in estimates for our ability to detect population declines and to prioritize species of conservation concern. I use three analysis methods: USNBS route regression (Geissler and Sauer 1990), CWS route regression (Erskine et al. 1992) and rank-trends analysis (Titus 1990) on a common set of BBS data from British Columbia. The first two methods were used by the two agencies in recent published reports (Erskine et al. 1992; Peterjohn et al. 1995) and are based on a similar approach (route regression), but differ in the way that this approach is implemented. The third method is a nonparametric technique that estimates the direction of trend (i.e., whether the population is increasing or decreasing), but not its magnitude. - 14- Methods Breeding Bird Survey The Breeding Bird Survey comprises approximately 3700 39.2km long survey routes selected at random along secondary roads within all states/provinces in North America. Each year, during late May or June, a volunteer observer drives the survey route, stopping every 0.8 km and recording all the birds seen or heard within a 400 meter radius in a three minute period (Robbins et al. 1986, Appendix B). Approximately 70% of routes are surveyed each year. Here, I used all routes in British Columbia that were surveyed between 1968 and 1992, excluding those from physiographic strata 29, 30 and 68 (Closed Boreal Forest, Aspen Parklands and Northern Rocky Mountains respectively, Butcher 1990), since they are rarely surveyed, and are not normally included in quantitative analyses of BBS data (J. R. Sauer, USNBS Patuxent Environmental Science Center, pers. comm.). The data set thus included 82 routes (803 individual surveys) within physiographic strata 63, 64 and 94 (Fraser Plateau, Central Rocky Mountains and Northern Pacific Rainforests). Data selection Population trend analysis from BBS data involves two stages: data selection and trend estimation. Data selection criteria are designed to screen out surveys performed under questionable conditions, and to ensure that there are sufficient data to meet the requirements of the trend estimation procedure. Criteria depend largely upon the judgment of the analyst and vary widely among published analyses. The selection procedures used by the USNBS and CWS are similar, but differ in the criteria for the exclusion of individual surveys and the minimum data requirements of the route regression procedures. Standard data selection protocols for ranktrends analysis have not been developed. Because I was chiefly interested in differences due to trend estimation method, I based my comparisons of the three methods on the common subset of data that passed both current USNBS and CWS pers. comm.) and B. T. Collins (CWS criteria, as provided by J. R. Sauer (USNBS, National Wildlife Research Centre, pers. comm.). Nevertheless, in order to determine the effect of data selection, I also compared the results from -15 - the two route regression methods produced using their respective agency's data selection protocols. USNBS criteria: (1) Only surveys performed within a regionally specified range of dates, times and weather conditions, and by qualified observers were included in the analysis (i.e., surveys with USNBS survey codes one, two or three). (2) Routes with too few surveys remaining to enable estimation of trend and trend variance were then excluded (i.e., where number of surveys minus number of observers was less than two). (3) Routes on which no birds of the species under analysis had been seen were also discarded. (4) If (after data selection) the species had fewer than 10 routes remaining or had a mean count of less than 1.0 birds per survey, trend estimates were not calculated for that species. CWS criteria: (1) Surveys that did not meet CWS survey criteria were discarded. These criteria are similar (but not identical) to those used in USNBS criterion 1. (2) Routes with too few surveys remaining to enable estimation of trend (but not trend variance) were then excluded (i.e., where the number of surveys minus number of subroutes was less than one. A subroute is a subset of the surveys done on a single route such that all surveys in a subroute were performed by the same observer within a span of 19 or fewer calendar days across years and under similar weather conditions. Subroute designations were provided by B.T. Collins). Criteria 3 and 4 were identical to those used by the USNBS. Trend estimation The three trend estimation methods are summarized in Table 2.1 and are described in detail below. USNBS route regression: In route regression, the overall trend for the region and species of interest is calculated using the weighted average of the trend on each route (see below). A route trend is the slope of a log-linear regression of birds counted (dependent variable) against time (independent variable). I performed the analysis described here using programs written in C and C++ (Appendix 1), following the methods of Geissler (1984), Geissler and Sauer (1990) and Peterjohn et al. (1995), with additional details provided by J. R. Sauer (pers. comms.). The programs included some numerical routines written by Press et al. (1992). - 16- For each route, I computed the natural logarithm of the number of birds counted (plus 0.5 to accommodate zero counts; Component 1, Table 2.1). Route trends were calculated on the logarithmic scale as the slope of the regression of log count on year of count (Component 2). To account for differences among observers I treated them as covariables in the regression (Component 3). I then back-transformed the route trends using the method developed by Bradu and Mundlac (1970), which is approximately equal to taking the exponent of the trend estimate minus one-half of its variance (Component 4). The overall species trend for British Columbia was calculated using the weighted mean of the back-transformed route trends (Component 5), in two stages as follows. I first calculated the weighted mean slope for each physiographic area, weighting back-transformed route trends by the marginal mean count on the route (Searle et al. 1980; Component 6), and a measure of the reliability of the route trend estimate, calculated as the variance of the slope divided by the mean square error for the route-specific regression (Component 7). The marginal mean is an estimate of the count that would have been recorded in the mid-year of the analysis period (i.e., 1980) by an average observer. Geometric mean count was used in place of marginal mean on routes where estimating the marginal mean required extrapolating beyond the years when counts took place. In this context geometric mean was defined as the arithmetic mean of the log-transformed counts, backtransformed onto the multiplicative scale. I then calculated the overall trend as the mean of the trends for each physiographic stratum, weighted by the stratum area (Component 8). I estimated the variance of the overall trend estimate from 400 bootstrapped subsamples of the routes within physiographic strata (Component 10). Bootstrapping is a technique that allows the estimation of the sampling distribution of a non-standard statistic (in this case the overall trend estimate). Many repeat subsamples are taken with replacement from the original data and the distribution of the statistic in the subsamples is taken as the best indication of the true distribution of the statistic (Efron 1982). I also used the mean of the bootstrapped trend estimates as the final species trend since it is in theory less biased than the original estimate. Statistical significance of the trend was assessed using a z-test (Component 11). This tested the -17- alternate hypothesis that routes show consistent log-linear trends with a slope different from zero across physiographic strata within the study area. CWS route regression: The CWS method also uses route regression to calculate regional trends from a weighted average of the trend on each route. However it differs from the USNBS method in the way the route trends are calculated, the method of averaging route trends and the weightings used to produce the regional trend. I performed this analysis using a F O R T R A N program supplied by B. T. Collins (Appendix 1), which implements the methodology of Erskine et al. (1992; see also Collins and Wendt 1989). As with the USNBS route regression, I calculated route trends as the slope of the regression of the natural logarithm of counts on year of count (Components 1, 2 and 3). In this case 0.23 (rather than 0.5) was added to each count before log transforming to accommodate zero counts, and subroutes were used as covariables in the regression (rather than observers; see CWS data selection criteria for a definition of subroutes). I calculated the species trend for British Columbia directly from the weighted mean of the route trends, still on the logarithmic scale (Component 5). The weightings used were the marginal mean count on the route (Searle et al. 1980; Component 6), a precision estimate given by the squared deviations of the counts from the subroute mean (Component 7), and an area weighting (Component 8). I constrained the marginal mean so that it could not exceed the maximum number of birds seen on all routes, or 200, whichever was the least. The area weightings were provided by B.T. Collins and these consider the number of routes and the proportion of land area in each degree-block. I then back-transformed the overall trend estimate to the arithmetic scale by simple exponentiation (Component 9), and calculated the variance of the estimate by jackknifing (Component 10). Jackknifing is similar to bootstrapping, except that in this case the population distribution of the statistic (overall trend estimate) is approximated using the set of subsamples generated when one data point (route) is removed from the sample in turn (Effon 1982). Unlike the USNBS method, no bias adjustment was performed. The statistical significance of the alternate hypothesis, that there are consistent non-zero log-linear trends across routes in British - 18- Columbia, was calculated using a /-test, with the degrees of freedom equal to the number of routes minus one (Component 11). Nonparametric rank trend analysis: With this method, trend statistics are calculated separately for each route using a nonparametric ranking procedure, and are summed to produce regional trend statistics. This enables calculation of the direction (increasing or decreasing) and statistical significance of the trend, but not its magnitude (the size of increase or decrease). The procedure described in Titus et al. (1990) was implemented fully with a C program (Appendix 1), and will be briefly summarized here. On each route, I arranged counts in ascending order and assigned ranks, giving ties the mean of the tied ranks. I calculated a trend statistic, D, for each route as the sum of (/?,• - i) , 2 where i?, is the rank of the fth yearly count (i.e. / = 1 in 1968 to 25 in 1992; Component 2). If counts tend to increase over time on the route then D will be small, and if counts tend to decrease over time then D will be large. Under the null hypothesis of random change in count over time, this statistic is normally distributed with an expected mean and variance that can be easily calculated (Lehman 1975). To test the alternate hypothesis, that counts on routes tend to increase or decrease over time throughout British Columbia, I summed the values of D (Component 5), the expected mean of D and the expected variance of D (Component 10) over routes and performed a z-test (Component 11). Comparison of trend estimates For each pairwise combination of methods, I quantified differences in trend estimates by calculating the proportion of species with trends in different directions (i.e., increasing in one method but decreasing in the other). I also calculated the proportion of differences in statistical significance (i.e., not significant (p<=0.05) in one method but significant (p>0.05) in the other). In order to determine whether the differences were related to sample size (number of routes) or species abundance (mean count per survey), I performed logistic regressions on the proportion of differences (dependent variable) against sample size and log species abundance (independent variables). -19- For the two route regression methods, I further examined the pattern of differences by plotting estimates of trend against one another and determining the slope of the principal axis (Model II regression, Sokal and Rohlf 1981). If the two methods produce similar trends, on average, then the principal axis will have a slope of 1.0 and run through the origin. I also calculated the median absolute difference between trends, and determined whether this was related to sample size or log species abundance using linear regression. Results Of the 265 species recorded at least once, 119 passed both USNBS and CWS data selection standards. For this common subset of the data, four species had estimated annual rates of population change of greater than 15% per year using CWS route regression (Tennessee warbler Vermivoraperegrina: +18.8%; Black tern Chlidonias niger: +22.5%; Magnolia warbler Dendroica magnolia: -33.6%; Black swift Cypseloides niger. +41.8%; see Discussion). These extreme trends are biologically implausible, so I excluded them from further analysis. For the remaining 115 species, the median sample size was 39 routes (range 10-61) and the median abundance was 4.76 birds per survey (range 1.09 - 62.86). Using these data, about the same number of species were estimated to be increasing or decreasing in all three methods (Fig. 2.1). However, there was considerable variation among methods in the species making up these totals: the proportion of species assigned different trends ranged between 0.21 and 0.33, depending upon the pair of methods being compared (Table 2.2). The proportion of differences did not appear to be strongly related to sample size or species abundance (Table 2.3). The estimates had quite large confidence intervals (Table 2.3), so I cannot rule out the possibility that there were too few samples (N = 115 species) to detect such an effect had one occurred. However, in many cases, the estimates were in the opposite direction to that predicted (i.e., proportion of differences increased with increasing sample size and abundance). There were large differences among methods in the overall number of species estimated to be undergoing statistically significant population changes, with USNBS route regression -20- assigning almost twice as many species significant trends as CWS route regression, and nonparametric rank trends assigning twice as many again as the USNBS method (Fig. 2.1). The proportion of species assigned trends of different significance levels was between 0.18 and 0.43 (Table 2.2), and again did not appear to be strongly related either to species abundance or sample size (Table 2.3). The comments about sample size are relevant here also. In all comparisons (i.e., direction and significance), the two route regression methods were the most similar, while the greatest differences were between CWS route regression and rank-trends analysis (Table 2.2). For the two route regression methods, there was little overall difference in trend estimates (median difference USNBS trend - CWS trend = +0.04% per year; 95% confidence interval (CI) estimated from 10 000 bootstrap subsamples of the 115 species = -0.60 to +0.44). However, the absolute magnitudes of trends were greater in CWS route regression trends than in the USNBS method: negative trends tended to be more negative and positive trends tended to be more positive (Fig. 2.2., slope of principal axis = 0.58, 95% parametric CI = 0.45 to 0.72). There was considerable scatter about this overall pattern (Fig. 2.2), with the median absolute difference in magnitude of trends being 1.19% per year, (range 0.00 - 13.06). Contrary to the results for the proportion of differences in the direction and significance of trends, the mean absolute difference in the magnitude of trends decreased with increasing sample size (slope of linear regression: -0.46; 95% parametric CI = -0.87 to -0.05) and with increasing log species abundance (slope: -0.044; 95% parametric CI = -0.067 to -0.020), although in both cases the proportion of variance explained by the relationship was small (r = 2 0.04 and 0.11 respectively). I speculated that this may be due to decreasing variance of the estimates, since variance decreases with increasing sample size, and sample size and log abundance were correlated (Pearson r = 0.44). In order to test this, I corrected for decreasing variance by dividing the difference in trend estimates by the square root of the sum of their variances, to yield z-statistics. The absolute magnitude of these z-statistics did not appear to vary with sample size or with log abundance in linear regressions (sample size: slope = -0.002; 95% parametric CI = -0.008 to +0.005; r = 0.00; log abundance: slope = 0.05; 95% parametric 2 CI = -0.06 to +0.16; r = 0.02), indicating that the increasing similarity between trend estimates 2 -21 - with sample size and abundance was indeed explained by decreasing variance. I conclude that trend estimates that have low variance when derived using route regression are less sensitive to the particular variant of route regression used. Finally, I evaluated the effect of data selection criteria. When only USNBS data criteria were applied to the data, the average species dataset contained 3.5 more routes (median sample size = 42.5 routes, range 11 - 64), and the median trend estimated using USNBS route regression was 0.45% per year higher than when both data selection criteria were applied (range 1.20 lower to 6.02 higher). For CWS data selection criteria alone, the average increase in sample size was 2 (median sample size = 41 routes, range 10 - 62), and the median trend was 0.07% per year lower (range 4.36 lower to 2.90 higher). Comparing USNBS route regression and CWS route regression using trend estimates derived from the datasets appropriate to each, these changes caused the USNBS method to produce on average more positive trends than the CWS method (median difference USNBS trend - CWS trend = +0.73% per year; 95% CI estimated from 10 000 bootstrap subsamples = +0.36 to +1.01). This is dissimilar to the result using the uniform data set, where both methods produced similar trends on average. I therefore conclude that data selection criteria can have an important effect on the results obtained. Discussion The Breeding Bird Survey is the most comprehensive source of data available for monitoring trends in North American landbird populations, and for prioritizing species of conservation concern. I have shown that the method of trend estimation can affect the magnitude, direction and statistical significance of population trends assigned to species. In addition, small changes in the way that data are selected for analysis can affect the overall magnitude of trends. In this section I consider possible causes of the differences, their implications for users of BBS trend information, and suggest some priorities for further work. C a u s e s of the differences This study was not designed to explicitly test which of the many differences between the methods caused the observed differences in trend estimates. Nevertheless, based on the -22- predicted effect of differences between the methods (Table 2.1) and some preliminary analyses, I can suggest which components of the methods could account for the results. In summary, the results were: (1) large differences in the number of statistically significant trends among the three methods; (2) greater absolute magnitude of trends in CWS route regression compared with the USNBS method; (3) slightly more positive trends in the USNBS method compared with CWS route regression when the respective agencies' data selection criteria were applied separately; and (4) large variation about these patterns in the direction, magnitude and significance of trends. I deal with each of these in turn. Differences in the number of statistically significant trends were consistent with differences in the way that regional trends and trend variances were calculated by the different methods, although many other factors may have been involved. Nonparametric rank trends analysis produced the greatest proportion of statistically significant results, probably because it did not consider between-route variance when calculating the expected variance of the regional trend (Component 10, Table 2.1). USNBS route regression produced more statistically significant results than CWS route regression. The USNBS method tested for consistency of trends nested within physiographic areas, while CWS route regression tested for consistency of trends over the whole region (Component 5). Hence in the USNBS method, the variance of trends among physiographic areas was removed from the statistical test (analogous to a nested analysis of variance), resulting in an increased chance that the trend was significant. The greater absolute magnitude of trends in CWS route regression relative to the USNBS method may have been due to differences in the constant added before log-transformation (Component 1). Addition of a constant biases the trend towards zero (no trend), and this effect is stronger for larger values of the constant (Geissler and Link 1988; Collins 1990). The constant was larger in the USNBS method than in the CWS method (0.5 vs. 0.23), so the potential bias in USNBS route regression was greater. The bias is also stronger for low abundance species (Geissler and Link 1988; Collins 1990). I thus predicted that i f differences in the constant were causing the observed differences in absolute magnitude of trend, the effect would be greater for low abundance species than those with high abundance. I tested this -23 - possibility by dividing the 115 species into two equal groups by abundance and calculating the slope of the principal axis separately for each group (one species was excluded from the test to make both groups of equal size). The median abundance in the low abundance group was 2.08 birds per survey (range 1.09 - 4.59), and the slope of the principal axis was 0.52 (95% parametric CI = 0.32 to 0.76). For the high abundance group, the median abundance was 8.29 birds per survey (range 4.96 - 62.86) and the slope of the principal axis was 0.63 (95% parametric CI = 0.49 to 0.79). The slopes vary in the anticipated direction (low abundance group smaller) but are clearly not significantly different. Simulations performed by Collins (1990) showed that the bias due to addition of a constant decreases rapidly with increasing abundance. I thus suspect that the minimum data requirement for a mean of at least 1.0 birds per survey was effective in controlling the difference in bias to below the level that I could detect statistically. Other explanations must account for most of the observed difference in the absolute magnitude of trends between route regression methods. The systematic change in USNBS trend estimates resulting from changes in data selection criteria currently defies explanation. I cannot imagine any scenario in which the small differences in data selected would cause systematic differences in trend. These results highlight the need for standardization of data selection, since many published analyses use data selection criteria even more divergent than the two I compared here (e.g., Bohning-Gaese et al. 1993; James etal. 1996). Differences in trend estimates among methods that are not explained by the above patterns are likely due to the additive effect of the many differences among the methods, each acting independently. Of these, I predict that the most important are the weightings used when combining route trends to produce trends for the region (Components 6, 7 and 8). These weightings are designed to correct for geographic variation in the number of routes and the abundance of birds, and to increase the precision of the regional trend estimate. The method of calculating the weightings differs between route regression methods, which can result in differences of up to three orders of magnitude in the weightings assigned individual routes. No system of weightings has been developed for use with nonparametric rank trends analysis. The -24- importance of weightings is demonstrated by the four species assigned improbably large trend estimates in CWS route regression and excluded from the comparisons. In all four of these species, one route with a large apparent change in the density of birds but few surveys was assigned a very large abundance weighting (due to the extrapolation of the route trend to calculate the least squares mean) and thus dominated the regional trend estimate. The abundance weighting was calculated differently in USNBS route regression (geometric mean was used where extrapolation was necessary), and these routes were assigned small weightings, leading to more credible regional trend estimates for the four species. Implications of the R e s u l t s To be effective, conservation programs must be based upon accurate information. Partners in Flight, the multi-agency coalition responsible for initiating the Neotropical Migratory Bird Conservation Program (Hagan 1992; Finch and Stangel 1993), has defined an effective monitoring scheme as one that has a 90% chance of detecting a 50% decline in species abundance over 25 years (Butcher et al. 1993). I have shown here that for British Columbia B B S data the direction, magnitude and significance of trends attributed to species shows considerable variation due to analysis method alone. For the two route regression methods, which used a similar approach, 15 of 115 species (17%) showed a difference in the magnitude of the trend due to analysis method of greater than 50% over the 25 years analyzed. For these species the uncertainty due to analysis method alone was larger than the absolute change in population size that Partners in Flight wishes to detect. Differences in the magnitude of route regression estimates were related to sample size (number of routes in the analysis); hence, I expect these differences to be smaller for estimates based on larger geographic areas or states/provinces with a greater density of routes (e.g., those in eastern North America). I do not know whether trend estimates from methods based on very different assumptions, such as nonlinear nonparametric route regression (James et al. 1992) or Mountford's method (Mountford 1982), will show similar convergence with increasing sample size. A number of species prioritization schemes have used the statistical significance of the trend estimate as one criterion for setting the conservation priority of a species (e.g., Hunter et al. -25- 1993). The results show large differences among methods in the statistical significance assigned to species (Table 2.2), with no consistent increase in similarity with increasing sample size or abundance (Table 2.3). Of the 31 species estimated to be showing significant declines (i.e., p<0.05) by one or more of the methods, 15 (47%) were estimated to be declining significantly by USNBS route regression, 8 (25%) by CWS route regression and 29 (91%) by rank-trends analysis. Hence species prioritization schemes based upon the statistical significance of trends are likely to highlight a very different set of species depending upon which analysis method is used. The differences reported here indicate that analysis method will also affect the results of geographic analyses of the significance of population trends (e.g., Sauer and Droege 1990a, d; Peterjohn et al. 1995). I do not know whether analysis method will influence the outcome of multispecies comparisons (e.g., Bohning-Gaese et al. 1993; Peterjohn and Sauer 1993; James et al. 1996), since differences in individual species trends may "average out" when species are grouped (Fig 2.1). I thus tentatively support Droege (1990:3) in his statement: "Statistical analyses of these data and their subsequent interpretation should dwell on the patterns of population change rather than on the magnitudes of calculated trends and variances." Priorities for Further Work In this chapter, I have demonstrated that the manner in which BBS data are selected and analyzed can affect the result obtained and inferences drawn. These results raise two important questions: 1) which components of the methods have the largest effect on the resulting estimates, and 2) which method is the most accurate? The first question requires a detailed comparison of many variants of each method applied to the same dataset. This type of study would act as a sensitivity analysis, enabling detailed methodological research (e.g., Sauer et al. 1994b) to focus on the aspects of the methods that are the most sensitive to the approach used. The second question is much harder to resolve. Comparisons of results among methods are useful for highlighting problems with the methods that lead to unfeasible estimates (e.g., the abundance weighting problem in CWS route regression). However, these comparisons cannot definitively tell which method is the most accurate because as there is no benchmark of known trends with -26- which to compare the estimates. True population trends in British Columbia, or for any other large geographic part of the continent, are unknown for any species monitored. One possibility is to compare results to those from other surveillance programs, such as migration monitoring from bird observatories or checklist programs (e.g., Hussell et al. 1992). However, this approach is flawed: because all programs measure trends with error, neither concordance nor lack of it can be used to infer which method is biased. Monte Carlo studies offer a possible solution, in which "known trends" are generated a priori using stochastic simulations, and are used to compare analysis methods. Trends could be simulated under a variety of realistic scenarios, such as those suggested by Wilcove and Terborgh (1984). In general, the most appropriate method will be the one that gives the most accurate results over the widest range of conditions. A n example of the approach is given in Chapter 3. The BBS has the potential to provide population trends for the majority of breeding birds in North America over a large part of their ranges. Although trend estimation from B B S data is a formidable statistical challenge, many of the causes of bias are predictable and can be minimized by the choice of a method appropriate to the hypothesis under examination. Given the current high levels of concern about population trends of North American landbirds it is clearly imperative that we determine how best to draw inference from the wealth of information that is collected by the BBS each year. Summary The purpose of this chapter was to determine whether analysis method has an important influence on the inferences drawn about population trends. I compared trend estimates produced by three different analysis methods applied to the same set of BBS data for 115 species in British Columbia. The methods were U.S. National Biological Service route regression, Canadian Wildlife Service route regression, and nonparametric rank trends analysis. Overall, the number of species estimated to be declining was similar among the three methods, but the number of statistically significant declines was not similar (15, 8 and 29 respectively). In addition, there -27- were many differences among methods in the trend estimates assigned to individual species. Comparing the two route regression methods, Canadian Wildlife Service estimates had a greater absolute magnitude on average than those of the U.S. National Biological Service method. U.S. National Biological Service estimates were on average more positive than Canadian Wildlife Service estimates when the respective agencies data selection criteria were applied separately. These results imply that our ability to detect population declines, and to prioritize species of conservation concern are strongly dependent upon the analysis method used. I discussed the possible causes of the observed differences, and suggested the use of simulation studies to evaluate the performance of the methods. -28- CO u « ^ ^a C/3 N 11 a~ GO CO <u o c S-i g S g cd ^ <+-i a a« a lo o o o 00 GO bo -2 ^ « *a !+-< m .5; .2 * CO <^ 2 ^ a CO <L> bo I-I p aa o t-l CD ES ° CO CO O O co <U w O > 8 &° .*-i CS 1 ^ PQ C5X) J +3 a C/D c0 o •2 U > T3 a ii a o CtH CO s-i 'to + <D r-i VH a bo o u CO J—i lop O (D o ~ I CO PQ g o + o o bo o o — bp u t) o N m co (-< o 8 co GO O .2 l-l x i co ^ St CO <4H o CO co I<L> O « 9 GO <L) I-i 3 O O O 1 u O ci a a .—i CO CO a Z3 T3 CO 1) o 3 CO -. P< o a O <u b0 <+< o 0 03 (D a^ 3 d> - 2 CO 3 ^—* no. u 1) CO o o egness co O ES o o H ubr co O bO CN o a o •s O t-i > O O CO -29- eg CO 8 C (U is o '5b 3 13 c o 'S CD N CD bO oo GO CD TO .2 o CD O Li ^ GO o g CD (L) a o CO CD a o c T3 o TO o o GO U a o 3 CO 8 1 * 3 fe O rv O S .° o Jf TO o '5b u M t-i > u CD cu g a> o ON ON II " PL, <D CD £ A o TO CD 43 r j 1/3 CD > <M TO T—I -Pi • • • M CO pi TO TOPQ _ <D 2 — I O GO *t3 O CO TO 4 co 2 ON OO ON • I—I PQ H ro CD g O CD +-» O O i ro oo t3 <z> CD ON ON CD 8 43 u a o S o U > © -30- GO o H-> GO PQ DO CD GO Si GO 33 CD CO CO co '5 o .5 <o U H CD CD <D CD CD CD GO GO GO 00 00 .2 £ 00 oo K> ra S W l-i TO S I-I bO »-• o • O^ o ©' o CN »/-> O O -t-» vo CN cd 00 o a oo > 00 £ ^ G3 S <u _o is 'oo co . , <D ^ *-< C bfl cd 1) I-i CN CO © O -*-> 00 o §1 I-I g PQ & o © o CN a . u• 8 £ 'S 00 CN o ^ ©' 2? a 12 Cd <D oo > .2 'oo 00 Sb C3 O oo 00 oo oo CN M a 2 "3 3 P 2 oo co PQ ^ £ u © o CN © in CN © o ©, oo © £ g 8 o •g ° C i 1 N? o> ?? '5b ' £ o _ 00 PQ "c3 c o .2 oo o o T3 a o o •a OT PQ of S OO cd bp bo -31 - CO CN cu J3 cd CN in o o oo cn o o o H-> o CN CN o O o c-<N v ' VO o o o i—i o H-> « 8 S3 VO m Js o o • TO 00 3 o 42 12 P « M« is" .2 CO P O C CD '55 +3 CO w cu 4^ o CU bo ' £3 CU o cu »-H CN 0 ^ O 0 o H—» cn cn 4-* O -»-» O -»-» OO 00 vo 0 i 0 ^1^ ON O 1 o a O T7i ^ . 1 vo 1—( r O 0 O 0 o CN r-l 'co CO CO cu CU cu £ o CO s 0 CN m >— in in «n 00 O O 0 1 O w O ,—1 <— 1 O 0 -*-» ,—I m O H-» CN CN O O ^Ls in VO CN O O O O O >—1 CO c^ CU 4 i cu cucPuH 7? c d cu .S cou -a CU o .a 60 O a 5 T5 00 o cou S3 cu cu p 3 .2 o cl-Hu T3 P 1) o Q > P 5 o o CU 4 3 -J-i 8 Q. CO T3 H-» CU 3 o CO T3 a CU 43 <+-< o u 4 ? o O PH < u o •g CL Pi O GO <U ' O S cu H CU PI Pj CU o •5 43 T3 CU CC! cu CU O cu CL £ • Pi 60 O T3 3 PH CO .a CD CO 3 CU a 8 ^5 c t3 J*-, CO 4c3u CCi C H co pa cu o o CO ^ 3 II S3 cu co o 4 3 4 3 > ' £ cu o T3 P J 43 S cu H3 "5 H—» 0 O u 0 O 0 43 ^ Pj o £ S co 3 ^3 o S CcO CU CU o .§ '3 c C+H CO cu ^ « n CO > M 43 ccj cou &0 co 3 CO O O 1-1 _ O r-l a cn 1—I 0 p & CN VO (N > w P cd - cPo ra o o o VO & in <N 43 oo .a a cu H ^ S^ u o o > cd ^ ts 13 3 c o c j x a CU io 4 = « o cu O N cd cu ' § ell cu t3 cu CO 3 > « -tl PJ a cu +^e cu PH 2 bO > _!3 CO pi CO c u S3 CU O PQ 43 cu CO ^ 2o Z 43 C O ro cu CQ P >49 -32- Figure Legends Figure 2.1. Number of species with increasing and decreasing trends from three analysis methods using the same Breeding Bird Survey data for 115 species from British Columbia. USNBS is U.S. National Biological Service and CWS is Canadian Wildlife Service. Figure 2.2. Plot of the magnitude of trends estimated from U.S. National Biological Service (USNBS) route regression against the magnitude of trends from Canadian Wildlife Service (CWS) route regression for 115 species. If the methods produced identical trends then all the species would be on the dashed line. The solid line is the principal axis. -33 - Figure 2.1 60 r CD w CD CO CD CD »- 8 = 40 20 1- Q. CO 0 CD -Q E 3 C CO CO CD i_ 20 h p< O p > 0 . 0 5 40 o CD Q 60 80 =0.05 _L USNBS route regression -34- CWS route regression Nonparametric rank trends Figure 2.2 1.10-j 1.08H 0.940.92- 0.90 0.92 0.94 0.96 0.98 1.00 1.02 1.04 1.06 CWS route regression Population change (% per year) -35- 1.08 1.10 1.12 Chapter 3 Issues in the Estimation of Long-term Population Change Introduction A number of issues make the analysis of long-term multi-site population monitoring data challenging and controversial. In this chapter I explore these issues, and suggest one possible method for resolving them. In the first part of the chapter I consider how the data could be analyzed in an idealized survey, and then discuss the effect of various practical constraints on the analysis. As a concise reference, I summarize the major features of analysis methods that have been used to date in Tables 3.1 and 3.2. In the second part I suggest the use of Monte Carlo simulations to empirically test the performance of these different methods under plausible scenarios of population change, and provide an example of the approach. Much of the emphasis is on the estimation of long-term population change in North American landbirds from the Breeding Bird Survey (BBS); however these methods and the discussion are applicable to any population where repeated counts are performed at a large number of sites. Analysis of survey data under ideal conditions Since it is not feasible to census the entire population of most bird species, population monitoring is almost always based upon surveys of a sample of the population. In an ideal sample survey, the population of interest is divided into non-overlapping units, a sample of these units is selected according to a pre-defined procedure and a complete census is taken of each unit. In bird surveys, the sample unit is usually all the birds within a specified survey site ("route" in BBS terminology). Because population change (i.e., change in abundance) is the main parameter of interest, the most efficient procedure involves repeated annual censuses of the same sample of sites. Two issues arise in the estimation of population change from the resulting data: (1) selection of the unit of analysis and (2) model of trend. -36- Selection of the unit of analysis Population change is usually portrayed using annual indices of abundance or as population trend (prevailing tendency underlying annual indices). In the ideal situation, annual indices are given by the average of the counts in each year, and population trend is the average of the trends at each sample site. If the sample design involves stratification then the average within each stratum should be weighted by the proportion of units in the stratum, usually calculated as the proportion of the total area covered that is in the stratum (area weighting). In addition, trends should be weighted by the abundance of individuals at the sample site (abundance weighting), since declines at a site containing many individuals represent a greater proportion of the population than declines at a site containing few individuals (see James et al. 1990 for an overview of weightings). Alternatively, population trends can be calculated from the estimates of annual indices, or annual indices could be calculated as the residuals from a trend analysis. In the ideal situation, these indirect estimators produce the same results as using counts as the unit of analysis. However, because the analysis methods used in practice are modified to deal with departures from the ideal, the estimates may diverge. Hence, analyses that present both annual indices and population trends (e.g., Greenwood et al. 1994, Peterjohn et al. 1995) tend to use estimates of one parameter to produce the other. The unit of analysis for each of the methods is given in Tables 3.1 and 3.2. Model of trend Repeated counts at a survey site form an ecological time series, and so may be thought of as being composed of four sources of pattern: trend (prevailing tendency), interventions (irregular perturbations such as environmental effects), autocorrelation (due to population processes such as density dependence) and sampling error (Barker and Sauer 1992). Estimating population trend thus involves separating the prevailing tendency from the other components. Because the relative contribution of these components is not known a priori, the model used depends upon the judgment of the analyst (Jassby and Powell 1990). Many of the major differences between trend-estimation methods can be attributed to differences in the way that -37- trend is modeled. Although there are a number of methods available for deriving trends from time series (Jassby and Powell 1990), analyses of bird counts have focused on regression-based approaches. Four regression models have been applied to count data (Table 3.2), each of which can produce a different picture of the prevailing tendency (Fig. 3.1), and each of which may be advantageous under certain conditions. (1) In linear-multiplicative models (e.g. linear route regression), population change is modeled as a multiplicative process in which the expected count in year.y+1 is a multiple of the expected count in yeary, leading to E(C } = E[C }B y y 0 where E((J ) is the expected count y in year y and p is the trend. This model is conceptually simple, involves estimating the minimum number of parameters and is unbiased in the face of random interventions and autocorrelation. However, it will be inaccurate if the prevailing tendency changes over time (i.e., non-linear trend). (2) Polynomial-multiplicative models (e.g. loglinear modeling) can allow for changes in trend over time, using polynomial parameters to estimate higher order effects, i.e., E{C ^ = E(C )p y y 0 p{p{.... The order of the model is usually determined a priori, but could be derived from the data using model selection techniques, such as likelihood ratio tests (Burnham and Anderson 1992). Although more general, it is also more susceptible to incorporating inappropriate patterns such as population cycles. (3) In additive models (e.g. smoothing route regression), no apiori model is postulated, and non-trend sources of variation in counts are removed through filters such as LOESS. The advantage of this model is that no explicit assumptions are made about the shape of the trend. One disadvantage is the requirement for specification of a smoothing parameter, which determines how much variation is filtered from the data. This must be assigned subjectively. A common criticism of both polynomial-multiplicative and additive models is that it is hard to draw inferences about overall patterns of population change since they estimate many parameters. Both can, however, provide a single parameter estimate i f -38- trend is defined as the difference in expected count between the beginning and the end of the time period of interest (e.g., James et al. 1996). (4) In rank models (e.g. rank-trends), trend is modeled as the tendency for counts to increase or decrease over time. Counts are ranked, and the distribution of ranks over time is used to estimate whether the population has shown a tendency to increase or decline (Fig. 3.1). Because ranks are used, only the direction of trend is estimated, not its magnitude. This model is only appropriate if the prevailing tendency is uni-directional, but does not require that the trend be linear. Ultimately, the most appropriate model of trend will depend upon the pattern of population change in the population being modeled. While the true pattern is not known for any species monitored, the robustness of alternative models given expected population changes over space and time can be tested using population simulations (e.g., Anganuzzi 1993). Note also that the accuracy of all models is conditional on the length of the time series and the validity of the starting point (Jassby and Powell 1990, Barker and Sauer 1992). For example, i f a species shows population cycles but the time series covers only part of one cycle then the observed changes in abundance will be interpreted as trends. A second example is the increasing trends recorded for some species in BBS analyses which may be due to recoveries from pesticide use or severe winters that occurred near the start of the survey (references in Barker and Sauer 1992). Additional factors that complicate the analysis of survey data Extensive bird surveys often depart from the ideal situation in three ways: (1) the population from which the sample is drawn is a subset of the population about which information is required (incomplete sampling frame), (2) the number of individuals in each sample site is not counted accurately (measurement error), and (3) most sites in the selected sample are not censused every year (missing data). Incomplete sampling frame If the sample frame is incomplete then the degree to which results obtained from the sample can be extrapolated to the population as a whole is limited. For example, in the BBS, -39- survey sites are allocated along secondary roads; hence limited inferences can be drawn about species in which a large proportion of individuals breed in habitats such as floodplains, mountains or forest interior that are not in close proximity to roads, or in the roadless northern part of Canada. Although this can be an important limitation to the utility of the survey program (Droege 1990; but see Bart 1995), it does not affect the way that estimates are derived from the population sampled, and so is not considered further here. Measurement error Measurement error is a feature of almost all bird counts. The probability that a bird is counted given that it is present at the survey site that year (its detectability) is a function of a large number of factors, such as survey method, species, observer, habitat, time of year, weather, bird density, and breeding phenology (references in Droege 1990). Hence, estimates of population size or absolute abundance are not possible from the count data, and indices of abundance are seldom comparable between species or survey programs. In addition, sources of error that change systematically over time, such as habitat, traffic noise, inter- and intra-observer ability can generate spurious estimates if they are not modeled in the analysis method (Barker and Sauer 1992). In particular, differences in detection rates between observers are often cited as a major potential source of bias (Barker and Sauer 1992, Sauer et al. 1994b, Peterjohn et al. 1995). Sauer et al. (1994b) used several approaches to show that the number of birds counted by observers has increased, independently of trend, for the majority of species in the BBS. This result was corroborated by James et al. (1996), who found that trends were usually more positive if observer differences were not modeled in their nonlinear trend-estimation technique. However, incorporating observer differences (or other systematic sources of error) into models may result in a drastic reduction in efficiency due to overfitting, and the resulting accuracy of the estimate may be decreased. Currently, some analyses model observer differences, while others do not (Tables 3.1 and 3.2). The problem of determining whether to include additional explanatory variables in a model has been encountered in other areas of ecology (e.g., mark-recapture studies: Lebreton et al. 1992), and the use of model selection techniques such as likelihood ratio tests and Akaike's -40- Information Criterion have proven useful in selecting the most parsimonious model (Burnham and Anderson 1992). One possible solution, therefore, is to incorporate these techniques into analyses of count data that use likelihood-based estimation techniques (e.g., loglinear modeling) to determine whether inter-observer differences should be modeled on a case-by-case basis. Another solution is offered by an extension to the estimating equations approach to trend estimation, in which observer differences may be modeled using relatively few parameters, with a corresponding increase in efficiency (Link and Sauer 1994). A third possibility is to apply the BBS field technique in areas where more intensive studies are taking place in an attempt to quantify observer differences and, more generally, to estimate the effect of variation in detectability. Lastly, simulation studies can be used to determine the robustness of the methods of analysis to inter-observer variation and other systematic sources of error (e.g., Geissler and Link 1988). Missing data Missing data occur when sites are not surveyed in a particular year, when sites are discontinued or new ones initiated, or when the location of a site changes. Such occurrences are a ubiquitous feature of extensive volunteer-based surveys: in the BBS approximately 35% of sites are missed each year. Complete time series are not required to estimate population trends, although missing data does decrease the precision of the estimate. For this reason many published trend analyses do not include sites that are rarely or irregularly surveyed, although criteria used to exclude sites varies widely (e.g., Geissler and Sauer 1990, Bohning-Gaese et al. 1993, James et al. 1996). A complimentary technique used to increase overall precision is to weight site trends by their precision when averaging them to produce regional trends ("precision weighting"). In this way, site trends estimated with low precision have little influence on the overall population trend. The calculation of the weightings should not be based upon the estimated variance of the site trends, because the autocorrelation of counts between years makes this estimate biased (James et al. 1990). Methods of analysis that use precision weightings (Tables 3.1 and 3.2) base the weightings on statistics that are proportional to the expected precision, calculated from the number and distribution of surveys. -41 - For annual indices, imbalance in the data structure caused by missing observations is a potential source of bias and must be accounted for (Geissler and Noon 1981, James et al. 1990, ter Braak et al. 1994). A number of approaches have been taken, most of which involve either implicitly or explicitly replacing missing observations with predicted values (Table 3.1). The chaining method is an exception in that it creates balance by excluding observations. This further decreases precision of chaining indices and is one reason why that method is not recommended (Geissler and Noon 1981, Mountford 1985, Peterjohn et al. 1995). Missing data and measurement error combine to obscure the true pattern of population change over space and time, and thus make the selection of an appropriate statistical model of counts difficult. A number of methods assume that counts are log-normally distributed (Tables 3.1 and 3.2). This simplifies the calculation of estimates, but has a number of disadvantages, including the implicit assumption that counts are continuous (i.e., non-integer), and the need to add a constant before log-transformation (Collins 1990, Link and Sauer 1994). These problems are avoided by loglinear modeling (Tables 3.1 and 3.2), which assumes that counts are a Poisson random variable. A l l of these methods may, however, be biased if the model of counts is incorrect (e.g., Geissler and Link 1988). Because of this a number of trend estimation methods use non-parametric estimators, which do not require specification of the distribution of counts (Table 3.2), but are usually less efficient where a parametric model would be justified. Evaluating the methods using Monte Carlo simulations From the above discussion it is clear that there are a number of unresolved issues that prevent a complete evaluation of the relative merits of the methods of analysis. One approach is to use Monte Carlo simulations (Manly 1991) to generate count data containing known population changes, and to use these data to evaluate the methods of analysis. While the true patterns of population change and measurement error over space and time are not known, a number of plausible hypotheses can be used to construct realistic count data. These data can then be used to determine the accuracy of the methods under a variety of models for population change. A n example of the approach is given below. -42- Partners in Flight, the multi-agency coalition responsible for initiating the Neotropical Migratory Bird Conservation Program (Finch and Stangel 1993), has suggested that continental monitoring programs such as the BBS should be able to detect a 50% decline in the population of a species over 25 years with 90% power (Butcher et al. 1993). The purpose of the simulation study described here was to evaluate the bias and efficiency of three BBS trend-estimation methods in the context of linear-multiplicative declines of this magnitude and under three different geographic scenarios of population change. The methods evaluated were the same as those compared in Chapter 2: the U.S. National Biological Service (USNBS) implementation of linear route regression (Geissler and Sauer 1990), the Canadian Wildlife Service (CWS) implementation of the same approach (Erskine et al. 1992) and rank-trends analysis (Titus 1990). The USNBS and CWS co-administer the BBS program, and results from their analyses are widely cited in secondary publications. The methods of analyses are very similar: both logtransform the count data, calculate trends for each site (route) as the slope of a linear regression of log count against time (hence "linear route regression"), and take a weighted average of these slopes to produce an overall trend for the population of interest. However the implementations differ in many details, such as the criteria for the exclusion of data, the constant that is added before log-transformation, the treatment of observer differences, the way that site trends are averaged to form the population trend (arithmetic mean in the USNBS implementation, geometric mean in the CWS version), the method of back-transforming the population trend from the log scale to the arithmetic scale, the method of calculating area, abundance and precision weightings, and the technique used to estimate confidence intervals. Rank-trends analysis is a nonparametric method that estimates the direction and consistency (but not the magnitude) of population change over time at each site, and averages these statistics over sites to produce an estimate of population change. No weightings are used, and no attempt is made to account for observer differences. In Chapter 2,1 found that these three methods produced different results when they were applied to the same subset of BBS data; however, because there was no baseline of known trends, I could not tell which method was the most accurate. -43- Wilcove and Terborgh (1984) have proposed that species undergoing population declines will show one of three broad geographic patterns of population change: (1) declines in abundance at the center of the species geographic range but little contraction in the size of the range, (2) range contractions with little change in abundance at the center, and (3) both declines at the center and range contraction. These scenarios were modeled by assuming that species abundance was proportional to the normal distribution with respect to distance from the center of the range (Brown 1984, Lacy and Bock 1986, Lawton 1993), as shown in Fig. 3.2. Population declines were simulated by reducing the overall abundance by 2.73% per year (i.e., 50% over the 25 years of the survey, Fig 3.2). Annual variation in addition to trend was modeled by multiplying the overall abundance by a random log-normal variable with mean (on the log scale) of zero and variance set by a model parameter. The actual count at a site was calculated as the abundance at that site, given its distance from the center of the species range and the year, multiplied by variables representing random within-site variation and measurement error due to inter-observer differences, and rounded to the nearest integer value. Within-site variation was modeled as a random variable in a manner similar to annual variation. Each observer was assigned an observer error from a random lognormal distribution with a mean that was proportional to the year that the observer began surveying the route and variance set by a model parameter. This allowed the modeling of increases or decreases in inter-observer quality over time (Sauer et al. 1994b). The locations of sample sites were selected at random from the circular sample space centered around the center of the species range and with a radius such that 99.9% of the initial population was contained in the space. The pattern of surveillance at each site was simulated using three parameters: the probability that a new observer adopts the site, the probability of the observer surveying the site in a given year and the probability of the observer never surveying the site again. In a full simulation study, data would be generated at a number of levels for each parameter. Here, for the purposes of illustration, results are presented for each geographic scenario of population change at one level of the other parameters (Table 3.3). These parameter values were chosen to represent median values obtained from an analysis of BBS data for British -44- Columbia. For each scenario, 100 sets of data were generated and analyzed by the three methods. The main features of the results were: (1) Both route regression methods were efficient in "detecting" the declines (i.e., high proportion of significantly declining trends) in all three scenarios (Fig. 3.3). O'Connor (1992b) suggested that these methods may be less sensitive to declines that are concentrated at the edge of species ranges, due to the abundance weightings used. No evidence for this was found in these simulations. (2) Rank-trends analysis produced inaccurate results when declines were concentrated at the center of the species range (Fig. 3.3). Observer differences were not accounted for in the rank-trends method, but increasing observer quality was a major component of this simulation, resulting in increasing counts on some routes. Rank-trend estimates for scenarios involving range contraction were "correct" due to the rapid decline in abundance at peripheral sites outweighing the tendency for increases due to observer differences. Since route trends are not weighted by abundance in this method, trends at peripheral sites were given equal weight to those at the center. These results demonstrate the potential importance of accounting for observer differences, and the effect of weightings on the overall estimates. (3) The CWS route regression method showed a small negative bias in all three scenarios (median -0.82% year" ), while the USNBS method was not significantly biased (Table 3.4). The 1 bias may be due to the use of the geometric mean in estimating overall population trends from route trends rather than the arithmetic mean, although there are many other differences between the methods (Chapter 2). (4) Both route regression methods underestimated the standard deviation of the trend estimates (Table 3.4). This implies that the frequency of type I errors (significant population trends detected when none are occurring) may be higher than the nominal level of a = 0.05, and suggests that estimates of the number of species adequately sampled by the BBS based upon power analysis (e.g., Peterjohn et al. 1995) may be overestimates. Simulations under the null hypothesis of no population change are required to evaluate this possibility. -45- In conclusion, this example has demonstrated the utility of simulation studies in highlighting potential areas of bias in the methods under a small number of "What if..." scenarios. A more complete study, over a wider range of conditions, is required to confirm the preliminary results presented here. Summary In this chapter I considered the issues involved in the estimation of long-term population change from multi-site count data. It is not feasible to census the entire population of most bird species, so monitoring programs are based on sample surveys. In an ideal sample survey the population is divided into sample units (sites), a random sample of these units is selected, and the selected units are censused at regular intervals. When analyzing the resulting data, two subjective decisions are required: (1) the primary unit of analysis (annual indices or population trend), and (2) the assumed model of trend (linear multiplicative, polynomial multiplicative, additive or rank). Real surveys often depart from the ideal in three ways: (1) the population from which the sample is drawn is a subset of the population about which information is required (incomplete sample frame), (2) the number of individuals at each sample site is not counted accurately (measurement error), and (3) most sites in the selected sample are not censused every year (missing data). The analysis methods take different approaches to deal with these problems and it in many cases it is not possible to determine which is the best on theoretical grounds alone. Simulation studies offer a method of evaluating the analysis methods under known conditions. I give an example where three analysis methods are tested using three geographic scenarios of population decline: declines at the center of the range, declines at the periphery and a mixture of the two. One of the methods, rank-trends analysis, produced inaccurate results in one scenario. Another, Canadian Wildlife Service route regression, was slightly negatively biased. Both this method and US National Biological Service route regression underestimated the variance of their estimates of trend. Although the example is simplistic, it demonstrates the potential for the use of simulation studies to test the methods under realistic scenarios of population change. -46- TO G a T3 o cfci Cd G + a 'oo O CJ G T3 <+H o td CJ TO c3 P S -3 > B o ;G k CJ TO C8 < H TO "El CO TO TO CO CJ H-» O a CJ CJ TO "TO OH T3 o a o o > \rj co Go cd <D O cd o + o 2 a | .5 .a T3 £ bO o CJ X id a g D 91 G S S B) > <+H -c» CO CJ CO <D +-» O CO CJ < « G o o CJ o TO TO td O H a O T3 CJ bB G VO *—I 1—I + + bO bfi o o o o G T3 O O CJ OH CO _ IH X CJ CO O o eco o o oo G O O bO G td c3 OH CJ G o G o %a CJ o o G CJ I G OH CJ <U. (D CJ CJ co co <D • £ * <L) G ^o t1o = 1 o -1 a T3 +-> O o" G m ~ cd cd ~G G H _G o CJ CJ T<XL3> O O -•-» 3 CO co o CO 'c>> o CJ d Id *3 _G c CJ cod Id G •G G rO OH <+H 'cd H-» •3 a cal[culat o OH O co cj id 8 «» CJCcdO CJ -o •a <D o O oun .2 o G T3 VH 53 X> cj T3 CJ o o O o o G CO CO £ <-s CJ I/IS <u . G S G CO O G <+H G b yee J B +-» -O o" G CO ccl CI > CJ _> CO ce CJ Cd & a 2S OH G ISB £82 >. 2 » -i cj G O co -.ts cd * a G _G ^ CJ bO T 3 O CJ JJ imp] cd b p <D G aCO co* ^ bO G T3 o G CO CO CJ o G -O CJ CJ t CJ CO o o -47- X5 o bo G td WH o & OH CO ^ « 'S X> td if G o o • £ ^^ .a K > CO CJ CD > T3 co CD cn 2 "2 .2 -3 CD O - £ 3 43 ° 2 P ft + 3 co co <U _ 'co CD cd cd ra cd ~c3 3 cd O o o cd ~P 42 S3 PH <-H bo a J3 50 CD P '43 43 CD 3 ° 2 cd a 2 ° S3 OH 3 CO T3 B CD o t3 * q=j o ft CD bo cd P P4 S3 p o o 1 co CD co l *2 co T3 CJ_j -S3 'o 2 « CO »H o O o l ft O a O CD ft <-s T3 T3 CD IS S3 i T3 . - t J-Cd .22 "co CO CD ft S3 O 43 O ID > CD > CD CD T3 S3 H-» O S3 o S3 CO CO O ft T3 cd -f-> O S3 ft CD S3 CD H CD co >,43 3 t> o o 'co 43 §•.& S3 P S3 s ^ rH in c > § cd ^ CD CD ft CO -*-H CD O S3 S3 O S3 O S3 4i a o S3 o S3 P o cd CD bfl CD S3 CO <+H g ° * GO 1 o CD > CD cd bO S3 p T3 o o cd S-l CD CO ^ P CD CD 3o T3 CD 'o ft CD CO H-» O S3 o S3 Cd 6 0 "> (-H O S3 co 2 •3 a 4H^ 03 M O S3-, id -a CD S3 9> > ' 5 O" S3 o bo rij <L> CO CD S3 cd CD > CO CD a h-1 CD o ft CD H-» 'co H-H o P CO S3 CD o cd o .CD CD T3 CD *£ CD S3 ft o co ^o .22 o S3 PH -a Cd -t-» O •i3 a CD CO T3 O T3 O 4) co CD o o CO CD 43 cn m CN M M " CD CO CO CD CD O o 43 S3 43 a 1? 8 *S co J3 '3 bo a ^ O O o "5 •a CD 42 42 >, |H CD 13 -P S3 Cd CO o o + + |H Ed S3 O S3 /—s S3 b o cd CD CD CD CO CD co H-> a H-> 8 43 I o O O • - H PH O PH -2 cd O CM -48- . P H-» ^ il > co CD Table 3.3. Parameters used in the simulation study of the performance of three BBS trendestimation methods in the context of overall population declines and under three different geographic scenarios of population change. Where three values are given, the first refers to the scenario declines at center of range, the second to range contraction and the third to both declines at the center and range contraction (see Fig. 2). Parameter Value Abundance in year 0 at center of range 50 Variance of range in year 0 (km ) 90000 2 Rate of change of abundance at center of range (% year" ) -2.73/0/-1.37 Rate of change of range variance (% year" ) 0/-2.73/-1.37 Year-to-year fluctuations about trend (variance, log scale) 0.16 Within-site variation in counts (variance, log scale) 0.04 Rate of change of inter-observer mean (log scale) 0.05 Inter-observer variation (variance, log scale) 0.16 Site density (sites km ) 0.0005 Probability that a new observer adopts a site 0.1 Probability of an observer surveying a site in a given year 0.8 Probability of the observer never surveying the site again 0.2 1 1 2 -49- Table 3.4. Bias and precision of two trend estimators (USNBS and CWS linear routeregression) from 100 simulations of Breeding Bird Survey data at each of three geographic scenarios of population decline. 1 Scenario Median bias 2,3 1 Median estimated std. dev. ' (Actual std. dev.) 2 3 2 USNBS 1 CWS 1 USNBS 1 CWS 1 Declines at center 0.27 n.s. 0.64 *** 0.54 *** (1.45) 0.60 *** (1.28) Range contraction -0.14 n.s. I i\ *** 0 27 *** (1.29) 0.22 *** (1.44) Both -0.21 n.s. 0.86 *** 0.60 *** (1.36) 0.68 *** (1.40) 1 2 3 U S N B S = U.S. National Biological Service; CWS = Canadian Wildlife Service Units are % rate of annual change Asterisks indicate the statistical significance of Wilcoxon signed-rank sum tests of the hypothesis that the median bias is 0 (Median bias column) or that the median estimated standard deviation is equal to the actual standard deviation calculated from the 100 estimates (Median estimated std. dev. column), n.s. = p > 0.05; *** = p < 0.001. -50- Figure legends Figure 3.1. Four models of trend applied to the same hypothetical count data. Figure 3.2. Three geographic scenarios of population decline. In each case the total population declines by 50% over 25 years. Curves show the distribution of abundance across a cross section of the species range at successive 5 year intervals, with the outer top line being year 0 and the inner bottom line being year 25. Figure 3.3. The direction and statistical significance of trends estimated by three methods of trend analysis from 100 simulations of BBS data. USNBS = U.S. National Biological Service linear route-regression method; CWS = Canadian Wildlife Service linear routeregression method; RT = rank-trends method. -51 - Figure 3.2 -53 - Figure 3.3 Chapter 4 Four Case Studies Introduction In Chapter 2,1 showed that the way in which Breeding Bird Survey (BBS) data are selected and analyzed can have important consequences for the inferences drawn about regional population trends. It is therefore important to choose the method of analysis carefully, and to evaluate the results cautiously. Here I provide a detailed empirical evaluation of the results when eight candidate analysis methods are applied to continent-wide BBS data for four species. These four case studies allow the performance of many methods to be discussed at a level of detail not possible when all species are being analyzed. Previous case studies (Geissler and Noon 1981; Mountford 1982; Sauer and Droege 1990c; Peach and Baillie 1994) and comparisons of results across methods based on many species (Link and Sauer 1994; Peterjohn, Sauer and Link 1996 Appendix 1; Chapter 2) have included only a small subset of the available methods. The analysis methods are listed in Table 4.1, and are described in Chapter 1 and Tables 3.1 and 3.2. The six route regression (RR) methods were chosen because they have all been advocated in recent analyses of BBS data (refs. in Table 4.1). The other two approaches, loglinear modeling and the Mountford moving windows method, have been applied to other extensive bird surveys, and have been recommended for general use. A secondary aim of this chapter, therefore, is to investigate the feasibility of applying them to the BBS. The four bird species were chosen to represent extremes with respect to abundance (average count per route) and sample size (number of routes). Trends for low abundance species may be poorly estimated by methods that log-transform the data before analysis (e.g., least squares RR, Table 4.1), because of the need to add a constant before log-transformation (Geissler and Link 1988; Collins 1990; Link and Sauer 1994; Chapters 1,2 and 3). High abundance species can show large variation in counts between sites and time periods and may therefore be -55- poorly modeled by some estimators. Differences among methods may be more apparent in species that are seen on few routes and therefore have a small sample size (Chapter 2), but widespread species that occur on many routes will likely show greater regional variation in trends. I therefore chose to analyze data for the following species: (1) the Cerulean Warbler (Dendroica ceruled), an eastern-forest breeding Neotropical migrant that is seen on relatively few routes and in low numbers, and is thought to be sensitive to forest fragmentation and be in decline (Ehrlich et. al. 1988; Robbins et al. 1992; but see Villard and Maurer 1996); (2) the Carolina Wren (Thyrothorus ludovicianus), a widespread and abundant species, resident in the eastern and central US, that appears to have undergone a population crash in the harsh winter of 1976 (Robbins et al. 1986) and may then have recovered (Barker and Sauer 1992); (3) the Chestnut-collared Longspur (Calcarius ornatus), a central prairie dwelling shortdistance migrant that occurs in relatively large numbers but on few BBS routes and may be declining along with many other grassland species due to overgrazing and modern farming practices (Johnson and Schwartz 1993); and (4) the Belted Kingfisher (Ceryle alcyon), a pan-continental but low abundance shortdistance migrant that inhabits riparian areas and is not thought to be of conservation concern. For all of these species we have a vague idea of the broad-scale patterns of population change, but the actual levels of change are unknown. Because the "truth" is unknown, it is not possible to tell which analysis method produces the most accurate results (simulations are • required for this, Chapter 3). Nevertheless, as we shall see, many of the methods give results for some species that are not feasible. Detailed case studies are a useful way to reveal these problems and so to make preliminary judgements about the reliability of the methods. I present the results of the Cerulean Warbler first and in some detail, and then treat the other species in successively less detail as common patterns emerge. I discuss the performance of the methods and, based on this and the previous chapters, suggest how they may be improved. -56- Methods Data selection A l l methods require some amount of data screening to block out poorly performed surveys and to ensure that the data are not too sparse to allow the estimator to work. Selection criteria vary widely among published analyses, from stringent criteria that remove all but the best surveys and exclude all routes except those with a nearly continuous record of surveys (e.g., James et al. 1996) to tolerant criteria that remove only the most questionable surveys and exclude only those routes with too few surveys to allow population change to be estimated (e.g., Link and Sauer 1994). Even small differences in data selection criteria can cause important differences in results (Chapter 2). In this analysis I used two sets of criteria: (1) stringent selection criteria that excluded many routes and surveys, but allowed the estimates from (almost) all of the methods be based on exactly the same data, and (2) tolerant criteria that excluded few routes and surveys, but resulted in some of the estimators not being able to use all of the data. The tolerant criteria were used to examine the effect of a more representative amount of missing data (-35% missing surveys per route on average, as opposed to only -10% after stringent data selection). However, differences among estimates based on tolerantly selected data will be partly caused by differences in the data that each method could use (see Results). The stringent criteria comprised choosing only data that passed the selection criteria of all analysis methods except Canadian Wildlife Service (CWS) R R and Mountford moving windows (Table 4.1). CWS R R selection criteria were not included because the appropriate information on survey quality (subroute designation, Chapter 2) is not available outside of Canada. Nevertheless, the CWS R R estimator was able to utilize all of the data that passed the stringent selection criteria, so results are directly comparable with those of the other estimators. The Mountford moving windows estimator requires that a route has at least two surveys in a "window" (block of successive years) for that window to be analyzed. With a window size of 4 (the smallest possible), this means that there must be almost no missing data for all windows on a route to be analyzed. Very few routes qualify under this criterion (2.8% of the total B B S -57- routes). However most of the routes that pass the stringent selection criteria (Table 4.1) have only one or two windows that could not be analyzed. Thus the Mountford estimator with window size 4 used slightly less data than the other estimators (approx. 3% fewer surveys). With window size 6 and greater, the difference in data used was trivial. The tolerant criteria were the same as those for the U.S. National Biological Service (USNBS) R R method (Table 4.1). This resulted in a much larger sample size and a more realistic amount of missing data, but also meant that some estimators (particularly nonlinear semiparametric R R (NSRR) and Mountford moving windows with small window sizes) could not use many of the routes in the dataset. For both criteria the initial dataset consisted of all BBS data collected between 1966 and 1995, excluding routes located in some peripheral state-strata and routes placed non-randomly or off roads (Appendix 2). For three species, data from the first one or two years could not be analyzed using T R I M model 5 because there were no surveys in some physiographic strata in those years (see Loglinear modeling in Data Analysis section). In these cases the dataset was reduced to 1967 or 1968 to 1995 before data selection. Data analysis For each species, I began with a simple exploratory analysis of the data. Many of the estimators assume that the data are log-normally distributed, implying that the standard deviation of counts is proportional to their expected value. Others assume a Poisson distribution, or they square-root transform the data, assuming that the variance of counts is equal to the expected value. As an initial check, I calculated the mean, standard deviation/mean (coefficient of variation) and variance/mean on each route. If the either measure of variation divided by the mean produces nearly constant values, that measure is proportional to the mean. These analyses are preliminary because variation among counts is also caused by covariables (such as time and observers; see Loglinear modeling subsection in Discussion). To get an initial impression of the change in counts over time I calculated the mean count in each year. I also calculated the area weighted mean count, which corrected for geographic variation in sample intensity by weighting -58- each count by a measure of the area represented by the route (area of the physiographic stratum within the state or province divided by the number of routes). I then estimated annual indices and population trends using each of the eight analysis methods. These methods are discussed in Chapters 1, 2 and 3, and references are given in Table 4.1. The programs used to implement the methods are listed in Appendix 1. Below I briefly describe how each method was applied to the data in this analysis. Descriptions of the loglinear modeling and Mountford moving window methods are given in more detail because I used slightly different methods from those previously applied to bird count data. Linear route regression: In linear RR, linear-multiplicative trends are calculated for each route, and a weighted average of these route trends is used to estimate the overall trend for the region of interest. For each of the three linear R R methods (Table 5.1)1 calculated the overall trend, and tested whether this was statistically different from zero trend at an a level of 0.05 using z-tests. (Note that the implementation of CWS R R used here was slightly different from that used in Chapter 2, reflecting recent changes in practice (Downes and Collins 1996). The abundance weight (Component 6, Table 2.1) in this case was the geometric mean of the route counts.) I calculated annual indices of abundance using the residual methods for USNBS R R (Sauer and Geissler 1990) and CWS R R (Collins 1994), but not for USNBS estimating equations R R because appropriate methods of annual index calculation have not been published for this method. For comparison with the non-linear R R estimators I also calculated the proportion of routes on which populations were estimated to be declining, and estimated confidence intervals (CIs) for this proportion using the normal approximation to the binomial distribution (Zar 1996). Nonlinear route regression: In both smoothing methods (nonlinear nonparametric R R (NNRR) and NSRR, Table 4.1), counts are first smoothed at the route level and then the smoothed counts are aggregated into strata. These stratum averages are smoothed again before being aggregated (with area weightings) to form overall trend estimates. As with the original publications, (James et al. 1990, 1992, 1996) I used smoothing parameters of 0.5 for the route smoothing and 0.25 for the stratum smoothing. To test statistically for overall change in -59- abundance, I used a z-test to compare the mean of the estimated indices for the fifth, sixth and seventh years with the mean estimated indices for the fifth, sixth and seventh from last years in each route (James et al. 1990, 1996; e.g., for 1966-95 data the comparison would be mean estimates from 1970, 71 and 72 with those from 1989, 90, and 91). For the rank trend method, I calculated a trend statistic D is calculated for each route, and summed these over all routes. I compared the total D with its expected value using a z-test, based on an expectation of no overall increase or decrease in counts over time (Chapter 2). Loglinear modeling: I implemented this approach using program T R I M (TRends and Indices for Monitoring data, Pannekoek and van Strien 1996). T R I M allows five models for population change: (1) no change over time, (2) loglinear trend, (3) annual indices, (4) separate loglinear trend for each level of a categorical covariable, and (5) separate set of annual indices for each level of a categorical covariable. A l l models allow for differences in expected counts among routes. In addition, there are options to incorporate extra-Poisson variance (overdispersion) and first-order serial correlation (autocorrelation) between counts on a route into the models. Because of the exploratory nature of this analysis, I attempted to fit each of the five models both with and without overdispersion and serial-correlation. Instead of using the estimating equation method for fitting overdispersed models provided in TRIM, I estimated overdispersion for each model as Pearson's chi-squared goodness-of-fit statistic divided by the degrees of freedom (McCullagh and Nelder 1989; Lebreton et al. 1992; Anderson et al. 1994) because this is computationally more efficient. I then incorporated overdispersion into the parameter estimates by multiplying their estimated variance by the overdispersion parameter. For models 4 and 5 I used a nested set of 3 covariables, representing increasingly fine geographic scales: BBS region (Eastern, Mid, Western), physiographic stratum (up to 99 strata, Fig 1.1) and route. (Note that model 5 could not be fit using route as a covariable, because of missing counts in many years on most routes.) In all models, I used area weighting (as described for mean counts, above) to correct for geographic variation in the density of routes. In addition to model-based indices of abundance, program T R I M also reports "imputed indices", which are similar in spirit to the residual method annual indices of the linear R R -60- methods. Imputed indices are calculated as the mean count in each year, with missing data replaced by the predicted count estimated from the model. For comparison with the linear R R indices, I calculated imputed indices for TRIM models 2 and 4. Variation among observers in the number of birds counted can be an important factor in influencing counts at a site (Sauer et al. 1994b; Chapter 2). I therefore repeated the above modeling in program T R I M but assigning a new route number to each observer on each route. In a few cases for rare species some observers on a route saw no birds during their surveys. Routes with all zero counts cannot be analyzed by TRIM, so I set these counts to 1 for this method. This had a negligible effect on the results, while ensuring that the "with observer effects" and "without observer effects" analyses used exactly the same number of observations. In some cases the numerical algorithm failed to converge, especially for models that included serial correlation (these use a different fitting algorithm based on estimating equations) and models that had a moderate number of estimable parameters relative to the number of observations. In these cases I used the last calculated deviance value for model selection if the algorithm appeared to be moving towards convergence (steadily decreasing difference in deviance over successive iterations) and the last deviance was less than 1% different from the previous value. I selected the best fitting model using Akaike's Information Criterion (AIC, Sakamoto et al. 1986) rather than the Wald tests produced by TRIM, because AIC values can be compared for non-nested models (Burnham and Anderson 1992). If the largest-parameter model contained significant overdispersion, as indicated by a low p-value in the Pearson goodness-of-fit test, then I selected models using the quasi-likelihood AIC (QAIC: Lebreton et al. 1992; Anderson et al. 1994). To calculate QAIC for each sub-model, I used the estimate of overdispersion from the largest-parameter model (McCullagh and Nelder 1989). In most cases (see Results), the final model that was selected was not dependent on whether AIC, QAIC or Wald tests (where applicable) were used. Mountford moving-window method: To implement the moving-window approach, I selected the appropriate window size using the goodness-of-fit test described by Mountford -61 - (1982). I started with a window size of 4 (the minimum possible) and increased window size in steps of 2 until one or more windows failed the test, at an alpha-level of 0.05. (This procedure is identical to that recommended by Peach and Baillie (1994), except that they started with window size 6). To quantify the uncertainty in the yearly indices at each window size, I conducted 1000 bootstrap replications (resampling routes), and calculated pointwise 95% CIs for each yearly index using the percentile method (Buckland 1984; Peach and Baillie 1994). I also calculated linear-multiplicative trend estimates, using linear regression of the Mountford indices on the log scale, and back-transforming by the method of Bradu and Mundlac (1970; see USNBS R R methods, Chapter 1). I calculated CIs for the trend estimate from the bootstrap replicates. As with the loglinear modeling, I calculated Mountford indices both with and without observer differences (i.e., coding each new observer as a new route). Comparison of estimates Quantitative comparisons among estimators are difficult because they estimate different parameters (annual population index, loglinear trend or smoothed trend) at different scales (logscale or untransformed). In addition, statistical tests of differences among estimates for each species are not appropriate because the estimates are based on the same (or similar) data and so are not independent, but their covariance is unknown. I therefore focus on graphical displays of the estimates to allow qualitative comparisons. I also compare loglinear trend estimates and their corresponding CIs, for those estimators where this is possible. For route-regression estimators I compare the proportion of routes with estimated declines. To enable comparison of indices among estimators, graphs of population change over time are scaled so that the index for the first year is 1.0. The only exceptions are annual indices of abundance from the linear R R estimators, which are properly displayed as residuals about the linear trend line. Because population change is thought to be a multiplicative process, results are presented using a log-scaled axis for the indices: this gives equal emphasis to proportional increases as decreases (Moses and Rabinowitz 1990). I used the log base 2 because this makes the axes easy to interpret: a change in index from 1 to 2 indicates an estimated doubling of -62- populations and a change from 1 to 0.5 indicates a halving; both are equally spaced along a log 2 axis. Here I focus exclusively on estimates of population change at the continental level. Regional estimates are likely to have similar properties to the continental estimates for the two species that were seen on few routes. For comparability I use the same models for each estimator for both the tolerant and strict data selection criteria. Results Cerulean Warbler Stringent data selection: Between 1966 and 1995, Cerulean Warblers were counted on 388 routes, but only 41 of these passed the stringent selection criteria (Table 4.2). The distribution of counts among these routes was highly skewed, with most routes having a median count of zero, while two, 82031 (Elk Valley, Tennessee) and 82032 (Smokey Junction, Tennessee), contained 31% and 11% of the total counts respectively (Fig. 4.1). The coefficient of variation (standard deviation of counts within routes/mean route count) decreased with increasing mean (Fig 4.2a, left plot), implying that counts are likely not approximately lognormally distributed. The pattern was reminiscent of a simple inverse relationship, but this would imply that the standard deviation was independent of the mean when in fact there was a significant positive relationship (least-squares linear regression slope = 0.49, 95% CI 0.52 to 0.46; the relationship was similar if the two high abundance routes were deleted). The variance/mean ratio was greater than 1.0 and increased with increasing mean count (Fig 4.2a, right plot; linear regression slope = 3.72, 95% CI 3.97 to 3.48; the slope was still significantly greater than zero if the two outlying high abundance routes were deleted), indicating possible departures from the Poisson distribution as well. The mean counts remained relatively stable until about 1977, then declined quite sharply until the mid 1980s, with some indication of a recovery in the last few years (Fig 4.3a). The area-weighted mean counts showed a similar pattern but were slightly more negative, with the decline perhaps starting after 1972 (Fig 4.3a). The means are clearly heavily influenced by -63- counts on routes 82031 and 82032: for example the increase in mean count in 1970 is an artifact of the missing data on route 82031 before that time (Fig. 4.1). A l l three linear R R methods estimated that Cerulean Warbler populations declined significantly over the period of the survey (Table 4.3). The magnitude of the estimated change differed among the methods (USNBS estimating equations < USNBS R R < CWS RR, Fig. 4.3b), although CIs for the estimates showed a large overlap (Table 4.3). The estimating equations estimate was much less precise than the other estimates, with a CI three times as wide as the comparable USNBS R R estimate (Table 4.3). Annual indices of abundance showed similar patterns for the USNBS R R and CWS R R methods (Fig 4.3b; annual indices for estimating equations approach have not been developed). In both cases the annual indices fluctuated around their loglinear trend line, which appeared to describe the overall pattern of population change well. The results from this methods imply that Cerulean Warblers have been declining at a reasonably constant rate for the time period of the survey. Although all three nonlinear R R methods also estimated overall declines, the estimates were not statistically significant for either smoothing R R method (Table 4.3) and there appeared to be some important changes in trend over time which reflect changes seen in the mean counts (compare Figs. 4.3c and d with 4.3a). Estimates that incorporated differences among observers (NSRR, Fig. 4.3c) were slightly more positive than those that did not (NNRR, Fig. 4.3d). In both methods, indices for the last 3 years were not joined smoothly to those from earlier years, indicating a problem with the methods. N S R R estimates were less precise than those from N N R R (compare CI widths in Table 4.3). The rank-trends method estimated highly significant declines, although the proportion of routes with estimated declines was considerably less than the other R R methods and was not significantly different from 0.5 (Table 4.4). The results of loglinear model fitting and selection using program T R I M (Table 4.5) can be summarized as follows: (1) when observer effects were included, the fit of the model to the data was adequate for all but model 1 (p for Pearson % non-significant, Table 4.5) with no 2 indication of overdispersion (Pearson % /d.f. close to 1.0, Table 4.5); (2) T R I M model 3 (separate 2 effects for each year) had the lowest AIC and is therefore the selected model; (3) modeling serial -64- correlation caused a small reduction in deviance and AIC so could be included in the final model, although the estimated correlation was small for the well-fitting models (Table 4.5). Index estimates from model 3 including observer effects (Figs. 4.3c and 4.4a) indicated that Cerulean Warblers were relatively stable or perhaps declining slowly until 1977, declined quite rapidly up to the late 1980s and then stabilized, or perhaps showed a small recovery. After 1980, all yearly indices were statistically significantly less than in the first year of the survey (upper 95% confidence limit less than the first year index of 1.0, Fig. 4.4a). Estimates with or without serial correlation were essentially identical (compare Figs. 4.4a and b), which is not surprising given the small correlation estimate (Table 4.5, upper half). If observer effects were not included the estimated decline stops in the early 1980s (Fig 4.4c), but the "without observers" model did not fit the data very well (Table 4.5). The lack of fit translated into an estimated overdispersion greater than 1.0 (Table 4.5), which resulted in slightly wider CIs when this was incorporated into the estimate of standard error (Figs. 4.4c and d). The estimated serial correlation was also nontrivial in this model (Table 4.5, lower half), which also produced wider CIs, although it had almost no effect on the estimate (compare Figs. 4.4c and d). The estimates from Model 3 "without observers" were very similar to the yearly mean counts and area weighted mean counts (compare Figs. 4.3a and d), but when observer differences were taken into consideration the indices were noticeably more negative (Fig. 4.3c). Although the data were better modeled in T R I M as annual indices, it is informative to compare linear trend estimates from program TRIM (model 2) with those from the linear R R methods (Table 4.3). When observer effects were included, the estimated magnitude of decline was between that of the USNBS estimating equations estimator and CWS RR, but with a smaller CI than for these methods. Without observer effects, the estimated decline was still statistically significant, but was of much smaller magnitude. The Mountford moving window method did not perform well with these data. When observer differences were incorporated, all windows passed the goodness-of-fit test for all years with window sizes 4 and 6 (Table 4.6, columns 2 and 3), but not with larger windows. Nevertheless, the estimated population index in 1991 with window size 6 (Fig. 4.5b) was -65- implausibly high and exceeded the CIs by a large margin. Estimates were outside of the CIs in several other years, indicating that they are probably highly biased. When observers were not included, the goodness-of-fit test gave negative % values for one window with window size 4 2 and another with window size 6 (Table 4.6). Negative % values are not admissible, and indicate 2 that one or more of the estimated counts in those windows were below zero. In addition, the estimates for window size 6 were clearly spurious for years 1978 and 1979 (Fig. 4.5d). The overall patterns of population change estimated with the Mountford movingwindows method did not match any of the other estimators (Figs 4.3c and d). The most credible estimates from the Mountford method ("with observers", window size 4: Fig 4.5a) showed population fluctuations but no consistent trend, and the Mountford trend estimates with this window size showed close to no trend with narrow CIs (Table 4.3). Tolerant data selection: When the tolerant data selection criteria were applied to the data, the number of routes available to the estimators was increased from 41 to 198 (Table 4.2). This was accompanied by an increase in the percentage of missing data, from 11% to 32% (Table 4.2). The annual mean counts showed a qualitatively similar pattern to those after stringent data selection but with smaller fluctuation (compare Figs. 4.3a and 4.6a), probably indicating a reduction in variation due to sampling error. Even with almost 200 routes, however, the indices were still noticeably influenced by a few routes with high counts: for example the increases in mean counts after 1991 disappear with the removal of two routes that were surveyed only in these later years (66080 Dell, Ohio and 90028 Strange Creek, West Virginia). Both the USNBS R R and CWS estimates were more positive using tolerantly selected data than with stringent data selection (Fig. 4.6b), and the residual annual indices showed much less variation. The USNBS estimating equations R R estimate, on the other hand, was more negative, and 95% CIs for the two USNBS estimates no longer overlapped (USNBS R R : estimate -2.16% per year, 95% CI -1.65 to -2.86; USNBS estimating equations RR: estimate 5.30% per year, 95% CI -3.20 to -7.49). No identifiable small group of routes seemed to be responsible for this difference between the methods. -66- The smoothing R R methods were not able to use all of the data that passed the tolerant selection criteria: N N R R used 177 of 198 routes for the smoothed index estimates and 125 routes for the /-tests, N S R R used 121 routes for the indices and 90 for the Mests. Both methods produced clearly spurious estimates of population change between 1966 and 1967 (Fig. 4.6c and d), because there were no suitable routes in three strata in 1966 (see Discussion). This year apart, the indices for both methods no longer showed the estimated stable period in the early years of the survey, but now estimated relatively constant declines over that period (compare Figs. 4.3c and d with 4.6c and d). The overall estimates of population change were still not statistically significant for either method (NNRR: estimate -2.42% per year, 95% CI +0.76 to 4.76; N S R R estimate -4.60% per year, 95% CI +1.26 to -8.07). The results of model fitting and selection in program T R I M differed slightly from the stringently selected dataset: models 3 and 4 (strata) "with observers" had almost the same AIC (4481.99 and 4481.68 respectively; the next closest was model 2 "with observers" 4494.79. Note that AIC values are only comparable across models fitted to the same data so these values should not be compared with those from the stringently selected data). Model 5 (strata) would likely have had an even lower AIC, but unfortunately it could not be fit because there were no surveys in several years for several strata. Index estimates from model 3 were more positive using the tolerant data selection criteria (Fig. 4.6c and d) than the stringent criteria (Fig. 4.3 c and d), particularly for the years 1979 1992 in the "with observers" model. Once again the "without observers" model estimates are very similar to the yearly mean counts (compare Figs 4.6a and d), except for the last few years. In addition both model estimates track the smoothed indices from N N R R and N S R R quite well (excluding the first year, Fig. 4.5c and d). Estimated trends using Model 2 were also more positive (-2.96% per year "with observers", -1.25% per year "without observers"; compare with Table 4.3). CIs with the larger dataset were narrower, although the estimated overdispersion in the "without observers" models was comparable to those reported in Table 4.5. The Mountford method was able to make use of some of the extra data (e.g., mean number of routes used per window, "without observers" models: 146.5 of a possible 198 for -67- window size 4, 156.9 for window size 6). The data passed the goodness-of-fit test for window sizes 4 and 6 for both "with observers" and "without observers" models, and the estimates exceeded the 95% consistency bounds on only two occasions, and then by only a small amount. Nevertheless, the estimated indices showed a similar overall pattern to those from the stringently selected data: seemingly random fluctuations that did not match estimates from any of the other methods. Summary of estimates: Excluding the Mountford method results, and the few clearly spurious estimates from the smoothing RR methods, the results for Cerulean Warbler showed overall declines of somewhere between 30% and 90% between 1966 and 1995 (Table 4.3 and tolerant selection criteria results above). Declines were either relatively continuous over that time period (Linear R R methods: Figs. 4.3b and 4.6b) or occurred largely in the late 1970s and early-to-mid 1980s (other methods: Figs 4.3c, d and 4.6c, d). Carolina Wren Stringent data selection: After stringent selection of Carolina Wren data, 268 routes remained from a possible 1435 (Table 4.2). Counts from 1966 were excluded because there were no counts in some strata in that year (see Data Selection Criteria in Methods section), leaving 29 years of data. Both the mean and median counts were much higher than for the Cerulean Warbler (Table 4.2), but the pattern of decreasing standard deviation/mean with increasing mean was similar (Fig. 4.2b, left plot). On average, the variance was appreciably larger than the mean (mean (variance/mean) = 4.33; median 3.54) indicating overdispersion relative to a Poisson distribution. The variance/mean ratio increased with increasing mean (Fig 4.2b, right plot; least-squares linear regression slope = 0.078, 95% CI 0.112 to 0.044), again indicating departure from the Poisson model. The yearly mean counts clearly demonstrated the main feature of this time series: the catastrophic effect of the cold winter in 1976 (Fig. 4.7a). Both before and after this, it appeared that counts of Carolina Wrens were increasing at a near-exponential rate (i.e., linear on the log scale). This pattern was mirrored in the residual indices from the USNBS and CWS R R methods (Fig. 4.7b). There was also some evidence of a second recent reduction in numbers -68- (1993 to 1995), although of much smaller magnitude. A l l three linear R R estimators gave almost identical estimates of trend over the whole time period (Table 4.7), although these are not particularly meaningful given the obviously non-linear changes. The smoothing R R methods poorly depicted the dynamics of this species, because the smooths obscured the interesting pattern. In both N N R R and NSRR, the decline was estimated to have been of relatively small magnitude, and to have occurred gradually over the period 1974 to 1981 (Fig 4.7c and d). The methods gave similar estimates, except that N N R R estimated a stronger recovery at the end of the time period. The nonparametric rank-trends estimate was also uninformative: it showed an overall statistically significant increase over the whole time period (Table 4.7), but provided no other information. None of the loglinear models fit the data very well (Pearson % /df = 2.41 to 4.32 and 2 goodness-of-fit p < 0.001 for all models, Table 4.8), indicating either that the counts were overdispersed relative to a Poisson distribution or that the models were incorrectly specified (see Discussion). Assuming for the moment that the lack of fit was due to overdispersion, I selected a model based on the QAIC criterion with the overdispersion parameter equal to 2.41, the estimate from the highest-dimension model for which an estimate was available (model 5 (region) "with observers"; the higher-dimension models 5 (stratum) and 4 (route) failed to converge). The model with the lowest QAIC also had the lowest AIC, so in this case the model selection process was robust to the choice of overdispersion parameter. Models that included observer effects fit considerably better than those that did not, and models including serial correlation fit slightly worse than those that did not (Table 4.8). Models with separate parameters for each time period fit better than those with one overall trend parameter, and, generally, the finer the spatial scale of the covariable the better the model fit. Model 5 (stratum) with observer effects had the lowest QAIC and AIC (Table 4.8). Unfortunately, parameter estimates were not available for this model because the fitting algorithm did not reach convergence. I therefore selected the next best model, Model 5 (region), as the model to use for parameter estimates (Figs. 4.7c and d). Population indices from this -69- method were generally similar to the residual linear R R indices, although the amplitude of the population fluctuations (peaks in 1976 and 1992 and troughs in 1977 and 1994) were less. The Mountford moving window estimates were again problematic. When observer differences were ignored all windows passed the goodness-of-fit test for window sizes 4, 6 and 8, and one window narrowly failed the test at window size 10. Although population indices for all windows showed the decline from 1976 to 1977, the estimated size of the decline increased with window size, as did the estimated rate of recovery (Fig. 4.8). The recovery rates were implausibly fast for window sizes 8 and 10. Bootstrap CIs were suspiciously narrow for all window sizes, indicating that important sources of variation were likely being ignored. When observer differences were included, one window narrowly failed the goodness-of-fit test at window sizes 4, 8, and 10, but the index estimates were very similar to the "with observer" indices. In both cases using window size 6 gave the most plausible index estimates, but these indices still showed a relatively flat trend after the mid 1980s when all other methods indicated an increase in abundance (Figs. 4.7c and d). Tolerant data selection: Population change estimates for all methods were almost identical to those described for the stringent data selection criteria (not shown). The weighted annual means were slightly more positive towards the end of the time series, and this was reflected in more positive overall trend estimates (approx. 0.5% per year) for all of the estimators except the Mountford method. Summary of estimates: The well-known crash in Carolina Wren populations between 1976 and 1977 was documented in all of the estimators except the loess smoothing methods, although estimates of the size of the decline varied (approximately 40% reduction for Mountford method, 50% for T R I M model 5 (region), and 70% for linear R R methods, Fig. 4.7). A l l of the estimators except the Mountford method indicated that Wren populations were increasing at a near-exponential rate both before and after the crash (Fig. 4.7); there was also some suggestion of a recent set-back in numbers between 1992 and 1994 (Fig. 4.7b). -70- Chestnut-collared Longspur Stringent selection criteria: For the Chestnut-collared Longspur, only 21 of a possible 190 routes remained after stringent data selection. There was a very large variation in the number of birds counted per survey, from 0 to 340, and the average counts were much higher than for the other species analyzed (Table 4.2). Nevertheless, the patterns of declining standard deviation/mean and increasing variance/mean with increasing mean per route appeared to be similar to those of the other species (Fig. 4.2; not investigated quantitatively because of the small sample size). The annual mean counts fluctuated around the initial value until 1982, after which there was some suggestion of a step to a lower level (Fig. 4.9a). This pattern was reflected to some extent in the residual linear R R indices (Fig. 4.9b), although an increase in indices in the mid 1970s gave the impression of little overall change. Trend estimates were almost identical for the three linear R R methods: all showed a non-significant decline of approximately 0.5% per year with relatively similar CIs (Table 4.9). The nonlinear R R estimates were more negative, with the N N R R and rank-trend estimates being statistically significant (Table 4.9). These were also the only R R methods where the proportion of routes with estimated declines was significantly different from 0.5 (Table 4.4). The estimated indices for N N R R (the "without observers" model) matched the overall pattern in the mean counts (compare Figs. 4.9a and d), while the estimates for N S R R were somewhat similar to the pattern in the residual linear R R indices (compare Figs. 4.9b and c). The loglinear models fit the data very poorly, with an estimated overdispersion for the highest parameter, best fitting model of 7.02 (all models failed the Pearson goodness-of-fit test withp < 0.001). Using AIC as the model selection criterion, the best model was model 5 (stratum) "with observers" (Deviance 2618.95, 139 parameters, AIC 2896.95, QAIC 651.24). However, when overdispersion was taken into account (i.e., using QAIC), the lower-dimension model 4 (route) "with observers" was favored (Deviance 3056.89, 76 parameters, AIC 3208.89, QAIC 587.65). As with previous species, models that included observer effects invariably had lower AIC and QAIC. Including serial correlation in the models worsened the fit in all cases. -71 - The index estimates for model 4 (route) "with observers" showed an initial decline, but then an increasing trend over the period of the survey (Fig 4.9c). This was quite different to the estimates from the linear R R methods, which also fit individual loglinear trends to each route (but used a different method to combine these trends into regional estimates). The corresponding "without observers" model estimates showed an increasingly shallow declining trend, closer to the estimates from the other methods (Fig. 4.9d). Imputed indices for model 4 (route) may be expected to show similar variation to the residual indices of abundance in the linear R R methods, especially for the "with observers" model. However, the method that I used to incorporate observer differences (coding each new observer as a new site) created many missing cells in the data matrix, and the use of model estimates in place of this "missing data" resulted in imputed indices that were very similar to the model indices (Fig 4.9c). The imputed indices in the "without observers" model (Fig. 4.9d) were similar to the linear R R residual indices and the mean counts in the latter part of the survey (after 1985), but not in the first part of the survey, where there were more missing data. The Mountford models did not pass the goodness-of-fit test for any window size from 4 to 8, giving negative % values in many cases. Population indices showed unfeasible estimated 2 population changes in some years (e.g., Fig. 4.9c and d), exceeding the 95% CIs which were themselves very wide (e.g., for window size 4 the index estimate for 1995 was 1.43 with 95% CI 0.28 to 28.12). Tolerant selection criteria: With tolerant data selection, the number of useable routes increased to 128. The annual means showed a somewhat different pattern, with an initial rise in counts followed by a large decline between 1975 and 1985 to about half the initial level (Fig. 4.10a). The CWS R R trend estimate was nearly unchanged, but the USNBS R R and estimating equations estimates were substantially higher (Fig. 4.10b) with much greater variance (see Discussion). Like the CWS R R indices, the N S R R smooth showed no overall trend (Fig 4.10c), while the N N R R smooth was similar to the annual means (Fig. 4.10d). In program T R I M , neither model 4 nor model 5 reached convergence, so parameter estimates were not -72- available. The Mountford moving window indices showed small fluctuations about the initial value, but they did not appear to correspond with either the annual means or CWS indices. Summary of estimates: It is difficult to summarize the findings for this species because the methods gave so many different estimates, and the effect of data selection appeared to be strong. Both the annual means and the N N R R smoothed indices indicated that the population declined during the mid 1980s, but neither of these took account of possible differences among observers. The other methods either indicate no overall change, or gave unrealistic estimates. Species like the Chestnut-collared Longspur that show a very wide variation in counts both spatially and temporally seem to be particularly difficult to analyze. Belted Kingfisher Stringent selection criteria: Belted Kingfishers were counted on over 2000 routes throughout the area of survey coverage, but the maximum count on any one route was 15 (far fewer than the other species, Table 4.2). Surveys on the western routes were begun in 1968 so only data from 1968 to 1995 were used in this analysis. After stringent data selection 529 routes remained, with a mean count per survey of 0.58. The patterns of declining standard deviation/mean and increasing variance/mean with increasing mean were similar to the other species (Fig. 4.2d; least-squares linear regression slope = 1.67, 95% CI 1.78 to 1.56). Compared with the other species, the annual means showed relatively small fluctuations, although (Fig. 4.1 la). Overall there was perhaps some evidence of a gradual decline in mean count, although this was largely due to several high count years between 1972 and 1977. The small overall decline was reflected in the linear R R estimates, all of which were statistically significant (Table 4.10). The estimating equations estimate was noticeably lower than the other two, but all three CIs overlapped. The pattern of fluctuations about the trend line mirrored those of the mean counts, but the amplitude was considerably less (compare Figs. 4.1 l a and b). The nonlinear smoothing R R methods both showed an increase in counts up to 1973 followed by a consistent decline (Fig. 4.1 l c and d). This was different from the other methods, but the difference was mainly due to an estimated decline between the first and second years in all other methods. Because the statistical tests for trend in these methods only considered the -73- period from 1972 to 1991 (see Methods) the initial increase was excluded and both N N R R and N S R R estimated an overall significant decline (Table 4.10). The rank-trends method also estimated an overall significant decline (Table 4.10). The proportion of routes with estimated declines was similar among the R R methods, and all were significantly greater than 0.5 (Table 4.4). The outcome of the loglinear modeling was different from previous species because models without observer effects had a lower AIC and QAIC than those that accounted for observer differences (Table 4.11). The models fit the data better than for the previous two species, but the large sample size meant that all models failed the Pearson goodness-of-fit test (p < 0.001) even though the estimated overdispersion was relatively small (1.11 to 1.24). The final model included parameters for each year, but surprisingly (given the large amount of data and the wide geographic range of the species) did not include region or stratum effects. Note, however, that neither model 4 (route) nor model 5 (stratum) reached convergence and these could conceivably provide a better fit. Models with or without serial correlation were very similar. The estimated indices were very similar to the annual mean counts (Figs. 4.1 l c and d). Unlike the previous species, all windows passed the Mountford goodness-of-fit tests for window sizes from 4 to 16 for the "with observers" data and 4 to 14 "without observers". I therefore used the window size 14 for both datasets. The annual indices were very similar to those from T R I M model 3 (Fig. 5.11c and d) except that the decline in abundance between 1977 and 1979 was larger and became the major feature of the dataset. The linear trend estimates were well within the range of the other methods (Table 4.10). As with the loglinear model estimates, there were only minor differences between the "with" and "without observers" indices. Tolerant selection criteria: When data were selected using the tolerant criteria the sample size nearly tripled (Table 4.2), but the estimates remained essentially the same (Fig. 4.12). The only exception was the NSRR method, which estimated strong declines at the end of the time series (Fig 4.12c). Summary of estimates: There was generally good agreement among the methods for this species, apart from estimates for the first two years in the smoothing R R methods, and -74- estimates for the period 1977 to 1979 in the Mountford moving window method. There appears to have been a small but significant decline (-15 - 35%) in the abundance of Belted Kingfishers that either occurred mainly in the years 1977-79 (Mountford moving window method) or was relatively continuous from the mid-1970s onward (other methods). Summary of performance of methods The methods produced quantitatively similar estimates in many cases, but there were also a number of obvious problems (Table 4.12). The linear R R methods produced seemingly reasonable estimates in all but one case (the tolerantly selected Chestnut-collared Longspur data, Fig 4.10b). There were some quite differences in results among the three implementations (e.g., Fig 4.6b, CIs given in text), although the 95% CIs often overlapped indicating that high variance of estimates was a contributing factor. Annual indices based on the residual method appeared to be biased toward the trend line in some cases (Figs. 4.3b and 4.10b). The smoothing R R methods produced seemingly realistic summaries of the patterns in some cases (Fig. 4.9c and d). However, they were plagued by problems at the beginning and end of the time series (Fig. 4.3c and d, Fig. 4.6c and d). The two different implementations, N N R R and NSRR, could produce quite different estimates (Fig. 4.10 c and d, Fig. 4.12 c and d), but the 95% CIs for the overall trend always overlapped. The rank-trend method apparently suffered no problems, but it does not give enough information in the results for lack of fit to be judged to the same extent as the other methods. Loglinear modeling gave seemingly reasonable results for three of the four species, although many of the more complex models did not reach convergence. Estimates for Chestnutcollared Longspur were clearly unreliable, and the fit of the models (judged by deviance/d.f) was very poor. There were significant problems with the Mountford moving-window method, including highly unrealistic estimates (Fig. 4.5), estimates that exceeded the 95% bootstrap CIs (Fig. 4.5b and d), and difficulties in selecting an appropriate window size (Fig. 4.8). The method gave seemingly reliable results in only one case, the Belted Kingfisher (Figs. 4.11 and 4.12). -75- Discussion Empirical comparisons of estimates cannot determine which estimation method is the most accurate, because the true parameter values are unknown. Nevertheless, case studies are useful for testing expectations about the differences among methods, for exposing situations under which they produce implausible estimates, and for suggesting areas that need further development. Here I discuss each group of estimators in turn. I place the results in the context of previous work, consider what inferences can be made about the reliability of the methods, and suggest some possible improvements to the methods. Linear route regression Linear R R methods are probably the most extensively studied of the population trend estimators. Several studies have investigated the bias caused by the addition of a constant before log-transformation in both CWS and USNBS RR methods (simulations by Geissler and Link 1988; Collins 1990; Link and Sauer 1994). Adding a constant biases the route trend estimates towards zero (i.e., no trend); this effect is greater for larger values of the constant (although the reverse is not always true, Collins 1990) and for routes that have few birds counted. In Chapter 2 I showed that there was a tendency for USNBS RR estimates (0.5 added) to be closer to zero than CWS trend estimates (0.23 added), and for both estimates to become more similar as mean abundance and sample size increased. The results in this chapter show a similar pattern. The CWS estimate was larger for the low abundance, low sample size species (Cerulean Warbler) under stringent data selection criteria (Table 4.3, Fig. 4.3b), but the two became more similar when the sample size was increased (Fig. 4.6b). The estimates were very similar for the low abundance, high sample size species (Belted Kingfisher, Table 4.10, Figs. 4.1 l b and 4.12b), but again the CWS estimate was slightly larger. The estimates were also similar for the high abundance species (except for Chestnut-collared Longspur with tolerant data selection, discussed presently). The same pattern is evident when looking at individual route trends within a species: the magnitude of the trend estimate tends to be closer to zero for lowabundance routes, and the effect is slightly stronger for the USNBS method than the CWS method (compare Figs. 4.13a and b). -76- The estimating equations method was introduced because it does not require the data to be log-transformed and so bypasses the issue of what constant to add and the accompanying bias (e.g., Fig 4.13c). It also makes no assumptions about the underlying distribution of the counts other than the expected value being a multiplicative function of an observer effect and an exponential (loglinear) trend (Link and Sauer 1994, equation 1). Simulations have shown that it is slightly positively biased (Link and Sauer 1994; Barker and Sauer 1992 have suggested that the other R R estimators are also positively biased). However, comparisons with the U S N B S R R estimator (using 348 species that occurred on 40 or more BBS routes after tolerant data selection criteria) showed that the estimating equations estimates were slightly more negative, especially for high abundance species (Link and Sauer 1994). Of the four species studied in this chapter, estimating equation estimates were more negative than both USNBS and CWS R R estimates in all cases, including the species with an increasing trend (Carolina Wren). The causes of the difference in behavior between the simulated and real data are unclear. Another difference between the linear RR methods is the weightings used in combining route trends to produce regional estimates. In Chapter 2 I found that trend estimates for four of the 115 species analyzed were unrealistically large using CWS RR because of the method of calculating the abundance weightings. In all four cases, one route with a large apparent change in density but few surveys was assigned a very large abundance weighting due to the extrapolation of the route trend to calculate the least-squares mean weighting value. In this chapter, I used the geometric mean as the abundance weight in CWS R R to reflect recent changes in practice (Downes and Collins 1996). However, both USNBS R R and the estimating equations estimator use the least-squares abundance estimate, rather than geometric mean, when this does not involve extrapolation (Chapter 2). When there are few surveys close to the midyear this can produce very imprecise estimates and lead to the same effect as that noted for the former CWS weighting. This is exactly what happened for the tolerantly-selected Chestnutcollared Longspur data (Fig 4.10b): route 04008 (Bow City, Alberta) had an abundance weighting of 5939, which is more than 10 times larger than the highest count for that species. The overall trend estimates were unrealistically high (+7.22% per year for USNBS R R and -77- +3.22% per year for the estimating equations estimator) and the variances were also large (SE 8.54 and 3.72 respectively). When this route is removed, the USNBS R R estimate becomes -0.32%) per year (SE 0.61) and the estimating equations estimate becomes -0.33% per year (SE 0.44); both are closer to the CWS estimate of -0.23. There are reasonable theoretical reasons for using the least-squares mean as an abundance weight; however even in the USNBS estimators it is currently used only part of the time (when there are surveys both before and after the midyear). To prevent the sabotage of the trend estimate by poor estimates on single routes it seems necessary to either constrain (Winsorize) the least-squares mean estimate (for example, to be less than the maximum observed count), or to use the geometric mean count throughout. Another possibility is to incorporate uncertainty in the estimate of abundance weighting into the precision weight, although the use of precision weights is controversial (James et al. 1990). The calculation of annual indices from route regression estimates has received less attention. The only simulation study (Collins 1994) compared CWS residual indices with two ANOVA-based methods and concluded that an index that incorporated both year effects and route trends (non-parallel A N O V A ) showed the least bias, although the results were preliminary. Despite this, residual indices are still the accepted approach for linear R R (e.g., Peterjohn et al. 1995), and are used as informal checks on the regional trend estimate. The CWS and USNBS implementations differ in three ways: (1) CWS indices are calculated entirely on the log-scale while USNBS indices are calculated mostly on the multiplicative scale; (2) CWS indices are calculated relative to the mean survey year (i.e., mean of the years on which surveys took place), while USNBS indices use the mid-year (i.e., mid point of the time period analyzed, regardless of when the surveys took place, e.g., 1980.5 for the time period 1966 to 1995); (3) CWS indices are based on residuals from the route trend estimates, while the USNBS estimator first re-fits the regional trend to each route before calculating residuals. Despite these differences, the estimates were very similar. Both appeared to accurately portray the non-linearity in the Carolina Wren series (Fig. 4.7b), but there was a suggestion that both were biased towards the straight trend line in the Cerulean Warbler series (Fig 4.3b). In addition, the USNBS residual indices did not strongly highlight the problem with route 04008 in the Chestnut-collared Longspur data (Fig -78- 4.10b). There is clearly a need to further investigate the use of residual indices before they can confidently be used as a diagnostic to judge the fit of the linear trend estimates or to draw conclusions about patterns of population change (see also discussion by Collins 1994). Nonlinear route regression At their best, the smoothing-based RR methods (NNRR and NSRR) can produce a plausible picture of the trend underlying noisy BBS data (e.g., Fig. 4.3d, Fig. 4.9d). However, in many cases it appeared that the estimates were unreliable near the beginning and the end of the time series. The problem arises because many routes have missing data in the first or last few years (even after data selection), and the smoothing-based methods do not extrapolate index estimates beyond the limits of the available data. When smoothed route indices are averaged to produce strata means, which are then smoothed again, the strata smooths at the beginning and end of the series will be based on a subset of the routes used in the middle of the series. This leads to obvious discontinuities such as at the end of the stringently selected Cerulean Warbler time series (Fig. 4.3d), and may account for the implausible estimates at the beginning of the Belted Kingfisher series (Figs. 4.1 l c and d, and 4.12c and d; also the end of 4.12c). When there are no surveys in a particular stratum, jumps may occur, such as at the beginning of the tolerantly selected Cerulean Warbler time series (Fig 4.6c and d). Missing index estimates also cause the statistical tests for overall trend to be based on fewer data points than the smooths, which is clearly undesirable. These problems are minimized in published smooth-based estimates (James et al. 1990, 1996) by using stringent data selection criteria. They are also less likely to be important with large sample sizes (although stringent data selection criteria lead to a drastic reduction in the amount of available data). Nevertheless, given the results of the case studies, it would be useful to investigate the consequences of extrapolating the smoothed estimates to the beginning and end of the time series. The loess smoother used in N N R R and N S R R is capable of extrapolation, as are other smoothers such as the cubic smoothing spline (which is linear beyond the bounds of the data, Hastie and Tibshirani 1990). Smoothing at the route level will, however, always be datahungry. It will probably be necessary to use stringent data selection criteria to ensure that there -79- are some surveys near the beginning and end of the time series in each route and so prevent overly long extrapolations. One aspect of the NSRR estimator not discussed by James et al. (1996) is the failure of the numerical algorithm for more than 30% of the routes that would otherwise pass the stringent data selection criteria. This failure occurs because the matrix of smoothed observer residuals (X'Xres in James et al. 1996 Appendix 1) is close to singular, and is not properly inverted by the LU-decomposition and back-substitution algorithm. Rescaling the initial matrix will sometimes solve such problems (Bard 1974), but in this case the poor condition of the matrix is likely due to the distribution of observers among years and not to large differences in the matrix elements. Alternative numerical methods, such as singular value decomposition (Press et al. 1992) and other direct decomposition methods (McCullagh and Nelder 1989) cope better with i l l conditioned matrices and are worth investigating. James et al. (1996) noted that NSRR estimates were, on average more negative than N N R R estimates for 26 species of warblers using stringently selected data from eastern and central North America. This result was consistent with a detailed analysis of the differences in counts among observers in the BBS based on linear R R analyses (Sauer et al. 1994b), which showed that for many species observers in recent years tended to count more birds than those in earlier years. Sauer et al. (1994b) speculated that this may be due to increasingly stringent screening by BBS coordinators of new observers, leading on average to increased observer quality over time. They conclude that ignoring observer differences will tend to positively bias trend estimates, although the effect varied among species. The results of one test were reported on a species-by-species basis, and can be compared with the differences observed for three of the four species analyzed in this chapter (Table 4.13; the fourth species showed non-significant results in Sauer et al. 1994b and so results were not reported). The effect of incorporating observer differences into the analysis is qualitatively similar in only one of the three species. Further comparisons based on more species are required to determine whether the conclusions of Sauer et al. (1994b) about observer differences are robust to the choice of analysis method. -80- The smoothing estimators were not flexible enough to accurately portray the sudden change in abundance of Carolina Wrens between 1976 and 1977 (Figs. 4.7c and d). The estimated low point using N N R R and NSRR was 1981, 4 years after the actual crash had occurred. As implemented, these methods are clearly unsuitable for detecting large, short term changes in abundance that may be of conservation concern (Peterjohn et al. 1995) and annual index analysis may be more appropriate if this is the goal. Alternatively, a smaller smoothing parameter may be used i f sudden changes in trend are suspected, although this would be highly subjective. A better solution would be to allow the smoothing parameter to be dependent on the year-to-year variability in the data (data-driven bandwidth selection, L i and Heckman 1996). Unfortunately, this is not straightforward for time series data because positive correlations between successive counts tends to produce over-smoothing (Hart and Wehrly 1986, L i and Heckman 1996). The problem of choosing an appropriate smoothing parameter deserves more study. The rank-trends estimates followed the pattern established in Chapter 2 of having higher statistical significance than the other estimates for each species. Titus (1990) suggested that the rank-trends method may not be suitable for species with many zero counts. In this analysis, the trend estimates did not appear to be any more different for the low-abundance species than the high-abundance species, although the proportion of routes with estimated declines was the most divergent from the other methods for Cerulean Warbler (Table 4.3). In general, the rank-trend method is quite unsatisfactory because it provides little information about population trends, and no way of testing for biologically interesting patterns of population change over space or time. Loglinear modeling The loglinear modeling approach appears to have much potential. It is flexible, allowing data-based selection of the appropriate scale for parameters in space and time. It appears to provide plausible estimates of population change and variance of estimates, except when the model is clearly inadequate (e.g., Chestnut-collared Longspur, see below). In addition, it allows for model checking through formal and informal analyses of residuals. -81 - There are, however, a number of practical problems in implementing the method. Firstly, it is computationally very expensive, especially when the number of sites or years is large. Previous attempts to apply linear modeling to BBS data have had problems with the large amount of memory needed to invert the design matrix (Sauer and Geissler 1990). Program T R I M bypasses these problems by taking advantage of the block-diagonal nature of the design matrix (Pannekoek and van Strien 1996), and so allows the analysis of data from a very large number of sites. Nevertheless, the program is still computationally intensive, especially for models that incorporate serial correlation (these use a different fitting algorithm). For example, fitting the model 5 (region) "with observers" to Belted Kingfisher data (1618 parameters) took 26 minutes on a 90 Mhz Pentium PC without inclusion of serial correlation, but 286 minutes when serial correlation was included (see Appendix 1 for computer, operating system and compiler specifications). A second problem is the failure of the estimates to converge for many of the more complex models, especially those that use the generalized estimating equations estimator (i.e., the serial correlation models). Two types of failure can be distinguished: (1) where the difference in deviance between successive iterations becomes small but some parameter estimates become increasingly large, and (2) where the deviance actually increases in some iterations, and there is no sign that it is reaching a minimum. In the first case the likelihood function for the troublesome parameters is very flat or is asymptotic. The estimated deviance value can be used in model selection, but the parameter estimates will be very inaccurate. The second type of failure is an indication of a complex likelihood function; in this case the deviance from the last iteration may not be close to the maximum-likelihood value so cannot be used. Rescaling may help in both cases, together with an increase in the accuracy of computations (to decrease problems with rounding errors and calculation of the gradient), and the use of alternative fitting algorithms (Bard 1974). It will, however, always be difficult to reliably fit large-dimension or poorly specified models. For three of the four species, the final model chosen did not adequately fit the data (Pearson x /d.f. from the final model was very close to the expected value of 1.0 for Cerulean 2 -82- Warbler, but was 2.41 for Carolina Wren, 8.03 for Chestnut-collared Longspur and 1.24 for Belted Kingfisher). There are two likely explanations for this: either the counts are not distributed as a Poisson random variable or there are some important covariables missing from the model (model incorrectly specified). Biological count data can show greater variability than allowed by the Poisson distribution for numerous ecological reasons (Eberhardt 1978; White and Bennetts 1995), such as clumping of birds in space due to habitat heterogeneity or social behaviors. In addition, variation in the probability of detection of individual birds both within and between surveys creates extra variability in the counts. Nevertheless, model fits as poor as those observed for the Chestnut-collared Longspur are unlikely to be due to overdispersion alone (Anderson et al. 1994), and incorrect model specification must be suspected. I do not pursue this possibility further here, although the BBS dataset contains some additional covariables (such as weather during the survey, time of day and date) that could prove useful. Instead, I pursue the issue of overdispersion, which is almost certainly responsible for at least some of the lack of fit. A common solution to dealing with extra-Poisson variability, and the one implemented by program T R I M , is to estimate the overdispersion parameter, a, based on the Pearson % 2 divided by the degrees of freedom, and to use this to inflate the variance estimates (McCullagh and Nelder 1989). However, this assumes that the variance is a linear function of the mean (variance is now assumed to be c ju, rather than n for a Poisson distribution). This assumption 2 can be checked with a plot of the estimated variance of each count divided by the predicted count against predicted count (Cameron and Trivedi 1990; Winkelmann 1994). Estimated variance is calculated as (observed count - predicted count) . Using the most parsimonious model for each 2 species, these plots show a great deal of scatter with little apparent pattern (Figure 4.14). However, linear regression slopes (Cameron and Trivedi 1990) are significantly greater than zero for all species except Cerulean Warbler (Cerulean Warbler 0.013 + 0.024 (slope + SE); Carolina Wren 0.051 + 0.004; Chestnut-collared Longspur 0.031 + 0.010; Belted Kingfisher 0.232 + 0.057), indicating the possibility of a higher-order relationship between the variance and mean. A n alternative distribution often suggested for overdispersed count data is the negative binomial (Winkelmann 1994; White and Bennetts 1996). In the negative binomial distribution, -83 - the variance is a quadratic function of the mean, and is specified by a dispersion parameter, k, such that variance = mean + k*(mean ) 2 In a preliminary investigation, I compared the fit of the overdispersed Poisson distribution with the negative binomial for the stringently selected Chestnut-collared Longspur data. I chose the Chestnut-collared Longspur because the sample size is relatively small and it is feasible to fit the models using ordinary statistical software. I used SAS PROC G E N M O D (SAS Institute Inc. 1993a) to repeat the Poisson regression from program T R I M and a SAS macro by Hilbe (1994) that used PROC G E N M O D to fit the negative binomial distribution with a log link-function. For model 4 (routes) "with observers", the negative binomial regression (deviance of fit 497.99 on 76 parameters, AIC 649.99) produced very similar parameter estimates to the Poisson regression, and the regression line of estimated variance/predicted count against count was almost identical (slope 0.032 + 0.010). Similar results were obtained when the other T R I M models were used. For these data the negative binomial distribution offered no apparent advantage over the overdispersed Poisson, although this question deserves further investigation. A n alternative distribution worthy of study is the Poisson-logistic (Winkelmann 1994), in which a binomial process, such as the detection of individual birds, is superimposed on a Poisson process, such as the occurrence of birds within the counting radius of each survey stop. This line of reasoning could be further extended to produce other biologically-motivated mixture distributions. The residual analysis presented here could be extended to include numerous other formal and informal checks for specification (McCullagh and Nelder 1989; Venables and Ripley 1994). One obvious improvement to program T R I M in this respect would be reporting of the hat (influence) matrix. Ecological time series are generally assumed to be serially correlated (Jassby and Powell 1990). However, including serial correlation in the models produced only a small improvement in fit for three of the four species, and resulted in a worse fit for the fourth (Carolina Wren). The estimated magnitude of serial correlation was small for all species except Carolina Wren, and so had little effect on the estimated variances or the statistical inferences drawn about the parameter -84- estimates. Little work has been done to date on serial correlation in annual songbird counts, so it is hard to generalize these findings. Fewster and Buckland (1996) incorporated serial correlation into their additive models of British Common Bird Census data using an autocovariate (the observed count in the previous year), and this had a reasonably strong effect on the estimates. Model selection using AIC and QAIC does not select the model that is the most faithful representation of the true data generating process (the "true" model); rather it attempts to select the most parsimonious model, i.e., the model with the optimum balance between bias (too few parameters) and imprecision (too many parameters, Burnham and Anderson 1992). One would therefore expect that, other factors being equal, higher sample sizes (number of routes) would lead to more complex (higher parameter number) models being selected. This trend was evident when comparing the selected models for the stringently selected (low sample size) and tolerantly selected (high sample size) datasets within species. It was not evident when comparing the selected models between high and low sample size species: for example the final model for the largest sample size species (Belted Kingfisher) contained only one set of geographic parameters. Thus, the amount of temporal and spatial variation in population change depends on the species being considered. This highlights the advantage of an estimator that has the flexibility to choose the appropriate scale based on the patterns evident in the data. For example, modeling of loglinear trends on a route-by-route basis (as in linear R R estimators) was the selected model in only one of the three species. The models of temporal change used here could be extended in four ways. Firstly, more complex models that combine trends and annual indices could be used to try and improve the fit of poorly modeled species such as the Chestnut-collared Longspur. For example, the linear model of Sauer and Geissler (1990), which included route-specific linear-multiplicative trends and continental indices, could be extended to the loglinear context. Secondly, models that are intermediate between annual indices (models 3 and 5 and linear trends (models 2 and 4) could be used to try and improve the parsimony of the models for species that are well modeled (e.g., Cerulean Warbler). This could be done by (1) fitting different loglinear trends in different time -85- periods ("broken stick" regression), e.g., in the Carolina Wren dataset fitting one trend for 19671976 and one for 1977 to 1995; (2) adding polynomial terms to the model; and (3) using Generalized Additive Models to estimate smoothed trends. This last option, was implemented on British Common Bird Census data by Fewster and Buckland (1996; see Chapter 5). Incorporating observer differences produced a more parsimonious fit for all species except the Belted Kingfisher. As with the loess smoothing R R estimators, this result is different from Sauer et al. (1994b), who found significant observer effects for all species except the Cerulean Warbler, although their analysis procedure was quite different. Interestingly, the effect of adding observer covariables on the overall trend estimate (Model 2) was qualitatively similar to the smoothing R R methods (Table 4.13). Differences among observers on a route could be caused by variation in observer ability (i.e., in the detection function), or by differences in the way the survey is performed (e.g., in the exact place used as a survey location), or by a combination of both factors. If observer ability is the more important factor, then a more parsimonious model could be fit by using one covariable for each observer by name, rather than treating the same observer on different routes as a new observer as I did in this analysis. A second extension of the observer effects modeling would be to incorporate "start-up" effects, whereby each observer tends to count fewer birds on the first year they survey a route (Kendall et al. 1996). This could be modeled using a binary covariable to indicate the first year by an observer, and tested using AIC or other measures of model fit. Moving-window Mountford's method The Mountford method is based on large-sample theory (Mountford 1982), so it is not surprising that it performed poorly on species observed on few routes. Peterjohn et al. (1995, discussing Mountford's method without moving windows) mention that "there is evidence that technical and practical aspects of model fitting [...] invalidate its use at smaller scales." It was, however, disappointing to find significant problems associated with the method at larger sample sizes as well (for Carolina Wren). Based on this limited sample of four species, it appears that larger window sizes lead to larger estimates of population change and smaller CIs. Larger window sizes give a greater overlap between successive windows, which could result in more -86- accurate "chaining" of successive models. Peach et al. (1994) suggested choosing the largest window size for which all windows that passed the Mountford goodness-of-fit test, but this strategy can apparently lead to implausible estimates of population change in some cases (e.g., for Carolina Wren, Fig. 4.8) but not others (e.g., for Belted Kingfisher, Fig. 4.11). Mountford's method was introduced because it is more efficient than the chaining method it replaced, and because it is less subject to the "random walk" that plagues chaining estimates at low sample sizes (Geissler and Noon 1981, Moss 1985). Because species with low sample sizes in this analysis could only be fit using short window sizes, we may expect to see evidence of random walk in their estimates. This could be tested by fitting an auto-regressive model to the indices and testing for any residual trend, but such a test would have low power with 30 time points. Informal inspection of the plotted indices (Figs 4.3c and d, and 4.9c and d) shows that apart from dramatic fluctuations in one poorly fitting model (1970 and 1971 estimates for Chestnut-collared Longspur) the indices tend to be clustered around the no-change line, and that this effect also occurs for large sample size species (e.g., Fig. 4.8a). I suspect that, rather than random walk, some other failures of the model assumptions are occurring. One possible reason for model failure would be heterogeneous population change between sites, although this is expected to be more prevalent with large sample and window sizes. Plots of counts in successive years provide an informal test of this assumption: the principal axis of the plot of counts in one year against the counts in another should be a straight line running through the origin (Mountford 1982). A n indirect piece of evidence for the importance of between-site heterogeneity is that the only species that produced feasible parameter estimates using Mountford's method was the Belted Kingfisher. Despite the large sample size for this species, the most parsimonious model using T R I M was model 3 (no regional variation), indicating that the population changes were likely very similar throughout the species range. By contrast, the other large sample size species, Carolina Wren, was poorly modeled using Mountford's method. This species was best fit with model 5 (strata) in TRIM, indicating significant geographic variation in patterns of population change. Because of its lack of -87- robustness, it appears that the Mountford method is not useful for estimating population change from B B S data. Generality of findings In this chapter I have estimated population change for four species, chosen from the B B S dataset to represent extremes of abundance and sample size. By using a small number of species I have been able to present and discuss the results in some detail. Nevertheless, this means that my findings are based on a small subset of the entire B B S data. Some conclusions, therefore, must be treated as preliminary. One reliable conclusion is that all of the estimators require further development to solve the obvious problems encountered (Table 4.12). However, based on these case studies alone I cannot reliably conclude when or how often these problems will occur in the whole dataset, or in data from other similar surveys. Further case studies, or less detailed comparisons using the entire B B S dataset, would therefore be useful. In the short term, however, I believe that a focus on methodological development to solve the obvious problems highlighted by these case studies would be more rewarding. S u m m a r y I analyzed continental BBS data for four case study species using eight analysis methods, and used the results to make inferences about the reliability of the methods. The species were Cerulean Warbler (low abundance, low sample size), Carolina Wren (high abundance, high sample size), Chestnut-collared Longspur (high abundance, low sample size), and Belted Kingfisher (low abundance, high sample size). I used three linear route-regression (RR) methods (USNBS and CWS R R and the estimating equations method), two nonlinear R R methods (NNRR, N S R R and rank-trends analysis), loglinear modeling and the moving-windows Mountford Method (see Tables 3.1 and 3.2 for summaries of the methods). I also compared the results under both stringent and tolerant data selection criteria. The linear route-regression methods produced reasonable estimates in all cases except for one, in which a single route-specific weighting caused poor estimates by the USNBS methods. -88- Weighting schemes are seemingly the bugbear of route regression methods and will require further research before an appropriate and robust scheme is found. The estimating equations estimator has, in theory at least, solved the problem of estimating loglinear trends at the route level. Annual indices calculated using the residual method appeared to be biased toward the linear trend estimates in some cases so are probably not a useful way of diagnosing nonlinearities in trend. The loess smoothing methods produced some useful summaries of the data, but also encountered a number of problems: (1) poor performance because of missing data at the beginning and end of the time series; (2) numerical problems with the N S R R estimator; (3) oversmoothing of large, sudden, biologically interesting population changes. Development is needed in all three of these areas. As in Chapter 2, the rank-trends method gave more statistically significant results than the other methods. Although there were no obvious problems with the method, it is hard to recommend because it summarizes the patterns of population change into one statistic which has no biological meaning. The loglinear modeling approach appeared promising, although there were numerical problems in program T R I M with higher-dimension models. A l l species except Belted Kingfisher showed significant lack of fit, indicating that the models were inadequate and/or that the counts are generally overdispersed relative to the Poisson distribution. This was particularly true for the Chestnut-collared Longspur, which showed the greatest variation in counts. The models used here could be extended to incorporate more survey- and observer-specific covariables and to give smoothed estimates of trend over time and/or space. Alternative overdispersed count distributions could also be modeled explicitly, rather than through quasilikelihood. The Mountford moving-window method performed poorly with the low sample-size species, and there were some questions regarding the appropriate window size to choose with one of the high-sample size species (Chestnut-collared Longspur). Belted Kingfisher was the only species that was well fit, and the estimates appeared similar to those from one of the T R I M models. Based on these results, this method cannot be recommended for use with BBS data. -89- oo do CD CD cof >2 ys. .2 <N CO 13 r- CD X CD »—' CD H-» rs vi o 3O Cd T—1 . „ o 1H C C D OI-I T—< *-H N /-v ON 13 CD the rr CN C N , ON o « cd oo d P ta ° -o T3 *5 o -22 H-H P H O •*-> CO cd CD (966 no. td £ cd d cd co cn o o >•» 00 1-1 .S CD H-» O d o CD- o ±j cd OO H—» +3 N 4^T" CO rote" TO O 43 a .2 CD C D C cD o Q . Id +H co C D co +do-» d ON C+H Pi cd -4-» O Q CD C CO D CD ro O dcr d d '3 .§ -+-» GO d O S 3 43 44 >-i 3 c3 CD w 3 ^ PH -90- d— bO O T3 43 cd co VO~ CD o cd <*O CD 43 2 GO u O CD CO CD 2 CD a ^3 3 CN ON ON d '> o c2 T3 CD PH a t-H d ON 1^ u ON ON ^H & § d_P, <u 2 S d H CO ^ S P d 43 CO . CD GO CO <D > CO cd PQ OO CD ^ 4-» .TH d d dO co co d 13 44 1 co <—i OO ^ V cd <+H ON ON D ro C O +-» 44 CD co CD co * PH H- PH O ON ON T3 Cd PH 42 -S ° g ON ON ^H > 3 43. cd vo d 00 ^ -d in 3 3 T3 CD CO g £ 8 u o x C d CD" « d --H co o d O o d d 'aCD 3 O ON ON PH CD u "K d d >^ 43 O o ON CD* PH l^-t s ° g>CD CDp+2 CD C CD D O CO l-l r-~ C C D CD vo vo £-3 cd CO CD cJ3 d 3 "in •B O o •n .2 c£ vi o d 4 3 {0 43° GO , 13 o § 4 3 cn d o C D CO — t "3 d CD C N a o % O C C D cD o O '3 fD T3 _do d 'r^ 00 c-- CD l S sa .—1 1 cd vo ON w VO co o 3 2a 2 CO •d * CD CD CD CD •a 3^ P CT 4 3 CD CD 1 X cd i g o G O 00 o >x. £ S G S o3 GG £ 8 00 G •G o ON cn tx cn 00 OO TT o r-H CN OO 0O ON vo vd cn CN in un © © xt 00 vo <—,• vo ocn NI i n CN 1—1 O r-x CN CN ON o i n oo i n •xf ON r—I I-H i n CN CN ON cn in in cn CN 00 oo rx_ oo o 1—1 CN CN •xT i—; cn r—i CN cn cn cn ON oo i n o oo ON vd © cn T T CN cn o cn VO •xt CN xt in oo c-~ r—i ON cn r-H © ON cn cn cn cn O VO CN CN CN cn cn vo cn cn C N C N C N C N C N C N r~H i n 00 CN VO O ON ON cn Ti- •xf en xj- cn ON o o r-x vo ^ ,-H © CN CN VO Ti- ON cn en cn oo "xt- oo cj oo x, .S3 £ r-x m vo i n cn CN G in in vi 222 Tt- © 1—1 in *—i © *— 1 x—> o Os VO o 1 o o TJ© vo £ 00 n ° ^3 td CN H-H O 00 6 m ^ in in xjvo 00 G O ^ CN encnen '5b o CJ 6 I G cn o 1 ON CN VO o o CN 1—1 T T r- 00 —' CN oo fx VO o CN "xt" in ON i-i oo in m CN r-x cn oo cn m oo ocn CN 00 i—i i—i cn 00 <+H tO I1 O cj in 00 CJ G ON CJ PH 00 in ON in ON VO VO VO VO VO ON - oo <-> 'o OO CN ON VO vo .2 cd 8 « 0C0J ^ ON ON a 11 «a .8 .a 85 m oo cd »H cj JJ in ON 5 cn in in ON VO VO ON r-l ON <U G S <H° « .9 <N 2 2 in in ON VO t-^ VO VO VO VO vo ON ON ON m ON ON ON ON "3 a> £ a r.a .9 is ^o aa oo H a 1 6 p oo H . oo ^ I00 CJ 0000 PH CJ rG u r§ *o o G O •J -91 - ON CN m in ON cn vo .-H vo cn o ^ CN in ON in ON oo OO VO vo VO VO ON G- ON CJ tn 00 G _G i i 'C *o •*-» H oo tH T3 CJ oo CJ « PQ 00 G 3 ON cd a Table 4.3. Estimates of overall population trend for the Cerulean Warbler between 1966 and 1995 (data selected using stringent data selection criteria). Estimation method Trend (% per year) SE 95% CI z P -2.76 -4.85 -3.80 0.48 1.51 _i -1.80 to-3.72 -1.83 to-7.87 -1.56 to-5.98 -5.88 -3.21 -3.38 O.001 0.001 0.001 _i _i +1.01 to-5.90 +0.67 to -4.86 - - -1.52 -1.72 -5.97 0.129 0.085 O.001 -4.06 -2.41 0.48 0.43 -3.10 to -5.02 -1.55 to -3.27 -8.46 -5.60 <0.001 <0.001 -0.13 -0.05 _i +0.40 to -0.30 +0.34 to -0.39 -0.01 -0.03 0.992 0.976 Linear route regression USNBS RR U S N B S est. eqns. R R CWSRR Nonlinear route regression NSRR NNRR Rank-trends -3.04 -2.58 2 2 - T R I M Model 2, incorporating serial correlation with observers without observers 3 Mountford window size 4 with observers without observers 1 2 3 _i Standard errors not calculated on the multiplicative scale and so not comparable. Population change in N N R R and N S R R is calculated as the weighted mean over all routes of the mean estimated value in the years 1970, 71 and 72 minus the mean estimated value in years 1989, 90 and 91 (see appendix 1). Therefore the estimates cover only the middle 22 years of the survey. Standard error (SE) and confidence intervals (CIs) calculated after adjusting for overdispersion (Table 4.4). -92- Ttin £ .s :5 a u SH CJ in CO bo /-s a ID ON • n CN in Td o CJ 53 Td cj ^— I* "I a bo to cj gO rC rJ u N U \° in ON CJ CJ Td PH CJ U cCJ m oli I? cd3 VO oo U ON TT in TT in o © o o o o o o o r-H 1—1 CN CN o VO VO VO VO © vq © O © o © r-H in © o ON in © CJ 3 53 o o vo in in © o in rin 00 © © © in in © *—i cn cn 00 VO 00 in in cn in oo © o o o •*-> o -*-> o cn cn TT cn -4-» 00 TT -t-» ON ON r-- oo © o © o © © © © © CN in © VO CN CN in VO © O © © o cn ON r- in © o cn © O © © o-> o o o oo -* ON VO CN -4-» -4-» CN -f-» Tf Tf TT -*-> cn cn CN ON CN VO CN O © © © vo cn Tien m cn cn © © © © VO © © CN © O -r-» © cn O © VO cn in CN © © CN o CJ Td U 0s in ON Cj-iTO = + co N 9? cn in w CJ -*-» "CJ 'co 'co ON CN in in cn m cn vq cn vq TT vo vq m © © © © © © o o o o o o o ON OO ON OO CN o cn © o © © © © oo VO 00 00 © © ON OO ON CJ <L) o CJ Td vo r- oo VO o © © in ©' Td O rC .5? t> -*-» CJ U es a a CO o Td G e CJ lH on C/5 w 93- c3 Table 4.5. Comparison of loglinear models for the analysis of population trends in Cerulean Warbler using B B S data from 1966-1995 (data selected stringent data selection criteria; sample sizes given in Table 4.2). Calculations were performed using program T R I M . Models are ordered according to the number of parameters estimated; for models 4 and 5 the covariable is given in parentheses. (Because this species only occurred in 1 B B S region, Models 4 (region) and 5 (region) are equivalent to models 2 and 3.) Definitions of the models are given in the text. The selected models (lowest AIC) are shown in boldface. Model No. of params. Deviance p for Pearson % Serial correlation - - - 0.98 1.02 1.02 1.02 1.10 0.627 0.347 0.330 0.364 0.016 0.01 0.05 0.07 0.07 0.11 Pearson X /d.f. 2 2 AIC 1 With observer effects, with serial correlation 5 (stratum) 4 (route) 3 4 (stratum) 2 1 2 228 153 141 116 113 112 776.48 857.61 868.38 929.99 930.12 996.66 1232.48 1163.61 1150.38 1161.99 1156.12 1220.66 With observer effects, without serial correlation 5 (stratum) 4 (route) 3 4 (stratum) 2 1 228 153 141 116 113 112 776.79 857.75 869.60 932.29 932.43 1002.62 0.99 0.98 1.02 1.02 1.01 1.10 0.604 0.627 0.354 0.335 0.372 0.015 1232.79 1163.75 1151.60 1164.29 1158.43 1226.62 Without observers effects, with serial correlation 5 (stratum) 3 4 (stratum) 2 1 2 157 70 45 42 41 980.20 1107.59 1168.40 1172.17 1224.62 - - - 1.30 1.30 1.30 1.32 <0.001 <0.001 <0.001 <0.001 0.15 0.16 0.17 0.20 1294.20 1247.59 1258.40 1256.17 1306.62 - 1294.20 1144.23 1248.52 1265.15 1263.87 1316.97 Without observers effects, without serial correlation 5 (stratum) 4 (route) 3 4 (stratum) 2 1 2 1 2 3 157 82 70 45 42 41 980.20 980.23 1108.52 1175.15 1179.87 1234.97 - - 1.09 1.30 1.31 1.31 1.33 0.019 O.001 <0.001 O.001 O.001 • • Akaike information criterion (AIC) = deviance + 2*no. free parameters. Models came close to convergence, but failed to find accurate parameter estimates. Model 4 (route) failed to converge. -94- Table 4.6. Goodness-of-fit tests for Mountford moving windows applied to BBS data for Cerulean Warbler 1966-1995 (n=41 routes, 112 observers; stringent data selection criteria). X values > x are shown in boldface. 2 2 crit Years With observers No . "routes" (i.e. observers) in model (max 112) Window size 4; d.f. = 3 ; 2 X crit 66-71 68-73 70-75 72-77 74-79 76-81 78-83 80-85 82-87 84-89 86-91 88-93 90-95 No. routes in model (max 41) x 5.19 1.31 0.10 1.54 4.83 3.07 4.12 2.99 2.89 3.26 1.52 1.67 1.78 0.53 34 38 41 41 41 41 40 40 40 40 39 39 39 36 -0.59 2.66 2.02 0.67 1.78 2.64 4.91 3.60 2.24 3.00 1.76 2.19 2.30 0.60 8.72 3.33 3.68 5.43 5.31 10.72 0.16 8.11 9.81 7.24 4.10 4.73 4.67 39 41 41 41 41 41 41 41 40 41 40 39 40 5.99 7.15 2.48 1.84 8.34 -18.13 3.79 9.49 2.54 5.31 2.74 0.78 2.59 x 2 2 = 7.82 (a = 0.05) 66-69 68-71 70-73 72-75 74-77 76-79 78-81 80-83 82-85 84-87 86-89 88-91 90-93 92-95 Window size 6; d.f. = 10; x Without observers 33 38 41 41 43 42 41 38 41 41 44 40 38 36 2 crit = 18-30 (a = 0.05) 41 45 44 45 46 46 45 43 43 52 48 42 42 -95 - Table 4.7. Estimates of overall population trend for the Carolina Wren between 1967 and 1995 (data selected using stringent data selection criteria). Estimation method Trend (% per year) SE 95% CI z P +1.04 +0.85 +1.27 0.25 0.25 _i +1.54 to+0.54 +1.34 to+0.35 +1.89 to+0.67 4.14 3.43 4.16 O.001 <0.001 O.001 _i +2.03 to +0.78 +2.09 to+1.17 - 4.80 7.70 19.31 <0.001 O.001 O.001 Linear route regression USNBS R R USNBS est. eqns. R R CWSRR Nonlinear route regression NSRR NNRR Rank-trends +1.38 +1.62 - 2 2 _i - T R I M Model 2, incorporating serial correlation and overdispersion with observers without observers +1.02 +1.95 0.15 0.14 +1.33 to +0.72 +2.08 to+1.68 6.80 13.93 O.001 O.001 +0.37 +0.40 - +0.47 to +0.26 +0.50 to +0.30 5.69 6.19 O.001 <0.001 Mountford window size 6 with observers without observers 1 2 1 i Standard errors not calculated on the multiplicative scale and so not comparable. Population change in N N R R and N S R R is calculated as the weighted mean over all routes of the mean estimated value in the years 1971, 72 and 73 minus the mean estimated value in years 1989, 90 and 91 (see appendix 1). Therefore the estimates cover only the middle 21 years of the survey. -96- Table 4.8. Comparison of loglinear models for the analysis of population trends in the Carolina Wren using BBS data from 1967-1995 (see Table 4.5 legend for further information). The model with the lowest AIC and QAIC is shown in boldface, but this model was not used for parameter estimation because the estimates were unreliable (see and text). A l l models failed the Pearson % goodness-of-fit test (p < 0.001 in all cases). 4 2 Model No. of params. With observer effects, with 4 (route) 1084 4 (stratum) 836 4 (region) 818 2 817 1 816 4 Deviance Pearson X /d.f 2 serial correlation 17516.25 3.14 19737.89 20045.21 3.18 20050.71 3.18 20275.97 3.20 Serial correlation AIC - 19684.25 21409.89 21681.21 21684.71 21907.97 9437.66 9863.69 9955.24 9955.51 10047.01 14029.19 19202.22 17108.89 17444.89 20973.55 21271.38 21275.91 21514.43 7432.30 9237.60 8120.75 8227.44 9682.60 9785.15 9758.86 9883.68 23584.87 26692.21 28021.00 28053.16 29481.45 10415.37 11414.84 11945.26 11957.44 12549.04 17227.06 22902.71 23880.37 23208.00 26992.21 27537.42 27579.64 29267.30 8661.21 10132.26 9871.64 9978.17 11414.84 11744.57 11760.92 12460.17 0.32 0.33 0.33 0.33 4 4 4 QAIC 2 3 With observer effects, without serial correlation 5 (stratum) 1376 11277.19 1084 17034.22 4 (route) 2.41 5 (region) 872 15364.78 844 15756.89 2.47 3 3.12 19301.55 4 (stratum) 836 19635.38 3.17 4 (region) 818 817 19641.91 3.17 2 19882.43 3.20 1 816 Without observers 4 (route) 4 (stratum) 4 (region) 2 1 1 effects, with serial correlation 536 22512.87 26116.21 288 4.08 270 27481.00 27515.16 4.09 269 4.32 28945.45 268 3 0.46 0.46 0.48 Without observer effects, without serial correlation 5 (stratum) 4 (route) 5 (region) 3 4 (stratum) 4 (region) 2 1 4 1 2 3 4 4 828 536 324 296 288 270 269 268 16879.06 21830.71 22224.37 22616.00 26116.21 26997.42 27041.64 28731.30 3.23 3.30 3.96 4.06 4.07 4.32 Akaike information criterion (AIC) = deviance + 2*no. free parameters. Quasi-likelihood AIC (QAIC) = deviance/dispersion + 2*no. free parameters. The dispersion value used was 2.41, from model 5 (region) with observers, the best fitting model for which an estimate of dispersion could be obtained. Models 5 (stratum), 5 (region) and 3 failed to converge. Models came close to convergence, but failed to find accurate parameter estimates. -97- Table 4.9. Estimates of overall population trend for the Chestnut-collared Longspur between 1967 and 1995 (data selected using stringent data selection criteria). Estimation method Trend (% per year) SE 95% CI z P -0.41 -0.48 -0.33 0.62 0.58 _i +0.83 to-1.65 +0.68 to-1.64 +1.02 to-1.82 -0.66 -0.83 -0.45 0.51 0.41 0.65 _i _i +2.28 to -3.44 -0.32 to -4.93 - - -0.71 -2.20 -5.17 0.48 0.03 <0.001 +0.99 to -2.65 -0.04 to -4.33 -0.92 -1.99 0.356 0.046 +2.58 to-1.64 +1.01 to-1.31 -0.03 +0.32 0.976 0.746 Linear route regression USNBS R R U S N B S est. eqns. R R CWSRR Nonlinear route regression NSRR NNRR Rank-trends -1.05 -2.94 2 2 - T R I M Model 2, incorporating serial correlation and overdispersion with observers without observers -0.84 -2.17 0.91 1.09 -1.56 +4.01 _i Mountford window size 4 with observers without observers 1 2 _i Standard errors not calculated on the multiplicative scale and so not comparable. Population change in N N R R and N S R R is calculated as the weighted mean over all routes of the mean estimated value in the years 1971, 72 and 73 minus the mean estimated value in years 1989, 90 and 91 (see appendix 1). Therefore the estimates cover only the middle 21 years of the survey. -98- Table 4.10. Estimates of overall population trend for the Belted Kingfisher between 1968 and 1995 (data selected using stringent data selection criteria). Estimation method Trend (% per year) SE 95% CI z P -0.70 -1.47 -0.81 0.13 0.41 _i -0.43 to -0.97 -0.66 to -2.28 -0.35 to-1.26 -5.21 -3.63 -3.50 O.001 <0.001 <0.001 -2.10 -1..41 2 _i 2 _l -0.39 to-3.53 -0.45 to -2.27 - - - -2.42 -2.88 -3.32 0.015 0.004 0.001 Linear route regression USNBS R R USNBS est. eqns. R R CWSRR Nonlinear route regression NSRR NNRR Rank-trends T R I M Model 2, incorporating serial correlation and overdispersion with observers without observers -1.08 -0.75 0.20 0.15 -0.66 to-1.49 -0.45 to-1.04 -5.14 -5.00 <0.001 O.001 -0.52 -0.42 _i +0.14 to-1.11 +0.15 to-1.03 -1.27 -1.13 0.204 0.258 Mountford window size 4 with observers without observers 1 2 - 1 Standard errors not calculated on the multiplicative scale and so not comparable. Population change in N N R R and N S R R is calculated as the weighted mean over all routes of the mean estimated value in the years 1972, 73 and 74 minus the mean estimated value in years 1989, 90 and 91 (see appendix 1). Therefore the estimates cover only the middle 20 years of the survey. -99- Table 4.11. Comparison of loglinear models for the analysis of population trends in the Belted Kingfisher using B B S data from 1968-1995 (see Table 4.4 legend for further information). A l l models failed the Pearson % goodness-of-fit test (p < 0.001 in all cases). 2 Model No. of params. Deviance Pearson X /d-f- Serial correlation 2 AIC 1 QAIC 2 With observer effects, with serial correlation 5 (region) 1618 11553.46 1.11 4 (stratum) 1570 11687.10 1.11 3 1564 11646.63 1.11 4 (region) 1540 11719.67 1.12 2 1538 11732.80 1.12 1 1537 11760.91 1.12 -0.05 -0.04 -0.04 -0.04 -0.04 -0.04 14789.46 14827.10 14774.63 14799.67 14808.80 14834.91 13673.01 13697.74 13649.18 13667.16 13675.02 13698.42 With observer effects, without serial correlation 5 (region) 1618 11555.99 1.11 4 (stratum) 1570 11691.19 1.11 3 1564 11648.46 1.11 4 (region) 1540 11723.62 1.12 2 1538 11736.63 1.11 1 1537 11764.47 1.12 - 14791.99 14831.19 14776.46 14803.62 14812.63 14838.47 13675.30 13701.43 13650.83 13670.73 13678.48 13701.63 - 14534.98 14569.86 14527.91 14647.00 14592.04 14616.41 13248.31 13270.54 13231.49 13286.14 13284.40 13306.22 14553.47 14566.69 14527.59 14590.12 14588.61 14612.61 13246.95 13267.68 13231.20 13283.05 13281.30 13302.79 3 4 3 4 Without observers 5 (region) 4 (stratum) 3 4 (region) 2 1 4 effects, with serial correlation 610 13314.98 562 13445.86 1.24 556 1.24 13415.91 532 13529.54 1.25 530 13532.04 1.25 13558.41 529 1.25 3 Without observer effects, without serial correlation 5 (region) 610 13313.47 1.23 4 (stratum) 562 1.24 13442.69 3 556 1.24 13415.59 4 (region) 532 1.24 13526.12 2 530 1.24 13528.61 1 13554.61 1.24 529 4 1 2 3 4 0.02 0.02 0.03 0.03 0.03 3 - - - Akaike information criterion (AIC) = deviance + 2*no. free parameters. Quasi-likelihood AIC (QAIC) = deviance/dispersion + 2*no. free parameters. The dispersion value used was 1.11, from model 5 (region) with observer effects and serial correlation, the best fitting model for which an estimate of dispersion could be obtained. Model 4 (route) and model 5 (stratum) filed to converge Models came close to convergence, but failed to find accurate parameter estimates. -100- C3 CD ,G < D < < < < U < Q < 00 G CD CD cci a CD (D IS OO G << < < O < Q < OH CO 00 T3 o o 1 < D CO u G PH < < < Q < O H 1 II u a T3 a PH CD 00 G 00 << < < «J <C Q PQ" T3 PH G PQ -G CD « "o ea J> CO — X G VH (D JD VH t— OOO O O O Q X 1 G 00 a JG C+H O CO O II G u -d »-7 bO T3 bo CD C -D > C .G CD || co c o < D O ^ T3 "GGCDco<D CD O G O CD CD -d S £ co CD i n 3 ON V O CD -Q •- rG O CD G bO CD •5G « tD 0 CO •*-» CD PQ h-1 D co C JO -g •H ^^ •§ HG 00 G O LG "G 00 2 || CD PH O CD d OH O ro ooo O OOQ ffi oo CD o "2 co +i G G CD — •H CD CD T3 rG -G «s C D CDCO O JCD D CD 1 H << < U U < Q O PQ JG CD O CO 3 Q CD s a u el <D (D 00 PQ PQ < O O < Q PH PQ oo L o CD CD *c •H S o ~o o a "CD CO G CD PQ iH •+JrH **-^ £ is o T3 G ~ CD cD o O C rt < ft .U H -O a o -a oo CD < DG OH .a O co o G .5 c o *H f | , 00 00 Ci o .3 CD > Tc 3 T 3 o G O -d o CD CD a 5G < S3 u a a ^ - 101 - 1 B bo o 3 T3 o g CO II CD° *c3 < -° CD CD CD -d G Table 4.13. Effect of incorporating observer differences on trend estimates for four species from three analyses (one prior analysis that used the linear R R method and two analyses in this chapter that used both "with observer" and "without observer" models). Species Cerulean Warbler Sauer etal. 1994b 1 Smoothing R R 2 Loglinear modeling not significant Decrease decrease Carolina Wren decrease decrease decrease Chestnut-collared Warbler decrease increase increase Belted Kingfisher increase decrease decrease 1 2 3 3 Analysis of the correlation between difference between observers (calculated from linear R R analysis) and time. Pearson correlations were reported for the three species with statistically significant results: r = 0.15, 0.23 and -0.05 respectively. Difference in overall trend estimates between N S R R and N N R R methods, from Tables 4.3, 4.7, 4.9 and 4.10. In all cases, 95% confidence intervals for N S R R and N N R R estimates overlap by a large margin, indicating that the differences are small in relation to the variability in trend estimates. Difference in overall trend estimates between model 2 "with observers" and "without observers", from Tables 4.3, 4.7, 4.9, and 4.10. Introducing observer effects into the models resulted in a lower AIC in all cases except Belted Kingfisher, but in all cases the 95% confidence intervals for the estimates overlap. - 102- Figure L e g e n d s Figure 4.1. Counts of Cerulean Warblers on each BBS route between 1966 and 1995. Only the 41 routes that passed the stringent data selection criteria are shown. The two routes that together contain 45 percent of the total counts are highlighted with dark lines (upper line is route 82031, lower line is route 82032). Figure 4.2. Relationship between measures of variation in counts on a BBS route and the mean count on that route for four species. If counts on route follow a log-normal distribution then the standard deviation would be proportional to the mean and the plots of standard deviation (stdev)Anean (i.e., coefficient of variation) against mean would show no trend. If counts were Poisson distributed then the variance would be equal to the mean and the plots of variance/mean against mean would show no trend, and have an average value of 1.0. Values higher than 1.0 in the variance/mean-mean plots indicate overdispersion relative to the Poisson distribution. Figure 4.3. Population indices for Cerulean Warblers between 1966 and 1995, calculated using eight different estimators. BBS data were screened using stringent selection criteria (see text). As with all subsequent figures, the population indices have been scaled so that their value in the first year is 1.0, and are plotted on a logarithmic axis with log-base 2. Figure 4.4. Population indices for the Cerulean Warbler between 1966 and 1995 calculated using model 3 of program TRIM. The dashed lines are 95% parametric confidence limits for the comparison between the index estimate and the index estimate in first year. In figures c) and d) the confidence limits have been scaled to account for the estimated overdispersion in these models (Table 4.4). Figure 4.5. Mountford moving-window estimates of population indices for Cerulean Warblers between 1966 and 1996. The dashed lines are 95% bootstrap consistency intervals. Estimates outside the range of the consistency intervals indicate failure of assumptions of the Mountford model. Figure 4.6. Population indices for the Cerulean Warbler between 1966 and 1995, calculated using eight different estimators, with data selected using tolerant data selection criteria. Figure 4.7. Population indices for the Carolina Wren between 1967 and 1995, calculated using eight different estimators, with data selected using stringent data selection criteria. Figure 4.8. Mountford moving-window estimates of population indices for Carolina Wren between 1966 and 1996 for "without observers" models at 4 window sizes. The dashed lines are 95% bootstrap confidence intervals. Figure 4.9. Population indices for the Chestnut-collared Longspur between 1967 and 1995, calculated using eight different estimators, with data selected using stringent data selection criteria. - 103 - Figure 4.10. Population indices for the Chestnut-collared Longspur between 1967 and 1995, calculated using seven different estimators (program T R I M failed to produce log-linear model estimates), with data selected using tolerant data selection criteria. Figure 4.11. Population indices for the Belted Kingfisher between 1968 and 1995, calculated using eight different estimators, with data selected using stringent data selection criteria. Figure 4.12. Population indices for the Belted Kingfisher between 1968 and 1995, calculated using eight different estimators, with data selected using tolerant data selection criteria. Figure 4.13. Relationship between mean abundance and estimated route trend for the linear route regression estimators using tolerantly selected data for Belted Kingfisher between 1968 and 1995 (n = 1603 routes). For both the USNBS and CWS methods, trend estimates for lowest abundance routes are biased towards zero. Figure 4.14. Residual analysis from loglinear models. These graphs are comparable with those on the right side of Fig. 4.2 (variance/mean against mean), except that these use the residuals and predicted values of each count using the best-fitting loglinear model for each species. If counts were distributed as Poisson random variables the graphs should show no trend and an average variance/predicted value of 1.0. In all cases, a linear regression line shows an increasing trend with increasing predicted value. -104- Figure 4.1 Figure 4.2 a) Cerulean Warbler ro 3- c 3ro (D E r "a? o-l c ro CO 1- | ~i—•—i— —r 1 6 8 - 10 i T 12 14 1 1 r 6 16 b) Carolina wren 5- 12 14 16 201 c ro is- 4c ro CD 3 E leo c ro JD 2 CO 10 Mean Mean e 8 s' H i 0 —i— —i—•—i—'—i—'—i— —r 10 20 30 40 50 60 Mean 1 —i—•—i—•—i—•—i—•—i—•—r 10 20 30 40 Mean 50 60 1 C) Chestnut-collared longspur 3 70c c 6 0 50H i * 8 30 20H > 10-1 0-1 3 ro o. CO o-i 100 200 300 100 200 Mean Mean d) Belted kingfisher 51 51 5 " ro 4 4- c ro CD o. <D <U 2 CD E E 3 3 0 1 CO .i 2 1H 04 1 2 Mean Mean -106- 300 Figure 4.3 a) Mean count b) Linear route regression 2.000- 2.000- i.oooH 1.000 CD CD o o TC3 T3 -§ 0.500 "§ 0.500 C CO C CO C x §5 •a c "D C 0.250 0.250 4 USNBS RR USNBS est eqns. RR -CWS RR — Mean count Area weighted mean count 0.125 I 65 70 I I 0.125 I • 75 80 85 Year 90 |—I—I—I—I—1—I—I—I—l—1—I—l—i—f— 70 75 80 85 Year -i—i—I—i—i—i—r~ 90 95 d) Other methods, without observers 2.000-1 c) Other methods, with observers 2.000- © o c I 65 95 i.oooH CD 1.000 H o c •a c CO •D "§ 0.500 •co 0.500 x x CO c TD CD C 0.250 0.250 — NSRR — TRIM Model 3 Mountford window size 4 0.125 i 65 1111 i 1 1 1 70 •i 1 75 •i • ' i ' i • • i 80 85 90 95 Year 11 11 1 11 - 107- 0.125-U r 65 -NNRR -TRIM Model 3 Mountford window size 4 - r" 90 r_ 70 75 80 85 Year r 95 Figure 4.4 a) with observers, with serial correlation 2.000 i 65 70 75 80 Year 85 90 95 -108- b) with observers, without serial correlation 2.000-1 65 70 75 80 Year 85 90 95 Figure 4.5 -109- Figure 4.6 a) Mean count 2.0001 CD b) Linear route regression 2.000- 1.000 1.000- CD O o c c CO TD CO TD C C 3 ro 0.500 3 0.500 TD C 0.250 0.250 H - USNBS RR USNBS est eqns. RR -CWS RR -Mean count •-Weighted mean count 0.125 65 n—i—| i i—i—i—pi—i—r—i—|—i—i—i i 70 75 80 Year | i—i i i—[ i 85 0.125 i—i—i—j- 90 95 c) Other methods, with observers 2.000 i <D 65 •i• 75 • 1 • 85 80 Year 1 70 1 1 1 ' i 90 1 95 1 1 1 1 d) Other methods, without observers 2.0001 i.oooH CD 1.000 o c CO •o c o c co TD C v\ 3 •§ 0.500 •§ 0.500 X CD x C TD TD /vyv--' V <D C 0.250 A 0.125A.i ' , , ,• i 1 1 65 0.250 NSRR TRIM Model 3 Mountford window size 4 1 1 1 1 70 1 >> 1 ' 1 75 80 85 Year 1 1 1 1 1 1 1 1 •i 90 1 1 1 1 1 95 - 110- 0.125 -NNRR -TRIM Model 3 Mountford window size 4 i 65 i i | i i i 70 i | i 75 i i i | i i i i | 80 Year i i 85 i i | i i 90 i i | 95 Figure 4.7 a) Mean count b) Linear route regression Mean count Area weighted mean count 4.0 4 0 1 CO 1 0 c 0 2.0- c 2.0" X) cd $ _, ~ USNBS RR USNBS est eqns. RR — CWS RR co 1 0 0.5 r 65 1 T 1 - 70 -i—|—r- 75 80 Year 85 90 0.54 T 95 70 i i 75 i—i—pT—i—i—i—|—i—i—i—i—|—i—i—i—i—r 80 Year 85 90 95 d) Other methods, without observers c) Other methods, with observers 4.0- i i i i i — | I-I-I-I | 65 -NSRR TRIM Model 5 (region) Mountford window size 6 4 0 . CO C OD CO2.0TJ c « 2.0H o c NNRR — TRIM Model 5 (region) Mountford window size 6 c XI CO x> CO H— o X 1.0C D TcS 0.565 §5 1 0 0.5 i i 65 70 1111 70 75 80 Year 85 90 95 - Ill- 1 75 80 85 Year 90 T 95 Figure 4.8 - 112- Figure 4.9 c) Other methods, with observers d) Other methods, without observers 4.00- 4.00- 0 NNRR TRIM Model 4 (route) Mountford window size 4 2.00 o c co T> c n 1.00 XI CO X CD C 0.50 — NSRR — TRIM Model 4 (route) Mountford window size 4 0.25-L r 65 ~V r r 70 75 80 Year 85 ~* i 90 0.25 95 -113 - -i 70 1 65 l— 75 80 Year 85 90 95 Figure 4.10 a) Mean count b) Unear route regression 8.001 8.001 Mean count Area weighted mean count —USNBS RR USNBS est. eqns. RR CWS RR 4.00 -\ CD o c -g 2.00 c Z2 n CO H— o x CD TJ C LOO H 0.50 0.254i • • • • i -• r " 75 65 70 ,- r ' i' '''i'''• i 80 85 90 Year 1111 i 95 -114- 0.2565 70 75 80 Year 85 90 95 Figure 4.11 a) Mean count b) Linear route regression 2.00- 2.00- CD O c CO 1.00T3 C 3 Si CO o X 0.50D TCD C CD O c 1.00CO T3 C 3 XI CO o X 0.50CD T3 C -Mean count Area weighted mean count ~ USNBS RR USNBS est.eqns. RR — CWS RR 0.25- 0.2565 75 70 80 85 90 95 • i 1 65 70 75 T 85 80 r 90 T 95 Year Year c) Other methods, with observers d) Other methods, without observers 2.00- 2.00- CD O c 1.00CO T3 C 3 _(0 Q <D o . § 1.00-1 c CO X " ~r- 11 H— 0.50 O X 0.50CD TJ C — NSRR — TRIM Model 3 Mountford window size 14 <D C 0.25 -1—i—i—I—i— 65 70 75 1 i' 111 80 0.25 i' ''i ' i 1 85 11 90 1 95 NNRR TRIM Model 3 Mountford window size 14 • • i• • • i • • i• • i ' • • i • 1 65 1 70 75 1 1 80 Year Year -115- 11 85 1 1 90 11 95 i Figure 4.12 a) Mean count b) Linear route regression 2.00- 2.00- CD O c CO 1.00T3 C Z3 xi CO M— o X 0.50CD •D C C OD .§ 1-00-1 c XI CO XD0.50 C C — Mean count • Area weighted mean count 0.2565 0.25 ]—i—i—i—i—|—i—i—i—i—| i I—i—i—|—i—i—i—r | l—i—i—i—]—r~ 70 75 80 Year 85 90 95 - USNBS RR USNBS est eqns. RR CWS RR -]—i—i—i—i—j—I—i—i—i—]—r—I—i—i—[—r 65 70 75 80 Year 'l 85 11 11 l 90 1 T 95 c) Other methods, with observers d) Other methods, without observers 2.00- 2.00- 0 o c: CO 1.00c XI CO H— o X 0.50CD T3 C 0.2565 o « 1.00 c XI CO XD0.50 C C — NSRR — TRIM Model 3 Mountford window size 14 1 1 i i i 70 75 1111 1 111 i• i i 80 85 90 Year 1 111 111 -NNRR -TRIM Model 3 " Mountford window size 14 0.25 11 1 1 95 -116- ~r 75 _T 65 70 r_ i i 80 85 Year 1 1 11 1 1 T 90 1 95 a) USNBS route regression 8.000TD 4.000c 2 2.000"S 1.000tS E 0.500V-• co 111 0.250 0.125- 2 - 1 0 1 log (mean count) b) CWS route regression 8.000TJ 4.000c 2 2.000 "S 1.000 E 0.500 + * C-O Ul 0.250 H 0.125-3 - 2 - 1 0 1 log (mean count) c) USNBS est. eqns. route regression 8.000 -O 4.000 c 2 2.000 "S 1.000•s E 0.500 to Ul 0.250 0.125 -2 log (mean count) -117- -r 2 Figure 4.14 a) Cerulean Warbler Model 3 with observers b) Carolina Wren Model 5 (region) with observers 130120- CD 110- co > 100H TD CD 90 ^ TJ CD 80 Q. 70H ^O D 60 C CO 50 H '§ TD 40 CD to 30 \— E 20 CO UJ 10 0 5 10 15 20 Predicted value 25 1 1 1 1 c) Chestnut-collared longspur Model 4 (route) with observers d) Belted kingfisher Model 3 without observers 160- 100- -1 1 90- (D •S f0 1 2 0 CD O TJ C D100 i— £0 1 0 10 20 30 40 50 60 70 80 90 Predicted value (D 140CO 8 C .cO T —i——i——i——i—•—i——i—'—r —i——r 80 o 70 60 g 50 c "S 8 0 CO 60 40 CO > TCD 3 30 40 |2 20 to £ 20 "•4—» I 50 100 150 200 Predicted value 250 CO w 10 1 Predicted value -118- T 4 Chapter General 5 C o n c l u s i o n s Introduction In this thesis I have investigated the statistical methods used to estimate long-term population change from extensive wildlife surveys like the North American Breeding Bird Survey (BBS). I have reviewed previous work (Chapter 1), demonstrated the importance of analysis method and data selection (Chapter 2), synthesized the issues involved and shown the potential of simulation studies to resolve some of these issues (Chapter 3), and compared the methods empirically using case studies (Chapter 4). Here, I bring this work together by providing an evaluation of the methods, with recommendations for when they should be used and how they may be further developed. The evaluations are summarized in Tables 5.1 (annual index methods) and 5.2 (trend estimation methods); these are intended to complement the summary descriptions of the methods given in Tables 3.1 and 3.2. I also provide a simplified guide to selecting among the methods (Fig 5.1). In the second part of the chapter, I discuss some future directions for research and suggest three promising new approaches to trend analysis. Lastly, I discuss the implications of my research for the use of extensive wildlife surveys like the BBS in estimating population change. Evaluation of analysis m e t h o d s Annual index methods Calculating annual indices is appropriate where the focus is on the actual yearly population changes, rather than the underlying trend. This is useful when the goal is to detect sudden changes in abundance (e.g., the 1976-77 population crash in Carolina Wren, Fig. 4.7), or when further modeling is to be undertaken (e.g., the effect of winter temperature on breeding population fluctuations in Song Thrush, Baillie 1990). Annual indices are also used forjudging the fit of trend estimates (e.g., Peterjohn et al. 1995; see Collins 1994 for further discussion of - 119- the uses of annual indices). In all of these cases the goal is to produce index estimates that accurately reflect the true population changes. Of the six annual index methods summarized in table 3.1, three (chaining, imputing and linear modeling) can be discounted from consideration on the basis of previous work (Chapter 1; Table 5.1) and were not included in the empirical evaluations in this thesis. This leaves three candidate methods: Mountford's method, the residual method, and loglinear modeling. Mountford method: Even with the moving-window extension, Mountford's method performed poorly in three out of four case study species (Chapter 4). Poor performance was probably caused by violation of the assumption that population change is the same at all sites. This assumption will likely only be met in small datasets or geographically homogeneous species, so the method cannot be recommended for general use. Residual method: This approach was developed to provide informal indices that complement linear regression trend estimates (Sauer and Geissler 1990). As the original authors noted, the indices are biased when the trend lines are invalid, and this was confirmed in at least one of the case studies (Fig 4.10b; see also Figs. 4.3b, 4.6b, 4.1 lb and Discussion, Chapter 4). The method is clearly not useful forjudging the fit of trend estimates, and cannot be generally recommended without further studies to determine the situations under which it works well. Loglinear modeling: Loglinear modeling, with separate effects for each year, can be recommended on both theoretical and empirical grounds. Theoretically, it allows a defensible model of the data to be constructed using the well-established tools of generalized linear modeling ( G L M , McCullagh and Nelder 1989). Empirically, the method worked well in the case studies, although the final model was not satisfactory for one species (see discussion, Chapter 4). One disadvantage of the method is the practical problems associated with implementation: complex models require a large amount of computer time and sometimes fail to converge, even using the specially written computer program T R I M (Chapter 4). Conclusion: The best annual index method appears to be loglinear modeling, although there may be problems implementing the method for large datasets and for species that show high variation in counts over time and space. - 120- Trend analysis methods Determining population trends is the principal goal of many long-term monitoring programs like the BBS. Unlike annual index analyses the focus is no longer on the actual population changes. Instead it is the underlying tendency in abundance that is of interest, and this is defined subjectively (Chapter 3). Thus, the assumed model of trend is a key factor that differentiates the analysis methods (Table 3.1). The methods evaluated in this thesis used three trend models: linear multiplicative, additive and rank (Fig 3.1; polynomial multiplicative models have been suggested but not implemented, and can be viewed as a type of additive model (Hastie and Tibshirani 1990)). The linear model requires estimation of fewest parameters (so is useful when there are many missing data), but is only realistic for most populations for relatively short time periods (perhaps 5-15 years for most species monitored by the BBS). Additive trend models are more realistic, but require more parameters and are also susceptible to incorporating changes in abundance that may not be regarded as trend. In rank models only the direction of trend is estimated, not its magnitude. Rank models are therefore more robust, but less informative. The second major difference between trend analysis methods is in the approach to trend estimation (Table 3.1). Route regression (RR) methods first estimate trends on each site (route) and then combine these to produce regional trend estimates (a form of multilevel or hierarchical modeling, Goldstein 1995). By contrast, loglinear modeling estimates trend directly from the counts at a level that can be determined by the geographic variation in the data. Both approaches are well established in the statistical literature and I have not been able to determine in this thesis whether one is generally superior. Both have practical problems associated with them: route regression methods have problems associated with the imprecision of route trend estimates and the methods for combining route trends to form regional estimates, while loglinear modeling can run into computational problems with large, complex models. I consider these issues in the specific context of each method in the following paragraphs. Linear route regression: This method has a number of practical advantages including low computational requirement and ability to deal with much missing data. Of the two methods -121 - for estimating route trends (Table 3.1), the estimating equations approach is preferred because it is not biased by low counts (Link and Sauer 1994; Chapter 4). However, a number of other decisions must be made when implementing the method (Table 2.1) and it is often not clear which is best. Different implementations usually produced reasonable (although divergent) estimates in the empirical comparisons, except in some cases where poorly estimated weightings (used in combining route trends) lead to obviously unfeasible results (Chapters 2 and 4). The example simulation study (Chapter 3) indicated that the variance estimates may be too small and that one implementation (CWS RR) may be negatively biased, although these results are preliminary. In conclusion, the method may be recommended where there are a priori reasons to wish to estimate linear trends, and it will be particularly useful in large surveys with many missing counts. Smoothing route regression: The primary advantage of the smoothing methods (NN and NSRR) is their ability to estimate nonlinear (additive) trends. This flexibility comes at a price: the methods cannot deal with much missing data. The case studies (Chapter 4) also showed that there can be significant problems with estimates at the beginning and end of the time series, and with fitting the observer-effects model (NSRR). Assuming that these problems can be resolved (see Discussion, Chapter 4), the smoothing R R methods can be recommended for most situations where there are few missing data, or where the analyst is willing to use stringent data selection to exclude sites with few counts. Rank trends route regression: As implemented, the rank trends method uses a naiv approach to combine route trends into regional estimates and does not take account of differences among observers. The empirical comparisons showed that it tends to give more statistically significant results than the other methods (Chapters 2 and 4), and the simulations showed that it may be biased by geographic variation in trend (Chapter 3). The method may be further developed, but given that rank models are intrinsically uninformative (because they only estimate the direction of trend, not its magnitude) this is not a high priority. Loglinear modeling: This method provides an appealing alternative to linear RR. modeling framework is well established and there are fewer arbitrary decisions to be made. The -122- The case studies did not focus on linear trend modeling using this approach because annual indices were found to give a better fit to the data in most cases (using Akaike's Information Criterion, Chapter 4). It would be useful to extend this method to allow fitting of additive trend models (see next section). Practical problems fitting large, complex models mean that this method is best suited to relatively small datasets. Regression of annual indices: This approach is a natural extension of annual index analyses, but relies on the indices being unbiased (e.g., poor performance of linear regression on Mountford indices, Chapter 4). The best annual index method is loglinear modeling, but I did not evaluate the use of loglinear model indices because the method can also be used to produce trend estimates directly. It is unlikely that this approach would offer any improvements over direct trend estimation. Conclusion: Unlike annual index estimation, it is not possible to recommend one single method for trend estimation in all situations. The choice of method depends on the assumed model of trend and the approach to trend estimation (Fig. 5.1). The length of time series, size of dataset and amount of missing data are all important considerations when making these decisions. Depending on the situation, one or more of the following may be appropriate: linear RR, smoothing R R or loglinear modeling. A l l of these methods suffer from practical drawbacks and could benefit from further development. Future research Development of current methods For linear RR, the biggest issue is the choice of weightings. Three weightings are used: area, abundance and precision weights (Table 2.1). Of these, it is the occasional very high abundance weight that causes the few unreliable trend estimates (Chapters 2, 4). This can be avoided by capping (Winsorizing) the abundance weights, by using another method to estimate the weights, or by including the uncertainty in abundance weight estimates into the precision weights (Chapter 4). Precision weighting is also controversial (James et al. 1990). One option worth investigating is the use of hierarchical Bayes methods (Gustafson 1997) to combine route - 123 - trends into regional estimates. This approach would automatically account for the precision of the route trend estimates and could also be used to constrain the estimates to within reasonable limits by choosing a suitable prior distribution for route trend. The issue of weightings is not so problematic for the smoothing R R methods because they do not require abundance weights (James et al. 1990). There are, however, two other improvements that must be made before the methods can be considered credible: (1) increased reliability of estimates at both ends of the time series, perhaps through extrapolation of the route trend estimates, and (2) solution of the numerical problems associated with fitting the N S R R model (Chapter 4). Another area of potential development is in the selection of the amount of smoothing to perform on the counts, which is currently set subjectively but could be set according to the amount of between-year variability in the data (Chapter 4). The loglinear models used in Chapter 4 (both annual index and trend) could be improved in several ways. Alternative count distributions should be investigated, such as the negative binomial and Poisson-logistic (Chapter 4). Models that incorporate site-specific linear trends and overall regional annual indices simultaneously (e.g., Sauer and Geissler 1990; Collins 1994) should be tried. More realistic models for trend should be considered, such as additive models (see below). Lastly, if the method is to be useful for large datasets then there must be improvements in the fitting algorithm (Chapter 4). Further evaluation of methods The empirical and theoretical work in this thesis has resulted in a substantial advance in our understanding of the relative merits of the analysis methods. Nevertheless, a number of important questions remain unanswered: (1) what is the relative accuracy of the R R approach compared with loglinear modeling? (2) which weighting schemes produce the most reliable results in R R methods? (3) how robust are the methods to mis-specification of trend? (4) are the variance estimates of the methods biased? I suggest that these questions can be answered by simulation studies similar to the example in Chapter 3. The model used to simulate data in Chapter 3 could be improved in a number of respects. Firstly it could be made more realistic by incorporating a simple stochastic population model. -124- For example, assume the actual population size in yearly at site s, n , is a Poisson random vs variable: n ys where N ~ Poisson (N ) ys is the expected population size. This is given by a discrete version of the simple y s logistic growth model: Ny = n(y.i) (\ +R(\ - n( -i) IKy )) S s y s S where R is the per capita growth rate and Ky is the carrying capacity (Krebs 1985). R could be S set at various values between 1 and, say, 1.5 to represent different species. (In a more complex model R could depend upon>> and/or s.) Ky could be set to change through time using an S additional function that depends upon the assumed local trend, regional trend, and distribution of random interventions. Trends could be linear-multiplicative or could change over time. The distribution of N()s and K()s among sites could be set by the geographic models proposed in Chapter 3 (Fig. 3.2) or by some other non-geographic model. A similar simulation model was proposed by Greenwood (1991). A second improvement would be to explicitly model observation error by assuming c ys ~ binomial (n ,py ) ys S where Cy is the count in yeary at site s, and pyS is the expected proportion of the S that are counted. p ys riy S individuals could change over time depending upon observer, and could also vary stochastically (giving Cy an overdispersed binomial distribution). S Simulation studies are a useful tool for investigating the properties of the methods under known conditions, but they can be contentious because the results depend upon the model and the range of parameter values used to generate the simulated data. One solution to this problem would be to invite all interested parties (developers of methods and users of results) to co-operate in planning and executing the simulation trials. In this approach, an initial workshop would be held to agree a general reference model, range of parameters, and evaluation criteria. Many simulated datasets would then be generated by an independent coordinator and sent to the developers of the methods for analysis. The analysts would not know the exact parameters used in the simulations. These would then be returned to the coordinator, analyzed, and the results -125- discussed at a follow-up workshop. Further simulations that are more specific may be agreed upon to address questions that arise. This approach (on a much larger scale) was successfully used by the International Whaling Commission to resolve the contentious issue of whale stock management procedures (International Whaling Commission 1992). Alternative approaches Generalized additive modeling (GAM): G A M is an extension of generalized li modeling where the independent variable is nonlinearly related to the covariables using smooth functions (Hastie and Tibshirani 1990). This approach provides an additive model counterpart to loglinear modeling, allowing smooth trends in abundance over time to be modeled. In a preliminary analysis, Fewster and Buckland (1996) investigated the feasibility of using G A M s to model British Common Bird Census Data. They compared G A M with loglinear modeling and found G A M s to be superior on both theoretical and empirical grounds (comparison of trend estimates for 14 species). They suggested the use of auto-regression (in which the predicted count from the previous year is included as a covariable) to model serial dependence of counts, and also provided a method of identifying significant changes in trend (by testing for non-zero approximate second derivative of predicted count). This approach could, in theory, readily be applied to other datasets. However, no specialist software is available, and the large memory requirements of complex models using standard software (e.g., function gam in S-plus, StatSci 1995) may severely limit the amount of data that can be analyzed. Nevertheless, G A M s often require fewer interaction terms than comparable linear models (Hastie and Tibshirani 1990) so this may not be a problem in practice. As with loglinear modeling, the approach could be extended to include alternative count distributions, such as negative binomial and Poisson-logistic. Another intriguing possibility is to fit a spatio-temporal G A M for counts, allowing more precise inferences to be made about the geographic patterns of population change. State-space modeling: The observed variation in counts comes from many sources: population trend, irregular perturbations, density dependent responses and measurement error (Chapter 3). However, all the current estimators assume a form for trend and put the remaining - 126- non-trend variation into a single catch-all error term. This is inefficient and uninformative. State-space models (Harvey 1989; Aoki 1990) attempt to model and estimate these components explicitly. The approach is theoretically equivalent to traditional A R I M A time-series analysis (Box and Jenkins 1970) but offers several practical advantages, such as easier model construction, fitting and interpretation, and better handing of multivariate (i.e., multi-site) data. In state-space modeling, the system (population in this case) is described in terms of parameters which change through time according to rules described in transition equations. The system is not observed directly but is sampled, so introducing noise (error). This is modeled by an observation equation. For example, a simple model of counts has the observation equation ys c = ys * ys n e (where Cy is the observed count at time y and site s, n S ys is the actual population size, and e ys is the measurement error) and the transition equations ys n = (y-l)s* b = b .j y b -i n y y *fy * g y (where by is the regional trend in year y,fy represents stochastic interventions and gy is the change in trend over time. The distribution family of the error terms e,f and g is assumed known). Estimates of b give the trend over time, and estimates of n can be used to calculate annual indices (area weighted mean, n .). y The model could be refined by allowing b,f and g to vary by geographic stratum and by the addition of other explanatory covariables (such as weather and observer) in the same way as loglinear modeling. One practical advantage over the G L M and G A M approaches is that the estimation procedure (Kahlman filtering) is much less computer-intensive and that updating the estimates to account for a new year of data is simple (Harvey 1989). Fitting ecological population models: As Solow (1995) has commented, ecolog fond of constructing theoretical population models, but rarely attempt to fit them to real data. However, ecologically-based (as opposed to purely empirical) modeling has two important advantages. Firstly, the parameters are biologically meaningful. It is far easier to interpret estimates of, say, changes in carrying capacity or per-capita growth rate than a smoothed - 127- population trend curve. Secondly, the models are likely to be a better approximation to reality, and may therefore provide a better (i.e., more parsimonious) fit to the data. A n example of a candidate population model is the simulation model proposed in the previous section. Such a model could be fit using maximum likelihood techniques or, with the addition of suitable priors, Bayesian methods. A n example of the former approach for intensive bird count data is given in Lebreton(1990). The models could be extended to incorporate information from other data sources. For example, in North America there are at least six regional or continental surveillance programs that provide information about relative bird abundance and three that give demographic information (Butcher et al. 1993). These could be combined using a simple demographic model in either a likelihood or Bayesian framework (Hillborn and Walters 1992; Rafferty et al. 1995 and subsequent papers). This approach has the potential to provide far more information about the dynamics of wildlife populations, and begin to explicitly address the second major goal of biological monitoring programs: generating preliminary hypotheses of causation. In reality it is likely that the quality of the data will seriously constrain the inferences that can be drawn from fitting population models, especially for the more ambitious demographic models. The exercise may serve instead to highlight the parameters for which more information is required. In either case, it will be informative and is worth attempting. Implications of research To be effective, population monitoring programs must provide reliable estimates of population change. However, as I have shown in this thesis, estimating population change from extensive survey data is a considerable statistical challenge. A wide diversity of analysis methods are available, but each is problematic in some areas. I have recommended that the use of some methods be discontinued, and have attempted to provide guidelines for choosing among those that remain (Fig 5.1). Nevertheless, there are a number of outstanding issues that make it difficult to select an appropriate method in many situations. I have suggested a consensus-based simulation approach to help resolve these remaining issues, but it seems likely that there will - 128- always be more than one seemingly appropriate method available in some circumstances. In these cases I recommend using all of the candidate methods and comparing the results to see how robust they are to the different approaches. If the different methods lead to substantially different results, the conclusions about population changes should be tempered accordingly. The field of population change estimation is evolving rapidly. A l l of the methods recommended here were newly developed or significantly improved during the period of this research, and I have suggested some further developments and alternative approaches. Improvements in current methods will, I hope, resolve the practical problems that constrain their use, while new approaches will allow for increased scope of inferences (e.g., estimating separate components of the time series with state-space modeling). Nevertheless, the accuracy of the results will always be constrained by unqualified biases in the survey methodology, such as an incomplete sampling frame or unknown but systematic changes in measurement error (Chapter 3). New surveys should be carefully designed to minimize these biases (e.g., the new British Breeding Bird Survey, Gregory and Carter 1994), while existing surveys must be validated as much as possible (e.g., Bart et al. 1995). Confronted with all of the problems outlined in this thesis, the reader may be tempted to ask whether reliable estimates of population change can ever be derived from large-scale monitoring schemes such as the BBS. To some extent they already are: trends from the B B S are positively correlated with those from a wide range of independent monitoring schemes, each of which has different survey biases and uses different methods of analysis (Butcher et al. 1990; Hagan et al. 1992, Hussell et al. 1992, Cyr and Larivee 1993). Long-term, extreme changes in the population status of birds well sampled by the BBS are almost certain to be detected using any of the recommended methods of analysis. However the degree to which these schemes can meet more stringent targets for accuracy, such as those laid out by Partners in Flight (Butcher et al. 1993, Chapters 2 and 3), is unknown. Determining and realizing the full potential of extensive monitoring programs will require strong commitments by survey biologists to investigate and minimize potential sources of bias and by statisticians to evaluate and improve their methods of analysis. - 129- a _o to -a c To 3I C s s <u o Pi o O "on s a co >< W T3 oo CO > •a CO CD 60 CO T3 00 CO T3 1 > 3 OH < s o o .3 .a _a o c -a a o o CD co o 2 « b To co 3 ss -a co to o £ 6 -9 •8 <2 § ° § J CD .3 CO T3 0 *-« OH OH O - CO oo 8. o -2 rS -S o on co a J - cu (13 o • CO t3 e <D CD ao J3 o .y § > > cS o, rS -n & o TS 5J -) o TO (jI Ii j>> CO _ a cd "ob-S 1 o c <u 6 to c/) CL> e O o •3 OH xt co 4-; LO 3 o Wl QH C? o J3 CO <J S .a •a o 00 £ 00 .3 .3 CO CJ c o £ I •§ ?11 S3 ?P .3 .3 *a> "5b u •J B -130- 60 .3 •*-» 3 OH a> .3 CO CO 04 U Tcd3 T3 CQ 60 CD <2 -2ffl v. s <« ^ Ui "73 CD.3 S2 csCD cn <B uX C D r* CD2 S =§ s c n 3 .§ c i>> CO <D 3 •a c _o IS -a c is CD - 3 •a 3 >—l n OH Ci TCD3 CD CD o 3 3 g J> a "S -a 8fe-S3 cn ft o .2 <" s ft a 3 cn CD 60 CD X w CD -a a C D C D CD Ui u O cn 2 -3 CD — iS CD CD T.*3S X cr s CD -a T3 c c3a G 3 oo a CDfeO •3 § 3 cn 60 o s cn U. CD 1 ) 1) h Ccvn s•O Xca 'e13 g S •M o "P 2 °h •£ cUo3 fo 23 O CU o X }-< a a CD O CD i . CD T3 3 CD O o 1 3 Pi a. CD 2 S ft ca u, >">> X 3 CO fH CD — o J. SO 5 ft 2cn 2 fe •S fe > cn U <UI o a 60 2 cn 5 s ca .a CD ca CD -2 O T3 U CD a ° 7 a CD CO CD 60 C3 •*-» d 2 .2 § <2 > a CD 60 C _o S 3 — o fe ^ 3 O • s i "P °° a O3 u cn a & -a « cr cn ca - 3 c*-i O 2 2 c a CD Ocn -*-» 3 s -2 3 ca SC. --u T3 3 \3 «? I ? g»8Uift M P - s CD •!= cO *£3 ca 60 3 3 ui ca C D H e a -3 3 cn u3> o 60 C D C D & ft a a -g a o .a •9 8 .S 1 CD cn o o '+3 u § -9 ^ 3ca C cn cn CD o -4-» C 'C 23 c a u a ft 3 oD C 60 T3 B o cn a a 3 Cn ft l-i , p> O ci c 3 iS CD s cn nP ,~ CD C — 60 « .H n h C D _ ca cn J3 ft cn c n CD & ft S CD _ <D cn ™ 8 .3 D ca.a £ .a C cS a 3 u i > P-3 6 O T3 O0 . 4-) a> o 111 o .a T3 CO CD c I CD ount abo ft C D 60 cou ea 1 > £ is fe 3 O 3 6 0 CD o '55 3 C D .n U oo b £ S a S .3 « P ca a 3 'x .3 CM on TO 3 T3 a CD 3 CD a T3 O a is T3 3 T3 ca is CD C D C< D cn cn CL cn o ro" CN" P .3 60 TO3 3 CD o Ui CD ia D CD C j_. cn .a In p CO fa 3 | .a - 131 - T3 O Ui 3 3 •'*O c-»a g. CD 3 S3 o 3 CD 3 6CD0 3 P 55 CD is I cn S 3 cd o "C cd »-< C O CD 2 ^3 c o "3. cd | ° j=i SP 2 2 a. <u > o is o. <& 2 "S > CD 3 T3 • cd c CD a a o o CD Q a, £ CD Pi .a o 2 % "co a 4^ (D O 1/5 2? i > CD p C T3 ° C cd 1TCD3. 34sS cn r. CD CD cd o O 43 3 + CD cn ^ O CH « • ' . .o. £ § .a on "> fe cd O g< J 6 CDx & S «.a a z .a •8 a cd ca > o stributions - inc e oth eathe h trei era lized addit od AM) c>d oo o 'cn (3 cn T3 <-< CD CD X w 'C O 00 00 •3 CD UIS - o o CD T3 O 00 8 cd O CD w <D CD cn H x > JD a cd oo a cd a 00 CD S ed CD m a CD -a o .5 • h a o cd CD CD CD00 DO .3 _3CD tJ o w sx CD sx i a— WJ S CH O0 -—i O i_, q — ,o cn cd > J3 cSX d i~ 00 o CD 00 T3 •a a o CD a .3a o cd > 42 cn cn >> g C D 00 C O -4-» O ~ 3 T3 CD T3 CD CD « -9 acn cd cd cd 5 '13 T3 o> > CD 6 o c Cd Q 3 o -2 & o O 42 00 oou <cD n CD £ .a .3 • cd 43 C D ccn n 42 c3 O C CD J a •a CD o ^ (3 CX CcT cd CD 2 cn 3 >< 2 cd 3 i cd .S3 o Q .a T3 O s 2 CD ° o cd M CD .g •3 "ob ° 2 i-j a •il cd 132 Figure L e g e n d s Figure 5.1. A flow-chart guide to choosing among the analysis methods. Only methods that recommended are included in the chart. -133 - Figure 5.1 Parameter of interest? Actual yearly population changes^ r Annual index methods Long-term underlying tendency in abundance Trend methods Both Separate annual index and trend analyses or Regression of annual indices (not evaluated) Loglinear modeling with year as a factor Approach to trend analysis? Route-based /route variationDirect estimation \flexiblemodeling of 'easily incorporated, patterns of population ' low computer k change, ' requirements, v but problems with 'but problems combining large, complex models ' route trends into overall trend: Route regression (RR) Direct methods Model of trend? Model of trend? Linear I unrealistic Additive \ can be more realistic, Linear unreal-Additivepotentially istic for more realfor long time but requires degree istic long time periods, but of smoothing to be set periods works with subjectively, and much missing cannot tolerate much ' data missing data Linear RR, estimating equations implementation Smoothing RR: N N R R or NSRR - 134- Loglinear modeling with year as a linear covariate Generalized additive modeling (new method - see text) Literature Cited Anderson, D. R., K . P. Burnham, and G. C. White. 1994. AIC model selection in overdispersed capture-recapture data. Ecology 75: 1780-1793. Anganuzzi, A . A . 1993. A comparison of tests for detecting trends in abundance indices of dolphins. Fishery Bulletin 91: 183-194. Aoki, M . 1990. State Space Modeling of Time Series, 2nd Edition. Springer-Verlag, Berlin, Germany. Askins, R. A . 1993. Population trends in grassland, shrubland, and forest birds in eastern North America. Pages 1-34 in D. M . Power, editor. Current Ornithology. Volume 11. Plenum Press, New York. Askins, R. A., J. F. Lynch, and R. S. Greenberg. 1990. Population declines in migratory birds in eastern North America. Pages 1-57 in D. M . Power, editor. Current Ornithology. Volume 7. Plenum Press, New York. Bailey, R. S. 1967. A n index of bird population changes on farmland. Bird Study 14: 195-209. Baillie, S. R. 1990. Integrated population monitoring of breeding birds in Britain and Ireland. Ibis 132: 151-166. Bard, Y . 1974. Nonlinear Parameter Estimation. Academic Press, New York, New York. Barker, R. J., and J. R. Sauer. 1992. Modeling population change from time series data. Pages 182-194 in D. R. McCullough and R. Barrett, editors. Wildlife 2001: Populations. Elsevier, New York, New York. Bart, J., M . Hofschen, and B. G. Peterjohn. 1995. Reliability of the Breeding Bird Survey: effects of restricting surveys to roads. The Auk 112: 758-761. Bohning-Gaese, K., M . L. Taper, and J. H. Brown. 1993. Are declines in North American insectivorous songbirds due to causes on the breeding range? Conservation Biology 7: 76-86. Borland International Inc. 1994. Borland C++for OS/2 User's Guide, Version 2.0. Borland International Inc., Scotts Valley, C A . Box, G. E. P., and G. M . Jenkins. 1970. Time Series Analysis, Forecasting and Control. Holden-Day, San Fransisco, CA. Bradu, D., and Y . Mundlac. 1970. Estimation in lognormal linear models. Journal of the American Statistical Association 65: 198-211. Brown, J. H . 1984. On the relationship between the abundance and distribution of species. The American Naturalist 124: 255-279. Buckland, S. T. 1984. Monte Carlo confidence intervals. Biometrics 40: 811-817. -135- Burnham, K . P., and D. R. Anderson. 1992. Data-based selection of an appropriate biological model: the key to modern data analysis. Pages 16-30 in D. R. McCullough and R. Barrett, editors. Wildlife 2001: Populations. Elsevier, New York, New York. Butcher, G. S. 1990. Audubon Christmas Bird Counts. Pages 5-14 in J. R. Sauer and S. Droege, editors. Survey Designs and Statistical Methods for the Estimation of Avian Population Trends. Biological Report 90(1). U.S. Fish and Wildlife Service, Washington D.C. Butcher, G. S., M . R. Fuller, L. S. McAllister, and P. H . Geissler. 1990. A n evaluation of the Christmas Bird Count for monitoring population trends for selected species. Wildlife Society Bulletin 18: 129-139. Butcher, G. S., B . G. Peterjohn, and C. J. Ralph. 1993. Overview of national bird population monitoring programs and databases. Pages 192-203 in D. M . Finch and P. W. Stangel, editors. Status and Management of Neotropical Migratory Birds; Proceedings of the 1992 Partners in Flight National Training Workshop; September 21-25; Estes Park, Colorado. General Technical Report RM-229. U.S.D.A. For. Serv. Rocky Mountain Forest and Range Experiment Station, Fort Collins, CO. Bystrak, D. 1981. The North American Breeding Bird Survey. Pages 34-41 in C. J. Ralph and J. M . Scott, editors. Estimating Numbers of Terrestrial Birds. Studies in Avian Biology 6. Cooper Ornithological Society, Lawrence, Kansas. Cameron, A . C , and P. K . Trivedi. 1990. Regression-based tests for overdispersion in the Poisson model. Journal of Econometrics 46: 347-364. Collins, B . T. 1990. Using rerandomizing tests in route-regression analysis of avian population trends. Pages 63-70 in J. R. Sauer and S. Droege, editors. Survey Designs and Statistical Methods for the Estimation of Avian Population Trends. Biological Report 90(1). U.S. Fish and Wildlife Service, Washington D.C. Collins, B . T. 1994. Annual population indices from BBS data. Unpublished Canadian Wildlife Service internal document, March 31 1994. Collins, B . T., and J. S. Wendt. 1989. The Breeding Bird Survey in Canada, 1966-1983: Analysis of Trends in Breeding Bird Populations. Technical Report 75. Canadian Wildlife Service, Ottawa, Ontario. Cyr, A., and J. Larivee. 1993. A checklist approach for monitoring neotropical migrant birds: twenty-year trends in birds of Quebec using EPOQ. Pages 229-237 in D. M . Finch and P. W. Stangel, editors. Status and Management of Neotropical Migratory Birds; Proceedings of the 1992 Partners in Flight National Training Workshop; September 21-25; Estes Park, Colorado. General Technical Report RM-229. U.S.D.A. For. Serv. Rocky Mountain Forest and Range Experiment Station, Fort Collins, CO. Downes, C. M . , and B. T. Collins. 1996. The Canadian Breeding Bird Survey, 1966-1994. Progress Note 210. Canadian Wildlife Service, Ottawa, Ontario. Droege, S. 1990. The North American Breeding Bird Survey. Pages 1-4 in J. R. Sauer and S. Droege, editors. Survey Designs and Statistical Methods for the Estimation of Avian Population Trends. Biological Report 90(1). U.S. Fish and Wildlife Service, Washington D.C. -136- Eberhardt, L . L . 1978. Appraising variability in population studies. Journal of Wildlife Management 42: 207-238. Efron, B . 1982. The Jackknife, the Bootstrap and Other Resampling Plans. Society for Industrial and Applied Mathematics, Philadelphia, Pennsylvania. Ehrlich, P. R., D. S. Dobkin, and D. Wheye. 1988. The Birder's Handbook. Simon & Schuster Inc., New York, New York. Erskine, A . J. 1978. The First Ten Years of the Co-operative Breeding Bird Survey in Canada. Canadian Wildlife Service Report 42. Canadian Wildlife Service, Ottawa, Ontario. Erskine, A . J., B . T. Collins, E. Hayakawa, and C. M . Downes. 1992. The Cooperative Breeding Bird Survey in Canada, 1989-91. Progress Note 199. Canadian Wildlife Service, Ottawa, Ontario. Fewster, R. M . , and S. T. Buckland. 1996. Quantifying population trends from C B C counts. Unpublished contract report to British Trust for Ornithology, M A F F contract CSA3109. Finch, D. M . , and P. W. Stangel. 1993. Introduction. Pages 1-4 in D. M . Finch and P. W. Stangel, editors. Status and Management of Neotropical Migratory Birds; Proceedings of the 1992 Partners in Flight National Training Workshop; September 21-25; Estes Park, Colorado. General Technical Report RM-229. U.S.D.A. For. Serv. Rocky Mountain Forest and Range Experiment Station, Fort Collins, Colorado. Flather, C. H., and J. R. Sauer. 1996. Using landscape ecology to test hypotheses about largescale abundance patterns in migratory birds. Ecology 77: 28-35. fPrint U K Ltd. 1996. Virtual Pascal Users Guide, Version 1.0. fPrint U K Ltd., London, U K . Geissler, P. H . 1984. Estimation of animal population trends and annual indices from a survey of call-counts or other indications. Proceedings of the American Statistical Association, Section on Survey and Research Methods 1984: 472-477. Geissler, P. H., and W. A . Link. 1988. Bias of animal population trend estimates. Pages 755759 in E. J. Wegman, D. T. Gantz and J. J. Miller, editors. Computing Science and Statistics; Proceedings of the 20th Symposium on the Interface; Fairfax, Virginia; April 1988. American Statistical Society, Alexandria, Louisiana. Geissler, P. H., and B. R. Noon. 1981. Estimates of avian population trends from the North American Breeding Bird Survey. Pages 42-51 in C. J. Ralph and J. M . Scott, editors. Estimating numbers of terrestrial birds. Studies in Avian Biology 6. Cooper Ornithological Society, Lawrence, Kansas. Geissler, P. H., and J. R. Sauer. 1990. Topics in route-regression analysis. Pages 54-57 in J. R. Sauer and S. Droege, editors. Survey Designs and Statistical Methods for the Estimation of Avian Population Trends. Biological Report 90(1). U.S. Fish and Wildlife Service, Washington D.C. Godambe, V . P. 1991. Estimating Functions. Clarendon Press, Oxford, U K . Goldstein, H . 1995. Multilevel Statistical Models. Edward Arnold, London, U K . -137- Greenwood, J. J. D. 1991. Notes on simulations to investigate 'random walk' in C B C and WBS. Unpublished British Trust for Ornithology internal document, June 13 1991. Greenwood, J. J. D., S. R. Baillie, H. Q. P. Crick, J. H. Marchant, and W. J. Peach. 1993. Integrated population monitoring: detecting the effects of diverse changes. Pages 267-342 in R. W. Furness and J. J. D. Greenwood, editors. Birds as Monitors of Environmental Change. Chapman and Hall, London, U.K. Greenwood, J. J. D., S. R. Baillie, R. D. Gregory, W. J. Peach, and R. J. Fuller. 1994. Some new approaches to conservation monitoring of British breeding birds. Ibis 137: S16-S28. Gregory, R. D., and S. Carter. 1994. Results from the Pilot Census Project. BTO News 194: 1213. Gustafson, P. 1997. Large hierarchical Bayesian analysis of multivariate survival data. Biometrics 53: 230-242. Hagan, J. M . 1992. Conservation biology when there is no crisis - yet. Conservation Biology 6: 475-476. Hagan, J. M . , T. L . Lloyd-Evans, J. L . Atwood, and D. S. Wood. 1992. Long-term changes in migratory landbirds in the north-eastern United States: evidence from migration capture data. Pages 115-130 in J. M . Hagan and D. W. Johnston, editors. Ecology and Conservation of Neotropical Migrant Landbirds. Smithsonian Institution Press, Washington D.C. Hart, J. D., and T. E. Wehrly. 1986. Kernel regression estimation using repeated measurements data. Journal of the American Statistical Association 81: 1080-1088. Harvey, A . C. 1989. Forecasting, structural time series models and the Kalman filter. Cambridge University Press, Cambridge, U K . Hastie, T. J., and R. J. Tibshirani. 1990. Generalized Additive Models. Chapman and Hall, London, U K . Heckert, J. R. 1997. Bobolink Dolichonyx oryzivorus population decline in agricultural landscapes in the midwestern USA. Biological Conservation 80: 107-112. Hellawell, J. M . 1991. Development of a rationale for monitoring. Pages 1-14 in F. B . Goldsmith, editor. Monitoring for Conservation and Ecology. Chapman and Hall, London, U.K. Hilbe, J. H . 1994. Log negative binomial regression using the G E N M O D procedure SAS/STAT software. Pages 1199-1204 in Proceedings of the Ninteenth Annual SAS Users Groups International Conference. SAS Institute, Inc., Cary, North Carolina. Hillborn, R., and C. J. Walters. 1992. Quantitative Fisheries Stock Assessment: Choice, Dynamics and Uncertainty. Chapman and Hall, London, U K . Holmes, R. T., and T. W. Sherry. 1988. Assessing population trends of New Hampshire forest birds: Local vs. regional patterns. The Auk 105: 756-768. -138- Hunter, W. C , M . F. Carter, D. N . Pashley, and K. Barker. 1993. The Partners in Flight species prioritization scheme. Pages 109-120 in D. M . Finch and P. W. Stangel, editors. Status and Management of Neotropical Migratory Birds; Proceedings of the 1992 Partners in Flight National Training Workshop; September 21-25; Estes Park, Colorado. General Technical Report Rm-229. U.S.D.A. For. Serv. Rocky Mountain Forest and Range Experiment Station, Fort Collins, Colorado. Hussell, D. J. T., M . H . Mather, and P. H . Sinclair. 1992. Trends in numbers of tropical- and temperate-wintering migrant landbirds in migration at Long Point, Ontario, 1961-1988. Pages 101-114 in J. M . Hagan and D. W. Johnston, editors. Ecology and Conservation of Neotropical Migrant Landbirds. Smithsonian Institution Press, Washington, D.C. International Whaling Commission. 1992. Report of the Sub-Committee on Management Procedures. Pages 87-136 in Forty-Second Report of the International Whaling Commission. International Whaling Commission, Cambridge, U K . James, F. C , and C. E. McCulloch. 1995. The strength of inferences about causes of trends in bird populations. Pages 40-51 in D. M . Finch and T. E. Martin, editors. Ecology and Management of Neotropical Migrant Birds: a Synthesis and Review of the Critical Issues. Oxford University Press, Oxford, U.K. James, F. C , C. E. McCulloch, and L. E. Wolfe. 1990. Methodological issues in the estimation of trends in bird populations with an example: the Pine Warbler. Pages 84-98 in J. R. Sauer and S. Droege, editors. Survey Designs and Statistical Methods for the Estimation of Avian Population Trends. Biological Report 90(1). U.S. Fish and Wildlife Service, Washington D.C. James, F. C , D. A . Wiedenfeld, and C. E. McCulloch. 1992. Trends in breeding populations of warblers: declines in the southern highlands and increases in the lowlands. Pages 43-56 in J. M . Hagan and D. W. Johnston, editors. Ecology and Conservation of Neotropical Migrant Landbirds. Smithsonian Institution Press, Washington D.C. James, F. C , C. E. McCulloch, and D. A. Wiedenfeld. 1996. New approaches to the analysis of population trends in land birds. Ecology 77(1): 13-27. Jassby, A . D., and T. M . Powell. 1990. Detecting changes in ecological time series. Ecology 71: 2044-2052. Johnson, D. H., and M . D. Schwartz. 1993. The Conservation Reserve Program and grassland birds. Conservation Biology 7: 934-937. Kendall, W. L., B . G. Peterjohn, and J. R. Sauer. 1996. First-time observer effects in the North American Breeding Bird Survey. The Auk 113: 823-829. Krebs, C. J. 1985. Ecology: The Experimental Analysis of Distribution and Abundance, 3rd Edition. Harper & Row, New York, New York. Lacy, R. C , and C. E. Bock. 1986. The correlation between range size and local abundance of some North American birds. Ecology 67: 258-260. Lawton, J. H . 1993. Range, population abundance and conservation. Trends in Evolution and Ecology 8: 409-413. - 139- Lebreton, Jean-D. 1990. Modelling density dependence, environmental variability, and demographic stochasticity from population counts: an example using Wytham Wood Great Tits. Pages 89-102 in J. Blondel, editor. Population Biology of Passerine Birds. SpringerVerlag, Berlin, Germany. Lebreton, Jean-D., K. P. Burnham, J. Clobert, and D. R. Anderson. 1992. Modeling survival and testing biological hypotheses using marked animals: a unified approach with case studies. Ecological Monographs 62: 67-118. Lehman, E. L . 1975. Nonparametrics - Statistical Methods Based on Ranks. Holden-Day Inc., San Francisco, California. L i , X . , and N . E. Heckman. 1996. Local linear forecasting. Technical Report 167. University of British Columbia, Department of Statistics, Vancouver, Canada. Link, W. A., and J. R. Sauer. 1994. Estimating equation estimates of trends. Bird Populations 2: 23-32. Link, W. A., and J. R. Sauer. 1994. Within-site variability in surveys of wildlife populations. Ecology 75: 1097-1108. Link, W. A., and J. R. Sauer. 1995. Estimation and confidence intervals for empirical mixing distributions. Biometrics 51: 810-821. Link, W. A., and J. R. Sauer. 1996. Extremes in ecology: avoiding the misleading effects of sampling variation in summary analyses. Ecology 77: 1633-1640. Manly, B . F. J. 1991. Randomization and Monte Carlo Methods in Biology. Chapman and Hall, London, U.K. Marchant, J. 1990. Common Bird Census report 1988-1989. Pages 52-61 in D. Stroud and D. Glue, editors. Britain's Birds in 1989-90: the Conservation and Monitoring Review. British Trust for Ornithology, Thetford, U . K . McCullagh, P., and J. A . Nelder. 1989. Generalized Linear Models, 2nd Edition. Chapman and Hall, London, U K . Morton, E. S. 1992. What do we know about the future of migrant landbirds? Pages 75-84 in J. M . Hagan and D. W. Johnston, editors. Ecology and Conservation of Neotropical Migrant Landbirds. Smithsonian Institution Press, Washington D.C. Moses, L . E., and D. Rabinowitz. 1990. Estimating (relative) species abundance from route counts of the Breeding Bird Survey. Pages 71-79 in J. R. Sauer and S. Droege, editors. Survey Designs and Statistical Methods for the Estimation of Avian Population Trends. Biological Report 90(1). U.S. Fish and Wildlife Service, Washington D.C. Moss, D. 1985. Some statistical checks on the BTO Common Birds Census index - 20 years on. Pages 175-179 in K . Taylor, R. J. Fuller and P. C. Lack, editors. Bird Census and Atlas Studies, Proceedings of the CHI International Conference on Bird Census and Atlas Work. British Trust for Ornithology, Tring, U K . Mountford, M . D. 1982. Estimation of population fluctuations with application to the Common Bird Census. Applied Statistics 31: 135-143. - 140- Mountford, M . D. 1985. A n index of population change with application to the Common Bird Census. Pages 121-132 in B. J. T. Morgan and P. M . North, editors. Statistics in Ornithology. Lecture Notes in Statistics 29. Springer-Verlag, Berlin, Germany. O'Connor, R. J. 1992a. Population variation in relation to migrancy status in some North American birds. Pages 64-74 in J. M . Hagan and D. W. Johnston, editors. Ecology and Conservation of Neotropical Migrant Landbirds. Smithsonian Institution Press, Washington D.C. O'Connor, R. J. 1992b. Trends in populations: introduction. Pages 23-25 in J. M . Hagan and D. W. Johnston, editors. Ecology and Conservation of Neotropical Migrant Landbirds. Smithsonian Institution Press, Washington. Ogilvie, M . A . 1967. Population changes and mortality of the mute swan in Britain. Wildfowl 18: 64-73. Pannekoek, J., and A . J. van Strien. 1996. T R I M (TRends & Indices for Monitoring data). Research paper 9634. Statistics Netherlands, Voorburg, The Netherlands. Peach, W. J., and S. R. Baillie. 1994. Implementation of the Mountford indexing method for the Common Bird Census. Pages 653-662 in E. J. M . Hagemeijer and T. J. Verstrael, editors. Bird Numbers 1992. Distribution, monitoring and ecological aspects. Proceedings of the 12th International Conference of International Bird Census Committee and European Ornithological Atlas Committee. Statistics Netherlands, Voorburg, The Netherlands. Peterjohn, B . G., and J. R. Sauer. 1993. North American Breeding Bird Survey annual summary 1990-1991. Bird Populations 1:1-15. Peterjohn, B . G., and J. R. Sauer. 1995. Population trends of the loggerhead shrike from the North American Breeding Bird Survey. Pages 117-121 in R. Yosef and F. Lohrer E, editors. Shrikes (Laniidae) of the world: biology and conservation. Proceedings of the First International Shrike Symposium. Proceedings of the Western Foundation of Vertebrate Zoology 6. Peterjohn, B. G., J. R. Sauer, and W. A . Link. 1994. The 1992-1993 summary of the North American Breeding Bird Survey. Bird Populations 2: 46-61. Peterjohn, B . G., J. R. Sauer, and W. A . Link. 1996. The 1994 and 1995 summary of the North American Breeding Bird Survey. Bird Populations 3: 48-66. Peterjohn, B . G., J. R. Sauer, and C. S. Robbins. 1995. Population trends from the North American Breeding Bird Survey. Pages 3-39 in D. M . Finch and T. E. Martin, editors. Ecology and Management of Neotropical Migrant Birds: a Synthesis and Review of the Critical Issues. Oxford University Press, Oxford, U . K . Press, W. H . , S. A . Teukolsky, W. T. Vetterling, and B. P. Flannery. 1992. Numerical Recipes in C: the Art of Scientific Computing, 2nd Edition. Cambridge University Press, Cambridge, U.K. Rafferty, A . E., G. H . Givens, and J. E. Zeh. 1995. Inference from a deterministic population dynamics model for Bowhead Whales (with discussion). Journal of the American Statistical Association 90: 402-430. -141 - Robbins, C. S., and W. T. van Velzen. 1970. Progress report on the North American Breeding Bird Survey. Pages 22-30 in S. Svensson, editor. Bird Census Work and Environmental Monitoring. Bulletin from the Ecological Research Committee 9. Swedish Natural Science Research Council, Stockholm, Sweden. Robbins, C. S., D. Bystrak, and P. H . Geissler. 1986. The Breeding Bird Survey: Its First Fifteen Years, 1965-1979. Resource Publication 157. U.S. Fish and Wildlife Service, Washington D.C. Robbins, C. S., J. R. Sauer, R. S. Greenberg, and S. Droege. 1989. Population declines in North American birds that migrate to the neotropics. Proceedings of the National Academy of Science, U.S.A. 86: 7658-7662. Robbins, C. S., J. W. Fitzpatrick, and P. B. Hamel. 1992. A warbler in trouble: Dendroica cerulea. Pages 549-578 in J. M . Hagan and D. W. Johnston, editors. Ecology and conservation of neotropical migrant landbirds. Smithsonian Institution Press, Washington, DC. Sakamoto, Y . , M . Ishiguro, and G. Kitagawa. 1986. Akaike Information Criterion Statistics. Mathematics and Its Applications (Japanese Series). K T K Scientific Publishers, Tokyo, Japan. SAS Institute Inc. 1993a. The G E N M O D Procedure, Release 6.09. Technical Report P-243. SAS Institute Inc., Cary, North Carolina. SAS Institute Inc. 1993b. SAS Companion for the OS/2 Environment, Version 6, 2nd Edition. SAS Institute Inc., Cary, North Carolina. Sauer, J. R., and S. Droege. 1990a. Recent population trends of the Eastern Bluebird. Wilson Bulletin 102: 239-252. Sauer, J. R., and S. Droege, editors. 1990b. Survey Designs and Statistical Methods for the Estimation of Animal Population Trends. Biological Report 90(1). U.S. Fish and Wildlife Service, Washington D.C. Sauer, J. R., and S. Droege, editors. 1990c. Scissor-tailed Flycatcher analysis. Pages 157-166 in J. R. Sauer and S. Droege, editors. Survey Designs and Statistical Methods for the Estimation of Avian Population Trends. Biological Report 90(1). U.S. Fish and Wildlife Service, Washington D.C. Sauer, J. R., and S. Droege. 1990d. Wood Duck population trends from the North American Breeding Bird Survey. Pages 225-231 in L. H. Fredrickson, G. V . Burger, S. P. Havera, D. A . Graber, R. E. Kirby and T. S. Taylor, editors. Proceedings of the 1988 North American Wood Duck Symposium, St. Louis, Missouri. Sauer, J. R., and S. Droege. 1992. Geographic patterns in population trends of neotropical migrants in North America. Pages 26-42 in J. M . Hagan and D. W. Johnston, editors. Ecology and Conservation of Neotropical Migrant Landbirds. Smithsonian Institution Press, Washington D.C. Sauer, J. R., and P. H . Geissler. 1990. Estimation of annual indices from roadside counts. Pages 58-62 in J. R. Sauer and S. Droege, editors. Survey Designs and Statistical Methods for the Estimation of Avian Population Trends. Biological Report 90(1). U.S. Fish and Wildlife Service, Washington D.C. - 142- Sauer, J. R., D. D. Dolton, and S. Droege. 1994a. Mourning dove population trend estimates from Call-Count and North American Breeding Bird Surveys. Journal of Wildlife Management 58: 506-515. Sauer, J. R., B . G. Peterjohn, and W. A . Link. 1994b. Observer differences in the North American Breeding Bird Survey. The Auk 111: 50-62. Sauer, J. R., G. W. Pendleton, and B. G. Peterjohn. 1996. Evaluating causes of popualtion changes in North American insectivorous songbirds. Ecology 10: 465-478. Searle, S. R., F. M . Speed, and G. A . Milliken. 1980. Population marginal means in the linear model: an alternative to least squares means. The American Statistician 34: 216-221. Sokal, R. R., and J. F. Rohlf. 1981. Biometry. W H Freeman and Co., New York, New York. Solow, A . R. 1995. Fitting population models to time series data. Pages 21-27 in T. M . Powell and J. H . Steele, editors. Ecological Time Series. Chapman and Hall, New York, New York. StatSci. 1995. S-PLUS Guide to Statistical and Mathematician Analysis, Version 3.3. StatSci, a division of Mathsoft, Inc., Seattle, WA. Taub, S. R. 1990. Smoothed scatterplot analysis of long-term Breeding Bird Census data. Pages 80-83 in J. R. Sauer and S. Droege, editors. Survey Designs and Statistical Methods for the Estimation of Avian Population Trends. Biological Report 90(1). U.S. Fish and Wildlife Service, Washington D.C. Terborgh, J. 1989. Where Have A l l the Birds Gone?: Essays on the Biology and Conservation of Birds That Migrate to the American Tropics. Princeton University Press, Princeton, New Jersey. Terborgh, J. 1992. Why American songbirds are vanishing. Scientific American 1992(May): 98-104. ter Braak, C. J. F., A . J. van Strien, R. Meijer, and T. J. Verstrael. 1994. Analysis of monitoring data with many missing values: which method? Pages 663-673 in E. J. M . Hagemeijer and T. J. Verstrael, editors. Bird Numbers 1992. Distribution, monitoring and ecological aspects. Proceedings of the 12th International Conference of International Bird Census Committee and European Ornithological Atlas Committee. Statistics Netherlands, Vooburg, The Netherlands. Titus, K . 1990. Trends in counts of scissor-tailed flycatchers based on a nonparametric ranktrend analysis. Pages 164-166 in J. R. Sauer and S. Droege, editors. Survey Designs and Statistical Methods for the Estimation of Avian Population Trends. Biological Report 90(1). U.S. Fish and Wildlife Service, Washington D.C. Titus, K., M . Fuller, and D. Jacobs. 1990. Detecting trends in hawk migration count data. Pages 105-113 in J. R. Sauer and S. Droege, editors. Survey Designs and Statistical Methods for the Estimation of Avian Population Trends. Biological Report 90(1). U.S. Fish and Wildlife Service, Washington D.C. Underhill, L. G., and R. P. Prys-Jones. 1994. Index numbers for waterbird populations. I: Review and methodology. Journal of Applied Ecology 31: 463-480. - 143- U.S. Fish and Wildlife Service. 1995. Nongame Birds of Management Concern - The 1995 List. U.S. Fish and Wildlife Service, Office of Migratory Bird Management, Washington, D.C. (Available on the internet at http://www.fws.gov/~r9mbmo/reports/speccon/tblconts.html) Venables, W. N . , and B. D. Ripley. 1994. Modern Applied Statistics with S-Plus. Statistics and Computing. Springer-Verlag, New York, New York. Villard, Marc-A., and B. A . Maurer. 1996. Geostatistics as a tool for examining hypothesized declines in migratory songbirds. Ecology 77: 59-68. White, G. C , and R. E. Bennetts. 1996. Analysis of frequency count data using the negative binomial distribution. Ecology 77: 2549-2557. Wilcove, D. S., and J. Terborgh. 1984. Patterns of population decline in birds. American Birds 38: 10-13. Winkelmann, R. 1994. Count Data Models. Econometric Theory and an Application to Labor Mobility. Lecture Notes in Econometric and Mathematical Systems 410. Springer-Verlag, Berlin, Germany. Witham, J. W., and M . L . J. Hunter. 1992. Population trends of neotropical migrant landbirds in northern coastal New England. Pages 85-95 in J. M . Hagan and D. W. Johnston, editors. Ecology and Conservation of Neotropical Migrant Landbirds. Smithsonian Institution Press, Washington D.C. -144- Appendix 1 Computer programs for analysis of population change In this appendix I list the computer programs that I used, and give details of how they may be obtained. Most of the analyses described in this thesis cannot be performed using standard statistical software without considerable macro programming. There is one exception: loglinear modeling can be performed by any package that allows generalized linear modeling, and I did some of the loglinear analyses in Chapter 4 using SAS (v6.11, SAS Institute Inc. 1993b; see Chapter 4, Discussion). In all other cases, however, I used custom-written computer programs that I either developed myself or adapted from programs supplied by the originators of the methods (Table A l . 1). I modified the programs that were sent to me for four reasons: (1) to enable them to work with the continent-wide B B S dataset (programs N N R R , NSRR, T R I M , M O U N T F ) ; (2) to update the analysis methods (MOUNTF); (3) to alter the input, output, and error handling capabilities (NNRR and NSRR); and (4) to enable them to work with my compilers (TRIM, MOUNTF). A l l programs were compiled and run on a 90 Mhz Pentium PC with 32MB R A M operating under I B M OS/2 (v2.11 for the analyses in Chapter 2 and 3, v4 for the analyses in Chapter 4). The compilers were Borland C++ v2.0 for OS/2 (Borland International Inc. 1994), Virtual Pascal vl.O (fPrint U K Ltd. 1996), and a port of G N U Fortran (g77-0.5.19) by Michael Holzapfel (unpublished) for the emx runtime environment (version 0.9c fix 2, Eberhard Mattes, unpublished; available from ftp://ftp.hobbes.nmsu.edu/pub/os2/dev/emx/). I will supply the source code or complied versions of the programs that I wrote upon request; however I cannot distribute source code for any programs that are based on the work of others without their prior consent. The names and addresses of the original authors are given in Table A L L - 145- cu Td t— 1 1 , 3 cu •CU 00 Ul cu t ui co ' S i o cu 13 Bu o CO CN CN o" o o CO CO cu l-l Td < -t-» in CD «" o oo 3 g o PQ •a o O c/f co cd o T3 8 cu cS CO 3 '-id CU CU CU to a CO SH (H UI CU CU CD cd* cu cd* cu cd* cu CU CU • r l CN PQ cu co r . H ^ 7 d CN o cu co ffl co £ Cd g '5b Td C+H CU CO O o C oO O -a T3 CO Cd a U o -a 13 -71 ca * £ £ r-H Q PQ Q J J CD cd o t-- cd cd o o i H H p., 00 w -a cd a o H HH UH UH CO 00 co a O g £ ^ 8 T3 ^ 00 w , cd M CU CD C PC !>• r Ji C cd co ~ CP CU Td (H co "cd 1 1 1 i ? l o o o .5 ° £ % £ £ GO T3 Jd Cd •a a S3 2 § CD TH c3 co r^H O cd u< PH PH Z CD t 3 g tcO g o o co i3 -2 >, ^ 8 S £ s CD co CH O cd a ^ Td --H ^ c 2 ^CD Ul CD a Ja s u cn SH PQ " W CO J ° "§ <+H CO X || | Jg t2 II O g •"2 S id Td § o o Td C CH f~~- 00 Cd o cd CO I a cn CH ^ H CD Td co ^ cd co av " o 5 3 * o -g •a > Ul + + o cd CO ^ H O fe + + PH -u. 2 b ^ 1 H - eg H ii CU CU O t3 S3 O ^ ^ >H & Vi .2 bo PH CU ^3 ^3 r£ Vi CO cd CU ^ cd td CD T d X ^ JH cu -S i » -*-> P H o tS co <H Td fe H CO PQ o CO ^ PH CO fe P H 5 S CO CO www co co - 146- * H PQ M PQ Z co < P II u <H 2P i & 3 UH CD T d C cd .c co ^PH O C Ul PH 1 O _ g ra Appendix 2 Breeding Bird Survey data The Breeding Bird Survey data are in the public domain, and are available upon request to the Canadian Wildlife Service (Canadian data) or the Biological Resources Division, U.S. Geological Survey (North American data). Addresses are given below. The data used in this thesis were taken from the 1966-1995 dataset, supplied on C D R O M by the Biological Resources Division, U.S. Geological Survey. However, not all of the 4191 routes in this dataset should be analyzed (B. G. Peterjohn, pers. comm.). A small number (219) have been established in non-random locations (such as within national parks) or use nonstandard methods (such as surveying away from roads). These are recognizable in the dataset because they are coded with route numbers of 800 or greater. In addition 48 routes were established in northern Mexico in 1994 and have not yet built up enough data to be useful for trend analysis. These routes were excluded from the dataset before any analysis was performed, leaving 3924 routes (Fig. 1.1). For the analyses in Chapter 4,1 excluded a further 260 routes that are located in states or provinces that contain physiographic strata where the density of routes is so low that a meaningful analysis of population change is not possible (Table A2.1). Another 3 routes are the only ones in a physiographic stratum within a state, which makes estimation of between-trend variance impossible using linear route regression methods (Table A2.1). This left 3664 routes that were submitted to the data selection criteria outlined in Table 4.1 (see also Table 4.2). Addresses for obtaining data Canadian Wildlife Service, Migratory Bird Surveys Division, National Wildlife Research Centre, Ottawa, Ontario K 1 A 0H3. Canada. Biological Resources Division, U.S. Geological Survey, Breeding Bird Survey, Patuxent Environmental Science Center, 12100 Beech Forest Road, Laured M D 20708. U S A . -147- Table A2.1. Number of Breeding Bird Survey routes and surveys (1966 - 1995) in physiographic strata excluded from the case studies of Chapter 4. State(s) / Province(s) n. routes n. surveys Tennessee 1 3 Rhode Island 1 8 14 - Upper Coastal Plain Virginia 1 11 25 - Open Boreal Forest All 44 170 29 - Closed Boreal Forest All 107 769 68 - Northern Rockies All 49 236 94 - N . Pacific Rainforest Alaska 18 94 96 - S. Alaska Coast Alaska 20 117 All 19 53 Total excluded 260 1461 Total remaining in BBS dataset after exclusion of these data 3664 54027 Physiographic stratum 3 - Coastal Flatwoods 9 - Glaciated Coastal Plain 99 - Tundra - 148-
- Library Home /
- Search Collections /
- Open Collections /
- Browse Collections /
- UBC Theses and Dissertations /
- Evaluation of statistical methods for estimating long-term...
Open Collections
UBC Theses and Dissertations
Featured Collection
UBC Theses and Dissertations
Evaluation of statistical methods for estimating long-term population change from extensive wildlife… Thomas, Len 1997
pdf
Page Metadata
Item Metadata
Title | Evaluation of statistical methods for estimating long-term population change from extensive wildlife surveys |
Creator |
Thomas, Len |
Date Issued | 1997 |
Description | Monitoring long-term changes in population abundance is an integral part of conservation-oriented research and management. Many extensive monitoring programs are based on annual counts at a large number of permanent survey sites. However, estimating population change from the resulting data is challenging and controversial. At least eight statistical methods have been proposed and used, and in this thesis I evaluated these methods using data from the North American Breeding Bird Survey (BBS). I began by demonstrating the importance of analysis method for the inferences drawn about species trends. I compared the results when three methods of trend analysis were applied to BBS data for 115 species in British Columbia, and found differences in the direction, magnitude and statistical significance of the estimates. I then explored the issues surrounding the analysis of these data. I showed that even under ideal survey conditions subjective decisions must be made about whether to estimate annual indices or population trends and, for trend methods, which model to use for trend. In reality the analysis is further complicated by measurement error and missing data. I suggested the use of simulation studies to determine the accuracy of the methods, and provided a simple example using three geographic scenarios of population decline. I next made a preliminary assessment of the reliability of eight candidate analysis methods by applying them to continental data for four case study species (Cerulean Warbler, Carolina Wren, Chestnut-collared Longspur, and Belted Kingfisher). All methods except one (Mountford moving windows) produced seemingly plausible results in most cases. However each method demonstrated obvious flaws in some cases, and I used these to suggest potential areas of further methodological development. I concluded that loglinear modeling is currently the best method for estimating annual indices. For trend analysis, I could not determine whether route regression or direct estimation (based on generalized modeling) is generally superior. Large, long-term population changes in well-surveyed species are almost certain to be detected by any of these methods. However, their accuracy and robustness remains uncertain, and further research is required both to develop the methods and to evaluate them through comprehensive simulation studies. |
Extent | 8730589 bytes |
Genre |
Thesis/Dissertation |
Type |
Text |
File Format | application/pdf |
Language | eng |
Date Available | 2009-04-17 |
Provider | Vancouver : University of British Columbia Library |
Rights | For non-commercial purposes only, such as research, private study and education. Additional conditions apply, see Terms of Use https://open.library.ubc.ca/terms_of_use. |
DOI | 10.14288/1.0088110 |
URI | http://hdl.handle.net/2429/7326 |
Degree |
Doctor of Philosophy - PhD |
Program |
Forestry |
Affiliation |
Forestry, Faculty of |
Degree Grantor | University of British Columbia |
Graduation Date | 1997-11 |
Campus |
UBCV |
Scholarly Level | Graduate |
Aggregated Source Repository | DSpace |
Download
- Media
- 831-ubc_1997-251748.pdf [ 8.33MB ]
- Metadata
- JSON: 831-1.0088110.json
- JSON-LD: 831-1.0088110-ld.json
- RDF/XML (Pretty): 831-1.0088110-rdf.xml
- RDF/JSON: 831-1.0088110-rdf.json
- Turtle: 831-1.0088110-turtle.txt
- N-Triples: 831-1.0088110-rdf-ntriples.txt
- Original Record: 831-1.0088110-source.json
- Full Text
- 831-1.0088110-fulltext.txt
- Citation
- 831-1.0088110.ris
Full Text
Cite
Citation Scheme:
Usage Statistics
Share
Embed
Customize your widget with the following options, then copy and paste the code below into the HTML
of your page to embed this item in your website.
<div id="ubcOpenCollectionsWidgetDisplay">
<script id="ubcOpenCollectionsWidget"
src="{[{embed.src}]}"
data-item="{[{embed.item}]}"
data-collection="{[{embed.collection}]}"
data-metadata="{[{embed.showMetadata}]}"
data-width="{[{embed.width}]}"
async >
</script>
</div>
Our image viewer uses the IIIF 2.0 standard.
To load this item in other compatible viewers, use this url:
http://iiif.library.ubc.ca/presentation/dsp.831.1-0088110/manifest