MT_spatial_avg.Rmd

---
title: "Species relationships in the extremes and their influence on community stability"
author:
 - Shyamolina Ghosh\textsuperscript{1}, Kathryn L. Cottingham\textsuperscript{2}, Daniel C. Reuman\textsuperscript{1,3,\dag}
 - \textsuperscript{1}Department of Ecology and Evolutionary Biology and Kansas Biological Survey, University of Kansas, Lawrence, KS, 66045, USA; \textsuperscript{2}Department of Biological Sciences, Dartmouth College, Hanover, NH, 03755, USA; \textsuperscript{3}Laboratory of Populations, Rockefeller University, 1230 York Ave., New York, NY 10065, USA
date: '\textsuperscript{\dag}Correspondence: Daniel C. Reuman, 2101 Constant Ave, Lawrence, KS, 66047, reuman@ku.edu, 626 560 7084'
fontsize: 11pt
geometry: "left=1 in,right=1 in,top=1 in,bottom=1 in"

output: 
  pdf_document:
    number_sections: no
    keep_tex: yes
    fig_caption: yes

header-includes:
      - \usepackage{xr} \externaldocument[SI-]{SI_spatial_avg}
      - \input{head_maintext.sty}

mainfont: Times New Roman 
tables: True
link-citations: True
urlcolor : blue
indent : True
csl: vancouver.csl
bibliography: REF_CSS.bib
---

<!--***DAN: store alternative titles here, decide on one closer to the end
1) Skewness ratio: a complimentary measure for variance ratio and community stability
2) A complete description of synchrony and compensatory dynamics yields a better understanding of community stability
3) Asymmetric tail association, skewness and community stability
4) Asymmetric tail association, skewness and community variability 
5) Asymmetric tail association affects community stability
6) Species relationships influence community stability in more ways than classical approaches reveal. This is OK, but
may make people think of it as a methods paper?
7) Through thick and thin: species relationships in extreme times and their specific influence on community stability 
8) Species relationships in extreme times and their specific influence on community stability
-->

```{r setup_Paper, echo=F}
knitr::opts_chunk$set(echo = TRUE, fig.pos = "H")
options(scipen = 1, digits = 4) #This option round all numbers appeared in the inline r code upto specified digit, default value=7
```

<!--Abstract-->
\newpage

\noindent \textbf{ABSTRACT.} 
<!--background-->
Synchrony among population fluctuations of multiple 
coexisting species has a major impact on community stability, i.e., on the relative temporal constancy 
of aggregate properties such as total community biomass. 
However, synchrony and its impacts are usually measured 
using covariance methods, which do not account for whether species abundances may be 
more correlated when species are relatively common than when they are scarce, or vice versa. 
Recent work showed that species commonly exhibit such “asymmetric tail associations”. 
<!--main goal-->
We here consider the influence of asymmetric tail associations on 
community stability. 
<!--new methods-->
We develop a “skewness 
ratio” which quantifies how much species relationships and tail associations 
modify stability. The skewness ratio complements 
the classic variance ratio and related metrics. 
Using multi-decadal grassland datasets, we 
show that accounting for tail associations gives new viewpoints on synchrony and stability. 
<!--results-->
E.g., species associations can alter community stability differentially for
community crashes or explosions to high values, a fact not previously detectable.
Species associations can mitigate explosions of community abundance to high values, increasing
one aspect of stability, while simultaneously exacerbating crashes to low values, decreasing another 
aspect of stability; or vice versa.
<!--summary-->
Our work initiates a new, more 
flexible paradigm for exploring species relationships and community stability. 

<!--keywords-->
\noindent \textbf{\textit{Keywords:}} variance ratio, community stability, community variability, 
synchrony, tail association, compensatory dynamics

\newpage

# Introduction\label{Introduction}

Understanding how the dynamics of individual species within an ecological community combine to determine 
the temporal stability of key aggregated properties of the system as a whole is a topic that has 
fascinated ecologists for decades [@macarthur1955fluctuations; @pimm1984complexity; @RN1; @gonzalez2009causes]. 
An important early insight was that an aggregate 
community property such as the total biomass of all species can be relatively stable through time, even 
while the dynamics of individual species are highly variable, so long as different species exhibit offsetting, asynchronous or 
partially asynchronous fluctuations [@peterson1975;@frost1995species].
Such fluctuations are often referred to as compensatory dynamics, because 
decreases in some species' abundances are compensated
for by simultaneous increases in other species [@peterson1975;@frost1995species;@bai2004ecosystem;@hallett2014]. 
However, when species instead exhibit 
dynamics that are positively correlated through time -- synchrony -- the aggregate community property tends to
have increased variability [@RN1;@keitt2008coherent;@loreau2008species;@ma2017climate]. 
Other fields, with which some readers may be more familiar, have very similar notions of synchrony. 
For instance, synchronous volume fluctuations of insect mating calls can produce large oscillations in the total volume of
an acoustic signal [@Sheppard2020]. Thus, species relationships, 
and specifically the degree of synchrony between the population dynamics 
of different species, are important contributors to the stability of an aggregate community property.
Here, the terminology "species relationships" encompasses direct species interactions as well as related 
responses to environmental and other drivers. In addition to total biomass, aggregated community
properties can include total abundance, net primary productivity, decomposition rate, and various nutrient cycling rates.

Synchronous versus compensatory dynamics, and their influence on community stability, have been studied 
in a wide variety of ecological communities [e.g., [@peterson1975; @RN13; @houlahan2007]].
<!--, including marine benthic invertebrates (e.g., @peterson1975); 
freshwater plankton (e.g., @frost1995species; @RN6; @RN7; @vasseur2007; @downing2008; @keitt2008coherent; @RN10;
@Vasseur2014; @brown2016); butterflies and moths [@houlahan2007; @RN13]; small mammals [@RN14; @RN15; @RN16];
birds (e.g., @RN17); and plants in deserts (e.g., @RN14; @RN15), grasslands (e.g., @bai2004ecosystem; @houlahan2007;
@hallett2014), and forests (e.g., @houlahan2007; @RN18). 
These and other studies suggest that synchrony is often -->
Synchrony is often due to common responses to environmental 
changes [e.g., [@peterson1975; @RN19; @RN6; @houlahan2007]], whereas 
compensatory dynamics result from multiple mechanisms, including competitive interactions and differing responses 
to environmental changes [e.g., [@RN19; @gonzalez2009causes; @houlahan2007; @RN21; @RN22]]. Severe environmental 
stressors that impact all populations are expected to synchronize dynamics, whereas stressors with differential 
impacts allow for compensatory dynamics [@RN23]. 
Although past work on the general concept of "community stability" has led to myriad precisely defined alternative measures
of "stability" [@Ives2003], we here focus on the common practice of examining the variability, through time,
of an aggregate community property such as the total biomass of all species, referring to this henceforth 
as "community variability."

Much past work on species relationships and their influence on community variability has examined covariances
between population time series, to characterize species relationships, and the temporal variance of an 
aggregate community property such as total biomass, to characterize community variability 
[@peterson1975; @schluter1984variance; @RN23; @frost1995species; @loreau2008species; @Gross2014].
For instance, if $x_i(t)$ are data representing the population abundance or biomass 
of species $i=1,\ldots,N$ at the sampling times
$t=1,\ldots,T$, and if $\mu_i=\mean(x_i)$, $v_{ii}=\var(x_i)$, and $v_{ij}=\cov(x_i,x_j)$ (here $j=1,\ldots,N$
is another species index), community variability
has commonly been quantified using the squared coefficient of variation of $\xtot(t)=\sum_i x_i(t)$, i.e., 
$\CVcomsq=\var(\xtot)/\left(\mean(\xtot)\right)^2=\left(\sum_{i,j} v_{ij}\right)/\left(\sum_i \mu_i\right)^2$.
This quantity obviously connects to species relationships measured via 
the covariances, $v_{ij}$, which appear directly in the formula.
Denoting by 
$\CVindsq=\left(\sum_i v_{ii}\right)/\left(\sum_i \mu_i\right)^2$ the value that $\CVcomsq$ would take if 
the dynamics of each species were independent of the dynamics of other species (so $v_{ij}=0$ for all $i \neq j$), the 
classic *variance ratio* was defined [@peterson1975] as 
$\phicv=\CVcomsq/\CVindsq=\var(\xtot)/\left(\sum_i v_{ii}\right)=\left(\sum_{i,j} v_{ij}\right)/\left(\sum_i v_{ii}\right)$. 
Thus $\CVcomsq=\phicv \CVindsq$.
Because $\CVindsq$ has been interpreted as what community variability would be, without species relationships,
$\phicv$ is the factor by which species relationships inflate or decrease community variability.
The variance ratio has therefore been used [@frost1995species;@hallett2014] as an index of whether dynamics 
are synchronous or compensatory. If $\phicv>1$, the interpretation has been that dynamics are, on balance,
synchronous because then $\CVcomsq>\CVindsq$, i.e., the aggregate property $\xtot$ is more variable than
it would be if species dynamics were independent. Conversely, if $\phicv<1$, the interpretation has been
that dynamics are compensatory [but see @loreau2008species, and below, for alternative perspectives].

