MainText.Rmd

---
title: "Tail associations in ecological variables and their impact on extinction risk"
author:
 - Shyamolina Ghosh\textsuperscript{1}, Lawrence W. Sheppard\textsuperscript{1}, Daniel C. Reuman\textsuperscript{1,2,\dag}
 - \textsuperscript{1}Department of Ecology and Evolutionary Biology and Kansas Biological Survey, University of Kansas, Lawrence, KS, 66045, USA; \textsuperscript{2}Laboratory of Populations, Rockefeller University, 1230 York Ave., New York, NY 10065, USA
date: '\textsuperscript{\dag}Correspondence: Daniel Reuman, 2101 Constant Ave, Lawrence, KS, 66047, reuman@ku.edu, 626 560 7084'
fontsize: 12pt
geometry: "left=1in,right=1in,top=1in,bottom=1in"
output: 
  pdf_document:
    number_sections: no
    keep_tex: yes
    fig_caption: yes
    includes:
      in_header: head_maintext.sty

mainfont: Times New Roman      
tables: True
link-citations: True
urlcolor : blue
indent : True

csl: ecosphere.csl
bibliography: Ref_extrisk.bib
---


```{r setup_MainText, echo=F}
library(rmarkdown)
knitr::opts_chunk$set(echo = TRUE, fig.pos = "H")
options(scipen = 1, digits = 3) #This option round all numbers appeared in the inline r code upto 5th digit
seed<-101
```

<!--\noindent \emph{Affiliations:}

\noindent Ghosh: Department of Ecology and Evolutionary Biology and Kansas Biological Survey, University of Kansas, Lawrence, KS, 66045, USA

\noindent Sheppard: Department of Ecology and Evolutionary Biology and Kansas Biological Survey, University of Kansas, Lawrence, KS, 66045, USA

\noindent Reuman: Department of Ecology and Evolutionary Biology and Kansas Biological Survey, University of Kansas, Lawrence, KS, 66045, USA; Laboratory of Populations, Rockefeller University, 1230 York Ave, New York, NY, 10065, USA

\noindent \emph{Correspondence:} Daniel Reuman, 2101 Constant Ave, Lawrence, KS, 66047, reuman@ku.edu, 626 560 7084. 
-->

\newpage

\noindent \textbf{ABSTRACT} 
<!--350 words limit : Ecosphere-->

\noindent Extreme climatic events are becoming more frequent and more intense
due to climate change. Furthermore, there is reason to believe extreme 
climatic events 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 modelling 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.

\vspace{12pt}
\noindent \textbf{\textit{Keywords:}} copula; density-dependent model; extinction risk; extreme events; mathematical ecology; tail association. 

\newpage

```{r read_res, echo=F}
srho_common<-readRDS("./Results/srho_common.RDS")
```

# Introduction\label{Introduction}

<!--ECEs have a lot of impacts, it is important to understand their impacts-->
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 
[@ummenhofer2017extreme]. ECEs are well known to produce severe
impacts at the population, community and ecosystem levels [@bragazza2008climatic; 
@gutschick2003extreme; @jiguet2011community; @mckechnie2009climate;
@felton2017integrating], and thus the goal of improving our understanding of
the variety of ecological impacts of ECEs is an important applied goal for 
modern ecology [@katz2005statistics; @smith2011ecological; @diez2012will; @moritz2013future; @bailey2016tackling]. 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 
[@hansen2012perception]. All such changes have the potential to influence 
populations, including through risks of species extinction or extirpation from an area 
[@GarciaCarerras2013; @Vasseuretal2014].
For instance, @frederiksen2008demographic used $43$ years 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.
@mcdermott2017sensitivity 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 [@thomas2004extinction; @edgar2010nino; @moreno2011extreme; @hoover2014resistance]. 
For example, @chiu2013impact showed that extreme floods affected 
the annual survival of a prey community, in turn influencing predators in 
aquatic as well as adjacent riparian habitats. 
@dreesen2014successive, 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.

<!--Earlier study on metapopulation extinction risk : Large-scale spatial synchrony 
may also be generated by other processes such as migration and predation, but regional
stochasticity in the form of spatially correlated weather conditions is
probably the dominant synchronizing mechanism (ref: Nature1998) but what about taildep. effect on extrisk?.-->

<!--ECEs have impacts in the spatial context as well-->
ECEs 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 [@Liebhold2004], is linked to metapopulation
extinction risk [@hanski1998metapopulation; @earn2000coherence].
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-year-long
study of @tack2015increasing 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.

