## New Research In

### Physical Sciences

### Social Sciences

#### Featured Portals

#### Articles by Topic

### Biological Sciences

#### Featured Portals

#### Articles by Topic

- Agricultural Sciences
- Anthropology
- Applied Biological Sciences
- Biochemistry
- Biophysics and Computational Biology
- Cell Biology
- Developmental Biology
- Ecology
- Environmental Sciences
- Evolution
- Genetics
- Immunology and Inflammation
- Medical Sciences
- Microbiology
- Neuroscience
- Pharmacology
- Physiology
- Plant Biology
- Population Biology
- Psychological and Cognitive Sciences
- Sustainability Science
- Systems Biology

# Black-swan events in animal populations

Edited by Stephen R. Carpenter, University of Wisconsin-Madison, Madison, WI, and approved January 18, 2017 (received for review August 3, 2016)

### This article has a reply. Please see:

## Significance

Black swans—statistically improbable events with profound consequences—happen more often than expected in financial, social, and natural systems. Our work demonstrates the rare but systematic presence of black-swan events in animal populations around the world (mostly birds, mammals, and insects). These events are predominantly downward, implying that unexpected population crashes occur more frequently than increases. Black-swan events are not driven by life history (e.g., lifespan) but by external causes such as extreme winters and disease. Ignoring the presence of downward black swans may severely underestimate extinction risk in animal populations, particularly under a changing climate, where such extreme events are expected to increase in frequency and magnitude.

## Abstract

Black swans are improbable events that nonetheless occur—often with profound consequences. Such events drive important transitions in social systems (e.g., banking collapses) and physical systems (e.g., earthquakes), and yet it remains unclear the extent to which ecological population numbers buffer or suffer from such extremes. Here, we estimate the prevalence and direction of black-swan events (heavy-tailed process noise) in 609 animal populations after accounting for population dynamics (productivity, density dependence, and typical stochasticity). We find strong evidence for black-swan events in

Major surprises (black-swan events) happen more often than expected in financial, social, and natural systems (1⇓⇓⇓–5). Massive unpredictable market swings are responsible for the majority of financial gains and losses (3), fatalities from the largest wars dwarf those from all others (6), and the frequency of the most damaging earthquakes has exceeded past expectations (5). In ecological systems, background rates of global extinction are punctuated by mass extinction (7), evolution is characterized by long periods of stasis interrupted by bursts of speciation (8), and billions of animals can die at once in mass mortality events (9).

Indeed, such die-offs may be the most important element affecting population persistence (10) and their importance is likely to increase given projected increases in the frequency and magnitude of climate-related extremes (11). However, the overwhelming majority of population model-fitting and risk-forecasting assumes that deviations from model predictions can be represented by a normal distribution [on a log scale (e.g., refs. 12 and 13)]. If black swans occur, however, a normal distribution would underestimate the probability of extreme events occurring (3).

Whereas there are many reports of black-swan events, only a flexible comparative approach consistently applied to a large number of time series can yield insights into the frequency, strength, and correlates of such events. We are unaware of such a comparative analysis. Previous comparative analyses fitting heavy-tailed distributions (distributions with higher probabilities of extreme events than the normal distribution) to time series have not accounted for an underlying population dynamics model (14⇓⇓–17). Alternatively, most examples of population die-offs come from identifying sudden changes in abundance that exceed a chosen threshold of decline (9, 18⇓–20). However, this approach does not distinguish events from expected dynamics. For instance, the range of natural variability can differ dramatically across taxa.

Here, we develop an approach to estimate the frequency and magnitude of black-swan dynamics across time series of 609 populations from a wide array of taxonomic groups—including many birds, mammals, insects, and fishes (*SI Appendix*, Table S1). We identify characteristics of time series or life-history traits associated with the detection of black-swan events and verify known causes. We develop a framework for identifying heavy-tailed (black-swan) process noise in population dynamics. We test whether the largest stochastic jumps in log abundance from one time step to the next are more extreme than typically seen with a normal distribution. Our framework allows for a range of population dynamic models, can incorporate observation uncertainty and skewness in process noise, and can be readily applied to abundance time series.

We fit population dynamics models in which the process noise is drawn from a Student

## Results

We detected black-swan dynamics most frequently for birds (7%), mammals (5%), and insects (3%) but almost never in fishes (Fig. 2). Black-swan events were taxonomically widespread, occurring in 38% of taxonomic orders. Accounting for time-series length and partially pooling inference across taxonomic class and order with a hierarchical model, we found stronger evidence for black swans in insect populations than these statistics suggest—two of five orders with the highest median probability of heavy tails were insect orders (Fig. 3*A*).

