Trend of the Yellowstone Grizzly Bear Population

Abstract

Yellowstone's grizzlies (Ursus arctos) have been studied for more than 40 years. Radiotelemetry has been used to obtain estimates of the rate of increase of the population, with results reported by Schwartz et al. (2006). Counts of females with cubs-of-the-year "unduplicated" also provide an index of abundance and are the primary subject of this report. An exponential model was fitted to n=24 such counts, using nonlinear leastsquares. Estimates of the rate of increase, r, were about 0.053. 95% confidence intervals, were obtained by several different methods, and all had lower limits substantially above zero, indicating that the population has been increasing steadily, in contrast to the results of Schwartz et al. (2006), which could not exclude a decreasing population. The grizzly data have been repeatedly mis-used in current literature for reasons explained here.

L. L. Eberhardt 1 and J. M. Breiwick 2

Recommended by Mats Olsson

1, 2528 W. Klamath Avenue, Kennewick, WA 99336, USA
2, National Marine Fisheries Service, National Marine Mammal Laboratory, 7600 Sand Point Way NE, Bldg. 4, Seattle, WA 98115, USA

Received 4 March 2010; Accepted 8 April 2010

1. Introduction

Yellowstone's grizzlies have been studied extensively for over 40 years. A recent monograph [1] provided a detailed analysis of much of the accumulated data. Two methods have been used to estimate trends of the population. One, the "unduplicated" counts of females with cubs-of-the-year, was developed by Dr. R. R. Knight and has been maintained since 1973. The second approach uses radiotelemetry data and was pioneered by Craighead [2]. The "unduplicated" counts along with some earlier observations were converted to an index described as a "minimum population" estimate [3]. That index used a three-year moving average to smooth the annual numbers of adult female bears and a calculated proportion of adult females in the population of 0.274. A long series of "study committees" [4] continued use of the index because a "minimum population" estimate appeared to be something readily understood by administrators and the public.

A weakness in the "minimum population" index is the presence of a strong serial correlation, induced by the fact that each annual estimate contains the counts from 2 other years. Difficulties caused by serial correlation have been assessed by Watt [5, 6], Chapman[7], and Freckleton et al. [8]. Eberhardt [9] indicated the sizable degree of correlation thus induced when the underlying data are constructed only from random numbers and showed that a correlation of -0.707 could be induced in analysis of two independent sets of random numbers. Royama [10] reported the same result. The "minimum population" index is not used here, and the Durbin-Watson test [11] does not show evidence of serial correlation in the data used here.

An alternative approach using radiotelemetry was applied by Eberhardt, et al. [12] using an approximation to Lotka's equation proposed by Eberhardt [13]. That analysis yielded an estimate of λ of 1.046 with bootstrap 95% confidence limits of 1.00 to 1.09. That study depended on bears caught and marked in a set of "backcountry" traps, located away from centers of human activity. The major, and essentially the only, cause of adult grizzly mortality comes from conflicts with humans. When such conflicts arise, efforts are made to capture and move bears. Much experience has shown that bears involved in such a situation have a low survival rate. Consequently, bears first captured in "conflict" situations were not used in the analysis, but bears first-captured in the so-called "research" trapping continued to be used in the analysis after they were involved in conflicts with humans. Subsequently, Pease and Mattson [14] proposed that there are "wary" and "unwary" bears so that the "unwary" bears first-captured in "conflict" situations should be included in the telemetry sample. This, of course, reduced the lower confidence limit on population trend below 1.0. The issues thus involved have been extensively explored in the monograph by Schwartz et al. [1] and need not be further considered here. However, probability distributions generated by stochastic simulations of Yellowstone grizzly numbers [15] have lower limits below λ =1.00, thus not excluding the prospect that the population has not increased. Again, these results have been amply examined in the monograph of Schwartz et al. [1]. The purpose of the present study is to utilize the original "unduplicated" counts as an independent measure of population trend and to correct some erroneous uses of the data.

2. Study Area

The study area constitutes Yellowstone and Grand Teton National Parks, 6 adjacent National Forests plus some state and private lands, and totals about 34,500 km² . It is known as "The Greater Yellowstone Ecosystem" (GYE). The GYE constitutes the Yellowstone Plateau, and surrounding mountain ranges above 1500 m. Long cold winters and short summers characterize the climate. Low elevations are covered by grasslands or shrub steppes. Douglas fir (Pseudotsuga menziesii ) and lodgepole pine (Pinus contorta ) are the dominant tree species. A detailed description of the study area is available in Schwartz et al. [1].