<!--concept for correlation in the extreme tail: tail association-->
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 (Figure \ref{fig_intro_cop}A, 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 
(Figure \ref{fig_intro_cop}C, D). A precise definition of a measure of tail 
association is given in Methods. Figure \ref{fig_intro_cop} shows 
that two variables can have a wide 
range of patterns of tail association, from extreme left- to 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_copula]. Asymmetric tail associations have been documented 
between environmental variables, and between the same environmental variable 
measured in different places [@Serinaldi2008; @li2013; @She2018; @Goswami2018]. 
Asymmetric tail associations have also been documented between numerous ecological
variables and between environmental and ecological variables [@Ghosh_copula].

Extinction risk is very commonly estimated for age-structured models 
such as stochastic matrix models [@MorrisDoak2002; @Caswell2000] and 
integral projection models
[@Merowetal2014]. 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. \ref{VitalRatePedagogFig}). 

<!--ECEs probably produce tail associations-->
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, i.e., 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 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_copula 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. \ref{MetapopPedagogFig}). The present study extends tests 
of this idea to metapopulation models incorporating density dependence. 

<!--objectives, last para of intro-->
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 age-structured 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_copula], using non-density-dependent models, that 
left-tail dependence in environmental variables 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 
[@alexander2008market; @li2013; @chen2011flood; @Goswami2018], to our
knowledge asymmetric tail associations have rarely been assessed for their
importance for extinction risks of populations and species (but see @Ghosh_copula).
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. 

<!--
How does environmental noises with asymmetric tail association across multiple patches
influence extinction risk of a density-dependent metapopulation model? Our speculation 
is that the general finding from density-independent metapopulation model [@Ghosh_copula]
(higher extinction risk if noises across patches were left-tail dependent than the one 
if noises had right-tail association) would hold true only in the under-compensatory regime, but 
should be reverse if the model shows over-compensatory dynamics. 
(Q3) We further hypothesize that when there is no dispersal, the above-said switching should occur
near the critical parameter value where the deterministic model shows transition from 
under-compensatory to over-compensatory dynamics.-->

<!--FIG : pedagogical fig showing 4 types of tail associations: 
extreme left, moderate left, moderate right, extreme right-->
\begin{figure*}[!h]
\begin{center}
\includegraphics[width=15 cm]{./Results/introfig_copula.pdf}
\caption{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 $\rho$ was the same, $\rho =$ `r srho_common`,
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.
\label{fig_intro_cop}}
\end{center}
\end{figure*}

\begin{figure}
\begin{center}
\includegraphics[width=25cm]{vitalrate_pedagog_fig.pdf}
\caption{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). \label{VitalRatePedagogFig}}
\end{center}
\end{figure}

\begin{figure}
\begin{center}
\includegraphics[width=25cm]{metapop_pedagog_fig.pdf}
\caption{Under extreme left-tail association of environments between two locations, catastrophes happen 
simultaneously in both locations (A). Whereas under extreme right-tail association, catastrophes do not happen
simultaneously, even when overall environmental association is the same in both scenarios (B). \label{MetapopPedagogFig}}
\end{center}
\end{figure}

# Methods\label{Methods}

<!--stochastic matrix model : age-structured-->
To investigate Q1 from the Introduction,
we considered a simple stochastic age-structured model
with two age classes, 
\begin{equation}
\left({\begin{array}{cc} E(t+1) \\ A(t+1) \end{array}}\right) = 
\left(\begin{array}{cc} 0 & f_a(t)\\ s_e(t) & 0 \end{array}\right)
\left(\begin{array}{c} E(t) \\ A(t) \end{array}\right), \label{eq_age_str}
\end{equation}
\noindent 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. \ref{eq_age_str} were selected, independently 
across times, as follows. The egg survival rate $s_e(t)$ was 
drawn from a beta distribution with parameters $\alpha=3$ and $\beta=2$. 
Adult fecundity was $f_a(t)=\tilde{f}_a(t)\exp(-A(t)/\Omega)$, where $\tilde{f}_a(t)$
was drawn from a gamma distribution with shape and scale parameters $k=2$
and $\theta=1$, respectively, 
and $\Omega$ is a density dependence parameter. Allowing $\Omega$ to tend to $\infty$,
$\exp(-A(t)/\Omega) \rightarrow 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 
$\Omega$ 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 @Caswell2000, and @MorrisDoak2002.

