\SweaveOpts{ eval=TRUE, prefix=FALSE, width=4.5, height=4.5, eps=FALSE, include=FALSE} \chapter{Community Composition and Diversity} <>= setwd("~/Documents/Projects/Book/draft/Chap10") rm(list=ls()) badfonts = FALSE .detach.stuff = function() { s1 = grep("package|Autoloads",search()) nattach = length(search()) xattach = search()[c(-1,-s1)] for (i in xattach) eval(substitute(detach(i),list(i=i))) } .detach.stuff() ## ps.options(family="NimbusSan") ps.options(colormodel="cmyk") options(width=60, SweaveHooks=list(fig=function() par(mar=c(5,4,.5,1.5))) ) ## palette(gray((0:8)/8)) library(primer) trellis.par.set(canonical.theme(color=FALSE)) @ <>= bsdiv <- read.csv('C3rank.csv') matplot(bsdiv[,1], bsdiv[,c(2:4, 7,9,11)], type='l', lty=1:3, lwd=rep(1:2, each=3), col=1, log="y", xlab="Abundance Rank", ylab="Relative Abundance") legend('topright', paste("Year ", c(1,2,4,15,25,35), sep=""), ncol=2, lty=1:3, lwd=rep(1:2, each=3), bty='n') ?moths @ \begin{figure} \centering \includegraphics[width=.6\textwidth]{BSsuc.pdf} \caption{Empirical rank--abundance distributions of successional plant communities (old-fields) within the temperate deciduous forest biome of North America. ``Year'' indicates the time since abandonment from agriculture. Data from the Buell-Small succession study (http://www.ecostudies.org/bss/)}\label{fig:BSsuc} \end{figure} It seems easy, or at least tractable, to compare the abundance of a single species in two samples. In this chapter, we introduce concepts that ecologists use to compare entire communities in two samples. We focus on two quantities: \emph{species composition}, and \emph{diversity}. We also discuss several issues related to this, including species--abundance distributions, ecological neutral theory, diversity partitioning, and species--area relations. Several packages in \R~include functions for dealing specifically with these topics. Please see the ``Environmetrics'' link within the ``Task Views'' link at any CRAN website for downloading \R packages. Perhaps the most comprehensive (including both diversity and composition) is the \texttt{vegan} package, but many others include important features as well. \section{Species Composition} \index{species composition}Species composition is merely the set of species in a site or a sample. Typically this includes some measure of abundance at each site, but it may also simply be a list of species at each site, where ``abundance'' is either presence or absence. Imagine we have four sites (A--D) from which we collect density data on two species, \emph{Salix whompii} and \emph{Fraxinus virga}. We can enter hypothetical data of the following abundances. <<>>= dens <- data.frame(Salwho=c(1,1,2,3), Fravir=c(21,8,13,5)) row.names(dens) <- LETTERS[1:4] dens @ Next, we plot the abundances of both species; the plotted points then are the sites (Fig. \ref{fig:dens}). <>= plot(dens, type='n'); text(dens, row.names(dens)) @ In Fig. \ref{fig:dens}, we see that the species composition in site A is most different from the composition of site D. That is, the distance between site A and D is greater than between any other sites. The next question, then, is \emph{how} far apart are any two sites? Clearly, this depends on the scale of the measurement (e.g., the values on the axes), and also on how we measure distance through multivariate space. \begin{figure} \centering \includegraphics[width=.5\textwidth]{dens.pdf} \caption{ Hypothetical species composition for four sites (A--D).} \label{fig:dens} \end{figure} \subsection{Measures of abundance} Above we pretended that the abundances were absolute densities (i.e., 1 = one stem per sample). We could of course represent all the abundances differently. For instance, we could calculate \emph{relative density}, where each species in a sample is represented by the \emph{proportion} of the sample comprised of that species. For site A, we divide each species by the sum of all species. <<>>= dens[1,]/sum(dens[1,]) @ We see that \emph{Salix} makes up about 5\% of the sample for Site A, and \emph{Fraxinus} makes up about 95\% of the sample. Once we calculate relative densities for each species at each site, this eliminates differences in total density at each site because all sites then total to 1. We could also calculate relative measures for any type of data, such as biomass or percent cover. In most instances, \emph{\index{relative density}relative density} refers to the density of a species \emph{relative} to the other species in a sample (above), but it can also be density in a sample relative to other samples. We would thus make each \emph{species} total equal 1, and then its abundance at each site reflects the proportion of a species total abundance comprised by that site. For instance, we can make all \emph{Salix} densities relative to each other. <<>>= dens[,1]/sum(dens[,1]) @ Here we see that sites A and B both have about 14\% of all \emph{Salix} stems, and site D has 43\%. Whether our measures of abundance are absolute or relative, we would like to know how different samples (or sites) are from each other. Perhaps the simplest way to describe the difference among the sites is to calculate the \emph{distances} between each pair of sites. \subsection{Distance} There are many ways to calculate a \emph{\index{distance}distance} between a pair of sites. One of the simplest, which we learned in primary school, is \emph{\index{Euclidean distance|see{distance}}\index{distance!Euclidean distance}Euclidean distance}. With two species, we have two dimensional space, which is Fig. \ref{fig:dens}. The Euclidean distance between two sites is merely the length of the vector connecting those sites. We calculate this as $\sqrt{x^2+y^2}$, where $x$ and $y$ are the $(x,y)$ distances between a pair of sites. The $x$ distance between sites B and C is difference in \emph{Salix} abundance between the two sites, <<>>= x <- dens[2,1] - dens[3,1] @ where \texttt{dens} is the data frame with sites in rows, and species in different columns. The $y$ distance between sites B and C is difference in \emph{Fraxinus} abundance between the two sites. <<>>= y <- dens[2,2] - dens[3,2] @ The Euclidean distance between these is therefore <<>>= sqrt(x^2 + y^2) @ Distance is as simple as that. We calculate all pairwise Euclidean distances between sites A--D based on 2 species using built-in functions in \R. <<>>= (alldists <- dist(dens)) @ We can generalize this to include any number of species, but it becomes increasingly harder to visualize. We can add a third species, \emph{Mandragora officinarum}, and recalculate pairwise distances between all sites, but now with three species. <>= dens[["Manoff"]] <- c(11,3,7,5) (spp3 <- dist(dens)) @ We can plot species abundances as we did above, and \texttt{pairs(dens)} would give us all the pairwise plots given three species. However, what we really want for species is a 3-D plot. Here we load another package\footnote{You can install this package from any R CRAN mirror.} and create a 3-D scatterplot. <>= pairs(dens)# not shown library(scatterplot3d) sc1 <- scatterplot3d(dens, type='h', pch="", xlim=c(0,5), ylim=c(0, 25), zlim=c(0,15)) text(sc1$xyz.convert(dens), labels=rownames(dens)) @ In three dimensions, Euclidean distances are calculated the same basic way, but we add a third species, and the calculation becomes $\sqrt{x^2 + y^2 +z^2}$. Note that we take the square root (as opposed to the cube root) because we originally \emph{squared} each distance. We can generalize this for two sites for $R$ species as \begin{gather} \label{eq:euc} D_E = \sqrt{\sum_{i=1}^R\left( x_{ai} - x_{bi}\right)^2} \end{gather} Of course, it is difficult (impossible?) to visualize arrangements of sites with more than three axes (i.e., $>3$ species), but we can always \emph{calculate} the distances between pairs of sites, regardless of how many species we have. There are many ways, in addition to Euclidean distances, to calculate distance. Among the most commonly used in ecology is \emph{Bray--Curtis} distance, which goes by other names, including \emph{\index{Sorenson distance@S{\o}orenson distance|see{distance}}\index{Sorenson distance@S{\o}orenson distance!distance}S{\o}orenson} distance. \begin{equation} \label{eq:sor} D_{BC} = \sum_{i=1}^R\frac{\left| x_{ai} - x_{bi}\right|}{x_{ai} + x_{bi}} \end{equation} where $R$ is the number of species in all samples. \index{Bray--Curtis distance|see{distance}}\index{distance!Bray--Curtis}Bray--Curtis distance is merely the total difference in species abundances between two sites, divided by the total abundances at each site. Bray--Curtis distance (and a couple others) tends to result in more intuitively pleasing distances in which both common and rare species have relatively similar weights, whereas Euclidean distance depends more strongly on the most abundant species. This happens because Euclidean distances are based on squared differences, whereas Bray--Curtis uses absolute differences. Squaring always amplifies the importance of larger values. Fig. \ref{fig:mds} compares graphs based on Euclidean and Bray--Curtis distances of the same raw data. \subsubsection{Displaying \index{multidimensional distance|see{distance}} \index{multidimensional distance!distance}multidimensional distances} A simple way to display distances for three or more species is to create a plot in two dimensions that attempts to arrange all sites so that they are \emph{approximately} the correct distances apart. In general this is impossible to achieve precisely, but distances can be approximately correct. One technique that tries to create an optimal (albiet approximate) arrangement is \emph{\index{non-metric multidimensional scaling}non-metric multidimensional scaling}. Here we add a fourth species (\emph{Aconitum lycoctonum}) to our data set before plotting the distances. <<>>= dens$Acolyc <- c(16, 0, 9,4) @ The non-metric multidimensional scaling function is in the \texttt{vegan} package. It calculates distances for us using the original data. Here we display Euclidean distances among sites (Fig. \ref{mds1}). <>= library(vegan) mdsE <- metaMDS(dens, distance="euc", autotransform=FALSE, trace=0) plot(mdsE, display="sites", type="text") @ Here we display Bray--Curtis distances among sites (Fig. \ref{mds2}). <>= mdsB <- metaMDS(dens, distance="bray", autotransform=FALSE, trace=0) plot(mdsB, display="sites", type="text") @ \begin{figure} \centering \subfloat[Euclidean distances] { \includegraphics[width=.47\textwidth]{mds1} \label{mds1}} \subfloat[Bray--Curits distances] { \includegraphics[width=.47\textwidth]{mds2} \label{mds2}} \caption{Nonmetric multidimensional (NMDS) plots showing approximate distances between sites. These two figures display the same raw data, but Euclidean distances tend to emphasize differences due to the more abundant species, whereas Bray-Curtis does not. Because NMDS provides iterative optimizations, it will find slightly different arrangements each time you run it.}\label{fig:mds} \end{figure} \subsection{Similarity} Sometimes, we would like to know how \emph{similar} two communities are. Here we describe two measures of \index{similarity}similarity, \emph{percent similarity}, and \emph{S{\o}rensen similarity} \cite{Magurran2004}. Percent similarity may be the simplest of these; it is simply the sum of the minimum percentages for each species in the community. Here we convert each species to its relative abundance; that is, its proportional abundance at each site. To do this, we treat each site (row) separately, and then divide the abundance of each species by the sum of the abundances at each site. <<>>= (dens.RA <- t( apply(dens, 1, function(sp.abun) sp.abun/sum(sp.abun) )) ) @ Next, to compare two sites, we find the minimum relative abundance for each species. Comparing sites A and B, we have, <<>>= (mins <- apply(dens.RA[1:2,], 2, min)) @ Finally, we sum these, and multiply by 100 to get percentages. <<>>= sum(mins) *100 @ The second measure of similarity we investigate is \index{Sorensen's similarity@S{\o}rensen's similarity}S{\o}rensen's similarity, \begin{equation} \label{eq:ss} S_s = \frac{2C}{A+B} \end{equation} where $C$ is the number of species two sites have in common, and $A$ and $B$ are the number of species at each site. This is equivalent to dividing the shared species by the average richness. To calculate this for sites A and B, we could find the species which have non-zero abundances both sites. <<>>= (shared <- apply(dens[1:2,], 2, function(abuns) all(abuns != 0) ) ) @ Next we find the richness of each. <<>>= (Rs <- apply(dens[1:2,], 1, function(x) sum(x>0)) ) @ Finally, we divide the shared species by the summed richnesses and multiply by 2. <<>>= 2*sum(shared)/sum(Rs) @ S{\o}rensen's index has also been used in the development of more sophisticated measures of similarity between sites \cite{Morlon:2008tx,Plotkin:2002eu}. \section{Diversity} To my mind, there is no more urgent or interesting goal in ecology and evolutionary biology than understanding the determinants of biodiversity. \index{biodiversity|see{diversity}}Biodiversity is many things to many people, but we typically think of it as a measure of the variety of biological life present, perhaps taking into account the relative abundances. For ecologists, we most often think of \emph{\index{diversity}species diversity} as some quantitative measure of the variety or number of different species. This has direct analogues to the genetic diversity within a population \cite{Crow1970,Vellend2005a}, and the connections between species and genetic diversity include both shared patterns and shared mechanisms. Here we confine ourselves entirely to a discussion of species diversity, the variety of different species present. <>= H <- function(x) {x <- x[x>0] p <- x/sum(x) - sum (p * log(p) ) } Sd <- function(x) {p <- x/sum(x) 1 - sum(p^2)} h1 <- H(x=rep(.25,4)); h2<- H(x=c(50,75,75));h3<- H(c(20,20,20,140)); h4<- H(x=200) sd1 <- Sd(x=rep(.25,4)); sd2<- Sd(x=c(50,75,75));sd3<- Sd(c(20,20,20,140)); @ Consider this example. We have four stream insect communities (Table \ref{tab:inverts}). Which has the highest ``biodiversity''? \begin{table}[ht] \caption{Four hypothetical stream invertebrate communities. Data are total numbers of individuals collected in ten samples (sums across samples). Diversity indices (Shannon-Wiener, Simpson's) explained below.} \centering \begin{tabular}[c]{lcccc} \hline \hline Species & Stream 1 & Stream 2 & Stream 3 & Stream 4\\ \hline \emph{Isoperla} & 20 & 50 & 20 & 0\\ \emph{Ceratopsyche} & 20 & 75 & 20 & 0\\ \emph{Ephemerella} & 20 & 75 & 20 & 0\\ \emph{Chironomus} & 20 & 0 & 140 & 200\\ \hline Number of species ($R$)&4& 3 & 4 & 1\\ Shannon-Wiener $H$ & \Sexpr{round(h1,2)} & \Sexpr{round(h2,2)} & \Sexpr{round(h3,2)} & 0 \\ Simpson's $S_D$ & \Sexpr{round(sd1,2)} & \Sexpr{round(sd2,2)} & \Sexpr{round(sd3,2)} & 0 \\ \hline \end{tabular} \label{tab:inverts} \end{table} We note that one stream has only one species --- clearly that can't be the most ``diverse'' (still undefined). Two streams have four species --- are they the most diverse? Stream 3 has more bugs in total (200 \emph{vs.} 80), but stream 1 has a more equal distribution among the four species. \subsection{Measurements of variety} So, how shall we define ``diversity'' and measure this variety? There are many mathematical expressions that try to summarize biodiversity \cite{Hill:1973it, Keylock:2005fy}. The inquisitive reader is referred to \cite{Magurran2004} for a practical and comprehensive text on measures of species diversity. Without defining it precisely (my pay scale precludes such a noble task), let us say that diversity indices attempt to quantify \begin{itemize} \item the probability of encountering different species at random, or, \item the uncertainty or multiplicity of possible community states (i.e., the \index{entropy}entropy of community composition), or, \item the \index{variance in species composition}variance in species composition, relative to a multivariate centroid. \end{itemize} For instance, a simple count of species (Table \ref{tab:inverts}) shows that we have 4, 3, 4, and 1 species collected from streams 1--4. The larger the number of species, the less certain we could be about the identity of an individual drawn blindly and at random from the community. To generalize this further, imagine that we have a \index{species pool}species pool\footnote{A \emph{species pool} is the entire, usually hypothetical, set of species from which a sample is drawn; it may be all of the species known to occur in a region.} of $R$ species, and we have a sample comprised of only one species. In a sample with only one species, then we know that the possible \emph{states} that sample can take is limited to one of only \emph{R} different possible states --- the abundance of one species is 100\% and all others are zero. On the other hand, if we have two species then the community could take on $R\left(R-1\right)$ different states --- the first species could be any one of $R$ species, and the second species could be any one of the other species, and all others are zero. Thus increasing diversity means increasing the possible states that the community could take, and thus increasing our uncertainty about community structure \cite{Keylock:2005fy}. This increasing lack of information about the system is a form of \emph{entropy}, and increasing diversity (i.e., increasing multiplicity of possible states) is increasing entropy. The jargon and topics of \index{statistical mechanics}statistical mechanics, such as \index{entropy}entropy, appear (in 2009) to be an increasingly important part of community ecology \cite{Harte:2008nq, Pueyo:2007ow}. Below we list three commonly used diversity indices: species \index{richness}richness, \index{Shannon-Wiener}Shannon-Wiener index, and \index{Simpson's diversity}Simpson's diversity index. \begin{itemize} \item Species richness, $R$, the count of the number of species in a sample or area; the most widely used measure of biodiversity \cite{Hurlbert1971}. \item Shannon-Wiener diversity.\footnote{Robert May stated that this index is connected by merely an ``ectoplasmic thread'' to information theory \cite{May1976}, but there seems be a bit more connection than that.} \begin{equation} \label{eq:H1} H'=-\sum_{i=1}^{R}p\ln\left(p\right) \end{equation} where $R$ is the number of species in the community, and $p_i$ is the relative abundance of species $i$. \item \index{Simpson's diversity}Simpson's diversity. This index is (i) the probability that two individuals drawn from a community at random will be different species \cite{Nei1987}, (ii) the initial slope of the species-individuals curve \cite{Lande2000} (e.g., Fig. \ref{fig:est}), and (iii) the expected variance of species composition (Fig. \ref{fig:Simp}) \cite{Lande1996,Stevens:2003nc}. \begin{gather} S_D=1-\sum_{i=1}^{R}p_i^{2} \end{gather} The summation $\sum_{i=1}^{R}p_i^{2}$ is the probability that two individuals drawn at random are the same species, and it is known as Simpson's ``\index{dominance}dominance.'' Lande found that this Simpon's index can be more precisely estimated, or estimated accurately with smaller sample sizes, than either richness or Shannon-Wiener \cite{Lande1996}. \end{itemize} These three indices are actually directly related to each other --- they comprise estimates of \emph{entropy}, the amount of disorder or the multiplicity of possible states of a system, that are directly related via a single constant \cite{Keylock:2005fy}. However, an important consequence of their differences is that richness depends most heavily on rare species, Simpson's depends most heavily on common species, and Shannon-Wiener stands somewhere between the two (Fig. \ref{fig:divs}). \begin{figure}[ht] \centering \subfloat[Equally abundant]{\includegraphics[width=.46\linewidth]{EqAb} \label{eqab}} \subfloat[Most common sp. = 90\%]{\includegraphics[width=.46\linewidth]{Ab90} \label{ab90}} \caption{Relations between richness, Shannon-Weiner, and Simpson's diversities (note difference in $y$-axis scale between the figures). Communities of 2--20 species are composed of either equally abundant species \subref{eqab} or with the most common species equal to 90\% of the community \subref{ab90}.} \label{fig:divs} \end{figure} \medskip \noindent \subsubsection{Relations between number of species, relative abundances, and diversity} \emph{This section relies heavily on code and merely generates Fig. \ref{fig:divs}.} Here we display diversities for communities with different numbers and relative abundances of species (Fig. \ref{fig:divs}). We first define functions for the diversity indices. <<>>= H <- function(x) {x <- x[x>0] p <- x/sum(x) - sum (p * log(p) ) } Sd <- function(x) {p <- x/sum(x) 1 - sum(p^2)} <>= h1 <- H(x=rep(.25,4)); h2<- H(x=c(50,75,75));h3<- H(c(1,1,1,197)); h4<- H(x=200) sd1 <- Sd(x=rep(.25,4)); sd2<- Sd(x=c(50,75,75));sd3<- Sd(c(1,1,1,197)); @ Next we create a \emph{list} of communities with from 1 to 20 \emph{equally} abundant species, and calculate $H$ and $S_D$ for each community. <>= Rs <- 2:20 ComsEq <- sapply(Rs, function(R) (1:R)/R) Hs <- sapply(ComsEq, H); Sds <- sapply(ComsEq, Sd) plot(Rs, Hs, type='l', ylab="Diversity", ylim=c(0,3)); lines(Rs, Sds, lty=2); legend("right", c(expression(italic("H")), expression(italic("S"["D"]))), lty=1:2, bty='n') @ Now we create a \emph{list} of communities with from 2 to 25 species, where one species always comprises 90\% of the community, and the remainder are equally abundant rare species. We then calculate $H$ and $S_D$ for each community. <>= Coms90 <- sapply(Rs, function(R) {p <- numeric(R) p[1] <- .9 p[2:R] <- .1/(R-1) p} ) Hs <- sapply(Coms90, H); Sds <- sapply(Coms90, Sd) plot(Rs, Hs, type='l', ylim=c(0,.7)); lines(Rs, Sds, lty=2) @ \subsubsection{\index{Simpson's diversity}Simpson's diversity, as a \index{variance in species composition}variance of composition} \emph{This section relies heavility on code.} Here we show how we would calculate the variance of species composition. First we create a pretend community of six individuals (rows) and 3 species (columns). Somewhat oddly, we identify the degree to which each individual is comprised of each species; in this case, individuals can be only one species.\footnote{Imagine the case where the columns are traits, or genes. In that case, individuals could be characterized by affiliation with multiple columns, whether traits or genes.} Here we let two individuals be each species. <<>>= s1 <- matrix(c(1,1,0,0,0,0,0,0,1,1,0,0,0,0,0,0,1,1), nr=6) colnames(s1) <- c("Sp.A", "Sp.B", "Sp.C");s1 @ We can plot these individuals in community space, if we like (Fig. \ref{simp1}). <>= library(scatterplot3d) s13d <- scatterplot3d(jitter(s1, .3), type='h', angle=60, pch=c(1,1,2,2,3,3), # main="Simpson's Div = 0.67", xlim=c(-.2, 1.4), ylim=c(-.2, 1.4), zlim=c(-.2, 1.4)) s13d$points3d(x=1/3, y=1/3, z=1/3, type='h', pch=19, cex=2) @ Next we can calculate a \emph{centroid}, or multivariate mean --- it is merely the vector of species means. <<>>= (centroid1 <- colMeans(s1)) @ Given this centroid, we begin to calculate a variance by (i) substracting each species vector (0s, 1s) from its mean, (ii) squaring each of these deviates, and (3) summing to get the sum of squares. <<>>= (SS <- sum( sapply(1:3, function(j) (s1[,j] - centroid1[j]))^2) ) @ We then divide this sum by the number of individuals that were in the community ($N$) <<>>= SS/6 @ We find that the calculation given above for Simpson's diversity returns exactly the same number. We would calculate the relative abundances, square them, add them, and subtract that value from 1. <<>>= p <- c(2,2,2)/6 1-sum(p^2) @ In addition to being the variance of species composition, this number is also the probability that two individuals drawn at random are different species. As we mentioned above, there are other motivations than these to derive this and other measures of species diversity, based on entropy and information theory \cite{Keylock:2005fy} --- and they are all typically correlated with each other. <>= s2 <- matrix(c(1,1,1,1,0,0,0,0,0,0,1,0,0,0,0,0,0,1), nr=6) colnames(s2) <- c("Sp.A", "Sp.B", "Sp.C");s2 centroid2 <- colMeans(s2) s23d <- scatterplot3d(jitter(s2, .2), type='h', angle=60, pch=c(1,1,1,1,2,3), # main="Simpson's Div = 0.50", xlim=c(-.2, 1.4), ylim=c(-.2, 1.4),zlim=c(-.2, 1.4)) s23d$points3d(x=4/6, y=1/6, z=1/6, type='h', pch=19, cex=2) <>= s3 <- matrix(c(1,1,1,1,1,0,0,0,0,0,0,1,0,0,0,0,0,0), nr=6) colnames(s3) <- c("Sp.A", "Sp.B", "Sp.C");s3 (centroid3 <- colMeans(s3)) s33d <- scatterplot3d(jitter(s3, .3), type='h', angle=60, pch=c(1,1,1,1,1,2), xlim=c(-.2, 1.4), ylim=c(-.2, 1.4), zlim=c(-.2, 1.4) # main="Simpson's Div = 0.28" ) s33d$points3d(x=5/6, y=1/6, z=0, type='h', pch=19, cex=2) @ \begin{figure}[ht] \centering \subfloat[$S_D = 0.67$]{% \includegraphics[width=.32\textwidth]{simpson1.pdf} \label{simp1}} \subfloat[$S_D = 0.50$]{% \includegraphics[width=.32\textwidth]{simpson2.pdf} \label{simp2}} \subfloat[$S_D = 0.28$]{% \includegraphics[width=.32\textwidth]{simpson3.pdf} \label{simp3}} \caption{Plotting three examples of species composition. The centroid of each composition is a solid black dot. The third example (on right) has zero abundances of species C. Simpson's diversity is the variance of these points around the centroid. Individual points are not plotted at precisely 0 or 1 --- they are plotted with a bit of jitter or noise so that they do not overlap entirely.\label{fig:Simp}} \end{figure} % \begin{table}[ht] % \caption{Simpson's Diversity --- an example using three communities (\textbf{A--C}) and three species (1--3); there are six different individuals in each community (total of 18 different individuals).} % \centering % \begin{tabular}[c]{lccc|c|ccc|c|ccc} % & & \textbf{A} & &~& & \textbf{B }& &~& & \textbf{C} & \\ \hline % Individual & Sp. 1 & Sp. 2 & Sp. 3 &~& Sp. 1 & Sp. 2 & Sp. 3 &~& Sp. 1 & Sp. 2 & Sp. 3 \\ \hline % 1 & 1 & 0 & 0& & 1 & 0 & 0&& & 1 & 0 & 0\\ % 2 & 1 & 0 & 0& & 1 & 0 & 0&& & 1 & 0 & 0\\ % 3 & 0 & 1 & 0& & 1 & 0 & 0&& & 1 & 0 & 0\\ % 4 & 0 & 1 & 0& & 1 & 0 & 0&& & 1 & 0 & 0\\ % 5 & 0 & 0 & 1& & 0 & 1 & 0&& & 1 & 0 & 0\\ % 6 & 0 & 0 & 1& & 0 & 0 & 1&& & 0 & 1 & 0\\ % \hline % \end{tabular} % \label{tab:simps} % \end{table} @ \subsection{Rarefaction and total species richness} \index{rarefaction}Rarefaction is the process of generating the relationship between the number of species \emph{vs.} the number of individuals in one or more samples. It is typically determined by randomly resampling individuals \cite{Gotelli2001a}, but could also be determined by resampling samples. Rarefaction allows direct comparison of the richness of two samples, corrected for numbers of individuals. This is particularly important because $R$ depends heavily on the number of individuals in a sample. Thus rarefaction finds the general relation between the number(s) of species \emph{vs.} number of individuals (Fig. \ref{fig:est}), and is limited to less than or equal to the number of species you actually observed. A related curve is the \emph{\index{species-accumulation curve}species-accumulation} curves, but this is simply a useful but haphazard accumulation of new species (a cumulative sum) as the investigator samples new individuals. Another issue that ecologists face is trying to estimate the \emph{true} number of species in an area, given the samples collected. This number of species would be larger than the number of species you observed, and is often referred to as \emph{\index{total species richness}total species richness} or the \emph{\index{asymptotic richness}asymptotic richness}. Samples almost always find only a subset of the species present in an area or region, but we might prefer to know how many species are really there, in the general area we sampled. There are many ways to do this, and while some are better than others, none is perfect. These methods estimate minimum numbers of species, and assume that the unsampled areas are homogeneous and similar to the sampled areas. Before using these methods seriously, the inquisitive reader should consult \cite{Gotelli2001a, Magurran2004} and references at http://viceroy.eeb.uconn.edu/EstimateS. Below, we briefly explore an example in \R. \subsubsection{An example of rarefaction and total species richness} Let us ``sample'' a seasonal tropical rainforest on Barro Colorado Island (BCI) http://ctfs.si.edu/datasets/bci/). Our goal will be to provide baseline data for later comparison to other such studies. We will use some of the data from a 50\,ha plot that is included in the \texttt{vegan} package \cite{Condit2002, Oksanen:2008}. We will pretend that we sampled every tree over 10\,cm dbh,\footnote{``dbh'' is diameter at 1.37\,m above the ground.} in each of 10 plots scattered throughout the 50\,ha area. What could we say about the forest within which the plots were located? We have to consider the scale of the sampling. Both the \emph{\index{experimental unit}experimental unit} and the \emph{\index{grain}grain} are the 1\,ha plots. Imagine that the plots were scattered throughout the 50\,ha plot, so that the extent of the sampling was a full 50.\,ha\footnote{See John Wiens' foundational paper on spatial \index{scale}scale in ecology \cite{Wiens1989} describing the meaning of grain, extent, and other spatial issues.} <>= ### Here I play around with stuff library(vegan) data(BCI); N <- sort(colSums(BCI), decr=TRUE) subs1 <- c( seq(50, 400, by=50), sum(BCI[1,])) rar1 <- rarefy(BCI[1,], MARGIN=1, sample=subs1, se=T) subs2 <- c( seq(50, 400, by=50), sum(BCI[2,])) rar2 <- rarefy(BCI[2,], MARGIN=1, sample=subs2, se=T) plot(subs1, rar1[1,], type='b', ylim=c(0, 100), ylab="Individual Rarefied Richness", xlab="No. of Individuals") lines(subs1, rar1[1,]+2*rar1[2,], lty=3) lines(subs1, rar1[1,]-2*rar1[2,], lty=3) lines(subs2, rar2[1,], type="b") @ First, let's pretend we have sampled 10 1\,ha plots by extracting the samples out of the larger dataset. <<>>= library(vegan) data(BCI); bci <- BCI[seq(5,50, by=5),] @ Next, for each species, I sum all the samples into one, upon which I will base rarefaction and total richness estimation. <>= N <- colSums(bci) NoOfInds <- c(seq(500, 4500, by=500), sum(N)) rar3 <- rarefy(N, sample=NoOfInds, se=T, MARG=2) plot(NoOfInds, rar3[1,], ylab="Species Richness", axes=FALSE, xlab="No. of Individuals", type='n', ylim=c(0,250), xlim=c(500,7000) ) axis(1, at=1:5*1000); axis(2); box(); text(2500, 200, "Individual-based rarefaction (10 plots)") lines(NoOfInds, rar3[1,], type='b') lines(NoOfInds, rar3[1,]+2*rar3[2,], lty=3) lines(NoOfInds, rar3[1,]-2*rar3[2,], lty=3) ace <- estimateR(N) segments(6000, ace['S.ACE']-2*ace['se.ACE'], 6000, ace['S.ACE']+2*ace['se.ACE'], lwd=3) text(6000, 150, "ACE estimate", srt=90, adj=c(1,.5) ) chaoF <- specpool(bci) segments(6300, chaoF[1,'Chao'] - 2*chaoF[1,'Chao.SE'], 6300, chaoF[1,'Chao'] + 2*chaoF[1,'Chao.SE'], lwd=3, col='grey') text(6300, 150, "Chao2 estimate", srt=90, adj=c(1,.5) ) points(6700, dim(BCI)[2], pch=19, cex=1.5) text(6700, 150, "Richness observed in 50 ha", srt=90, adj=c(1,.5) ) @ Next we combine all the plots into one sample (a single summed count for each species present), select numbers of individuals for which I want rarefied samples (multiples of 500), and then perform rarefaction for each of those numbers. <<>>= N <- colSums(bci) subs3 <- c(seq(500, 4500, by=500), sum(N)) rar3 <- rarefy(N, sample=subs3, se=T, MARG=2) @ Next we want to graph it, with a few bells and whistles. We set up the graph of the 10 plot, individual-based rarefaction, and leave room to graph total richness estimators as well (Fig. \ref{fig:est}). <<>>= plot(subs3, rar3[1,], ylab="Species Richness", axes=FALSE, xlab="No. of Individuals", type='n', ylim=c(0,250), xlim=c(500,7000) ) axis(1, at=1:5*1000); axis(2); box(); text(2500, 200, "Individual-based rarefaction (10 plots)") @ \begin{figure}[ht] \centering \includegraphics[width=.6\linewidth]{est} \caption{Baseline tree species richness estimation based on ten 1\,ha plots, using individual-based rarefaction, and two different total richness estimators, ACE and Chao 2. The true total tree richness in the 50\,ha plot is present for comparison.} \label{fig:est} \end{figure} Here we plot the expected values and also $\pm$ 2\,SE. <<>>= lines(subs3, rar3[1,], type='b') lines(subs3, rar3[1,]+2*rar3[2,], lty=3) lines(subs3, rar3[1,]-2*rar3[2,], lty=3) @ Next we hope to estimate the minimum total number of species (asymptotic richness) we might observe in the area around (and in) our 10 plots, if we can assume that the surrounding forest is homogeneous (it would probably be best to extrapolate only to the 50\,ha plot). First, we use an \emph{abundance-based \index{coverage estimator}coverage estimator}, \emph{\index{ACE}ACE}, that appears to give reasonable estimates \cite{Magurran2004}. We plot intervals, the expected values $\pm$ 2\,SE (Fig. \ref{fig:est}). <<>>= ace <- estimateR(N) segments(6000, ace['S.ACE']-2*ace['se.ACE'], 6000, ace['S.ACE']+2*ace['se.ACE'], lwd=3) text(6000, 150, "ACE estimate", srt=90, adj=c(1,.5) ) @ Next we use a frequency-based estimator, \index{Chao 2}Chao 2, where the data only need to be presence/absence, but for which we also need multiple sample plots. <<>>= chaoF <- specpool(bci) segments(6300, chaoF[1,'Chao'] - 2*chaoF[1,'Chao.SE'], 6300, chaoF[1,'Chao'] + 2*chaoF[1,'Chao.SE'], lwd=3, col='grey') text(6300, 150, "Chao2 estimate", srt=90, adj=c(1,.5) ) @ Last we add the observed number of tree species (over 10\,cm dbh) found in the entire 50\,ha plot. <<>>= points(6700, dim(BCI)[2], pch=19, cex=1.5) text(6700, 150, "Richness observed in 50 ha", srt=90, adj=c(1,.5) ) @ This shows us that the total richness estimators did not overestimate the total number of species within the extent of this relatively homogenous sample area (Fig. \ref{fig:est}). If we wanted to, we could then use any of these three estimators to compare the richness in this area to the richness of another area. \section{Distributions} In addition to plotting species in multidimensional space (Fig. \ref{fig:mds}), or estimating a measure of diversity or richness, we can also examine the \emph{distributions} of species abundances. Like any other vector of numbers, we can make a histogram of species abundances. As an example, here we make a histogram of tree densities, where each species has its own density (Fig. \ref{bci1}). This shows us what is patently true for nearly all ecological communities --- \emph{most species are rare}. \subsection{Log-normal distribution} Given general empirical patterns, that most species are rare, Frank Preston \cite{Preston1948,Preston1962} proposed that we describe communities using the logarithms of species abundances (Fig. \ref{fig:bci}).\footnote{Preston used base 2 logs to make his histogram bins, and his practice remains standard; we use the natural log.} This often reveals that a community can be described approximately with the normal distribution applied to the log-transformed data, or the \emph{\index{log-normal ditribution|see{species--abundance distribution}}\index{log-normal!species--abundance distribution}log-normal ditribution}. We can also display this as a \emph{\index{rank--abundance distribution}rank--abundance distribution} (Fig. \ref{bci3}). To do this, we assign the most abundant species as rank = 1, and the least abundant has rank = $R$, in a sample of $R$ species, and plot log-abundance \emph{vs.} rank. \index{rank--abundance distribution|see{species--abundance distribution}} May \cite{May1975} described several ways in which common processes may drive log-normal distributions, and cause them to be common in data sets. Most commonly cited is to note that log-normal distributions arise when each observation (i.e., each random variable) results from the product of independent factors. That is, if each species' density is determined by largely independent factors which act multiplicatively on each species, the resulting densities would be log-normally distributed. \begin{figure}[ht] \centering \subfloat[Raw data]{% \includegraphics[width=.325\textwidth]{bci1.pdf} \label{bci1}} \subfloat[Species--abundance dist.]{% \includegraphics[width=.325\textwidth]{bci2.pdf} \label{bci2}} \subfloat[Rank--abundance dist.]{% \includegraphics[width=.325\textwidth]{bci3.pdf} \label{bci3}} \caption{Three related types of distributions of tree species densities from Barro Colorado Island \cite{Condit2002}. \subref{bci1} Histogram of raw data, \subref{bci2} histogram of log-transformed data; typically referred to as the ``species--abundance distribution,'' accompanied here with the normal probability density function, \subref{bci3} the ``rank--abundance distribution,'' as typically presented with the log-transformed data, with the complement of the cumulative probability density function (1-pdf) \cite{May1975}. Normal distributions were applied using the mean and standard deviation from the log-transformed data, times the total number of species.\label{fig:bci}} \end{figure} \medskip \noindent \begin{boxedminipage}{\linewidth} {\footnotesize \paragraph{Log-normal abundance distributions (Fig. \ref{fig:bci})} We can plot tree species densities from Barro Colorado Island \cite{Condit2002}, which is aailable online, or in the \texttt{vegan} package. First we make the data available to us, then make a simple histogram. <>= data(BCI); N <- sort(colSums(BCI), decr=TRUE) hist(N, main=NULL) @ Next we make a \emph{species--abundance distribution}, which is merely a histogram of the log-abundances (classically, base 2 logs, but we use base $e$). In addition, we add the normal probability density function, getting the mean and standard deviation from the data, and plotting the expected number of species by multiplying the densities by the total number of species. <>= hist(log(N), xlab="Log Density Classes", main=NULL) m.spp <- mean(log(N)); sd.spp <- sd(log(N)); R <- length(N) curve(dnorm(x,m.spp,sd.spp)*R, 0, 8, add=T) @ Next we create the \emph{rank--abundance distribution}, which is just a plot of log-abundances \emph{vs.} ranks. <>= plot(log(N), type='b', ylim=c(0,8), main=NULL, xlab="Species Rank", ylab="Log Density") ranks.lognormal <- R * (1-pnorm( log(N), m.spp, sd.spp)) lines(ranks.lognormal, log(N)) @ We can think of the rank of species $i$ as the total number of species that are more abundant than species $i$. This is essentially the opposite (or complement) of the integral of the species--abundance distribution \cite{May1975}. That means that if we can describe the species abundance distribution with the normal density function, then 1-cumulative probability function is the rank. } \end{boxedminipage} \medskip \subsection{Other distributions} Well over a dozen other types of abundance distributions exist to describe abundance patterns, other than the log-normal \cite{Magurran2004}. They can all be represented as \emph{rank--abundance distributions}. The \emph{\index{geometric!species--abundance distribution}geometric distribution}\footnote{This probability mass function, $ \mathrm{P}_i = d\left(1-d\right)^{i-1}$, is the probability distribution of the number of attempts, $i$, needed for one success, if the independent probability of success on one trial is $d$.} (or pre-emption niche distribution) reflects a simple idea, where each species pre-empts a constant fraction of the remaining niche space \cite{May1975, Motomura1932}. For instance, if the first species uses 20\% of the niche space, the second species uses 20\% of the remaining 80\%, etc. The frequency of the $i$th most abundant species is \begin{equation} \label{eq:geomf} N_i = \frac{N_T}{C}d\left(1-d\right)^{i-1} \end{equation} where $d$ is the abundance of the most common species, and $C$ is just a constant to make $\sum N_i = N_T$, where $C=1-\left(1-d\right)^{S_T}$. Thus this describes the geometric rank--abundance distribution. The \emph{log-series distribution} \cite{Fisher1943} describes the frequency of species with $n$ individuals, \begin{equation} \label{eq:fisher} F\left(S_n\right) = \frac{\alpha x^n}{n} \end{equation} where $\alpha$ is a constant that represents diversity (greater $\alpha$ means greater diversity); the $\alpha$ for a diverse rainforest might be 30--100. The constant $x$ is a fitted, and it is always true that $0.9 < x < 1.0$ and $x$ increases toward 0.99 as $N/S \to 20$ \cite{Magurran2004}. $x$ can be estimated from $S/N=[(1-x)/x] \cdot [-\ln(1-x)]$. Note that this is not described as a rank--abundance distribution, but species abundances can nonetheless be plotted in that manner \cite{May1975}. The \index{log-series!species--abundance distribution}log-series rank--abundance distribution is a bit of a pain, relying on the standard exponential integral \cite{May1975}, $E_1(s) = \int_s^\infty exp(-t)/t\,dt$. Given a range of $N$, we calculate ranks as \begin{equation} \label{eq:1} F\left(N\right) = \alpha \int_s^\infty exp(-t)/t\,dt \end{equation} where we can let $t=1$ and $s = N \log \left(1 + \alpha/N_T\right)$. The log-series distribution has the interesting property that the total number of species in a sample of $N$ individuals would be $S_T = \alpha \log(1+N/\alpha)$. The parameter $\alpha$ is sometimes used as a measure of diversity. If your data are log-series distributed, then $\alpha$ is approximately the number of species for which you expect 1 individual, because $x\approx 1$. Two very general theories predict a log-series distribution, including neutral theory, and maximum entropy. Oddly, these two theories both predict a log-series distirbution, but make opposite assumptions about niches and individuals (see next section). MacArthur's \emph{\index{MacArthur's broken stick!species--abundance distribution}broken stick distribution} is a classic distribution that results in a very even distribution of species abundances \cite{MacArthur1957}. The number of individuals of each species $i$ is \begin{equation} \label{eq:brokenstick} N_i = \frac{N_T}{S_T}\sum_{n=i}^{S_T}\frac{1}{n} \end{equation} where $N_T$ and $S_T$ are the total number of individuals and species in the sample, respectively. MacArthur described this as resulting from the simultaneous breakage of a stick at random points along the stick. The resulting size fragments are the $N_i$ above. MacArthur's broken stick model is thus both a \emph{stochastic} and a \emph{deterministic} model. It has a simulation (stick breakage) that is the direct analogue of the deterministic analytical expression. Other similarly tactile stick-breaking distributions create a host of different rank--abundance patterns \cite{Magurran2004, Tokeshi:1999fv}. In particular, the stick can be broken\index{sequential broken stick!species--abundance distribution} \emph{sequentially}, first at one random point, then at a random point along one of two newly broken fragments, then at an additional point along any one of the \emph{three} broken fragments, \emph{etc.}, with $S_T-1$ breaks creating $S_T$ species. The critical difference between the models then becomes \emph{how each subsequent fragment is selected}. If the probability of selecting each fragment is related directly to its size, then this becomes identical to MacArthur's broken stick model. On the other hand, if each subsequent piece is selected randomly, regardless of its size, then this results in something very similar to the \index{log-normal!species--abundance distribution}log-normal distribution \cite{Sugihara1980, Tokeshi1990}. Other variations on fragment selection generate other patterns \cite{Tokeshi:1999fv}. \begin{figure}[ht] \centering \includegraphics[width=.67\linewidth]{other.pdf} \caption{A few common rank--abundance distributions, along with the \index{BCI}BCI data \cite{Condit2002}. The log-normal curve is fit to the data, and the broken stick distribution is always determined by the number of species. Here we let the geometric distribution be determined by the abundance of the most common species. The log-series was plotted so that it matched qualitatively the most abundant species.} \label{fig:other} \end{figure} \medskip \noindent \begin{boxedminipage}{\linewidth} {\footnotesize \paragraph{Generating other rank--abundance distributions)} We can illustrate the above rank--abundance distributions as they might relate to the BCI tree data (see previous code). We start with MacArthur's broken stick model. We use cumulative summation backwards to add all $1/n_i$, and then re-sort it by rank (cf. eq. \ref{eq:brokenstick}). <<>>= N <- sort(colSums(BCI), decr=TRUE) f1 <- sort(cumsum(1/(R:1)), decr=TRUE) Nt <- sum(N) NMac <- Nt*f1/sum(f1) @ Next, we create the geomtric rank--abundance distribution, where we let the BCI data tell us $d$, the density of the most abundant species; therefore we can multiply these by $N_T$ to get expected abundances. <<>>= d <- N[1]/Nt Ngeo.f <- d*(1-d)^(0:(R-1)) Ngeo <- Nt * Ngeo.f @ Last, we generate a log-series relevant to the BCI data. First, we use the \texttt{optimal.theta} function in the \texttt{untb} package to find a relevant value for Fisher's $\alpha$. (See more and $\theta$ and $\alpha$ below under neutral theory). <>= library(untb) alpha <- optimal.theta(N) @ To calculate the rank abundance distribution for the log-series, we first need a function for the ``standard exponential integral'' which we then integrate for each population size. <<>>= sei <- function( t=1 ) exp(-t)/t ### standard expo integral alpha <- optimal.theta(N) ranks.logseries <- sapply(N, function(x) { n <- x * log(1 + alpha/Nt ) f <- integrate(sei, n, Inf) fv <- f[["value"]] alpha * fv } ) @ } \end{boxedminipage} \medskip @ \medskip \noindent \begin{boxedminipage}{\linewidth} {\footnotesize \paragraph{Plotting other rank--abundance distributions (Fig. \ref{fig:other})} Now we can plot the BCI data, and all the distributions, which we generated above. Note that for the log-normal and the log-series, we calculated ranks, based on the species--abundance distributions, whereas in the standard form of the geometric and broken stick distributions, the expected abundances are calculated, in part, from the ranks. <>= plot(1:R, N, ylim=c(.1,2000), xlim=c(1,301), axes=FALSE, log='y') axis(2); axis(1, 1 + seq(0, 300, by=50));box() lines(1:R, NMac, lty=4) lines(1:R, Ngeo, lty=3) lines(ranks.logseries, N, lty=2) lines(ranks.lognormal, N, lty=1) legend("topright", c("Broken Stick", "Log-series", "Log-normal", "Geometric"), lty=c(4,2,1,3), bty='n' ) @ Note that we have not \emph{fit} the the log-series or geometric distributions to the data, but rather, placed them in for comparison. Properly fitting curves to distributions is a picky business \cite{Magurran2004, Oksanen:2008}, especially when it comes to species abundance distributions. } \end{boxedminipage} \medskip \subsection{Pattern \emph{vs.} process} \index{pattern vs. process@pattern \emph{vs.} process} \index{pattern vs. process@pattern \emph{vs.} process!species--abundance distribution} Note the resemblance between stick-breaking and niche evolution --- if we envision the whole stick as all of the available niche space, or energy, or limiting resources, then each fragment represents a portion of the total occupied by each species. Thus, various patterns of breakage act as models for niche partitioning and relative abundance patterns. Other biological and stochastic processes create specific distributions. For instance, completely random births, deaths, migration, and speciation will create the log-series distribution and the log-normal-like distributions (see \emph{neutral theory} below). We noted above that independent, multiplicatively interacting factors can create the log-normal distribution. On the other hand, all of these abundance distributions should probably be used primarily to describe \emph{patterns} of commonness, rarity, and not to infer anything about the processes creating the patterns. These graphical descriptions are merely attractive and transparent ways to illustrate the abundances of the species in your sample/site/experiment. The crux of the issue is that different \emph{processes} can cause the same abundance distribution, and so, sadly, \emph{we cannot usually infer the underlying processes from the patterns we observe}. That is, correlation is different than causation. Abundances of real species, in nature and in the lab, are the result of \emph{mechanistic processes}, including as those described in models of abundance distributions. However, we cannot say that a particular pattern was the result of a particular process, based on the pattern alone. Nonetheless, they are good summaries of ecological communities, and they show us, in a glance, a lot about \emph{diversity}. Nonetheless, interesting questions remain about the relations between patterns and processes. Many ecologists would argue that describing patterns and predicting them are the most important goals of ecology \cite{Peters1988}, while others argue that process is all important. However, both of these camps would agree about the fallacy of drawing conclusions about processes based on pattern --- nowhere in ecology has this fallacy been more prevalent than with abundance distributions \cite{Hairston:1991sf}. \section{Neutral Theory of Biodiversity and Biogeography} One model of abundance distributions is particularly important, and we elaborate on it here. It is referred to variously as the \emph{\index{unified neutral theory of biodiversity and biogeography|see{neutral theory}}unified neutral theory of biodiversity and biogeography} \cite{Hubbell2001}, or often ``neutral theory'' for short. Neutral theory is important because it does much more than most other models of abundance distributions. It is a testable theory that makes quantitative predictions across several levels of organizations, for both evolutionary and ecological processes. Just as evolutionary biology has neutral theories of gene and allele frequencies, ecology has neutral theories of population and community structure and dynamics \cite{Bell2000, Caswell1976, Hubbell2001}. Neutral communities of species are computationally and mathematically related and often identical to models of genes and alleles \cite{Ewens:2004zp} (Table \ref{tab:jargon}). Thus, lessons you have learned about genetic drift often apply to neutral models of communities. Indeed, neutral ecological dynamics at the local scale are often referred to as \emph{ecological drift}, and populations change via demographic stochasticity.\footnote{Some have called neutral theory a \emph{null model}, but others disagree \cite{Gotelli:2006bv}, describing distinctions between dynamical, process-based neutral models with fitted parameters, \emph{vs.} static data-based null models \cite{Gotelli1996}. Both can be used as benchmarks against which to measure other phenomena. Under those circumstances, I suppose they might both be null hypotheses.} Stephen Hubbell proposed his ``unified neutral theory of biodiversity and biogeography'' (hereafter NTBB, \cite{Hubbell2001}) as a null model of the dynamics of individuals, to help explain relative abundances of tropical trees. Hubbell describes it as a direct descendant of MacArthur and Wilson's theory of island biogeography \cite{MacArthur1963, MacArthur1967} (see below, species--area relations). Hubbell proposed it both as a null hypothesis and also --- and this is the controversial part --- as a model of community dynamics that closely approximates reality. <>= library(untb) par(mar=c(1,1,1,1)) display.untb(start=rep(10, 900), prob=0.1, gens=10000) #box() polygon(c(6.5,13.5,13.5,6.5), y=c(6.5,6.5,13.5,13.5), lwd=6) <>= out <- untb(rep(1:10, 90), prob=0, gens=1000, D=9, keep=TRUE) matplot(1:1000,species.table(out), type='l', lwd=1.5, ylab="Population Density", xlab="Generation") @ \begin{figure}[ht] \centering \includegraphics[width=.6\linewidth]{metacomm} \caption{A cartoon of a local community of forest canopy trees (small box) nested inside part of the metacommunity of a tropical forest. The true metaccommunity would extend far beyond the boundaries of this figure to include the true potential source of propagules. Shades of grey indicate different species. The local community is a sample of the larger community (such as the 50\,ha forest dynamics plot on BCI) and receives migrants from the metacommunity. Mutation gives rise to new species in the metacommunity. For a particular local community, such as a 50\,ha plot on an island in the Panama canal, the metacommunity will include not only the surrounding forest on the island, but also Panama, and perhaps much of the neotropics \cite{Jabot:2009mi}.} \label{fig:meta} \end{figure} The relevant ``world'' of the NTBB is a \index{metacommunity|see{neutral theory}} \index{metacommunity!neutral theory}metacommunity (Fig. \ref{fig:meta}), that is, a collection of similar local communities connected by dispersal \cite{Leibold:2004fe}.\footnote{The metacommunity concept is quite general, and a neutral metacommunity is but one caricature \cite{Leibold:2004fe}.} The metacommunity is populated entirely by individuals that are functionally identical. The NTBB is a theory of the \emph{dynamics of individuals}, modeling individual births, deaths, migration and mutation. It assumes that within a guild, such as late successional tropical trees, species are essentially neutral with respect to their fitness, that is, they exhibit \emph{\index{fitness equivalence|see{neutral theory}}\index{fitness equivalence!neutral theory}fitness equivalence}. This means that the probabilities of birth, death, mutation and migration are \emph{identical for all individuals}. Therefore, changes in population sizes occur via \emph{random walks}, that is, via stochastic increases and decreases with time step (Fig. \ref{fig:drift}). Random walks do not imply an absence of competition or other indirect enemy mediated negative density dependence. Rather, competition is thought to be diffuse, and equal among individuals. We discuss details of this in the context of simulation. Negative density dependence arises either through a specific constraint on the total number of individuals in a community \cite{Hubbell2001}, or as traits of individuals related to the probabilities of births, deaths, and speciation \cite{Volkov:2005rt}. \begin{figure}[ht] \centering \includegraphics[width=.6\linewidth]{driftfig} \caption{Neutral ecological drift. Here we start with 10 species, each with 90 individuals, and let their abundances undergo random walks within a finite local community, with no immigration. Here, one generation is equal to nine deaths and nine births. Note the slow decline in unevennes --- after 1000 deaths, no species has become extinct.} \label{fig:drift} \end{figure} A basic paradox of the NTBB, is that in the absence of migration or mutation, diversity gradually declines to zero, or monodominance. A random walk due to fitness equivalence will eventually result in the loss of all species except one.\footnote{This is identical to allele frequency in \index{parallels with genetic drift!neutral theory}genetic drift.} However, the loss of diversity in any single local community is predicted to be very, very slow, and is countered by immigration and speciation (we discuss more details below). Thus, species do not increase deterministically when rare --- this makes the concept of coexistence different than the stable coexistance criteria discussed in previous chapters. Coexistance here \emph{is not stable} but rather only a stochastic event with a limited time horizon which is balanced by the emergence of new species. If all communities are thought to undergo random walks toward monodominance, how is diversity maintained in any particular region? Two factors maintain species in any local community. First, immigration into the local community from the metacommunity can bring in new species. Even though each local community is undergoing a random walk toward monodominance, each local community may become dominated by any one of the species in the pool because all species have equal fitness. Thus separate local communities are predicted to become dominated by different species, and these differences among local communities help maintain diversity in the metacommunity landscape.\footnote{This is the same as how genetic drift operates across subpopulations connected by low rates of gene exchange.} Second, at longer time scales and larger spatial scales, speciation (i.e., mutation and lineage-splitting) within the entire metacommunity maintains biodiversity. Mutation and the consequent speciation provide the ultimate source of variation. Random walks toward extinction in large communities are \emph{so} lengthy that the extinctions are balanced by speciation. Introducing new symbols and jargon for ecological neutral theory, we state that the diversity of the metacommunity, $\theta$, is a function of the number of individuals in the metacommunity, \index{Jm@$J_M$}\index{Jm@$J_M$!neutral theory}$J_M$, and the per capita rate at which new species arise (via mutation) \index{$\nu$|see{neutral theory}}\index{$\nu$}$\nu$ ($\theta = 2J_M \nu$; Table \ref{tab:jargon}). A local community undergoes \index{drift|see{neutral theory}}ecological drift; drift causes the slow loss of diversity, which is balanced by a per capita ($J_L$) immigration rate $m$. \begin{table} \centering \caption{Comparison of properties and jargon used in ecological and population genetic neutral theory (after Alonso et al. \cite{Alonso:2006pi}). Here $x$ is a continuous variable for relative abundance of a species or allele ($0 \le x \ge 1$, $ x = n/J$). Note this is potentially confused with Fisher's log-series $\langle \phi_n \rangle = \theta x^n /n$, which is a discrete species abundance distribution in terms of numbers of individuals, $n$ (not relative abundance), and where $x=b/d$.} \label{tab:jargon} \begin{tabular}[c]{p{.36\linewidth}p{.315\linewidth}p{.315\linewidth}} \hline \hline \textbf{Property} & \textbf{Ecology} & \textbf{Population Genetics} \\ \hline Entire System (size) & Metacommunity ($J_M$) & Population ($N$)\\ Nested subsystem (size) & Local community ($J_L$) & Subpop. or Deme ($N$)\\ Smallest neutral system unit & Individual organism & Individual gene \\ Diversity unit & Species & Allele \\ Stochastic process & Ecological drift & Genetic drift \\ Generator of diversity (rate symbol) & Speciation ($\nu$) & Mutation ($\mu$) \\ Fundamental diversity number & $\theta \approx 2J_M\nu$ & $\theta \approx 4N\mu$\\ Fundamental dispersal number & $I \approx 2J_Lm$ & $\theta \approx 4Nm$ \\ Relative abundance distribution ($\Phi \left(x\right)$) & $\frac{\theta}{x}\left(1-x\right)^{\theta - 1}$ & $\frac{\theta}{x}\left(1-x\right)^{\theta - 1}$\\ Time to common ancestor & $\frac{-J_M x}{1-x}\log x$ &$\frac{-N x}{1-x}\log x$ \\ \end{tabular} \end{table} It turns out that the abundance distribution of an infinitely large metacommunity is Fisher's \index{log-series!species--abundance distribution}log-series distribution, and that \index{$\theta$, diversity|see{neutral theory}}$\theta$ of neutral theory is $\alpha$ of Fisher's log-series \cite{Alonso2004, Leigh:1999qr, Volkov:2003ly}. However, in any one local community, random walks of rare species are likely to include zero, and thus become extinct in the local community by chance alone. This causes a deviation from Fisher's log-series in any \emph{local community} by reducing the number of the the rarest species below that predicted by Fisher's log-series. In particular, it tends to create a log-normal--like distribution, much as we often see in real data (Fig. \ref{fig:other}). These theoretical findings and their match with observed data are thus consistent with the hypothesis that communities may be governed, in some substantive part, by neutral drift and migration. Both theory and empirical data show that species which coexist may be more similar than predicted by chance alone \cite{Leibold:2006ey}, and that similarity (i.e., fitness equality) among species helps maintain higher diversity than would otherwise be possible \cite{Chesson2000}. Chesson makes an important distinction between \emph{\index{stabilizing vs. equalizing mechanisms!neutral theory}stabilizing mechanisms}, which create attractors, and \emph{equalizing mechanisms}, which reduce differences among species, slow competitive exclusion and facilitate stabilization \cite{Chesson2000}. The NTBB spurred tremendous debate about the roles of chance and determinism, of \index{dispersal-assembly and niche-assembly!neutral theory}dispersal-assembly and niche-assembly, of evolutionary processes in ecology, and how we go about ``doing community ecology'' (see, e.g., articles in \emph{Functional Ecology}, 19(1), 2005; \emph{Ecology}, 87(6), 2006). This theory in its narrowest sense has been falsified with a few field experiments and observation studies \cite{Clark2003c, Wootton2005}. However, the \emph{degree} to which stochasticity and dispersal \emph{versus} niche partitioning structure communities remains generally unknown. Continued refinements and elaboration (e.g., \cite{Alonso2004, Etienne:2007fe, Green:2007nx, Jabot:2009mi, Plotkin:2002eu, Volkov:2007sy}) seem promising, continuing to intrigue scientists with different perspectives on natural communities \cite{Jabot:2009mi, Latimer:2005ft}. Even if communities turn out to be completely non-neutral, NTBB provides a baseline for community dynamics and patterns that has increased the rigor of evidence required to demonstrate mechanisms controlling coexistence and diversity. As Alonso \emph{et al.} state, ``\ldots good theory has more predictions per free parameter than bad theory. By this yeardstick, neutral theory fares fairly well'' \cite{Alonso:2006pi}. \subsection{Different flavors of neutral communities} Neutral dynamics in a local community can be simulated in slightly different ways, but they are all envisioned some type of \emph{random walk}. A random walk occurs when individuals reproduce or die at random, with the consequence that each population increases or decreases by chance. The simplest version of a \index{random walk}\index{random walk!neutral theory}random walk assumes that births and deaths and consequent increases and decreases in population size are equally likely and equal in magnitude. A problem with this type of random walk is that a community can increase in size (number of individuals) without upper bound, or can disappear entirely, by chance. We know this doesn't happen, but it is a critically important first step in conceptualizing a neutral dynamic \cite{Caswell1976}. Hubbell added another level of biological reality by fixing the total number of individuals in a local community, $J_L$, as constant. When an individual dies, it is replaced with another individual, thus keeping the population size constant. Replacements come from within the local community with probability $1-m$, and replacements come from the greater metacommunity with probability $m$. The dynamics of the metacommunity are so slow compared to the local community that we can safely pretend that it is fixed, unchanging. Volkov \emph{et al.} \cite{Volkov:2003ly} took a different approach by assuming that each species undergoes independent \emph{\index{random walks, biased !neutral theory}biased random walks}. We imagine that each species undergoes its own completely independent random walk, \emph{as if} it is not interacting with any other species. The key is that the birth rate, $b$, is slightly \emph{less} than the death rate, $d$ --- this bias toward death gives us the name \emph{biased random walk}. In a deterministic model with no immigration or speciation, this would result in a slow population decline to zero. In a stochastic model, however, some populations will increase in abundance by chance alone. Slow random walks toward extinctions are balanced by speciation in the metacommunity (with probability~$\nu$). If species all undergo independent biased random walks, does this mean species don't compete and otherwise struggle for existance? No. The \emph{reason} that $b < d$ is precisely because all species struggle for existance, and only those with sufficiently high fitness, \emph{and which are lucky}, survive. Neutral theory predicts that it is these species that we observe in nature --- those which are lucky, and also have sufficiently high fitness. In the metacommunity, the average number of species, $\langle \phi_n^M \rangle$, with population size $n$ is \begin{equation} \label{eq:meta} \langle \phi_n^M \rangle = \theta \frac{x^n}{x} \end{equation} where $x = b/d$, and $\theta = 2J_M\nu$ \cite{Volkov:2003ly}. The $M$ superscript refers to the metacommunity, and the angle brackets indicate merely that this is the average. Here $b/d$ is barely less than one, because it a biased random walk which is then offset by speciation, $\nu$. Now we see that this is exactly \index{Fisher's log-series|see{species--abundance distribution, log-series}}Fisher's log-series distribution (eq. \ref{eq:fisher}), where that $x=b/d$ and $\theta = \alpha$. Volkov \emph{et al.} thus show that in a neutral metacommunity, $x$ has a biological interpretation. The expected size of the entire metacommunity is simply the sum of all of the average species' $n$ \cite{Volkov:2003ly}. \begin{equation} \label{eq:Jm} J_M = \sum_{n=1}^{\infty} n \langle \phi_n^M \rangle = \theta \frac{x}{1-x} \end{equation} Thus the size of the metacommunity is an emergent property of the dynamics, rather than an external constraint. To my mind, it seems that the number of individuals in the metacommunity must result from responses of individuals \emph{to} their external environment. Volkov \emph{et al.} went on to derive expressions for births, deaths, and average relative abundances in the local community \cite{McKane2000, Volkov:2003ly}. Given that each individual has an equal probability of dying and reproducing, and that replacements can also come from the metacommunity with a probability proportional to their abundance in the metacommunity, one can specify rules for populations in local communities of a fixed size. These are the probability of increase, $b_{n,k}$, or decrease, $d_{n,k}$, in a species, $k$, of population size $n$. \begin{align} b_{n,k} & = \left(1-m\right)\frac{n}{J_L}\frac{J_L-n}{J_L-1} + m\frac{\mu_k}{J_M}\left(1-\frac{n}{J_L}\right) \label{eq:b}\\ d_{n,k}&=(1-m)\frac{n}{J_L}\frac{J_L-n}{J_L-1} + m\left(1-\frac{mu_k}{J_M}\right)\frac{n}{J_L} \label{eq:d} \end{align} These expressions are the sum of two joint probabilities, each of which is comprised of several independent events. These events are immigration, and birth and death of individuals of different species. Here we describe these probabilities. For a population of size $n$ of species $k$, we can indicate per capita, per death probabilities including \begin{itemize} \item $m$, the probability that a replacement is immigrant from the metacommunity, and $1-m$, the probability that the replacement is from the local community. \item $n/J_L$, the probability that an individual selected randomly from the local community belongs to species $k$, and $1 - n/J_L$ or $(J-n)/(J_L)$, the probability that an individual selected randomly from the local community belongs to any species other than $k$. \item $(J-n)/(J_L-1)$, the conditional probability that, given that an individual of species $k$ has already been drawn from the population, an individual selected randomly from the local community belongs to any species other than to species $k$. \item $\mu_k/J_M$, the probability that an individual randomly selected from the metacommunity belongs to species $k$, and $1 - n/J_M$, the probability that an individual randomly selected from the metacommunity belongs to any species other than $k$. \end{itemize} Each of these probabilities is the probability of some unspecified event --- that event might be birth, death, or immigration. Before we translate eqs. \ref{eq:b}, \ref{eq:d} literally, we note that $b$ and $d$ each have two terms. The first term is for dynamics related to the local community, which happen with probability $1-m$. The second is related to immigration from the metacommunity which occurs with probability $m$. Consider also that if a death is replaced by a birth of the same species, or a birth is countered by a death of the same species, they cancel each other out, as if nothing ever happened. Therefore each term requires a probability related to species $k$ and to non-$k$. Eq. \ref{eq:b}, $b_{n,k}$, is the probability that an individual will be added to the population of species $k$. The first term is the joint probability that an addition to the population comes from within the local community ($1-m$) \emph{and} the birth comes from species $k$ ($n/J_L$) \emph{and} there is a death of an individual of any other species ($(J_L-n)/(J_L-1)$).\footnote{The denominator of the death probability is $J_L-1$ instead of $J_L$ because we have already selected the first individual who will do the reproduction, so the total number of remaining individuals is $J_L-1$ rather than $J_L$; the number of non-$k$ individuals remains $J_L-n$} The second term is the joint probability that the addition to the population comes from the metacommunity via immigration ($m$) \emph{and} that the immigrant is of species $k$ ($\mu_k/J_M$) \emph{and} is not accompanied by a death of an individual of its species ($n/J_L$). An individual may be subtracted from the population following similar logic. Eq. \ref{eq:d}, $d_{n.k}$, is the probability that a death will remove an individual from the population of species $k$. The first term is the joint probability that the death occurs in species $k$ ($n/J_L$) and the replacement comes from the local community ($1-m$) and is some species other than $k$ ($(J_L-n)/(J_L-1)$). The second term is the joint probability that the death occurs in species $k$ ($n/J_L$), and that it is replaced by an immigrant ($m$) \emph{and} the immigrant is any species in the metacommunity other than $k$ ($1-\mu_k/J_M$). \subsection{Investigating neutral communities} Here we exxplore netural communities using the \index{untb@\texttt{untb}}\texttt{untb} package, which contains a variety of functions for teaching and research on neutral theory. \subsubsection{Pure drift} After loading the package and setting graphical parameters, let's run a simulation of drift. Recall that drift results in the slow extinction of all but one species (Fig. \ref{fig:drift}). We start with a local community with 20 species, each with 25 individuals\footnote{This happens to be the approximate tree density (450 trees\,ha$^{-1}$, for trees $>10$\,cm DBH) on BCI.} The simulation then runs for 1000 generations (where 9/900 individuals die per generation).\footnote{Note that \texttt{display.untb} is great for pretty pictures, whereas \texttt{untb} is better for more serious simulations.} <<>>= library(untb) a <- untb(start=rep(1:20, 25), prob=0, gens=2500, D=9, keep=TRUE) @ We \texttt{keep}, in a matrix, all 450 trees from each of the 1000 time steps so that we can investigate the properties of the community. The output is a matrix where each element is an integer whose value represents a species ID. Rows are time steps and we can think of columns as spots on the ground occupied by trees. Thus a particular spot of ground may be occupied by an individual of one species for a long time, and suddenly switch identity, because it dies and is replaced by an individual of another species. Thus the community always has 450 individuals (columns), but the identities of those 450 change through time, according to the rules laid out in eqs. \ref{eq:b}, \ref{eq:d}. Each different species is represented by a different integer; here we show the identitites of ten individuals (columns) for generations 901--3. <<>>= ( a2 <- a[901:903, 1:10 ] ) @ Thus, in generation 901, tree no. 1 is species \Sexpr{a2[1,1]} and in generation 902, tree no. 3 is species \Sexpr{a2[2,3]}. We can make pretty pictures of the communities at time steps 1, 100, and 2000, by having a single point for each tree, and coding species identity by shades of grey.\footnote{We could use a nice color palette, \texttt{hcl}, based on hue, chroma, and luminance, for instance \texttt{hcl(a[i,]*30+50)}} <<>>= times <- c(1,50,2500) sppcolors <- sapply(times, function(i) grey((a[i,]-1)/20) ) @ This function applies to the community, at each time step, the \texttt{grey} function. Recall that species identity is an integer; we use that integer to characterize each species' shade of grey. Next we create the three graphs at three time points, with the appropriate data, colors, and titles. <>= layout(matrix(1:3, nr=1)) par(mar=c(1,1,3,1)) for(j in 1:3){ plot(c(1,20), c(1,25), type="n", axes=FALSE) points(rep(1:20,25), rep(1:25,each=20), pch=19, cex=2, col=sppcolors[,j]) title(paste("Time = ", times[j], sep="")) } <>= layout(matrix(1:9, nr=3, byrow=TRUE)) par(mar=c(1,1,3,1)) for(j in 1:3){ plot(c(1,20), c(1,25), type="n", axes=FALSE) points(rep(1:20,25), rep(1:25,each=20), pch=19, cex=1.5, col=sppcolors[,j]) title(paste("Time = ", times[j], sep="")) } par(mar=c(5,4,1,1)) for(i in times){ plot(as.count(a[i,]), ylim=c(1,max(as.count(a[2000,]))), xlim=c(0,20)) } par(mar=c(6,4,1,3)) out <- lapply(times, function(i) preston(a[i,]) ) bins <- matrix(0, nrow=3, ncol=length(out[[3]]) ) for(i in 1:3) bins[i,1:length(out[[i]])] <- out[[i]] bins colnames(bins) <- names(preston(a[times[3],])) for(i in 1:3){ par(las=2) barplot(bins[i,], ylim=c(0,20), xlim=c(0,8), ylab="No. of Species", xlab="Abundance Category" ) } @ \begin{figure} \centering \includegraphics[width=\linewidth]{comms} \caption{Three snapshots of one community, drifting through time. Shades of grey represent different species. Second row contains rank abundance distributions; third row contains species abundance distributions. Drift results in the slow loss of diversity.} \label{fig:comms} \end{figure} From these graphs (Fig. \ref{fig:comms}), we see that, indeed, the species identity of the 450 trees changes through time. Remember that this is \emph{not spatially explicit} --- we are not stating that individual $A$ is next, or far away from, individual $B$. Rather, this is a spatially implicit representation --- all individuals are characterized by the same probaibllities. Next let's graph the rank abundance distribution of these communities (Fig. \ref{fig:comms}). For each time point of interest, we first coerce each ``ecosystem'' into a \texttt{count} object, and then plot it. Note that we plot them all on a common scale to facilitate comparison. <<>>= layout(matrix(1:3, nr=1)) for(i in times){ plot(as.count(a[i,]), ylim=c(1,max(as.count(a[times[3],]))), xlim=c(0,20)) title(paste("Time = ", i, sep="")) } @ Next we create species abundance distributions, which are histograms of species' abundances (Fig. \ref{fig:comms}). If we want to plot them on a common scale, it takes a tad bit more effort. We first create a matrix of zeroes, with enough columns for the last community, use \texttt{preston} to create the counts of species whose abundances fall into logarithmic bins, plug those into the matrix, and label the matrix columns with the logarithmic bins. <<>>= out <- lapply(times, function(i) preston(a[i,]) ) bins <- matrix(0, nrow=3, ncol=length(out[[3]]) ) for(i in 1:3) bins[i,1:length(out[[i]])] <- out[[i]] bins colnames(bins) <- names(preston(a[times[3],])) @ Finally, we plot the species--abundance distributions. <<>>= layout(matrix(1:3, nr=1)) for(i in 1:3){ par(las=2) barplot(bins[i,], ylim=c(0,20), xlim=c(0,8), ylab="No. of Species", xlab="Abundance Category" ) } @ Bottom line: \emph{drift causes the slow loss of species from local communities} (Fig. \ref{fig:comms}). What is not illustrated here is that without dispersal, drift will cause different species to become abundant in different places because the variation is random. In that way, drift maintains diversity at large scales, in the metacommunity. Last, low rates of dispersal among local communities maintains some diversity in local communities without changing the entire metacommunity into a single large local community. Thus dispersal limitiation, but not its absence, maintains diversity. Next, let's examine the dynamics through time. We will plot individual species trajectories (Fig. \ref{spptable}). <>= sppmat <- species.table(a) matplot(1:times[3], sppmat[1:times[3],], type='l', ylab="Population Density") @ The trajectories all start at the same abundance, but they need not have. The trajectories would still have the same drifting quality. \begin{figure}[ht] \centering \subfloat[20 populations]{\includegraphics[width=.45\linewidth]{spptable} \label{spptable}} \subfloat[Average variability]{\includegraphics[width=.45\linewidth]{vartime} \label{vartime}} \caption{Dynamics and average variation within populations. In random walks, average variation (measured with the coefficient of variation) increases with time. } \label{fig:walks} \end{figure} For random walks in general, the observed variance and \index{coefficient of variation}coefficient of variation ($CV=\hat{\sigma}/\bar{x}$) of a population will grow over time \cite{Clark2003c}. Here we calculate the average population $CV$ of cumulative observations (note the convenient use of nested \texttt{(s)apply} functions). Let's calculate the $CV$'s for every tenth time step. <<>>= cvtimes <- seq(2, 2500, by=10) CV.within <- sapply(cvtimes, function(i) {# Use sapply with i as a row counter ## For each row i, use apply to caluculate the variance for each ## population from 1:i cvs <- apply(sppmat[1:i,],2, function(x) sd(x)/mean(x) ) mean(cvs) # get the mean among populations (for each i). } ) #popvar.among <- apply(sppmat, 1, var) @ Now plot the average $CV$ through time. The average observed $CV$ should increase (Fig. \ref{vartime}). <>= plot( cvtimes, CV.within, type='l' ) @ This shows us that the populations never approach an equilibrium, but wander aimlessly. \subsubsection{Real data} Last, we examine a \index{BCI}BCI data set \cite{Condit2002}. We load the data (data from 50 1\,ha plots $\times$ 225 species, from the \texttt{vegan} package), and sum species abundances to get each species total for the entire 50\,h plot (Fig. \ref{bcipreston}). <<>>= library(vegan) data(BCI) n <- colSums(BCI) par(las=2, mar=c(7,4,1,1)) xs <- plot(preston(rep(1:225,n), n=12, original=TRUE)) @ We would like to estimate $\theta$ and $m$ from these data, but that requires specialized software for any data set with a realistic number of individuals. Specialized software would provide maximum likelihood estimates (in a reasonable amount of computing time) for $m$ and $\theta$ for large communities \cite{Etienne2005, Hankin:2007zm, Jabot:2009mi}. The BCI data have been used repeatedly, so we rely on estimates from the literature ($\theta \approx 48$, $m \approx 0.1$) \cite{Etienne2005, Volkov:2003ly}. We use the approach of Volkov et al. \cite{Volkov:2003ly} to generate expected species abundances (Fig. \ref{bcipreston}). <<>>= v1 <- volkov(sum(n), c(48, 0.1), bins=TRUE) points(xs, v1[1:12], type='b') @ More recently, Jabot and Chave \cite{Jabot:2009mi} arrived at estimates that differed from previous estimates by orders of magnitude (Fig. \ref{bciJabot}). Their novel approach estimated $\theta$ and $m$ were derived from \emph{both} species abundance data and from phylogenetic data ($\theta \approx 571,\,724$, $m \approx 0.002$). This is possible because neutral theory makes a rich array of predictions, based on both ecological and evolutionary processes. Their data were only a little bit different (due to a later census), but their analyses revealed radically different estimates, with a much greater diversity and larger metacommunity (greater $\theta$), and much lower immigration rates (smaller $m$). Here we derive expected species abundance distributions. The first is based soley on census data, and is similar to previous expections (Fig. \ref{bciJabot}). <<>>= v2 <- volkov(sum(n), c(48, 0.14), bins=TRUE) xs <- plot(preston(rep(1:225,n), n=12, original=TRUE), col=0, border=0) axis(1, at=xs, labels=FALSE) points(xs, v2[1:12], type='b', lty=2, col=2) @ However, when they also included phylogenetic data, they found very different expected species abundance distributions (Fig. \ref{bciJabot}). <<>>= v4 <- volkov(sum(n), c(571, 0.002), bins=TRUE) points(xs, v4[1:12], type='b', lty=4, col=2) @ <>= par(las=2, mar=c(5,4,1,1)) xs <- plot(preston(rep(1:225,n), n=12, original=TRUE)) points(xs, v1[1:12], type='b') <>= par(las=2, mar=c(5,4,1,1)) xs <- plot(preston(rep(1:225,n), n=12, original=TRUE), col=0, border=0) axis(1, at=xs, labels=FALSE) points(xs, v2[1:12], type='b') points(xs, v4[1:12], type='b', lty=2) @ If nothing else, these illustrate the effects of increasing $\theta$ and reducing $m$. \begin{figure}[ht] \centering \subfloat[Original estimates]{\includegraphics[width=.48\linewidth]{bcipreston.pdf} \label{bcipreston}} \subfloat[Jabot and Chave estimates]{\includegraphics[width=.48\linewidth]{bciJabot.pdf} \label{bciJabot}} \caption{Species abundance distributions for BCI trees. \subref{bcipreston} Data for histogram from \cite{Condit2002}, overlain with expected abundances with $\theta$ and $m$ values fitted to the data \cite{Etienne2005, Volkov:2003ly}. \subref{bciJabot} Jabot and Chave found that when they used only species abundances (as did previous investigators) their pattern was similar to previous findings (solid line). However, adding phylogenetic information led to very different expectations (dashed line).} \label{fig:Jabot} \end{figure} \subsection{Symmetry and the rare species advantage} An important advance in neutral theory is the quantitfication of a \emph{symmetric rare species advantage}. The \index{symmetry|see{neutral theory}}\index{symmetry!neutral theory}symmetry hypothesis posits that all species have a \emph{symmetrical} \emph{rare species advantage} \cite{Volkov:2005rt, Volkov:2007sy}. That is, all species increase when rare to the same degree (equal negative density dependence). In a strict sense, all individuals remain the same in that their birth and death probabilities change with population size in the same manner for all species. This obviously reduces the chance of random walks to extinction, but is nonetheless the same among all species. Estimation of the magnitude of the rare species advantage is interesting addition to stepwise increasing complexity. To sum up: it is safe to say that neutral theory has already made our thinking about community structure and dynamics more sophisticated and subtle, by extending island biogeography to individuals. The theory is providing quantitative, pattern-generating models, that are analogous to null hypotheses. With the advent of specialized software, theory is now becoming more useful in our analysis of data \cite{Etienne2005, Hankin:2007zm, Jabot:2009mi}. @ \section{\index{diversity partitioning}Diversity Partitioning} \index{partitioning|see{diversity partitioning}} We frequently refer to biodiversity (i.e., richness, Simpson's, and Shannon-Wiener diversity) at different spatial scales as \index{$\alpha$-diversity}$\alpha$, \index{$\beta$-diversity}$\beta$, and \index{$\gamma$-diversity}$\gamma$ diversity (Fig. \ref{fig:moths}). \begin{itemize} \item Alpha diversity, $\alpha$, is the diversity of a point location or of a single sample. \item Beta diversity, $\beta$, is the diversity due to multiple localities; $\beta$ diversity is sometimes thought of as \emph{turnover} in species composition among sites, or alternatively as the number of species in a region that are \emph{not} observed in a sample. \item Gamma diversity, $\gamma$, is the diversity of a region, or at least the diversity of all the species in a set of samples collected over a large area (with large extent relatve to a single sample). \end{itemize} Diversity across spatial scales can be further be \emph{partitioned} in one of two ways, either using \emph{additive} or \emph{multiplicative} partitioning. \emph{Additive partitioning} \cite{Crist:2006lr,Crist:2003ym,Lande1996} is represented as \begin{equation} \label{eq:add} \bar{\alpha} + \beta = \gamma \end{equation} where $\bar{\alpha}$ is the average diversity of a sample, $\gamma$ is typically the diversity of the pooled samples, and $\beta$ is found by difference ($\beta = \gamma - \bar{\alpha}$). We can think of $\beta$ as the \emph{average} number of species \emph{not} found in a sample, but which we know to be in the region. \index{additive partitioning|see{diversity partitioning}}Additive partitioning allows direct comparison of average richness among samples at any hierarchical level of organization because all three measures of diversity ($\alpha$, $\beta$, and $\gamma$) are \emph{expressed in the same units}. This makes it analogous to partitioning variance in ANOVA. This is not the case for multiplicative partitioning diversity. Partitioning can also be \emph{multiplicative} \cite{Whittaker:1960it}, \begin{equation} \label{eq:mult} \bar{\alpha}\beta = \gamma \end{equation} where $\beta$ is a conversion factor that describes the relative change in species composition among samples. \index{multiplicative partitioning|see{diversity partitioning}}Sometimes this type of $\beta$ diversity is thought of as the number of different community types in a set of samples. However, one must use this interpretation with great caution, as either meaning of $\beta$ diversity depends completely on the sizes or extent of the samples used for $\alpha$ diversity. << echo=false, results=hide>>= library(maps) library(ellipse) N <- 1000 crd <- map(, xlim=c(12,42), ylim=c(-5,5)) coords <- matrix(NA, nr=N, ncol=2) for(i in 1:N) coords[i,] <- { xr <- range(crd$x[!is.na(crd$x)]); yr <- range(crd$y[!is.na(crd$y)]) x <- runif(1, xr[1], xr[2]); y <- runif(1,yr[1], yr[2]) c(x,y) } ### Pick letters relative to latitude a <- (coords[,1] + .1 - min(coords[,1] ) )*2 b <- max(a) - a +.1 inds <- sapply(1:N, function(i) { sample(1:26, 1, replace=T, prob=dbeta(seq(0.001,.999,length=26), a[i], b[i])) }) inds <- sapply(1:N, function(i) { sample(LETTERS, 1, replace=T, prob=dbeta(seq(0.001,.999,length=26), a[i], b[i])) }) rws <- rep(c(-2,2), each=6) cls <- rep(seq(15,40,5), 2) x1s <- coords[1,] - rws[1] pts <- ellipse(0, centre=c(0,0), level=.5) sqrt(apply(pts^2,1,sum)) #quartz(height=7, width=7) <>= library(ellipse) #install.packages("maps", dep=TRUE) library(maps) map.text('world', xlim=c(12,42), ylim=c(-5,5)) map.axes() for(i in 1:12) { pts <- ellipse(0, centre=c(cls[i], rws[i]), level=.5, type="l" ) lines(pts) } points(coords, pch=inds, cex = .5) @ \begin{figure}[ht] \centering \includegraphics[width=.8\linewidth]{moths} \caption{Hierarchical sampling of \index{moths}moth species richness in forest patches in Indiana and Ohio, USA \cite{Summerville:2004ym}. $\alpha$-diversity is the diversity of a single site (richness indicated by numbers). $\gamma$-diversity is the total number of species found in any of the samples (here $\gamma=230$\,spp.). Additive $\beta$-diversity is the difference, $\gamma - \bar{\alpha}$, or the average number of species \emph{not} observed in a single sample. Diversity partitioning can be done at two levels, sites within \index{ecoregions}ecoregions and ecoregions within the geographic region (see example in text for details). } \label{fig:moths} \end{figure} Let us examine the the limits of $\beta$ diversity in extremely small and extremely large samples. Imagine that our sample units each contain, on average, one individual (and therefore one species) and we have 10$^6$ samples. If richness is our measure of diversity, then $\bar{\alpha}=1$. Now imagine that in all of our samples we find a total of 100 species, or $\gamma=100$. Our additive and multiplicative partitions would then be $\beta_{A}=99$, and $\beta_{M}=100$, respectively. If the size of the sample unit increases, each sample will include more and more individuals and therefore more of the diversity, and by definition, $\beta$ will decline. If each sample gets large enough, then each sample will capture more and more of the species until a sample gets so large that it includes all of the species (i.e., $\bar{\alpha} \to \gamma$). At this point, $\beta_A \to 0$ and $\beta_M \to 1$. Note that $\beta_A$ and $\beta_M$ do not change at the same rates (Fig. \ref{fig:beta}). When we increase sample size so that each sample includes an average of two species ($\bar{\alpha} = 2$), then $\beta_A = 98$ and $\beta_M =50$. If each sample were big enough to have on average 50 species ($\bar{\alpha} = 50$), then $\beta_A = 50$ and $\beta_M =2$. So, the $\beta$'s do not change at the same rate (Fig. \ref{fig:beta}). \begin{figure}[ht] \centering \includegraphics[width=.5\linewidth]{beta.pdf} \caption{Relations of $\beta_A$ (with additive partitioning) and $\beta_M$ (with multiplicative partitioning) to $\bar{\alpha}$, for a fixed $\gamma=500$ species. In our example, we defined diversity as species richness, so the units of $\beta_A$ and $\alpha$ are number of species per sample, and $\bar{\alpha}$ is the mean number of species in a sample.} \label{fig:beta} \end{figure} <>= a <- seq(1,100, by=.1) g <- max(a) bA <- g-a bM <- g/a matplot(a, cbind(bA, bM), type='l', xlab=quote(bar(italic(alpha))), ylab=quote(italic(beta)), col=1) @ Multiplicative $\beta_M$ is sometimes thought of as the number of independent ``communities'' in a set of samples. This would make sense if our sampling regime were designed to capture representative parts of different communities. For example, if we sampled an elevational gradient, or a productivity gradient, and our smallest sample was sufficiently large so as to be representative of that point along the gradient\footnote{One type of sample that attempts to do this is a relev\'e.} then $\beta_M$ could provide a measure of the relative turnover in composition or ``number of different communities.'' However, we know that composition is also predicted to vary randomly across the landscape \cite{Condit2002}. Therefore, if each sample is small, and not really representative of a ``community,'' then the small size of the samples will inflate $\beta_M$ and change the interpretation. As an example, consider the BCI data, which consists of 50 contiguous 1\,ha plots. First, we find $\gamma$ (all species in the data set, or the number of columns), and $\bar{\alpha}$ (mean species number per 1\,ha plot). <<>>= (gamma <- dim(BCI)[2]) (alpha.bar <- mean( specnumber(BCI) )) @ Next we find additive $\beta$-diversity and multiplicative $\beta$-diversity. <>= (beta.A <- gamma - alpha.bar) (beta.M <- gamma/alpha.bar) @ Now we interpret them. These plots are located in a relatively uniform tropical rainforest. Therefore, they each are samples drawn from a single community type. However, the samples are small. Therefore, each 1\,ha plot ($10^4$\,m$^2$ in size) misses more species than it finds, on average ($\beta_A > \bar{\alpha}$). In addition, $\beta_M=2.48$, indicating a great deal of turnover in species composition. We could mistakenly interpret this as indicating something like $\sim 2.5$ independent community types in our samples. Here, however, we have a single community type --- additive partitioning is a little simpler and transparent in its interpretation. For other meanings of $\beta$-diversity, linked to neutral theory, see \cite{Condit2002, Morlon:2008tx}. \subsection{An example of diversity partitioning} Let us consider a study of moth diversity by Keith Summerville and Thomas Crist \cite{Summerville2003, Summerville:2004ym}. The subset of their data presented here consists of woody plant feeding moths collected in southwest Ohio, USA. Thousands of individuals were trapped in 21 forest patches, distributed in two adjacent ecoregions (12 sites - \index{North Central Tillplain}North Central Tillplain [NCT], and 9 sites - \index{Western Allegheny Plateau}Western Allegheny Plateau [WAP], Fig. \ref{fig:moths}). This data set includes a total of 230 species, with 179 species present in the NCT ecoregion and 173 species present in the WAP ecoregion. From these subtotals, we can already see that each ecoregion had most of the combined total species ($\gamma$). <>= moths <- read.table("Summervillle_CristMothsALL.txt", header=TRUE) rowsnct <- c(4,7,9,11,6,3,5,1,10,12,8,2) rowswap<- c(16, 20,21,14,17,19,15,18,13) #identify mothsnct <- moths[rowsnct,] mothswap <- moths[rowswap,] library(maps) #wap <- map.text("state",xlim=c(-83, -82.3), ylim=c(39.4,39.7)) #map.text("state",xlim=c(-85.1, -84.5), ylim=c(39.4,39.6)) #text(moths[['long']], moths[['lat']], moths[['site']], cex=.5) #identify(wap) #locator() map.text("state",xlim=c(-85.5, -82), ylim=c(38.9,39.95)) text(-85.5, 39.9, labels="Indiana", adj=c(0,1)) text(-83, 39.9, labels="Ohio", adj=c(1,1)) text(-84.5, 38.5, labels="Kentucky", adj=c(0,0)) map.axes() points(moths[['long']], moths[['lat']], pch=1) n <- dim(mothsnct)[1] x <- cos( seq(pi, 1.25*2*pi, length=n) ) y <- sin(seq(pi, 1.25*2*pi, length=n) ) x2 <- .3*x + mean(mothsnct[['long']]) y2 <- .25*y + mean(mothsnct[['lat']]) segments(mothsnct[['long']], mothsnct[['lat']], x2, y2) x3 <- .35*x + mean(mothsnct[['long']]) y3 <- .3*y + mean(mothsnct[['lat']]) text(x3,y3, mothsnct[['spp']]) n <- dim(mothswap)[1] x <- cos( seq(pi, 1.25*2*pi, length=n) ) y <- sin(seq(pi, 1.25*2*pi, length=n) ) x2 <- .3*x + mean(mothswap[['long']]) y2 <- .25*y + mean(mothswap[['lat']]) x2 <- .3*x + -82.531 y2 <- .25*y + 39.584 segments(mothswap[['long']], mothswap[['lat']], x2, y2) x3 <- .35*x + -82.531 y3 <- .3*y + 39.584 text(x3,y3, mothswap[['spp']], adj=c(.75,0.5)) @ We will partition richness at three spatial scales: sites within ecoregions ($\bar{\alpha}_1$), ecoregions ($\bar{\alpha}_2$), and overall ($\gamma$). This will result in two $\beta$-diversities: $\beta_1$ among sites within each ecoregion, and $\beta_2$ between ecoregions. The relations among these are straightforward. \begin{gather} \label{eq:part2} \bar{\alpha}_2 = \bar{\alpha}_1 + \beta_1\\ \gamma = \bar{\alpha}_2 + \beta_2 \\ \gamma = \bar{\alpha}_1 + \beta_1 + \beta_2 \end{gather} To do this in \R, we merely implement the above equations using the data in Fig. \ref{fig:moths} \cite{Summerville:2004ym}. First, we get the average site richness, $\bar{\alpha}_1$. Because we have different numbers of individuals from different site, and richness depends strongly on the number of individuals in our sample, we may want to weight the sites by the number of individuals. However, I will make the perhaps questionable argument for the moment that because trapping effort was similar at all sites, we will not adjust for numbers of individuals. We will assume that different numbers of individuals reflect different population sizes, and let number of individuals be one of the local determinants of richness. <>= data(moths) <<>>= a1 <- mean(moths[['spp']]) @ Next we calculate average richness richness for the ecoregions. Because we had 12 sites in NCT, and only nine sites in WAP for what might be argued are landscape constraints, we will use the weighted average richness, adjusted for the number of sites.\footnote{The arithmetic mean is $\sum{a_i Y_i}$, where all $a_i = 1/n$, and $n$ is the total number of observations. A weighted average is the case where the $a_i$ represent unequal \emph{weights}, often the fraction of $n$ on which each $Y_i$ is based. In both cases, $\sum{a}=1$.} We also create an object for $\gamma=230$. <<>>= a2 <- sum(c(NCT=179, WAP=173) * c(12,9)/21) g <- 230 @ Next, we get the remaining quantities of interest, and show that the partition is consistent. <<>>= b1 <- a2-a1 b2 <- g-a2 abg <- c(a1=a1, b1=b1, a2=a2, b2=b2, g=g) abg a1 + b1 + b2 == g <>= par(las=2) bp <- barplot(abg[c(1,2,4,5)], hor=T, xlim=c(0, 250), space=1, names=expression(alpha[1],beta[1],beta[2],gamma)) text(abg[1]/2, bp[1], c("Mean site\nrichness"), cex=.8) text(abg[2]/2, bp[2], c("Species not in sites,\nbut within ecoregion"), cex=.8) text(abg[4]/2, bp[3], c("Species not\nin ecoregions"), cex=.8) text(abg[5]/2, bp[4], c("Regional richness (all sampled species)"), cex=.8) @ \begin{figure}[ht] \centering \includegraphics[width=.7\linewidth]{partplot} \caption{Hierarchical partitioning of moth species richness in forest patches \cite{Summerville:2004ym}. See Fig. \ref{fig:moths} for geographical locations.} \label{fig:partplot} \end{figure} \index{hierarchical partitioning|see{diversity partitioning}}\index{hierarchical partitioning!diversity partitioning} The partitioning reveals that $\beta_1$ is the largest fraction of overall $\gamma$-richness (Fig. \ref{fig:partplot}). This indicates that in spite of the large distance between sampling areas located in different ecoregions, and the different soil types and associated flora, most of the variation occurs \emph{among sites within regions}. If there had been a greater difference in moth community composition among ecoregions, then $\beta_2$-richness would have made up a greater proportion of the total. These calculations show us how simple this additive partition can be, although more complicated analyses are certainly possible. It can be very important to weight appropriately the various measures of diversity (e.g., the number of individuals in each sample, or number of samples per hierarchical level). The number of individuals in particular has a tremendous influence on richness, but has less influence on Simpson's diversity partitioning. The freely available \index{PARTITION software}PARTITION software will perform this additive partitioning (with sample sizes weights) and perform statistical tests \cite{Veech:2007cr}. @ \subsection{Species--area relations} The relation between the number of species found in samples of different area has a long tradition \cite{Arrhenius1921, Bell2001, MacArthur1972, MacArthur1963, MacArthur1967, Preston:1960uq, Rosenzweig1995}, and is now an important part of the metastasizing subdiscipline of \emph{macroecology} \cite{Brown1995, Harte:2008nq}. Most generally, the \index{species--area relation}species--area relation (\index{SAR|see{species--area relation}}SAR) is simply an empirical pattern of the number of species found in patches of different size, plotted as a function of the sizes of the respective patches (Fig. \ref{fig:sars}). These patches may be isolated from each other, as islands in the South Pacific \cite{MacArthur1963}, or mountaintops covered in coniferous forest surrounded by a sea of desert \cite{Brown1977}, or calcareous grasslands in an agricultural landscape \cite{Hamback:2007yj}. On the other hand, these patches might be nested sets, where each larger patch contains all others \cite{Crawley2001a, Plotkin:2000qe}. Quantitatively, the relation is most often proposed as a simple power law, \begin{equation} \label{eq:area1} R = cA^z \end{equation} where $R$ is the number of species in a patch of size $A$, and $c$ and $z$ are fitted constants. This is most often plotted as a log--log relation, which makes it linear. \begin{equation} \label{eq:area2} \log\left(R\right) = b + zA \end{equation} where $b$ is the intercept (equal to $\log c$) and $z$ is the slope. \begin{figure}[ht] \centering \subfloat[Arithmetic scale]{% \includegraphics[width=.48\linewidth]{SARs1.pdf} } \subfloat[Logarithmic scale]{% \includegraphics[width=.48\linewidth]{SARs2.pdf} } \caption{Power law species--area relations.} \label{fig:sars} \end{figure} \medskip \noindent \begin{boxedminipage}{\linewidth} {\footnotesize \paragraph{Drawing power law species--area relations (Fig. \ref{fig:sars})} Here we simply draw some species area curves. <>= A <- 10^1:10; c <- 1.5; z <- 0.25 curve(c*x^z, 10, 10^10, n=500, ylab="No. of Species", xlab="Area (ha)") <>= A <- 10^1:10; c <- 1.5; z <- 0.25 curve(log(c,10) + z*x, 1, 10, ylab=quote(log[10]("No. of Species")), xlab=quote(log[10]("Area (ha)")) ) @ } \end{boxedminipage} \medskip \begin{figure}[ht] \centering \includegraphics[width=.6\linewidth]{SARmoths.pdf} \caption{Fitted power law species--area relations.} \label{fig:sarsfit} \end{figure} \medskip \noindent \begin{boxedminipage}{\linewidth} {\footnotesize \paragraph{Fitting a species--area relation (Fig. \ref{fig:sarsfit})} Here we fit a species--area curve to data, and examine the slope. We could fit a nonlinear power relation ($y=cA^z$); this would be appropriate, for instance, if the residual noise around the line were of the same magnitude for all patch areas. We could use reduced major axis regression, which is appropriate when there is equivalent uncertainty or error on both $x$ and $y$. Last (and most often), we could use a simple linear \index{regression}regression on the log-transformed data, which is appropriate when we know $x$ to a high degree of accuracy, but measure $y$ with some error, and the transformation causes the residual errors to be of similar magnitude at all areas. We start with the last (log-transformed). Here we plot the data, and fit a linear model to the common log-transformed data. <>= moths <- read.table("Summervillle_CristMothsALL.txt", header=TRUE) names(moths) plot(log(spp,10) ~ log(area,10), moths) mod <- lm(log(spp,10) ~log(area,10), data=moths) abline(mod) try(rm("a")); try(rm("z")) mod.nonlin <- nls(spp ~ a*area^z, start=list(a=1, z=.2),data=moths) curve(log(coef(mod.nonlin)[1],10) + x*coef(mod.nonlin)[2], 0, 3, add=TRUE, lty=2) <<>>= plot(log(spp,10) ~ log(area,10), moths) mod <- lm(log(spp,10) ~log(area,10), data=moths) abline(mod) @ Next we fit the nonlinear curve to the raw data, and overlay that fit (on the log scale). <<>>= mod.nonlin <- nls(spp ~ a*area^z, start=list(a=1, z=.2),data=moths) curve(log(coef(mod.nonlin)[1],10) + x*coef(mod.nonlin)[2], 0, 3, add=TRUE, lty=2) @ } \end{boxedminipage} \medskip \medskip \noindent \begin{boxedminipage}{\linewidth} {\footnotesize \paragraph{Assessing species--area relations} Note that in Figure \ref{fig:sarsfit}, the fits differ slightly between the two methods. Let's compare the estimates of the slope --- we certainly expect them to be similar, given the picture we just drew.\index{confint@\texttt{confint}} <<>>= confint(mod) confint(mod.nonlin) @ We note that the estimates of the slopes are quite similar. Determining the better of the two methods (or others) is beyond the scope of this book, but be aware that methods can matter. } \end{boxedminipage} \medskip The major impetus for the species--area relation came from (i) Preston's work on connections between the species--area relation and the log-normal species abundance distribution \cite{Preston:1960uq, Preston1962}, and (ii) MacArthur and Wilson's theory of \index{island biogeography}island biogeography\footnote{(Originally proposed in a paper entitled ``An Equilibrium Theory of Insular Zoogeography'' \cite{MacArthur1963})} \cite{MacArthur1967}. Preston posited that, given the log-normal species abundance distributions (see above), then increasingly large samples should accumulate species at particular rates. Direct extensions of this work, linking neutral theory and \index{maximum entropy theory}maximum entropy theory to species abundances and species--area relations continues today \cite{Bell2001, Harte:2008nq, Hubbell2001} \subsubsection{Island biogeography} MacArthur and Wilson proposed a simple theory wherein the number of species on an oceanic island was a function of the immigration rate of new species, and extinction rate of existing species (Fig. \ref{fig:IE}). The number of species at any one time was \emph{a dynamic equilibrium}, resulting from both slow inevitable extinction and slow continual arrival of replacements. Thus species composition on the island was predicted to change over time, that is, to undergo turnover. Let us imagine immigration rate, $y$, as a function of the number of species already on an island, $x$ (Fig. \ref{fig:IE}). This relation will have a negative slope, because as the number of species rises, that chance that a new individual actually represents a new species will decline. The immigration rate will be highest when there are no species on the island, $x=0$, and will fall to zero when every conceivable species is already there. In addition, the slope should be decelerating (concave up) because some species will be much more likely than others to immigrate. This means that the immigration rate drops steeply as the most likely immigrants arrive, and only the unlikely immigrants are missing. Immigrants may colonize quickly for two reasons. First, as Preston noted, some species are simply much more common than others. Second, some species are much better dispersers than others. Now let us imagine extinction rate, $y$, as a function of the number of species, $x$ (Fig. \ref{fig:IE}). This relation will have a positive slope, such that the probability of extinction increases with the number of species. This is predicted to have an accelerating slope (concave-up), for essentially the same sorts of reasons governing the shape of the immigration curve: Some species are less common than others, and therefore more likely to become extinct do to demographic and environmental stochasticity, and second, some species will have lower fitness for any number of reasons. As the number of species accumulates, the more likely it will become that these extinction-prone species (rare and/or lower fitness) will be present, and therefore able to become extinct. The \emph{rate of change} of the number of species on the island, $\Delta R$, will be the difference between immimigration, $I$, and extinction, $E$, or \begin{equation} \label{eq:deltaR} \Delta R = I - E. \end{equation} When $\Delta R = 0$, we have an equilibrium. If we knew the quantitative form of immigration and extinction, we could solve for the equilibrium. That equilibrium would be the point on the $x$ axis, $R$, where the two rates cross (Fig. \ref{fig:IE}). In MacArthur and Wilson's theory of island biogeography, these rates could be driven by the \emph{sizes} of the islands, where \begin{itemize} \item larger islands had lower extinction rates because of larger average population sizes. \item larger islands had higher colonization rates because they were larger targets for dispersing species. \end{itemize} The distance between an island and sources of propagules was also predicted to influence these rates. \begin{itemize} \item Islands closer to mainlands had higher colonization rates of new species because more propagules would be more likely to arrive there, \item The effect of area would be more important for islands far from mainlands than for islands close to mainlands. \end{itemize} Like much good theory, these were simple ideas, but had profound effects on the way ecologists thought about communities. Now these ideas, of dispersal mediated coexistence and landscape structure, continue to influence community ecologists \cite{Brown1995, Hubbell2001, Leibold:2004fe}. \begin{figure}[ht] \centering \includegraphics[width=.7\linewidth]{IE} \caption{Immigration and extinction curves for the theory of island biogeography. The declining curves represent immigration rates as functions of the number of species present on an island. The increasing curves represent extinction rates, also as functions of island richness. See text for discussion of the heights of these curves, i.e., controls on these rates. Here the dashed lines represent an island that is shrinking in size.} \label{fig:IE} \end{figure} \subsubsection{Drawing immigration and extinction curves} It would be fun to derive a model of immigration and extinction rates from first principles \cite{MacArthur1963}, but here we can illustrate these relations with some simple phenomenological graphs. We will assume that immmigration rate, $I$, can be represented as a simple negative exponential function $\exp (I_0 - iR)$, where $I_0$ is the rate of immigration to an empty island, and $i$ is the per species negative effect on immigration rate. <>= I0 <- log(1) b <- .1 curve(exp(I0 - b*x),0, 50, xlab='No. of Species (R)', ylab="Rate (I or E)") d <- .01 curve(exp(d*x)-1, 0, 50, add=TRUE) deltaR <- function(R) { (exp(I0-b*R) - (exp(d*R)-1))^2 } S1 <- optimize(f=deltaR, interval=c(1,50))$minimum I1 <- exp(I0-b*S1) segments(S1, -.1, S1, I1) I0 <- log(1/2) curve(exp(I0-b*x), 0,50, add=TRUE, lty=2) d <- .014 curve(exp(d*x)-1, 0, 50, add=TRUE, lty=2) S2 <- optimize(f=deltaR, interval=c(1,50))$minimum I2 <- exp(I0-b*S2) segments(S2, -.1, S2, I2, lty=2) <<>>= I0 <- log(1) b <- .1 curve(exp(I0 - b*x),0, 50, xlab='No. of Species (R)', ylab="Rate (I or E)") @ Note that extinction rate, $E$, must be zero if there are no species present. Imagine that extinction rate is a function of density and that average density declines as the number of species increases, or $\bar{N} = 1/R$.\footnote{Why would this make sense ecologically?}. <<>>= d <- .01 curve(exp(d*x)-1, 0, 50, add=TRUE) @ We subtract 1 merely to make sure that $E=0$ when $R=0$. The number of species, $R$, will result from $\Delta R = 0= I-E$, the point at which the lines cross. \begin{gather} I = e^{I_0 - bR}\\ E = e^{dR}-1 \\ \delta R = 0 = I-E \end{gather} Here we find this empricially by creating a function of $R$ to minimize --- we will minimize $(I-E)^2$; squaring the difference gives the quantity the convenient property that the minimum will be approached from either positive or negative values. <<>>= deltaR <- function(R) { (exp(I0-b*R) - (exp(d*R)-1))^2 } @ We feed this into an optimizer for one-parameter functions, and specify that we know the optimum will be achieved somewhere in the interval between 1 and 50. <<>>= optimize(f=deltaR, interval=c(1,50))[["minimum"]] @ The output tells us that the minimum was achieved when $R \approx 16.9$. Now imagine that rising sea level causes island area to shrink. What is this predicted to do? It could \begin{enumerate} \item reduce the base immigration rate because the island is a smaller target, \item increase extinction rate because of reduced densities.\footnote{We might also predict an increased extinction rate because of reduced rescue effect (Chapter 4).} \end{enumerate} Let us represent reduced immigration rate by reducing $I_0$. <<>>= I0 <- log(1/2) curve(exp(I0 - b*x), 0,50, add=TRUE, lty=2) @ Next we increase extinction rate by increasing the per species rate. <<>>= d <- .014 curve(exp(d*x)-1, 0, 50, add=TRUE, lty=2) @ If we note where the rates cross, we find that the number of predicted species has declined. With these new immigration and extinciton rates the predicted number of species is <<>>= optimize(f=deltaR, interval=c(1,50))[["minimum"]] @ or 11 species, roughly a 35\% decline ($(17-11)/17 = 0.35$). The beauty of this theory is that it focuses our attention on \emph{\index{landscape level processes}landscape level processes}, often outside the spatial and temporal limits of our sampling regimes. It specifies that any factor which helps determine the immigration rate or extinction rate, including island area or proximity to a source of propagules, is predicted to alter the equilibrium number of species at any point in time. We should further emphasize that the \emph{identity} of species should change over time, that is, undergo turnover, because new species arrive and old species become extinct. The rate of turnover, however, is likely to be slow, because the species that are most likely to immigrate and least likely to become extinct will be the same species from year to year. \subsection{\index{species--area relations!diversity partitioning}Partitioning species--area relations} You may already be wondering if there is a link island biogeography and $\beta$-diversity. After all, as we move from island to island, and as we move from small islands to large islands, we typically encounter additional species, and that is what we mean by $\beta$-diversity. Sure enough, there are connections \cite{Crist:2006lr}.\index{diversity partitioning!species--area relations} Let us consider the moth data we used above (Fig. \ref{fig:moths}). The total number of species in all of the patches is, as before, $\gamma$. The average richness of these patches is $\bar{\alpha}$, and also note that \emph{part of what determines that average is the area of the patch}. That is, when a species is missing from a patch, part of the reason might be that the patch is smaller than it could be. We will therefore partition $\beta$ into yet one more category: species missing due to patch size, $\beta_{area}$. This new quantity is the average difference between $\bar{\alpha}$ and the diversity \emph{predicted} for the largest patch (Fig. \ref{fig:SARpart}). In general then, \begin{equation} \label{eq:partSAR} \beta = \beta_{area} + \beta_{replace} \end{equation} where $\beta_{replace}$ is the average number of species missing that are \emph{not} explained by patch size. In the context of these data (Fig. \ref{fig:SARpart}), we now realize that $\beta_1 = \beta_{area} + \beta_{ecoregion}$, so the full partition becomes \begin{equation} \label{eq:partfull} \gamma = \bar{\alpha}_1 + \beta_{area} + \beta_{ecoregion} + \beta_{geogr.region} \end{equation} where $\beta_{replace} = \beta_{ecoregion} + \beta_{geogr.region}$. Note that earlier in the chapter, we did not explore the effect of area. In that case, $\beta_{ecoregion}$ included both the effect of area and the effect of ecoregion; here we have further partitioned this variation into variation due to patch size, as well as variation due to ecoregion. This reduces the amount of unexplained variation among sites within each ecoregion. Let's calculate those differences now. We will use quantities we calculated above for $\bar{\alpha}_1$, $\bar{\alpha}_2$, $\gamma$, and a nonlinear species--area model from above. We can start to create a graph similar to Fig. \ref{fig:SARpart}. <<>>= plot(spp ~ area, data=moths, ylim=c(30,230), #log='y', xlab="Area (ha)", ylab="No. of Species (R)") curve(coef(mod.nonlin)[1] * x^coef(mod.nonlin)[2], 0, max(moths[['area']]), add=TRUE, lty=2, lwd=2) #text(150, 200, # bquote(italic("R") == .(round(coef(mod.nonlin)[1],1)) * # italic("A")^.(round(coef(mod.nonlin)[2],2)) ) # ) abline(h=g, lty=3) text(275, g, quote(gamma), adj=c(.5,1.5), cex=1.5) @ Next we need to find the \emph{predicted richness} for the maximum area. We use our statistical model to find that. <>= (MaxR <- predict(mod.nonlin, list(area=max(moths[['area']]) )) ) @ We can now find $\beta_{area}$, $\beta_{eco}$ and $\beta_{geo}$. <<>>= b.area <- MaxR - a1 b.eco <- a2-(b.area+a1) b.geo <- g - a2 @ Now we have partitioned $\gamma$ a little bit more finely with a beastiary of $\beta$'s, where \begin{itemize} \item $\bar{\alpha}_1$ is the average site richness. \item $\beta_{area}$ is the average number of species not observed, due to different patch sizes. \item $\beta_{eco}$ is the average number of species not observed at a site, is not missing due to patch size, but is in the ecoregion. \item $\beta_{geo}$ is the average number of species not found in the samples from different ecoregions. \end{itemize} Finally, we add lines to our graph to show the partitions. <<>>= abline(h=g, lty=3) abline(h=b.eco+b.area+a1, lty=3) abline(h=b.area+a1, lty=3) abline(h=a1, lty=3) @ Now we have further quantified how forest fragment area explains moth species richness. Such understanding of the spatial distribution of biodiversity provides a way to better quantify patterns governed by both dispersal and habitat preference, and allows us to better describe and manage biodiversity in human-dominated landscapes. \begin{figure}[ht] \centering \includegraphics[width=.6\linewidth]{SARpart} \caption{Combining species--area relations with additive diversity partitioning. Forest fragment area explains relatively little of the diversity which accumulates in isolated patches distributed in space. However, it is likely that area associated with the collection of samples (i.e., the distances among fragments) contributes to $\beta_{eco}$ and $\beta_{geo}$.} \label{fig:SARpart} \end{figure} <>= par(mar=c(5,4,1,3), las=1) plot(spp ~ area, data=moths, ylim=c(0,250), #log='y', xlab="Area (ha)", ylab="No. of Species (R)") cols <- hcl(c(150), alpha=c(1,.5,.3,.1)) polygon(c(-20,320,320,-200), c(0,0,a1,a1), border=FALSE, col=cols[1]) polygon(c(-20,320,320,-200), c(0,0,a1+b.area,a1+b.area), border=FALSE, col=cols[2]) polygon(c(-20,320,320,-200), c(0,0,a1+b.area+b.eco,a1+b.area+b.eco), border=FALSE, col=cols[3]) polygon(c(-20,320,320,-200), c(0,0,g,g), border=FALSE, col=cols[4]) box() points(moths[['area']], moths[['spp']]) curve(coef(mod.nonlin)[1] * x^coef(mod.nonlin)[2], min(moths[['area']])/2, max(moths[['area']]), add=TRUE) abline(h=g, lty=3) abline(h=b.eco+b.area+a1, lty=3) abline(h=b.area+a1, lty=3) abline(h=a1, lty=3) text(225, g, quote(gamma), adj=c(-.5,0), cex=1.5) arrows(225, 0, 225, g, code=3, length=.05) #title(sub=quote(beta["replace"]==beta['eco']+beta['geo'])) #mtext(quote(beta["replace"]==beta['eco']+beta['geo'])) #text(175, 300, quote(beta["replace"]==beta['eco']+beta['geo']), cex=1.2) #arrows(200, MaxR, 200, g, code=3, length=.05) mtext(quote(alpha["1"]), side=4, at=a1, line=.5, cex=1.2) mtext(quote(alpha["2"]), side=4, at=a2, line=.5, cex=1.2) mtext(quote(italic("R"["max"])), side=4, at=a1+b.area, line=.5) #text(280, a1, quote(alpha["1"]), adj=c(0,0), cex=1.2) #text(280, a2, quote(alpha["2"]), adj=c(0,0), cex=1.2) text(150, a2+b2/2, quote(beta["geo"]), adj=c(1.1,0.5), cex=1.2) arrows(150, a1+b.area+b.eco, 150, g, code=3, length=.05) text(100, a1+b1/2, quote(beta["eco"]), adj=c(1.1,0.5), cex=1.2) arrows(100, a1+b.area, 100, a1+b.area+b.eco, code=3, length=.05) text(25, a1+b.area/2, quote(beta["area"]), adj=c(1.1,.5), cex=1.2) arrows(32, a1, 32, a1 + b.area, code=3, length=.05) text(32, a1/2, quote(bar(alpha)), adj=c(1.2,.5), cex=1.2) arrows(32, 0, 32, a1, code=3, length=.05) @ \section{Summary} We have examined communities as multivariate entities which we can describe and compare in a variety of ways. \begin{compactitem} \item Composition includes all species (multivariate data), whereas species diversity is a univariate description of the variety of species present. \item There are many ways to quantify species diversity, and they tend to be correlated. The simplest of these is richness (the number of species present), whereas other statistics take species' relative abundances into account. \item Species abundance distributions and rank abundance distributions are analogous to probability distributions, and provide more thorough ways to describe the patterns of abundance and variety in communities. These all illustrate a basic law of community ecology: most species are rare. Null models of community structure and processes make predictions about the shape of these distributions. \item Ecological neutral theory provides a dynamical model, not unlike a null model, which allows quantitative predictions relating demographic, immigration, and speciation rates, species abundance distributions, and patterns of variation in space and time. \item Another law of community ecology is that the number of species increases with sample area and appears to be influenced by immigration and extinction rates. \item We can partition diversity at different spatial scales to understand the structure of communities in landscapes. \end{compactitem} \section*{Problems} \addcontentsline{toc}{section}{Problems} \begin{table}[h] \centering \caption{Hypothetical data for Problem 1.} \label{tab:prob1} \begin{tabular}[c]{lccc} \hline \hline Site~ & ~Sp. A~ & ~Sp. B~ & ~Sp. C~ \\ \hline Site 1 & 0 & 1 & 10 \\ Site 2 & 5 & 9 & 10 \\ Site 3 & 25 & 20 & 10 \\ \hline \end{tabular} \end{table} \begin{prob} How different are the communities in Table \ref{tab:prob1}?\\ (a) Convert all data to relative abundance, where the relative abundance of each site sum to 1.\\ (b) Calculate the Euclidean and Bray-Curtis (S{\o}rensen) distances between each pair of sites for both relative and absolute abundances.\\ (c) Calculate richness, Simpson's and Shannon-Wiener diversity for each site. \end{prob} \begin{prob} Use rarefaction to compare the tree richness in two 1\,ha plots from the BCI data in the \texttt{vegan} package. Provide code, and a single graph of the expectations for different numbers of individuals; include in the graph some indication of the uncertainty. \end{prob} \begin{prob} Select one of the 1\,ha BCI plots (from the \texttt{vegan} package), and fit three different rank abundance distributions to the data. Compare and contrast their fits. \end{prob} \begin{prob} Simulate a neutral community of 1000 individuals, selecting the various criteria on yur own. Describe the change through time. Relate the species abundance distributions that you observe through time to the parameters you choose for the simulation. \end{prob} \begin{prob} Using the \texttt{dune} species data (\texttt{vegan} package), partition species richness into $\bar{\alpha}$, $\beta$, and $\gamma$ richness, where rows are separate sites. Do the same thing using Simpson's diversity. \end{prob} % \begin{figure}[ht] % \centering % \includegraphics[width=.5\linewidth]{sugi.pdf} % \caption{[Coming soon: simulated broken stick distributions].} % \label{fig:sugi} % \end{figure} % \medskip \noindent % \begin{boxedminipage}{\linewidth} % {\footnotesize % \paragraph{[Coming soon: Simulating Sugihari's minimal model(Fig. \ref{fig:sugi})]} % We can illustrate a stick-breaking niche-structured community. Let us let the stick be the interval from zero to one. Once we have broken the stick, each species abundance be a fraction of this segement, so the sum of species aunbances will equal 1. % Breaking this sequentially requires that we use a \emph{for-loop} for $R$ species. % <<>>= % R2 <- 25; spp <- numeric(R2) % brkpnt <- runif(1) % old.frag <- 1 % new.frags <- c(old.frag*brkpnt, old.frag*(1-brkpnt)) % <>= % plot(1,1) % @ % To be continued. % } % \end{boxedminipage} \medskip