Extensions, improvements and alternatives to the original variance ratio have been proposed, notably the synchrony metric of 
Loreau and de Mazancourt @loreau2008species, here denoted $\phiLdM$. 
Whereas newer metrics such as $\phiLdM$ have interpretive and other advantages over the 
original variance ratio [@loreau2008species; @Thibaut2013], in this study
we show that the ways in which species relationships influence community variability are likely to be 
more nuanced than is revealed by either the variance ratio approach or
by any existing extensions or alternatives. Whereas the variance ratio 
compares $\CVcomsq$ to the baseline value $\CVindsq$ that 
would be the community $\text{CV}^2$ under an assumption of species independence (by considering the ratio 
$\phicv=\CVcomsq/\CVindsq$), $\phiLdM$
instead compares $\CVcomsq$ to the baseline value $\CVsynsq$ that would be the community $\text{CV}^2$ under an assumption of 
perfect linear correlation between species dynamics (again by considering a ratio): $\phiLdM=\CVcomsq/\CVsynsq$. 
It is straightforward to show [@Thibaut2013] that $\CVsynsq=\left(  \sum_i \sqrt{v_{ii}} \right)^2/\left(\sum_i \mu_i\right)^2$,
i.e., that this expression is what community $\text{CV}^2$ would be under perfect linear correlation between species
[@Thibaut2013]; and therefore that $\phiLdM=\left( \sum_{i,j} v_{ij}  \right)/\left( \sum_i \sqrt{v_{ii}}  \right)^2$. Timescale-specific alternatives and extensions of the original variance ratio approach have also been developed
[@vasseur2007; @brown2016;@zhao2020]. Although each alternative metric improves upon the original variance ratio in 
at least some respects, we nevertheless construct the specific statistical
tools of this paper as extensions of the original variance ratio. This is for mathematical simplicity, and also
because the choice makes little difference: the concepts of this study can be applied as 
improvements to all the existing frameworks of which we are aware.  We return to this topic in the Discussion.

<!--Old version follows, adventurously trialing new versions, so keeping the old one here: 
Much past work on community variability and the influence of species relationships on it has proceeded by examining 
the coefficient of variation through time of an aggregate community property, 
together with covariances between population time series, to characterize species 
relationships [@peterson1975; @schluter1984variance; @RN23; @frost1995species].
For instance, if $x_i(t)$ are data representing the population abundance or biomass 
of species $i=1,\ldots,N$ at the sampling times
$t=1,\ldots,T$, and if $\mu_i=\mean(x_i)$, $v_{ii}=\var(x_i)$, and $v_{ij}=\cov(x_i,x_j)$ (here $j=1,\ldots,N$
is another species index), community variability
has commonly been quantified using the squared coefficient of variation of $\xtot(t)=\sum_i x_i(t)$, i.e., 
$\CVcomsq=\var(\xtot)/\left(\mean(\xtot)\right)^2=\left(\sum_{i,j} v_{ij}\right)/\left(\sum_i \mu_i\right)^2$.
This quantity obviously connects to species relationships measured via 
the covariances, $v_{ij}$, which appear directly in the formula. Denoting by 
$\CVindsq=\left(\sum_i v_{ii}\right)/\left(\sum_i \mu_i\right)^2$ the value that $\CVcomsq$ would take if 
the dynamics of each species were independent of the dynamics of other species (so $v_{ij}=0$ for all $i \neq j$), the 
classic *variance ratio* was defined [@peterson1975] as 
$\phicv=\CVcomsq/\CVindsq=\var(\xtot)/\left(\sum_i v_{ii}\right)=\left(\sum_{i,j} v_{ij}\right)/\left(\sum_i v_{ii}\right)$. 
Thus $\CVcomsq=\phicv \CVindsq$. Because $\CVindsq$ is interpreted as what community variability would be, without species relationships,
$\phicv$ is the factor by which species relationships inflate or decrease community variability.
The variance ratio has therefore commonly been 
used [@frost1995species;@houlahan2007;@hallett2014;@ma2017climate] as an index of whether dynamics 
are synchronous or compensatory. If $\phicv>1$, the interpretation has been that dynamics are, on balance,
synchronous because then $\CVcomsq>\CVindsq$, i.e., the aggregate property $\xtot$ is more variable than
it would be if species dynamics were independent. Conversely, if $\phicv<1$, the interpretation has been
that dynamics are compensatory.

Several extensions of the variance ratio approach have been proposed, but the original variance ratio 
is still commonly used. Timescale-specific extensions [@vasseur2007;
@brown2016;@zhao2020] have revealed that it is possible for synchronous and compensatory 
dynamics to occur simultaneously on different timescales. Other alternatives to $\phicv$ 
have also been proposed [@loreau2008species]. 
Here, we show that the ways in which species relationships influence community variability are likely to be 
more nuanced than is revealed by either the classic variance ratio approach or
by existing extensions. -->

We begin by giving an intuitive introduction to the main concepts of this study; 
we do so by considering two simulated community dynamics datasets 
(Fig. \ref{fig_pedag_cv2_skw}A, B; data generated as in Supporting Information (SI) 
section \ref{SI-PedagogFigMethods}, see also SI Fig. \ref{SI-fig_pedag_cv2_skw_SM}).
Both communities consist of $N=20$ species. Data were generated so that species marginal distributions 
(and therefore the means, $\mu_i$, and variances, $v_{ii}$) were the same, up to sampling variation, 
in both scenarios. In other words, the distribution of abundances, through time, of any given species $i$ was the same for 
both scenarios. Likewise, data were 
generated so that species covariances, $v_{ij}$ for $i \neq j$,
were essentially the same for Fig. \ref{fig_pedag_cv2_skw}A and B. Thus both $\CVcomsq$ and $\CVindsq$, which 
depend only on the $\mu_i$ and $v_{ij}$,
were the same in both scenarios, and so $\phicv$ was the same. Because $\CVsynsq$ also depends only on 
the $\mu_i$ and $v_{ii}$, $\phiLdM$ was also the same in our two scenarios, up to sampling variation
(Fig. \ref{fig_pedag_cv2_skw}A, B). Thus classical approaches do not distinguish between the scenarios, neither in the nature
of species relationships nor in the effects of relationships on community aggregate dynamics $\xtot(t)$.
Any approach that considers only species marginal distributions through time
and the covariances $v_{ij}$ will not detect the differences between our two scenarios, by construction. The 
metric of @Gross2014 is another such.

Nevertheless, species relationships and aggregate dynamics differed in important and 
related ways in our scenarios. In both scenarios,
the 20 species could be separated into two groups of 10 species (the red lines on Fig. \ref{fig_pedag_cv2_skw}A, B are one 
group and the black lines are the other), 
and the dynamics of species in the different groups were compensatory. However, in scenario 1 (Fig. \ref{fig_pedag_cv2_skw}A), 
species in the same group exhibited strongly synchronous dynamics when those species were scarce (below about 5.5 on the 
$y$-axis of Fig. \ref{fig_pedag_cv2_skw}A), 
and much less synchronous dynamics when those species were common (above about 5.5); whereas in scenario 
2 (Fig. \ref{fig_pedag_cv2_skw}B) species in the same group exhibited strongly synchronous dynamics 
when common and much less synchronous dynamics when scarce. Dynamics like scenario 1 could occur ecologically
if red species in Fig. \ref{fig_pedag_cv2_skw}A were sensitive to an environmental factor, $F$, but only when it 
goes below a threshold. So all the red species are controlled by the same factor, $F$, when it is below the threshold,
and hence are synchronous when scarce. 
When $F$ is above the threshold, the red species are each sensitive, instead, to other, 
distinct factors, rendering them asynchronous when common. Black species 
of Fig. \ref{fig_pedag_cv2_skw}A may also be sensitive to $F$, but in a reverse manner, being all adversely impacted 
by $F$ when it is above the threshold, and unaffected by $F$ below the threshold. 
Resulting community dynamics 
(Fig. \ref{fig_pedag_cv2_skw}C-F) differed between scenarios because of the distinct species relationships.
The aggregate quantity, $\xtot$, was more likely to crash to low values in the first scenario (Fig. \ref{fig_pedag_cv2_skw}C, E)
than in the second scenario (Fig. \ref{fig_pedag_cv2_skw}D, F; compare the probabilities on E and F listed to the left of the dashed lines), 
whereas $\xtot(t)$ was more likely to explode to high values 
in the second scenario than in the first (compare the probabilities on E, F to the right of the dotted lines).
Thus total abundance was left-skewed, frequently crashing to low values, for scenario 1 (C, E); 
but was right skewed, frequently exploding to high values, in scenario 2 (D, F). 

It is not difficult to imagine ecological ramifications or 
applied significance of the distinction illustrated here. For instance, grazers subsisting on a grass mixture will 
have different growth prospects if that mixture exhibits occasional crashes (which may harm the grazers) 
compared to if the mixture shows occasional explosions of abundance (unlikely to harm them). 
Thus species relationships and their effects on community aggregate dynamics may differ in potentially ecologically
meaningful ways not detected by classical approaches. 
The reliance of earlier approaches on variances, covariances and related linear tools
renders those approaches unable to register the differences illustrated here.

This study addresses three gaps in the research literature. First, commonly used past descriptions of species 
relationships have been oversimplified, based primarily on covariances or correlations of species time series.
The example above suggests that more nuanced measures should be developed that identify whether species 
abundances are more related to each other when the species are common, or scarce (Fig. \ref{fig_pedag_cv2_skw}). 
Our first goal (G1) is to develop 
an approach to studying relationships
between species that complements classic approaches relying on covariances. 
Second, past measures of community variability have typically focused on the variance or 
coefficient of variation of an aggregate property like  $\xtot$. The example above suggests
that more nuanced or complementary measures should be developed, taking skewness into account 
(Fig. \ref{fig_pedag_cv2_skw}). 
Our goal (G2) is to develop an
approach to studying community variability that complements the coefficient of variation.
Finally, the example above suggests that means by which species relationships
influence community variability should be studied using the new measures of G1 and G2.
Our goal (G3) is to do so, using data from long-term studies from mixed-grass prairie in Hays, Kansas and 
from the Konza Prairie Long-Term Ecological Research (LTER) site, also in Kansas.