The stochastic quantities $s_e(t)$ and $\tilde{f}_a(t)$ were generated 
to have the same overall association (Spearman's $\rho$ 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,\ldots,T$ with
$\rho=0.875$ (up to sampling variance), but with either extreme left-,
moderate left-, moderate right-, or extreme right-tail association 
(e.g., Fig. \ref{fig_intro_cop}). 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
$\alpha=3$ and $\beta=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 $\theta=1$, and we used those values for the $\tilde{f}_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)$ and $\tilde{f}_a(t)$ was the same as that of 
$u(t)$ and $v(t)$. The quantities $(u(t),v(t))$ for $t=1,\ldots,T$ were
generated as described in Appendix, section \ref{SM-makenoise}.

<!--calculation of ext risk with age-structured model-->
To calculate extinction risk for a given $\Omega$ and choice of tail dependence, 
we simulated the model for $25$ time steps, $10000$ 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 $10000$
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. \ref{fig_intro_cop}),
but we have instead used the transformed variables $s_e(t)$ and $\tilde{f}_a(t)$
in our model. In what sense are $s_e(t)$ and $\tilde{f}_a(t)$ asymmetrically
tail associated? Given two boundary quantities $l_b$ and $u_b$ such
that $0 \leq l_b < u_b \leq 1$, and letting $x(t)$ be the rank of 
$u(t)$ in the set $\{u(t) | t=1,\ldots,T\}$, divided by $T+1$, and letting 
$y(t)$ be the rank of 
$v(t)$ in the set $\{v(t) | t=1,\ldots,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. \ref{PartialSpearmanFig}). 
We can then define the *partial Spearman correlation*,
\begin{equation}
\cor_{l_b,u_b}(u,v) = \frac{\sum (x(t)-\mean(x))(y(t)-\mean(y))}{(T-1)\sqrt{\var(x) \var(y)}},\label{PartialSpearmanEq}
\end{equation}
\noindent 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 $x(t)+y(t) > 2l_b$ and $x(t)+y(t)<2u_b$. 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 of $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. \ref{PartialSpearmanFig}). 
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,\tilde{f}_a)-\cor_{0.5,1}(s_e,\tilde{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)$ and $\tilde{f}_a(t)$. Thus, in the sense captured
by the statistic $\cor_{0,0.5}-\cor_{0.5,1}$, $s_e(t)$
and $\tilde{f}_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 (e.g., @Genest2007). 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.

<!--FIG : schematic diagram for partial spearman correlation, adapted from BIVAN-->
\begin{figure*}[!h]
\begin{center}
\includegraphics[width=6 cm]{./Results/PartialSpearmanFig.pdf}
\caption{The partial Spearman correlation, $\cor_{l_b,u_b}(u, v)$, within a band can be computed
for any band to describe how the strength of association between $u$ and $v$
varies from one part of the two distributions to another. Diagonal lines 
show two bands, the data in the right/upper band showing
stronger association than those in the left/lower band. Each band is bounded by
the lines $x+y=2l_b$ and $x+y=2u_b$ described in Methods (for values of $l_b$ and $u_b$
that depend on the band to be used). This figure was reproduced
largely following Fig. 7 of Ghosh \emph{et al.} (2020).
\label{PartialSpearmanFig}}
\end{center}
\end{figure*}

<!--density-dependence metapopulation model-->
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 @cohen1995unexpected. 
The models are also listed in Table \ref{tab_model_summary}, and 
are analyzed in detail in Appendix, section \ref{SM-model_stability}, where 
original references for the models are also provided.
All models took the general form $\vec{P} (t+1) = D \lambda(t) \vec{P}(t)$,
where $\vec{P}(t)$ is a vector of population densities, 
the $i^{\text{th}}$ component $P_i(t)$ representing
the population density in the $i^{\text{th}}$ habitat patch ($i=1,2,\ldots,n$) at time $t$. 
The $n \times n$ matrix $\lambda(t)$ was diagonal with $i^{\text{th}}$ diagonal entry 
$\lambda_i(t)$ depending on the model used (Table \ref{tab_model_summary}). The matrix $D$
was an $n \times 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 @cohen1995unexpected. 
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, section \ref{SM-model_stability}
summarizes stability analyses of the models.
Whereas @cohen1995unexpected
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 @cohen1995unexpected 
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 [@cohen1995unexpected] 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 $\varepsilon_i(t)$ (Table \ref{tab_model_summary}), 
independent through time, that were generated to have the same overall associations 
between all patches and in all simulations (Spearman's $\rho$ again `r srho_common`, 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),\ldots,u_n(t))$
for $t=1,\ldots,T$ with $\rho=$ `r srho_common` between $u_i(t)$ and $u_j(t)$ for all $i \neq 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. \ref{fig_intro_cop}). 
Pairs $u_i(t)$ and $u_j(t)$ looked like Fig. \ref{fig_intro_cop} 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 $\varepsilon_i(t)$, the $u_i(t)$ were 
transformed so that the resulting numbers were
normally distributed with mean $0$ and standard deviation $\sigma$ (Table \ref{tab_model_summary}), 
for all $i$. This was 
accomplished with the inverse cumulative distribution function of the normal distribution
with mean $0$ and standard deviation $\sigma$. Because this was a monotonic transformation,
the Spearman correlation between $\varepsilon_i(t)$ and $\varepsilon_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),\ldots,u_n(t))$
for $t=1,\ldots,T$ were generated as described in Appendix, section \ref{SM-makenoise}.

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, $10000$ times, using $n=5$ habitat patches.
The equilibrium value of the 
one-patch, deterministic version of the model, $P_e$ (Table \ref{tab_model_summary})
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 10000 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.

