Paper.Rmd

---
title: "Copulas and their potential for ecology"
author: "Shyamolina Ghosh, Lawrence W. Sheppard, Mark T. Holder, Terrance D. Loecke, Philip C. Reid, James D. Bever, Daniel C. Reuman"
date: ""
fontsize: 12pt
geometry: "left=2.54cm,right=2.54cm,top=2.54cm,bottom=2.54cm"

output: 
  pdf_document:
    number_sections: yes
    keep_tex: yes
    fig_caption: yes
    includes:
      in_header: head_maintext.sty
      
tables: True
link-citations: True
urlcolor : blue
indent : True

csl: ecology-letters.csl
bibliography: REF_ALL.bib
---

```{r setup_Paper, echo=F}
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 specified digit
#seed<-101
#source("mtime.R") #A function needed for caching
```

\raggedright
\setlength\parindent{2em}
\setlength{\parskip}{6pt}

\noindent \emph{Affiliations:}

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

\noindent Holder: Department of Ecology and Evolutionary Biology and Biodiversity Institute, University of Kansas, Lawrence, KS, 66045, USA

\noindent Loecke: Environmental Studies Program and Kansas Biological Survey, University of Kansas, Lawrence, KS, 66047, USA

\noindent Reid: School of Biological & Marine Sciences, University of Plymouth, Drake Circus, Plymouth PL4 8AA, UK; Continuous Plankton Recorder Survey, The Marine Biological Association, The Laboratory, Citadel Hill, Plymouth PL1 2PB, UK

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

<!--\noindent \emph{Possible alternative titles:}

Complex and informative dependencies propagate throughout ecology and can be revealed using copulas

Copulas reveal complex and informative dependencies propagating throughout ecology

Copulas and their importance for ecology

Copulas and their importance in ecology

Copulas and their importance throughout ecology

Copulas and their potential throughout ecology

Copulas and their importance

Copulas in ecology

Copulas

The importance of copulas

The importance of copulas in ecology-->

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

\noindent \emph{The number of words in the abstract:} 241

\newpage

\noindent \textbf{Abstract}

\noindent All branches of ecology study relationships among and between environmental and biological variables. However, 
standard approaches to studying such relationships, based on correlation and regression, 
provide only some of the complex information contained in the relationships. Other 
statistical approaches exist that provide a complete description of relationships between variables, 
based on the concept of the *copula*; they are applied in finance, neuroscience and elsewhere, 
but rarely in ecology. We explore the concepts that underpin copulas and the potential 
for those concepts to improve our understanding of ecology. We find that informative copula 
structure in dependencies between variables is common across all the environmental, 
species-trait, phenological, population, community, and ecosystem functioning datasets we 
considered. Many datasets exhibited asymmetric tail associations, whereby two variables were more 
strongly related in their left compared to right tails, or *vice versa*. We describe mechanisms by which 
observed copula structure and tail associations can arise in ecological data, including a Moran-like 
effect whereby dependence structures are inherited from environmental variables; and asymmetric or nonlinear influences of environments on ecological variables, such as 
under Liebig's law of the minimum. We also describe consequences of copula structure for 
ecological phenomena, including impacts 
on extinction risk, Taylor's law, and the temporal stability of ecosystem services. By documenting the 
importance of a complete description of dependence between variables, advancing conceptual frameworks, and 
demonstrating a powerful approach, we encourage widespread use of copulas in ecology, which we
believe can benefit the discipline.

\noindent \emph{Keywords:} dependence, regression, correlation, copula, population, community, ecosystem functioning

\newpage

# Introduction\label{Intro}

All branches of ecology 
study relationships among biological variables and relationships between 
environmental and biological 
variables. However, commonly used correlation and regression 
approaches to studying 
such relationships are limited, and 
provide only a partial description of the relationship. 
For instance, datasets showing different 
relationships may have the same correlation coefficient
(Fig. \ref{fig_cop_pedag1}A-B). The variables of Fig. \ref{fig_cop_pedag1}A 
(respectively, Fig. \ref{fig_cop_pedag1}B) are principally 
related in the left (respectively, right) portions of their distributions, 
an asymmetric pattern of association that can
have ecological significance, as discussed below, but that is not captured by 
correlation. Although correlation is only one way to 
study relationships between variables, other common metrics also 
provide only partial information.

\begin{figure}[!h]
\begin{center}
\includegraphics[width=3.5in]{./Results/PedagogFig1.pdf}
\caption{(A-D) Bivariate datasets showing diverse relationships between 
variables. Some of these datasets nevertheless have the same Pearson (P)
or Spearman (S) correlation coefficients. (C) and (D) show normalized ranks of (A), (B), 
respectively. (E, F) Reflections
of (C), (D), respectively, about a vertical axis. Subscripts h and v 
in the axis labels stand for "horizontal" and "vertical".\label{fig_cop_pedag1}}
\end{center}
\end{figure}

Well-developed approaches do exist, however, and are 
applied widely in other fields, 
that provide a description of the association 
between two or more variables that is mathematically provably 
complete [@nelsen2006_copula; @joe2014_dependence; @MaiScherer2017]. 
The approaches are based on 
the concept of the *copula*. Copula approaches separate information in a bivariate 
random variable into two parts: the information in the 
marginal distributions (which says nothing about the association between the 
variables), and the remaining information (which is solely about
the association). We here introduce some simple concepts, based on ranks,
that relate to copulas and give a conceptual 
flavor of copula ideas to be introduced 
formally in the "Background" section below (section \ref{Background}). Given 
a sample $(x_i,y_i)$, $i=1,2,\ldots,n$ (e.g., 
Fig. \ref{fig_cop_pedag1}A-B), information about the structure of the association
between $x$ and $y$ can be separated from information contained 
in the marginal distributions by 
considering the plot of $u_i$ versus $v_i$, where
$u_i$ is the rank of $x_i$ in the set of $x_j$ ($j=1,\ldots,n$), divided by
$n+1$; and $v_i$ is defined in the same way but using the $y_i$. The rank
of the smallest element in a set is $1$. The $u_i$ and 
$v_i$ are called *normalized ranks* of the $x_i$ and $y_i$; they relate to 
the empirical cumulative distribution function of the $x_i$ and $y_i$,
respectively. We henceforth abbreviate "cumulative distribution function" as "cdf". 
See, for instance,
Fig. \ref{fig_cop_pedag1}C-D, which show the normalized ranks of 
A-B. Ranking makes the marginal distributions of 
the component datasets uniform, isolating the 
dependence structure. Dependence structure and marginals can then be studied 
separately. Normalized rank plots such as Fig. \ref{fig_cop_pedag1}C-D
relate to copulas in a way to be described 
(section \ref{Background}). Note that Pearson correlation, even though it is the most commonly
used measure of association, is modified by monotonic
transformation of the component variables (Appendix 
\ref{SM-pedagog1_details}, Fig. \ref{SM-fig_transf_pearson_spearman}),
and therefore reflects not only dependence information, but also
information on the marginals [@Genest2007]. 
Spearman and Kendall correlations, being rank-based, 
are not influenced by monotonic transformations of the variables. 
They provide information solely on the dependence between those variables.

A main benefit of a copula approach is 
that it can detect associations in the tails of distributions, and asymmetries of such associations. 
Tail association (introduced formally in section 
\ref{Background}) is association between extreme values of variables. 
If smaller values of two positively associated variables are more 
strongly associated than are larger values, the variables are said to
have stronger lower- or left-tail association than right- or upper-tail association
(Fig. \ref{fig_cop_pedag1}C);
and *vice versa* if larger values are more strongly associated 
(Fig. \ref{fig_cop_pedag1}D). Datasets of the same 
Spearman or Kendall correlation can have a range of tail associations 
(Fig. \ref{SM-fig_cop_pedag}), which can be quantified (see sections on Background and
methods for Q1 below, sections \ref{Background} and \ref{Methods}). 
Although copulas can be used to model the entire
complex dependence structure of data, we have found tail associations to
be a useful component of that structure, so tail association 
is a main focus of this paper. We give four examples, below, which
illustrate some of the ecological meaning of tail associations.

The goal of this paper is to explore the potential for 
applications of copulas in ecology and to 
estimate to what extent ecological understanding may benefit from
using copulas. We provide an introduction
to copula concepts in section \ref{Background}. 
Readers can find additional abbreviated [@Anderson2018; @Genest2007] 
and comprehensive [@nelsen2006_copula; @joe2014_dependence; @MaiScherer2017]
introductions elsewhere. Copulas have already been used effectively 
in a few studies to illuminate ecological phenomena 
[@Valpine2014; @Anderson2018; @Popovic2019], but such usage is
still rare. Whereas @Anderson2018 provide a specific method, based
on copulas, for important problems of the 
analysis of multivariate population count data,
we instead advocate more broadly for the application of copulas across ecology.
Our results suggest that environmental, 
ecological and evolutionary processes 
may commonly generate complex dependence structures, including asymmetric tail 
associations, and that greater use of copulas can illuminate underlying processes.
We believe copula approaches are among the tools all ecologists 
should be considering for analysis of their data in the 21st century.

\begin{figure}[!h]
\begin{center}
\includegraphics[width=7in]{./MasterFigure.pdf}
\caption{Summary and guide for analyses presented in the text. Middle boxes
correspond to ecological datasets for which Q1 was examined (in sections
\ref{Data}-\ref{Results}, see also Table \ref{tab_data_info}). Upper boxes 
correspond to potential causes of non-normal (non-N) copula structure that
were examined (in sections \ref{Causes}-\ref{Results_Q2} of the text). Lower boxes 
correspond to potential consequences we examined of non-normal copula structure for
ecological understanding and for applications (in sections 
\ref{Consequences}-\ref{Results_Q3} of the text). Arrow labels (A-D for analyses pertaining 
to causes, Q2; and X-Z for analyses pertaining to consequences, Q3) 
link to locations in the text.\label{MasterFigure}}
\end{center}
\end{figure}

One way we may expect, *a priori*, a study of tail associations and 
copula dependence structures 
to better reflect ecological relationships relates to Liebig's law of the
minimum. Liebig's law is 
the idea that growth is controlled not by total resources but by the resource which
is scarcest relative to organism needs. 
If, for instance, the growth of a plant depends on soil 
nitrogen, N, and other factors, a plot of growth rates versus soil N may 
look more like Fig. \ref{fig_cop_pedag1}A or C than like 
Fig. \ref{fig_cop_pedag1}B or D, i.e., the two variables may show 
left-tail association: N 
controls plant growth, producing a clear
relationship, only when it is limiting. 
This aspect of the association would be visible in a rank-based
plot, but may not be revealed by correlation. 

A second reason why we may expect detailed understanding of dependence structures
to benefit ecological
research is that prior work demonstrated, using copulas, complex 
structure in the spatial dependence of environmental 
variables [@Serinaldi2008; @She2018; @Goswami2018;
@li2013]. An environmental variable 
measured through time in two locations may show strong tail associations between
the locations if
intense meteorological events are also widespread, as seems 
frequently to be the case: extreme values are associated with intense events,
so happen in both places at the same time, whereas moderate values of the 
environmental variable may be associated with local phenomena, which differ
between the locations.
Spatial dependence in an environmental variable 
tends to beget spatial dependence 
in fluctuations of populations or other ecosystem variables 
(this is called \emph{spatial synchrony}) influenced by the environmental 
variable. This is called the \emph{Moran effect}. 
If asymmetric tail associations, or other complex dependence structure, is transmitted 
from environmental to ecological variables, then we would expect complex dependence 
structure and tail associations to be a common feature of the 
spatial synchrony of population, community, 
biogeochemical and other environmentally influenced ecological variables.
Synchrony attracts major interest in 
ecology [@Liebhold2004; @Sheppard2016; @walteretal2017].  

A third reason why we may expect an understanding of complex 
dependence structures and tail associations to benefit 
ecology is that such an approach may help illuminate 
causal mechanisms between variables. If two species,
Sp1 and Sp2, are strong competitors, abundances of the two species
across quadrats should be negatively related, as in 
Fig. \ref{fig_cop_pedag1}E-F. If Sp1 is plotted on the horizontal axis
and Sp2 on the vertical, then Fig. \ref{fig_cop_pedag1}E 
may suggest Sp1 is the dominant 
competitor: when Sp1 is abundant, Sp2 is necessarily rare because it is suppressed; 
whereas when Sp1 is rare, Sp2
may be abundant, or may also be relatively rare due to 
limiting factors other than Sp1. 
Alternatively, Fig. \ref{fig_cop_pedag1}F may suggest Sp2
is the dominant competitor. Other causal hypotheses may produce 
similar dependence structure, and it is usually impossible to obtain
complete information on causal pathways from analyses which are
fundamentally based on associations. However, 
examination of tail associations may
be combined with biological reasoning to rule out
some causal hypotheses which could not previously be eliminated. 
@Popovic2019 used copulas to illuminate causal relationships
between species. 

A fourth reason why we may expect an understanding of tail associations and complex 
dependence structures to be useful
for ecology has to do with spatially aggregated or averaged quantities. 
Many ecological variables of applied importance 
depend on the spatial average or sum of local 
quantities. For instance, regional methane and CO$_2$ 
fluxes are the sum of local fluxes; 
and the total economic value of a fishery is the sum of local catches. 
We will explore how tail associations between fluctuations of local
quantities can influence fluctuations of 
the spatial mean or sum, and how this may influence 
higher organizational levels in ecology and human concerns.
To illustrate this idea we cite @li2013, who demonstrated that the overall 
reliability of wind-generated electricity depends sensitively on details
of the dependence between wind speeds at multiple generator sites. 
Spatially aggregated ecological variables may be subject to a similar 
effect. For instance, if populations of a 
pest species in different locations are all positively associated and are also more strongly
related to each other in their right tails, then local outbreaks will tend 
to occur together, creating regional epidemics. 
Stronger left-tail associations in a pest species, even if overall 
correlation were the same, would have more benign effects. 

We approached our overall goal of exploring the potential of copula approaches 
for ecology by addressing the following specific questions. (Q1) Do datasets
in ecology have dependence structure distinct from that of a multivariate 
Gaussian/normal
distribution (here called non-normal copula structure)? Do positively associated 
ecological variables show tail 
associations distinct from those of a normal distribution, and in particular
do they show asymmetric tail associations? Normal copula structure is assumed
by standard approaches that use multivariate 
normal distributions or distributions obtained by transforming
the marginals of a normal distribution. So Q1 asks whether ecological
data contain dependency information distinct from a standard or default case. 
(Q2) If the answer to
Q1 is "yes", then what are some possible causes/mechanisms of non-normal
copula structure and asymmetric tail associations in ecology? And (Q3), what are the 
consequences of non-normal copula
structure and asymmetric tail associations for ecological understanding and
applications? After section \ref{Background}, 'Background on copulas', we address Q1 
(sections \ref{Data} through \ref{Results} of the paper) 
by analyzing several datasets, including 
environmental, species-trait, phenological, population, community, and 
ecosystem functioning data, selected to span multiple 
sub-fields and organizational levels of ecology. 
We address Q2 (sections \ref{Causes} through \ref{Results_Q2}) using 
simple models. We address Q3 (sections \ref{Consequences} through \ref{Results_Q3}) using 
both data and models. Section \ref{Discussion} is the Discussion.
Multiple analyses are brought to bear on each question, and these
are summarized diagrammatically in Fig. \ref{MasterFigure}, which can 
serve as a *post hoc* guide to the paper. 
Copula approaches have been used to great effect in neuroscience [@onken2009], 
bioinformatics [@Kim2008], medical research [@Emura2016], direct study of environmental 
variables [@Serinaldi2008; @She2018; @Goswami2018; @li2013], and finance [@Li2000],
and have also been used effectively, but rarely so far, in ecology [@Valpine2014; 
@Anderson2018; @Popovic2019]. 
We argue that benefits of wider use of
copula descriptions of dependence in ecology will be substantial.

# Background on copulas\label{Background}

We give a brief introduction to copulas, focusing, for simplicity, on bivariate 
copulas and on concepts needed
for the rest of the paper. We assume familiarity with the basic 
language of probability distributions and random variables.
See @nelsen2006_copula, @Genest2007, @joe2014_dependence, and
@MaiScherer2017 for general background on copulas, 
and see also @Anderson2018, part of 
which is a short introduction to 
copulas for ecologists. A bivariate *copula* can be defined as a bivariate cdf 
with both margins uniform on $(0,1)$ (@joe2014_dependence, p. 7). 
This will be the 
cdf of a bivariate random variable with uniform marginals, and 
the terminology *copula* is sometimes also applied to this random vector. 
A foundational theorem of
@sklar1959_theorem (see also @MaiScherer2017, p. 16) 
says that if $F$ is the cdf
of a random vector $(X,Y)$, with margins $F_X$ and $F_Y$, then there 
exists a copula $C$
such that for all $(x,y)$ in the Euclidean plane, $F(x,y)=C(F_X(x),F_Y(y))$.
The theorem also says $C$ is unique if $F_X$ and $F_Y$ are continuous.
Thus $C$ couples the bivariate cdf, $F$, with the 
cdfs of the marginals. 
Finally, Sklar's theorem says that if $D$ is any bivariate copula and 
$G_A$ and $G_B$ are univariate cdfs, then 
$D(G_A(x),G_B(y))$ is a valid cdf 
of some bivariate random variable.
The applications of this study fall within the case in which univariate marginal 
distributions have continuous, strictly 
monotonic cdfs, and this case is simpler. So we 
henceforth make such an assumption. Table \ref{SM-tab_summary_notation} provides a summary
of notation.
<!--***DAN: Shyamolina, once I email you an edited version of the
table you produced, please make a table in the sup mat Rmd and cause it to be
cited here in place of the SX above. Please place the table in such
a way that SX is replaced with S1,  since I believe this is now going
to be the first supplementary table we are citing.-->