The core ideas of this study relate to extreme events, and to how relationships between constituents of 
an aggregate quantity change during such events. There are strong connections to
risk pricing in finance. A portfolio of investments should reduce risk if constituent investments 
are negatively correlated or uncorrelated. However, if constituents become positively correlated, risk can be accentuated. 
Some observers have blamed the 2008 financial crisis on analysts estimating correlations between housing default
risks based on data from normal times, not taking into account that in a crisis, default risks become more correlated.
In crisis, financial products based on aggregating mortgages therefore became riskier. 
A single modelling approach [@li2000] was apparently widely used [@salmon2009]. The relationship between two 
real quantities, be they investment risks or species populations, can probably never be completely captured by a single number
such as a covariance [@cheung2013], nor considered to be constant through time. Just as for mortgage defaults risks, 
the abundances of two species may
be more or less correlated under different circumstances. Classical approaches that characterize relationships
between species using covariances may therefore neglect something important. For instance, correlations between species 
may change as global climate change progresses, or when extreme climatic events occur [@Ghosh_extinction]. 
Financial analysts 
have learned these lessons and have begun developing more appropriate models [@cheung2013]. Ecologists 
can benefit from application of the same methods. Overall, this study initiates a new and more flexible paradigm for
conceptualizing synchronous and compensatory dynamics in ecosystems, and their influence on community variability.

<!--Synchronous versus compensatory dynamics in communities, and their influence on community stability, has most commonly
been studied in plankton systems [@downing2008;@brown2016] and in grasslands 
[@bai2004ecosystem;@houlahan2007;@hallett2014]. Evidence suggests that synchronous dynamics may be dominant,
and compensatory dynamics rare, in plankton systems [@Vasseur2014]. On the other hand,
compensatory dynamics may be at least somewhat more common
in grasslands [@hallett2014;@zhao2020]. Though some studies have also found compensatory dynamics rare in
grasslands [@houlahan2007;@valone2008], others have found general evidence for compensatory dynamics [@bai2004ecosystem;@hector2010] 
or have found compensatory dynamics selectively, in some contexts [@grman2010;@hallett2014]. 
One outcome of our empirical work is a suggestion that compensatory dynamics in grasslands may be less 
common than was previously appreciated, when one adopts a broader conceptual view of what constitutes synchronous versus 
compensatory dynamics.--> 

<!--***DAN: Not sure how relevant this stuff is, but keeping it around for now, partly for the references.
There are reasons to focus on skewness to better understand the community-stability 
in this rapidly changing world. Earlier studies [@smith2011ecological; @hoover2014resistance; 
@felton2017integrating] showed climatic events that are extreme and 
sometimes spatially well-spread impact the net productivity of community [@ciais2005europe],
community assembly [@thibault2008impact] and even induce a shift in 
ecotone boundaries [@allen1998drought]. Not only the mean and variance, the symmetry of 
temporal distribution for such climatic signals is also changing [@ummenhofer2017extreme]. This 
change in symmetry results into change in the skewness that can act on the community as well 
as on ecosystem level. Another reason is like variance changing skewness had been 
used earlier as an indicator for changing ecosystem (i.e. for an early 
warning signal, regime shift, etc.) [@guttal2008changing; @guttal2009spatial; 
@dakos2010spatial; @dakos2012methods]. But to our knowledge, there is a gap to incorporate 
this skewness into direct community-stability metric.-->

# Theory \label{Theory}

\noindent We now illustrate theoretically how relationships between
species, beyond covariances, can influence variability of the 
total community. The theory here relates to, but goes substantially beyond, the example of Fig. \ref{fig_pedag_cv2_skw}. 
As in the Introduction, let $x_i(t)$ be an abundance 
measure for species $i=1,\ldots,N$ at time $t=1,\ldots,T$. We use $N=20$ and $T=100000$ in the theoretical examples here.
This large value of $T$ was used to ensure that sampling effects do not obscure our theoretical results.
We assume the multivariate random variables $x(t)=(x_1(t),\ldots,x_N(t))$ are independent
and identically distributed for distinct times, $t$. It should be straightforward to replace
this assumption with an appropriate assumption of stationarity and ergodicity, with no real change
to results, but we adopt the stronger assumption for simplicity. 
To describe species relationships, 
classic approaches use the covariances
$v_{ij}=\cov(x_i,x_j)$, and related quantities. We will instead consider how 
$x_i$ and $x_j$ are related in their distribution tails.
To describe community variability, classic
approaches use the variance of 
$\xtot=\sum_i x_i$, or its coefficient of variation. But these are only summaries of the 
distribution of $\xtot$. We will consider the 
whole distribution, including its skewness. 

We present four scenarios (scenarios 1-2 below, and 3-4 in SI section \ref{SI-TheoryMethods} and 
SI Fig. \ref{SI-fig_theory_SM}). For all 
scenarios, the marginal distributions $x_i$ and
the covariances $v_{ij}$ are the same, as are the quantities $\var(\xtot)$, $\CVcomsq$, $\CVindsq$, $\CVsynsq$, 
$\phicv$, and $\phiLdM$. Thus, as in the example of Fig. \ref{fig_pedag_cv2_skw}, 
our scenarios are indistinguishable to classic approaches but are
nevertheless ecologically distinct: the total abundance $\xtot$ achieves 
more extreme values under some scenarios than others. Hence our scenarios represent communities showing 
different characteristics of variability. 
Scenarios are intentionally simplified, with complexities such as autocorrelation excluded, for clarity, 
but the ideas also apply to real communities (see below). 
We additionally compare each scenario to a reference *comonotonic* scenario (called scenario C), for which 
the $x_i$ are related via perfect positive relationships. 
This is a scenario of maximal covariance between species, 
and therefore maximal values of $\var(\xtot)$ [@cheung2013] and of community variability, given
fixed species marginal distributions, $x_i$.

In scenario 1, which is a baseline, 
$(x_1,\ldots,x_N)$ is a multivariate normal distribution with mean $(0,\ldots,0)$ and covariance matrix
having $1$s along the diagonal and off-diagonal entries $0.6$. A sample from $(x_1,x_2)$ (Fig. \ref{fig_theory}A)
and marginal histograms (side panels of Fig. \ref{fig_theory}A) help illustrate the distribution. 
The standard deviation of $\xtot$ was $`r round(sqrt(vartotX1),1)`$, and its skewness
was close to $0$ (Fig. \ref{fig_theory}B). The probability of $\xtot$ exceeding a high threshold is given 
in Fig. \ref{fig_theory}C for a range of high thresholds. Likewise the probability of $\xtot$ falling below a
low threshold is in Fig. \ref{fig_theory}D for a range of low thresholds. These probabilities represent the 
chances that $\xtot$ will surpass or fall below a threshold considered disastrous or unacceptable from an
applied or ecological viewpoint, and we refer to these as *disaster thresholds*. It is possible to imagine 
applied scenarios in which exceeding a high threshold is disastrous, and other scenarios in which disaster
instead occurs when $\xtot$ falls below a low threshold. Our theory can be considered in relation to 
either of these applied scenarios, so we keep them both in mind while developing the theory.

For the comonotonic scenario, C, we assume $x_i=x_1$ for all $i$. 
This is a specific form of comonotonicity [@cheung2013], with the relationship between $x_1$ and $x_2$
pictured in Fig.  \ref{fig_theory}E.
Because synchrony between the dynamics of individual species is perfect, $\sd(\xtot)$ is much larger
than in scenario 1 (Fig.  \ref{fig_theory}F; the dashed black lines on that panel and the panels below it are 
reproduced from 
Fig. \ref{fig_theory}B, to facilitate comparisons). 
Probabilities of $\xtot$ exceeding a high disaster threshold or falling
below a low one are likewise larger (red lines in Fig.  \ref{fig_theory}G-H; dashed black lines on those panels
and the panels below them are reproduced from panels C and D, respectively).

Scenarios 2-4 illustrate the ideas, previously little recognized in ecology, 
of tail comonotonicity [@cheung2013] and tail associations [@Ghosh_extinction;@Ghosh_copula;@Ghosh_aphidplankton], and the consequences of tail
associations for $\xtot$ and its extreme values. 
In scenario 2, as in all our scenarios, each $x_i$ is standard-normally distributed (Fig. \ref{fig_theory}I side panels),
but $(x_1,\ldots,x_N)$ is not a multivariate normal distribution. Instead, $(x_1,\ldots,x_N)$ was engineered so
that $x_i$ and $x_j$ are perfectly related when they exceed some threshold, but imperfectly related 
below the threshold (Fig. \ref{fig_theory}I; see SI section \ref{SI-TheoryMethods} for 
details of each scenario). The strength of association below the threshold was chosen 
to make $\cov(x_i,x_j)$ equal to the same value as in scenario 1, so $\sd(\xtot)$ was also
the same, up to sampling variation, as in scenario 1 (Fig. \ref{fig_theory}J). (Recall that 
$\sd(\xtot)=\sqrt{\var(\xtot)}$ and $\var(\xtot)=\sum_{i,j} \cov(x_i,x_j)$.)
Comonotonicity in the upper tails produces right/positive skewness in $\xtot$ (Fig. \ref{fig_theory}J).
Skewness translates to elevated probabilities of $\xtot$ exceeding high disaster thresholds, compared to 
scenario 1 (blue lines in Fig. \ref{fig_theory}K; the dashed red lines on panels K and L are copied from the solid red lines of panels G and H, respectively, to facilitate comparisons). 
The probability of exceeding sufficiently high disaster thresholds is actually the
same as scenario C (Fig. \ref{fig_theory}K). So although scenarios 1 and 2 are 
indistinguishable to classic approaches, the probability 
of disaster in scenario 2, for high disaster thresholds, is much elevated, and is 
as high as it can possibly be given the species marginal
distributions used. Although it does not appear unstable to classic approaches, the scenario 2 community
is maximally unstable with respect to high disaster thresholds.

Scenario 3 parallels scenario 2, but with comonotonicity in the left tails instead of the right, and concomitant
elevated probabilities of $\xtot$ falling below a low disaster threshold. 
Scenario 4 shows comonotonicity in both tails, and thus elevated probabilities both of $\xtot$ exceeding high
disaster thresholds and falling below low ones (SI section \ref{SI-TheoryMethods} and Fig.  \ref{SI-fig_theory_SM}). 

