EMPIRICAL PERSPECTIVES FROM MATHEMATICAL ECOLOGY Tail associations in ecological variables and their impact on extinction risk

. Extreme climatic events (ECEs) are becoming more frequent and more intense due to climate change. Furthermore, there is reason to believe ECEs may modify "tail associations" between distinct population vital rates, or between values of an environmental variable measured in different locations. "Tail associations" between two variables are associations that occur between values in the left or right tails of the distributions of the variables. Two positively associated variables can be principally "left-tail associated" (i.e., more correlated when they take low values than when they take high values) or "right-tail associated" (more correlated when they take high than low values), even with the same overall correlation coef ﬁ cient in both cases. We tested, in the context of non-spatial stage-structured matrix models, whether tail associations between stage-speci ﬁ c vital rates may in ﬂ uence extinction risk. We also tested whether the nature of spatial tail associations of environmental variables can in ﬂ uence metapopulation extinction risk. For instance, if low values of an environmental variable reduce the growth rates of local populations, one may expect that left-tail associations increase metapopulation extinction risks because then environmental "catastrophes" are spatially synchronized, presumably reducing the potential for rescue effects. For the non-spatial, stage-structured models we considered, left-tail associations between vital rates did accentuate extinction risk compared to right-tail associations, but the effect was small. In contrast, we showed that density dependence interacts with tail associations to in ﬂ uence metapopulation extinction risk substantially: For population models showing undercompensatory density dependence, left-tail associations in environmental variables often strongly accentuated and right-tail associations mitigated extinction risk, whereas the reverse was usually true for models showing overcompensatory density dependence. Tail associations and their asymmetries are taken into account in assessing risks in ﬁ nance and other ﬁ elds, but to our knowledge, our study is one of the ﬁ rst to consider how tail associations in ﬂ uence population extinction risk. Our modeling results provide an initial demonstration of a new mechanism in ﬂ uencing extinction risks and, in our view, should help motivate more comprehensive study of the mechanism and its importance for real populations in future work.

change. Furthermore, there is reason to believe ECEs may modify "tail associations" between distinct population vital rates, or between values of an environmental variable measured in different locations. "Tail associations" between two variables are associations that occur between values in the left or right tails of the distributions of the variables. Two positively associated variables can be principally "left-tail associated" (i.e., more correlated when they take low values than when they take high values) or "right-tail associated" (more correlated when they take high than low values), even with the same overall correlation coefficient in both cases. We tested, in the context of non-spatial stage-structured matrix models, whether tail associations between stage-specific vital rates may influence extinction risk. We also tested whether the nature of spatial tail associations of environmental variables can influence metapopulation extinction risk. For instance, if low values of an environmental variable reduce the growth rates of local populations, one may expect that left-tail associations increase metapopulation extinction risks because then environmental "catastrophes" are spatially synchronized, presumably reducing the potential for rescue effects. For the non-spatial, stage-structured models we considered, left-tail associations between vital rates did accentuate extinction risk compared to right-tail associations, but the effect was small. In contrast, we showed that density dependence interacts with tail associations to influence metapopulation extinction risk substantially: For population models showing undercompensatory density dependence, left-tail associations in environmental variables often strongly accentuated and right-tail associations mitigated extinction risk, whereas the reverse was usually true for models showing overcompensatory density dependence. Tail associations and their asymmetries are taken into account in assessing risks in finance and other fields, but to our knowledge, our study is one of the first to consider how tail associations influence population extinction risk. Our modeling results provide an initial demonstration of a new mechanism influencing extinction risks and, in our view, should help motivate more comprehensive study of the mechanism and its importance for real populations in future work.