Sklar's theorem implies that
any random vector can be constructed from a unique copula and marginal 
cdfs, and, furthermore,
any copula and any univariate cdfs give rise to a random vector;
so Sklar's theorem provides a correspondence between random vectors
and pairs consisting of a copula and two univariate cdfs,
thereby making it possible to study a random vector by studying these two
constituents.
Marginals contain no information about the
dependence structure of a random vector, so the copula contains
all such information - it is a complete description of the dependence
between variables.
<!--***-->
The univariate cdfs associated with a random vector 
$(X,Y)$ are its margins, $F_X$ and $F_Y$, and the 
associated copula is the copula, $C$, guaranteed by Sklar's theorem
(and guaranteed unique, thanks to the assumptions of continuity 
and monotonicity made above). In fact, and making Sklar's correspondence
more concrete, we show in Appendix \ref{SM-ContStrictMonoton}
that $C$ is the bivariate distribution
function of the random variable $(F_X(X),F_Y(Y))$. 
<!--***-->
Conversely, given a copula
$D$ and univariate cdfs $G_A$ and $G_B$,
Sklar's theorem guarantees that $D(G_A(x),G_B(y))$ is a 
cdf, so the random variable with this 
cdf is the random variable corresponding
to $D$, $G_A$ and $G_B$ under Sklar's correspondence. 
Again making Sklar's correspondence
concrete, we show in Appendix \ref{SM-ContStrictMonoton}
that this random variable is $(G_A^{-1}(U),G_B^{-1}(V))$, 
where $G_A^{-1}$ and $G_B^{-1}$ are the inverses of $G_A$ and $G_B$
and $(U,V)$ is the random variable with cdf $D$
(which has uniform marginals).
If we conflate a copula with the random vector of 
which the copula is the cdf, then Sklar's correspondence
is between random variables, $(U,V)$, with uniform marginals 
and random variables, $(X,Y)$, with arbitrary continuous, strictly monotonic
marginals. The correspondence is simply via the application of the
univariate cdfs, or their inverses: $U=F_X(X)$,
$V=F_Y(Y)$, $X=F^{-1}_X(U)$, $Y=F^{-1}_Y(V)$.

Thus Sklar's theorem makes it possible to construct
bivariate distributions in two separate steps: by specifying the marginal 
distributions and by specifying the dependence structure.
Construction of several bivariate distributions is carried out in this way in
Fig. \ref{PedagogFigCopThry}. The contrast between Fig. \ref{PedagogFigCopThry}C and D
illustrates in particular how a bivariate distribution can be changed, while retaining 
the same marginals, by changing the copula. Fig. \ref{PedagogFigCopThry}C is a 
bivariate normal distribution, but Fig. \ref{PedagogFigCopThry}D clearly 
is not, although its marginals
are the same. Fig. \ref{PedagogFigCopThry}D shows stronger association of the two variables
in the left than in the right tails. Sklar's theorem
also makes it possible to study and model dependence separately 
from marginal distributions. For instance, in a pedagogical 
example, @Anderson2018 modelled the dependence between abundance
data for two fish species, red moki and black angelfish, by first 
modelling the marginal distributions and then modelling 
the dependence using Gaussian copulas (their section 2).
See the Discussion for a few words on multivariate copulas, which can
be developed in the same way as above [@nelsen2006_copula; @joe2014_dependence;
@MaiScherer2017; @Anderson2018].

\begin{figure}[!h]
\begin{center}
\includegraphics[width=6.5in]{./Results/PedagogFigCopThry.pdf}
\caption{Normal (A) and Clayton (B) copulas were combined with standard normal 
marginals (C, D) and gamma marginals (E, F) via Sklar's theorem to produce 
bivariate distributions. Note that a normal copula is distinct from a normal marginal 
distribution: a normal copula is the copula of a bivariate normal distribution.
See section \ref{Background} for more information on normal, Clayton,
and other copulas. Each copula was used with both sets of marginals and each
set of marginals was used with both copulas, to demonstrate that: both copula and 
marginals contribute fundamentally to the resulting distribution, the copula 
contributing the information on the dependence between the variables; and that 
one can select the copula and marginals independently. Bivariate
distributions, including the copulas, are depicted via their log-scale probability density functions (pdfs), 
and marginals are depicted via their linear-scale pdfs. Grey dots are
50 random samples from each distribution. The parameter $0.7$ was used for the normal 
copula, the parameter $2$ for the Clayton 
copula, and shape and scale parameters $5$ and $1$ for the gamma marginals. 
Variables $u$ and $v$ were used for copulas
and $x$ and $y$ for distributions created by pairing a copula with
(non-uniform) marginals.
\label{PedagogFigCopThry}}
\end{center}
\end{figure}

We introduce some standard copula families, as examples, and also
because we will fit these families to data. The bivariate normal
copula family, which consists of copulas of bivariate normal
distributions, is a one-parameter family. The parameter, $p$,
corresponds to the correlation
of the related bivariate normal distribution, and controls the degree of association
between the variables. A normal copula with $p=0.7$ was already introduced (Fig.
\ref{PedagogFigCopThry}A). The formula for the
copula is $\Phi_{2,p}(\Phi_{1}^{-1}(u),\Phi_{1}^{-1}(v))$, where $\Phi_1^{-1}$ is the inverse of the
cdf of a (univariate) standard 
normal and $\Phi_{2,p}$ is the cdf of a
bivariate normal distribution with mean $(0,0)$ and covariance matrix having $1$s on the diagonal
and $p$ in the off diagonal positions. Formulas for all the 
copula families we use were provided by @Brechmann2013. 
However, formulas for copulas are often not as 
directly informative about copula properties as probability density function (pdf) 
pictures; Fig. \ref{figped_N} has 
example pdfs corresponding to bivariate normal copulas for various values of $p$.
Note that the pdfs are symmetric
around the diagonal line $v=-u+1$ (Fig. \ref{figped_N}), so normal copulas have symmetric 
associations between the two variables in the left and right tails. 
The Clayton copula family,
of which we already pictured an example (Fig.
\ref{PedagogFigCopThry}B), is another one-parameter family. In contrast to normal copulas, Clayton copulas
have stronger left- than right-tail association. The formula is 
$\left[  \reumax(u^{-p}+v^{-p}-1,0)\right]^{-1/p}$, for parameter $p$; though this again provides 
less direct intuition than do example pdfs (Fig. \ref{figped_C}). Higher values of the parameter,
$p$ produce copulas with higher Kendall or Spearman correlation. The survival
Clayton family is the symmetric opposite of the Clayton family, showing
stronger right- than left-tail association (Fig. \ref{SM-figped_SC}). The BB1 copula 
family is a two-parameter family, thereby providing more flexibility: for some parameters
BB1 copulas have stronger left- than right-tail association, and for others the reverse
(Fig. \ref{figped_BB1}). See @nelsen2006_copula, @joe2014_dependence and 
@Brechmann2013 for further
information on these and other families.

Having introduced tail association conceptually and referred to it
in examples, we now need a precise definition of a measure of tail association.
We use the term *tail association* to describe the general idea of the strength of
association between two variables in the tails of their distributions. One 
way this is measured is called *tail dependence* (@joe2014_dependence, section 2.13; 
@nelsen2006_copula, section 5.4), defined here. Given a random vector $(X,Y)$
with margins $F_X$, $F_Y$ and defining $U=F_X(X)$, $V=F_Y(Y)$, the upper-tail
dependence of $X$ and $Y$ is defined as 
$\lambda_U = \lim_{u \rightarrow 1^-} P[Y>F_Y^{-1}(u) | X>F_X^{-1}(u)]$. This equals
$\lim_{u \rightarrow 1^-} P[U>u | V>u]$, which in turn equals
$\lim_{u \rightarrow 1^-} P[U>u, V>u]/P[V>u]=\lim_{u \rightarrow 1^-} P[U>u, V>u]/(1-u)$,
which shows upper-tail dependence is defined symmetrically in the two variables.
All the variables we consider were positively associated when they were significantly
associated, so we think of positively associated variables when conceptualizing the
definitions here. Lower tail dependence is defined analogously, as
$\lambda_L = \lim_{u \rightarrow 0^+} P[Y \leq F_Y^{-1}(u) | X \leq F_X^{-1}(u)]$.
Tail dependence is a property of $(U,V)$, so depends only
on the copula and not on the marginals combined with it. It can be shown 
(@nelsen2006_copula, section 5.4) that $\lambda_L=\lambda_U=0$ for the normal copula,
$\lambda_L=2^{-1/p}$ and $\lambda_U=0$ for a Clayton copula with parameter $p>0$,
and $\lambda_L=0$ and $\lambda_U=2^{-1/p}$ for a survival Clayton  copula with parameter $p>0$.
In section \ref{Methods} we will introduce
ways tail dependence can be applied to data, by first fitting copulas to the data. We will also 
introduce nonparametric measures of tail association.

Plots using normalized ranks were proposed in the Introduction as being conceptually
similar to copulas; we explain the connection. As stated above, the copula associated
with a random variable $(X,Y)$ is the cdf of the random variable
$(F_X(X),F_Y(Y))$. Given a sample $(x_i,y_i)$, $i=1,\ldots,n$ from $(X,Y)$, let $\hat{F}_X$ and
$\hat{F}_Y$ be the empirical cdfs associated with the $x_i$ and $y_i$, respectively.
These are step functions of $x$ and $y$, respectively, that start at $0$ for low values of 
$x$ and $y$ and jump by $1/n$
at each of the data points. Therefore $\hat{F}_X(x_j)$ is the rank of $x_j$ in
the set $\{x_i : i=1,\ldots,n\}$, divided by $n$, 
and analogously $\hat{F}_Y(y_j)$ is the rank of $y_j$ in the set 
$\{y_i : i=1,\dots,n\}$, divided
by $n$. For large $n$, $\hat{F}_X(x_i)$ and $\hat{F}_Y(y_i)$ 
approximate the cdfs $F_X$ and $F_Y$, respectively. 
But $\hat{F}_X(x_i)$ equals the normalized rank of $x_i$ times
$(n+1)/n$, and likewise $\hat{F}_Y(y_i)$ equals the normalized rank of $y_i$ times
$(n+1)/n$. So the normalized rank in turn approximates the empirical cdf and therefore 
the cdf. Thus the normalized rank pairs $(u_i,v_i)$ for $i=1,\ldots,n$ can be 
regarded as approximate samples from the random variable $(F_X(X),F_Y(Y))$, 
which is the random variable associated with the copula of $(X,Y)$. The scatterplot 
of $v_i$ versus $u_i$ can be used to infer aspects of copula 
structure (see section \ref{Methods} below).

\begin{figure}[!h]
\begin{center}
\includegraphics[width=\textwidth]{./Results/PedagogSuppMat_Normal.pdf}
\caption[Example normal copulas]{Log-transformed pdfs (A-E) and samples (F-J) from example normal copulas. 
K is Kendall correlation; $p$ is the value of the parameter for the normal
family (it is a 
one-parameter family); and LT and UT are the measures of lower- and upper-tail dependence, respectively (for the definitions of these, see the section on Background on copulas, 
section \ref{Background}). The 
parameter range for the family is $p \in [-1,1]$, and lower- and upper-tail
dependence are always $0$.\label{figped_N}}
\end{center}
\end{figure}

\begin{figure}[!h]
\begin{center}
\includegraphics[width=\textwidth]{./Results/PedagogSuppMat_Clayton.pdf}
\caption[Example Clayton copulas]{Log-transformed pdfs (A-E) and samples (F-J) from example Clayton  
copulas. K is Kendall correlation; $p$ is the value of the parameter for the Clayton family (it 
is a one-parameter family); and LT and UT are the measures of lower- and 
upper-tail dependence, 
respectively (for the definitions of these, see the section on Background on copulas, 
section \ref{Background}). The parameter range for the family 
is $p \in (0,\infty)$, lower-tail dependence is $2^{-1/p}$ and upper-tail
dependence is $0$.\label{figped_C}}
\end{center}
\end{figure}

\begin{figure}[!h]
\begin{center}
\includegraphics[width=\textwidth]{./Results/PedagogSuppMat_BB1.pdf}
\caption[Example BB1 copulas]{Log-transformed pdfs for example BB1 copulas. K is Kendall correlation; $p_1$ and $p_2$
denote the two parameters of the family (it is a 
two-parameter family); and LT and UT are our measures of lower- and upper-tail dependence, 
respectively (for the definitions of these, see the section on Background on copulas, 
section \ref{Background}). The parameter ranges for 
the family are $p_1 \in (0,\infty)$ and $p_2 \in [1,\infty)$, 
lower-tail dependence is $2^{-1/(p_1 p_2)}$ and upper-tail
dependence is $2-2^{1/p_2}$.\label{figped_BB1}}
\end{center}
\end{figure}

# Data \label{Data} 

```{r read_somedata, echo=F, results='hide'}
d_soilCN<-read.csv("./Data/SoilCN/SoilCNdata.csv")
rownames(d_soilCN)<-d_soilCN$X
d_soilCN<-d_soilCN[,-1]
d_soilCN<-na.omit(d_soilCN)
d_birdsBMR<-readRDS("./Data/BMR/BirdBodyMassesMetabolicRates/my_birdsBMR_data.RDS")
d_mammalsBMR<-readRDS("./Data/BMR/MammalianBodyMassesMetabolicRates/my_mammalsBMR_data.RDS")
d_cc<-readRDS("./Results/Ceder_creek_results/cop_HB_all_yr.RDS")
```

The datasets we used included environmental, species-trait, phenological, 
population, community, and ecosystem functioning data (Table \ref{tab_data_info}), 
selected to span multiple fields and levels
of organization within ecology. Copula structure 
between atmospheric weather variables such as rainfall or wind speed, 
as measured in multiple locations through time, has been 
examined previously in the meteorological literature [e.g., @Serinaldi2008; @li2013]. We therefore 
examined environmental variables from the soil instead, using the Rapid Carbon 
Assessment database [RaCA; @Wills2014].
<!--***DAN: Add Wills and Loecke pending reference if published in time--> 
The database comprises measurements of soil organic carbon and total soil 
nitrogen (megagrams C or N per hectare of soil surface) from `r nrow(d_soilCN)` 
locations across the coterminous United States (Fig. \ref{SM-fig_SoilCNmapUSA} for locations).
Species-trait data were average species basal metabolic rate (BMR, KJ per hour) and 
body mass (grams) for `r nrow(d_birdsBMR)`  species of birds [@Mcnab_2009] 
and `r nrow(d_mammalsBMR)`  species of mammals [@Mcnab_2008].
These data have been much studied, but to our knowledge the copula structure 
of the association has not been examined. 
Species-trait data such as these reflect the coevolution of the two traits.
Phenological data were aphid first-flight dates from 10 locations (Fig. \ref{SM-fig_aphid_map}) 
across the United Kingdom (UK) for 20 aphid species (Table \ref{SM-tab_aphid_info}), for the 
35 years 1976 to 2010. These time series were computed from the Rothamsted Insect 
Survey (RIS) suction-trap dataset [@Harrington2014; @Bell2015]. The first of our two 
population-level datasets was also derived from the RIS suction-trap data, and comprised 
total counts of aphids trapped for the same locations, species, and years. 
Our second population dataset comprised average annual plankton abundance estimates 
for 14 locations (Fig. \ref{SM-fig_plankton_map})
in the North Sea and British seas for 22 taxa (Table \ref{SM-tab_plankton_info}) for the 
56 years 1958 to 2013. The locations are $2^\circ$ by $2^\circ$ grid cells.
These data were computed from the Continuous Plankton Recorder survey of the 
UK Marine Biological Association.
Community-level data, obtained from the Cedar Creek Ecosystem Science Reserve, were 
plant aboveground biomass [@ccdataple120] and Shannon's diversity index (computed from plant 
species percent cover data [@ccdatapce120]) for `r nrow(d_cc[[5]])` plots [@ccdataplotse120], 
each 9m by 9m, as sampled in the years 1996-2000 and 2007, each year 
analyzed separately [@Tilman2001; @Tilman2006]. 
Finally, ecosystem functioning data were methane (CH$_{4}$) fluxes between the soil or 
water surface and the troposphere, measured at 13 locations at daily to weekly intervals 
from September 2015 to September 2016 at the Great Miami Wetland mitigation bank in Trotwood, 
Ohio [@Holland1999; @Jarecke2016; @smyth2019using]. Each included location was measured on at least 50 dates. See Appendix \ref{SM-Data_details} for additional data details. 

Our environmental, species-trait, and community datasets happen to be 
"bivariate" datasets in the sense that they comprise two quantities measured at 
different locations or for different species (Table \ref{tab_data_info}). 
Our phenology, population, and ecosystem functioning datasets are 
"multivariate" in that they comprise, for each taxon, measurements through 
time at multiple locations; 
copula structure was studied for each location pair.

<!--summary table with info of all data used in Paper-->
```{r tab_summary_data_info, echo=F, results="asis",message=F}
library(kableExtra)
library(dplyr)
tab_data_info<-as.data.frame(matrix(NA,nrow=8,ncol=4))
colnames(tab_data_info)<-c("Data","Category","No. of measurements we used", "References")
tab_data_info$Data<-c("Soil C and N","Bird body masses and BMR", "Mammal body masses and BMR", 
                      "Green spruce and other aphid species abundances", "Leaf-curling plum and other aphid first flight dates","\\textit{Ceratium furca} and other plankton taxa abundances", "Plant diversity and aboveground biomass","Methane-flux")
tab_data_info$Category<-c("Environmental","Species-trait","Species-trait",
                          "Population","Phenological","Population","Community level","Ecosystem functioning")
tab_data_info$`No. of measurements we used`<-c(paste(nrow(d_soilCN)," locations",sep=""),
                                               paste(nrow(d_birdsBMR)," birds",sep=""),
                                               paste(nrow(d_mammalsBMR)," mammals",sep=""),
                                               "10 locations with at least 30-yr. timeseries for each of 20 sp.",
                                               "10 locations with at least 30-yr. timeseries for each of 20 sp.",
                                               "14 locations with at least 45-yr. timeseries for each of 22 taxa",
                                               paste(nrow(d_cc[[5]])," plots",sep=""),
                                               "13 locations with at least 50 dates of data")
tab_data_info$References<-c("US Rapid Carbon Assessment database (RaCA), Wills et al. (2014)", "McNab (2009)","McNab (2008)",
                            "Rothamsted Insect Survey","Rothamsted Insect Survey","Continuous Plankton Recorder Survey",
                            "Cedar Creek Ecosystem Science Reserve, biodiversity experiment e120",
                            "Great Miami Wetland Mitigation Bank, Smyth et al. (2019)")

knitr::kable(tab_data_info, 
             #format="latex", align="l",
             format="latex", align="l",escape=F,
             caption = "Summary table for the data we used. Bold entries are the multivariate data, the rest are bivariate datasets (see section \\ref{Data}). Basal metabolic rate=BMR. \\label{tab_data_info}", 
             booktabs = T,linesep = "\\addlinespace") %>%
             column_spec(1:4,width="3.4 cm")%>%
             row_spec(c(4,5,6,8), bold=T)
```