3. Material and Methods

The model used here is the exponential [figure omitted; refer to PDF] r is the annual rate of increase and t is time in years, while a denotes an initial value. This equation was fit using nonlinear least-squares [16], as implemented in the R-language [17], and ordinary least-squares regression after a logtransformation of the trend index. Utility of the log-transformation was demonstrated by Eberhardt [13]. The data were also examined using a smoothing technique, "lowess" [18, 19]. Confidence intervals on r were obtained by bootstrapping [19, 20] for the nonlinear model and from the usual linear regression model for the log-transformed data. Several different confidence intervals were calculated in the R-language program boot.ci, as expounded by [20, chapter 5] A generalized linear model (GLM) [21] with Poisson errors was also fit to the data using program glm in the R-language. All bootstraps used 5,000 calculations.

The Durbin-Watson test [11] was used to test for serial correlation. This test depends on the fact that the squared difference between successive deviations can be used to approximate the variance of the deviations if the pattern of deviations is random. The test is [figure omitted; refer to PDF] where the observations are deviations from a fitted model. Values of d lie between 0 and 4, with zero correlation at d=2 . Bias can be estimated in bootstrapping, using the following: [figure omitted; refer to PDF]

where ^θ* is the mean of the bootstraps and t^{F *} is the original estimate [19].