The majority of our heavy-tailed estimates were robust to alternative population models, observation error, and choice of Bayesian priors (*SI Appendix*). Our conclusions were not systematically altered when we included autocorrelation in the residuals, modeled population growth rates with or without density dependence, or modeled the population dynamics as Ricker-logistic instead of Gompertz (*SI Appendix*, Fig. S1). Similarly, introducing moderate observation error (coefficient of variation = 0.2) only slightly decreased the estimated prevalence of black-swan events (*SI Appendix*, Fig. S1), and the strength of the prior on *SI Appendix*, Fig. S2). Finally, our simulation testing shows that, if anything, our models underpredict the true magnitude and probability of heavy-tailed events—especially given that the time series are relatively short, with a median of only 26 y in our analysis (*SI Appendix*, Fig. S3 and S4).

For model fits to the population data, the probability of detecting black-swan dynamics was positively related to time-series length and negatively related to magnitude of process noise but not clearly related to population growth rate, density dependence, or maximum lifespan (Fig. 3*B* and *SI Appendix*, Fig. S5). Longer time-series length was the strongest covariate of observing black-swan dynamics: there is a 1.5-times greater probability of detecting a black-swan event in populations with 60 time steps than in one with 30 time steps (*SI Appendix*, Fig. S5).

The majority of black-swan events (86%) were downward (die-offs) rather than upward (unexpectedly rapid abundance increases). Of the black-swan events with published explanations (*SI Appendix*, Table S2), the majority involved a combination of factors. For example, a synchronization of environmental- and predation-mediated population cycles is thought to have caused a downward black-swan event for a water vole (*Arvicola amphibius*) population (21). Other black swans were the result of a sequence of extreme climate events on their own. For instance, severe winters in 1929, 1940–1942, and 1962–1963 were associated with black-swan downswings in gray herons (*Ardea cinerea*) in the United Kingdom (22) (Fig. 4*C*). Our analysis finds that the last event was a combination of two heavy-tailed events in a row and that the population took three times longer to recover than predicted (22). Downward black swans were sometimes followed by upward black swans. For example, during a period of population crowding and nest shortages, a population of European shag cormorants (*Phalacrocorax aristotelis*) on the Farne Islands, United Kingdom, declined suddenly following a red tide event in 1968 (23). This population decline freed up quality nest sites for first-time breeders, productivity increased, and the population experienced a rapid upswing in abundance (23).

Given the prevalence of downward events, we refit our heavy-tailed models to measure the degree of skewness *A*) and used these models to make near-term risk forecasts. Aggregated across populations with strong evidence of heavy tails (median *B* and *SI Appendix*, Fig. S6). In contrast, populations that did not have heavy tails (median *B*; 89% of 95% credible intervals overlapped 1). Projecting these heavy-tailed populations forward 5 y revealed that assuming standard normal process noise underestimated risk (99% lower credible interval of abundance) by 1.2– to 1.9-fold (interquartile range; Fig. 4*C* and *SI Appendix*, Fig. S7).

## Discussion

We systematically evaluated the prevalence of black-swan events in hundreds of animal populations. We find strong evidence for the occurrence of black-swan events in animal populations. Black-swan population crashes are substantially more frequent than black-swan population increases and are usually driven by external events such as weather and disease. Thus, our analysis provides strong evidence for downward-skewed heavy-tailed events in abundance time series of higher taxa, and ignoring these events will tend to underestimate the risk of population declines. Next, we consider (*i*) the possible mechanisms underlying black-swan events, (*ii*) caveats to our findings, (*iii*) the consequences of our findings for making risk forecasts, and (*iv*) for managing natural resources.

There are many possible causes of black-swan events, including unmodeled intrinsic properties of populations or extrinsic forces acting on populations. For example, we could have observed black-swan dynamics if we missed an underlying mixture of processes, because a mixture of normal distributions with different variances can generate a

There are a number of caveats when considering the generality of our results. First, the Global Population Dynamics Database (GPDD) [Version 2 (2010); Natural Environment Research Council Centre for Population Biology, Imperial College London] that we used represents a taxonomically and geographically biased sample of populations—the longer time series we focused on are dominated by commercially and recreationally important species and a disproportionate number of populations are located in the United Kingdom. Although we would expect to find qualitatively similar evidence for black swans in other large taxonomic or geographic samples of populations, the common forces driving those black swans would likely differ. Second, some apparent black-swan events could be recording mistakes, although conversely, some extreme observations may have been discarded or altered if they were erroneously suspected of being recording mistakes. Indeed, we discarded three of the populations that our method initially identified as heavy-tailed because they turned out to be data-entry errors (*SI Appendix*). Third, the temporal scales of observation and population dynamics vary considerably across populations in the GPDD, and these likely influence the detection of heavy tails. For example, if we make frequent observations relative to generation time (e.g., for many large-bodied mammals), we will average across generations and perhaps miss black-swan events. Conversely, if we census populations infrequently relative to generation time (e.g., many insects in the GPDD), the recorded data may average across extreme and less-extreme events and also dampen black-swan dynamics. Finally, our models assume that the fitted parameters (e.g., productivity and density dependence) do not vary over time. Allowing for time-varying parameters could be an important development for understanding the dynamics of longer time series.