```{r tab_summary_copula_info, echo=F, results="asis",message=F}
library(kableExtra)
library(dplyr)
tab_copula_info<-as.data.frame(matrix(NA,nrow=16,ncol=4))
colnames(tab_copula_info)<-c("Copula family","No. of parameters","Tails with stronger association", "Reference")

tab_copula_info$`Copula family`<-c("Normal",
                                   "Clayton","Survival Clayton",
                                   "Gumbel","Survival Gumbel",
                                   "Joe","Survival Joe",
                                   "Frank", 
                                   "BB1","Survival BB1",
                                   "BB6","Survival BB6", 
                                   "BB7","Survival BB7", 
                                   "BB8","Survival BB8")

tab_copula_info$`No. of parameters`<-c("1",
                                       "1","1",
                                       "1","1",
                                       "1","1",
                                       "1",
                                       "2","2",
                                       "2","2",
                                       "2","2",
                                       "2","2")

tab_copula_info$`Tails with stronger association`<-c("Neither (symmetric)",
                                             "Lower","Upper",
                                             "Upper","Lower",
                                             "Upper","Lower",
                                             "Neither (symmetric)",
                                             "Either (depending on params)","Either (depending on params)",
                                             "Upper","Lower",
                                             "Either (depending on params)","Either (depending on params)",
                                             "Upper or symmetric (depending on paras)","Lower or symmetric (depending on params)")
  
# caution: this fig. numbers are not auto-linked
tab_copula_info$Reference<-c("Fig. 4",
                              "Fig. 5","Fig. S3",
                              "Fig. S7","Fig. S8",
                              "Fig. S9","Fig. S10",
                              "Fig. S11",
                              "Fig. 6","Fig. S12",
                              "Fig. S13","Fig. S14",
                              "Fig. S15","Fig. S16",
                              "Fig. S17","Fig. S18")

knitr::kable(tab_copula_info, 
             format="latex", align="l",escape=F,
             caption = "Summary table for the sixteen copula families we used (see section \\ref{Background}). For the parameters we consider, all copula families model positive association between variables. \\label{tab_copula_info}", 
             booktabs = T,linesep = "\\addlinespace") %>%
             column_spec(2,width="2 cm")
```


# Concepts and methods for Q1 \label{Methods} 

Recall that Q1 has three parts: Do datasets
in ecology have non-normal copula structure? Do positively associated 
ecological variables show tail 
associations distinct from those of a normal distribution? And in particular
do they show asymmetric tail associations? 
We addressed Q1 first via a model selection procedure in which several
families of copulas were fitted to 
our ecological datasets and fits were compared via the Akaike and Bayesian 
Information Criteria (AIC and BIC). One of the fitted copulas was the normal
copula, making possible comparisons of the degree to which the normal versus
other copula families were good descriptions of data. 

For bivariate datasets 
$(x_i,y_i)$ for $i=1,...,n$, model selection involved several steps. 
First, we produced normalized ranks 
$u_i$ and $v_i$ as in the Introduction. 
Second, we tested the independence of the $u_i$ and $v_i$
using the statistic $\sqrt{(9n(n-1))/(4n+10)}|\tau|$, where $\tau$
is Kendall's tau for the data. @Genest2007 argue that this statistic is approximately standard normally
distributed. We used the 
implementation of this test in `BiCopIndTest` in the `VineCopula` package
in R. We tested for independence because our model selection 
algorithms were ineffective if data could not 
be distinguished from independent data, 
since many copula families include the independent copula. 
If independence could be rejected ($0.05$ level), model 
selection proceeded. Third, we fit 16 bivariate copula families (see below) to the 
normalized ranks via maximum likelihood. The approach of fitting
copula families to normalized ranks was recommended 
by @Genest1995 and @Shih1995. Their estimator of copula
parameters, which we use, is consistent, asymptotically normal,
and fully efficient at indepencence [@Genest1995]. @Genest2007
recommend carrying out inferences of dependence structures (which was our
goal here) using normalized ranks.
We used the fitting implementation given in `BiCopEst` in `VineCopula`.
Fourth, we obtained AIC and BIC values and accompanying model weights 
$\text{AIC}_{\text{w}}$ and 
$\text{BIC}_{\text{w}}$ [@burnham2003_modelselection] for each fitted copula. 
`BiCopEst` also 
provided lower- and 
upper-tail dependence of the best-fitting member of each family. 
$\text{AIC}_{\text{w}}$ values were used to get
model-averaged lower- and upper-tail dependence values using standard model averaging 
formulas [@burnham2003_modelselection]; likewise for BIC. 

Thus for each bivariate dataset, 
the end products of our procedure were threefold, corresponding to the
three parts of Q1 listed in the Introduction: A) the AIC, BIC, $\text{AIC}_{\text{w}}$ and 
$\text{BIC}_{\text{w}}$ values for each of our 16 copula families (see below for list),
including the normal family, providing an inference as to whether the normal
copula or an alternative was better supported by 
data; B) lower- and upper-tail dependence measures for each fitted copula and 
model-averaged tail dependence measures, providing information on whether,
and to what extent, tail dependence differed from the tail dependence of 
a normal copula (i.e., 0); and C) the difference of lower- and upper-tail
dependence measures for each fitted copula and model averages of those 
quantities, providing information on whether tail
dependence was asymmetric.

Model selection methods give the *relative* support of several models, but do not 
indicate whether any of the models are an objectively good fit. To test that, 
we tested the goodness of fit of our AIC-best copula family using a 
bootstrapping procedure of @Wang2000 and @genest2006goodness,
implemented as `BiCopGofTest` in `VineCopula`. The procedure performed 
one test based on a Cramer-von Mises statistic and another based on a 
Kolmogorov-Smirnov statistic. To keep computation times reasonable, 
a run using $100$ bootstraps was performed. If
the $p$-value from either test was $<0.2$, tests were re-run with 
$1000$ bootstraps.

We fit 16 bivariate copula families, exhibiting a variety of 
lower- and upper-tail dependence characteristics, with bivariate datasets. The purpose
of using a large collection of families was to include a 
variety of alternative dependence structures to have a robust model
selection and multi-model inference procedure. For that purpose, it 
is not important for the reader to
understand the details of these copulas,
and, additionally, these copulas have been 
described in detail elsewhere [@joe1997multivariate; @Brechmann2013]. So 
we only briefly identify each family, say a few words about its tail 
dependence, and provide pictorial descriptions (Figs 
\ref{figped_N}-\ref{figped_BB1}, \ref{SM-figped_SC} and \ref{SM-figped_G}-\ref{SM-figped_SBB8}). 
The pictorial descriptions can also be used to aid quick comparisons of copula
families when alternative families are being considered for future applications. 
We used several families that can exhibit positive lower-tail dependence 
(of strength depending on parameters) and 
zero upper-tail dependence: the Clayton, survival Gumbel, survival Joe
and survival BB6 families. We used several families that can exhibit positive
upper-tail dependence (strength depends on parameters) and zero lower-tail
dependence: the survival Clayton, Gumbel, Joe and BB6 families. We used
families that show zero upper- and lower-tail dependence: the normal
and Frank families. These families have pdfs symmetric about the line
$v=-u+1$. We also used several families that can show both upper- and
lower-tail dependence, in relative amounts depending on parameter
values: the BB1, survival BB1, BB7 and survival BB7 families. 
We also used the BB8 family, which shows zero lower-tail dependence,
and zero upper-tail dependence except for a boundary case for the parameters.
And we used the survival BB8 copula, which shows zero upper-tail dependence, 
and zero lower-tail dependence except for a boundary case.
"Survival" families are rotations of the copula with the similar name by
180 degrees. 
See @joe1997multivariate and @Brechmann2013 for details on all families. We used the implementations
provided in the `VineCopula` package for the R programming
language. See Figs \ref{figped_N}-\ref{figped_BB1}, \ref{SM-figped_SC} and \ref{SM-figped_G}-\ref{SM-figped_SBB8} 
for visual depictions of the pdfs of these copulas and how
pdfs and tail dependence are influenced by the parameters. See 
Table \ref{tab_copula_info} for a summary.

<!--***DAN: Info on families kept here for convenience, see also the file
CopulaFamilyParamRangesAndTailDependence.jpg, stored in the BIVAN folder:
lower not upper dependence:
Clayton (3), survival Gumbel (14), survival Joe (16),
SBB6 (18)
zero both, and symmetric:
  normal (1), Frank (5)
upper not lower dependence:
  SC (13), Gumbel (4), Joe (6)
  BB6 (8)
both, in relative amounts that depend on parameters:
  BB1 (7) - maybe need to look at formulas for this one to understand what the   
    limits are, since I only chose some random values of parameters in my pdf 
    pictures
  SBB1 (17)
  BB7 (9), SBB7 (19)
BB8 (10) - 0 LT dep, UT dep 0 except in a boundary case for the parameters
SBB8 (20) - reverse of BB8-->

For multivariate datasets, we performed the bivariate analysis
for all possible pairwise combinations of distinct locations. We carried out
pairwise bivariate analyses instead of trying to fit a multivariate copula,
for simplicity and because that approach was
sufficient to answer our research questions; but see the Discussion
for a few words on multivariate copulas. 
We present
the number of pairs for which a non-normal copula
was the AIC-best copula, and we characterize AIC differences between
the normal and AIC-best copulas across location pairs. 
We also computed the model-averaged lower- and
upper-tail statistics, and differences
between these, for each pair of locations, and we characterize
the distributions of these values across location pairs.  

In addition to our model selection approach, we also used a
nonparametric approach, to provide greater 
confidence in our answers to Q1. We used three 
statistics which quantify the extent to which 
the normalized ranks $u_i$ and $v_i$ are 
related in any part of their distributions. We here describe the 
statistics, with additional details in Appendix \ref{SM-nonparam_stats}.
The statistics are defined with positively associated
variables in mind. All our variables were positively associated
when they were significantly associated.
Given two bounds 
$0 \leq l_b < u_b \leq 1$, define the lines $u+v=2l_b$ and 
$u+v=2u_b$, which intersect the unit square (Fig. \ref{NonparamStatsFig}). 
Our statistics quantify the association
between $u_i$ and $v_i$ in the region bounded by these lines.
Using $l_b=0$ and $u_b\leq 0.5$ gives information about 
association in the left parts of the distributions of $u$ and $v$, and 
using $u_b=1$ and $l_b \geq 0.5$ gives information about association in 
the right parts. The first statistic, a partial Spearman correlation,
\begin{equation}
\cor_{l_b,u_b}(u,v)=\frac{\sum 
   (u_i-\mean(u)) (v_i-\mean(v))}{(n-1)\sqrt{\var(u)\var(v)}},\label{eq:partialspearman}
\end{equation}
\noindent is the portion of the Spearman 
correlation of $u_i$ and $v_i$ that is attributable to the points 
between the bounds. 
Here sample means and variances are computed using all $n$ data 
points, but the
sum is over only the indices $i$ for which $u_i+v_i > 2l_b$ and 
$u_i+v_i < 2u_b$.
Larger values of the partial Spearman correlation indicate stronger positive association.
The sum of 
$\cor_{0,0.5}(u,v)$ and $\cor_{0.5,1}(u,v)$ (or some other choice of 
$\cor_{l_{b_k},u_{b_k}}(u,v)$ for bounds $l_{b_k},u_{b_k}$ that partition the
interval $(0,1)$) equals the standard Spearman correlation
as long as no points fall exactly on the bounds.
We also defined a statistic $\Ps_{l_b,u_b}$ (Appendix \ref{SM-nonparam_stats}), 
which has a similar interpretation to $\cor_{l_{b},u_{b}}$.
Our third statistic, $\Dsq_{l_b,u_b}$, is the average squared 
distance between points satisfying $u_i+v_i>2l_b$ and 
$u_i+v_i<2u_b$ and the line $v=u$.
Unlike $\cor_{l_b,u_b}$ and 
$\Ps_{l_b,u_b}$, for which large values indicate strong association
between the bounds, small values of 
$\Dsq_{l_b,u_b}$ indicate strong association. These statistics are not estimators of the 
tail dependence quantities
defined previously, but rather are conceptually similar measures of associations in
the tail portions of the distributions when appropriate values of $l_b$ and $u_b$ are
selected. 

\begin{figure}[!h]
\begin{center}
\includegraphics[width=3.5in]{./Results/NonparamStatsFig.pdf}
\caption{The partial Spearman correlation, $\cor_{l_{b},u_{b}}(u,v)$, 
within a band can be computed
for any band (section \ref{Methods}) to describe how the strength of association
between $u$ and $v$
varies from one part of the two distributions to another, 
as can the statistics $\Ps_{l_b,u_b}(u,v)$ and $\Dsq_{l_b,u_b}(u,v)$. Diagonal lines
show two bands, the data in the left/lower band showing stronger
association than those in the right/upper band.\label{NonparamStatsFig}}
\end{center}
\end{figure}

For large datasets (large $n$),
we used $l_b$ and $u_b$ close together without incurring undue
sampling variation in our statistics, and we considered 
multiple bands $(l_b,u_b)$ to understand how association
varies across different parts of the distributions. But for datasets with 
smaller $n$ we considered only $l_b=0$, $u_b=0.5$ and $l_b=0.5$, $u_b=1$,
abbreviating $\cor_l=\cor_{0,0.5}$ ($l$ is for "lower") 
and $\cor_u=\cor_{0.5,1}$ ($u$ is for "upper"). Likewise $\Ps_l=\Ps_{0,0.5}$,
$\Ps_u=\Ps_{0.5,1}$, $\Dsq_l=\Dsq_{0,0.5}$, $\Dsq_u=\Dsq_{0.5,1}$.
To test for asymmetry of association in upper and lower portions of distributions,
we used differences $\cor_l-\cor_u$, $\Ps_l-\Ps_u$ and $\Dsq_u-\Dsq_l$
for smaller datasets (note the opposite order
in the last of these); and $\cor_{0,u_b}-\cor_{1-u_b,1}$, 
$\Ps_{0,u_b}-\Ps_{1-u_b,1}$ and
$\Dsq_{1-u_b,1}-\Dsq_{0,u_b}$ with $u_b$ close 
to $0$ for large datasets.
We tested our statistics in Appendix \ref{SM-Test_npa_stat} (see also 
Figs \ref{SM-fig_stat_testing35}-\ref{SM-fig_stat_testing1000}). 

For bivariate datasets, we compared values 
of $\cor_{l_b,u_b}$,
$\Ps_{l_b,u_b}$, $\Dsq_{l_b,u_b}$, $\cor_{0,u_b}-\cor_{1-u_b,1}$, 
$\Ps_{0,u_b}-\Ps_{1-u_b,1}$ and
$\Dsq_{1-u_b,1}-\Dsq_{0,u_b}$ to 
distributions of the same statistics computed on surrogate datasets that were produced from 
the empirical data by randomizing it in a special way to have no 
tail dependence (Appendix \ref{SM-surrog_test}). 
The surrogate/randomized datasets had exactly 
the same marginal distributions as the empirical data and had very similar Kendall or Spearman 
correlation (the surrogate algorithm had two versions, one for preserving each correlation coefficient), 
but had normal copula structure. Thus our comparisons tested the null hypothesis 
that our statistics took values on the empirical data no different from what would have been expected 
if the copula structure of the data were normal, but the data were otherwise 
statistically unchanged. The comparison of one of the
statistics $\cor_{l_b,u_b}$, $\Ps_{l_b,u_b}$, $\Dsq_{l_b,u_b}$, as computed on the 
empirical data, to the distribution
of its values computed on surrogate datasets 
provides a test of whether association between the variables in the 
part of the distributions specified 
by $l_b$ and $u_b$ is different from
what would be expected from a null hypothesis of normal copula structure. 
Thus the comparison addresses the first two parts of Q1
for bivariate datasets. 
Significant deviations correspond to deviations from normal 
copula structure. In particular, deviations 
using $l_b=0$ and $u_b$ small (say, $0.1$) correspond to lower-tail associations 
different from that of a normal copula;
likewise, using $l_b=0.9$ and $u_b=1$ tests for upper-tail associations different 
from a normal copula.
The comparison of one of the statistics 
$\cor_{0,u_b}-\cor_{1-u_b,1}$, 
$\Ps_{0,u_b}-\Ps_{1-u_b,1}$ and
$\Dsq_{1-u_b,1}-\Dsq_{0,u_b}$,
as computed on the empirical data, to the distribution of its values on surrogate datasets 
provides a test of asymmetry of tail associations. Thus this comparison helps address the
third part of Q1. 

For multivariate datasets, we calculated Spearman and Kendall correlations
and $\cor_l$, $\cor_u$, $\Ps_l$, $\Ps_u$, $\Dsq_l$, $\Dsq_u$,
$\cor_l-\cor_u$, $\Ps_l-\Ps_u$, and $\Dsq_u-\Dsq_l$
for all pairs of sampling locations, and we then computed means across the pairs
for each statistic. We used a spatial resampling scheme 
(Appendix \ref{SM-spatial_resamp}) to calculate confidence intervals of these means.
The scheme is identical to that proposed by @bjornstadfalck2001 and also 
used by @walteretal2017. Code and data for this
project are archived at \url{www.github.com/sghosh89/BIVAN}.

# Results for Q1: Ecological datasets have non-normal copula structure and asymmetric tail associations \label{Results}

```{r read_BivMS_soilCN_data,echo=F}
d<-read.csv("./Data/SoilCN/SoilCNdata.csv")
rownames(d)<-d$X
d<-d[,-1]
d<-d[,c("SOCstock100","TSNstock100")] #raw soilCN data
BivMS_res_Loecke<-readRDS("./Results/fitting_results/BivMS_soilCN.RDS")
```

```{r read_RES_aphid_count, echo=F}
RES_aphid_count<-readRDS("./Results/fitting_results/AphidCopulaFit_selecloc_count_species_10.RDS")
#------------ also reading some extra stuffs -------
stat_aphid_count<-readRDS("./Results/stat_results/stat_aphid_count/stat_aphid_count_sp_10.RDS")
stat_aphid_ff<-readRDS("./Results/stat_results/stat_aphid_ff/stat_aphid_ff_sp_11.RDS")
stat_plankton_north_sea<-readRDS("./Results/stat_results/stat_plankton_north_sea/stat_plankton_north_sea_sp_16.RDS")
stat_methane<-readRDS("./Results/stat_results/stat_methane/stat_methane.RDS")
tab_multivar_summary<-readRDS("./Results/tab_multivar_summary.RDS")
```

