In atredennick/community_synchrony: Calculates syunchrony and stability metrics from time series data
\renewcommand\linenumberfont{\normalfont\tiny\sffamily\color{gray}}
library("methods")
library("knitr")
basename <- "ms_communitySynchrony"
opts_chunk$set(fig.path = paste("components/figure/", basename, "-", sep=""),
cache.path = paste("components/cache/", basename, "/", sep=""))
opts_chunk$set(cache = 2)
opts_chunk$set(tidy=FALSE, warning=FALSE, message=FALSE,
comment = NA, verbose = TRUE, echo=FALSE)
# PDF-based figures
opts_chunk$set(dev='pdf')
\begin{center}
\textbf{\Large{Online Supporting Information: Appendix S1}} \
A.T. Tredennick, C. de Mazancourt, M. Loreau, \& P.B. Adler \
``Environmental responses, not species interactions, determine synchrony...'' \
\emph{Ecology}
\end{center}
\renewcommand{\theequation}{S\arabic{equation}}
\renewcommand{\thetable}{S\arabic{table}}
\renewcommand{\thefigure}{S\arabic{figure}}
\tableofcontents
\listoffigures
\listoftables
\newpage{}
Derivation of limiting case predictions for synchrony of per capita growth rates (Equations 3-4)
Following @Loreau2013 and @DeMazancourt2013a, we define population growth, ignoring observation error, as
\begin{align}
r_i(t) = r_{mi} \left[ 1- \frac{N_i(t)+\sum_{j \neq i} \alpha_{ij}N_j(t)} {K_i} + \sigma_{ei}u_{ei}(t) + \frac{\sigma_{di}u_{di}(t)}{\sqrt{N_i(t)}} \right]
\end{align}
\noindent where $N_i(t)$ is the biomass of species i in year t, and $r_i(t)$ is its growth rate in year t.
$r_{mi}$ is species i's intrinsic rate of increase, $K_i$ is its carrying capacity, and $\alpha_{ij}$ is the interspecific competition coefficient representing the effect of species j on species i.
Environmental stochasticity is incorporated as $\sigma_{ei}u_{ei}(t)$, where $\sigma_{ei}^2$ is the environmental variance and $u_{ei}$ are normal random variables with zero mean and unit variance that are independent through time but may be correlated between species.
Demographic stochasticity arises from variations in births and deaths among individuals (e.g., same states, different fates), and is included in the model as a first order, normal approximation [@Lande2003; @DeMazancourt2013a].
$\sigma_{di}^2$ is the demographic variance and $u_{di}(t)$ are independent normal variables with zero mean and unit variance.
First-order approximations of the temporal variance of total community biomass are obtained as follows [@Ives1995; @Hughes2000; @Ives2002; @Loreau2008a; @DeMazancourt2013a]. Let $\delta N_i(t) = N_i(t) - N^_i$ denote the deviation of observed species i’s biomass from its equilibrium value in the community, $N^i$, in the absence of stochasticity. Equation (S1) can be Taylor expanded around $\delta N_i(t) = u{ei}(t) = u_{di}(t) = u^S_{oi}(t) = 0$ to yield, after dropping terms of order two and higher,
\begin{equation}
\boldsymbol{\delta}\textbf{N}(t+1) = \textbf{A}\boldsymbol{\delta}\textbf{N}(t) + \textbf{z}(t),
\end{equation}
\noindent where $\boldsymbol{\delta}\textbf{N}(t)$ is the vector of deviations of species biomasses from their deterministic equilibrium value, A is the community matrix, also known as the Jacobian matrix around the equilibrium, with elements $(a_{ij})_{1<i,j<S}$:
\begin{align}
A_{ij} =
\begin{cases}
1-r_{mi} \frac{N_{i}^}{K_i}, & i=j \
-r_{mi} \frac{N_{i}^}{K_i} \alpha_{ij}, & i \neq j
\end{cases}
\end{align}
\noindent and $\textbf{z}(t)$ is a vector that encapsulates the effects of environmental and demographic stochasticity whose elements are
\begin{equation}
z_i(t) = N^i \sigma{ei} u_{ei}(t) + \sqrt{N^i} \sigma{di} u_{di}(t)
\end{equation}
When the system reaches a stationary distribution, the variances and covariances between species biomass time series are:
\begin{equation}
\big\langle \boldsymbol{\delta}\textbf{N}(t) \boldsymbol{\delta}\textbf{N}(t)^T \big\rangle = \textbf{C}^\infty = (\text{cov}(N_i,N_j)){1<i,j<S} = (\text{cov}(\delta N_i, \delta N_j)){1<i,j<S}
\end{equation}
\noindent where $\boldsymbol{\delta}\textbf{N}^T$ is the transpose of vector $\boldsymbol{\delta}\textbf{N}$, i.e. a row vector.
Our assumptions listed above lead to the following correlation structure of z:
\begin{equation}
\big\langle \textbf{z}(t) \textbf{z}(t)^T \big\rangle = \textbf{B}
\end{equation}
\begin{equation}
\big\langle \textbf{z}(t-s) \textbf{z}(t)^T \big\rangle = 0 \qquad \text{for} \: s>0
\end{equation}
We use Equation S2 to write a dynamical equation for the covariance C:
\begin{footnotesize}
\begin{align}
\textbf{C}(t+1) =& \big\langle \boldsymbol{\delta}\textbf{N}(t+1) \boldsymbol{\delta}\textbf{N}(t+1)^T \big\rangle \
=& \textbf{A}\big\langle \boldsymbol{\delta}\textbf{N}(t) \boldsymbol{\delta}\textbf{N}(t)^T \big\rangle \textbf{A}^T + \textbf{A}\big\langle \boldsymbol{\delta}\textbf{N}(t) \boldsymbol{\delta}\textbf{N}(t)^T \big\rangle + \big\langle \textbf{z}(t) \boldsymbol{\delta}\textbf{N}(t)^T \big\rangle \textbf{A}^T + \big\langle \textbf{z}(t) \textbf{z}(t)^T \big\rangle \
=& \textbf{A} \textbf{C}(t) \textbf{A}^T + \textbf{Z}_0
\end{align}
\end{footnotesize}
\noindent Taking the limit $t \rightarrow \infty$ on both sides, we get:
\begin{equation}
\textbf{C}^\infty = \textbf{A} \textbf{C}^\infty \textbf{A}^T + \textbf{B}
\end{equation}
\noindent where A is as in Equation S3 and B is:
\begin{equation}
B_{ij} = N_{i}^ N_{j}^ \sigma_{ei} \sigma_{ej} \text{cov}(u_{ei},u_{ej}) + \sqrt{N_{i}^ N_{j}^} \sigma_{di}\sigma_{dj} \text{cov}(u_{di},u_{dj})
\end{equation}
\noindent Similarly, R is the variance-covariance matrix for population growth rates at steady state
\begin{equation}
R_{ij} = \frac{r_{mi}r_{mj}}{K_i K_j} \sum_{k,l} \alpha_{ik} \alpha_{jl} C_{kl} + \sigma_{ei} \sigma_{ej} \text{cov}(u_{ei},u_{ej}) + \frac{\sigma_{di}\sigma_{dj} \text{cov}(u_{di},u_{dj})}{\sqrt{N_{i}^ N_{j}^}}.
\end{equation}
\noindent Then, following the synchrony metric of @Loreau2008a, the sychrony of population biomasses is
\begin{equation}
\phi_N = \frac{\sum_{i,j}C_{ij}}{\left( \sum_i \sqrt{C_{ii}} \right)^2}
\end{equation}
\noindent and the synchrony of per capita growth rates is
\begin{equation}
\phi_R = \frac{\sum_{i,j}R_{ij}}{\left( \sum_i \sqrt{R_{ii}} \right)^2}.
\end{equation}
Synchrony of population biomasses and growth rates emerge from complex interactions among species' intrinsic growth rates, interspecific interactions, environmental stochasticity, and demographic stochasticity.
Given these complexities, it is impossible to determine expected effects of different parameters in a multi-species case.
Thus, we analyze a simplified case where interspecific interactions are zero.
The variance-covariance matrix of population biomasses at steady state then becomes
\begin{equation}
c_{ij} = \frac{B_{ij}}{1-A_ii A_jj} = \frac{N_{i}^ N_{j}^ \sigma_{ei} \sigma_{ej} \text{cov}(u_{ei},u_{ej}) + \sqrt{N_{i}^ N_{j}^} \sigma_{di}\sigma_{dj} \text{cov}(u_{di},u_{dj})}{1 - (1-r_{mi})(1-r_{mj})}
\end{equation}
and the variance-covariance matrix of population growth rates simplifies to
\begin{equation}
R_{ij} = \frac{1}{1- \frac{r_{mi}r_{mj}}{r_{mi}+r_{mj}}} \left( \sigma_{ei} \sigma_{ej} \text{cov}(u_{ei},u_{ej}) + \frac{\sigma_{di}\sigma_{dj} \text{cov}(u_{di},u_{dj})}{\sqrt{N_{i}^ N_{j}^}} \right)
\end{equation}
\noindent Note that synchrony in population sizes and synchrony in growth rates use weighted factors of the environmental variances and covariances with species-specific parameters.
So, we further assume that, along with interspecific interactions being zero, all species have identical growth rates, environmental stochasticity is absent, and all species have identical demographic variance.
This represents a theoretical limiting case where the community consists of identical species coexisting in a constant environment where only demographic stochasticity causes temporal fluctuations.
Under such conditions, synchrony of population biomasses is
\begin{equation}
\phi_N = \frac{1}{\left(\sum_i p_i^{1/2} \right)^2}
\end{equation}
\noindent and synchrony of growth rates is
\begin{equation}
\phi_R = \frac{\sum_i p_i^{-1}}{\left(\sum_i p_i^{-1/2} \right)^2}
\end{equation}
\noindent where $p_i$ is the average frequency of species i, $p_i = N_i/N_T$.
When all species have identical abundances and $p_i = 1/S$, where S is species richness, the both synchrony values equal 1/S [@Loreau2008a]. The prediction represented in Equation S19 is Equation 3 in the main text, which we refer to as $\phi_{R,\mathcal{M}_D}$.
Another limiting case is where only environmental stochasticity is operating: no interspecific interactions, no demographic stochasticity, identical intrinsic growth rates, and environmental stochasticity with same standard deviation for all species.
Under such constraints, the synchrony of of species' biomasses is
\begin{equation}
\phi_N = \sum_{i,j} p_i p_j \text{cov}(u_{ei}, u_{ej})
\end{equation}
\noindent{} and the synchrony of species' growth rates is
\begin{equation}
\phi_R = \frac{\sum_{i,j} \text{cov}(u_{ei}, u_{ej})}{S^2}
\end{equation}
\noindent The prediction represented in Equation S21 is Equation 4 in the main text, which we refer to as $\phi_{R,\mathcal{M}_E}$. If we know the covariance matrix of the species environmental responses, we can calculate the above expectations directly. Note that the environmental responses are normalized in the above equations, therefore the covariances are correlations in the above equations.
\newpage{}
Materials and methods details
Vital rate statistical models
We modeled survival probability and growth on individual genets as a function of genet size, the crowding experienced by the focal genet from both heterospecific and conspecific genets in its neighborhood (described below), temporal varation among years, and spatial variation among quadrat groups. Groups are sets of quadrats located in close proximity within a pasture or grazing exclosure.
We follow the approach of @Chu2015 to estimate crowding, assuming that the crowding experienced by a focal genet depends on distance to each neighbor genet and the neighbor's size, u:
\begin{equation}
w_{ijm,t} = \sum_k e^{-\delta_{jm}d_{ijkm,t}^{2}}u_{km,t}.
\end{equation}
In the above, $w_{ijm,t}$ is the crowding that genet i of species j in year t experiences from neighbors of species m. The spatial scale over which species m neighbors exert influence on any genet of species j is determined by $\delta_{jm}$. The function is applied for all k genets of species m that neighbor the focal genet at time t, and $d_{ijkm,t}$ is the distance between genet i in species j and genet k in species m. When $k=m$, the effect is intraspecific crowding. We use regression-specific (survival and growth) $\delta$ values estimated by Chu and Adler (2015).
We used logistic regression to model survival probability ($S$) of genet $i$ from species $j$ in quadrat group $g$ from time $t$ to $t+1$:
\begin{align}
\text{logit}(S_{ijg,t}) &= \gamma^{S}{j,t} + \phi^{S}{jg} + \beta^{S}{j,t}x{ij,t} + \boldsymbol{\omega}^{S}{j} \textbf{w}{ij,t}
\end{align}
where $x_{ij,t}$ is the log of genet size, $\gamma^{S}{j,t}$ is a year-specific intercept, $\beta^{S}{j,t}$ is the year-specific slope parameter for size, $\phi^{S}_{jg}$ is the random effect of quadrat group location, and $\omega$ is a vector of per capita interaction coefficients which determine the impact of crowding, \textbf{w}, by each species on the focal species.
We modeled genet growth, conditional on survival, in a similar manner:
\begin{align}
\mu_{ijg,t+1} &= \gamma^{G}{j,t} + \phi^{G}{jg} + \beta^{G}{j,t}\mu{ij,t} + \boldsymbol{\omega}^{G}{j} \textbf{w}{ij,t}
\end{align}
where $\mu$ is log genet size and all other parameters are as described for the survival regression. We capture non-constant error variance in growth by modeling the variance around the growth regression ($\varepsilon$) as a nonlinear function of predicted genet size:
\begin{align}
\varepsilon_{ij,t} = a e^{b \mu_{ijg,t+1}}
\end{align}
Our data allows us to track new recruits, but we cannot assign a specific parent to new genets. THerefore, we model recruitment at the quadrat level: the number of new individuals of species $j$ in quadrat $q$ recruiting at time $t+1$ as a function of quadrat "effective cover" ($A'$) in the previous year ($t$). Effective cover is a mixture of observed cover ($A$) in the focal quadrat ($q$) and the mean cover across the entire group ($\bar{A}$) of $Q$ quadrats in which $q$ is located:
\begin{equation}
A'{jq,t} = p{j}A_{jq,t} + (1-p_{j})\bar{A}_{jQ,t}
\end{equation}
where $p$ is a mixing fraction between 0 and 1 that is estimated within the model.
We assume the number of individuals, $y^{R}$, recruiting at time $t+1$ follows a negative binomial distribution:
\begin{equation}
y^{R}{jq,t+1} \sim \text{NegBin}(\lambda{jq,t+1},\zeta)
\end{equation}
where $\lambda$ is the mean intensity and $\zeta$ is the size parameter. We define $\lambda$ as:
\begin{equation}
\lambda_{jq,t+1} = A'{jq,t}e^{(\gamma^{R}{j,t} + \phi^{R}{jQ} + \theta^{R}{jk}C_{k,t} + \omega^{R}\sqrt{A'_{q,t}})}
\end{equation}
where $A'$ is effective cover ($\text{cm}^2$) of species $j$ in quadrat $q$ and all other terms are as in the survival and growth regressions.
Multi-species populations models
The individually based model (IBM) is straighforward and described in the main text, so here focus on the structure of the integral projeciton model (IPM).
We built an environmentally stochastic IPM. Our IPM follows the specification of @Chu2015 where the population of species j is a density function $n(u_{j},t)$ giving the density of sized-u genets at time t. Genet size is on the natural log scale, so that $n(u_{j},t)du$ is the number of genets whose area (on the arithmetic scale) is between $e^{u_{j}}$ and $e^{u_{j}+du}$. So, the density function for any size v at time $t+1$ is
\begin{equation}
n(v_{j},t+1) = \int_{L_{j}}^{U_{j}} k_{j}(v_{j},u_{j},\bar{\bold{w_{j}}}(u_{j}))n(u_{j},t)
\end{equation}
where $k_{j}(v_{j},u_{j},\bar{\bold{w_{j}}})$ is the population kernal that describes all possible transitions from size $u$ to $v$ and $\bar{\bold{w_{j}}}$ is a vector of estimates of average crowding experienced from all other species by a genet of size $u_j$ and species $j$. The integral is evaluated over all possible sizes between predefined lower (L) and upper (U) size limits that extend beyond the range of observed genet sizes.
The population kernal is defined as the joint contributions of survival (S), growth (G), and recruitment (R):
\begin{equation}
k_{j}(v_{j},u_{j},\bar{\bold{w_{j}}}) = S_j(u_j, \bar{\bold{w_{j}}}(u_{j}))G_j(v_{j},u_{j},\bar{\bold{w_{j}}}(u_{j})) + R_j(v_{j},u_{j},\bar{\bold{w_{j}}}),
\end{equation}
\noindent{} which means we are calculating growth (G) for individuals that survive (S) from time t to t+1 and adding in newly recruited (R) individuals of an average sized one-year-old genet for the focal species. Our stastical model for recruitment (R, described below) returns the number of new recruit produced per quadrat. Following previous work, we assume that fecundity increases linearly with size ($R_j(v_{j},u_{j},\bar{\bold{w_{j}}}) = e^{u_j}R_j(v_{j},\bar{\bold{w_{j}}})$) to incorporate the recruitment function in the spatially-implicit IPM.
\newpage{}
Results for synchrony of percent cover
Synchrony of percent cover from population model simulations did not compare well with observed synchrony of cover or analytical predictions for synchrony cover (Fig. S2, S3).
We did not expect synchrony of percent cover from simulations to compare well with observed synchrony of cover because our models represent equilibrium dynamics rather than a specific set of observed years.
Also, unlike growth rates, synchrony in cover is highly influenced by drift and legacy effects [@Loreau2008a].
More interesting is that our simulation results did not compare well with analytical predictions.
One issue is that some species are not well-regulated around their equilibrium cover, so that our linear approximation required for anayltical predictions is likely to fail.
Analytical predictions for synchrony of growth rates are more informative than synchrony in population sizes because the linear approximation is more proximate.
Another issue is the influence of demographic stochasticity in some communities, for example, in Idaho where our simulations with only demographic stochasticity yielded results far different from analytical predicitons (Fig. S2).
Our analytical predictions require the unrealistic assumption that demographic variances among species are equal.
However, in Idaho, A. tripartita has much higher demographic stochasticity than the other species, so that variation in its abundance dominates.
Thus, the assumption that all species have similar stochasticity fails.
In combination, our results from analyzing synchrony of percent cover indicate that the best way to decipher the mechanisms contributing to community synchrony is to analyze growth rates, in agreement with the theoretical arguments of @Loreau2008a.
\newpage{}
Supporting Figures
\begin{figure}[!ht]
\centering
\includegraphics[width=6in]{./components/figureS1.png}
\caption{Observed (vertical dashed lines) and simulated (solid density curves) synchrony of species per capita growth rates at each site from the IBM (top panels) and the IPM (bottom panels). IPM density curves come from 100 random contiguous sections from 2,000 iteration IPM runs, where the length of each randomly selected section is equal to the number of observation years for each data set. IBM density curves come from 100 replicate simulations of 75 iterations each. Synchrony values come from simulations where environmental stochasticity and interspecific interactions are present. The IBM was run on a 5 by 5 meter landscape to reduce the effect of demographic stochasticity.}
\end{figure}
\pagebreak{}
\begin{figure}[!ht]
\centering
\includegraphics[width=6in]{./components/formatted_figures/formatted_figureS2.png}
\caption{Community-wide synchrony of species percent cover from model simulation experiments. Synchrony of species' percent cover for each study area are from simulation experiments with demographic stochasticity, environmental stochasticity, and interspecific competition present (All Drivers''), demographic stochasticity removed (No D.S.''), environmental stochasticity removed (No E.S.''), interspecific competition removed (No Comp.''), interspecific competition and demographic stochasticity removed (No Comp. + No D.S.''), and interspecific competition and environmental stochasticity removed (No Comp. + No E.S.''). Abbreviations within the bars for the New Mexico site indicate whether the IBM or IPM was used for a particular simulation. Error bars represent the 2.5\% and 97.5\% quantiles from model simulations.}
\end{figure}
\pagebreak{}
\begin{figure}[!ht]
\centering
\includegraphics[width=6in]{./components/formatted_figures/formatted_figureS3.png}
\caption{Synchrony of species' percent cover for each study area from IBM simulations across different landscape sizes when only demographic stochastcity is present (D.S. Only'') and when environmental stochasticity is also present removed (D.S. + E.S.''). The horizontal lines show the analytical predictions $\mathcal{M}_D$ (dashed line) and $\mathcal{M}_E$ (dotted line). The strength of demographic stochasticity decreases as landscape size increases because population sizes also increase. Error bars represent the 2.5\% and 97.5\% quantiles from model simulations.}
\end{figure}
\pagebreak{}
\begin{figure}[!ht]
\centering
\includegraphics[width=6in]{./components/formatted_figures/environmental_variances_cover.png}
\caption{Variance of percent cover for each species (along x-axes) in each site through time from simulations with only environmental stochasticity operating (IPM with no species interactions).}
\end{figure}
\begin{figure}[!ht]
\centering
\includegraphics[width=6in]{./components/formatted_figures/environmental_variances_pgr.png}
\caption{Variance of per capita growth rates for each species (along x-axes) in each site through time from simulations with only environmental stochasticity operating (IPM with no species interactions).}
\end{figure}
\newpage{}
\begin{figure}[!ht]
\centering
\includegraphics[width=6in]{./components/formatted_figures/demographic_variances_cover.png}
\caption{Variance of percent cover for each species (along x-axes) in each site through time from simulations with only demographic stochasticity operating (IBM simulated on 5 by 5 meter landscape).}
\end{figure}
\begin{figure}[!ht]
\centering
\includegraphics[width=6in]{./components/formatted_figures/demographic_variances_pgr.png}
\caption{Variance of per capita growth rates for each species (along x-axes) in each site through time from simulations with only demographic stochasticity operating (IBM simulated on 5 by 5 meter landscape).}
\end{figure}
\pagebreak{}
\begin{figure}[!ht]
\centering
\includegraphics[width=3in]{./components/formatted_figures/formatted_figureS4.png}
\caption{Community-wide synchrony of species' growth rates from model simulation experiments for the Idaho community with \textit{Artemisia tripartita} removed. Synchrony of species' growth rates are from simulation experiments with demographic stochasticity, environmental stochasticity, and interspecific competition present (All Drivers''), demographic stochasticity removed (No D.S.''), environmental stochasticity removed (No E.S.''), interspecific competition removed (No Comp.''), interspecific competition and demographic stochasticity removed (No Comp. + No D.S.''), and interspecific competition and environmental stochasticity removed (No Comp. + No E.S.''). Error bars represent the 2.5\% and 97.5\% quantiles from model simulations.}
\end{figure}
\newpage{}
Supporting Tables
In all tables, species codes are as follows (also see Table 1 in the main text):
- Bouteloua eriopoda = BOER
- Sporobolus flexuosus = SPFL
- Bouteloua rothrockii = BORO
- Bouteloua curtipendula = BOCU
- Bouteloua hirsuta = BOHI
- Schizachyrium scoparium = SCSC
- Bouteloua gracilis = BOGR
- Hesperostipa comata = HECO
- Pascopyrum smithii = PASM
- Poa secunda = POSE
- Artemisia tripartita = ARTR
- Pseudoroegneria spicata = PSSP
library(communitySynchrony)
library(plyr)
library(reshape2)
library(synchrony)
library(ggplot2)
library(xtable)
library(gridExtra)
library(RColorBrewer)
site <- "Kansas"
spp_list <- c("BOCU","BOHI","SCSC")
num_spp <- length(spp_list)
ks_data <- data.frame(quad=NA, year=NA, totCover=NA, species=NA)
for(dospp in 1:num_spp){ #loop through species to read in data
spp_now <- spp_list[dospp]
quad_file <- paste("../data/", site,"/",spp_now,"/quadratCover.csv",sep="")
spp_data <- read.csv(quad_file)
spp_data$species <- spp_now
ks_data <- rbind(ks_data, spp_data)
} #end species looping for raw data
ks_data <- ks_data[2:nrow(ks_data),] #remove first NA row
tmp1<-which(ks_data$quad_data=="q25" & (ks_data$year<35 | ks_data$year>62))
tmp2<-which(ks_data$quad_data=="q27")
tmp3<-which(ks_data$quad=="q28")
tmp4<-which(ks_data$quad=="q30")
tmp5<-which(ks_data$quad=="q31" & (ks_data$year<35 | ks_data$year>39))
tmp6<-which(ks_data$quad=="q32" & (ks_data$year<35 | ks_data$year>41))
tmp<-c(tmp1,tmp2,tmp3,tmp4,tmp5,tmp6)
ks_data<-ks_data[-tmp,]
# exclude the records later than 1968, to keep the same random year effect...
ks_data<-subset(ks_data,year<68)
site <- "Idaho"
spp_list <- sort(c("PSSP","HECO","POSE","ARTR"))
num_spp <- length(spp_list)
id_data <- data.frame(quad=NA, year=NA, totCover=NA, species=NA)
for(dospp in 1:num_spp){ #loop through species to read in data
spp_now <- spp_list[dospp]
quad_file <- paste("../data/", site,"/",spp_now,"/quadratCover.csv",sep="")
spp_data <- read.csv(quad_file)
spp_data$species <- spp_now
id_data <- rbind(id_data, spp_data)
} #end species looping for raw data
id_data <- id_data[2:nrow(id_data),] #remove first NA row
# id_data <- subset(id_data, species!="ARTR") #take out the shrub
site <- "Montana"
spp_list <- sort(c("BOGR","HECO","PASM","POSE"))
num_spp <- length(spp_list)
mt_data <- data.frame(quad=NA, year=NA, totCover=NA, species=NA)
for(dospp in 1:num_spp){ #loop through species to read in data
spp_now <- spp_list[dospp]
quad_file <- paste("../data/", site,"/",spp_now,"/quadratCover.csv",sep="")
spp_data <- read.csv(quad_file)
spp_data$species <- spp_now
mt_data <- rbind(mt_data, spp_data)
} #end species looping for raw data
mt_data <- mt_data[2:nrow(mt_data),] #remove first NA row
site <- "NewMexico"
spp_list <- sort(c("BOER","SPFL"))
num_spp <- length(spp_list)
nm_data <- data.frame(quad=NA, year=NA, totCover=NA, species=NA)
for(dospp in 1:num_spp){ #loop through species to read in data
spp_now <- spp_list[dospp]
quad_file <- paste("../data/", site,"/",spp_now,"/quadratCover.csv",sep="")
spp_data <- read.csv(quad_file)
spp_data$species <- spp_now
nm_data <- rbind(nm_data, spp_data)
} #end species looping for raw data
nm_data <- nm_data[2:nrow(nm_data),] #remove first NA row
site <- "Arizona"
spp_list <- sort(c("BOER","BORO"))
num_spp <- length(spp_list)
az_data <- data.frame(quad=NA, year=NA, totCover=NA, species=NA)
for(dospp in 1:num_spp){ #loop through species to read in data
spp_now <- spp_list[dospp]
quad_file <- paste("../data/", site,"/",spp_now,"/quadratCover.csv",sep="")
spp_data <- read.csv(quad_file)
spp_data$species <- spp_now
az_data <- rbind(az_data, spp_data)
} #end species looping for raw data
az_data <- az_data[2:nrow(az_data),] #remove first NA row
out_ks <- get_comm_synchrony(ts_data = ks_data)
out_id <- get_comm_synchrony(ts_data = id_data)
out_mt <- get_comm_synchrony(ts_data = mt_data)
out_nm <- get_comm_synchrony(ts_data = nm_data)
out_az <- get_comm_synchrony(ts_data = az_data)
get_table_metrics <- function(output){
tmp <- output
exp_tmp <- round(tmp$pgr_expected_synch_ind_flucts, 2)
obs_tmp <- round(as.numeric(tmp$pgr_synchrony[1]), 2)
demo_tmp <- obs_tmp - exp_tmp
exp_cover <- round(as.numeric(tmp$cover_expected_synch_ind_flucts), 2)
obs_cover <- round(as.numeric(tmp$abund_synchrony[1]), 2)
demo_cover <- obs_cover - exp_cover
plant_size <- tmp$avg_size
cover_var <- tmp$variability
return(c(obs_tmp, exp_tmp, demo_tmp, obs_cover, exp_cover,
demo_cover, cover_var))
}
table_data <- as.data.frame(rbind(get_table_metrics(out_az),
get_table_metrics(out_id),
get_table_metrics(out_ks),
get_table_metrics(out_mt),
get_table_metrics(out_nm)))
table_data$Site <- c("Arizona", "Idaho", "Kansas", "Montana", "New Mexico")
alltable_data <- table_data
exp_synchrony_env <- readRDS("../results/envstoch_predictions.RDS")
exp_synchrony_env <- exp_synchrony_env[complete.cases(exp_synchrony_env),]
exp_synchrony_demo <- alltable_data[,c("V2","V5")]
exp_synchrony_demo$site <- table_data$Site
colnames(exp_synchrony_demo) <- c("pgr", "cover", "site")
sims_output <- readRDS("../results/allsims_plot_data.RDS")
pred_table <- data.frame(site = table_data$Site,
exp_demo = exp_synchrony_demo$pgr,
sim_demo = round(subset(sims_output, simulation=="6No Comp. + No E.S.")["synchrony"],2),
exp_env = round(exp_synchrony_env$growthrate_prediction,2),
sim_env = round(subset(sims_output, simulation=="2No Comp. + No D.S.")["synchrony"],2))
colnames(pred_table) <- c("Site",
"Predicted $\\phi_{R,\\mathcal{M}_D}$",
"Simulated $\\phi_{R,\\mathcal{M}_D}$",
"Predicted $\\phi_{R,\\mathcal{M}_E}$",
"Simulated $\\phi_{R,\\mathcal{M}_E}$")
myorder <- c(5,1,3,4,2)
xpred_table <- pred_table[myorder,]
pred_cap <- "Comparisons between our analytical predictions and simulation results for synchrony of species' per capita growth rates. Analytical predictions represent two limiting cases where only demographic stochasticity is operating ($\\phi_{R,\\mathcal{M}_D}$) and where only environmental stochasticity is operating ($\\phi_{R,\\mathcal{M}_E}$). Simulated synchrony values come from our empirically-based, multi-species population models when simulated under conditions that match the limiting case conditions (e.g., environmental stochasticity and competition removed for $\\mathcal{M}_D$)."
print(xtable(xpred_table, caption = pred_cap),
caption.placement="top",
include.rownames = F,
sanitize.colnames.function = identity,
comment=FALSE,
size="normalsize")
library(reshape2)
library(plyr)
library(ggplot2)
library(xtable)
perc_changes <- read.csv("../results/synchsims_percent_diffs.csv")
perc_changes <- perc_changes[which(perc_changes$simulation!="3All Drivers"), c("site", "simulation", "percent_diff")]
table_diff <- dcast(perc_changes, simulation~site, value.var = "percent_diff")
diff_cap <- "Percent differences of synchrony of per capita growth rates between each removal simulation experiment and the 'All Drivers' simulation."
print(xtable(table_diff, caption = diff_cap),
caption.placement="top",
include.rownames = FALSE,
sanitize.colnames.function = identity,
comment=FALSE)
library(reshape2)
library(plyr)
library(ggplot2)
library(xtable)
####
#### Read in regression fits --------------------------------------------------
####
all_files <- list.files("../results/")
param_files1 <- all_files[grep("param*", all_files)]
param_files <- param_files1[grep("RDS", param_files1)]
num_vitals <- length(param_files)
list_yr_effects <- list()
for(j in 1:num_vitals){
param_list <- readRDS(paste0("../results/",param_files[j]))
vital_name <- unlist(strsplit(param_files[j], "_"))[1]
if(vital_name != "recruit"){
site_names <- names(param_list)
site_year_list <- list()
for(do_site in site_names){
tmp_site <- param_list[[do_site]]
num_spp <- length(tmp_site)
tmp_yr_int <- matrix(ncol = num_spp, nrow=nrow(tmp_site[[1]]))
tmp_yr_slope <- matrix(ncol = num_spp, nrow=nrow(tmp_site[[1]]))
for(i in 1:num_spp){
spp_df <- tmp_site[[i]]
tmp_yr_int[,i] <- spp_df$Intercept.yr
col2use <- grep("logarea", colnames(spp_df))
col2use2 <- grep(".yr", colnames(spp_df[col2use]))
tmp_yr_slope[,i] <- spp_df[,col2use[col2use2]]
} # end species within site loop
site_year_list[[do_site]][["intercept"]] <- tmp_yr_int
site_year_list[[do_site]][["slope"]] <- tmp_yr_slope
} # end site within vital rate loop
} # end NOT recruit loops
if(vital_name == "recruit"){
site_names <- names(param_list)
site_year_list <- list()
for(do_site in site_names){
tmp_site <- param_list[[do_site]]
tmp_yrs <- tmp_site[grep("intcpt.yr", rownames(tmp_site)),"Mean"]
stringlist <- strsplit(names(tmp_yrs), ",")
num_spp <- length(unique(sapply(stringlist, "[[", 2)))
mat_yrs <- t(matrix(tmp_yrs, nrow=num_spp, byrow = TRUE))
site_year_list[[do_site]] <- mat_yrs
} # end site within vital rate loop
} # end YES recruit loops
list_yr_effects[[vital_name]] <- site_year_list
} # end vital rate loop
# ## Plot histograms of year effects by vital rate
# vital_rates <- names(list_yr_effects)
# par(mfrow=c(3,5))
# for(do_vital in vital_rates){
# tmp_vital <- list_yr_effects[[do_vital]]
# site_names <- names(tmp_vital)
# for(do_site in site_names){
# tmp_site <- tmp_vital[[do_site]]
# hist(tmp_site, main=paste(do_site,do_vital), xlab="Year Effect", xlim=c(-3,8))
# }
# }
## Get standard deviation of year effect posterior
####
#### Calculate average correlation of year effects ----------------------------
####
vital_rates <- names(list_yr_effects)
out_df <- data.frame(site=NA, vital_rate=NA, term=NA, corr=NA)
for(do_vital in vital_rates){
tmp_vital <- list_yr_effects[[do_vital]]
site_names <- names(tmp_vital)
for(do_site in site_names){
tmp_site <- tmp_vital[[do_site]]
if(do_vital != "recruit"){
for(i in 1:length(tmp_site)){
tmp_cor <- cor(tmp_site[[i]])
avg_cor <- mean(tmp_cor[upper.tri(tmp_cor)])
tmp_out <- data.frame(site=do_site, vital_rate=do_vital, term=names(tmp_site)[i], corr=avg_cor)
out_df <- rbind(out_df, tmp_out)
}
}
if(do_vital == "recruit"){
tmp_cor <- cor(tmp_site)
avg_cor <- mean(tmp_cor[upper.tri(tmp_cor)])
tmp_out <- data.frame(site=do_site, vital_rate=do_vital, term="intercept", corr=avg_cor)
out_df <- rbind(out_df, tmp_out)
}
}
}
cor_df <- out_df[2:nrow(out_df),]
cor_cast <- dcast(cor_df, site+term~vital_rate)
cor.cap <- "Correlations of species' year random effects for each site by term, where term refers to the random effect on the slope or the intercept."
print(xtable(cor_cast, caption = cor.cap),
caption.placement="top",
include.rownames = FALSE,
sanitize.colnames.function = identity,
comment=FALSE)
library(reshape2)
library(plyr)
library(ggplot2)
library(xtable)
####
#### Get Vital Rate Satistical Results Files ----------------------------------
####
vital_rates <- c("growth", "surv", "recruit")
result_files <- list.files("../results/")
vr_files1 <-result_files[grep(paste(vital_rates,collapse="|"),
list.files("../results/"))]
vr_files <- vr_files1[grep("RDS", vr_files1)]
####
#### Loop Over Files and Extract Crowding Effects -----------------------------
####
# Growth and survival first
vital_mat_list <- list()
for(do_vital in vital_rates[c(1,2)]){
vr_do <- vr_files[grep(do_vital,vr_files)]
torms <- grep("yr", vr_do)
vr_do <- vr_do[-torms]
tmp_vital <- readRDS(paste0("../results/",vr_do))
sites <- names(tmp_vital)
mat_list <- list()
for(do_site in sites){
tmp <- tmp_vital[[do_site]]
spp <- names(tmp)
tmpmat <- matrix(0,nrow = length(spp), ncol = length(spp))
for(i in 1:length(spp)){
do_spp <- spp[i]
tmpspp <- tmp[[do_spp]]
tmpw <- tmpspp[1,grep("crowd*", colnames(tmpspp))]
tmpmat[i,] <- as.numeric(tmpw)
}
rownames(tmpmat) <- spp
colnames(tmpmat) <- spp
mat_list[[do_site]] <- tmpmat
}
vital_mat_list[[do_vital]] <- mat_list
}
# Recruitment
vr_do <- vr_files[grep("recruit",vr_files)]
torms <- c(grep("yr", vr_do), grep("MONT", vr_do))
vr_do <- vr_do[-torms]
tmp_vital <- readRDS(paste0("../results/",vr_do))
names_vital <- readRDS(paste0("../results/",vr_files[1]))
sites <- names(tmp_vital)
mat_list <- list()
for(do_site in sites){
tmp <- tmp_vital[[do_site]]
spp <- names(names_vital[[do_site]])
tmpdd <- tmp[grep("dd", rownames(tmp)),"Mean"]
tmpmat <- matrix(tmpdd,length(spp),length(spp), byrow = FALSE)
rownames(tmpmat) <- spp
colnames(tmpmat) <- spp
mat_list[[do_site]] <- tmpmat
}
vital_mat_list[["recruit"]] <- mat_list
####
#### Calculate Average Inter and Intraspecific Interactions -------------------
####
vital_comp_list <- list()
for(do_vital in vital_rates){
tmp <- vital_mat_list[[do_vital]]
sites <- names(tmp)
comp_list <- list()
for(do_site in sites){
tmpsite <- tmp[[do_site]]
avgintra <- mean(diag(tmpsite))
avginter <- mean(c(tmpsite[upper.tri(tmpsite)], tmpsite[lower.tri(tmpsite)]))
avgcomp <- matrix(c(avginter, avgintra),1,2)
colnames(avgcomp) <- c("inter", "intra")
comp_list[[do_site]] <- avgcomp
}
vital_comp_list[[do_vital]] <- comp_list
}
df <- melt(unlist(vital_comp_list))
pieces <- strsplit(rownames(df), split = "[.]")
df$vitalrate <- sapply(pieces, "[", 1)
sites <- sapply(pieces, "[", 2)
df$site <- substr(sites, start = 1, stop = nchar(sites)-1)
df$interintra <- substr(sites, start = nchar(sites), stop = nchar(sites))
dftable <- dcast(df, site+interintra~vitalrate, value.var = "value")
colnames(dftable) <- c("Site", "Interaction Type", "Growth", "Recruitment", "Survival")
dftable[which(dftable$`Interaction Type`==1),"Interaction Type"] <- "Interspecific"
dftable[which(dftable$`Interaction Type`==2),"Interaction Type"] <- "Intraspecific"
# avg_interintra <- ddply(df, .(site, interintra), summarise,
# avg_val = mean(value))
# avg_interintra[which(avg_interintra$interintra==1),"interintra"] <- "Interspecific Competition"
# avg_interintra[which(avg_interintra$interintra==2),"interintra"] <- "Intraspecific Competition"
# interintra_table <- dcast(avg_interintra, site~interintra, value.var = "avg_val")
comp_cap <- "Average interaction coefficients for each vital rate for each community."
print(xtable(dftable, caption = comp_cap, digits=4),
caption.placement="top",
include.rownames = FALSE,
sanitize.colnames.function = identity,
comment=FALSE)
####
#### Get Vital Rate Satistical Results Files ----------------------------------
####
vital_rates <- c("growth", "surv", "recruit")
result_files <- list.files("../results/")
vr_files <-result_files[grep(paste(vital_rates,collapse="|"),
list.files("../results/"))]
####
#### Loop Over Files and Extract Crowding Effects -----------------------------
####
# Growth and survival first
vital_mat_list <- list()
for(do_vital in vital_rates[c(1,2)]){
vr_do <- vr_files[grep(do_vital,vr_files)]
torms <- grep("yr", vr_do)
vr_do <- vr_do[-torms]
tmp_vital <- readRDS(paste0("../results/",vr_do))
sites <- names(tmp_vital)
mat_list <- list()
for(do_site in sites){
tmp <- tmp_vital[[do_site]]
spp <- names(tmp)
tmpmat <- matrix(0,nrow = length(spp), ncol = length(spp))
for(i in 1:length(spp)){
do_spp <- spp[i]
tmpspp <- tmp[[do_spp]]
tmpw <- tmpspp[1,grep("crowd*", colnames(tmpspp))]
tmpmat[i,] <- as.numeric(tmpw)
}
rownames(tmpmat) <- spp
colnames(tmpmat) <- spp
mat_list[[do_site]] <- tmpmat
}
vital_mat_list[[do_vital]] <- mat_list
}
# Recruitment
vr_do <- vr_files[grep("recruit",vr_files)]
vr_do <- vr_do[which(vr_do=="recruit_parameters.RDS")]
tmp_vital <- readRDS(paste0("../results/",vr_do))
names_vital <- readRDS(paste0("../results/",vr_files[1]))
sites <- names(tmp_vital)
mat_list <- list()
for(do_site in sites){
tmp <- tmp_vital[[do_site]]
spp <- names(names_vital[[do_site]])
tmpdd <- tmp[grep("dd", rownames(tmp)),"Mean"]
tmpmat <- matrix(tmpdd,length(spp),length(spp), byrow = FALSE)
rownames(tmpmat) <- spp
colnames(tmpmat) <- spp
mat_list[[do_site]] <- tmpmat
}
vital_mat_list[["recruit"]] <- mat_list
names(vital_mat_list) <- c("growth", "survival", "recruitment")
for(do_vital in names(vital_mat_list)){
tmp_list <- vital_mat_list[[do_vital]]
for(do_site in names(tmp_list)){
site_list <- tmp_list[[do_site]]
synch_cap <- paste0("Interaction coefficients for ", do_vital, " regressions in ", do_site, ".")
# digitvec <- rep(3,ncol(site_list))
print(xtable(site_list, caption = synch_cap, digits = 4),
caption.placement="top",
include.rownames = TRUE,
sanitize.colnames.function = identity,
comment=FALSE)
}
}
\newpage{}
References
atredennick/community_synchrony documentation built on May 10, 2019, 2:10 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
\renewcommand\linenumberfont{\normalfont\tiny\sffamily\color{gray}}
library("methods") library("knitr") basename <- "ms_communitySynchrony" opts_chunk$set(fig.path = paste("components/figure/", basename, "-", sep=""), cache.path = paste("components/cache/", basename, "/", sep="")) opts_chunk$set(cache = 2) opts_chunk$set(tidy=FALSE, warning=FALSE, message=FALSE, comment = NA, verbose = TRUE, echo=FALSE) # PDF-based figures opts_chunk$set(dev='pdf')
\begin{center} \textbf{\Large{Online Supporting Information: Appendix S1}} \ A.T. Tredennick, C. de Mazancourt, M. Loreau, \& P.B. Adler \ ``Environmental responses, not species interactions, determine synchrony...'' \ \emph{Ecology} \end{center}
\renewcommand{\theequation}{S\arabic{equation}} \renewcommand{\thetable}{S\arabic{table}} \renewcommand{\thefigure}{S\arabic{figure}}
\tableofcontents \listoffigures \listoftables \newpage{}
Derivation of limiting case predictions for synchrony of per capita growth rates (Equations 3-4)
Following @Loreau2013 and @DeMazancourt2013a, we define population growth, ignoring observation error, as
\begin{align} r_i(t) = r_{mi} \left[ 1- \frac{N_i(t)+\sum_{j \neq i} \alpha_{ij}N_j(t)} {K_i} + \sigma_{ei}u_{ei}(t) + \frac{\sigma_{di}u_{di}(t)}{\sqrt{N_i(t)}} \right] \end{align}
\noindent where $N_i(t)$ is the biomass of species i in year t, and $r_i(t)$ is its growth rate in year t. $r_{mi}$ is species i's intrinsic rate of increase, $K_i$ is its carrying capacity, and $\alpha_{ij}$ is the interspecific competition coefficient representing the effect of species j on species i. Environmental stochasticity is incorporated as $\sigma_{ei}u_{ei}(t)$, where $\sigma_{ei}^2$ is the environmental variance and $u_{ei}$ are normal random variables with zero mean and unit variance that are independent through time but may be correlated between species. Demographic stochasticity arises from variations in births and deaths among individuals (e.g., same states, different fates), and is included in the model as a first order, normal approximation [@Lande2003; @DeMazancourt2013a]. $\sigma_{di}^2$ is the demographic variance and $u_{di}(t)$ are independent normal variables with zero mean and unit variance.
First-order approximations of the temporal variance of total community biomass are obtained as follows [@Ives1995; @Hughes2000; @Ives2002; @Loreau2008a; @DeMazancourt2013a]. Let $\delta N_i(t) = N_i(t) - N^_i$ denote the deviation of observed species i’s biomass from its equilibrium value in the community, $N^i$, in the absence of stochasticity. Equation (S1) can be Taylor expanded around $\delta N_i(t) = u{ei}(t) = u_{di}(t) = u^S_{oi}(t) = 0$ to yield, after dropping terms of order two and higher,
\begin{equation} \boldsymbol{\delta}\textbf{N}(t+1) = \textbf{A}\boldsymbol{\delta}\textbf{N}(t) + \textbf{z}(t), \end{equation}
\noindent where $\boldsymbol{\delta}\textbf{N}(t)$ is the vector of deviations of species biomasses from their deterministic equilibrium value, A is the community matrix, also known as the Jacobian matrix around the equilibrium, with elements $(a_{ij})_{1<i,j<S}$:
\begin{align} A_{ij} = \begin{cases} 1-r_{mi} \frac{N_{i}^}{K_i}, & i=j \ -r_{mi} \frac{N_{i}^}{K_i} \alpha_{ij}, & i \neq j \end{cases} \end{align}
\noindent and $\textbf{z}(t)$ is a vector that encapsulates the effects of environmental and demographic stochasticity whose elements are
\begin{equation} z_i(t) = N^i \sigma{ei} u_{ei}(t) + \sqrt{N^i} \sigma{di} u_{di}(t) \end{equation}
When the system reaches a stationary distribution, the variances and covariances between species biomass time series are:
\begin{equation} \big\langle \boldsymbol{\delta}\textbf{N}(t) \boldsymbol{\delta}\textbf{N}(t)^T \big\rangle = \textbf{C}^\infty = (\text{cov}(N_i,N_j)){1<i,j<S} = (\text{cov}(\delta N_i, \delta N_j)){1<i,j<S} \end{equation}
\noindent where $\boldsymbol{\delta}\textbf{N}^T$ is the transpose of vector $\boldsymbol{\delta}\textbf{N}$, i.e. a row vector.
Our assumptions listed above lead to the following correlation structure of z:
\begin{equation} \big\langle \textbf{z}(t) \textbf{z}(t)^T \big\rangle = \textbf{B} \end{equation}
\begin{equation} \big\langle \textbf{z}(t-s) \textbf{z}(t)^T \big\rangle = 0 \qquad \text{for} \: s>0 \end{equation}
We use Equation S2 to write a dynamical equation for the covariance C:
\begin{footnotesize} \begin{align} \textbf{C}(t+1) =& \big\langle \boldsymbol{\delta}\textbf{N}(t+1) \boldsymbol{\delta}\textbf{N}(t+1)^T \big\rangle \ =& \textbf{A}\big\langle \boldsymbol{\delta}\textbf{N}(t) \boldsymbol{\delta}\textbf{N}(t)^T \big\rangle \textbf{A}^T + \textbf{A}\big\langle \boldsymbol{\delta}\textbf{N}(t) \boldsymbol{\delta}\textbf{N}(t)^T \big\rangle + \big\langle \textbf{z}(t) \boldsymbol{\delta}\textbf{N}(t)^T \big\rangle \textbf{A}^T + \big\langle \textbf{z}(t) \textbf{z}(t)^T \big\rangle \ =& \textbf{A} \textbf{C}(t) \textbf{A}^T + \textbf{Z}_0 \end{align} \end{footnotesize}
\noindent Taking the limit $t \rightarrow \infty$ on both sides, we get:
\begin{equation} \textbf{C}^\infty = \textbf{A} \textbf{C}^\infty \textbf{A}^T + \textbf{B} \end{equation}
\noindent where A is as in Equation S3 and B is:
\begin{equation} B_{ij} = N_{i}^ N_{j}^ \sigma_{ei} \sigma_{ej} \text{cov}(u_{ei},u_{ej}) + \sqrt{N_{i}^ N_{j}^} \sigma_{di}\sigma_{dj} \text{cov}(u_{di},u_{dj}) \end{equation}
\noindent Similarly, R is the variance-covariance matrix for population growth rates at steady state
\begin{equation} R_{ij} = \frac{r_{mi}r_{mj}}{K_i K_j} \sum_{k,l} \alpha_{ik} \alpha_{jl} C_{kl} + \sigma_{ei} \sigma_{ej} \text{cov}(u_{ei},u_{ej}) + \frac{\sigma_{di}\sigma_{dj} \text{cov}(u_{di},u_{dj})}{\sqrt{N_{i}^ N_{j}^}}. \end{equation}
\noindent Then, following the synchrony metric of @Loreau2008a, the sychrony of population biomasses is
\begin{equation} \phi_N = \frac{\sum_{i,j}C_{ij}}{\left( \sum_i \sqrt{C_{ii}} \right)^2} \end{equation}
\noindent and the synchrony of per capita growth rates is
\begin{equation} \phi_R = \frac{\sum_{i,j}R_{ij}}{\left( \sum_i \sqrt{R_{ii}} \right)^2}. \end{equation}
Synchrony of population biomasses and growth rates emerge from complex interactions among species' intrinsic growth rates, interspecific interactions, environmental stochasticity, and demographic stochasticity. Given these complexities, it is impossible to determine expected effects of different parameters in a multi-species case. Thus, we analyze a simplified case where interspecific interactions are zero. The variance-covariance matrix of population biomasses at steady state then becomes
\begin{equation} c_{ij} = \frac{B_{ij}}{1-A_ii A_jj} = \frac{N_{i}^ N_{j}^ \sigma_{ei} \sigma_{ej} \text{cov}(u_{ei},u_{ej}) + \sqrt{N_{i}^ N_{j}^} \sigma_{di}\sigma_{dj} \text{cov}(u_{di},u_{dj})}{1 - (1-r_{mi})(1-r_{mj})} \end{equation}
and the variance-covariance matrix of population growth rates simplifies to
\begin{equation} R_{ij} = \frac{1}{1- \frac{r_{mi}r_{mj}}{r_{mi}+r_{mj}}} \left( \sigma_{ei} \sigma_{ej} \text{cov}(u_{ei},u_{ej}) + \frac{\sigma_{di}\sigma_{dj} \text{cov}(u_{di},u_{dj})}{\sqrt{N_{i}^ N_{j}^}} \right) \end{equation}
\noindent Note that synchrony in population sizes and synchrony in growth rates use weighted factors of the environmental variances and covariances with species-specific parameters. So, we further assume that, along with interspecific interactions being zero, all species have identical growth rates, environmental stochasticity is absent, and all species have identical demographic variance. This represents a theoretical limiting case where the community consists of identical species coexisting in a constant environment where only demographic stochasticity causes temporal fluctuations. Under such conditions, synchrony of population biomasses is
\begin{equation} \phi_N = \frac{1}{\left(\sum_i p_i^{1/2} \right)^2} \end{equation}
\noindent and synchrony of growth rates is
\begin{equation} \phi_R = \frac{\sum_i p_i^{-1}}{\left(\sum_i p_i^{-1/2} \right)^2} \end{equation}
\noindent where $p_i$ is the average frequency of species i, $p_i = N_i/N_T$. When all species have identical abundances and $p_i = 1/S$, where S is species richness, the both synchrony values equal 1/S [@Loreau2008a]. The prediction represented in Equation S19 is Equation 3 in the main text, which we refer to as $\phi_{R,\mathcal{M}_D}$.
Another limiting case is where only environmental stochasticity is operating: no interspecific interactions, no demographic stochasticity, identical intrinsic growth rates, and environmental stochasticity with same standard deviation for all species. Under such constraints, the synchrony of of species' biomasses is
\begin{equation} \phi_N = \sum_{i,j} p_i p_j \text{cov}(u_{ei}, u_{ej}) \end{equation}
\noindent{} and the synchrony of species' growth rates is
\begin{equation} \phi_R = \frac{\sum_{i,j} \text{cov}(u_{ei}, u_{ej})}{S^2} \end{equation}
\noindent The prediction represented in Equation S21 is Equation 4 in the main text, which we refer to as $\phi_{R,\mathcal{M}_E}$. If we know the covariance matrix of the species environmental responses, we can calculate the above expectations directly. Note that the environmental responses are normalized in the above equations, therefore the covariances are correlations in the above equations.
\newpage{}
Materials and methods details
Vital rate statistical models
We modeled survival probability and growth on individual genets as a function of genet size, the crowding experienced by the focal genet from both heterospecific and conspecific genets in its neighborhood (described below), temporal varation among years, and spatial variation among quadrat groups. Groups are sets of quadrats located in close proximity within a pasture or grazing exclosure.
We follow the approach of @Chu2015 to estimate crowding, assuming that the crowding experienced by a focal genet depends on distance to each neighbor genet and the neighbor's size, u:
\begin{equation} w_{ijm,t} = \sum_k e^{-\delta_{jm}d_{ijkm,t}^{2}}u_{km,t}. \end{equation}
In the above, $w_{ijm,t}$ is the crowding that genet i of species j in year t experiences from neighbors of species m. The spatial scale over which species m neighbors exert influence on any genet of species j is determined by $\delta_{jm}$. The function is applied for all k genets of species m that neighbor the focal genet at time t, and $d_{ijkm,t}$ is the distance between genet i in species j and genet k in species m. When $k=m$, the effect is intraspecific crowding. We use regression-specific (survival and growth) $\delta$ values estimated by Chu and Adler (2015).
We used logistic regression to model survival probability ($S$) of genet $i$ from species $j$ in quadrat group $g$ from time $t$ to $t+1$:
\begin{align} \text{logit}(S_{ijg,t}) &= \gamma^{S}{j,t} + \phi^{S}{jg} + \beta^{S}{j,t}x{ij,t} + \boldsymbol{\omega}^{S}{j} \textbf{w}{ij,t} \end{align}
where $x_{ij,t}$ is the log of genet size, $\gamma^{S}{j,t}$ is a year-specific intercept, $\beta^{S}{j,t}$ is the year-specific slope parameter for size, $\phi^{S}_{jg}$ is the random effect of quadrat group location, and $\omega$ is a vector of per capita interaction coefficients which determine the impact of crowding, \textbf{w}, by each species on the focal species.
We modeled genet growth, conditional on survival, in a similar manner:
\begin{align} \mu_{ijg,t+1} &= \gamma^{G}{j,t} + \phi^{G}{jg} + \beta^{G}{j,t}\mu{ij,t} + \boldsymbol{\omega}^{G}{j} \textbf{w}{ij,t} \end{align}
where $\mu$ is log genet size and all other parameters are as described for the survival regression. We capture non-constant error variance in growth by modeling the variance around the growth regression ($\varepsilon$) as a nonlinear function of predicted genet size:
\begin{align} \varepsilon_{ij,t} = a e^{b \mu_{ijg,t+1}} \end{align}
Our data allows us to track new recruits, but we cannot assign a specific parent to new genets. THerefore, we model recruitment at the quadrat level: the number of new individuals of species $j$ in quadrat $q$ recruiting at time $t+1$ as a function of quadrat "effective cover" ($A'$) in the previous year ($t$). Effective cover is a mixture of observed cover ($A$) in the focal quadrat ($q$) and the mean cover across the entire group ($\bar{A}$) of $Q$ quadrats in which $q$ is located:
\begin{equation} A'{jq,t} = p{j}A_{jq,t} + (1-p_{j})\bar{A}_{jQ,t} \end{equation}
where $p$ is a mixing fraction between 0 and 1 that is estimated within the model.
We assume the number of individuals, $y^{R}$, recruiting at time $t+1$ follows a negative binomial distribution:
\begin{equation} y^{R}{jq,t+1} \sim \text{NegBin}(\lambda{jq,t+1},\zeta) \end{equation}
where $\lambda$ is the mean intensity and $\zeta$ is the size parameter. We define $\lambda$ as:
\begin{equation} \lambda_{jq,t+1} = A'{jq,t}e^{(\gamma^{R}{j,t} + \phi^{R}{jQ} + \theta^{R}{jk}C_{k,t} + \omega^{R}\sqrt{A'_{q,t}})} \end{equation}
where $A'$ is effective cover ($\text{cm}^2$) of species $j$ in quadrat $q$ and all other terms are as in the survival and growth regressions.
Multi-species populations models
The individually based model (IBM) is straighforward and described in the main text, so here focus on the structure of the integral projeciton model (IPM). We built an environmentally stochastic IPM. Our IPM follows the specification of @Chu2015 where the population of species j is a density function $n(u_{j},t)$ giving the density of sized-u genets at time t. Genet size is on the natural log scale, so that $n(u_{j},t)du$ is the number of genets whose area (on the arithmetic scale) is between $e^{u_{j}}$ and $e^{u_{j}+du}$. So, the density function for any size v at time $t+1$ is
\begin{equation} n(v_{j},t+1) = \int_{L_{j}}^{U_{j}} k_{j}(v_{j},u_{j},\bar{\bold{w_{j}}}(u_{j}))n(u_{j},t) \end{equation}
where $k_{j}(v_{j},u_{j},\bar{\bold{w_{j}}})$ is the population kernal that describes all possible transitions from size $u$ to $v$ and $\bar{\bold{w_{j}}}$ is a vector of estimates of average crowding experienced from all other species by a genet of size $u_j$ and species $j$. The integral is evaluated over all possible sizes between predefined lower (L) and upper (U) size limits that extend beyond the range of observed genet sizes.
The population kernal is defined as the joint contributions of survival (S), growth (G), and recruitment (R):
\begin{equation} k_{j}(v_{j},u_{j},\bar{\bold{w_{j}}}) = S_j(u_j, \bar{\bold{w_{j}}}(u_{j}))G_j(v_{j},u_{j},\bar{\bold{w_{j}}}(u_{j})) + R_j(v_{j},u_{j},\bar{\bold{w_{j}}}), \end{equation}
\noindent{} which means we are calculating growth (G) for individuals that survive (S) from time t to t+1 and adding in newly recruited (R) individuals of an average sized one-year-old genet for the focal species. Our stastical model for recruitment (R, described below) returns the number of new recruit produced per quadrat. Following previous work, we assume that fecundity increases linearly with size ($R_j(v_{j},u_{j},\bar{\bold{w_{j}}}) = e^{u_j}R_j(v_{j},\bar{\bold{w_{j}}})$) to incorporate the recruitment function in the spatially-implicit IPM.
\newpage{}
Results for synchrony of percent cover
Synchrony of percent cover from population model simulations did not compare well with observed synchrony of cover or analytical predictions for synchrony cover (Fig. S2, S3). We did not expect synchrony of percent cover from simulations to compare well with observed synchrony of cover because our models represent equilibrium dynamics rather than a specific set of observed years. Also, unlike growth rates, synchrony in cover is highly influenced by drift and legacy effects [@Loreau2008a].
More interesting is that our simulation results did not compare well with analytical predictions. One issue is that some species are not well-regulated around their equilibrium cover, so that our linear approximation required for anayltical predictions is likely to fail. Analytical predictions for synchrony of growth rates are more informative than synchrony in population sizes because the linear approximation is more proximate.
Another issue is the influence of demographic stochasticity in some communities, for example, in Idaho where our simulations with only demographic stochasticity yielded results far different from analytical predicitons (Fig. S2). Our analytical predictions require the unrealistic assumption that demographic variances among species are equal. However, in Idaho, A. tripartita has much higher demographic stochasticity than the other species, so that variation in its abundance dominates. Thus, the assumption that all species have similar stochasticity fails. In combination, our results from analyzing synchrony of percent cover indicate that the best way to decipher the mechanisms contributing to community synchrony is to analyze growth rates, in agreement with the theoretical arguments of @Loreau2008a.
\newpage{}
Supporting Figures
\begin{figure}[!ht] \centering \includegraphics[width=6in]{./components/figureS1.png} \caption{Observed (vertical dashed lines) and simulated (solid density curves) synchrony of species per capita growth rates at each site from the IBM (top panels) and the IPM (bottom panels). IPM density curves come from 100 random contiguous sections from 2,000 iteration IPM runs, where the length of each randomly selected section is equal to the number of observation years for each data set. IBM density curves come from 100 replicate simulations of 75 iterations each. Synchrony values come from simulations where environmental stochasticity and interspecific interactions are present. The IBM was run on a 5 by 5 meter landscape to reduce the effect of demographic stochasticity.} \end{figure}
\pagebreak{}
\begin{figure}[!ht]
\centering
\includegraphics[width=6in]{./components/formatted_figures/formatted_figureS2.png}
\caption{Community-wide synchrony of species percent cover from model simulation experiments. Synchrony of species' percent cover for each study area are from simulation experiments with demographic stochasticity, environmental stochasticity, and interspecific competition present (All Drivers''), demographic stochasticity removed (No D.S.''), environmental stochasticity removed (No E.S.''), interspecific competition removed (No Comp.''), interspecific competition and demographic stochasticity removed (No Comp. + No D.S.''), and interspecific competition and environmental stochasticity removed (No Comp. + No E.S.''). Abbreviations within the bars for the New Mexico site indicate whether the IBM or IPM was used for a particular simulation. Error bars represent the 2.5\% and 97.5\% quantiles from model simulations.}
\end{figure}
\pagebreak{}
\begin{figure}[!ht]
\centering
\includegraphics[width=6in]{./components/formatted_figures/formatted_figureS3.png}
\caption{Synchrony of species' percent cover for each study area from IBM simulations across different landscape sizes when only demographic stochastcity is present (D.S. Only'') and when environmental stochasticity is also present removed (D.S. + E.S.''). The horizontal lines show the analytical predictions $\mathcal{M}_D$ (dashed line) and $\mathcal{M}_E$ (dotted line). The strength of demographic stochasticity decreases as landscape size increases because population sizes also increase. Error bars represent the 2.5\% and 97.5\% quantiles from model simulations.}
\end{figure}
\pagebreak{}
\begin{figure}[!ht] \centering \includegraphics[width=6in]{./components/formatted_figures/environmental_variances_cover.png} \caption{Variance of percent cover for each species (along x-axes) in each site through time from simulations with only environmental stochasticity operating (IPM with no species interactions).} \end{figure}
\begin{figure}[!ht] \centering \includegraphics[width=6in]{./components/formatted_figures/environmental_variances_pgr.png} \caption{Variance of per capita growth rates for each species (along x-axes) in each site through time from simulations with only environmental stochasticity operating (IPM with no species interactions).} \end{figure}
\newpage{}
\begin{figure}[!ht] \centering \includegraphics[width=6in]{./components/formatted_figures/demographic_variances_cover.png} \caption{Variance of percent cover for each species (along x-axes) in each site through time from simulations with only demographic stochasticity operating (IBM simulated on 5 by 5 meter landscape).} \end{figure}
\begin{figure}[!ht] \centering \includegraphics[width=6in]{./components/formatted_figures/demographic_variances_pgr.png} \caption{Variance of per capita growth rates for each species (along x-axes) in each site through time from simulations with only demographic stochasticity operating (IBM simulated on 5 by 5 meter landscape).} \end{figure}
\pagebreak{}
\begin{figure}[!ht]
\centering
\includegraphics[width=3in]{./components/formatted_figures/formatted_figureS4.png}
\caption{Community-wide synchrony of species' growth rates from model simulation experiments for the Idaho community with \textit{Artemisia tripartita} removed. Synchrony of species' growth rates are from simulation experiments with demographic stochasticity, environmental stochasticity, and interspecific competition present (All Drivers''), demographic stochasticity removed (No D.S.''), environmental stochasticity removed (No E.S.''), interspecific competition removed (No Comp.''), interspecific competition and demographic stochasticity removed (No Comp. + No D.S.''), and interspecific competition and environmental stochasticity removed (No Comp. + No E.S.''). Error bars represent the 2.5\% and 97.5\% quantiles from model simulations.}
\end{figure}
\newpage{}
Supporting Tables
In all tables, species codes are as follows (also see Table 1 in the main text):
- Bouteloua eriopoda = BOER
- Sporobolus flexuosus = SPFL
- Bouteloua rothrockii = BORO
- Bouteloua curtipendula = BOCU
- Bouteloua hirsuta = BOHI
- Schizachyrium scoparium = SCSC
- Bouteloua gracilis = BOGR
- Hesperostipa comata = HECO
- Pascopyrum smithii = PASM
- Poa secunda = POSE
- Artemisia tripartita = ARTR
- Pseudoroegneria spicata = PSSP
library(communitySynchrony) library(plyr) library(reshape2) library(synchrony) library(ggplot2) library(xtable) library(gridExtra) library(RColorBrewer) site <- "Kansas" spp_list <- c("BOCU","BOHI","SCSC") num_spp <- length(spp_list) ks_data <- data.frame(quad=NA, year=NA, totCover=NA, species=NA) for(dospp in 1:num_spp){ #loop through species to read in data spp_now <- spp_list[dospp] quad_file <- paste("../data/", site,"/",spp_now,"/quadratCover.csv",sep="") spp_data <- read.csv(quad_file) spp_data$species <- spp_now ks_data <- rbind(ks_data, spp_data) } #end species looping for raw data ks_data <- ks_data[2:nrow(ks_data),] #remove first NA row tmp1<-which(ks_data$quad_data=="q25" & (ks_data$year<35 | ks_data$year>62)) tmp2<-which(ks_data$quad_data=="q27") tmp3<-which(ks_data$quad=="q28") tmp4<-which(ks_data$quad=="q30") tmp5<-which(ks_data$quad=="q31" & (ks_data$year<35 | ks_data$year>39)) tmp6<-which(ks_data$quad=="q32" & (ks_data$year<35 | ks_data$year>41)) tmp<-c(tmp1,tmp2,tmp3,tmp4,tmp5,tmp6) ks_data<-ks_data[-tmp,] # exclude the records later than 1968, to keep the same random year effect... ks_data<-subset(ks_data,year<68) site <- "Idaho" spp_list <- sort(c("PSSP","HECO","POSE","ARTR")) num_spp <- length(spp_list) id_data <- data.frame(quad=NA, year=NA, totCover=NA, species=NA) for(dospp in 1:num_spp){ #loop through species to read in data spp_now <- spp_list[dospp] quad_file <- paste("../data/", site,"/",spp_now,"/quadratCover.csv",sep="") spp_data <- read.csv(quad_file) spp_data$species <- spp_now id_data <- rbind(id_data, spp_data) } #end species looping for raw data id_data <- id_data[2:nrow(id_data),] #remove first NA row # id_data <- subset(id_data, species!="ARTR") #take out the shrub site <- "Montana" spp_list <- sort(c("BOGR","HECO","PASM","POSE")) num_spp <- length(spp_list) mt_data <- data.frame(quad=NA, year=NA, totCover=NA, species=NA) for(dospp in 1:num_spp){ #loop through species to read in data spp_now <- spp_list[dospp] quad_file <- paste("../data/", site,"/",spp_now,"/quadratCover.csv",sep="") spp_data <- read.csv(quad_file) spp_data$species <- spp_now mt_data <- rbind(mt_data, spp_data) } #end species looping for raw data mt_data <- mt_data[2:nrow(mt_data),] #remove first NA row site <- "NewMexico" spp_list <- sort(c("BOER","SPFL")) num_spp <- length(spp_list) nm_data <- data.frame(quad=NA, year=NA, totCover=NA, species=NA) for(dospp in 1:num_spp){ #loop through species to read in data spp_now <- spp_list[dospp] quad_file <- paste("../data/", site,"/",spp_now,"/quadratCover.csv",sep="") spp_data <- read.csv(quad_file) spp_data$species <- spp_now nm_data <- rbind(nm_data, spp_data) } #end species looping for raw data nm_data <- nm_data[2:nrow(nm_data),] #remove first NA row site <- "Arizona" spp_list <- sort(c("BOER","BORO")) num_spp <- length(spp_list) az_data <- data.frame(quad=NA, year=NA, totCover=NA, species=NA) for(dospp in 1:num_spp){ #loop through species to read in data spp_now <- spp_list[dospp] quad_file <- paste("../data/", site,"/",spp_now,"/quadratCover.csv",sep="") spp_data <- read.csv(quad_file) spp_data$species <- spp_now az_data <- rbind(az_data, spp_data) } #end species looping for raw data az_data <- az_data[2:nrow(az_data),] #remove first NA row out_ks <- get_comm_synchrony(ts_data = ks_data) out_id <- get_comm_synchrony(ts_data = id_data) out_mt <- get_comm_synchrony(ts_data = mt_data) out_nm <- get_comm_synchrony(ts_data = nm_data) out_az <- get_comm_synchrony(ts_data = az_data) get_table_metrics <- function(output){ tmp <- output exp_tmp <- round(tmp$pgr_expected_synch_ind_flucts, 2) obs_tmp <- round(as.numeric(tmp$pgr_synchrony[1]), 2) demo_tmp <- obs_tmp - exp_tmp exp_cover <- round(as.numeric(tmp$cover_expected_synch_ind_flucts), 2) obs_cover <- round(as.numeric(tmp$abund_synchrony[1]), 2) demo_cover <- obs_cover - exp_cover plant_size <- tmp$avg_size cover_var <- tmp$variability return(c(obs_tmp, exp_tmp, demo_tmp, obs_cover, exp_cover, demo_cover, cover_var)) } table_data <- as.data.frame(rbind(get_table_metrics(out_az), get_table_metrics(out_id), get_table_metrics(out_ks), get_table_metrics(out_mt), get_table_metrics(out_nm))) table_data$Site <- c("Arizona", "Idaho", "Kansas", "Montana", "New Mexico") alltable_data <- table_data exp_synchrony_env <- readRDS("../results/envstoch_predictions.RDS") exp_synchrony_env <- exp_synchrony_env[complete.cases(exp_synchrony_env),] exp_synchrony_demo <- alltable_data[,c("V2","V5")] exp_synchrony_demo$site <- table_data$Site colnames(exp_synchrony_demo) <- c("pgr", "cover", "site") sims_output <- readRDS("../results/allsims_plot_data.RDS") pred_table <- data.frame(site = table_data$Site, exp_demo = exp_synchrony_demo$pgr, sim_demo = round(subset(sims_output, simulation=="6No Comp. + No E.S.")["synchrony"],2), exp_env = round(exp_synchrony_env$growthrate_prediction,2), sim_env = round(subset(sims_output, simulation=="2No Comp. + No D.S.")["synchrony"],2)) colnames(pred_table) <- c("Site", "Predicted $\\phi_{R,\\mathcal{M}_D}$", "Simulated $\\phi_{R,\\mathcal{M}_D}$", "Predicted $\\phi_{R,\\mathcal{M}_E}$", "Simulated $\\phi_{R,\\mathcal{M}_E}$") myorder <- c(5,1,3,4,2) xpred_table <- pred_table[myorder,] pred_cap <- "Comparisons between our analytical predictions and simulation results for synchrony of species' per capita growth rates. Analytical predictions represent two limiting cases where only demographic stochasticity is operating ($\\phi_{R,\\mathcal{M}_D}$) and where only environmental stochasticity is operating ($\\phi_{R,\\mathcal{M}_E}$). Simulated synchrony values come from our empirically-based, multi-species population models when simulated under conditions that match the limiting case conditions (e.g., environmental stochasticity and competition removed for $\\mathcal{M}_D$)." print(xtable(xpred_table, caption = pred_cap), caption.placement="top", include.rownames = F, sanitize.colnames.function = identity, comment=FALSE, size="normalsize")
library(reshape2) library(plyr) library(ggplot2) library(xtable) perc_changes <- read.csv("../results/synchsims_percent_diffs.csv") perc_changes <- perc_changes[which(perc_changes$simulation!="3All Drivers"), c("site", "simulation", "percent_diff")] table_diff <- dcast(perc_changes, simulation~site, value.var = "percent_diff") diff_cap <- "Percent differences of synchrony of per capita growth rates between each removal simulation experiment and the 'All Drivers' simulation." print(xtable(table_diff, caption = diff_cap), caption.placement="top", include.rownames = FALSE, sanitize.colnames.function = identity, comment=FALSE)
library(reshape2) library(plyr) library(ggplot2) library(xtable) #### #### Read in regression fits -------------------------------------------------- #### all_files <- list.files("../results/") param_files1 <- all_files[grep("param*", all_files)] param_files <- param_files1[grep("RDS", param_files1)] num_vitals <- length(param_files) list_yr_effects <- list() for(j in 1:num_vitals){ param_list <- readRDS(paste0("../results/",param_files[j])) vital_name <- unlist(strsplit(param_files[j], "_"))[1] if(vital_name != "recruit"){ site_names <- names(param_list) site_year_list <- list() for(do_site in site_names){ tmp_site <- param_list[[do_site]] num_spp <- length(tmp_site) tmp_yr_int <- matrix(ncol = num_spp, nrow=nrow(tmp_site[[1]])) tmp_yr_slope <- matrix(ncol = num_spp, nrow=nrow(tmp_site[[1]])) for(i in 1:num_spp){ spp_df <- tmp_site[[i]] tmp_yr_int[,i] <- spp_df$Intercept.yr col2use <- grep("logarea", colnames(spp_df)) col2use2 <- grep(".yr", colnames(spp_df[col2use])) tmp_yr_slope[,i] <- spp_df[,col2use[col2use2]] } # end species within site loop site_year_list[[do_site]][["intercept"]] <- tmp_yr_int site_year_list[[do_site]][["slope"]] <- tmp_yr_slope } # end site within vital rate loop } # end NOT recruit loops if(vital_name == "recruit"){ site_names <- names(param_list) site_year_list <- list() for(do_site in site_names){ tmp_site <- param_list[[do_site]] tmp_yrs <- tmp_site[grep("intcpt.yr", rownames(tmp_site)),"Mean"] stringlist <- strsplit(names(tmp_yrs), ",") num_spp <- length(unique(sapply(stringlist, "[[", 2))) mat_yrs <- t(matrix(tmp_yrs, nrow=num_spp, byrow = TRUE)) site_year_list[[do_site]] <- mat_yrs } # end site within vital rate loop } # end YES recruit loops list_yr_effects[[vital_name]] <- site_year_list } # end vital rate loop # ## Plot histograms of year effects by vital rate # vital_rates <- names(list_yr_effects) # par(mfrow=c(3,5)) # for(do_vital in vital_rates){ # tmp_vital <- list_yr_effects[[do_vital]] # site_names <- names(tmp_vital) # for(do_site in site_names){ # tmp_site <- tmp_vital[[do_site]] # hist(tmp_site, main=paste(do_site,do_vital), xlab="Year Effect", xlim=c(-3,8)) # } # } ## Get standard deviation of year effect posterior #### #### Calculate average correlation of year effects ---------------------------- #### vital_rates <- names(list_yr_effects) out_df <- data.frame(site=NA, vital_rate=NA, term=NA, corr=NA) for(do_vital in vital_rates){ tmp_vital <- list_yr_effects[[do_vital]] site_names <- names(tmp_vital) for(do_site in site_names){ tmp_site <- tmp_vital[[do_site]] if(do_vital != "recruit"){ for(i in 1:length(tmp_site)){ tmp_cor <- cor(tmp_site[[i]]) avg_cor <- mean(tmp_cor[upper.tri(tmp_cor)]) tmp_out <- data.frame(site=do_site, vital_rate=do_vital, term=names(tmp_site)[i], corr=avg_cor) out_df <- rbind(out_df, tmp_out) } } if(do_vital == "recruit"){ tmp_cor <- cor(tmp_site) avg_cor <- mean(tmp_cor[upper.tri(tmp_cor)]) tmp_out <- data.frame(site=do_site, vital_rate=do_vital, term="intercept", corr=avg_cor) out_df <- rbind(out_df, tmp_out) } } } cor_df <- out_df[2:nrow(out_df),] cor_cast <- dcast(cor_df, site+term~vital_rate) cor.cap <- "Correlations of species' year random effects for each site by term, where term refers to the random effect on the slope or the intercept." print(xtable(cor_cast, caption = cor.cap), caption.placement="top", include.rownames = FALSE, sanitize.colnames.function = identity, comment=FALSE)
library(reshape2) library(plyr) library(ggplot2) library(xtable) #### #### Get Vital Rate Satistical Results Files ---------------------------------- #### vital_rates <- c("growth", "surv", "recruit") result_files <- list.files("../results/") vr_files1 <-result_files[grep(paste(vital_rates,collapse="|"), list.files("../results/"))] vr_files <- vr_files1[grep("RDS", vr_files1)] #### #### Loop Over Files and Extract Crowding Effects ----------------------------- #### # Growth and survival first vital_mat_list <- list() for(do_vital in vital_rates[c(1,2)]){ vr_do <- vr_files[grep(do_vital,vr_files)] torms <- grep("yr", vr_do) vr_do <- vr_do[-torms] tmp_vital <- readRDS(paste0("../results/",vr_do)) sites <- names(tmp_vital) mat_list <- list() for(do_site in sites){ tmp <- tmp_vital[[do_site]] spp <- names(tmp) tmpmat <- matrix(0,nrow = length(spp), ncol = length(spp)) for(i in 1:length(spp)){ do_spp <- spp[i] tmpspp <- tmp[[do_spp]] tmpw <- tmpspp[1,grep("crowd*", colnames(tmpspp))] tmpmat[i,] <- as.numeric(tmpw) } rownames(tmpmat) <- spp colnames(tmpmat) <- spp mat_list[[do_site]] <- tmpmat } vital_mat_list[[do_vital]] <- mat_list } # Recruitment vr_do <- vr_files[grep("recruit",vr_files)] torms <- c(grep("yr", vr_do), grep("MONT", vr_do)) vr_do <- vr_do[-torms] tmp_vital <- readRDS(paste0("../results/",vr_do)) names_vital <- readRDS(paste0("../results/",vr_files[1])) sites <- names(tmp_vital) mat_list <- list() for(do_site in sites){ tmp <- tmp_vital[[do_site]] spp <- names(names_vital[[do_site]]) tmpdd <- tmp[grep("dd", rownames(tmp)),"Mean"] tmpmat <- matrix(tmpdd,length(spp),length(spp), byrow = FALSE) rownames(tmpmat) <- spp colnames(tmpmat) <- spp mat_list[[do_site]] <- tmpmat } vital_mat_list[["recruit"]] <- mat_list #### #### Calculate Average Inter and Intraspecific Interactions ------------------- #### vital_comp_list <- list() for(do_vital in vital_rates){ tmp <- vital_mat_list[[do_vital]] sites <- names(tmp) comp_list <- list() for(do_site in sites){ tmpsite <- tmp[[do_site]] avgintra <- mean(diag(tmpsite)) avginter <- mean(c(tmpsite[upper.tri(tmpsite)], tmpsite[lower.tri(tmpsite)])) avgcomp <- matrix(c(avginter, avgintra),1,2) colnames(avgcomp) <- c("inter", "intra") comp_list[[do_site]] <- avgcomp } vital_comp_list[[do_vital]] <- comp_list } df <- melt(unlist(vital_comp_list)) pieces <- strsplit(rownames(df), split = "[.]") df$vitalrate <- sapply(pieces, "[", 1) sites <- sapply(pieces, "[", 2) df$site <- substr(sites, start = 1, stop = nchar(sites)-1) df$interintra <- substr(sites, start = nchar(sites), stop = nchar(sites)) dftable <- dcast(df, site+interintra~vitalrate, value.var = "value") colnames(dftable) <- c("Site", "Interaction Type", "Growth", "Recruitment", "Survival") dftable[which(dftable$`Interaction Type`==1),"Interaction Type"] <- "Interspecific" dftable[which(dftable$`Interaction Type`==2),"Interaction Type"] <- "Intraspecific" # avg_interintra <- ddply(df, .(site, interintra), summarise, # avg_val = mean(value)) # avg_interintra[which(avg_interintra$interintra==1),"interintra"] <- "Interspecific Competition" # avg_interintra[which(avg_interintra$interintra==2),"interintra"] <- "Intraspecific Competition" # interintra_table <- dcast(avg_interintra, site~interintra, value.var = "avg_val") comp_cap <- "Average interaction coefficients for each vital rate for each community." print(xtable(dftable, caption = comp_cap, digits=4), caption.placement="top", include.rownames = FALSE, sanitize.colnames.function = identity, comment=FALSE)
#### #### Get Vital Rate Satistical Results Files ---------------------------------- #### vital_rates <- c("growth", "surv", "recruit") result_files <- list.files("../results/") vr_files <-result_files[grep(paste(vital_rates,collapse="|"), list.files("../results/"))] #### #### Loop Over Files and Extract Crowding Effects ----------------------------- #### # Growth and survival first vital_mat_list <- list() for(do_vital in vital_rates[c(1,2)]){ vr_do <- vr_files[grep(do_vital,vr_files)] torms <- grep("yr", vr_do) vr_do <- vr_do[-torms] tmp_vital <- readRDS(paste0("../results/",vr_do)) sites <- names(tmp_vital) mat_list <- list() for(do_site in sites){ tmp <- tmp_vital[[do_site]] spp <- names(tmp) tmpmat <- matrix(0,nrow = length(spp), ncol = length(spp)) for(i in 1:length(spp)){ do_spp <- spp[i] tmpspp <- tmp[[do_spp]] tmpw <- tmpspp[1,grep("crowd*", colnames(tmpspp))] tmpmat[i,] <- as.numeric(tmpw) } rownames(tmpmat) <- spp colnames(tmpmat) <- spp mat_list[[do_site]] <- tmpmat } vital_mat_list[[do_vital]] <- mat_list } # Recruitment vr_do <- vr_files[grep("recruit",vr_files)] vr_do <- vr_do[which(vr_do=="recruit_parameters.RDS")] tmp_vital <- readRDS(paste0("../results/",vr_do)) names_vital <- readRDS(paste0("../results/",vr_files[1])) sites <- names(tmp_vital) mat_list <- list() for(do_site in sites){ tmp <- tmp_vital[[do_site]] spp <- names(names_vital[[do_site]]) tmpdd <- tmp[grep("dd", rownames(tmp)),"Mean"] tmpmat <- matrix(tmpdd,length(spp),length(spp), byrow = FALSE) rownames(tmpmat) <- spp colnames(tmpmat) <- spp mat_list[[do_site]] <- tmpmat } vital_mat_list[["recruit"]] <- mat_list
names(vital_mat_list) <- c("growth", "survival", "recruitment") for(do_vital in names(vital_mat_list)){ tmp_list <- vital_mat_list[[do_vital]] for(do_site in names(tmp_list)){ site_list <- tmp_list[[do_site]] synch_cap <- paste0("Interaction coefficients for ", do_vital, " regressions in ", do_site, ".") # digitvec <- rep(3,ncol(site_list)) print(xtable(site_list, caption = synch_cap, digits = 4), caption.placement="top", include.rownames = TRUE, sanitize.colnames.function = identity, comment=FALSE) } }
\newpage{}
References
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.