INTRODUCTION
Over the past few decades, extreme climatic events (ECEs) such as temperature or precipitation extremes, heat waves, and extreme wildfires have become more common or more extreme due to anthropogenic climate change (Ummenhofer and Meehl 2017). Extreme climatic events are well known to produce severe impacts at the population, community, and ecosystem levels (Gutschick and BassiriRad 2003, Bragazza 2008, McKechnie and Wolf 2009, Jiguet et al. 2011, Felton and Smith 2017, and thus, the goal of improving our understanding of the variety of ecological impacts of ECEs is an important applied goal for modern ecology (Katz et al. 2005, Smith 2011, Diez et al. 2012, Moritz and Agudo 2013, Bailey and van de Pol 2016. Frequent extreme events can help alter the mean temperature or precipitation of an area, can be associated with changes through time in the overall variance of climate signals, and can directly contribute to changes in the symmetry or the skewness of distributions of climate signals through time (Hansen et al. 2012). All such changes have the potential to influence populations, including through risks of species extinction or extirpation from an area (Garcia-Carreras and Reuman 2013, Vasseur et al. 2014). For instance, Frederiksen et al. (2008) used 43 yr of ring-recovery data to examine the demographic impact of ECEs on European shags via their effects on temporal variability in juvenile, immature and adult survival rates. Long et al. (2017) investigated the sensitivity of 41 UK butterfly species to four different types of ECEs. At the community level, research has shown that ECEs can reorder the dominance hierarchy of species in an area, sometimes leading to selective local extinctions, and can also produce effects that cascade across trophic levels (Thomas et al. 2004, Edgar et al. 2010, Moreno and Møller 2011, Hoover et al. 2014. For example, Chiu et al. (2013) showed that extreme floods affected the annual survival of a prey community, in turn influencing predators in aquatic as well as adjacent riparian habitats. Dreesen et al. (2014), in an experimental study, showed that ECEs can reduce the ability of plant communities to withstand subsequent ECEs that occur after a short time interval.
Extreme climatic events also have impacts in ecological contexts for which spatial dynamics are important. In the metapopulation context, spatial synchrony of population dynamics, which can be caused by spatial synchrony of climatic fluctuations (Liebhold et al. 2004), is linked to metapopulation extinction risk (Hanski 1998, Earn et al. 2000. If ECEs are also spatially more extensive than moderate events, as seems likely for some types of events, extreme events may accentuate population synchrony and thereby increase extinction risk. For instance, a 21-yr-long study of Tack et al. (2015) on the Glanville fritillary butterfly demonstrated that frequency of drought during early larval instars increased the strength and extent of spatial synchrony over a metapopulation network of 4000 dry meadows, and could thereby influence the long-term viability of the species in the area of study.
Though earlier literature considered the influence of spatial synchrony in environmental variables on metapopulation extinction risk, to our knowledge the distinct influence of "asymmetric tail associations" between environmental time series in different locations has hardly been considered. "Tail associations" are associations between two random quantities in their extreme values. For instance, if two positively associated variables are more closely associated in their smaller values than in their larger values, they are said to have stronger "left-tail" or "lower-tail" association ( Fig. 1A, B). If two positively associated variables are more associated in their larger values than in their smaller values, they are said to have stronger "right-tail" or "upper-tail" association (Fig. 1C, D). A precise definition of a measure of tail association is given in Methods. Fig. 1 shows that two variables can have a wide range of patterns of tail association, from extreme leftto extreme right-tail association, with the same overall correlation. So tail association is a distinct concept from correlation itself and may have distinct ecological ramifications (Ghosh et al. 2020). Asymmetric tail associations have been documented between environmental variables, and between the same environmental variable measured in different places (Serinaldi 2008, Li et al. 2013, Goswami et al. 2018, She and Xia 2018. Asymmetric tail associations have also been documented between numerous ecological variables and between environmental and ecological variables (Ghosh et al. 2020).
Extinction risk is very commonly estimated for age-structured models such as stochastic matrix models (Caswell 2000, Morris andDoak 2002) and integral projection models (Merow et al. 2014). For such models, extinction risk depends, in part, on the long-term population growth rate, a population parameter which in turn depends on the means, variances, and covariances of population vital rates such as the survival and fecundity rates of age classes. Extinction risks may, a priori, also depend on possible asymmetries of tail association between population vital rates. For instance, it seems reasonable to imagine a case in which both adult fecundity and juvenile survival rates in a population depend on an environmental variable such as spring rainfall in similar ways, but only when that variable takes low values. When rainfall takes typical or high values, fecundity and survival may be limited instead by different, unrelated factors. This could produce left-tail dependence between these vital rates, possibly resulting in different extinction risk for the model compared to a case in which vital rates are symmetrically associated in their tails or right-tail associated but otherwise statistically similar (Fig. 2).
There are reasons to hypothesize that ECEs may produce asymmetric tail associations between climatic variables measured in different locations. The idea, here, is based on the reasonable assumption that ECEs are both more intense and spatially more extensive than non-extreme climatic events. For clarity, we explain the idea using the example case of extreme rain events. Suppose a single extreme rain event is associated with spatially widespread large precipitation values; that is, not only does it rain a lot in each locale, but it does so for many locales across a large area. Also suppose that smaller local precipitation measurements are associated, instead, with local, non-extreme weather which can differ from location to location, rather than with a single extreme rain event. The result could be stronger right-tail association, through time, between measurements of precipitation made in different locations: Large precipitation measurements are coherent across space, since they are all attributable to a single event, whereas small values are less spatially correlated because they come from a multiplicity of distinct local events.
There are also reasons to hypothesize that asymmetry of tail association in environmental variables may influence extinction risks for metapopulations. For instance, assuming for the sake of a simple example that low values of an environmental variable are "bad" for populations of a given species and high values are "good", greater left-tail associations between measurements of the variable in different habitat patches may cause higher metapopulation extinction risk because then very bad years for component populations would occur at the same time in many patches, reducing the potential for rescue effects. Greater right-tail associations would probably not have the same effect, even if overall correlations between  Two stochastic quantities can have the same Spearman correlation, but different degrees of tail association: (A) extreme left-tail association, (B) moderate left-tail association, (C) moderate right-tail association, and (D) extreme right-tail association. In all cases, Spearman's q was the same, q = 0.875, up to sampling variance (the sample correlation for each dataset is displayed on its panel). The extreme left-tail-(respectively, right-tail-) dependent case is perfectly correlated below (respectively, above) the threshold of 1/2. See Methods for information on how bivariate noise was generated with asymmetric tail associations. locations were the same, because in that case, it would be very good years for the component populations that would occur at the same time in many patches. Ghosh et al. (2020) substantiated this intuition by showing that it held true for a modeling setup using a spatial version of a very simple non-density-dependent population growth model (the Lewontin-Cohen model). In essence, if environmental "catastrophes" (i.e., extremely bad years for a population) are widely spatially synchronized it should create much greater metapopulation extinction risk than if "bonanzas" (i.e., extremely good years) are widely synchronized, even if overall (non-tail-specific) levels of environmental synchrony are the same in both these scenarios (Fig. 3). The present study extends tests of this idea to metapopulation models incorporating density dependence.
In this paper, we will explore the potential for asymmetric tail associations to influence extinction risk in populations and metapopulations. Specifically, we will address the following sets of questions. (Q1) Can the nature of tail associations between vital rates in an agestructured population model influence extinction risk for the model, and to what extent? (Q2) Does the nature of tail associations between measurements of an environmental variable made in different locations influence the extinction risk of a density-dependent metapopulation model that depends on the environmental variable? How does the nature of the density dependence mediate this influence? Will density dependence reverse or reproduce the earlier result (Ghosh et al. 2020), using non-density-dependent models, that lefttail dependence in environmental variables Year 1 Year 2 Year 1 Fig. 2. Left-tail association between vital rates makes it more likely that low values of fecundity and survival will happen simultaneously (A), causing both life stages to crash simultaneously. This seems likely to increase extinction risk. In contrast, right-tail association between vital rates instead makes it more likely that high values of fecundity and survival will happen simultaneously (B).
accentuates and right-tail dependence mitigates extinction risks? For Q1, we simulate a simple 2-stage model. We hypothesize that left-tail dependence between juvenile survival and adult fecundity rates increases extinction risk relative to right-tail dependence, even when the overall correlation is kept fixed. For Q2, we simulate spatial extensions of several simple and commonly used population models. We intentionally use a simulation approach and simple models for both Q1 and Q2 instead of attempting mathematical analysis and/or using models that are more realistic and complex, because our intention here is to provide an initial exploration only of the potential mechanisms outlined above. By so doing, we hope to open these potentially practically important questions to wider examination by other researchers, using a variety of models and perspectives. Though tail associations are used to assess risks in fields such as finance and hydrology (Alexander 2008, Chen et al. 2011, Li et al. 2013, Goswami et al. 2018, to our knowledge asymmetric tail associations have rarely been assessed for their importance for extinction risks of populations and species (but see Ghosh et al. 2020). Our essential goal is to advance the idea that tail associations may also be important in this context. Our questions are also an important step for understanding the potential influence of ECEs on extinction risks if ECEs are related to tail associations in environmental variables.

METHODS
To investigate Q1 from the Introduction, we considered a simple stochastic age-structured model with two age classes, where E(t) and A(t) are the numbers or population densities of eggs and adults, respectively, at time t. We chose a very simple model because our goal was to carry out an initial evaluation, only, of the realism of a potential mechanism. The stochastic vital rates s e (t) and f a (t) in Eq. 1 were selected, independently across times, as follows.
The egg survival rate s e (t) was drawn from a beta distribution with parameters a = 3 and b = 2.
Adult fecundity was f a ðtÞ ¼f a ðtÞ exp ðÀAðtÞ=XÞ, wheref a ðtÞ was drawn from a gamma distribution with shape and scale parameters k = 2 and h = 1, respectively, and Ω is a density dependence parameter. Allowing Ω to tend to ∞, exp (ÀA(t)/Ω) ? 1, so that, in that limit, the model did not include density dependence (it was then a standard stochastic matrix model). Finite, positive values of Ω made the model nonlinear, with density dependence in the adult fecundity rate. Stochastic matrix models are very commonly used in evolutionary ecology and quantitative conservation biology, so our model is an example of a very widely used category of models. See, for instance, the classic books of Caswell (2000) and Morris and Doak (2002).
The stochastic quantities s e (t) andf a ðtÞ were generated to have the same overall association (Spearman's q equal to 0.875, up to sampling variance), but different degrees and types of asymmetry of tail association in different simulations, as follows. For a simulation of length T, we used methods that can generate T independent pairs (u(t), v(t)) for t = 1, . . ., T with q = 0.875 (up to sampling variance), but with either extreme left-, moderate left-, moderate right-, or extreme right-tail association (Fig. 1). The methods we used generated (u(t), v(t)) such that u(t) was uniformly distributed on (0, 1), and likewise for v(t). We then applied a monotonic transformation to the u(t) so that the resulting numbers were distributed according to a beta distribution with parameters a = 3 and b = 2, and we used those values for the s e (t). We also applied a (different) monotonic transformation to the v(t) so that the resulting numbers were distributed according to a gamma distribution with parameters k = 2 and h = 1, and we used those values for thef a ðtÞ. Because the u(t) and v(t) were uniformly distributed, the inverse cumulative distribution functions of the desired beta and gamma distributions were the transformations. Because these were monotonic transformations, the Spearman correlation of s e (t) andf a ðtÞ was the same as that of u(t) and v(t). The quantities (u(t), v(t)) for t = 1, . . ., T were generated as described in Appendix S1: Section 1.
To calculate extinction risk for a given Ω and choice of tail dependence, we simulated the model for 25 time steps, 10,000 times, each simulation starting from E(t) = 5 and A(t) = 10. The population for a simulation was considered to have gone extinct if the total population E(t) + A(t) ever went below an extinction threshold (we used 2) within the 25-time-step time horizon. Extinction risk was the fraction of 10,000 simulations that went extinct. Confidence intervals (95%) of extinction risk were also computed based on a binomial distribution.
The variables u(t) and v(t) are visibly asymmetrically tail associated (Fig. 1), but we have instead used the transformed variables s e (t) and f a ðtÞ in our model. In what sense are s e (t) and f a ðtÞ asymmetrically tail associated? Given two boundary quantities l b and u b such that 0 ≤ l b < u b ≤ 1, and letting x(t) be the rank of u(t) in the set {u(t)|t = 1, . . ., T}, divided by T + 1, and letting y(t) be the rank of v(t) in the set {v(t)|t = 1, . . ., T}, divided by T + 1 (smallest elements are here defined to have rank 1), consider the portions of the lines x + y = 2l b and x + y = 2u b , which lie within the unit square (Fig. 4). We can then define the partial Spearman correlation, ( 2) where sample means and variances in this equation are computed using all the data, but the sum is computed only over the t such that The variance, var, here and throughout, is the sample variance, with Bessel's correction (denominator T À 1). The partial Spearman correlation is the component of the Spearman correlation u and v which can be attributed to the points for which the normalized ranks, x and y, lie between the boundary lines we defined in the unit square (Fig. 4). The quantity cor 0,0.5 (u, v)cor 0.5,1 (u, v) is a way to measure asymmetry of tail associations between u and v. And this quantity is the same as cor 0;0:5 ðs e ;f a Þ À cor 0:5;1 ðs e ;f a Þ, because the partial Spearman correlation, like the Spearman correlation itself, is based on ranks. Ranks are the same for u(t) and for s e (t) because these quantities are related by an increasing monotonic transformation; and likewise for v(t) andf a ðtÞ. Thus, in the sense captured by the statistic cor 0,0.5cor 0.5,1 , s e (t) andf a ðtÞ have the same asymmetry of tail association as the values u(t) and v(t) from which they were constructed. Rank-based approaches to understanding the nature of association between variables have been recommended, as such approaches are not influenced by the marginal distributions of the quantities (Genest and Favre 2007). Marginal distributions contain no information about association between variables, as each marginal pertains to only one of the variables. These ideas are revisited in the Discussion.
To investigate Q2 from the Introduction, we considered spatialized, stochastic versions of six density-dependent population models, the Ricker, Hassell, Maynard Smith, Pennycuick, Verhulst, and Austin-Brewer models. These models were all described by Cohen (1995). The models are also listed in Table 1 and are analyzed in detail in Appendix S1: Section 2, where original references for the models are also provided. All models took the general form Pðt þ 1Þ ¼ DkðtÞPðtÞ, where PðtÞ is a vector of population densities, the ith component P i (t) representing the population density in the ith habitat patch (i = 1, 2, . . ., n) at time t. The n 9 n matrix k(t) was diagonal with ith diagonal entry k i (t) depending on the model used (Table 1). The matrix D was an n 9 n dispersal matrix. For simplicity, we assumed the n habitat patches were arranged, evenly spaced, in a line. We considered both "local" and "global" dispersal among patches in different simulations. For local dispersal, a fraction d of each population dispersed during each time step, equally distributed to the two or one nearest neighbors. For global dispersal the same fraction d dispersed equally to the other n À 1 patches. When d was 0, dispersal did not occur and rescue effects were not possible.
We obtained our list of classic population models (which we then spatialized and made stochastic, as above) from Cohen (1995). Model parameters were selected so that varying the growth parameter, r, across a range, while keeping the other parameters fixed at selected values, caused the deterministic, one-patch versions of the models to transition from dynamics showing an undercompensatory (i.e., monotonic) approach to an equilibrium, to dynamics showing an overcompensatory (i.e., oscillatory) approach to an equilibrium. The transition value of r was denoted r c . Appendix S1: Section 2 summarizes stability analyses of the models. Whereas Cohen (1995) was interested in studying chaotic model behavior and therefore chose parameters for which the models exhibited chaotic dynamics, we were interested in non-chaotic behaviors, which seem likely to be more common in real populations. Two models that Cohen (1995) used were excluded from our analyses because their functional forms complicated the analyses we wished to perform. First, the equilibrium of the Malthus-Condorcet-Mill model (Cohen 1995) could not be determined analytically, making downstream analyses difficult for that model. Second, for the Varley model and for the portions of parameter space we explored, varying r produced no transition from under-to overcompensatory dynamics.
We simulated our metapopulation models using environmental noises e i (t) ( Table 1), independent through time, that were generated to have the same overall associations between all patches and in all simulations (Spearman's q again 0.875, up to sampling variance), but with different degrees and types of asymmetry of tail associations between patches in different simulations, as follows. For a simulation of length T, we used methods that can generate T independent n-tuples (u 1 (t), . . ., u n (t)) for t = 1, . . ., T with q = 0.875 between u i (t) and u j (t) for all i 6 ¼ j (up to sampling variance), but with either extreme left-, moderate left-, moderate right-, or extreme right-tail association between all pairs of locations (as in Fig. 1). Pairs u i (t) and u j (t) looked like Fig. 1 and were structured, statistically, in the same way as the stochastic quantities u(t) and v(t) used for the non-spatial simulations described above. For any i, u i (t) was again uniformly distributed on the interval (0, 1). To obtain the values e i (t), the u i (t) were transformed so that the resulting numbers were normally distributed with mean 0 and standard deviation r (Table 1), for all i. This was accomplished with the inverse cumulative distribution function of the normal distribution with mean 0 and standard deviation r. Because this was a monotonic transformation, the Spearman correlation between e i (t) and e j (t) was the same as that between u i (t) and u j (t); likewise partial Spearman correlations were the same. The quantities (u 1 (t), . . ., u n (t)) for t = 1, . . ., T were generated as described in Appendix S1: Section 1.
To calculate extinction risk for a given metapopulation model, given parameters, and a given choice of tail association, we simulated the model for 25 time steps, 10,000 times, using n = 5 habitat patches. The equilibrium value of the one-patch, deterministic version of the model, P e (Table 1) was used as the initial condition in all patches and simulations. After each time step, if P i (t) was less than P e /10, then it was set to 0. The whole metapopulation was considered to have gone extinct if all patches were 0 at time step 25. Extinction risk was the fraction of 10,000 metapopulations that went extinct. Estimates of extinction risk were plotted together with 95% confidence intervals based on a binomial distribution. Though these confidence intervals have limited utility because they can be made arbitrarily narrow by increasing the number of simulations, they are nonetheless presented to illustrate the confidence achieved with the number of simulations we performed.
To answer (Q2) of the Introduction, we estimated extinction risk for each metapopulation model, using each form of tail association for the environmental noise, using values of r that included under-and overcompensatory dynamics, as indicated above. We chose values of r within a range, the range specified as follows. The second part of Q2 (about how density dependence mediates the influence of asymmetric tail associations on extinction risk) was addressed by comparing results for under-and overcompensatory model regimes. Defining r min as the minimum value of r for which P e was non- Table 1. Summary of the six density-dependent stochastic metapopulation models we used.

Models
k i (t) Parameters used P e r min r c r bf Ricker exp(e i (t)) exp[r(1 -P i (t)/K)] K = 50, r = 1 K 0 1.00 2.00 Hassell exp(e i (t))r/(1 + aP i (t)) b a = 0.5, b = 100, r = 1 ( r 1/b À 1)/a 1 2.73 7.54 Maynard Smith Notes: The population growth rate in habitat patch i at time t, k i (t) (second column), is P i (t + 1)/P i (t), and defines how the model works within each patch. Here, P i (t) is the population in patch i at time t. The growth rate depends on model parameters, P i (t), and the environment e i (t) in patch i at time t. e i (t) was normally distributed with mean 0 and standard deviation r, and values were generated independently across times. P e is the equilibrium of the deterministic, one-patch version of the model. r min is the smallest r, given fixed values of the other parameters, for which P e was non-negative. r c is the value of r for which a transition occurs from under-to overcompensatory dynamics (see main text). r bf is the value of r for which a transition occurs from a stable to an unstable equilibrium (a bifurcation; see main text). For extinction-risk simulations, we used r such that r min < r < r bf . Dispersal between patches was also implemented (see main text). Formulas are given for P e , but specific values which pertain for the listed parameters are given for r min , r c , and r bf . See Appendix S1: Section 2 for the formulas that gave these values.
negative, and defining r bf as the value of r at which a bifurcation occurs, the range of r used spanned 95% of the interval (r min , r bf ). For details, see Appendix S1: Section 2. The range of r encompassed r c , so included a transition from under-to overcompensatory dynamics. For each value of r used within the range, we calculated extinction risk as indicated above. Dispersal parameter values d = 0, 0.1, . . ., 0.5 were also used in different simulations. We plotted extinction risk as a function of d for several fixed values of r, some producing under-and some overcompensatory dynamics. For d = 0, we also plotted extinction risk as a function of r. All code for this project is available at https://github.com/ sghosh89/ERC.

RESULTS
For our age-structured, non-spatial model, extinction risks were higher when associations between the vital rates s e and f a were principally for low values of these vital rates (left-tail association), and were lower when associations were principally for higher values, but the effect was minor for this model (Fig. 5). Though extinction risks were overall higher for the density-dependent version of our model compared to the density-independent version, extinction-risk differences across tail association scenarios were similar for the two versions. Because the same Spearman correlation (up to sampling variance) was used in all simulations, the differences between lines in the panels of Fig. 5 are attributable solely to the changes we made in the asymmetry of tail associations, and answer Q1. Differences were small, corresponding only to a few percentage points of additional extinction risk after 25 yr. A potential reason for this outcome is elaborated in the Discussion.
For our metapopulation simulations, for low values of r (close to r min ), left-tail associations of environmental noises accentuated extinction risks, and right-tail associations mitigated them; but the reverse was true for high values of r (close to r bf ). Thus, density-dependent models can either produce similar results, regarding the influence of tail associations on extinction risk, as the density-independent models of Ghosh et al. (2020), or opposite results, depending on the nature of density dependence. For each model of Table 1, using a value of r on the low end of the range (r min , r bf ), extinction risk was highest for extreme left-tail-associated noise, was lowest for extreme right-tail-associated noise, and was intermediate for the moderate tail association scenarios, for both local (Fig. 6, left panels) and global (Appendix S1: Fig. S1, left panels) dispersal, and for essentially all the dispersal rates, d, that we considered. However, using a value of r on the high end of the range (r min , r bf ), extinction risks were ordered oppositely (right panels of Fig. 6 and Appendix S1: Fig. S1). These results help address Q2 from the Introduction. Intermediate values of r are shown in Appendix S1: Figs. S2-S7 for local dispersal and Figs. S8-S13 for global dispersal among habitat patches.
Because the deterministic, one-patch versions of our models showed undercompensatory dynamics for r < r c (and therefore for r on the left end of the range [r min , r bf ]), and showed overcompensatory dynamics for r > r c , Fig. 6 and Appendix S1: Fig. S1 suggest the hypothesis that the opposite effects of tail association on extinction risk correspond to the domains in which models exhibit under-vs. overcompensatory dynamics. This was approximately, but not exactly, true. Careful inspection of Appendix S1: Figs. S6, S7 and S12, S13 and the r c values of Table 1 reveals that for the Verhulst and Austin-Brewer models, the switch occurs at a value of r less than r c . For the example case of d = 0, plots of extinction risk vs. r (Fig. 7) seem to indicate that for the Ricker, Hassell, and Pennycuick models, left-tail association accentuates and right-tail association mitigates extinction risk for undercompensatory dynamics (r < r c ), whereas the reverse occurs for overcompensatory dynamics (r > r c ). Fig. 7 reveals a similar reversal for the Maynard Smith, Verhulst, and Austin-Brewer models, but it occurs at a value of r close to but distinguishable from r c . Inspection of Appendix S1: Figs. S2-S13 shows that a qualitatively similar conclusion holds for nonzero d, for both local and global dispersal.

DISCUSSION
Our results suggest that asymmetry of tail associations between vital rates in a stage-structured, non-spatial stochastic model may be of limited importance for extinction risk. In our stage-structured, non-spatial model, extinction risk was only weakly influenced by asymmetry ❖ www.esajournals.org of tail associations between vital rates. Our results provide suggestive, rather than definitive evidence, here, because we considered only one model (but see below). In contrast, our metapopulation results show that asymmetries in tail associations between environmental variables in different locations can have a substantial influence on metapopulation extinction risk, for density-dependent as well as density-independent models. These results extend the result of Ghosh et al. (2020), which was for a density-independent model. For a density-independent model, Ghosh et al. (2020) corroborated the intuition that when "harmful" environmental events (i.e., events which reduce the growth rates of local populations) are correlated across space, extinction risk of metapopulations will be accentuated compared to when "beneficial" events (i.e., events which increase the growth of local populations) are spatially correlated, even if overall correlation is the same in both scenarios. This was presumably because the first of these Fig. 5. Extinction risk vs. time for our stochastic, age-structured model. Risks were higher when associations between vital rates were stronger in the left compared to the right tails of the rates, but only slightly so. Results are for the linear version of the model (Ω ? ∞; A, C) and for the nonlinear model with density dependence in the fecundity rate (Ω = 100; B, D). Results show cases of moderate (A, B) and extreme (C, D) asymmetry of tail association. The solid line is the estimated extinction risk, and the area of shading, which is barely visible because it is very narrow, shows 95% confidence intervals for extinction risk based on the number of simulations done. Abbreviations are LT, left tail; and RT, right tail.  , Table 1). Dispersal was local, see Appendix S1: Fig. S1 for global dispersal. Values of r used are in the panel headers. Extinction risks are for a 25-yr time horizon. Confidence intervals for extinction risk are not plotted, to provide a clearer plot. But because 10,000 simulations were done for each extinction-risk estimate, confidence intervals were very narrow, as in Fig. 5. Abbreviations are LT, left tail; and RT, right tail.   scenarios reduces the potential for rescue effects. For the density-dependent metapopulation models we considered, the same result held, typically, but not always, for undercompensatory model dynamics. However, for overcompensatory dynamics, the opposite result held, again typically, but not always: When "harmful" events were correlated across space, extinction risk of metapopulations was actually mitigated compared to when "beneficial" environmental events were spatially correlated by the same amount. This was probably because the correlated beneficial events caused populations to exceed carrying capacity, and to subsequently crash, across all or most habitat patches.
Our exploration for stage-structured models suggested limited influence of tail associations between vital rates on extinction risk, but some additional lines of inquiry might be worth exploring as part of an effort to understand how general our results are. First, we considered a model with positive association between egg and adult vital rates, corresponding to similar impacts of an environmental variable on eggs and adults. Negative associations between vital rates can occur as well, and the impacts of tail associations should be considered in that context as well. Negative associations relate to life-history trade-offs, for example, between survival and reproduction (Lande et al. 2003, p. 59). Asymmetric tail associations may occur in such a context if, for instance, there is a three-way trade-off among vital rates r 1 , r 2 , and r 3 . Values of r 1 near an upper limit may be strongly associated with r 2 , whereas lower values of r 1 may not be strongly associated with r 2 because r 1 trades off against both r 2 and r 3 . It is possible that tail associations in negatively related vital rates produce different outcomes for extinction risk from what we observed by studying positively associated vital rates.
A second topic of research that our work suggests has to do with Tuljapurkar's well-known approximation of the long-term stochastic growth rate (Tuljapurkar 1982(Tuljapurkar , 1990. It may be possible to obtain a better understanding, using Tuljapurkar's approximation, of our result that asymmetries of tail association may have limited influence on the extinction risk of stage-structured, non-spatial models. The long-term stochastic growth rate is an important quantity in demography (Caswell 2000) and population viability analysis for conservation (Morris and Doak 2002), and it factors prominently into extinction-risk estimates for density-independent age-or stage-structured population models (Lande and Orzack 1988), Pðt þ 1Þ ¼ MðtÞPðtÞ. Here, PðtÞ is the age-or stage-structured population, and M is a matrix-valued stationary stochastic process. The matrix entries of M can correspond, for instance, to the fluctuating fecundity and survival rates of the life stages of the population. For such a model, Tuljapurkar's approximation expresses the long-term stochastic growth rate as a function of the variances and covariances of the entries of M and of the eigenvalues and eigenvectors of the expected value of M. None of these quantities depends on tail associations between matrix entries. Thus, whenever Tuljapurkar's approximation is a good approximation, the longterm stochastic growth rate should not depend sensitively on tail association. Lande et al. (2003) provide reasoning supporting a claim that Tuljapurkar's approximation is "likely to be accurate in many cases" (p. 60), though we are unaware of the extent to which accuracy of the approximation has been tested for scenarios with asymmetric tail associations. It seems potentially worth exploring the accuracy of Tuljapurkar's approximation under scenarios of tail association in the vital rates, and then relating results to the effects of tail associations on the long-term stochastic growth rate and extinction risks. We note that the above reasoning pertains only to density-independent models and says nothing about our results for the density-dependent version of our stage-structured model.
Extensive evidence now exists that ECEs have become more frequent and more intense over the last several decades, and substantial progress has recently been made so that some event categories and even specific events can now be attributed to human-induced climate change (Ummenhofer and Meehl 2017). We now provide an explanation of how such changes may modify the asymmetry of tail association of weather variables measured in different places. Such modification may thereby influence metapopulation extinction risk through the mechanisms explored earlier in this paper. Let X i (t) and X j (t) denote the values of some weather variable as measured in fixed locations i and j at time t. Although changes in the correlation, through time, between X i (t) and X j (t) are increasingly documented, and may also be a consequence of climate change (Black et al. 2018), we are unaware of any studies examining the extent to which tail associations between X i (t) and X j (t), and asymmetries of tail association, may be changing. If ECEs are becoming more spatially extensive, as seems likely given that they are becoming more intense (Ummenhofer and Meehl 2017), then the conditional probability that both of X i (t) and X j (t) are influenced by an ECE, given that one of them is so influenced, should be increasing. That is, since ECEs are spatially bigger, it is now more likely that a single weather event influences both locations i and j given that it influences one of them. This should contribute to increased association of the two variables in one of their tails (which tail depends on the type of ECE and what the variable X is). A study of changes through time in tail associations of weather variables would be straightforward and potentially important given the mechanisms we have explored.
It is known that the distributions X i (t) are changing, in mean, variance, and possibly other ways, depending on what environmental variable X represents. Studying associations between X i (t) and X j (t) using standard methods such as Pearson correlation may be complicated by these changes; rankbased measures of association such as we have used may be more appropriate for some applications. If f is a strictly monotonic, increasing transformation, then the Spearman (or Kendall) correlation of f(X i (t)) and f(X j (t)) is the same as the Spearman (or Kendall) correlation of X i (t) and X j (t). However, the Pearson correlation of f(X i (t)) and f(X j (t)) may differ substantially from the Pearson correlation of X i (t) and X j (t). Relatedly, if the marginal distributions of X i (t) and X j (t) are changing through time, it may also cause changes through time in the Pearson correlation of these quantities, even while the values of rank-based correlations such as Spearman and Kendall may be unchanging. Essentially, in such a case, associations between X i (t) and X j (t) would be unchanging through time even though Pearson correlation appears to indicate a change, which is instead a change in marginal distributions. Perhaps unfortunately, the Pearson correlation, which is the most commonly used measure of association of two variables, depends not only on association information, but also on the marginals. For this reason, it is possible that recently observed changes in environmental (Black et al. 2018) and biological (Koenig and Liebhold 2016) synchrony are, at least in part, just another feature of alreadystudied changes in distributional characteristics of variables, rather than reflecting true changes in association as would be measured using rankbased correlation methods.
Rank-based methods are linked to a well-developed suite of statistical approaches related to the copula concept of statistics. Copulas and related ideas can be used to formally separate the information in the joint distribution of X i and X j into the information that pertains to the marginal distributions of these variables, and the information that pertains solely to their association. Ghosh et al. (2020) argued generally for the potential of copulas to improve our understanding of ecology (see also Valpine et al. 2014, Anderson et al. 2018, Popovic et al. 2019. Rank-based methods such as the Spearman and partial Spearman correlations we used are closely related to copulas, and Ghosh et al. (2020), Genest and Favre (2007), and others have recommended copula and rank approaches when the goal is to understand associations between variables. We recognize that the Pearson correlation of X i (t) and X j (t), and changes through time in that correlation, will still be useful for many applications. But to understand potential changes in the association of X i (t) and X j (t), unconfounded by changes in their marginal distributions, it is necessary to use rank-based or copula methods (Genest andFavre 2007, Ghosh et al. 2020).
Finally, our results suggest that combining ideas about tail association with recently developed ideas about the geography of synchrony may provide fertile territory for future research. Our metapopulation models assumed that habitat patches were equispaced on a transect or were all equally mutually accessible through dispersal. Furthermore, associations between environmental variables in different patches were assumed to be spatially homogeneous. In real metapopulations, neither of these assumptions holds. Several recent studies, reviewed and synthesized by Walter et al. (2017), have explored geographic variation in the spatial synchrony of environmental and population variables. Walter et al. (2017) coined such variation the "geography of synchrony" and explored its mechanisms, consequences, and ecological importance. For instance, synchrony in several systems was found to have prominent geographic variation, which can help provide inferences about the causes of synchrony and about organism ecology. Two among many research questions that could be explored that blend the geography of synchrony with tail associations include the following: (1) What is the nature of the geography of tail associations of environmental and ecological variables? and (2) Do asymmetries of tail associations of environmental variables interact with geographies of synchrony to influence metapopulation extinction risk? Connecting the ideas of tail association that we have developed with the ideas on the geography of synchrony of Walter et al. (2017) seems like a promising future direction.