\begin{figure}[!h]
\begin{center}
\includegraphics[width=6.5in]{./Results/BivarDataPlot.pdf}
\caption{Upper panels show raw data plots for (A) $\log_{10}$(soil N) vs. $\log_{10}$(soil C) data, (B) $\log_{10}$(basal metabolic rate (BMR)) vs. $\log_{10}$(body mass) 
for birds, (C) the same for mammals, and (D) above-ground plant biomass vs. Shannon's index from Cedar Creek. Bottom panels show the corresponding normalized rank plots. See section \ref{Data} for units used in upper panels. \label{fig_biv_multi_raw_cop}}
\end{center}
\end{figure}

Copula structures were non-normal and 
showed asymmetric tail 
associations for most, but not all datasets we examined, 
answering Q1 in the affirmative. 
To make results easier to absorb, we
present results first for an example bivariate dataset, then in summary
for all bivariate datasets, then for an example multivariate dataset, then 
in summary for all multivariate datasets.

<!--bivariate example, model selection-->
For the soil C and N data (section \ref{Data}, Table \ref{tab_data_info}, 
Fig. \ref{fig_biv_multi_raw_cop}A), variables were 
non-independent ($p=$ `r BivMS_res_Loecke$IndepTestRes`, to within the precision available 
from `BiCopIndTest` in the `VineCopula` package). Non-independence is also visually apparent.
The Kendall correlation was `r round(BivMS_res_Loecke$TauVal,2)`. We fitted our 16 copulas
to the normalized ranks (Fig. \ref{fig_biv_multi_raw_cop}E) 
and computed AIC and BIC weights and corresponding
lower- and upper-tail statistics for the fitted distributions (Table \ref{tab_soilCNfit}).
`r BivMS_res_Loecke$InfCritRes$copname[which.min(BivMS_res_Loecke$InfCritRes$AIC)]` was the 
minimum-AIC copula, with AIC `r round(min(BivMS_res_Loecke$InfCritRes$AIC),2)`, whereas the 
normal copula had much higher AIC, `r round(BivMS_res_Loecke$InfCritRes$AIC[1],2)`, answering the 
first part of Q1 for these data: the data have a dependence structure markedly distinct 
from a normal copula. Recall that AIC differences of $2$
or $3$ are considered meaningful and differences of $8$ or greater
are definitive, so that data provide essentially no support
for the higher-AIC model in that case [@burnham2003_modelselection].
Model-averaged lower- and upper-tail dependence statistics were `r round(BivMS_res_Loecke$relLTdep_AICw,3)` 
and `r round(BivMS_res_Loecke$relUTdep_AICw,3)`, respectively (AIC weights were used for averaging). 
These values are distinct from what a normal copula would give, namely, 0, helping to answer 
the second part of
Q1 for these data. The values also differed substantially 
from each other, helping answer the third part of Q1. 
The numbers reflect stronger upper- than lower-tail
dependence, and this is also visible in the extreme upper-right
and lower-left corners of the copula plot (Fig. \ref{fig_biv_multi_raw_cop}E).
Apparently C and N are more strongly associated with each other in soils
for which they are both plentiful than in soils for which they are 
both scarce, and C and N are more strongly associated in both
tails than would be modelled using a normal copula model or methods that
assume such a model.

However, our model selection procedures do not reveal whether tail dependence parameters are 
*significantly* different from $0$ and from each other. 
Furthermore, we caution that, while our model selection results do convincingly show 
non-normal copula structure, model-averaged tail-dependence statistics may have been biased 
because even the lowest-AIC copula 
(`r BivMS_res_Loecke$InfCritRes$copname[which.max(BivMS_res_Loecke$InfCritRes$AICw)]`) 
was a poor fit ($p=$ `r BivMS_res_Loecke$GofRes_CvM` and `r BivMS_res_Loecke$GofRes_KS` for the 
Cramer-von Mises and Kolmogorov-Smirnov goodness-of-fit tests, respectively, to within the precision 
available from `BiCopGofTest`).
Our nonparametric results, detailed next, provide information about significance of tail
associations and asymmetries of tail associations.

The values of $\cor_{l_b,u_b}$, $\Ps_{l_b,u_b}$ and $\Dsq_{l_b,u_b}$ for the soil
C and N data
were compared to distributions of their values on 1000 Kendall- or Spearman-preserving normal surrogates,
separately in two comparisons for each of the ranges $(l_b,u_b)=(0,0.1),\ldots,(0.9,1)$ 
(Fig. \ref{fig_soilCN_nonparam}). Results confirmed that tail associations of data are stronger in 
both the lower and upper tails than would be expected under a normal-copula null hypothesis (Fig. \ref{fig_soilCN_nonparam}B-F).
The values of the statistics $\cor_{0,0.1}-\cor_{0.9,1}$, $\Ps_{0,0.1}-\Ps_{0.9,1}$ and
$\Dsq_{0.9,1}-\Dsq_{0,0.1}$ for the soil C and N data were also compared to distributions of their values
on surrogates (Table \ref{tab_soilCN_asym}). Results confirmed that
upper-tail association was significantly stronger than lower-tail association; i.e., C and N values are
more related when high than when low. These results
are consistent with the model selection results, but go beyond them by providing information
about significance. Thus our results provide an
affirmative answer to Q1 for the soil C and N data; this answer is 
represented in Fig. \ref{MasterFigure} as a solid outline around "soil C and N"
in the right-most box in the middle row.

<!--bivariate summary-->
Table \ref{tab_bivar_summary} summarizes Q1 results 
for the bivariate datasets, with some details in the Appendices. The two 
variables were significantly related for all datasets 
(Table \ref{tab_bivar_summary}, row 1). A non-normal copula 
always emerged as the lowest-AIC copula 
(Table \ref{tab_bivar_summary}, row 2). 
The normal copula was a poor fit compared to the lowest-AIC copula 
(Table \ref{tab_bivar_summary}, rows 3-4), except for Cedar Creek, for which the AIC
difference between these fits was marginal.
Thus we answer the first part of Q1 in the affirmative for 3 of 4 of our bivariate datasets. 
Often, either lower- or upper-tail dependence statistics differed substantially from 0 
(Table \ref{tab_bivar_summary}, rows 7-8), and/or these statistics differed from each other (Table \ref{tab_bivar_summary}, row 9), 
helping to answer the second and third parts of Q1. Though again the model
selection results do not provide 
information on significance of these differences, and are subject to the
caveat that, for some datasets, even the lowest-AIC copula was not an objectively good fit
(Table \ref{tab_bivar_summary}, rows 5 and 6), nonparametric results 
(Table \ref{tab_bivar_summary}, row 10) showed that tail associations
deviated significantly from what would be expected from a normal-copula
null model, and were also significantly asymmetric, except for the Cedar Creek data.
See Figs \ref{SM-fig_bmr_birds_nonparam}-\ref{SM-fig_cc2000_nonparam} for analogues to 
Fig. \ref{fig_soilCN_nonparam} for the bivariate datasets other than the soil C and N
data, see Tables \ref{SM-tab_bmr_birds_fit}-\ref{SM-tab_cc2000fit} for analogues to Table \ref{tab_soilCNfit},
and see Table \ref{SM-tab_3biv_npa_summary} for analogues to Table \ref{tab_soilCN_asym}.
The first three datasets had stronger upper- 
than lower-tail dependence, and the last dataset had the reverse, though non-significantly. 
Apparently species body masses and metabolic rates, for both birds and mammals, are more strongly 
associated with each other when both these quantities are large than when they are 
small. Our generally affirmative (except for Cedar Creek), 
empirically based answer to Q1 for the bivariate datasets is 
represented in Fig. \ref{MasterFigure} as solid outlines around the names of 
the soil C and N and bird and mammal body mass versus metabolic rate
datasets, in the middle row of boxes.
We did the same analyses 
with Cedar Creek data from other available years, 1996-1999 and 2007. For 1996 and 
2007, independence of biomass and Shannon's index could not be rejected; 
for 1997-1999, results were similar to 2000. 

We present green spruce aphid abundance (section \ref{Data}; Table \ref{tab_data_info}) results
as an example multivariate analysis.
Independence was rejected for each pair of the 10 sampling locations. 
Best-fitting (lowest-AIC) copulas were non-normal
for the large majority of location pairs (`r tab_multivar_summary$Green_spruce_aphid_count$num.non_N` of `r tab_multivar_summary$Green_spruce_aphid_count$num.pairs` pairs), 
and the AIC
of the normal copula minus the minimum AIC, averaged across 
location pairs, was `r tab_multivar_summary$Green_spruce_aphid_count$Avg.N_minus_bestfit.AIC`.
Whereas `r tab_multivar_summary$Green_spruce_aphid_count$Avg.N_minus_bestfit.AIC`
would be only a marginally meaningful AIC difference for a single
location pair, it is more meaningful as an average across many
pairs. Relatedly, and illustrating the concept here, 
the chances of getting `r tab_multivar_summary$Green_spruce_aphid_count$num.non_N` or more 
non-normal location pairs if the chances were equal of getting
a normal versus a non-normal result (and taking into account that
the location pair (A,B) will necessarily 
produce the same result as the pair
(B,A)) is low, 
$`r 1-pbinom(tab_multivar_summary$Green_spruce_aphid_count$num.non_N/2-1,45,0.5)`$. 
The average AIC difference `r tab_multivar_summary$Green_spruce_aphid_count$Avg.N_minus_bestfit.AIC` was
much less than the AIC difference between `r round(min(BivMS_res_Loecke$InfCritRes$AIC),2)` and 
`r round(BivMS_res_Loecke$InfCritRes$AIC[1],2)` that was obtained for the soil C and N data
at least in part because the soil C and N data were much more numerous (Table \ref{tab_data_info}).
These results help answer the first part of Q1: non-normal
copula structure appears to be meaningfully common in
these data, even if not universal.
Goodness of fit tests in every case 
failed to reject the hypothesis that the AIC-best copula 
family was also 
an objectively adequate description of 
the data; i.e., the collection of 
copula families we used was sufficiently broad to characterize 
these data. Model-averaged lower- and upper-tail
dependence statistics had $2.5^{th}$ and $97.5^{th}$ quantiles (across
location pairs)
(`r tab_multivar_summary$Green_spruce_aphid_count$LT_0.025_CI`, `r tab_multivar_summary$Green_spruce_aphid_count$LT_0.975_CI`) and 
(`r tab_multivar_summary$Green_spruce_aphid_count$UT_0.025_CI`, `r tab_multivar_summary$Green_spruce_aphid_count$UT_0.975_CI`), respectively,
thereby commonly differing 
from what a normal copula would give (i.e., 0), and helping to
answer the second
part of Q1: these data have greater tail dependence 
(lower and upper) than would be expected from a normal-copula null hypothesis. We note however that our model selection 
procedures again do not reveal whether tail-dependence 
parameters are *significantly* different from $0$, and we refer 
the reader to nonparametric results below for that information. 
Model-averaged lower- minus upper-tail dependence statistics 
were positive for all but a few location pairs (`r sum((RES_aphid_count$LTdep_AICw-RES_aphid_count$UTdep_AICw)>0,na.rm=T)` 
out of `r (dim(RES_aphid_count$gfc_p_CvM)[1]^2-dim(RES_aphid_count$gfc_p_CvM)[1])-2*RES_aphid_count$num_indep_loc_pair`). The chances
of getting `r sum((RES_aphid_count$LTdep_AICw-RES_aphid_count$UTdep_AICw)>0,na.rm=T)` 
or more positive values, here, under a null hypothesis of equal chances 
for positive and negative values (and again accounting for the fact that location pairs (A,B) and (B,A) will show the same result) was
again low, 
$`r 1-pbinom(sum((RES_aphid_count$LTdep_AICw-RES_aphid_count$UTdep_AICw)>0,na.rm=T)/2-1,45,0.5)`$.
Thus the spatial synchrony of rarity in the green spruce aphid is stronger 
than the spatial synchrony of outbreaks. This answers the third part of Q1.

Nonparametric statistics verified that tail associations were asymmetric for the green 
spruce aphid abundance data. Mean values across all pairs of 
locations of the Spearman and Kendall correlations and the statistics 
$\cor_l$, $\cor_u$, $\Ps_l$, $\Ps_u$, $\Dsq_l$, and $\Dsq_u$ were positive and
confidence intervals excluded $0$ (Table \ref{tab_count_resamp}). 
Mean values of the statistics $\cor_l-\cor_u$, $\Ps_l-\Ps_u$, and $\Dsq_u-\Dsq_l$
also were always positive and had confidence intervals that 
excluded $0$ (Table \ref{tab_count_resamp}).

<!--multivariate summary-->
Table \ref{tab_multivar_summary} summarizes Q1 results for
multivariate data. Results supported the 
conclusions that non-normal copula structure, non-normal tail dependence, and asymmetric tail 
associations were common, though not universal, answering 
Q1 in the affirmative for these data. 
Most location pairs were non-independent 
(Table \ref{tab_multivar_summary}, row 2). The large majority of non-independent location 
pairs had best-fitting copulas that were not the normal copula (Table \ref{tab_multivar_summary}, 
row 3), and AIC values of best-fitting copulas were, on average
across location pairs, between $2.714$ (*Ceratium furca* abundance data)
and $5.764$ (leaf-curling plum aphid first-flight data) lower than
AIC values for the normal copula (Table \ref{tab_multivar_summary}, 
row 4). Best-fitting copulas were nearly always considered an adequate fit 
(Table \ref{tab_multivar_summary}, row 5). Some of the datasets (green spruce aphid
abundance, \textit{Ceratium furca} abundance, methane-flux) showed stronger 
lower- than upper-tail dependence 
(Table \ref{tab_multivar_summary}, rows 6-8), whereas leaf-curling plum aphid first flight
data showed the reverse. 
Thus the spatial synchrony of rarity was stronger than that of outbreaks for 
green spruce aphid and \textit{Ceratium furca} abundance, for instance, whereas leaf-curling aphid
first flights were more spatially related to each other when these flights happened
late than when they happened early.
All of the values in Table \ref{tab_multivar_summary}, rows 8 deviated highly significantly
from what would have been expected under a null hypothesis 
of equal chances for positive and negative values. 

Asymmetry results were generally verified by nonparametric approaches,
with the exception of the methane-flux data. 
For instance, the $95\%$ confidence intervals of the
mean over pairs of sampling locations of the statistic $\cor_l-\cor_u$ were 
(`r unname(round(stat_aphid_ff$numericdf[9,c(3,4)],3))`)
for the leaf-curling plum aphid first flight data
(Table \ref{tab_multivar_summary}, row 9), 
indicating greater upper-tail dependence, and consistent with
the results of Table \ref{tab_multivar_summary}, row 8. 
For the *Ceratium furca* abundance data, confidence intervals were
(`r unname(round(stat_plankton_north_sea$numericdf[9,c(3,4)],3))`),
indicating greater lower-tail dependence, and again 
consistent with the results of Table \ref{tab_multivar_summary}, row 8.
For the methane-flux data, confidence intervals were 
(`r unname(round(stat_methane$numericdf[9,c(3,4)],3))`);
the asymmetry in tail dependence revealed by the model selection
results for the methane-flux data was apparently not strong enough to
also be revealed by the nonparametric analyses. 
Analogues to Table \ref{tab_count_resamp} are in Tables
\ref{SM-tab_ff_resamp}-\ref{SM-tab_methane_resamp}.
Our generally
affirmative, empirically based answer to Q1 for the multivariate datasets is 
represented in Fig. \ref{MasterFigure} as solid outlines around the names
of those datasets, in the middle row of boxes.
We also carried out the same analyses for abundance and first-flight data for the 18 
aphid species for which we had data other than the green spruce and leaf-curling plum aphids,
as well as for the 21 plankton taxa for which we had abundance data other than
\textit{Ceratium furca} (results not shown). Results supported the conclusion that 
non-normal copula structure and tail dependence, and asymmetry of tail associations, are 
common.

<!--bivariate results for soil CN data-->

<!--Display model selection results for soilCN data in a table-->
```{r tab_soilCNfit, echo=F, results="asis"}
#knitr::opts_knit$set(kable.force.latex = TRUE)
BivMS_res_Loecke<-readRDS("./Results/fitting_results/BivMS_soilCN.RDS")
h<-BivMS_res_Loecke$InfCritRes[,c("copname","AIC","AICw","BIC","BICw","LTdep","UTdep")]
colnames(h)<-c("Copula","AIC","AICw","BIC","BICw","LT","UT")
h_ord <- h[order(h$AIC),] #sorting data: min AIC to max AIC

#---- This line forces the formating which is not always true by using digits arg in kable function ----
h_ord$AIC<-formatC(h_ord$AIC,2,format="f")
h_ord$AICw<-formatC(h_ord$AICw,2,format="f")
h_ord$BIC<-formatC(h_ord$BIC,2,format="f")
h_ord$BICw<-formatC(h_ord$BICw,2,format="f")
h_ord$LT<-formatC(h_ord$LT,4,format="f")
h_ord$UT<-formatC(h_ord$UT,4,format="f")
#-------------------------------------------------

knitr::kable(h_ord, #longtable = T, 
             format = "latex",linesep = "",row.names = F, #digits = c(4,2,2,2,2,4,4), 
      caption.short = "Model fitting results for soil C and N dataset",
      caption = "Model fitting results for soil C and N dataset using 16 copula families: normal (N), Frank (F), Clayton (C), survival Clayton (SC), Gumbel (G), survival Gumbel (SG), Joe (J), survival Joe (SJ), BB1, survival BB1 (SBB1), BB6, survival BB6 (SBB6), BB7, survival BB7 (SBB7), BB8, and survival BB8 (SBB8). Table rows are sorted by AIC. AICw = AIC weight; BICw = BIC weight; LT = the lower-tail dependence statistic for the indicated copula family with fitted parameters; UT = the same for upper-tail dependence. \\label{tab_soilCNfit}",
      booktabs = TRUE)
```

<!--Displaying a table with npa summary results for soilCN data-->
```{r tab_soilCN_asym,echo=F,results="asis"}
tab_biv_npa_summary<-readRDS("./Results/stat_results/tab_biv_npa_summary.RDS")
tab_soilCN_asym<-tab_biv_npa_summary[1:3,2:4]
knitr::kable(tab_soilCN_asym, booktabs = T, escape=F,
             format = "latex",linesep = "",
      caption.short = "Nonparametric asymmetry tests against a normal-copula null, soil C and N",
      caption = "Soil C and N results for nonparametric tests of whether asymmetry of tail association was significant, compared to a normal-copula null hypothesis. The values of the listed statistics (cor$_{FS}$-cor$_{LS}$ = $\\cor_{0,0.1}-\\cor_{0.9,1}$, P$_{FS}$-P$_{LS}$ = $\\Ps_{0,0.1}-\\Ps_{0.9,1}$, D$_{LS}^2$-D$_{FS}^2$ = $\\Dsq_{0.9,1}-\\Dsq_{0,0.1}$; here FS stands for `first slice' and LS for `last slice') for real data were compared to their values for 1000 Kendall- or Spearman-preserving normal surrogates, in separate comparisons. A table entry $<X$ indicates the value of the given statistic on the data was less than its value on $X$ of the surrogates, so entries of the form $<X$ for $X$ equal to $975$ or above indicate that upper-tail dependence was significantly stronger than lower-tail dependence. Results were significant only for the $P$ and $D^2$ statistics.\\label{tab_soilCN_asym}")