Given our results, we suggest modelers should routinely consider representing population dynamic process noise with a heavy-tailed, and possibly downward-skewed, distribution, especially when making forecasts to evaluate risk. However, what heavy-tailed distribution should one pick and how heavy-tailed should it be? This is an open research question. The exact shape of the tails has proven important in the field of dispersal biology (32,33). Answers for the field of population dynamics might come from simulation analyses exploring the implications of tail shapes and from comparing the relative fit of heavy-tailed population models to time series with suspected mass mortality events [e.g., the *A*) if a modeler wishes to allow for downward heavy-tailed events. These are midrange values from the populations identified as heavy-tailed here. Above

In light of our findings, we suggest that natural resource management can learn from disciplines that focus on heavy tails. For example, earthquake preparedness and response is focused on black-swan events. Similarly to ecological black swans, we can rarely predict the specific timing of large earthquakes. However, earthquake preparedness involves spatial planning based on forecast probabilities to focus early detection efforts and develop disaster-response plans (35). The presence of ecological black swans also suggests that we develop management policy that is robust to heavy tails and encourages general resilience (36). For instance, setting target population abundances far back from critical limits may buffer against black-swan events (37), and maintaining genetic, phenotypic, and behavioral diversity may allow some components of populations to persist when others are affected by disease or extreme environmental forces (25). Finally, surprising, or counterintuitive, ecological dynamics offer a tremendous opportunity to learn about ecological systems, evaluate when models break down, and adjust future management policy (38, 39).

Rare catastrophes can have a profound influence on population persistence (10). In recent decades, ecology has moved toward focusing on aspects of variance in addition to mean responses (40). Our results suggest that an added focus on ecological extremes represents the next frontier, particularly in the face of increased climate extremes (11, 28, 40). Financial analysts are concerned with the shape of downward tails in financial returns because these directly impact estimates of risk—the probability of a market crash occurring. Similarly, ecologists should focus more on estimating and predicting downward tails of population abundance, because these may increase true extinction risk.

## Materials and Methods

### Data.

We selected abundance time series from the GPDD, which contains nearly 5,000 time series of abundance from *SI Appendix*) to remove populations from less reliable data sources, and those without sufficient data for our models, and then interpolated remaining missing values (sensu ref. 12). Our interpolation affected only *SI Appendix*, Table S1) and none of the data points that were later considered black-swan events. Our final dataset contained 609 populations across 39 taxonomic orders and 7 taxonomic classes, with a median of 26 time steps (range, 20–117) (*SI Appendix*, Table S1).

### Modeling Framework.

We fit heavy-tailed Gompertz population dynamics models to data from the GPDD. The Gompertz model represents population growth as a linear function in log space (here and throughout we use log to refer to the natural logarithm). If *t*, then*A* and *B*).

We chose the

One alternative approach would involve fitting the generalized extreme value (GEV) distribution, which represents the limit distribution of a series of maxima or minima. Although the GEV is well-suited for environmental variables [e.g., sedimentation rates or wind speeds (31, 41)], the GEV requires maxima per time block. Therefore, the GEV requires longer or higher frequency time series than typically available for wild animal population abundance. Furthermore, to our knowledge, the GEV cannot be easily integrated into population dynamics models.

For the Gompertz model, we chose weakly informative priors for all parameters (*SI Appendix*, Fig. S8). For the degrees of freedom parameter *E* and *H*), and the population model is otherwise a reasonable representation of the dynamics, the posterior for

We fit our models with Stan 2.14.1 (43⇓–45) and R 3.3.1 (46). We began with four chains and 2,000 iterations, discarding the first 1,000 as warmup. If

### Alternative Models and Simulation Testing.

We fit alternative population models to test if four key phenomena systematically changed our conclusions. Autocorrelation has been suggested as a reason for increased observed variability of abundance time series through time, which could create apparent heavy tails (47); therefore, we fit a model that included serial correlation in the residuals. Additionally, previous work has modeled abundance or growth rates without accounting for density dependence (16, 17); therefore, we fit a simpler model in which we assumed density independence. This model is equivalent to a random walk with drift on log abundance and therefore does not assume that the time series is stationary. Third, observation error could bias parameter estimates (13) or mask our ability to detect heavy tails (34); therefore, we fit a model where we allowed for a fixed quantity of observation error (*SI Appendix*).