<!--test hypothesis from Q2 : ext risk with density-dependence metapopulation model with r from different regime-->
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_\text{min}$ as the minimum value of $r$ 
for which $P_e$ was non-negative; and defining
$r_\text{bf}$ as the value of $r$ at which a bifurcation occurs, the range of
$r$ used spanned 95\% of the interval ($r_\text{min},r_\text{bf}$). For details, see Appendix, section 
\ref{SM-model_stability}. 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.2,\ldots,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\label{Results}

<!--Q1 : age-structured model 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. \ref{fig_extrisk_agestr}).
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. \ref{fig_extrisk_agestr} 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 years. A 
potential reason for this outcome is elaborated in the Discussion.

<!--Q2 : metapopulation results-->
For our metapopulation simulations, for low values of $r$ (close to $r_\text{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_\text{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_copula, or
opposite results, depending on the nature of density dependence. 
For each model of Table \ref{tab_model_summary}, using a value of $r$ on
the low end of the range $(r_\text{min},r_\text{bf})$, extinction risk was highest
for extreme left-tail associated noise, lowest for extreme right-tail associated
noise, and was intermediate for the moderate tail association scenarios, for
both local (Fig. \ref{fig_extrisk_allmodel_localdisp}, left panels) and global 
(Fig. \ref{SM-fig_extrisk_allmodel_globaldisp}, 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_\text{min},r_\text{bf})$,
extinction risks were ordered oppositely (right panels of Figs 
\ref{fig_extrisk_allmodel_localdisp} and \ref{SM-fig_extrisk_allmodel_globaldisp}). 
These results help address Q2 from the Introduction.
Intermediate values of $r$ are shown in Figs. \ref{SM-fig_ricker_risk_local_disp_scl_1} - 
\ref{SM-fig_risk_local_disp_abrewer_r_vary_scl_0.6} for local dispersal and 
Figs. \ref{SM-fig_ricker_risk_global_disp_scl_1} - 
\ref{SM-fig_risk_global_disp_abrewer_r_vary_scl_0.6} 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_\text{min},r_\text{bf})$),
and showed overcompensatory dynamics for $r>r_c$, Figs \ref{fig_extrisk_allmodel_localdisp}
and \ref{SM-fig_extrisk_allmodel_globaldisp} suggest the hypothesis that the
opposite effects of tail association on extinction risk correspond to the domains
in which models exhibit under- versus overcompensatory dynamics. This was 
approximately, but not exactly, true. Careful inspection of Figs 
\ref{SM-fig_risk_local_disp_verhulst_r_vary_scl_0.6} - \ref{SM-fig_risk_local_disp_abrewer_r_vary_scl_0.6} and Figs \ref{SM-fig_risk_global_disp_verhulst_r_vary_scl_0.6} - \ref{SM-fig_risk_global_disp_abrewer_r_vary_scl_0.6} 
and the $r_c$ values of Table \ref{tab_model_summary} 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 versus $r$ (Fig. \ref{fig_extrisk_d0}) 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. \ref{fig_extrisk_d0}
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 Figs
\ref{SM-fig_ricker_risk_local_disp_scl_1}-\ref{SM-fig_risk_global_disp_abrewer_r_vary_scl_0.6}
shows that a qualitatively similar conclusion holds for nonzero $d$, for both local and global
dispersal.

# Discussion\label{Discussion}

<!--summary results:-->
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 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 for density-independent 
models. These results
extend the result of @Ghosh_copula, which was for a density-independent model. For a
density-independent model, @Ghosh_copula 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 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. 

<!--future ideas start here-->
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, e.g. between survival and reproduction 
(e.g., @Lande2003, 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.

<!--Tuljapurkar approximation-->
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 [@tulja1982; @tulja1990].
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
[@Caswell2000] and population viability analysis for conservation 
[@MorrisDoak2002], and it factors prominently into
extinction risk estimates for density-independent age- or 
stage-structured population models 
[@Lande1988], $\vec{P}(t+1)=M(t)\vec{P}(t)$. Here $\vec{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 long-term stochastic growth rate should not depend sensitively
on tail association. @Lande2003 provides 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. 

<!--linking back to ECEs - this para about the metapop case-->
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 [@ummenhofer2017extreme]. 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 (e.g., @Black2018),
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 [@ummenhofer2017extreme], 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. 
<!--***DAN: Shyamolina, another aspect of an exploration you may do on the 
the geography of tail association may be to look at how tail association 
may be changing. This relates to ECEs, so may help better justify the work
to a funder. You had proposed to look at geography of tail association in 
NDVI, and that makes good sense. You could also look at how tail association
is changing in NDVI (the records are about 35 years long, so you could just
about do this). Even better (because data are better), you could get weather
station records (lots of these can be downloaded, over 1000 stations
have been operating for over a century continuously in North America) and 
look at how tail associations are changing in those data.-->

<!--copulas and rank methods-->
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; rank-based 
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
(e.g., @Black2018) and biological (e.g., @Koenig2016) synchrony 
are, at least in part, just another feature of 
already-studied changes in distributional characteristics of 
variables, rather than reflecting true changes in association as would be measured
using rank-based 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_copula argued generally for the potential of copulas 
to improve our understanding of ecology 
(see also @Popovic2019; @Valpine2014; @Anderson2018).
Rank-based methods such as the Spearman and partial Spearman correlations
we used are closely related to copulas, and @Ghosh_copula,
@Genest2007 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 [@Ghosh_copula;
@Genest2007].

<!--relate things to the geography of synchrony-->
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 @walteretal2017, have explored geographic variation in the spatial synchrony of environmental
and population variables. @walteretal2017 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: 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 @walteretal2017 seems like a promising future direction.

<!--future study: worthy for PVA, IUCN red list:
1. Identify which life stage or demographic process as management target?
2. Old study (Doak1994, Ferguson1995) showed including correlation btw vital rates for PVA makes a difference, 
   what if we include TD through copula? Multisite age str. model with dispersal, directional change
   in spatio-temporal variation for vital rates.
3. Tdep. with -ve correlation between vital rates?
4. extinction risk how propagated through a food web?
5. Not sure : see McCarthy2000: Most previous models of metapopulation viability have either ignored
correlated extinction of patches, or have included uniform correlations, with the correlation between 
nearby patches being the same as the correlation between distant
patches (Hanski 1994; Tilman et al. 1994; Lindenmayer et al. 1995; McCarthy
1996). 
-->

# Acknowledgments

\noindent The authors were partly supported by US National Science Foundation
grant 1714195 and the James S. McDonnell Foundation, and the California Department of Fish 
and Wildlife Delta Science Program. 

 <!--FIG : age structure model results : linear, nonlinear-->
\begin{figure*}[!h]
\textbf{ \hspace{5.4 cm} (A) \hspace{4 cm} (B)} \\
\vspace{-0.9 cm}
\begin{center}
\includegraphics[width=5 cm]{./Results/age_str_results/extrisk_age_str_C_SC_linear.pdf}
\includegraphics[width=5 cm]{./Results/age_str_results/extrisk_age_str_C_SC_nonlinear_omg_100.pdf}\\
\end{center}
\textbf{ \hspace{5.4 cm} (C) \hspace{4 cm} (D)} \\
\vspace{-0.9 cm}
\begin{center}
\includegraphics[width=5 cm]{./Results/age_str_results/extrisk_age_str_extremetail_linear.pdf}
\includegraphics[width=5 cm]{./Results/age_str_results/extrisk_age_str_extremetail_nonlinear_omg_100.pdf}
\caption{Extinction risk versus 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 ($\Omega \rightarrow \infty$; A, C) and for the nonlinear model with density
dependence in the fecundity rate ($\Omega = 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. LT=left tail, RT=right tail.\label{fig_extrisk_agestr}}
\end{center}
\end{figure*}

<!--FIG :  summary : all nonlinear model results for local dispersion : extrisk vs. d-->
\begin{figure*}[!h]
\begin{center}
\includegraphics[width=14 cm]{./Results/allmodel_ext_risk_local_disp_extreme_and_moderate_noise.pdf}
\caption{Extinction risk versus dispersal rate for different scenarios of tail association
for our six metapopulation models (Methods, Table \ref{tab_model_summary}). Dispersal was
local, see Fig. \ref{SM-fig_extrisk_allmodel_globaldisp} for global dispersal. Values of $r$
used are in the panel headers.
Extinction risks are for a 25-year time horizon. Confidence intervals for extinction 
risk are not plotted, to provide a clearer plot. But 
because 10000 simulations were
done for each extinction risk estimate, confidence intervals
were very narrow, as in Fig. \ref{fig_extrisk_agestr}. LT=left tail, RT=right 
tail. \label{fig_extrisk_allmodel_localdisp}}
\end{center}
\end{figure*}


<!--FIG : all nonlinear model results, d=0, extrisk vs. r-->
\begin{figure*}[!h]
\begin{center}
\includegraphics[width=12 cm]{./Results/common_legend_extreme_and_moderate_noise.pdf}\\
\includegraphics[width=6 cm]{./Results/ricker_scl_1/ricker_extrisk_d0_extreme_and_moderate_noise.pdf}
\includegraphics[width=6 cm]{./Results/hassell_scl_1/hassell_extrisk_d0_extreme_and_moderate_noise.pdf}\\
\vspace{-0.6 cm}
\includegraphics[width=6 cm]{./Results/msmith_scl_0.8/msmith_extrisk_d0_extreme_and_moderate_noise.pdf}
\includegraphics[width=6 cm]{./Results/pennycuick_scl_1/pennycuick_extrisk_d0_extreme_and_moderate_noise.pdf}\\
\vspace{-0.6 cm}
\includegraphics[width=6 cm]{./Results/verhulst_scl_0.6/verhulst_extrisk_d0_extreme_and_moderate_noise.pdf}
\includegraphics[width=6 cm]{./Results/austinbrewer_scl_0.6/abrewer_extrisk_d0_extreme_and_moderate_noise.pdf}\\
\caption{Extinction risk versus $r$, with $d=0$, for different scenarios of tail association for 
our six metapopulation models (Methods, Table \ref{tab_model_summary}). 
Extinction risks are for a 25-year time horizon. Vertical lines show $r_c$,
the transition from under- to overcompensatory dynamics of the deterministic, one-patch
version of the model. Error bars show 95\% confidence intervals of extinction risk.
LT=left tail, RT=right tail. \label{fig_extrisk_d0}}
\end{center}
\end{figure*}

<!--Summary tables for 6 nonlinear models used-->
```{r read_tab_model_summary,echo=F, results='asis',message=F}
library(kableExtra)
library(dplyr)
dt<-readRDS("./Results/tab_model_summary.RDS")
knitr::kable(dt, escape=F,
             format="latex",digits=3, align="c",table.env='table*',
             caption = "Summary of the six density-dependent stochastic metapopulation models we used. The population growth rate in habitat patch $i$ at time $t$, $\\lambda_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 $\\varepsilon_i(t)$ in patch $i$ at time $t$. $\\varepsilon_i(t)$ was normally distributed with mean $0$ and standard deviation $\\sigma$, and values were generated independently across times. $P_e$ is the equilibrium of the deterministic, one-patch version of the model. $r_\\text{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_\\text{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_\\text{min}<r<r_\\text{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_\\text{min}$, $r_c$ and $r_\\text{bf}$. See Appendix, section \\ref{SM-model_stability} for the formulas that gave these values.\\label{tab_model_summary}",
             booktabs=T, linesep = "\\addlinespace"
             )%>%column_spec(c(1,3,5,6,7),width="1.4 cm")%>%kable_styling(latex_options = "basic",position = "center")%>%column_spec(2, width = "3.8 cm")
```

\clearpage

# Literature cited

\setlength{\parindent}{-0.2in}
\setlength{\leftskip}{0.2in}
\setlength{\parskip}{1pt}
\noindent