```

\begin{figure}[!h]
\begin{center}
\includegraphics[width=13 cm]{./Results/stat_results/stat_soilCN/soilCN_bivfunctionplot.pdf}
\caption[Nonparametric tests against a normal-copula null, soil C and N data]{Nonparametric tests for tail association and other deviations from normal copula structure, for soil C and N data. As described in the main text, the values of the statistics $\cor_{l_b,u_b}$, $\Ps_{l_b,u_b}$ and $\Dsq_{l_b,u_b}$ for real data (crosses) were compared to distributions of their values on 1000 Kendall- or Spearman-preserving normal surrogates of the data (dark and light grey x's show 0.025 and 0.975 quantiles), separately in two comparisons for each of the ranges $(l_b,u_b)=(0,0.1),\ldots,(0.9,1)$. Whenever the cross was outside the range given by the x's, text at the top of panels indicates the number of surrogate values the real-data value was greater than or less than. For instance, a value $>N$ (respectively, $<N$) means the value of the statistic on real data was greater than (respectively, less than) its value on $N$ surrogates. When the statistic was greater than 975 or less than 975 surrogate values, it indicates significance ($95\%$ confidence level). When $\cor$ or $\Ps$ values (respectively, $\Dsq$ values) were greater than surrogates, it means association in that part of the distributions was stronger than (respectively, weaker than) expected from a normal-copula null hypothesis. \label{fig_soilCN_nonparam}}
\end{center}
\end{figure}

<!--summary table : results with bivariate data-->
```{r tab_bivar_summaryl, echo=F, results='asis',message=F}
tab_bivar_summary<-readRDS("./Results/tab_bivar_summary.RDS")
tab_biv_npa_summary<-readRDS("./Results/stat_results/tab_biv_npa_summary.RDS")

# adding npa stat summary res: spearman cor. preserving corl-coru for 4 bivariate data
my_mat<- matrix(tab_biv_npa_summary[c(1,4,7,10),4],nrow=1,ncol=4)
colnames(my_mat)<-colnames(tab_bivar_summary) # this line is needed to apply rbind later: colnames should be matched

tab_bivar_summary<-rbind(tab_bivar_summary,my_mat)
row.names(tab_bivar_summary)<-NULL

library(kableExtra)
library(dplyr)
library(tibble)

dt <- tibble(Description = c("1. p-value, independence test",
           "2. Minimum-AIC copula",
           "3. AIC for best copula",
           "4. AIC for normal copula",
           "5. p-value, Cramer-von Mises goodness of fit test",
           "6. p-value, Kolmogorov-Smirnov goodness of fit test",
           "7. Model averaged lower-tail dependence",
           "8. Model averaged upper-tail dependence",
           "9. Model averaged lower- minus upper-tail dependence",
           "10. Asymmetry rank in surrogates", # this row will be commented out, that's why the numbering starts from 10 in the next line
           "10. Placement of cor$_{FS}$-cor$_{LS}$ in surrogate distribution (see caption)")
)
tab_bivar_summary<-cbind(dt,tab_bivar_summary)
tab_bivar_summary<-tab_bivar_summary[c(1:9,11),]

knitr::kable(tab_bivar_summary, 
             format="latex", escape=F,row.names=F,
             caption = "Summary of Q1 results for bivariate datasets. The $p$-values (rows 5-6) are for the minimum-AIC copula. Model averaging used for rows 7-9 was based on AIC weights. Row 10 shows values as in Table \\ref{tab_soilCN_asym}, upper right table entry. Although the result shown in row 10 for the soil C and N data is non-significant, see Table \\ref{tab_soilCN_asym} for significant results using the $P$ and $D^2$ statistics. The first 3 datasets use $(l_b,u_b)=(0,0.1)$ for the first slice (FS) and $(0.9,1)$ for the last slice (LS), whereas $(0, 0.25)$ and $(0.75,1)$ were used for Cedar Creek due to fewer data for that system. Results for Cedar Creek for the year 2000 are shown. See Table \\ref{tab_soilCNfit} for copula family abbreviations.\\label{tab_bivar_summary}",
             booktabs=T, linesep = "\\addlinespace",
             col.names = c(" ", "Soil C and N", "Bird masses and BMR", "Mammal masses and BMR", "Cedar Creek data")
             )%>%column_spec(2:5,width="1.5cm")%>%
             column_spec(1, width = "8.4cm")%>%kable_styling(latex_options = "basic",position = "center")
```

```{r tab_count_resamp,echo=F,results="markup"}
s<-stat_aphid_count$numericdf[,c(3,2,4)]
sc<-as.data.frame(c("Spearman","Kendall","cor$_{l}$","cor$_{u}$","P$_{l}$","P$_{u}$","D$_{l}^2$","D$_{u}^2$",
                    "cor$_{l}$ - cor$_{u}$","P$_{l}$ - P$_{u}$","D$_{u}^2$ - D$_{l}^2$"))
s<-cbind(sc,s)
knitr::kable(s, digits = 3,
             col.names=c("","2.5$^{th}$ quantile", "Mean", "97.5$^{th}$ quantile"),
             format = "latex",linesep = "", escape = F,align="c",
             caption.short="Results for nonparametric statistics, green spruce aphid abundance data",
             caption = "Average values of statistics across all locations pairs and confidence intervals based on spatial resampling for green spruce aphid abundance data.\\label{tab_count_resamp}",
             booktabs=TRUE)
```

<!--summary table : results with multi-variate data-->
```{r tab_multivar_summaryl, echo=F, results='asis',message=F}
tab_multivar_summaryl<-tab_multivar_summary[-c(4,12:15),] #exclude extra rows for AICw, Corl, Coru

tab1<-tab_multivar_summaryl[c(6,8,12),] # see carefully if you rearrange tab_multivar_summary's row
                                                   #in suppmat this will change
tab2<-tab_multivar_summaryl[c(7,9,13),]  # see carefully if you rearrange tab_multivar_summary's row
                                                   #in suppmat this will change

two_tab_in1 <- as.data.frame(do.call(cbind, lapply(1:ncol(tab1), function(i) paste0(" (",tab1[,i], ", ", tab2[ , i], ")") )))   # see carefully if you rearrange tab_multivar_summary's row in suppmat this will change

two_tab_in1_again <- as.data.frame(do.call(cbind, lapply(1:ncol(two_tab_in1), 
                                                         function(i) paste0(tab_multivar_summaryl[11,i], ", ", two_tab_in1[3,i]) ))) 
  
colnames(two_tab_in1)<-colnames(two_tab_in1_again)<-colnames(tab_multivar_summaryl)


compiled_tab<-rbind(two_tab_in1[1:2,],two_tab_in1_again)
rownames(compiled_tab)<-c("CI_LT","CI_UT","Avg_CI_Corl-Coru")

tab_multivar_summary_ed<-rbind(tab_multivar_summaryl[c(1:5),],
                            compiled_tab[1:2,], #CI of LT, CI of UT
                            tab_multivar_summaryl[10,],
                            compiled_tab[3,] #avg and CI of Corl-Coru
                           )

row.names(tab_multivar_summary_ed)<-NULL
library(kableExtra)
library(dplyr)
library(tibble)
dt <- tibble(
            Description = c("1. Location pairs (excluding self comparisons)",
           "2. Number of non-independent pairs",
           "3. Number of pairs with non-normal copula as best fit",
           "4. Average of AIC differences (normal AIC minus min AIC) across location pairs",
           "5. Number of well-fitting location pairs",
           "6. 2.5$^{th}$ and 97.5$^{th}$ quantiles of model-avg. lower-tail dependence",
           "7. 2.5$^{th}$ and 97.5$^{th}$ quantiles of model-avg. upper-tail dependence",
           "8. Percent pairs showing stronger lower tail dependence",
           "9. Mean with 2.5$^{th}$ and 97.5$^{th}$ quantiles of cor$_{l}$ - cor${_u}$ across all location pairs")
           )

tab_multivar_summary_ed<-cbind(dt,tab_multivar_summary_ed)

knitr::kable(tab_multivar_summary_ed, 
             format="latex", escape = F,
             caption.short = "Summary of results for multivariate datasets",
             caption = "Summary of Q1 results for multivariate datasets. Rows 3-8 of the table were computed for the non-independent pairs (row 2) only. A well-fitting location pair (row 5) was one for which the best-fitting copula had $p$-values $>0.01$ for both the Cramer-von Mises and Kolmogorov-Smirnov goodness of fit tests. \\label{tab_multivar_summary}",
             booktabs=T, linesep = "\\addlinespace",
             col.names = c(" ","Green spruce aphid abundance",
                        "Leaf-curling plum aphid first flight",
                        "\\textit{Ceratium furca} abundance","Methane-flux")
             )%>%column_spec(2:5,width="2.5cm")%>%
             column_spec(1, width = "6cm")%>%kable_styling(latex_options = "basic",position = "center")
```

<!--Question 2 analyses-->
# Concepts and methods for Q2 \label{Causes}

Recall that Q2 is: what are some possible causes/mechanisms of non-normal
copula structure and asymmetric tail associations in ecology?
Having demonstrated that non-normal copula structure and asymmetric tail associations are 
common, we addressed Q2 by exploring, using models, three possible mechanisms
for these phenomena. 
Our models presented are initial explorations, only, of whether the proposed 
mechanisms may explain observed patterns. As such, simple models were used. 
Comprehensive explorations of
model parameter space and alternative model formulations were left for future work.

The first mechanism relates to the ideas in the Introduction about 
Liebig's law of the minimum, and to nonlinear influences of environmental
variables on ecological variables.
If an environmental variable influences an ecological variable 
disproportionately in one of its tails, we explored whether the ecological
variable could then exhibit asymmetric tail associations across space.
Let $E_i(t)$ be an ecological variable 
measured at location $i$ ($i=1,2$) at time $t$.
Assume the dynamics 
$E_i(t+1)=bE_i(t)+g(\varepsilon_i(t))+a \delta_i(t)$, where 
the $\delta_i(t)$ are standard-normally distributed 
and independent across time and locations, $a=0.2$,
$b=0.1, -0.1, 0.5, -0.5$ in different simulations, and 
the $(\varepsilon_1(t),\varepsilon_2(t))$ were drawn, independently through
time, from a bivariate normal distribution with $\var(\varepsilon_i)=1$ and 
$\cov(\varepsilon_1,\varepsilon_2)=0.8$. Thus basic ecological dynamics follow the
very simple AR(1) formulation, influenced by a "regional" environmental factor which
is correlated across locations ($\varepsilon$) and by a local factor ($\delta$). 
The function $g$ describes how $\varepsilon_i(t)$ influences
$E_i(t)$. We used $g$ equal to $g_1$ or $g_2$ in different simulations:
$g_1(\varepsilon)$ equals $\varepsilon$ if $\varepsilon<0$ and $0$ otherwise;
$g_2(\varepsilon)$ equals $\varepsilon$ if $\varepsilon>0$ and $0$ otherwise. Thus
$g_1$ represents environmental effects that negatively impact populations,
but only below the threshold $\varepsilon=0$; 
and $g_2$ represents effects that positively impact populations, but 
only above $\varepsilon=0$.
The values of $b$ provide a modicum of exploration of whether ecological 
dynamics may also influence how complex dependencies between
the $E_i(t)$ arise: 
negative $b$ corresponds to overcompensating dynamics 
and positive $b$ to undercompensating dynamics; larger $|b|$
means slower return to equilibrium after a disturbance. 

For each $b$ and $g_i$, we simulated the model for $25000$ time
steps and retained the $E_i(t)$ for the final $2500$ time 
steps. We applied our 
nonparametric statistic $\cor_{0,0.1}-\cor_{0.9,1}$ 
and our Spearman- or Kendall-preserving normal-copula 
surrogate comparison methods (see section \ref{Methods}) to these 
outputs to discover if the model could
produce asymmetric tail associations (and therefore non-normal copula structure)
between the $E_i(t)$. Because $(\varepsilon_1,\varepsilon_2)$
and $(\delta_1,\delta_2)$ have normal copula structure (they were drawn
from bivariate normal distributions), the Moran 
mechanism analyzed below does not operate here.

Our second mechanism is an extension of the well-known Moran
effect, and was summarized conceptually in the Introduction. 
We consider a linear model, as well as 
two parameterizations of a nonlinear population model which includes density dependence
of population growth. For the linear model, let $E_i(t)$ 
again be an ecological variable, $i=1,2$.
We use AR(1) dynamics, 
$E_i(t+1)=\beta E_{i}(t)+ \sqrt{1-\beta^2} \varepsilon_i(t)$, with $\beta=0.5$. 
The environmental noises $\varepsilon_i(t)$ were 
standard-normal random variables that were independent for distinct times, $t$, but 
exhibited different kinds of 
dependence across locations in different simulations (see below).
The variable $E$ is general. It could represent deviations
of a population density from a carrying capacity, deviations of
total plant community biomass from an average value,
flux of a biogeochemical variable such as methane, or other quantities.
The nonlinear population model was a stochastic, multi-habitat-patch version 
of the Ricker model, 
$P_i(t+1)=P_i(t)\exp \left[r\left(1-P_i(t)/K  \right)+\sigma\varepsilon_i(t)\right]$ for $i=1,2$, using $r=0.5$, 
$K=100$, $\varepsilon_i(t)$
as above, and $\sigma=0.1$ or 
$\sigma=1$ in different simulations. 
When $\sigma=0.1$, population dynamics stay close to the carrying capacity, 
$K$, where the model equation can be well approximated by a linear equation,
and the nonlinearities of the model therefore have limited influence on dynamics. When
$\sigma=1$, model dynamics are strongly nonlinear because the stochastic 
component of the model causes populations to stray far from the carrying capacity. 
We refer to these as
the weak-noise and strong-noise cases, though the importance of the noise here
is that, when strong, it brings the nonlinearities of the model into play. 

For each of the three model setups above (the linear model and the nonlinear model
with weak and strong stochasticity), for the Clayton and survival Clayton copula families, 
for each $\tau=0.1,0.2,\ldots,0.9$, and for each of
$50$ replicate simulations, we generated 5000 noise pairs 
$(\varepsilon_1(t),\varepsilon_2(t))$ from the
bivariate random variable with standard-normal marginals and with 
the given copula family and the given Kendall
correlation $\tau$.
We then used this noise to drive the model, and retained both the
noise and population values for the final 500 
time steps. For each simulation, the following
statistics were then computed for noises and populations: 
Pearson, Spearman and Kendall correlations, $\cor_l$, $\cor_u$, $\Ps_l$, $\Ps_u$, 
$\Dsq_l$, $\Dsq_u$, $\cor_l - \cor_u$, $\Ps_l - \Ps_u$, and $\Dsq_u - \Dsq_l$. 
Values were plotted against $\tau$ for noises and populations. If the 
hypothesis from the Introduction was reasonable that characteristics of 
the copula structure of spatial dependence in an
ecological variable may be inherited from characteristics of spatial 
dependence in an environmental variable through a Moran-like effect, then 
plots should be similar for populations and noises. 

The next mechanism we investigated is evolutionary, and pertains to bivariate 
trait data across species, e.g., our bird and mammal 
data. This mechanism is a hypothetical explanation for the bias toward right-tail association
observed in those data (Fig. \ref{fig_biv_multi_raw_cop}F,G). 
The hypothesis is that asymmetric tail association occurs in evolutionary changes
in bivariate characters, and gives rise to asymmetric tail association between the 
two character values across extant species.
We simulated bivariate character evolution on an estimate of 
the phylogeny, taken from @GenoudIM2018, of $817$ mammal species.
The root character state and change 
across each branch were randomly chosen from matrices of one million 
independent draws from bivariate distributions showing one of five distinct types of copula 
structure: 1) extreme or 2) moderate left-tail dependence,
3) symmetric tail dependence, or 4) moderate or 5) extreme
right-tail dependence (Appendix \ref{SM-evol_mech}). 
All distributions had standard-normal marginals and Spearman 
correlation $0.875$ between components, so our simulations assess 
the impact of copula structure only. For 
each of the five copulas, mammalian character
evolution was simulated 100 times. For each simulation, symmetry of tail associations 
of the two characters across phylogeny tips was assessed using our nonparametric
statistics.
We hypothesized that cases 1 and 2 above would yield stronger left- than right-tail 
associations in tip characters, and cases 4 and 5 would yield the reverse. 
The simulator was written in Python and used version 4.4 of the DendroPy 
package [@DendroPy]. 

# Results for Q2: Moran effects and asymmetric dependencies produce non-normal copula structure\label{Results_Q2}

Our model with asymmetric environmental effects 
produced outputs with visually apparent asymmetry 
of tail associations between the ecological variables $E_i$ in the two locations, 
to an extent that depended on $b$.
For $b=0.1$ (Fig. \ref{fig_ms_asym_b_0.1}) and $b=-0.1$ 
(Fig. \ref{SM-fig_ms_asym_b_-0.1_0.5_-0.5}A,B), 
for both $g_1$ and $g_2$,
asymmetry of tail association was strong; for 
$b=\pm 0.5$, asymmetry was weaker but still apparent
(Fig. \ref{SM-fig_ms_asym_b_-0.1_0.5_-0.5}C-F).
Lower-tail (respectively, upper-tail) spatial associations in the 
effects of noise ($g_1$, respectively, $g_2$) 
produced lower-tail (respectively, upper-tail) spatial 
associations in the ecological variable, $E_i$.
Results using our statistic $\cor_{0,0.1}-\cor_{0.9,1}$ 
strongly reflected the visually apparent asymmetry (Table \ref{tab_stat_npa_asym}). 
Thus asymmetry of environmental effects is a mechanism that may be partly 
responsible for non-normal copula structure and
asymmetric tail associations across space in ecological variables:
when populations or other ecological quantities are influenced principally by
low (respectively, high) values of an environmental variable, one may expect 
left-tail (respectively, right-tail) associations across space for the ecological variable. 
This result is represented in Fig. \ref{MasterFigure} as the solid box around 
"Nonlinear environmental effects, Liebig's law" and the arrows labelled "A".
It is explained in the Discussion why the box and some of the arrows 
are solid, instead of dashed, although the results here are theoretical.

\begin{figure}[!h]
\begin{center}
\includegraphics[width=12 cm]{./Results/asym_sens_results/a_0.2_b_0.1_r_0.8_asym_sens.pdf}
\caption{If environmental effects operate asymmetrically in their tails on
ecological variables, it can result in non-normal copula structure and asymmetric 
tail dependence across space in the ecological variables. Shown are the last
$500$ points for (A) $g=g_1$ and (B) $g=g_2$ from simulations described in the text, $b=0.1$. 
Asymmetric tail associations are visually apparent, but were also elaborated statistically in 
Table \ref{tab_stat_npa_asym}. \label{fig_ms_asym_b_0.1}}
\end{center}
\end{figure}

```{r tab_stat_npa_asym_simple,echo=F,results="markup"}
library(kableExtra)
source("./make_table_stat_fn.R")
tab_asym<-readRDS("./Results/asym_sens_results/stat_npa_asym_b_0.1_b_0.5.RDS")
tab_stat_npa_asym<-c()
for(i in c(1:8)){
  temp<-make_table_stat_fn(data=tab_asym[[i]])
  tab_stat_npa_asym<-rbind(tab_stat_npa_asym,temp)
}