Although scenarios 1-4 are indistinguishable to classic approaches, 
the probabilities of disaster are elevated in scenarios 2-4, and are maximal, given fixed species marginal distributions, 
for some disaster types (high or low disaster thresholds), depending on the scenario. 
Perfect comontonicity in
the tails of population distributions $x_i$ simplifies our presentation of the theoretical ideas, but is not 
necessary to generate skewness in $x_{\text{tot}}$ and elevated disaster probabilities. Similar results pertain 
for essentially any case where associations between species are asymmetric in the tails.

Our scenarios were constructed using mathematical results of @cheung2008 and @cheung2013. We refer the reader
to those papers for mathematical specifics, while here summarizing aspects important for ecological
applications. Given one-dimensional cumulative distribution functions (CDFs) $F_1,\ldots,F_N$, the 
\emph{Fr\'{e}chet space} $\mathcal{R}(F_1,\ldots,F_N)$ is the set of all $N$-dimensional random vectors 
$(x_1,\ldots,x_N)$ with the $F_1,\ldots,F_N$ as their marginal distribution functions. Thus the $\text{Fr{\'e}chet}$ space
represents, in our modelling context, all possible interspecific relationships (different degrees and kinds of 
synchrony and compensatory dynamics between species) that can pertain, given fixed
species marginal distributions. So the question of what forms the distribution of $\xtot$ can take across the $\text{Fr\'{e}chet}$ space
is a precisely formulated version of the classic question of how species interrelationships 
such as synchrony and compensatory dynamics can influence community variability.
It is well known to statisticians that the sum
$\xtot$ exhibits maximal variance when the $x_i$ are comonotonic [@cheung2013]; this is the maximal-synchrony case.
Cheung and Vanduffel @cheung2013 showed the converse, i.e., if $\var(\xtot)$ is maximal in the $\text{Fr\'{e}chet space}$, then 
$(x_1,\ldots,x_N)$ must be comonotonic (subject to some mild regularity assumptions about the $F_i$). 
The classical variance ratio and Loreau-de Mazancourt approaches essentially make sole use of
$\var(\xtot)$ to quantify community variance ($\CVcomsq$ is used, actually, but this is equivalent, since species 
marginals are fixed within the $\text{Fr\'{e}chet space}$). But @cheung2013 showed that, 
even if we consider only the sub-$\text{Fr\'{e}chet}$ space
consisting of random vectors $(x_1,\ldots,x_N)$ in $\mathcal{R}(F_1,\ldots,F_N)$ that additionally have 
$\var(\xtot)=S$ for some fixed $S$, there are random vectors that achieve extreme values just as extreme as the
comonotonic case (e.g., scenarios 2-4 in Fig. \ref{fig_theory} and SI Fig. \ref{SI-fig_theory_SM}). Existing approaches 
in the ecological literature tell us nothing
about true community variability, in the sense that if species relationships are tail associated,
the community can exhibit the same extreme behavior as the comonotonic case, while the classically important
quantity $\var(\xtot)$ gives the impression of a much more stable community. We next transition to describing the data
and the statistics used to assess to what extent the above theory applies to those data.

# Methods\label{Methods} 

## Data

<!--Data prep, and common/intermediate/rare species-->
\noindent We used data from two grassland sites, both monitored for decades. 
After initial processing, data from the Hays, Kansas site consisted of annual estimates of basal cover, 
averaged over 36 quadrats from which livestock were excluded, for each species present, for the years 1932-1972
[@adler2007long]. 
Data from Konza Prairie, after initial processing, consisted of annual canopy percent cover data by species for the years
1983-2018, averaged over 20 annually burned plots on tully soils, ungrazed by livestock [@Hartnett2019]. 
Additional details are in SI section \ref{SI-Data}.

Community data often include numerous rare species. Rare species complicate analyses of interspecific relationships 
and their effect on community variability because the relationship of a rare species with another species is difficult to accurately
assess. We therefore categorized species as "common" (present in the system for at least 35 years),
"rare" (present for at most two years), or "intermediate" (other species). Rare and intermediate species were combined
into a single pseudo-species for analyses. Common species (Tables
\ref{SI-tab_spinfo_hays} and \ref{SI-tab_spinfo_knz}) made up the overwhelming majority of both the Hays and Konza 
communities (Fig. \ref{SI-fig_spaceavg_timeseries}), so "intermediate" species were actually quite uncommon, and 
it is the variability of common species, and their
relationships, that mainly determine community variability. Henceforth the terminology "species" signifies
both true species and the pseudo-species formed by combining rare and intermediate species.

## Statistical methods: quantifying tail associations

<!--Normalized rank plots and partial spearman-->
\noindent Fig. \ref{fig_pedag_cv2_skw} and the Theory section
suggest that to accomplish goal G1 of the Introduction, we could apply measures of tail association between species, 
and asymmetries of tail association. We here define such measures, based on the *partial Spearman correlation*
of @Ghosh_copula. Notation is summarized in Table \ref{SI-fig_notation}. Given $x_i(t)$ for $t=1,\ldots,T$, we begin by defining the
*normalized rank* $u_i(t)$, equal to the rank of $x_i(t)$ in the set $\{x_i(1),\ldots,x_i(T)\}$, 
divided by $T+1$. The smallest element of this set is here considered to have rank 1. 
Plotting the normalized ranks $u_i(t)$ and $u_j(t)$ against each other for two positively associated species 
$x_i$ and $x_j$ provides a visual
indication of tail association and asymmetries in tail association (e.g., the points of Fig. \ref{fig_pedag_taildep},
see also SI Fig. \ref{SI-fig_pedag_taildep_SM}), 
as described by @Ghosh_copula.
Given two bounds $0\leq b_l < b_u \leq 1$ ($l$ stands for "lower" and $u$ for "upper"), we define the lines $u_i + u_j=2b_l$ and $u_i+u_j=2b_u$,
which intersect the unit square $[0,1]\times[0,1]$ of the normalized rank plot (Fig. \ref{fig_pedag_taildep} shows
these bounds for $b_l=0.1$ and $b_u=0.25$). The partial Spearman
correlation associated with the band circumscribed by these bounds is
\begin{equation}
\cor_{b_l,b_u}(x_i,x_j) = \frac{\sum (u_i(t)-\mean(u_i))(u_j(t)-\mean(u_j))}{(T-1)\sqrt{\var(u_i) \var(u_j)}},\label{PartialSpearmanEq}
\end{equation}
\noindent where sample means and variances are computed over $t=1,\ldots,T$ but the summation in the 
above formula is only computed over $t$ such that $u_i(t) + u_j(t)>2b_l$ and $u_i(t) + u_j(t)<2b_u$
(e.g., the points in the band delineated by the two diagonal lines on Fig. \ref{fig_pedag_taildep}).
The partial Spearman correlation is the component of the Spearman correlation which can be attributed to
the points in the band. For two positively associated variables, $\cor_{0,0.5}$ measures 
lower-tail association, and $\cor_{0.5,1}$ measures upper-tail association. 
The difference $\cor_{0,0.5}-\cor_{0.5,1}$ is a way to measure asymmetry of tail
associations. Positive values of this statistic indicate stronger lower-tail association, and negative values 
indicate stronger upper-tail association. We henceforth use the shorthand $\cor_l$ for $\cor_{0,0.5}$ and 
$\cor_u$ for $\cor_{0.5,1}$. 
Genest and Favre @Genest2007 recommend making inferences about dependence structures such as tail associations using 
rank-based approaches such as the tools we introduced here. 
Work of Ghosh and colleagues [@Ghosh_copula; @Ghosh_extinction; @Ghosh_aphidplankton] provides additional information on
why rank-based approaches and the partial Spearman correlation are an appropriate statistical choice 
for purposes such as ours. 

<!--Aggregating to the community-->
We computed $\cor_l-\cor_u$ for all positively associated species pairs (those with positive overall Spearman correlation); 
negatively associated pairs were ignored because theory suggested the importance of positively associated species pairs
(but see SI section \ref{SI-stats_details} for thoughts on future work). 
Fig. \ref{fig_pedag_cv2_skw} and the Theory section
suggest that communities with a preponderance of upper-tail association between positively
associated species will have $\xtot$ with distinct distributional properties from communities
with more lower-tail association. Community-aggregated tail association between positively 
associated species was measured by counting the number, $n_L$, of values
$\cor_l(x_i,x_j)-\cor_u(x_i,x_j)$, for $i \neq j$, that were positive, representing stronger lower- than
upper-tail association between positively associated species; as well as the  number, $n_U$, of values 
that were negative, representing stronger upper- than lower-tail association.
We also computed the sum, $\Atot$, of $\cor_l(x_i,x_j)-\cor_u(x_i,x_j)$
across all positively associated species pairs, $i,j$. This was positive if lower-tail 
association was stronger, in aggregate, than upper-tail association, and
negative if the reverse. We call $\Atot$ the *total community tail association*. 
Note that the word "positive" applies in three senses, here, all distinct: positive association 
of species pairs; positivity of $\cor_l(x_i,x_j)-\cor_u(x_i,x_j)$; and positivity of $\Atot$.
<!--***DAN: Do some search and standardize the use of lower v. left-tail and upper v. right-tail
association, and related quantities
***Shya: I did, it is always lower-tail and upper-tail association. Left- and right- tail for skewness.-->
<!--***DAN: Having defined this new piece of terminology, make sure to use it.-->

<!--omitted as of now: Significance tests for individual and overall asymmetry, tailsignif.R, binomial_sigtest.R
At present not shown, because we want to consider all interactions: weekest to strongest.

Tail-asymmetry could be weaker on individual species-pair, but on an aggregate level can make 
the community strongly tail-asymmetric. To test significantly stronger 
association in one of its extreme tail for a given species-pair, we simulate $10000$ made-up 
species pair with same Spearman correlation but removing their tail-asymmetry. If 
estimated tail-asymmetry $\cor_{l}-\cor_{u}$ on given data lies outside 
the $95\%$ confidence interval of the ones estimated with $10000$ surrogated data, then 
we consider the tail-asymmetric relationship is significant for that species-pair. 
On a community level, we calculated the fraction of significant species-pair among all correlated 
pairs considered (overall significance on the community level, $C_s$). We also perform a binomial 
test to check whether the community shows significantly stronger lower-tail (or upper-tail) association.
-->

<!--2. For Q2:
      - brief intro about the usual way to measure community variability (variance ratio method)
      - then introduce skewness ratio in analogy with variance ratio
      - separate the contributions from tail association and usual correlation (Pearson) 
      - need schematic diagram to show the relation (this diagram will again help to interpret results for Q2 later)
      - somewhere, in results or in discussion we should talk about the aggregated effect of weak interactions
-->

## Statistical methods: quantifying variability, and the skewness ratio 

\noindent To accomplish goal G2 of the Introduction, we need a measure of the distribution $\xtot$ that
characterizes community variability in a way that complements $\var(\xtot)$ or $\CVcomsq$. We use
the skewness, $\sk(\xtot)$ (see SI section \ref{SI-stats_details} for details). 
As illustrated in Fig. \ref{fig_pedag_cv2_skw}, skewness complements use of 
$\var(\xtot)$ or $\CVcomsq$ to characterize community variability, because strongly negative values indicate that $\xtot$
undergoes large downward departures from typical values (crashes), and large positive values
indicate that $\xtot$ undergoes upward departures (explosions); either of these dynamical behaviors 
can happen independently of the 
values of $\var(\xtot)$ and $\CVcomsq$. Because $\sk(\xtot)$ relates to upward
or downward departures of $\xtot$ from typical values, it also relates to probabilities that $\xtot$
will surpass a high disaster threshold (see Theory) or fall below a low one. 

Just as the variance ratio $\phicv$ is the quotient $\CVcomsq/\CVindsq$, we 
analogously define a *skewness ratio*. Let $\Scom=\sk(\xtot)$, and consider the value that 
$\Scom$ would take if species dynamics were independent, i.e., 
$\Sind=\frac{\sum_i m_3(x_i)}{(\sum_i v_{ii})^{3/2}}$, where $m_3(x_i)$ is the third central moment of $x_i$ 
(SI section \ref{SI-stats_details}).
The *skewness ratio* is $\phi_S=\Scom/\Sind$, so that $\Scom=\phi_S \Sind$. It is also possible to define 
a skewness metric in a manner that attempts to parallel the choices made in constructing $\phiLdM$. This is 
elaborated in the Discussion.

<!--phi_S > 1, anmd assuming \phicv \approx 1-->
Values of the skewness ratio provide a summary of how species relationships influence community variability that 
complements information provided by the variance ratio. The skewness ratio can be positive or negative because 
$\Scom$ and $\Sind$ can be positive or negative, unlike $\CVcomsq$ and $\CVindsq$. If $\phi_S>1$, 
$\Scom$ has greater magnitude, but the same sign, as what its value would be in 
the absence of species relationships ($\Sind$). 
Assuming, for simplicity $\phicv \approx 1$ (see SI section \ref{SI-phiSinterp} and Table \ref{SI-further_discussion_table} for more on this assumption),
if $\Sind$ is positive, then $\phi_S>1$ means $\Scom>\Sind>0$, i.e., species relationships magnify the tendency
for the community aggregate quantity $\xtot$ to explode to high values;
and if $\Sind$ is negative, then $\phi_S>1$ means $\Scom<\Sind<0$, i.e., species relationships magnify the tendency
for $\xtot$ to crash to low values.
This can be viewed as a new type of synchrony, in the sense that species relationships
inflate community variability by accentuating the probability that $\xtot$ will exceed large 
disaster threshold, or fall below small ones. 

<!-- 0<phi_S<1, and assuming \phicv \approx 1-->
If $0<\phi_S<1$, then $\Scom$ has lower magnitude, but the same sign, as $\Sind$. 
Again assuming, for simplicity, that $\phicv \approx 1$ (see SI section \ref{SI-phiSinterp} for more on this assumption),
if $\Sind$ is positive, then $0<\phi_S<1$ means $0<\Scom<\Sind$, i.e., species relationships
mitigate the tendency for $\xtot$ to explode to high values;
and if $\Sind$ is negative, then $0<\phi_S<1$ means $\Sind<\Scom<0$, i.e., species relationships mitigate the 
tendency for $\xtot$ to crash to low values. 
This can be viewed as a new type of compensatory dynamics, because species 
relationships reduce a kind of community variability. 

<!--\phi_S<0, and assuming \phicv \approx 1-->
If $\phi_S<0$, what happens? If $\Sind<0$, so that without species relationships $\xtot$ would have exhibited occasional
crashes, we instead have $\Scom>0$: with species relationships $\xtot$ instead
exhibits occasional explosions. If $\Sind>0$, so that without species relationships $\xtot$ 
would have exhibited occasional
explosions, we instead have $\Scom<0$: with species relationships $\xtot$ instead
exhibits occasional crashes. Thus the tendency of $\xtot$ to crash or explode is reversed by 
species relationships. We again assume, here, for simplicity, that $\phicv \approx 1$ -- see SI section \ref{SI-phiSinterp}.

Our skewness ratio does not supplant the variance ratio (or any of its alternatives, such as the Loreau-de Mazancourt
metric), but rather
complements it. Just as one understands more about a distribution by knowing both its variance and skewness, so
one understands more about community dynamics by combining the variance and skewness
ratio approaches. The values of $\phicv$ and $\phi_S$ provide complementary
information about how the distribution of $\xtot$ is influenced by species relationships (see SI section \ref{SI-phiSinterp}).
See Fig. \ref{fig_method_box}A,B for a summary of the two approaches.

## Statistical methods: the link from tail associations to variability

\noindent The skewness ratio $\phiS=\Scom/\Sind$ provides information 
about the influence of species relationships, of all kinds,
on community skewness, $\Scom=\sk(\xtot)$, because $\Scom$ is influenced by species relationships
of all kinds and $\Sind$ is influenced by none. To understand the specific influence of species tail associations
and asymmetries of tail associations on $\Scom$ (G3 of the Introduction), we developed a randomization procedure
that rendered symmetric the tail associations between species, while keeping other statistical properties of 
community dynamics approximately fixed. Given the data $x_i(t)$ ($i=1,\ldots,N$, $t=1,\ldots,T$), the randomization 
procedure produced any desired number $M$ of surrogate datasets $x_i^{(m)}(t)$, $m=1,\ldots,M$, 
with properties detailed in SI section \ref{SI-surrogs} (see also SI Figs \ref{SI-fig:varphiexample}-\ref{SI-fig:check_pairwisecor}). 
Since these surrogate datasets
had symmetric tail associations between species but were otherwise statistically similar to the original dataset,
$x_i(t)$, $\Scom=\sk(\xtot)$ was computed for the original data, and a surrogate skewness value, $\sk(\sum_i x_i^{(m)})$,
was computed for each of the surrogate datasets. 
A test of whether asymmetry of tail associations between species significantly influenced
$\Scom$ was obtained by examining whether $\Scom$ fell sufficiently in the tails of the distribution 
of surrogate values to meet a desired statistical confidence level. 
This approach produces a test of the null hypothesis that $\Scom$ was no more extreme than would have been expected
by chance if species relationships were symmetric in their tail associations, against the alternative hypothesis that 
asymmetries of species relationships contribute meaningfully to community skewness.
We used $M=10000$ to ensure the $p$-values produced by this method were very precise.

The quantity $\Snta$, which represents what community skewness would be without asymmetries of tail association
between species, was defined as the median of the values $\sk(\sum_i x_i^{(m)})$;
it is $\Snta$, and associated quantities which we will now define, that are used to accomplish goal G3 of the Introduction. 
Here and below, "ta" stands for "tail association" and "nta" stands for "no tail association".
We also defined the quantities $\phiSta=\Scom/\Snta$ and $\phiScor=\Snta/\Sind$.
The quantity $\phiScor$ satisfies $\Snta=\phiScor\Sind$, and therefore
is the factor by which correlations between species, as distinct from asymmetries of tail association,
influence community skewness. The subscript "cor" is a reference to the effects of correlation.
The quantity $\phiSta$ satisfies $\Scom=\phiSta\Snta$, and therefore
is the factor by which asymmetric tail associations between species further influence community skewness. 
So a value of $\phiSta$ which differs from $1$ indicates that asymmetric tail associations have a role in
determining $\Scom$. Statistical significance, here, is judged as in the previous paragraph.
It is straightforward to show that $\phiS=\phiSta\phiScor$, so
$\Scom=\phiSta \phiScor\Sind$.
Thus $\phiSta$ and $\phiScor$ separate the influence of 
species relationships on community skewness into factors due to asymmetries of tail associations and correlations
with symmetric tail associations. 
For comparison, $\CVcomsq$ was also computed for the data and for surrogates, 
defining quantities $\CVntasq$, $\phicvcor$ and $\phicvta$
which are related to each other in an analogous way. Fig. \ref{fig_method_box}A-D summarize the new 
quantities and their relationships. Because $\CVcomsq$ depends only 
on species means and covariances, which are preserved by surrogates, we know, *a priori*, that
$\CVntasq \approx \CVcomsq$ and $\phicvta \approx 1$.

Having defined methods, we can now frame hypotheses suggested by Fig. \ref{fig_pedag_cv2_skw} and the 
Theory section using the new methods; tests of the hypotheses will address goal G3 of the Introduction. 
Fig. \ref{fig_pedag_cv2_skw}-\ref{fig_theory} and SI Fig. \ref{SI-fig_theory_SM} 
suggest that communities exhibiting a preponderance of lower-tail association 
should have total community dynamics that are more left-skewed, or less right-skewed, than would have been the case
without these tail associations. This is the same as saying 
that for communities for which $n_L>n_U$ and $\Atot>0$, we should have 
$\Scom<\Snta$. Further, Fig. \ref{fig_pedag_cv2_skw}-\ref{fig_theory} and SI Fig. \ref{SI-fig_theory_SM} suggest that 
communities exhibiting mostly upper-tail association 
should have total community dynamics that are more right-skewed, or less left-skewed, than would have been the case
with symmetric tail associations. This is the same as saying that for 
communities for which $n_L<n_U$ and $\Atot<0$, we should have $\Scom>\Snta$. We test these hypothesis
with the two grassland datasets in Results and Discussion below.

All analyses were done in R. Complete computer code for this project is at https://github.com/sghosh89/Copula_spaceavg. 

```{r pval_skw, echo=F}
hays_skwres<-readRDS("./Results/hays_results/skewness_results/stability_CP_hays.RDS")
pval_skw_hays<-sum(hays_skwres$skw_surrogs<hays_skwres$df_stability$skw_real)/length(hays_skwres$skw_surrogs)
hays_stab<-hays_skwres$df_stability

knz_skwres<-readRDS("./Results/knz_results/skewness_results/stability_CP_knz.RDS")
pval_skw_knz<-sum(knz_skwres$skw_surrogs>knz_skwres$df_stability$skw_real)/length(knz_skwres$skw_surrogs)
knz_stab<-knz_skwres$df_stability
```

# Results and Discussion\label{Results}

<!--***DAN: Shyamolina, I have thought about how I am going to edit this Results section, and I've made a plan, but I need
a few things from you to implement that plan. 

Please compute the fraction of surrogate values on the current Fig. 6F (Dan, this is the new Fig. S5B) that are to the *left* of the black dot. This is the 
p-value for a one-tailed test. It looks like it should be a very small number. Please type this number in a reponse to this comment, 
but *also* please make the number available as an R variable, so when I edit the results I can render the number using an autolink. 
Please say what the R variable name is in a response to this comment. 
***Shya: 0.0027, variable name you can use: pval_skw_hays 

Please compute the fraction of surrogate values on the current Fig. 6H that are to the *right* of the black dot (Dan, this is the current Fig. S5D). This is the 
p-value for a one-tailed test. It looks like it should be a small number. Please type this number in a reponse to this comment, 
but *also* please make the number available as an R variable, so when I edit the results I can render the number using an autolink. 
Please say what the R variable name is in a response to this comment. 

***Shya: 0.0252, variable name you can use: pval_skw_knz 
-->

Addressing the goal G1 of the Introduction, of applying a new way of quantifying species relationships that complements
covariance approaches, the quantities $\cor_l-\cor_u$ for all positively associated 
pairs of species within the Hays and Konza datasets were computed (SI Fig. \ref{SI-fig_CorlmCoru}), as were the
community aggregate quantities $n_L$ and $n_U$ and the total community tail association values $\Atot$ (Methods). 
Results revealed marked differences between Hays and Konza in species' tail associations. 
Although both datasets had pairs of positively associated species with more 
lower-tail association and pairs with more upper-tail association, $n_L=70$ was greater than $n_U=40$ and $\Atot=6.1$ was positive
for Hays, and $n_L=70$ was less than $n_U=82$ and $\Atot=-1.3$ was negative for Konza. Thus 
lower-tail association was dominant at Hays and upper-tail association was dominant at Konza.
<!--***DAN: Note that the numeric results in this para are not auto-linked. This was done because I moved
these values to the text after Shyamolina left the lab, and close to the end of the project, so I wanted to get it done
and did not want to mess around with the autolinks. But if anything changes you'll have to change these values
manually, or else set up autolinks.-->

We earlier framed the hypotheses that communities exhibiting more lower- than upper-tail association 
should also have $\xtot$ more left-skewed, or less right-skewed, than if tail associations were symmetric; 
and communities exhibiting more upper- than lower-tail association 
should have $\xtot$ more right-skewed, or less left-skewed, than if tail associations were symmetric 
(Methods). We tested these hypotheses for Hays and Konza, thereby addressing 
goals G2 and G3 of the Introduction. The quantities of Fig. \ref{fig_method_box} are displayed for our 
data in Fig. \ref{fig_result_box}.
Because $n_L>n_U$ and $\Atot>0$ for Hays (SI Fig. \ref{SI-fig_CorlmCoru}A), $\Scom$ should be less than $\Snta$.
This was confirmed (Fig. \ref{fig_result_box}B). Both $\Scom$ and $\Snta$ were positive and the component 
of the skewness ratio due to tail association, 
$\phiSta=$ `r round(hays_stab$skw_real/hays_stab$skw_ntd_median,1)`, was less than 1. 
Becuse $n_L<n_U$ and $\Atot<0$ for Konza (SI Fig. \ref{SI-fig_CorlmCoru}B), $\Scom$ should be greater than 
$\Snta$. This was also confirmed (Fig. \ref{fig_result_box}D). $\Scom$ and $\Snta$ were again positive,
so for Konza $\phiSta=$ `r round(knz_stab$skw_real/knz_stab$skw_ntd_median,1)` was greater than 1. 
These results were statistically significant: for Hays, $\Scom$ was less than a fraction `r 1-pval_skw_hays`
of <!--the quantities $\sk(\sum_i x_i^{(m)})$--> the analogous quantities computed using the Hays surrogates
datasets (corresponding to $p=1-$ `r 1-pval_skw_hays` $=$ `r pval_skw_hays` for a one-tailed test);
whereas for Konza, $\Scom$ was greater than a fraction `r 1-pval_skw_knz` of the analogous surrogate 
values for Konza ($p=$ `r pval_skw_knz`, one-tailed test; Fig. \ref{SI-fig_result_box}). 
Thus asymmetric tail associations between species 
significantly modified the skewness of the total community abundance, $\xtot$, lowering it for Hays
and raising it for Konza and thereby influencing community variability and probabilities of passing
disaster thresholds. Fig. \ref{fig_empparallel} makes concrete some of the conclusions of this paragraph.

As anticipated (Methods), the variance ratio approach did not register the importance of 
tail associations for community variability and thus provided an incomplete picture of the influence of species
relationships on community variability, i.e., $\CVntasq \approx \CVcomsq$ and $\phicvta \approx 1$ for both 
Hays and Konza (Fig. \ref{fig_result_box}A, C). This result was expected because $\CVcomsq$ 
depends only on species means and covariances, which were preserved 
by the randomized, surrogate datasets from which $\CVntasq$ is computed (Methods). $\CVntasq$ and $\CVcomsq$ were 
not statistically distinguishable for the grassland data: the distribution, 
across surrogates, of the squared coefficients of variation of the quantities 
$\sum_i x_i^{(m)}$ (recall that $\CVntasq$ was defined to be the median of this distribution), was centered very 
close to the value of $\CVcomsq$, for both Hays and Konza (Fig. \ref{SI-fig_result_box}A, C).

We now compare conclusions about Hays that would typically be drawn by using
only the variance ratio approach, to conclusions drawn using both the variance and skewness ratios.
Because $\phicv<1$ (Fig. \ref{fig_result_box}), dynamics in Hays would classically be interpreted as
compensatory, i.e., species relationships buffer 
community variability. However, the combined approach yields more nuanced information.
The skewness ratio, $\phi_S$, was substantially less than $1$, with $\Scom=$ `r round(hays_stab$skw_real,1)` less 
than $\Sind=$ `r round(hays_stab$skw_indep,1)` (Fig. \ref{fig_result_box}),
suggesting that the tendency for the total community abundance $\xtot$ to explode to high values (right skewness)
and surpass high disaster thresholds was mitigated by species relationships. This finding 
supports the variance-ratio-based
conclusion that species relationships are generally stabilizing at Hays, but 
also extends the classic result by revealing 
the greater importance of community relationships for mitigating explosions of $\xtot$ to high values,
compared to their lesser role in mitigating crashes. Apparently, species relationships can mitigate or accentuate
extreme values of $\xtot$ differentially for high versus low extremes.
See SI section \ref{SI-ecdf_approach} for related results
and discussion.

We next compare conclusions about Konza drawn by using
only the variance ratio approach to conclusions drawn using both the variance and skewness ratios.
Because $\phicv$ was slightly less than $1$ (Fig. \ref{fig_result_box}), dynamics in Konza 
would classically be interpreted as slightly compensatory. However, using both 
the variance and skewness ratios together yields conclusions that differ 
substantially. The skewness ratio, $\phiS$, was much greater than $1$, 
with $\Scom=$ `r round(knz_stab$skw_real,1)` correspondingly much greater than $\Sind=$ 
`r round(knz_stab$skw_indep,1)` (Fig. \ref{fig_result_box}). 
So although species relationships slightly reduce the variance of $\xtot$, they also dramatically increase right skew
of $\xtot$. This in turn leads to an overall greater propensity for explosions of $\xtot$ to high values, and 
a correspondingly accentuated probability that $\xtot$ will exceed high disaster thresholds. Crashes of $\xtot$
to low values are reduced by species relationships. The variance ratio obscures the fact that species relationships
modify extreme behavior of $\xtot$ differently in the two tails of its distribution. Although the variance 
ratio suggests that species dynamics at Konza are compensatory, reducing community 
variability, in fact this is true only for mitigating crashes of $\xtot$ to low values; species relationships 
instead accentuate the tendency of $\xtot$ to explode to high values. 
See SI section \ref{SI-ecdf_approach} for related results and discussion. 

Thus the skewness ratio newly reveals differences between Hays and Konza.
For both communities, $\phicv$ is slightly or moderately less than $1$,
suggesting compensatory, stabilizing species relationships. However, $\phi_S$ was less than 1 for Hays 
and was markedly greater than 1 for Konza due partly to differing influences of interspecific tail 
associations on the distributions of $\xtot$. Hays can probably still be labeled as a community showing
stabilizing, compensatory dynamics. In contrast,
species relationships in Konza only reduce community variability in that they mitigate crashes of $\xtot$. 
Species relationships simultaneously accentuate explosions of $\xtot$ to high values, a feature that should probably 
not be considered as stabilizing. 
For some applications, it may indeed be more important to avoid crashes than to avoid explosions
of $\xtot$ to high values, and for those applications it may still be reasonable to characterize species
relationships at Konza as stabilizing; but this possibility depends on the applied context, and does not
reduce the importance of our larger point that species relationships can be differently "stabilizing"
or "destabilizing" (i.e., reducing versus increasing community variability) in the two tails 
of $\xtot$, a feature of empirical reality not captured
by the variance ratio and other standard approaches.

These results highlight the previously unrecognized importance of tail associations for community dynamics. 
The influence of tail associations on the distribution of $\xtot$ can be
as great or greater than the influence of species correlations with symmetric tail association. For Konza, 
$\phiSta$ and $\phiScor$ were similar (Fig. \ref{fig_result_box}), indicating that both tail associations and
correlations were comparable in their effects on the skewness
of $\xtot$ (recall that $\Scom=\phiSta \phiScor \Sind$). For Hays, however, $\phiScor$
was close to 1 and $\phiSta$ differed substantially from 1 (Fig. \ref{fig_result_box}), indicating that essentially all of the 
influence of species relationships on $\sk(\xtot)$ was due to asymmetric tail associations between species. 
The influence of tail associations on the CV of $\xtot$ was negligible, as indicated 
previously, but that reflects
the inability of the CV to detect the impact of tail associations, rather than indicating that
tail associations were unimportant. 

Fig. \ref{fig_empparallel} makes more concrete some of the above conclusions through comparison of $\xtot$ time series
in Hays and Konza to total abundance time series computed for randomized, surrogate datasets for which tail associations
were rendered symmetric. The figure also helps illuminate future
applied possibilities of the ideas of this study. For Hays, $\xtot$ would have risen to more extreme high values
were it not for the dominance of lower-tail associations between species in that system. For Konza, $\xtot$
would have fallen to more extreme low values were it not for the dominance of upper-tail assocations 
in that system. Because grazers seem more likely to be harmed by extreme low values of $\xtot$ than by extreme high
values, the differences between these systems illuminate some of the
applied possibilities of our ideas. Whereas we emphasize that our data were from 
livestock-excluded plots, so applications to grazing are only suggestive, 
our overall results nevertheless show that details of species interactions can mitigate or
accentuate extreme values of aggregate ecosystem functioning differently for high versus low extremes.
Since human systems that depend on ecosystem functions (e.g., livestock grazing, though there are many examples)
will react differently to high versus low extremes (and this will depend on the ecosystem function in question),
the details of tail associations between species should be of clear applied interest in many circumstances. In essence, 
we showed that tail associations between species can have important effects on extreme behavior of ecosystem 
functioning variables; so there is applied importance whenever extreme values would be 
a concern.

As described in the Introduction, the variance ratio $\phicv$ compares $\CVcomsq$ to the independence benchmark
$\CVindsq$ ($\phicv=\CVcomsq/\CVindsq$), whereas the Loreau-de Mazancourt metric 
$\phiLdM$ compares $\CVcomsq$ to the benchmark value $\CVsynsq$
that community $\text{CV}^2$ would take if species dynamics were perfectly linearly associated ($\phiLdM=\CVcomsq/\CVsynsq$). This
choice leads to substantial interpretive and other advantages of $\phiLdM$ over $\phicv$ [@loreau2008species; @Thibaut2013],
and as a result $\phiLdM$ is increasingly widely used.
An approach to community skewness using a paradigm similar to that of Loreau and de Mazancourt is available, and reveals the 
same main ideas we elaborated in this study using $\phiS$, but does not appear to have the same advantages over $\phiS$ that
$\phiLdM$ has over $\phicv$; nevertheless, such a metric may be worth considering in future work. 
<!--An approach to community skewness using a paradigm similar to that of Loreau and de Mazancourt is available-->
Because $\phiLdM=\CVcomsq/\CVsynsq$, whereas $\phicv=\CVcomsq/\CVindsq$, 
an alternative to $\phiS=\Scom/\Sind$ which follows the same paradigm may be $\phiSLdM=\Scom/\Ssyn$, where 
$\Ssyn$ is the benchmark value that community skewness would take in the case of perfect synchrony between species
dynamics. Whereas $\CVsynsq$ was community $\text{CV}^2$ computed under an assumption of perfect linear associations 
between the dynamics of different species [@Thibaut2013], under our framework it may make more sense to consider $\Ssyn$ to be
community skewness under an assumption of species comonotonicity. Perfect linear associations between species cannot 
generally be realized without altering the marginal distributions of the species, and such alteration seems, to us, to be
inconsistent with using a benchmark scenario in which species relationships are altered but other aspects of dynamics are 
kept constant. An alternative approach using $\phiSLdM=\Scom/\Ssyn$ would 
reveal essentially the same main ideas that we have revealed with
the current approach of this study, because $\Snta$ and $\CVntasq$ would probably have 
to be defined in the same way. The main ideas of this study center on the importance of tail associations for
the skewness of $\xtot$, and therefore for extreme values of $\xtot$. Those ideas can be explored regardless of whether 
independence or perfect synchrony is used as a benchmark for comparison.
<!--but does not appear to have the same advantages over $\phiS$ that $\phiLdM$ has over $\phicv$-->
One of the advantages
of $\phiLdM$ over $\phicv$ is that $\phiLdM$ always takes values between $0$ and $1$, facilitating comparisons 
across communities. But the same is not necessarily true for $\phiSLdM$, as a careful examination of Figs
\ref{fig_theory} and \ref{SI-fig_theory_SM} will reveal ($\sk(\xtot)$ for scenario 2 divided by the same quantity for 
scenario C is much greater than $1$, whereas computing the analogous quotient for scenario 3 gives a quantity much less
than $-1$).
Another advantage of $\phiLdM$ over $\phicv$ is that $\sqrt{\CVsynsq}$ equals a weighted average of the CVs of
individual species, and hence can be interpreted as population-level variability. Thus $\phiLdM$ can be elegantly interpreted
as a scaling factor converting population- to community-level variability, whereas an analogous interpretation for 
$\phicv$ is not available. But such an interpretive advantage does not appear to extend to $\phiSLdM$. 
Decades passed between the formulations of the original variance ratio and the 
Loreau-de Mazancourt and other improvements. Our goal with this study was to introduce the 
main ideas of tail association and its importance for community stability. We welcome the possibility that future research
will reveal a metric that captures these main ideas, but with additional advantages.

<!--New version-->
Our work may also provide lessons for other fields. More detailed statistical 
descriptions of synchrony have become increasingly useful in recent years. 
Tail association metrics and related "copula" statistics, though apparently 
seldom used so far in most biological 
research areas, may turn out to provide fruitful additional ways to 
make synchrony statistics more detailed and more useful across many fields.
In many fields, techniques for quantifying the
synchrony of signals have historically been crude, 
but have recently become more detailed, thereby enhancing progress. 
For instance, one paper in this special issue develops newly detailed 
frequency-specific descriptions of a kind of synchrony, and then uses them
to make new inferences about musical and rhythmic 
behavior in humans [@DotovInReview].
The examples of [@Reus_etal_InPress] on the interactive vocalizations of primates 
are also enhanced by a frequency-specific perspective: Fig. 1C of @Reus_etal_InPress shows coordination
of not only the timings but also the pitches of the calls of a pair of indris, a feature which would not
be revealed without the additional statistical detail provided by the frequency-specific techniques 
of that paper (spectrograms).
Sheppard, Zhao and colleagues [@sheppard2016; @sheppard2019; @zhao2020] developed detailed frequency-specific measures 
of synchrony in population ecology which provided major inferential benefits. 
One common feature of all these studies is that progress was achieved through use of 
frequency-specific statistical
descriptors of synchrony that were more detailed than earlier correlation and related basic tools. 
But spectral tools are not the only statistical approach by which additional 
information can be extracted from synchronous or asynchronous pairs of signals. 
Mathematically, the nature of the relationship between two quantities is an
infinite-dimensional object, even when neither quantity exhibits spectral structure 
of any kind [@Ghosh_copula]. Most synchrony measures extract only a small portion of the
available information. But the field of copula statistics [@nelsen2006_copula]
provides well-developed and general tools for assessing the complete information content. Tail 
association tools such as used in this study assess one aspect of copula structure. 
The potential of tail association and copula tools to reveal previously unsuspected
useful information in patterns of synchrony across multiple fields of biology may be great,
as demonstrated by our successes in this study. For instance, @LenseInReview indicates, in this special issue, 
that many childhood neurodevelopmental disorders are associated with deficits in 
rhythm, timing and synchrony. One wonders if diagnosis of disorders
may be improved by applying copula statistics to data on synchrony
tasks attempted by patients. Whereas we are not child psychologists and cannot
critique work in that field or in many of the other disparate fields represented in this special issue,
the power of copula statistics and the general scarcity of prior application of such tools to 
synchrony compel us to recommend these approaches to synchrony researchers in all fields.

<!--Examples of synchrony and aggregation in many fields-->
<!--In many hierarchical dynamical systems,
synchrony between multiple signals is more important than the individual signals themselves, because
what really matters to the next level of functioning in the system is the total value, or another kind of 
aggregated value, of the signals. And the aggregated value is crucially impacted by the nature of 
the synchrony of the signals. This is
true of both financial and population systems. 
But it seems likely to be a common feature of many systems. 
For instance, a neuron may fire only when its input neurons fire simultaneously, so that total
inputs pass a threshold. 
Or the electrical grid may crash only when the demands of a multitude of users become synchronized,
producing total usage spikes.--> 
<!--Examples of where improvements in stats has led to improved results. Refram so that looking 
at frequency specificity is a clear way that more detailed stats have already benefitted us-->

# Conclusion \label{Discussion}

It has been known for decades that species relationships can accentuate or mitigate the 
temporal variability of an aggregate community property such as the total abundance of 
all species, through synchronous or compensatory dynamics of species. However, our
results show that consequences of species relationships for community stability are more 
complex than previously recognized. Synchronous relationships between 
the dynamics of different species are not necessarily well characterized by standard, commonly 
used covariance measures: species can exhibit contrasting patterns of tail association which 
are not detected by standard approaches. Moreover, the temporal stability of 
total abundance is also not necessarily well characterized by commonly used measures 
based on the coefficient of variation. Total abundance can exhibit patterns of skewness 
which correspond to increased probabilities of crashes, or explosions to high values, 
that are readily interpretable as instability. Thus the original 
variance ratio approach, whereby patterns of species’ synchronous or compensatory dynamics 
are assessed using covariance tools and then related to community variability measured using 
a coefficient of variation, is probably inadequate, by itself. Species relationships can actually 
mitigate or accentuate community variability differentially for community crashes or explosions 
to high values, a fact which was apparently not previously recognized. 
In other words, species relationships can mitigate explosions of 
total abundance to high values, while simultaneously exacerbating the tendency 
to crash to low values; or vice versa. Our new skewness ratio approach complements the variance ratio 
to reveal a more multidimensional nature of species relationships and their effects on community variability. 
Overall, our work has begun a new and more flexible approach to synchrony, 
compensatory dynamics, and their influence on community stability.

# Acknowledgments
\noindent The authors thank Jim Bever, Bryan Foster, Lauren Hallett and Michel Loreau for helpful discussions and 
constructive criticism. The authors were partly supported 
by U.S. National Science Foundation grants 1714195 and 2023474, the James S. McDonnell Foundation, and the 
California Department of Fish and Wildlife Delta Science Program. 

<!--Pedagogical figure-->
\begin{figure}[!h]
\includegraphics[width=18cm]{./Results/pedagog_figs/pedagog_cv2_skw.pdf}
\caption{Pedagogical example showing two communities (A, B) of 20 species each that are not distinguishable 
via the classic variance ratio approach and related approaches, but that nonetheless differ markedly in 
potentially important ways. See Introduction for details and interpretation. See 
SI Fig. \ref{SI-fig_pedag_cv2_skw_SM} for larger versions of panels A, B.
Dashed and dotted lines on C-F represent thresholds. If the lower threshold (dashed line) were
a disastrous threshold for $\xtot(t)$ to cross, then the first scenario (C, E) would involve much
more risk than the second scenario (see probabilities of crossing the threshold listed to the left
of the dashed lines on E, F - the probability for scenario 1 is an order of magnitude higher).
On the other hand, if the upper threshold (dotted line) were a disastrous threshold for 
$\xtot(t)$, the second scenario would involve more risk (see probabilities to the right of the 
dotted lines on E, F - the probability for scenario 2 is much higher).
Time series were generated and statistics computed for 100000 years, but are plotted on A, B
for 60 years only, for clarity. \label{fig_pedag_cv2_skw}}
\end{figure} 
<!--\caption{Pedagogical example showing two communities (A, B) of twenty species each that are not distinguishable 
via the classic variance ratio approach and related approaches, but that nonetheless differ markedly in 
potentially important ways. Species marginal distributions, and therefore species means and variances
over time, were the same in both scenarios, as were species covariances, up to sampling variation. 
Thus the quantities $\CVcomsq$, $\CVindsq$ and $\phi_{\text{cv}}$ of the classic variance ratio approach
were essentially the same in the two scenarios (A, B), as was the Loreau-de Mazancourt (Loreau \emph{et al.} 2008)
variance ratio, denoted $\phi_{\text{LdM}}$. Nevertheless, differences between the scenarios in interspecific relationships, not 
captured by covariances, resulted in differently skewed total time series, $\xtot(t)$, in the two scenarios
(C-F). The total was more likely to crash to low values in the first scenario (C, E)
than in the second scenario (D, F; compare the probabilities on E and F listed to the left of the dashed line), 
whereas $\xtot(t)$ was more likely to explode to high values 
in the second scenario than in the first (probabilities on E, F to the right of the dotted line). 
Dashed and dotted lines show lower and upper thresholds which may be considered
disastrous for $\xtot(t)$ to cross. The probability of crossing each threshold was negligible for one scenario, but
non-negligible for the other.
In both scenarios, there were two groups of species (red and black lines on A, B; see SI 
Fig. \ref{SI-fig_pedag_cv2_skw_SM} for a larger version).
Species in different groups were compensatory, but for scenario 1, within-group synchrony occurred principally
when populations were low, and for scenario 2 when populations were high. See the Introduction and 
SI section \ref{SI-PedagogFigMethods} for notation and 
additional details. Time series were generated and statistics computed for 100000 years, but are plotted on A, B
for 60 years only, for clarity. \label{fig_pedag_cv2_skw}}
-->

<!--Theory figure-->
\begin{figure}
\includegraphics[width=.5\textwidth]{./Results/TheoryFig_MT.pdf}
\caption{Theoretical scenarios demonstrating how species relationships can influence a 
community aggregate
quantity in ways that are not detectable by classic approaches. The abbreviation ``sd'' stands 
for ``standard deviation'', ``sk'' stands for ``skewness'', ``th'' stands
for ``threshold'', and ``fr'' stands for ``frequency''. Scenario C is the reference, comonotonic scenario of
maximal community variability. See the Theory section for other notation and 
interpretations, and SI section \ref{SI-TheoryMethods} and Fig. \ref{SI-fig_theory_SM} 
for scenarios 3 and 4 and for mathematical details.\label{fig_theory}}
\end{figure}

<!--Methods figure for asymmetric tail association and partial Spearman correlation-->
\begin{figure}[!h]
\includegraphics[width=4 cm]{./Results/pedagog_figs/LTUT_rho_0p8_LTdep.pdf}
\caption{A normalized rank plot (see Methods) can visually help reveal that 
an association between variables is stronger in the lower tails of the variables' distributions than it is in the upper tails
(as in this example), or \emph{vice versa} (as in SI Fig. \ref{SI-fig_pedag_taildep_SM}). In this example, points are clustered
closer to the $u_i=u_j$ diagonal in the lower left than in the upper right, revealing stronger 
association between the variables in the lower tails of their distributions. 
The partial Spearman correlation, $\cor_{b_l,b_u}(x_i,x_j)$ (Methods), within a band defined by two bounds $0 \leq b_l < b_u \leq 1$
can be computed for any band to quantify the component of the overall Spearman correlation
due to that band. Solid diagonal lines show the band associated with $b_l=0.1$ and $b_u=0.25$.
Up to sampling variation, the (total) Spearman correlation is 0.8 on both this plot and 
SI Fig. \ref{SI-fig_pedag_taildep_SM}, though patterns of
tail association differ. This figure was slightly modified from Fig. 4 of (30)
and Fig. 7 of (31). \label{fig_pedag_taildep}}
\end{figure}
<!--***DAN: This is a reminder that the citations in the above caption are not autolinks, so need to be
changed manually if anything changes. They are supposed to be the extinction risk paper and then BIVAN. Do
not delete this reminder.-->

<!--Pedagogical figure to explain method-->
\begin{figure}[!h]
\includegraphics[width=6.5in]{./Results/pedagog_figs/method_fig_MT.pdf}\\
\caption{Summary of approach. Arrows in A-D imply an equation
where the quantity in the box at the arrow tail, times the quantity labelled on 
the arrow itself, equals the quantity in the box at the arrow head. For instance, the 
arrow in A represents the equation $\CVcomsq=\phicv \CVindsq$. \label{fig_method_box}}
\end{figure}

<!--Skewness results for both datasets-->
\begin{figure}[!h]
\includegraphics[width=14cm]{./Results/pedagog_figs/Result_box.pdf}
\caption{The quantities of Fig. \ref{fig_method_box} computed for the Hays (A, B) and Konza (C, D)
datasets. See SI Fig. \ref{SI-fig_result_box} for related results. \label{fig_result_box}}
\end{figure}

<!--Figure that is parallel in some ways to Fig. 1, but uses data-->
\begin{figure}[!h]
\includegraphics[width=4.5in]{./Results/EmpiricalParallelOfFig1.pdf}
\caption{Example, using data, of the effect of tail association on dynamics of aggregate community abundance,
$\xtot$. Emprical $\xtot$ time series for Hays (A) and Konza (D) are solid lines, with distributions of values 
across time shown in B, E, respectively. Dashed lines on A, D show community 
abundance sums for ``typical'' surrogate datasets produced via a randomization procedure that rendered tail associations
symmetric, but left unchanged other statistical aspects of the data (Methods). Distributions across time for the dashed
lines are in C, F, respectively. Thus comparing dashed to solid lines, or comparing B to C (for Hays) or 
E to F (for Konza) reveals the effects of tail associations between species on community dynamics. Though the 
effect appears subtle, extreme high values are mitigated
by tail associations in Hays (e.g., values on C exceed 8000, but values on B do not, and see also the 
dotted-line threshold on A); whereas extreme low
values are mitigated in Konza (e.g., values on F go just below 120 but values on E do not, and see also the dotted-line
threshold on D). If extreme values spell disaster in an applied scenario, these differences are important. This figure is an 
empirical parallel to aspects of Fig. 1. "Abundance" was basal cover (cm$^2$) for Hays (A-C) and percent cover
for Konza (D-F). Vertical ticks
between axes and histograms show the location of individual data points.
Skewness, ``sk'', is displayed on each histogram panel. Surrogates were selected for display which were median
with respect to skewness of the total abundance, so that values of skewness on C and F approximate $\Snta$.
\label{fig_empparallel}}
\end{figure}

\clearpage

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