We investigated the sensitivity of our results to weaker and stronger priors (exponential rate parameter*SI Appendix*, Fig. S8). Furthermore, we used simulated data to test how easily we could detect *SI Appendix*).

### Covariates of Heavy-Tailed Dynamics.

We fit a hierarchical beta regression model to the predicted probability of heavy tails, Pr(

We log-transformed *SI Appendix*, Fig. S5. All input variables were standardized by subtracting their mean and dividing by two SDs to make their coefficients comparable in magnitude (49). We excluded body length as a covariate because it was highly correlated with lifespan, and lifespan exhibited more overlap across taxonomy than body length. Lifespan is also more directly related to time and potential mechanisms driving black-swan dynamics.

We incorporated weakly informative priors into our model:

To derive taxonomic-order-level estimates of the probability of heavy tails accounting for time-series length (Fig. 3*A*), we fit a separate hierarchical model with the same structure but with only

### Skewed Student t Forecasts.

To evaluate the apparent skewness of heavy-tailed process noise, we fit Gompertz models with skewed Student

To generate 5-y forecasts of abundance, we combined the posterior parameter samples from the skewed *t* Gompertz models with stochastically generated process noise. We compared these forecasts to those generated by a standard Gompertz model with normally distributed process noise fit to the same data. We calculated the ratio of the lower 99% quantile credible interval between the two model projections for all populations in which Pr(

## Acknowledgments

We thank J. W. Moore, A. O. Mooers, L. R. Gerber, J. D. Yeakel, C. Minto, two anonymous reviewers, and members of the Earth to Ocean Group for helpful discussions and comments. We are grateful to the contributors and maintainers of the Global Population Dynamics Database and to Compute Canada’s WestGrid high-performance computing resources. Funding was provided by a Simon Fraser University Graduate Fellowship and David H. Smith Conservation Research Fellowship (to S.C.A.), the Natural Sciences and Engineering Research Council of Canada (N.K.D. and A.B.C.), the Canada Research Chairs Program (N.K.D.).

## Footnotes

- ↵
^{1}To whom correspondence should be addressed. Email: sandrsn{at}uw.edu.

Author contributions: S.C.A., T.A.B., A.B.C., and N.K.D. designed research; S.C.A. and T.A.B. conceived the project; S.C.A. performed research; S.C.A. analyzed data; and S.C.A., T.A.B., A.B.C., and N.K.D. wrote the paper.

The authors declare no conflict of interest.

This article is a PNAS Direct Submission.

Data deposition: The code and data reported in this paper have been deposited in GitHub, https://github.com/seananderson/heavy-tails, and Zenodo, http://doi.org/10.5281/zenodo.321930.

This article contains supporting information online at www.pnas.org/lookup/suppl/doi:10.1073/pnas.1611525114/-/DCSupplemental.

Freely available online through the PNAS open access option.

## References

- ↵.
- Sornette D

- ↵.
- Albeverio S,
- Jentsch V,
- Kantz H

- ↵.
- Taleb NN

- ↵
- ↵.
- Sornette D

- ↵
- ↵
- ↵
- ↵.
- Fey SB, et al.

- ↵
- ↵.
- Intergovernmental Panel on Climate Change

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵.
- Segura AM,
- Calliari D,
- Fort H,
- Lan BL

- ↵
- ↵
- ↵.
- Reed DH,
- O’Grady JJ,
- Ballou JD,
- Frankham R

- ↵.
- Saucy F

*Arvicola terrestris*. Oikos 71(3):381–392. - ↵.
- Stafford J

- ↵
- ↵.
- Gelman A, et al.

- ↵
- ↵
- ↵.
- Nuñez MA,
- Logares R

- ↵.
- Meehl GA,
- Tebaldi C

- ↵
- ↵
- ↵
- ↵
- ↵
- ↵.
- Ward EJ,
- Hilborn R,
- Towell RG,
- Gerber L

- ↵.
- NRC

- ↵
- ↵.
- Caddy JF,
- McGarvey R

- ↵
- ↵.
- Lindenmayer DB,
- Likens GE,
- Krebs CJ,
- Hobbs RJ

- ↵
- ↵
- ↵
- ↵.
- Hoffman MD,
- Gelman A

- ↵.
- Stan Development Team

*Stan Modeling Language User’s Guide and Reference Manual*Version 2.14.0. Available at mc-stan.org/documentation/. Accessed February 28, 2017. - ↵.
- Carpenter B, et al.

- ↵.
- R Core Team

- ↵.
- Inchausti P,
- Halley J

- ↵
- ↵
- ↵
- ↵
- ↵.
- Skaug H,
- Fournier D,
- Bolker B,
- Magnusson A,
- Nielsen A