ind_corstat<-which(tab_stat_npa_asym$Statistic %in% c("cor$_{FS}$-cor$_{LS}$"))
tab_corstat_asym<-tab_stat_npa_asym[ind_corstat,]
tab_corstat_asym<-tab_corstat_asym[,-1]

#gtext<-c("b=0.1, g$_{1}$","b=0.1, g$_{2}$","b=-0.1, g$_{1}$","b=-0.1, g$_{2}$","b=0.5, g$_{1}$","b=0.5, g$_{2}$","b=-0.5, #g$_{1}$","b=-0.5, g$_{2}$")
#tab_stat_npa_asym<-cbind(gtext,tab_corstat_asym)

tab_corstat_asym<-cbind(b=c(0.1,0.1,-0.1,-0.1,0.5,0.5,-0.5,-0.5),
  g=c("g$_{1}$","g$_{2}$","g$_{1}$","g$_{2}$","g$_{1}$","g$_{2}$","g$_{1}$","g$_{2}$"),tab_corstat_asym)

rownames(tab_corstat_asym)<-NULL

knitr::kable(tab_corstat_asym, booktabs = T, 
             format = "latex",linesep = "",escape=F, 
             caption.short = "Asymmetric sensitivity model, nonparametric statistics results",
             caption ="Asymmetric sensitivity model, nonparametric statistics results. 
Outputs of simulations from the first model detailed in 
section \\ref{Causes} (i.e., Figs \\ref{fig_ms_asym_b_0.1} and 
\\ref{SM-fig_ms_asym_b_-0.1_0.5_-0.5}) were subjected to nonparametric statistical 
analyses described in section \\ref{Methods}. The statistic cor$_{FS}$ - cor$_{LS}$ was computed for the
output of the model with the indicated b and g, and the value was compared to 1000 values
of the same statistic computed on surrogate time series randomized to have normal copula structure
(see section \\ref{Methods}). A table entry $<X$ (respectively, $>X$) indicates the value of the 
given statistic on the data was less (respectively, more) than its value on $X$ of the 
surrogates, so entries of the form $<X$ (respectively, $>X$) for $X$ equal to $975$ or above 
indicate that upper-tail (respectively, lower-tail) dependence was significantly stronger than 
lower-tail (respectively, upper-tail) associations compared to a 
normal-copula null hypothesis. `FS' stands for `first slice'
and refers to the bounds 0 to 0.1; `LS' stands for `last slice' and refers to the bounds
0.9 to 1.\\label{tab_stat_npa_asym}")
```

<!--Second mechanism, the Moran one-->
The hypothesis that characteristics of 
the copula structure of spatial dependence 
in an ecological variable can be inherited from an environmental 
variable was found to be reasonable, because it held for our models -
the Moran effect seems to extend to copula structure and tail associations - 
however, similarities between environmental- and 
ecological-variable copula structure were reduced when 
dynamics were strongly nonlinear.
For simulations using AR(1) models, our correlation and tail asymmetry statistics 
were always similar for noise and populations 
(Fig. \ref{fig_Cause4copula_C_kend_ar1} for the Clayton copula, 
Fig. \ref{SM-fig_Cause4copula_SC_kend_ar1} for the survival Clayton copula).
Though there were significant differences for many statistics
and simulations, these were small compared to the overall tendency for
larger (respectively, smaller) values of our statistics, as applied to
noise, to be paired with larger (respectively, smaller) values of the same
statistics for model outputs. 
For our nonlinear model with weak noise, values of the statistics
were again quite similar for noise and model outputs 
(Fig. \ref{SM-fig_Cause4copula_C_kend_sRicker_nsd_0.1} for Clayton, 
Fig. \ref{SM-fig_Cause4copula_SC_kend_sRicker_nsd_0.1} for survival Clayton).
Though there were again significant differences for many statistics
and simulations, these were again small relative to overall variation of 
values of the statistics.
Since many ecological models are nonlinear, this result provides the 
reasonable expectation for a Moran-effect-like correspondence 
between noise and model-output dependence structure across space 
for typical ecological dynamics,
as long as environmental noise is small enough that dynamics stay relatively close
to the model equilibrium. Theoretical results that hold for "weak noise" in this sense
are common in ecology.
For strong noise and using our nonlinear model, 
our correlation and tail-association-asymmetry statistics, generally 
speaking, were approximately similar between noise and model
outputs; however, similarity was reduced, 
and for a few simulations, asymmetry statistics
had opposite signs for noise and model outputs 
(Figs \ref{SM-fig_Cause4copula_C_kend_sRicker_nsd_1} and \ref{SM-fig_Cause4copula_SC_kend_sRicker_nsd_1}). For instance, 
using a Clayton copula with a large Kendall $\tau$,
$\cor_l - \cor_u$ was slightly positive for noise, but slightly negative for 
population outputs (Fig. \ref{SM-fig_Cause4copula_C_kend_sRicker_nsd_1}).

We repeated our analyses using the nonlinear model with $r=1.3$. The deterministic 
one-habitat-patch Ricker model exhibits a monotonic approach to a stable equilibrium
at $K$ when $r<1$ (undercompensating dynamics, e.g., the value $r=0.5$ used initially), 
but exhibits an oscillatory approach when $r>1$ (overcompensating
dynamics, e.g., $r=1.3$). For weak noise, similarities were again dominant between
values of our statistics on noise and population time series. 
For strong noise, however, discrepancies were often glaring. Apparently 
noise of standard deviation 1 interacted especially strongly with 
model nonlinearities when the model was in its overcompensatory 
regime. We repeated all analyses using the Gumbel, survival Gumbel, Joe, and 
survival Joe copulas, with substantially similar main conclusions (not shown).
 
Thus it is reasonable to hypothesize that a Moran-effect-like 
mechanism may produce non-normal copula structure and
asymmetric tail associations across space in ecological variables.
This result is represented in Fig. \ref{MasterFigure} as the dashed box around 
"Moran effects" and the arrows labelled "B".

<!--***DAN: Decided to remove this simplified version of the figure for AER
and replace by the bigger one which was in the supp mat, now below. Keeping the
code the simpler one here in case we ever need it.
\begin{figure}[!h]
\begin{center}
\includegraphics[width=5 cm]{./Results/Cause4copula_results/Cause4copula_stat_results_ar1_beta_0.5_Clayton/Clayton_scatter_CorlmCoru_Kendall's Tau.pdf}
\includegraphics[width=5 cm]{./Results/Cause4copula_results/Cause4copula_stat_results_ar1_beta_0.5_Clayton/Clayton_scatter_PlmPu_Kendall's Tau.pdf}
\includegraphics[width=5 cm]{./Results/Cause4copula_results/Cause4copula_stat_results_ar1_beta_0.5_Clayton/Clayton_scatter_D2umD2l_Kendall's Tau.pdf}
\caption{Example results showing similarity of copula statistics for 
environmental-noise inputs and ecological-variable outputs of a dynamical model.
Asymmetric tail dependence was similar in both noise inputs and model outputs
for this model. The 
AR(1) model and a Clayton copula were used (see section \\ref{Causes}).
Each point is the mean across 50 replicate simulations for which the same 
Kendall $\tau$ value was used. Error bars are standard errors 
and panel headers give Pearson correlation results for
the points. The regression line through the points (black line) was similar
to the 1-1 line (green line). 
See Fig. \ref{SM-fig_Cause4copula_C_kend_ar1} for additional results for the AR(1) model and Clayton copula, 
Figs \ref{SM-fig_Cause4copula_SC_kend_ar1}-\ref{SM-fig_Cause4copula_SJ_kend_ar1} for AR(1) results with other copula families, and 
Figs \ref{SM-fig_Cause4copula_C_kend_sRicker_nsd_0.1}-\ref{SM-fig_Cause4copula_SJ_kend_sRicker_nsd_1} for
the nonlinear models described in the methods section for Q2.\label{fig_Cause4copula_scatter_C_kend_ar1}}
\end{center}
\end{figure}
-->

<!--Plots for AR(1) with Clayton-->
\begin{figure}[!h]
\begin{center}
\includegraphics[width=13 cm]{./Results/Cause4copula_results/Cause4copula_stat_results_ar1_beta_0.5_Clayton/common_legend_cause4copula_stat.pdf}\\
\vspace{-0.4 cm}
\includegraphics[width=4.3 cm]{./Results/Cause4copula_results/Cause4copula_stat_results_ar1_beta_0.5_Clayton/Clayton_Spearman_vs_Kendall's Tau.pdf}
\includegraphics[width=4.3 cm]{./Results/Cause4copula_results/Cause4copula_stat_results_ar1_beta_0.5_Clayton/Clayton_Kendall_vs_Kendall's Tau.pdf}
\includegraphics[width=4.3 cm]{./Results/Cause4copula_results/Cause4copula_stat_results_ar1_beta_0.5_Clayton/Clayton_Pearson_vs_Kendall's Tau.pdf}\\
\vspace{-0.2 cm}
\includegraphics[width=4.3 cm]{./Results/Cause4copula_results/Cause4copula_stat_results_ar1_beta_0.5_Clayton/Clayton_Corl_vs_Kendall's Tau.pdf}
\includegraphics[width=4.3 cm]{./Results/Cause4copula_results/Cause4copula_stat_results_ar1_beta_0.5_Clayton/Clayton_Coru_vs_Kendall's Tau.pdf}
\includegraphics[width=4.3 cm]{./Results/Cause4copula_results/Cause4copula_stat_results_ar1_beta_0.5_Clayton/Clayton_Corl-Coru_vs_Kendall's Tau.pdf}\\
\vspace{-0.2 cm}
\includegraphics[width=4.3 cm]{./Results/Cause4copula_results/Cause4copula_stat_results_ar1_beta_0.5_Clayton/Clayton_Pl_vs_Kendall's Tau.pdf}
\includegraphics[width=4.3 cm]{./Results/Cause4copula_results/Cause4copula_stat_results_ar1_beta_0.5_Clayton/Clayton_Pu_vs_Kendall's Tau.pdf}
\includegraphics[width=4.3 cm]{./Results/Cause4copula_results/Cause4copula_stat_results_ar1_beta_0.5_Clayton/Clayton_Pl-Pu_vs_Kendall's Tau.pdf}\\
\vspace{-0.2 cm}
\includegraphics[width=4.3 cm]{./Results/Cause4copula_results/Cause4copula_stat_results_ar1_beta_0.5_Clayton/Clayton_D2u_vs_Kendall's Tau.pdf}
\includegraphics[width=4.3 cm]{./Results/Cause4copula_results/Cause4copula_stat_results_ar1_beta_0.5_Clayton/Clayton_D2l_vs_Kendall's Tau.pdf}
\includegraphics[width=4.3 cm]{./Results/Cause4copula_results/Cause4copula_stat_results_ar1_beta_0.5_Clayton/Clayton_D2u-D2l_vs_Kendall's Tau.pdf} \\
\hspace{-0.6 cm}
\includegraphics[width=4.3 cm]{./Results/Cause4copula_results/Cause4copula_stat_results_ar1_beta_0.5_Clayton/Clayton_scatter_CorlmCoru_Kendall's Tau.pdf}
\includegraphics[width=4.3 cm]{./Results/Cause4copula_results/Cause4copula_stat_results_ar1_beta_0.5_Clayton/Clayton_scatter_PlmPu_Kendall's Tau.pdf}
\includegraphics[width=4.3 cm]{./Results/Cause4copula_results/Cause4copula_stat_results_ar1_beta_0.5_Clayton/Clayton_scatter_D2umD2l_Kendall's Tau.pdf}
\caption{Comparison of correlation and tail association statistics between 
environmental-noise inputs and ecological-variable 
outputs, AR(1) model, Clayton copula. Grey and black points give 
means over 50 replicate simulations of the listed statistics, 
grey points are for noise inputs and black points are for
model outputs. Error bars are standard errors. 
$p$-values (triangles, right axis) are for a 
paired $t$-test of the null hypothesis that the distributions have the same mean. 
For bottom panels, headers give Pearson correlations for
the points. The regression line through the points (solid) was similar
to the 1-1 line (dashed). 
See section \ref{Causes} for details.\label{fig_Cause4copula_C_kend_ar1}}
\end{center}
\end{figure}

For our evolutionary model, the hypothesis was 
correct, for the models we employed, that asymmetric tail association in evolutionary changes 
can produce similarly asymmetric tail association between characters across
phylogeny tips.
Once character evolution was simulated 100 times for each of the
dependence structures we considered for evolutionary changes, 
we had 817 bivariate characters for each of 500 simulations (see section 
\ref{Causes}). We 
computed our nonparametric asymmetry statistics 
$\cor_{0,0.2}-\cor_{0.8,1}$, $\Ps_{0,0.2}-\Ps_{0.8,1}$,
and $\Dsq_{0.8,1}-\Dsq_{0,0.2}$ for each simulation output and 
produced a histogram for each statistic and for each
dependence structure (Fig. \ref{fig_cop_evo}). Results showed that 
asymmetries of tail associations in evolutionary changes
were associated with similar asymmetries of tail associations in extant characters.
Thus we cannot reject the possibility that
this is a proximate mechanism behind observed asymmetric tail associations
in our bird and mammal data.
This result is represented in Fig. \ref{MasterFigure} as the dashed box around 
"If character evolution has copula structure" and the arrow labelled "D".
This topic is revisited in the Discussion, as is the box "Asymmetric species
interactions" and the arrow labelled "C" in Fig. \ref{MasterFigure}.

<!--results for third mechanism-->
\begin{figure}[!h]
\begin{center}
\includegraphics[width=15 cm]{./Results/BMR_results/evo_cop/singlepage_hist_taildep.pdf}\\
\caption[Non-parametric statistics on a mechanism of tail dependence in 
character distributions]{Three measures of asymmetry of tail 
associations ($\cor_{0,0.2}-\cor_{0.8,1}$, A-E; 
$\Ps_{0,0.2}-\Ps_{0.8,1}$, F-J; $\Dsq_{0.8,1}-\Dsq_{0,0.2}$, K-O) between two characters, across extant species,  were 
computed for each of 100 simulations of mammalian character evolution for each of five types of tail dependence between 
evolutionary changes in the characters (extreme left-tail dependence, A, F, K; moderate left-tail dependence, B, G, L; 
symmetry of tail dependence, C, H, M; moderate right-tail dependence, D, I, N; and extreme right-tail dependence, 
E, J, O). See section \ref{Causes} for details. The number above each panel indicates how many statistics out of $100$ were less than or equal to $0$. Values substantially lower (respectively, higher) than $50$ indicate greater left (respectively, right) tail-dependence between characters.\label{fig_cop_evo}}
\end{center}
\end{figure}

<!--Question 3 analyses-->
# Concepts and methods for Q3 \label{Consequences}

Recall that Q3 is: what are the consequences of non-normal copula
structure and asymmetric tail associations for ecological understanding and
applications?
We addressed Q3 by exploring, using both data and models, three hypothesized
consequences. As for Q2, the models we used are intended as initial explorations, only. 
Therefore only simple models are used, and comprehensive explorations of the 
sensitivity of results to model structure and parameters are left for future work.

The hypothesis was presented in the Introduction that
the distribution (through time) of a spatially 
averaged quantity should be influenced by 
dependencies between the local quantities being averaged, including their
copula structure and tail associations. 
To make this hypothesis more precise, suppose an ecological variable 
$E_i(t)$ is measured at locations $i=1,\ldots,N$ and times $t=1,\ldots,T$, and
the spatial mean $\sum_i E_i(t)/N$ is of
interest. The $E_i(t)$ could be, for instance, local abundances of a pest or 
exploited species, or local fluxes of a greenhouse gas. If  $E_i$ and $E_j$ are 
associated primarily in their right tails for most location pairs $i$ and $j$, then
exceptionally large values tend to occur at the same time in most locations.
We hypothesize that this can increase the skewness of the distribution 
of the spatial mean. Similarly, left-tail associations between local variables
should decrease skewness. Strong positive skewness of the 
spatial-mean time series corresponds to "spikiness" of that time series,
i.e., occasional very large values, 
which corresponds to instability through time. 
The spatial-mean time series and its skewness may be quantities of 
principal importance for 
pest or resource abundance, for which extreme values (spikes) 
in the spatial mean may have large effects.

We tested the above hypothesis
using our multivariate datasets (Table \ref{tab_data_info}).
For each dataset we calculated the spatial mean time series, and then 
the skewness through time of that
mean. Then, for each dataset, we compared the value obtained 
to a distribution of values 
of the same quantity for each of 10000 surrogate datasets. Surrogate 
datasets were produced by randomizing the empirical data
in a special way to have the 
copula structure of a multivariate-normal distribution, but 
to retain exactly the same
distributions of values for each sampling location as the original data
(Appendix \ref{SM-Surrog_multi}). 
Surrogates also had very similar Spearman correlations 
between pairs of sampling locations
as the data. 
Our comparisons therefore tested the null 
hypothesis that the skewness of the spatial 
mean took values on the empirical data no different 
from what one would expect if the copula
structure of the data were multivariate normal (i.e., the same copula as a 
multivariate normal distribution), but the data were otherwise statistically
unchanged. Significant differences indicate that non-normal copula structure in the data 
contributed to the skewness of the spatial mean time series, 
i.e., to its instability and "spikiness" through time.

For green spruce aphid abundance data, *C. furca* abundance data, and methane-flux
data, because these datasets exhibited stronger lower- than upper-tail associations 
(Table \ref{tab_multivar_summary}), we compared empirical and surrogate skewness values
via a one-tailed test in the left tail: the $p$-value was the fraction
of surrogate skewnesses less than the skewness for the empirical data. The test
examines whether stronger lower-tail associations between local time series 
caused the spatial average to have significantly lower
skewness than would have been expected with symmetric tail associations. For
leaf-curling plum aphid data, because that dataset exhibited stronger upper-
than lower-tail 
associations (Table \ref{tab_multivar_summary}), we did the analogous one-tailed test
in the right tail. The test examines whether upper-tail associations caused the spatial
average to have significantly higher skewness than 
would have been expected with symmetric tail associations.

We also examined a hypothesis that 
asymmetric tail associations across space of an environmental variable can influence the 
extinction risk of a metapopulation. We hypothesized that environmental noises
exhibiting greater left-tail associations across habitat patches would 
cause higher metapopulation extinction risks because then very bad years for the 
component populations occur simultaneously in many patches,
reducing rescue effects. Here we assume, for simplicity, 
that low values of the environmental
variable are "bad" for the populations and high values are "good".
We tested the reasonableness of this hypothesis using a metapopulation
extension of the Lewontin-Cohen model,
$\vec{P}(t+1) = D \lambda(t) \vec{P}(t)$,
where the $i^{th}$ component of the length-$N$ vector 
$\vec{P}(t)$ represents
population density in the $i^{th}$ habitat patch at time $t$. 
The $N \times N$ matrix $\lambda(t)$ was diagonal with $i^{th}$ diagonal entry 
$\exp(r+\varepsilon_i(t))$. Here $r$ is a growth rate; we used $r=0$. The
$\varepsilon_i(t)$ represent environmental noises. They were standard-normally 
distributed, were independent through time, and showed the same 
spatial correlations for every simulation, 
but were made to exhibit stronger right- or left-tail associations between patches in
different simulations (Appendix \ref{SM-ext_risk_model}). The $N \times N$
matrix $D$ was a dispersal matrix
modelling local or global dispersal at rate $d$, in different simulations
(Appendix \ref{SM-ext_risk_model}).
After each step, 
if the density in a patch was $<1$, it was
set to $0$. We simulated the model 10000 times for each combination of parameters,
starting from $p_0 = 50$ in each patch, 
and calculated extinction risk after $25$ time steps.

Finally, and pursuing ideas from @Reuman2017, we 
tested whether the copula structure of the dependences between population variables 
measured in different locations has consequences for the spatial version of Taylor's
law. Taylor's law is a commonly observed and widely studied 
[@Taylor1961; @Taylor1988; @Cohen2013; @Xu2016; @Reuman2017] 
empirical pattern that relates the 
variances of groups, $g$, of population measurements to the means of the 
groups via a power law, $v_g = a \times m_g^b$, or equivalently, 
$\log(v_g)=\log(a)+b \times \log(m_g)$. Here $b$ is called 
the Taylor's law *exponent* or 
*slope*, and $\log(a)$ is the Taylor's law intercept. 
There are several versions of Taylor's law. For spatial Taylor's law,
given population density or abundance data $x_i(t)$ measured in locations
$i=1,\ldots,N$ and times $t=1,\ldots,T$, a group $g$ consists of 
all the measurements $x_i(t)$, $i=1,\ldots,N$ made in different locations at the same time.
So means $m_g$ and variances $v_g$ are computed across space. We refer to the matrix
with $x_i(t)$ in the $t^{\text{th}}$ row and the $i^{\text{th}}$ 
column as a *population matrix*. 
We consider that 
Taylor's law holds true for a dataset if the 
$\log(v_g)$ versus $\log(m_g)$ scatter plot 
for the dataset shows a linear, homoskedastic pattern.
Taylor's law has been 
verified empirically for hundreds of taxa and has been applied in numerous fields
including fisheries management, estimation of species persistence times, and
agriculture [@Cohen2015; @Reuman2017], so the question of whether copula
structure of spatial dependence influences Taylor's law may have applied
significance.

To explore the influence of copula structure and tail associations on spatial Taylor's
law, we carried out a series of simulations that generated sets of population 
matrices that had 
distinct copula structure between locations but that 
were otherwise statistically similar. For $k=1,\ldots,1000$
and X representing a Clayton, normal, Frank, or a survival Clayton copula family,
we generated $N \times n$ population matrices $m^{(\text{X},k)}$ with the following properties:
1) for any given values of $j$ and $k$, the values 
$\{ m_{t,j}^{(\text{X},k)} : t=1,\ldots,N  \}$ were exactly the same, as unordered
sets, for all X, i.e., the same actual population values were used for location $j$ 
in simulation $k$, regardless of X;
2) for any given $k$, and $j_1$ and $j_2$ such that $j_1 \neq j_2$, the 
Spearman correlations (computed through time)
$\cor_t(m_{t,j_1}^{(\text{X},k)},m_{t,j_2}^{(\text{X},k)})$ were the same, to within sampling
variation, i.e., correlations through time of populations in two locations were the same,
up to sampling variation, regardless of X;
3) the copula for the dependence between $m_{t,j_1}^{(\text{X},k)}$ and $m_{t,j_2}^{(\text{X},k)}$ was  from the family X.
We used $N=50$ and $n=25$. Details of how these matrices were generated are in 
Appendices \ref{SM-Consequence_TL_simulation} and \ref{SM-sect:copsurrognd}.
Thus for each $k$, the population matrices $m^{(\text{X},k)}$
for X taking the values Clayton, normal, Frank, and survival Clayton
were statistically similar except for the copula structure between locations, which was
X. Thus comparing how the $\log(v_g)$ versus $\log(m_g)$ relationship may
manifest differently for $m^{(\text{X},k)}$
for different values of X constitutes a test of the influence of copula
structure on Taylor's law. We selected X equal to the Clayton, normal, Frank, and survival 
Clayton copulas to explore a range of tail association patterns.  

A variety of Taylor's law statistics were computed for each population matrix
$m^{(\text{X},k)}$. First, the $N=50$ 
spatial means,
$m_g$, and variances, $v_g$, were computed, and $\log(v_g)$ vs. $\log(m_g)$ plots 
were considered. Linearity of the $\log(v_g)$ vs. $\log(m_g)$ relationship was 
tested for each simulation by comparing the linear regression
through these $50$ points to a quadratic alternative via an $F$-test, 
producing a $p$-value result. If 
these $p$-values were uniformly distributed across the unit interval for the 1000 replicate 
simulations which were generated for X, it supported the 
linearity assumption of Taylor's law for X; whereas if 
they were clustered toward smaller values
the test tended to reject that assumption. We also tested, for each simulation, 
the assumption of Taylor's 
law that the $\log(v_g)$ vs. $\log(m_g)$ plot was homoskedastic: we regressed the absolute 
residuals of the 
linear regression of $\log(v_g)$ versus $\log(m_g)$ against the predictions of that 
regression. A significant $p$-value
result of this test indicates heteroskedasticity. For each simulation we also computed the
root mean squared error of $\log(v_g)$ vs. $\log(m_g)$ data from the linear regression
of $\log(v_g)$ versus $\log(m_g)$, as well as the 
intercept and slope of the linear regression.
Finally we recorded the quadratic coefficient of the regression of $\log(v_g)$ against
$\log(m_g)$ and $(\log(m_g))^2$, and the mean curvature of 
the quadratic regression equation 
across the values $\log(m_g)$. The distributions of all these 
statistics across the 1000 replicate 
simulations were compared for X the Clayton, normal, Frank, and survival Clayton copulas, 
to determine if copula structure influences Taylor's law.

# Results for Q3: Tail dependence influences skewness of the spatial average, extinction risk, and Taylor's law \label{Results_Q3}

\noindent Empirical results were consistent with the hypothesis that the skewness
of a spatial-average time series is influenced by tail associations
between the local quantities being averaged. 
For datasets that
exhibited highly asymmetric tail associations in earlier analyses (green spruce 
aphid abundance data had stronger lower-tail association, and, 
respectively,  leaf-curling plum 
aphid first flight data had stronger upper-tail association, Table 
\ref{tab_multivar_summary}), skewness of the spatial average 
was less than (respectively, greater than) a significant fraction of 
surrogate skewness values 
(Fig. \ref{fig_skewness_spearman}A, B). For datasets with moderately stronger 
lower than upper-tail 
dependence (\emph{Ceratium furca} abundance and methane data), 
skewness of the spatial
average showed a non-significant or marginally significant tendency 
toward being less than surrogate skewnesses (Fig. \ref{fig_skewness_spearman}C, D). 
Thus copula 
structure and asymmetric tail associations are important for spatially averaged 
quantities and can influence the temporal stability of those quantities. 
This result is represented in Fig. \ref{MasterFigure}
as the solid box around "Instability/skewness of mean or total 
time series" and the solid arrows labeled "X". 

\begin{figure}[!h]
\textbf{ \hspace{5.6 cm} (A) \hspace{4 cm} (B)} \\
\vspace{-0.6 cm}
\begin{center}
\includegraphics[width=5 cm]{./Results/skewness_results/skewness_aphid_count/skewness_aphid_count_sp_10_spearman.pdf}
\includegraphics[width=5 cm]{./Results/skewness_results/skewness_aphid_ff/skewness_aphid_ff_sp_11_spearman.pdf}\\
\end{center}
\vspace{0.2 cm}
\textbf{ \hspace{5.6 cm} (C) \hspace{4 cm} (D)} \\
\vspace{-0.6 cm}
\begin{center}
\includegraphics[width=5 cm]{./Results/skewness_results/skewness_plankton_north_sea/skewness_plankton_north_sea_sp_16_spearman.pdf}
\includegraphics[width=5 cm]{./Results/skewness_results/skewness_methane/skewness_methane_spearman.pdf}
\caption{Skewness of spatially averaged green spruce aphid abundance (A), leaf-curling plum 
aphid first flight dates (B), \emph{Ceratium furca} abundance (C), and methane-flux
(D) time series compared to a multivariate 
normal-copula null hypothesis. Black dots are empirical 
skewnesses; see text 
and Appendix \ref{SM-Surrog_multi} for details of the null hypothesis. Results show a 
tendency for skewness of the spatial average to be affected as hypothesized 
by asymmetric tail associations. \label{fig_skewness_spearman}}
\end{center}
\end{figure}

Consistent with our extinction risk hypothesis, 
left-tail-associated environmental fluctuations increased metapopulation
extinction risk for the spatial 
Lewontin-Cohen model, 
for $N=5$ and $N=25$ habitat patches, and for local 
and global
dispersal (Fig. \ref{fig_extrisk_lcohen_local}). This result 
is represented in Fig. \ref{MasterFigure}
as the dashed box around "Extinction risk" and the dashed arrow labelled "Y".

Also consistent with hypothesis, copula structure had a substantial effect on
Taylor's law for the models we considered. 
Taylor's law was strongly influenced, and was often even invalidated in its assumptions of 
linearity and homoskedasticity, by non-normal copula structure. For normal copula structure 
(i.e., for simulations that gave normal-copula dependence between populations
in different locations), $p$-values for the linearity
and homoskedasticity tests were roughly uniformly distributed across replicate 
simulations and Taylor's law appeared visually to be a reasonable approximation
of the $\log(v_g)$ versus $\log(m_g)$ relationship (Fig. \ref{fig_spTL}B, E, F). 
Furthermore, quadratic coefficients
and curvature values were close to $0$ (Fig. \ref{fig_spTL}J, K), and root mean squared
errors from the linear regression were relatively small (Fig. \ref{fig_spTL}G).
But linearity or homoskedasticity were violated more frequently for non-normal 
copula structure
(Fig. \ref{fig_spTL}A, C, D, E, F); quadratic coefficients and curvatures were frequently 
non-zero (Fig. \ref{fig_spTL}J, K); and root mean squared
errors from the linear regression were much higher (Fig. \ref{fig_spTL}G). Slopes 
and intercepts of the linear regression were also strongly affected by copula
structure (Fig. \ref{fig_spTL}H, I), though some of the effect here was because
linear regressions do not always adequately represent the $\log(v_g)$ versus $\log(m_g)$
relationship when copula structure was not normal.
Thus our results substantiated the hypothesis, at least for the models we used,
that Taylor's law can be influenced by copula structure and asymmetric
tail associations. This is represented in Fig. \ref{MasterFigure} by the 
arrow labelled "Z".

\begin{figure}[!h]
\textbf{ \hspace{5.6 cm} (A) \hspace{4 cm} (B)} \\
\vspace{-1 cm}
\begin{center}
\includegraphics[width=10 cm]{./Results/ext_risk_copula/lcohen_ext_risk_local_disp.pdf}
\end{center}
\textbf{ \hspace{5.6 cm} (C) \hspace{4 cm} (D)} \\
\vspace{-1 cm}
\begin{center}
\includegraphics[width=10 cm]{./Results/ext_risk_copula/lcohen_ext_risk_global_disp.pdf}
\end{center}
\caption{Extinction risk for the metapopulation extension of the Lewontin-Cohen model, after 
$25$ time steps, was higher for environmental noise with stronger left-tail association
across space. Dispersal was local (A, B) or global (C, D),
with dispersal rate $d$ (Appendix \ref{SM-ext_risk_model}). Simulations used $N$ patches
for $N=5$ (A, C) and $N=25$ (B, D).\label{fig_extrisk_lcohen_local}}
\end{figure}

<!--adding plots showing influences of non-normal copula on spatial Taylors law-->
\begin{figure}[!h]
\textbf{ \hspace{2cm} (A) \hspace{2.92 cm} (B) \hspace{2.92 cm} (C) \hspace{2.92 cm} (D)} \\
\vspace{-0.8 cm} 
\begin{center}
\includegraphics[width=16 cm]{./Results/taylorslaw_results/SpatialTL_CNFSC_spearman.png}\\
\end{center}
\textbf{ \hspace{1 cm} (E) \hspace{1.4 cm} (F) \hspace{1.2 cm} (G) \hspace{1.2 cm} (H) \hspace{1.4 cm} (I) \hspace{1.4 cm} (J) \hspace{1.2 cm} (K)}\\
\vspace{-0.8 cm}
\begin{center}
\includegraphics[width=16 cm]{./Results/taylorslaw_results/SpatialTL_7stats_spearman.pdf}
\caption[Spatial Taylor's law results]{Spatial Taylor's law results. Log(spatial variance) vs. 
log(spatial mean) relationships for $1000$ simulations over $50$ time steps for $25$ locations 
using Clayton (A), normal (B), Frank (C), and survival Clayton (D) copula structures between
locations. Results from all $1000$ 
simulations were plotted on the same axes for each panel, but two example simulations are 
shown in solid and hollow black circles to help assess to what extent variation on the plots was between 
simulations or within simulations. For each simulation, we tested the linearity 
(E) and homoskedasticity (F) of the log(variance) vs. log(mean) plot for that simulation, 
and quantified the root mean squared error of points from the linear regression line (G).
We also calculated the linear regression intercept (H) and slope (I), the quadratic
term of a quadratic regression through the points on the log(variance) vs. log(mean) plot
(J), and the curvature of that quadratic regression (K). Distributions of values 
across all 1000 simulations are displayed. Spearman-preserving surrogates were used, 
though results using Kendall-preserving surrogates were similar. See section \ref{Consequences}
and Appendices \ref{SM-Consequence_TL_simulation} and 
\ref{SM-sect:copsurrognd} for details.\label{fig_spTL}}
\end{center}
\end{figure}

# Discussion \label{Discussion}

<!--what we did and why is it important to study?-->
We showed that non-normal copula structure and asymmetric tail associations are 
common across multiple sub-disciplines
in ecology, although these facets of data are only occasionally 
accounted for [@Valpine2014; @Anderson2018; @Popovic2019]. We hypothesized 
mechanisms that may cause non-normal copula structures and asymmetric tail associations;
we discuss below how commonly some of our
mechanisms may operate. 
We also demonstrated important consequences of 
non-normal copula structure and asymmetric tail associations for ecology. 
For instance, the skewness
of a spatial-average time series is influenced by asymmetric tail associations between its 
constituent time series: predominantly right-tail-associated local time series 
can lead to "spiky" spatially averaged time series, with large outbreaks; and 
predominantly left-tail-associated local time series can lead to spatially
averaged time series showing accentuated "crashes".
Thus tail associations could have implications for pests and exploited species.
Extinction risk and Taylor's law can also be altered by tail association 
patterns across space. 
In our view, our results make it reasonable to suggest 
that a more comprehensive understanding 
of many ecological phenomena may be possible if a
complete copula characterization of associations between variables is employed. 
Copula statistics are
well developed [@nelsen2006_copula; @joe2014_dependence; @MaiScherer2017], 
and have been introduced in accessible formats [@Anderson2018; @Genest2007]. Furthermore, 
open-source computer implementations exist (e.g., the `copula` and 
`VineCopula` packages in R). Ecologists can apply these tools
immediately. We created several interrelated randomization procedures
(Appendices \ref{SM-surrog_test}, \ref{SM-Surrog_multi}, \ref{SM-sect:copsurrognd}) 
that built upon existing copula methods. 

<!--tropical ecosystems-->
The approaches we demonstrated should apply equally well to 
data from tropical or temperate ecosystems.
There is no reason to expect that non-normal copula structure and 
asymmetric tail associations should be special properties of datasets from temperate 
regions. The mechanisms we proposed of non-normal copula structure and asymmetric tail 
associations seem equally likely to apply anywhere. 
For instance, the Moran effect, which underpins 
one of our proposed mechanisms (Fig. \ref{MasterFigure}, B), is a standard mechanism
that occurs whenever environmental variables influence populations. And Liebig's
law and nonlinear environmental influences on ecosystems, which underpin 
another one of our mechanisms (Fig. \ref{MasterFigure}, A), are widely
demonstrated phenomena. 

<!--Stuff about the mechanisms-->

<!--a para on how likely one is to get copula structure via the Liebig mechanism-->
Our first proposed causal mechanism (Fig. \ref{MasterFigure}, A) may well operate commonly,
for two reasons.
First, Liebig's law and the idea of limiting 
nutrients are dominant paradigms in ecology, and 
many studies have documented nonlinear or 
threshold influences of environmental
variables on ecological quantities. 
Second, fluctuations in environmental variables through time
are very commonly correlated across space. 
Because these factors, which are the essential ingredients of
the mechanism, are common, it is reasonable to hypothesize that the mechanism
may operate commonly and may be a dominant cause of 
asymmetric tail associations and non-normal copula structure of ecological 
dependencies across space.
We provide further support for the mechanism 
in our discussion of green spruce aphids and
winter temperature below.

<!--para on how likely one is to get copula structure via the Moran mechanism-->
There are also reasons to hypothesize that 
our Moran mechanism (Fig. \ref{MasterFigure}, B) 
may operate commonly: Moran effects are common [@Liebhold2004; @Sheppard2016; 
@Defriez2017_ocean; @Defriez2017_land], and 
non-normal copula structures and asymmetric tail associations are often found in 
environmental variables. If intense meteorological events are 
also widespread, then
environmental variables associated with these events should take extreme
values simultaneously across large spatial areas, producing tail associations in measurements
made through time at different locations. Non-extreme values may instead
be associated with local phenomena, and therefore may be
less correlated across large areas.
@Serinaldi2008 examined the spatial dependence 
of rainfall in Central Italy. Gumbel or Student 2-copulas 
were candidates for modelling dependence, and 
neither of these is a normal copula.
A long-term study (1950–2014) in the Loess Plateau of China [@She2018]
showed that a Gumbel copula effectively modeled 
the spatial dependence of drought variables. The Gumbel copula has 
asymmetric tail association.
Bivariate copula analysis was also used in forecasting the 
co-occurrence of extreme events (flood or drought) 
over the North Sikkim Himalayas using spatial datasets [@Goswami2018]. 

<!--Causal mechanism C, which we also did not test-->
We suggested in the Introduction that asymmetric competitive relationships
between species could yield asymmetric tail associations between abundance measurements for the 
species. This is another theoretical mechanism for non-normal copula
structure and asymmetric tail associations, represented in Fig. \ref{MasterFigure} 
by the box "Asymmetric species interactions" and the arrow labelled "C".
It could be tested by analyzing copulas of abundances of competing 
species, sampled across space or time. We note that, whereas all the datasets
we studied here have been positively associated when they were significantly
associated, for negatively associated variables such as abundances of competing species,
the definitions of left- and
right-tail association no longer apply, strictly speaking: the left tail
of one distribution corresponds to the right tail of the other. One must
be careful with terminology, but it is still possible to study asymmetries of
association.

<!--Holder stuff-->
Our simulations of character evolution suggest the hypothesis that
changes through evolutionary time in bird and 
mammal BMR and body size may exhibit greater right- than left-tail association,
contrary to standard normality assumptions of character evolution models. This is
a hypothesis only, because the greater right-tail association shown in 
Fig. \ref{fig_biv_multi_raw_cop}F, G
could have come about in another, unknown way rather than via the mechanism we suggested
which implicated asymmetric tail associations in evolutionary change. 
Our simulations show that asymmetric tail associations 
in evolutionary changes are sufficient, but may not be necessary, to produce 
the observed asymmetric tail
associations in characters of extant species. For instance (see below, and Appendix
\ref{SM-sect:missingdata} and Fig. \ref{SM-fig:missingdata}), 
systematically missing data can also produce tail dependence and may
have influenced results for the BMR-body mass datasets. Even if the hypothesized 
evolutionary mechanism (Fig. \ref{MasterFigure}, D)
is correct, our results only replace one question, i.e., why do we see greater right-
than left-tail
association between BMR and body mass, with another,
i.e., why might we see greater right- than left-tail association in evolutionary 
changes in these traits?

Performing statistical tests for associations between continuous traits across different species
    was a primary motivating example of phylogenetic comparative methods.
For example, Felsenstein's method of phylogenetically independent contrasts 
    [@Felsenstein1985] is currently the second most-cited paper in the history
    of The American Naturalist journal [@HueyGT2019].
Frustratingly, the field still lacks a dependable procedure for dealing with branch lengths, 
    which are a crucial input to the method.
Ideally, the branch lengths used to correct for phylogenetic effects would represent
    the expected amounts of change for the characters that are being analyzed.
Because researchers almost never have a reliable method for providing such branch lengths,
    most researchers rely on ultrametric trees -- those for which 
    the duration of the branch in time can be treated as a proxy for the branch length.
Frequently these branch lengths are transformed to assess sensitivity to different 
    assumptions about the degree of phylogenetic inertia displayed by the traits
    under study [@Ives2018; @Harmon2018].
Even if one were able to simply use a time-based set of branch lengths, assigning
    dates to nodes in phylogenies is difficult.
DNA sequence data can provide estimates of branch lengths, but these estimates
    are dependent on the adequacy of models which correct sequences for multiple
    substitutions occurring at the same location.
Biases in estimating the evolutionary distance can affect downstream analyses [@Phillips2009].
Additionally, changes
    in the rate of molecular evolution make the estimation of dates difficult [@HeathM2014] even when branch lengths in time are accurately estimated.

Without reliable branch length estimates, it is difficult to interpret the significance
    of the magnitude of changes in traits across a tree.
Developing tests of association based on copula structure may possibly lead the way to more
    robust methods for studying associations when we lack defensible estimates of
    branch lengths. 
We note that our phylogenetic analyses here consisted merely of simulations to assess
    whether interesting copula structure in a simple evolutionary process could 
    leave a detectable signal on the trait data for extant species.
Substantial work remains to be done before we have a copula-based method for 
    analyzing data on a phylogenetic tree. 
But the nonparametric nature of approaches
    based solely on data ranks, as many of our methods are, may be a promising avenue for 
    avoiding inaccuracies that arise via the highly structured model assumptions implicit
    in the method of phylogenetically independent contrasts and related methods.

Additionally, character evolution was simulated using one random draw from the relevant matrix
per phylogeny branch. This was principally because branch lengths are often hugely
uncertain. If branch lengths were well known, alternative simulation strategies may
include selecting one draw from the matrix of evolutionary changes per unit time, 
or selecting one draw for the whole branch but rescaling the variances of the selected
character changes according to the 
length of the branch. These choices amount to the same thing for normal copula
structure, but not otherwise. Modelling choices such as these may have influenced
our results. Additional research seems warranted testing the realism 
of our hypothesized mechanism and simulations.

<!--Discuss stuff here on tradeoff and life history that Bever told me about-->
Relationships between BMR and body mass relate
to a trade-off between mass-specific BMR (BMR per unit body mass) and body
mass itself. Copulas probably interrelate with 
life-history trade-off theory in additional ways beyond what we demonstrated.
For instance, it is well known that energy allocation to a life function, F (e.g., reproduction)
will reduce the energy that can be allocated to other functions, $G_1$, $G_2$, $G_3$ 
(e.g., growth, predation avoidance). This is the principle of 
allocation. But F can trade off against any or all of 
the $G_i$. Therefore, for large F, approaching absolute 
limitations, there may be a strong 
association between F and $G_1$, for instance. For small F, there may be
little association because resources not allocated to F can instead be allocated 
to any combination of the $G_i$. This constitutes 
asymmetric tail association between F and $G_1$. @Winemiller1992
described a three-way trade-off in fishes between age of reproductive
maturity, juvenile survivorship, and fecundity. The trade-off should, in theory,
produce a tight association between age of maturity and fecundity for 
fishes with low age of maturity, but little such association for later-maturing
fishes because those species may invest the resources not invested in maturing quickly into
either fecundity or juvenile survival.
These ideas suggest that copulas may be 
useful for studying multi-dimensional life-history trade-offs. But applications will 
require careful attention to the possible consequences of biased sampling:
if the degree of completeness of a dataset is associated with one or
more of the characters, then statistical artefacts
can bias conclusions (Appendix \ref{SM-sect:missingdata}, 
Fig. \ref{SM-fig:missingdata}). Multi-dimensional copulas may also be the 
appropriate copula approach
to studying multi-dimensional trade-offs. We used only bivariate 
analyses in this study solely for simplicity.
But statistical theory on multi-dimensional copulas is also well developed
[@nelsen2006_copula; @joe2014_dependence; @MaiScherer2017], 
and open-source computer implementations exist for multivariate as well as 
bivariate copula methods (e.g., the `VineCopula`
package for R). Such approaches may be a 
useful next step in life-history theory. 

<!-- We have only given a few mechanisms, but there are likely to be others, and discoverig
additional mechanisms of copula structure is certainly worthy of future work. -->
Additional mechanisms of non-normal copula structure and asymmetric tail associations
probably also operate. For instance, 
measurement error may modify copula structure. 
Our models investigating potential causes of copula structure 
were intentionally simple in other respects, 
too, not including factors such as
delayed density dependence, dispersal, population stage structure, trophic
interactions, etc.; and we did not comprehensively 
explore parameter space for our models. We re-emphasize 
that fuller explorations, in future work,
of some of our models and of variant models may be informative. 
We hope by enumerating a few potential mechanisms of copula structure
we will inspire additional research on the
potentially numerous mechanisms that may operate in diverse datasets,
and their relative importance under different circumstances.

<!--some consequences paragraphs next-->

<!--a para on stability of the total biomass time series in a community, and skewness - sets up Shya's later paper as future work-->
We also elaborated potential consequences of copula structure for ecological 
phenomena and understanding. 
Our results showed that the skewness of the 
spatial average of local time series is influenced 
by their tail associations. But the same logic should also apply to any collection of time
series, whether associated with locations in space or not. 
Another potential application is time series of abundances of 
all species from a single community, e.g., all
plants in a quadrat surveyed repeatedly over time. A large 
literature has focused on synchrony versus 
compensatory dynamics between such time series, and 
the influence of interspecific relationships on the 
variability of community
properties such as total biomass (e.g., @Doak1998, @Tilman1999, @Tilman2006,
@Gonzalez2009). Typically, variability of community biomass 
is measured with the coefficient of variation, but skewness
may also be of interest because it can help characterize "spikiness" 
through time. Future work on copula 
structure of interspecific relationships in communities and its implications for community 
variability seems likely to be valuable.

<!--A para that sets up Shya's extinction risk paper-->
Although we demonstrated that tail associations between environmental variables
can influence extinction risk, substantial work remains to determine the importance 
of this effect. First, we used 
a non-density-dependent model. Do similar results pertain
when density dependence
is involved? Second, we considered metapopulation extinction risk, but the large field 
of population viability analysis (PVA) via stochastic 
matrix modelling [@Caswell2000; @MorrisDoak2002]
uses a framework in which a single population's 
vital rates (e.g., life-stage-specific 
fecundity and survival rates) are considered to
vary stochastically through time due to environmental variation. Do relationships
between different vital rates exhibit asymmetric tail associations, 
and do tail associations influence extinction risk in this context?
Finally, is the copula structure of environmental variables or vital rates
changing through time, and, if so, how do such changes influence extinction risks?
Climate change is known to amplify the factors that lead to extreme weather events
[@Hansen2012] and hence may alter spatial tail associations for weather variables.

<!--green spruce aphid all the way from causes to consequences-->
Our hypotheses and results cover the presence, causes, and consequences
of non-normal copula structure and asymmetric tail associations in ecological 
systems (Fig. \ref{MasterFigure}), but this was done using a variety of 
datasets and models.
We here take a closer look at the green spruce aphid, because it simultaneously 
illustrates causes and consequences within one system.
Green spruce aphid abundance, as measured in the data we use,
is strongly positively associated with the temperature of the previous winter (@Sheppard2016,
their supplementary figure 6).
Winter temperature for year $t$ was here taken to be an average for December of year $t-1$ through
March of year $t$, was available for the locations of aphid sampling, 
and was preprocessed in the same way as @Sheppard2016.
For each of our 10 sampling locations, we therefore 
examined the copula of winter temperature
and aphid abundance time series for the location, finding 
stronger left- than right-tail associations in 7 of the 8 locations for which 
independence of winter temperature and aphid abundance could be rejected, according to 
the data specific to the location (Table \ref{tab:greenspruce}).
Apparently winter temperature has an asymmetric
influence on aphid abundance in that cold winters generally produce low abundances
but warm winters often do not yield higher abundances than moderate winters.
One of our hypothesized mechanisms (Fig. \ref{MasterFigure}, A), 
which our modelling results
supported (section \ref{Results_Q2}), therefore suggests that spatial 
dependence between green spruce aphid counts
in different locations should show stronger left- than right-tail associations. This is 
exactly what was observed (Table \ref{tab_multivar_summary}), providing 
empirical evidence supporting the mechanism (this is why the box around "Nonlinear 
environmental effects, Liebig's law" and some of the arrows labelled "A" are solid
instead of dashed in Fig. \ref{MasterFigure}). The consequences of such tail 
dependence for the skewness of spatially averaged aphid counts was 
described previously 
(Fig. \ref{fig_skewness_spearman}A). Thus the asymmetric influence of winter 
temperature ultimately causes spatially averaged aphid abundance time series 
to have lower skewness (i.e., less spikiness, and greater stability through time) 
than they would otherwise. It seems likely that asymmetric influence of winter 
temperature on populations may be a common phenomenon, so effects such as 
we have documented for green spruce aphid may be common.
 
```{r tab_greenspruceaphid_temp_summary, echo=F, results='asis',message=F}
library(kableExtra)
ms_npa_RES_greenspruceaphid_count_temp<-readRDS("./Results/aphid_count_winter_temp_results/ms_npa_RES_greenspruceaphid_count_temp.RDS")
res<-ms_npa_RES_greenspruceaphid_count_temp$sp10_Green_spruce_aphid$ms_summary_allloc
#row.names(res)<-NULL
#res<-cbind(c(2:11),res)
colnames(res)<-c("p-value, independence test","Best-fit copula",
                 "Best-fit AIC","Normal AIC","p, CvM","p, KS","Model-avg. LT","Model-avg. UT","Model-avg. LT minus UT","cor$_{l}$ - cor$_{u}$","P$_{l}$ - P$_{u}$","D$_{u}^2$ - D$_{l}^2$")
res<-res[,-1]
res<-cbind("Site"=c(2:11),res)
knitr::kable(res,"latex", booktabs = T, linesep = "",digits=3, escape=F,longtable=T,row.names = F,
             caption.short="Summary results for copulas for winter temperature versus green spruce aphid abundance",
      caption = "Summary results for analysis of copulas for winter temperature versus green spruce 
aphid abundance for each of the 10 sampling locations, using the same model selection
and nonparametric methods detailed in section \\ref{Methods}. The column labeled p, CvM is the p-value 
result for the goodness of fit test using the Cramer-von Mises statistic. The 
column labeled p, KS is the p-value 
result for the goodness of fit test using the Kolmogorov-Smirnov statistic. Entries with an NA are because model selection and other statistics
were only employed for copulas for which independence was rejected (5 percent significance level).\\label{tab:greenspruce}") %>% column_spec(1:12,width="1.1 cm")#%>%
#kable_styling(latex_options = c("scale_down"))
```

<!--Alternative methods-->

<!--Partial Person correlation-->
We carried out most of our analyses on data ranks, $u_i$ and $v_i$, for good reasons 
mentioned in sections \ref{Intro}, \ref{Background} and \ref{Methods} and 
reviewed here; but under some circumstances it may also be appropriate to use 
techniques which parallel our nonparametric statistics but that use unranked data, 
$x_i$, $y_i$. As mentioned previously, @Genest2007 and others
recommend carrying out inferences about dependence structures (which was our
goal) using normalized ranks, stating explicitly that "statistical inference concerning
dependence structures should always be based on ranks". The reason for this is
that measures of dependence which use unranked data conflate information about
the marginal distributions of variables with information about the association
between variables. For instance, suppose the true population densities $p_{1,t}$
and $p_{2,t}$ of two species of fish in a lake in year $t$ are unknown,
but are assessed using catch per unit effort (CPUE), a standard approach in fisheries
science [@Zale2013]. The CPUE measurements, which we denote $c_{1,t}$ and $c_{2,t}$, may be 
differently correlated through time than the true densities would be, if they were known,
if the function $f$ relating $p_{i,t}$ and $c_{i,t}$ is nonlinear and if Pearson
correlation is used. The function $f$ modifies the marginal distributions, which modifies
measures such as Pearson correlation that conflate dependence and marginal information.
Rank-based measures such as Spearman or Kendall correlation, or any of our copula approaches,
will be unaffected by nonlinear monotonic functions such as $f$. The same difficulty will pertain in any
case for which a measurement of an ecological quantity is a nonlinear index of the true quantity
of interest. Nevertheless, if one is certain that measurements are linearly related to 
the true quantities being measured, statistics based on unranked data may be more appropriate 
under some circumstances, in part because ecological mechanisms which produce relationships 
between variables are influenced not only by ranks but also 
by the relative-size information in the unranked data. Judgements of what statistics
are appropriate depend, of course, on the purpose of the analysis. The goal of this 
study was to study dependence in isolation, so we used ranks. If one is interested
in tail associations using unranked data, a partial Pearson correlation can be developed
straightforwardly, by replacing $u$ and $v$ in equation \ref{eq:partialspearman}
by $x$ and $y$, but still computing the sum in that equation over data points
with ranks $u_i$ and $v_i$ constrained by the bounds $u_i+v_i>2l_b$ and 
$u_i+v_i<2u_b$. This would result in a different approach to
tail associations that would conflate dependence and marginal information but 
that may also have its own utility for some applications.

<!--referee-requested para comparing to regression with heteroskedastic error term-->
Non-copula methods may be useful under some circumstances for exploring
asymmetries of tail association, but copula approaches have theoretical
and practical advantages. It was suggested to us that asymmetries of
tail association can also be explored using a regression 
model with heteroskedastic error term, $y=mx+b+\varepsilon(x)$. Whereas this 
approach seems capable of yielding insight for some datasets, it implicitly assumes
a particular causal relationship between $x$ and $y$, an untenable assumption for 
many applications. Furthermore, no equivalent
of Sklar's theorem exists under this approach, so information on marginal 
distributions cannot be separated from association information, and the two types
of information will both influence regression results. Sklar's theorem
means that the copula associated with two variables contains all and only the information
about the association. Such statements of 
mathematical completeness will not be available for other methods.

<!--para on data requirements of the methods, requested by the AER referees-->
Data requirements for copula methods will vary depending on multiple factors,
but our results show that requirements are not beyond what is common 
in ecology. With copula methods as with any statistical analysis, the data required to detect 
an effect or a phenomenon depends on the strength of the effect. For instance,
detection of mild asymmetry of tail association will require more data than
detection of strong asymmetry. Different copula or copula-related methods will
also require different amounts of data. For instance, our nonparametric methods 
such as the partial Spearman correlation should be effective with datasets 
only slightly bigger than typical guidelines for standard correlation 
methods; whereas our model-selection approaches may be most effective with more data.
Some of our datasets were large, but others were much smaller, e.g., aphid datasets
comprised 30 points, and our analyses still provided robust results. Thus
many copula methods will be suitable for datasets of sizes that are common in ecology.
Partial Spearman correlation seems particularly widely applicable.

<!--additional copula methods-->
Copula methods are numerous, and go well beyond the cases we have considered
in this introductory paper. We hope our work helps inspire applications of copulas
in ecology both of and 
beyond the specific tools we employed. For instance, multivariate copulas are useful for 
studying relationships between multiple interacting variables, and have been used in finance
to model tail risk for portfolio optimization problems [@Czado2019; @Joe2011]. 
Analogous ecological applications
exist that may be amenable to the same approach, e.g., ecosystem functioning variables
such as total primary productivity or carbon flux are the sum of the contributions of
multiple species in the same way that the value of an investment portfolio is the sum 
of the values of its constituent assets. Ecologists may likewise be interested in
managing the risk that an ecosystem functioning variable will take an extreme value. 
As another example, we assumed for this paper that univariate marginal 
distributions have continuous, strictly 
monotonic cdfs. This was for simplicity and because the simpler
approach was sufficient for our research questions. But count data in ecology have been analyzed
with approaches not making this assumption [@Anderson2018].

# Acknowledgments

We thank Joel E. Cohen for making us aware of copulas in the 
first place. We thank the many contributors to the large datasets 
we used; D. Stevens and P. Verrier for data 
extraction; and David Tilman, Lauren Hallet, Jonathan Walter, Thomas 
Anderson, Lei Zhao, Andrew Rypel, and editors and anonymous reviewers for helpful suggestions. 
We thank James Bell of the Rothamsted Insect Survey (RIS). 
The RIS, a UK Capability, is funded by the Biotechnology and Biological Sciences Research Council 
under the Core Capability Grant BBS/E/C/000J0200. SG, LWS and DCR were partly funded by US National 
Science Foundation grants 1714195 and 1442595, the James S McDonnell Foundation, and the 
California Delta Science Program. 

\clearpage
\newpage

# References

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

<!--<div id="refs"></div>-->