Data used here on the "unduplicated" females with cubs in the GYE are those given in Table 4 of the 2006 Annual Report of the Interagency Grizzly Bear Study Team, C. C. Schwartz, M. A. Haroldson, and K. West, Editors (available at http://www.nrmsc.usgs.gov/research/igbst-home.htm).

4. Results

The unduplicated data set was initially assessed by the lowess program, which is essentially a data smoothing program that calculates values from a weighted linear regression of a subset of points adjacent to the current location. It was used to produce Figure 1.

Figure 1: Plot of "unduplicated" data with lowess smoothed line.

[figure omitted; refer to PDF]

The J-shaped nature of the plot is in consequence that measures to protect and enhance grizzly populations in Yellowstone did not take effect until about 1983 [1, 3]. All subsequent analyses given here use the data from 1983 onwards. Fits of the data from 1983 onwards using nonlinear least-squares with (1) appear in Figure 2.

Figure 2: Plot of "unduplicated" data with fitted exponential function.

[figure omitted; refer to PDF]

The estimate of rate of increase (r) for the exponential fit (1) was 0.052 with a standard error of 0.0071. The Durbin-Watson test when applied to the exponential fit yielded a D-statistic of 2.32, a nonsignificant result, indicating no evidence of serial correlation. Bias calculated from (3) was negligible in all instances. Results of the various analyses appear in Table 1.

Table 1: Estimates of rate of change (r ), bias (3), and 95% confidence intervals for various methods of estimation using the exponential model.

			95% Confidence
METHOD	Estimate of r	Bias equation (4)	Interval on r
			Lower	Upper
Regression (log transformed counts)	0.0526	--	0.0422	0.07
Bootstrap (log transformed counts)	0.0563	2.00E-04	0.041	0.071
Basic bootstrap c.i.			0.04	0.071
Percentile method c.i.			0.042	0.072
BCA method c.i.			0.041	0.072
Nonlinear Least-squares exponential model	0.0526	--	0.039	0.066
Bootstrap of Nonlinear L.S. Model	0.0526	0.0007	0.037	0.066
Basic bootstrap c.i.			0.036	0.066
Percentile method c.i.			0.039	0.069
BCA method c.i.			0.036	0.066
Generalized linear model (Poisson errors)	0.054	--	0.043	0.065

5. Discussion

The fact that the difference between r from the exponential fit and 0 is over 7 times the standard error provides strong evidence that it is significantly different from zero. This result is clearly supported by the several confidence interval calculations (Table 1), so it seems clear that the "unduplicated" counts show a significant upwards trend. Two different models are considered, one the linear model with additive errors resulting from the logarithmic transformation of the data and the second a multiplicative error structure for the nonlinear least-squares fit to the exponential model. Confidence intervals from both models demonstrate that the rate of increase, r , is substantially greater than zero. We use r , the "intrinsic rate of increase" here, rather than λ = exp(r) as calculated from reproductive and survival rates, but, for small r , the difference is trivial.

The data indicate that the GYE bear population does not yet show signs of approaching an asymptote. This is a bit surprising, inasmuch as Schwartz et al. [1]. detected signs of decreasing survival in the population surrounding the primary inhabited area. One possible explanation is that the expansion zone, being crudely an annulus of substantial area, can still provide a substantial increment to the population in spite of having lower overall survival rates. Clearly, the population will ultimately level off, and its numbers may then be fit with a logistic curve. We have favored a "generalized logistic" model in which the rate of increase remains high until the asymptote (K) is approached [22], but believe the present data set is too variable for that model to be useful until the population begins to approach an asymptotic value.

Harris et al. [15]. proposed use of the Chao2 estimator [23] to "estimate the total number of females with cubs present from the estimated number observed" [15, page 17]. Using their estimates [24, Table 4 ] we fitted (1) by nonlinear least-squares giving the results shown in Figure 3. This fit had a residual standard error of 8.36 as contrasted to a value of 6.11 from the data of Figure 2, suggesting that the Chao2 estimator gives substantially more variable results than using the unadjusted estimates of females with cubs of the year.

Figure 3: Plot of Chao2 data with fitted exponential function.

[figure omitted; refer to PDF]

Conducting live-trapping, marking and radiotelemetry operations on grizzlies are very expensive operations. Hence, the main lesson from the present study is that index studies should also be attempted as a check on results. Transforming observed index data with the Chao2 estimator produced slightly more variable results, so the method may not be useful unless some improvements are possible. The radiotelemetry results provide essential data for management so should also be conducted whenever possible.

The time span of the reported data (approaching 40 years) has led to its use in a long list of applications [25-32]. All of these ignore the fact that the initial data in the series were constructed from different data sources not compatible with the 1973 development of the "unduplicated" index. As noted above, the population was evidently decreasing until about 1983. This lead to various forecasts of ultimate extinction [27, 28] a somewhat decreased risk [31,Table 1 ], a lower limit on lambda of 0.99 [26], and a very small probability of extinction [29]. Staudenmayer and Buonaccorsi [31] indicate that some of their estimates indicate a rate of increase less than zero. They also used the data to estimate extinction rates. Brook and Bradshaw [25] used the data to study densitydependence and did not consider extinction of the population.

Acknowledgments

The senior author is grateful to Dr. Richard Knight for the opportunity to work on data for the Yellowstone grizzly bears for over 20 years and to Bonnie M. Blanchard for her continual and cheerful efforts at tabulation and analysis of the data. Dave Stradley of the Gallatin Flying Service provided expert piloting for many hours of "looking for bears". Support from the Marine Mammal Commission and the National Marine Mammal Laboratory is acknowledged.

References

[1] C. C. Schwartz, M. A. Haroldson, G. C. White Temporal, Spatial and Environmental Influences on the Demographics of Grizzly Bears in the Greater Yellowstone Ecosystem , Wildlife Monographs, Washington, DC, USA, 1992.

[2] J. J. Craighead, J. Varney, F. C. Craighead A Population Analysis of the Yellowstone Grizzly Bears , Montana Forest and Conservation Experiment Station, Bulletin 40, School of Forestry, University of Montana, Missoula, Mont, USA, 1974.

[3] L. L. Eberhardt, R. R. Knight, "How many grizzlies in Yellowstone?," Journal of Wildlife Management , vol. 60, no. 2, pp. 416-421, 1996.

[4] C. Servheen, R. Schoemaker, "Delisting the Yellowstone grizzly bear," Yellowstone Science , vol. 16, pp. 25-29, 2008.

[5] K. E. F. Watt, "Density dependence in population fluctuations," Canadian Entomologist , vol. 96, pp. 1147-1148, 1964.

[6] K. E. F. Watt Ecology and Resource Management , McGraw-Hill, New York, NY, USA, 1968.

[7] D. G. Chapman, "Statistical problems in dynamics of exploited fisheries populations," in Proceedings of the 4th Berkeley Symposium on Mathematical Statistics and Probability, University of California Press, Berkley, Calif, USA, 1961.

[8] R. P. Freckleton, A. R. Watkinson, R. E. Green, W. J. Sutherland, "Census error and the detection of density dependence," Journal of Animal Ecology , vol. 75, no. 4, pp. 832-851, 2006.

[9] L. L. Eberhardt, "Correlation, regression, and density dependence," Ecology , vol. 51, pp. 306-310, 1970.

[10] T. Royama Analytical Population Dynamics , Chapman and Hall, London, UK, 1992.

[11] N. R. Draper, H. Smith Applied Regression Analysis , John Wiley & Sons, New York, NY, USA, 1998.

[12] L. L. Eberhardt, B. M. Blanchard, R. R. Knight, "Population trend of the Yellowstone grizzly bear as estimated from reproductive and survival rates," Canadian Journal of Zoology , vol. 72, no. 2, pp. 360-363, 1994.

[13] L. L. Eberhardt, "Assessing the dynamics of wild populations," Journal of Wildlife Management , vol. 49, pp. 997-1012, 1985.

[14] C. M. Pease, D. J. Mattson, "Demography of the Yellowstone grizzly bears," Ecology , vol. 80, no. 3, pp. 957-975, 1999., [email protected]

[15] R. B. Harris, C. C. Schwartz, M. A. Haroldson, G. C. White, C. C. Schwartz, M. A. Haroldson, G. C. White, "Trajectory of the yellowstone grizzly bear population under alternative survival rates," Temporal, Spatial and Environmental Influences on the Demographics of Grizzly Bears in the Greater Yellowstone Ecosystem , pp. 44-55, Wildlife Monographs, Washington, DC, USA, 2006.

[16] D. M. Bates, D. G. Watts Nonlinear Regression Analysis and Its Applications , John Wiley & Sons, New York, NY, USA, 1988.

[17] R Development Core Team, "R: A language and environment for statistical computing," R Foundation for Statistical Computing, Vienna, Austria, 2004, http://www.R-project.org

[18] W. S. Cleveland, "Robust locally-weighted regression and smoothing scatterplots," Journal of the American Statistical Association , vol. 74, no. 2, pp. 829-836, 1979.

[19] B. Efron, R. J. Tibishirani An Introduction to the Bootstrap , Chapman and Hall, New York, NY, USA, 1993.

[20] A. C. Davison, D. V. Hinkley Bootstrap Methods and Their Application , Cambridge University Press, New York, NY, USA, 1997.

[21] P. McCullagh, J. A. Nelder Generalized Linear Models , Chapman and Hall, New York, NY, USA, 1989., 2nd.

[22] L. L. Eberhardt, J. M. Breiwick, D. P. DeMaster, "Analyzing population growth curves," Oikos , vol. 117, no. 8, pp. 1240-1246, 2008.

[23] R. M. Wilson, M. F. Collins, "Capture-recapture estimation with samples of one using frequency data," Biometrika , vol. 79, pp. 543-553, 1992.

[24] R. B. Harris, G. C. White, C. C. Schwartz, M. A. Haroldson, "Population growth of Yellowstone grizzly bears: uncertainty and future monitoring," Ursus , vol. 18, no. 2, pp. 168-178, 2007.

[25] B. W. Brook, C. J. A. Bradshaw, "Strength of evidence for density dependence in abundance time series of 1198 species," Ecology , vol. 87, no. 6, pp. 1445-1451, 2006., [email protected]

[26] J. P. Buonaccorsi, J. Staudenmayer, "Statistical methods to correct for observation error in a density-independent population model," Ecological Monographs , vol. 79, no. 2, pp. 299-324, 2009.

[27] B. Dennis, P. L. Munholland, J. M. Scott, "Estimation of growth and extinction parameters for endangered species," Ecological Monographs , vol. 61, no. 2, pp. 115-143, 1991.

[28] B. Dennis, M. L. Taper, "Density dependence in time series observations of natural populations: estimation and testing," Ecological Monographs , vol. 64, no. 2, pp. 205-224, 1994.

[29] S. T. Lindley, "Estimation of population growth and extinction parameters from noisy data," Ecological Applications , vol. 13, no. 3, pp. 806-813, 2003.

[30] W. F. Morris, D. F. Doak Quantitative Conservation Biology: Theory and Practice of Population Variability Analysis , Sinauer Associates, Sunderland, Mass, USA, 2002.

[31] J. Staudenmayer, J. P. Buonaccorsi, "Measurement error in a random walk model with applications to population dynamics," Biometrics , vol. 62, no. 4, pp. 1178-1189, 2006., [email protected]; [email protected]

[32] D. F. Staples, M. L. Taper, B. Dennis, "Estimating population trend and process variation for PVA in the presence of sampling error," Ecology , vol. 85, no. 4, pp. 923-929, 2004., [email protected]