\SweaveOpts{ eval=TRUE, prefix=FALSE, width=4, height=4, eps=FALSE, include=FALSE} % Use this file as a template for yo \chapter{Density-independent Demography} \label{DIDemo} % Always give a unique label % use \chaptermark{} % to alter or adjust the chapter heading in the running head <>= setwd("~/Documents/Projects/Book/draft/Chap02") 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=75) ## palette(gray((0:8)/8)) library(lattice) trellis.par.set(canonical.theme(color=FALSE)) library(primer) @ Different populations have different numbers of individuals of different ages. Consider the human populations of Mexico and Sweden in 1990. Mexico had more individuals in total than Sweden, and a larger fraction of their population was of child bearing age or younger (Figs. \ref{fig:AgeDistM}, \ref{fig:AgeDistS}). \begin{figure}[H] \centering \subfloat[Mexico]{\includegraphics[width=.3\textwidth]{AgeDistM.pdf} \label{fig:AgeDistM}} \subfloat[Sweden]{\includegraphics[width=.3\textwidth]{AgeDistS.pdf} \label{fig:AgeDistS}} \subfloat[Mexico \emph{vs}. Sweden]{\includegraphics[width=.35\textwidth]{AgeSpecFert.pdf} \label{fig:AgeSF}} \caption{Demography of human populations of Mexico and Sweden. Based on 1990 data from US Census Bureau, Population Division, International Programs Center.} \end{figure} In addition, the \index{age-specific fertility}age-specific fertility rate is higher in Mexico, especially for younger women (Fig. \ref{fig:AgeSF}). How did this happen, and why does Mexico have so many young people? What are the consequences of this for their culture, their use of resources, their domestic and foreign policies, and their future population growth? How about Sweden? Demography is the study of populations with special attention to age or stage structure \cite{Lincoln:1998ta}. Originally, age-based human \index{demography}demography was the provenance of actuaries who helped governments keep track of the number citizens of different ages and thus, for instance, know how many would be available for conscription into the military.% \footnote{In his chapter entitled ``Interesting Ways to Think about Death'' G.E. Hutchinson \cite{Hutchinson1978} cites C.~F. Trenerry, E.~L. Gover and A.~S. Paul (\emph{The Origins and Early History of Insurance}, London, P.~S. King \& Sons, Ltd.) for description of early Roman actuarial tables.}% The demography of a population is the age (or stage) structure and the survival, fertility, and other demographic rates associated with those ages or life history stages. \emph{\index{age structure}Age structure} is the number or relative abundance of individuals of different ages or age classes. \emph{\index{stage structure}Stage structure} is the number or relative abundance of individuals of different stages. Stages are merely useful categories of individuals, such as size classes (e.g. diameters of tropical trees) or \index{life history stages}life history stages (e.g. egg, larvae, and adult anurans). Stages are particularly useful when (i) age is difficult to determine, and/or (ii) when stage is a better predictor of demographic rates (e.g. birth, death, survival) than is age. Demography is, in part, the study of how demographic rates vary among ages or stages, and the consequences of those differences. There are a few ways to study a population's demography, and all ecology text books can provide examples. \emph{Life tables} are lists of important demographic parameters such as survivorship, birth and death rates each age or age class. Commonly, both age and stage based demography now take advantage of matrix algebra to simplify and synthesize age and stage specific demography \cite{Caswell:2001gu}. This approach is essential when individuals don't proceed through stages in a simple sequential manner, perhaps reverting to an ``earlier'' stage. When used with age-based demography, these matrices are referred to as Leslie matrices \cite{Leslie:1945la}. L.~P. \index{Lefkovitch}Lefkovitch \cite{Lefkovitch:1965mq} generalized this approach to allow for complex demography. This could include, for instance, regression from a large size class to a smaller size class (e.g. a two-leaved woodland perennial to a one-leaved stage). Using matrices to represent a population's demography allows us to use the huge workshop of linear algebra tools to understand and predict growth in structured populations. Let's start with a hypothetical example. \section{A Hypothetical Example} Pretend you are managing a small nature reserve and you notice that a new invasive species, spotted mallwort (\emph{Capitalia globifera}),\footnote{Not a real species.} is popping up everywhere. You think you may need to institute some control measures, and to begin, you would like to understand its life cycle and population growth. Careful examination of the flowers reveals perfect flowers,\footnote{Individual flowers possess both female and male reproductive structures.} and you find from the literature that mallwort is a perennial capable of self-fertilizing. The seeds germinate in early fall, grow the following spring and early summer to a small adult that has 2--3 leaves and which sometimes produce seeds. In the second year and beyond, if the plants survive, they grow to large adults which have four or more leaves and only then do they produce the majority of their seeds. The seeds do not seem to survive more than one year, so there is probably no seed bank. You summarize what you have learned in the form of a \emph{\index{life cycle graph}life cycle graph} (Fig. \ref{fig:malldemo}). Demographers use a life cycle graph to summarize stages that may be observed at a single repeated point in time (e.g., when you go out to explore in June). It also can include the probabilities that an individual makes the \index{transition}transition from one stage to another over one time step (e.g. one year), as well as the fecundities. \begin{figure}[ht] \centering \includegraphics[width=.5\textwidth]{MallDemo2} \label{fig:malldemo} \caption{Life cycle graph of the imaginary spotted mallwort (\emph{Capitalia globifera}). $P_{ij}$ is the probability that an individual in stage $j$ transitions to stage $i$ over a single fixed time interval between samples. $F_{i}$ the number of progeny (transitioning into stage 1) produced by an individual in stage $j$. Thus, for mallwort, $P_{21}$ is the probability that a seed (Seeds) makes it into the small adult stage (A-Small). $P_{32}$ is the probability that a small adult shows up as a large adult the next year. $F_2$ is the average fertility for individuals in the small adult stage and $F_3$ is the average fertility for individuals in the large adult stage.} \end{figure} As the manager responsible for this small reserve, you decide to keep track of this new exotic species. After identifying a general area where the plant seems to have obtained a foothold, you established 50 permanent 1\,m$^2$ sample plots, located randomly throughout the invasion area. Each year, for two years, you sample in early summer when the fruits are on the plants (when the weather is pleasant and you can find interns and volunteers to help). In all plots you tag and count all first year plants (2--3 leaves), and all older plants (4+ leaves). You also are able to count fruits and have determined that there is one seed per fruit. Now that you have your data for two years, you would like to figure out how quickly the population growing. You could simply keep track of the total population size, $N$, or just the large adults. You realize, however, that different stages may contribute very differently to growth, and different stages may be better for focused control efforts. A description, or model, of the population that includes different stages will provide this. We call such a model a \emph{\index{demographic model}demographic model}, and it consists of a \emph{\index{projection!population projection matrix}population projection matrix}. The population projection matrix is a matrix that represents the life cycle graph. We use the projection matrix to calculate all kinds of fun and useful stuff including \begin{itemize} \item The finite rate of increase, $\lambda$ (the asymptotic population growth rate). \item The stable stage distribution (the population structure that would emerge if the demographic rates ($P$, $F$) do not change). \item Elasticty, the relative importance of each transition to $\lambda$. \end{itemize} \subsection{The population projection matrix} The population projection matrix (a.k.a. the transition matrix) is simply an organized collection of the per capita contribution of each stage $j$ to the next stage $i$ in the specified time interval (often one year). These contributions, or transitions, consist of (i) the probabilities that an individual in stage $j$ in one year makes it into stage $i$ the next year, and (ii) the per capita fecundities for reproductive stages (eq. \ref{Lefk1}). Each element of the projection matrix (eq. \ref{Lefk1}) relates its column to its row. Thus $P_{21}$ in our matrix, eq. \ref{Lefk1} is the probability that an individual in stage 1 (seeds; respresented by the column 1) makes it to the next census period and shows up in stage 2 (1 year old small adult; represented by row 2). Similarly, $P_{32}$ is the probability that an individual in stage 2 (a small one year old adult) has made it to the large adult stage at the next census period. The \index{fecundities}fecundities are not probabilities, of course, but are the per capita contribution of propagules from the adult stage to the seed stage. The population projection matrix allows us to multiply all of these transition elements by the abundances in each age class in one year to predict, or \emph{project}, the abundances of all age classes in the following year. \begin{equation} \label{Lefk1} \left( \begin{array}{ccc} 0 & F_{2} & F_{3}\\ P_{21}& 0 & 0\\ 0 & P_{32} & P_{33} \end{array} \right) \end{equation} \subsection{A brief primer on matrices} We refer to matrices by their rows and columns. A \index{matrix algebra}matrix with three rows and one column is a $3 \times 1$ matrix (a ``three by one'' matrix); we \emph{always} state the number of rows first. Matrices are composed of \emph{elements}; an element of a matrix is signified by its row and column. The element in the second row and first column is $a_{21}$. To multiply matrices, we multiply and then sum each row by each column (eq. \ref{exM}). More specifically, we multiply each row element of matrix \vec{A} times each column element of matrix \vec{B}, sum them, and place this sum in the respective element of the final matrix. Consider the matrix multiplication in eq. \ref{exM}. We first multiply each element of row 1 of \vec{A} ($a\quad b$), times the corresponding elements of column 1 of \vec{B} ($m\quad n$), sum these products and place the sum in the first row, first column of the resulting matrix. We then repeat this for each row of \vec{A} and each column of \vec{B} \begin{align} \label{exM} \vec{A} & = \left( \begin{array}{cc} a & b \\ c& d \end{array} \right); \; \vec{B} = \left(\begin{array}{cc} m & o\\ n & p\end{array}\right)\\ \vec{AB} &= \left( \begin{array}{cc} \left( am + bn \right) & \left(ao+bp \right)\\ \left(cm + dn \right) & \left(co + dp \right) \end{array}\right) \end{align} This requires that the number of columns in the first matrix must be the same as the number of rows in the second matrix. It also means that the resulting matrix will have the same number of rows as the first matrix, and the same number of columns as the second matrix. \medskip \noindent \begin{boxedminipage}{\linewidth} {\footnotesize \paragraph{Matrices in \R} Let's define two $2 \times 2$ matrices, filling in one by rows, and the other by columns. <<>>= M <- matrix( 1:4 , nr=2, byrow=T) M N <- matrix( c(10, 20, 30, 40), nr=2) N @ Following our rules above, we would multiply and then sum the first row of $M$ by the first column of $N$, and make this element $a_{11}$ of the resulting matrix product. <<>>= 1 * 10 + 2 * 20 @ We multiply matrices using \texttt{\%*\%} to signify that we mean \emph{matrix} multiplication. <<>>= M %*% N @ } \end{boxedminipage} \medskip \subsection{Population projection} With our spotted mallwort we could multiply our projection matrix by the observed abundances (seeds=$Sd$, small adults - $SA$, large adults - $LA$) to \emph{project} the abundances of all age classes in subsequent years. \begin{equation} \label{Lefk} \left( \begin{array}{ccc} 0 & F_{2} & F_{3}\\ P_{21}& 0 & 0\\ 0 & P_{32} & P_{33} \end{array} \right) \left(\begin{array}{l} N_{Sd}\\ N_{SA}\\ N_{LA}\end{array}\right) =\\ \left( \begin{array}{c} \left(0\times N_{Sd} + F_{2}\times N_{SA} + F_{3}\times N_{LA}\right)\\ \left(P_{21}\times N_{Sd} + 0\times N_{SA} + 0\times N_{LA}\right)\\ \left(0\times N_{Sd} + P_{32}\times N_{SA} + 0\times N_{LA}\right) \end{array}\right) \end{equation} The next step is to create the projection matrix. Let's pretend that over the two years of collecting these data, you found that of the small adults we tagged, about half (50\%) survived to become large adults the following year. This means that the transition from stage 2 (small adults) to stage 3 (large adults) is $P_{32}=0.50$. Of the large adults that we tagged, about 90\% of those survived to the next year, thus $P_{33}=0.90$. We estimated that, on average, each small adult produces 0.5 seeds (i.e. $F_2=0.50$) and each large adult produces 20 seeds (i.e. $F_{3}=20$). Last, we found that, on average, for every 100 seeds (fruits) we counted, we found about 30 small adults (one year olds), meaning that $P_{21}=0.30$. Note that this requires that seeds survive until germination, germinate, and then survive until we census them the following summer. We can now fill in our population projection matrix, \textbf{A}. \begin{equation} \label{Lefk2} \mathbf{A} = \left( \begin{array}{ccc} 0 & F_{2} & F_{3}\\ P_{21}& 0 & 0\\ 0 & P_{32} & P_{33} \end{array} \right) = \left(\begin{array}{ccc} 0 & 0.5 & 20\\ 0.30 & 0 & 0\\ 0 & 0.50 & 0.90 \end{array} \right) \end{equation} Next we can multiply it the projection matrix, \textbf{A}, by the last year for which we have data. \begin{equation} \label{Lefk2} \left(\begin{array}{ccc} 0 & 0.5 & 20\\ 0.3 & 0 & 0\\ 0 & 0.5 & 0.9 \end{array} \right) \left(\begin{array}{c} 100\\ 250\\ 50\end{array}\right) = \left( \begin{array}{c} \left(0\times 100 + 0.5\times 250 + 20\times 50\right)\\ \left(0.3\times 100 + 0\times250 + 0\times 50\right)\\ \left(0\times 100 + 0.5\times250 + 0.9 \times 50\right)\\ \end{array} \right) = \left(\begin{array}{c} 1125\\ 30\\ 170\end{array}\right) \end{equation} If we wanted more years, we could continue to multiply the projection matrix by each year's projected population. We will observe that, at first, each stage increases or decreases in its own fashion (Fig. \ref{fig:MallGrow}), and that over time, they tend to increase in a more similar fashion. This is typical for demographic models. It is one reason why it is important to examine \emph{\index{demography!stage-structured growth}stage-structured growth} rather than trying to lump all the stages together --- we have a much richer description of how the population is changing. \medskip \noindent \begin{boxedminipage}{\linewidth} {\footnotesize \paragraph{Stage structured growth - one step} First, we create a population projection matrix, and a vector of stage class abundances for year zero. <<>>= A <- matrix(c(0, 0.5, 20, 0.3, 0, 0, 0, 0.5, 0.9), nr=3, byrow=TRUE) N0 <- matrix(c(100, 250, 50), ncol=1) @ Now we perform matrix multiplication between the projection matrix and $N_0$. <<>>= N1 <- A%*%N0 N1 @ Note that the first stage declined, while the second and third stages increased. } \end{boxedminipage} \medskip \medskip \noindent \begin{boxedminipage}{\linewidth} {\footnotesize \paragraph{Stage structured growth - multiple steps} Now we project our population over six years, using a for-loop. We use a for-loop, rather than \texttt{sapply}, because each year depends on the previous year (see the Appendix, sec. \ref{sec:iter-acti-apply}). First, we set the number of years we want to project, and then create a matrix to hold the results. We put $N_0$ in the first column. <<>>= years <- 6 N.projections <- matrix( 0, nrow=nrow(A), ncol=years+1) N.projections[,1] <- N0 @ Now we perform the iteration with the for-loop. <<>>= for(i in 1:years) N.projections[,i+1] <- A%*%N.projections[,i] @ Last, we graph the results for each stage (Fig. \ref{fig:MallGrow}). To graph a matrix, \R~is expecting that the data will be in columns, not rows, and so we need to {\em t}ranspose the projection matrix. <>= matplot(0:years, t(N.projections), type="l", lty=1:3, col=1, ylab="Stage Abundance", xlab="Year") legend('topleft', legend=c("Seeds", "Small Adult","Large Adult"), lty=1:3, col=1, bty='n') @ } \end{boxedminipage} \medskip \begin{figure}[ht] \centering \subfloat[Projected Population]{\includegraphics[width=0.48\textwidth]{MallGrow.pdf} \label{fig:MallGrow}} \subfloat[Annual~Growth]{\includegraphics[width=0.48\textwidth]{MallR.pdf}\label{fig:lambda1}} \caption{Population dynamics and annual growth ($R=N_{t+1}/N_t$) of spotted mallwort. Note that stage abundance is on a log-scale.} \end{figure} \subsection{Population growth} We have projected the stages for six years --- what is its observed rate of increase, $R_t=N_{t+1}/N_t$? How do we even think about $R$ and $N$ in stage structured growth? The way we think about and calculate these is to add all the individuals in all stages to get a total $N$, and calculate $R$ with that, as we did in Chapter~1. \begin{equation} \label{eq:L} R_t = N_{t+1}/N_t. \end{equation} If we do that for our mallwort, we can see that $R_t$ changes with time (Fig. \ref{fig:lambda1}). We can summarize the projection as $\vec{n}(t)=\vec{A}^t\vec{n}_{0}$, where $\vec{A}^t$ is $\vec{A}$ multiplied by itself $t$ times. \medskip \noindent \begin{boxedminipage}{\linewidth} {\footnotesize \paragraph{Annual growth rate} Now let's calculate $R_t=N_{t+1}/N_t$ for each year $t$. We first need to sum all the stages, by {\em applying} the {\tt sum} function to each column. <<>>= N.totals <- apply(N.projections, 2, sum) @ Now we get each $R_t$ by dividing all the $N_{t+1}$ by each $N_t$. Using negative indices cause \R~to drop that element. <<>>= Rs <- N.totals[-1] / N.totals[-(years+1)] @ We have one fewer $R$s than we do years, so let's plot each $R$ in each year $t$, rather than each year $t+1$ (Fig. \ref{fig:lambda1}). <>= plot(0:(years-1), Rs, type="b", xlab="Year", ylab="R") @ } \end{boxedminipage} \medskip \section{Analyzing the Projection Matrix} You seem to have a problem on your hands (Fig. \ref{fig:MallGrow}). Being a well-trained scientist and resource manager, several questions come to mind: What the heck do I do now? What is this population likely to do in the future? Can these data provide insight into a control strategy? How confident can I be in these projections? After you get over the shock, you do a little more research on demographic models; Caswell \cite{Caswell:2001gu} is the definitive treatise. You find that, indeed, there is a lot you can do get more information about this population that might be helpful. Once you have obtained the projection matrix, $\vec{A}$, you can analysis it using \emph{eigenanalysis} to estimate \begin{itemize} \item $\lambda$, the finite rate of increase, \item stable stage structure, \item reproductive value, and \item sensitivities and elasticities. \end{itemize} Below, we explain each of these quantities. These quantities will help you determine which stages of spotted mallwort on which to focus eradication efforts. \subsection{Eigenanalysis} Eigenanalysis is a mathematical technique that summarizes multivariate data. Ecologists use \index{eigenanalysis!demographic matrix}eigenanalysis frequently, for (i) multivariate statistics such as ordination, (ii) stability analyses with two or more species, and (iii) analyzing population projection matrices. Eigenanalysis is simply a method to transform a square matrix into independent, orthogonal, and useful chunks --- the eigenvectors and their eigenvalues. In demography, the most useful piece is the dominant eigenvalue and its corresponding vector. Eigenanalysis is a technique that finds all the solutions for $\lambda$ and $\vec{w}$ of \begin{equation} \label{eigsol} \vec{Aw}=\lambda\vec{w}, \end{equation} where $\vec{A}$ is a particular summary of our data. \footnote{For ordination, we analyze a correlation or covariance matrix, and for stability analyses, we use the matrix of pairwise partial differential equations between each pair of species. In these eigenanalyses of a square $i \times j$ matrix $\vec{A}$, we can think of the elements of $\vec{A}$ describing the ``effect'' of stage (or species) $j$ on stage (or species) $i$, where $j$ is a column and $i$ is a row.} With projection matrix analysis, $\vec{A}$ is the projection matrix. $\lambda$ is an \emph{eigenvalue} and $\vec{w}$ is an \emph{eigenvector}. If we write out eq. \ref{eigsol} for a $3 \times 3$ matrix, we would have \begin{equation} \label{eigsollong} \left( \begin{array}{ccc} a_{11} & a_{12} & a_{13}\\ a_{21} & a_{22} & a_{33}\\ a_{31} & a_{32} & a_{33}\\ \end{array} \right) \left( \begin{array}[c]{c} w_{11}\\ w_{21}\\ w_{31} \end{array} \right) =\lambda \left( \begin{array}[c]{c} w_{11}\\ w_{21}\\ w_{31} \end{array} \right) \end{equation} There are typically an infinite number of solutions to this equation, and what eigenanalysis does is find set of solutions that are all independent of each other, and which capture all of the information in \vec{A} in a particularly useful way.\footnote{The number of solutions is infinite because they are just simple multiples of the set found with eigenanalysis.} Typically, the first solution captures the most important features of the projection matrix. We call this the \index{eigenvalue!dominant}dominant eigenvalue, \index{$\lambda_1$|see{eigenvalue, dominant}}$\lambda_1$ and its corresponding eigenvector, $w_1$. The first solution does not capture all of the information; the second solution captures much of the important remaining information. To capture all of the information in \vec{A} requires as many solutions as there are columns of \vec{A}. Nonetheless, the first solution is usually the most useful. \medskip \noindent \begin{boxedminipage}{\linewidth} {\footnotesize \paragraph{Eigenanalysis in \R} Here we perform eigenanalysis on $\vec{A}$. <<>>= eigs.A <- eigen(A) eigs.A @ Each eigenvalue and its corresponding eigenvector provides a solution to eq. \ref{eigsol}. } \end{boxedminipage} \medskip The first, or dominant, eigenvalue is the long term asymptotic finite rate of increase $\lambda$. Its corresponding eigenvector provides the \emph{stable stage distribution}. We can also use eigenanalysis get the \emph{reproductive values} of each stage out of \vec{A}. To be a little more specific, $\vec{w}$ we described above are \emph{right} eigenvectors, so-called because we solve for them with $\vec{w}$ on the right side of $\vec{A}$. We will also generate \emph{left} eigenvectors $\vec{v}$ (and their corresponding eigenvalues), where $\vec{vA}=\lambda\vec(v)$. The dominant left eigenvector provides the reproductive values (section \ref{sec:repval}). \subsection{Finite rate of increase -- $\lambda$} \index{$\lambda$|see{lambda}} The asymptotic annual growth rate finite rate of increase\index{lambda!dominant eigenvalue} is the dominant \emph{eigenvalue} of the projection matrix. Eigenvalues are always referred to with the Greek symbol $\lambda$, and provides a solution to eq. \eqref{eigsol}. The dominant eigenvalue of any matrix, $\lambda_{1}$, is the eigenvalue with the largest absolute value, and it is frequently a complex number.\footnote{When you perform eigenanalysis, it is common to get complex numbers, with real and imaginary parts. Eigenanalysis is, essentially, solving for the roots of the matrix, and, just like when you solved quadratic equations by hand in high school, it is possible to get complex numbers.} With projection matrices, $\lambda_{1}$ will always be positive and real. We use eigenanalysis to solve eq. \ref{eigsol} and give us the answers --- like magic. Another way to find $\lambda_1$ is to simply iterate population growth a very large number of times, that is, let $t$ be very large. As $t$ grows, the annual growth rate, $N_{t+1}/N_t$, approaches $\lambda_1$ (Fig. \ref{fig:converge}). \medskip \noindent \begin{boxedminipage}{\linewidth} {\footnotesize \paragraph{Finding $\lambda$} Next we explicitly find the index position of the largest absolute value of the eigenvalues. In most cases, it is the first eigenvalue. <<>>= dom.pos <- which.max( eigs.A[["values"]] ) @ We use that index to extract the largest eigenvalue. We keep the real part, using \texttt{Re}, dropping the imaginary part. (Note that although the dominant eigenvalue will be real, \R~will include an imaginary part equal to zero ($0i$) as a place holder if \emph{any} of the eigenvalues have a non-zero imaginary part). <<>>= L1 <- Re(eigs.A[["values"]][dom.pos]) L1 @ \texttt{L1} is $\lambda_1$, the aysmptotic finite rate of increase. } \end{boxedminipage} \medskip \medskip \noindent \begin{boxedminipage}{\linewidth} {\footnotesize \paragraph{Power iteration method of eigenanalysis} Because growth is an exponential process, we can figure out what is most important in a projection matrix by multiplying it by the stage structure many times. This is actually one way of performing eigenanalysis, and it is called the {\em \index{lambda!power iteration}power iteration method}. It is not terribly efficient, but it works well in some specific applications. (This method is \emph{not} used by modern computational languages such as \R.) The population size will grow toward infinity, or shrink toward zero, so we keep rescaling $N$, dividing the stages by the total $N$, just to keep things manageable. Let $t$ be big, and rescale $N$. <<>>= t <- 20 Nt <- N0/sum(N0) @ We then create a for-loop that re-uses $N_t$ for each time step, making sure we have an empty numeric vector to hold the output. <<>>= R.t <- numeric(t) for(i in 1:t) R.t[i] <- {Nt1 <- A%*%Nt R <- sum(Nt1)/sum(Nt) Nt <- Nt1/sum(Nt1) R } @ Let's compare the result directly to the point estimate of $\lambda_1$ (Fig. \ref{fig:converge}). <>= par(mar=c(5,4,3,2)) plot(1:t, R.t, type='b', main=quote("Convergence Toward "*lambda)) points(t, L1, pch=19, cex=1.5) @ } \end{boxedminipage} \medskip \begin{figure}[ht] \centering \includegraphics[width=0.5\textwidth]{Converge.pdf} \caption{Iterating the population, and recalculating $R_t=N_{t+1}/N_t$ at each time step converges eventually at the dominant eigenvalue, indicated as a solid point. It is possible to use the same power iteration method to get the other eigenvalues, but it is not worth the trouble.} \label{fig:converge} \end{figure} @ \subsection{Stable stage distribution} The relative abundance of the different life history stages is called the \emph{stage distribution}, that is, the distribution of individuals among the stages. A property of a stage structured population is that, if all the demographic rates (elements of the population projection matrix) remain constant, its stage structure will approach a \emph{\index{stage distribution!stable}stable stage distribution}, a stage distribution in which the \textbf{relative} number of individuals in each stage is constant. Note that a population can grow, so that the absolute number of individuals increases, but the relative abundances of the stages is constant; this is the stable stage distribution. If the population is not actually growing ($\lambda=1$) and demographic parameters remain constant, then the population is \emph{stationary} and will achieve a \emph{\index{stage distribution!stationary}stationary stage distribution}, where neither absolute nor relative abundances change. How do we find the stable stage distribution? It also turns out that $w_1$, which is the corresponding eigenvector of $\lambda_1$ (eq. \eqref{eigsol}), provides the necessary information. We scale the eigenvector $w_1$ by the sum of its elements because we are interested in the \emph{distribution}, where all the stages should sum to one. \footnote{Eigenvectors can only be specified up to a constant, arbitrary multiplier.} Therefore the stable stage distribution is \begin{equation} \label{eq:3} SSD = \frac{w_1}{\sum_{i=1}^S{w_1}} \end{equation} where $S$ is the number of stages. Once a population reaches its stable stage distribution it grows exponentially, \begin{gather*} \mathbf{N_t = A^{t} N_0}\\ N_t = \lambda^t N_0 \end{gather*} represented either in the matrix notation (for all stages), or simple scalar notation (for total $N$ only). \medskip \noindent \begin{boxedminipage}{\linewidth} {\footnotesize \paragraph{Calculating the stable stage distribution} The dominant eigenvector, \vec{w}, is in the same position as the dominant eigenvalue. We extract \vec{w}, keeping just the real part, and divide it by its sum to get the stable stage distribution. <<>>= w <- Re(eigs.A[["vectors"]][,dom.pos]) ssd <- w/sum(w) round(ssd, 3) @ This shows us that if the projection matrix does not change over time, the population will eventually be composed of 80\% seeds, 13\% small adults, and 7\% large adults. Iterating the population projection will also eventually provide the stable stage distribution (e.g., Fig. \ref{fig:MallGrow}). } \end{boxedminipage} \medskip \subsection{Reproductive value} \label{sec:repval} If the stage structure gives us one measure of the importance of a stage (its abundance), then the \emph{reproductive value} gives us one measure of the importance of an \emph{individual} in each stage. Reproductive value is the expected contribution of each individual to present and future reproduction. We characterize all individuals in a stage using the same expected \index{reproductive value}reproductive value. We find each stage's reproductive value by solving for the dominant \emph{left} eigenvector $\vec{v}$, where \begin{equation} \label{eq:4} \vec{vA}=\lambda\vec{v} \end{equation} Like the relation between the dominant right eigenvector and the stable stage distribution, this vector is actually \emph{proportional} to the reproductive values. We typically scale it for $v_{0}=1$, so that all reproductive values are relative to that of the first stage class (e.g. newborns or seeds). \begin{equation} \label{eq:5} RV = \frac{v_1}{\sum_{i=1}^S{v_1}} \end{equation} \medskip \noindent \begin{boxedminipage}{\linewidth} {\footnotesize \paragraph{Calculating reproductive value} We get the left eigenvalues and -vectors by performing eigenanalysis on the transpose of the projection matrix. The positions of the dominant right and left eigenvalues are the same, and typically they are the first. We perform eigenanalysis, extracting just the the dominant left eigenvector; we then scale it, so the stage 1 has a reproductive value of 1.0. <<>>= M <- eigen(t(A)) v <- Re(M$vectors[,which.max(Re(M$values))]) RV <- v/v[1] RV @ Here we see a common pattern, that reproductive value, $v$, increases with age. In general, reproductive value of individuals in a stage increases with increasing probability of reaching fecund stages. } \end{boxedminipage} \medskip \subsection{Sensitivity and elasticity} \index{sensitivity}Sensitivity and elasticity tell us the relative importance of each transition (i.e. each arrow of the life cycle graph or element of the matrix) in determining $\lambda$. They do so by combining information on the stable stage structure and reproductive values. The stage structure and reproductive values each in their own way contribute to the importance of each stage in determining $\lambda$. The stable stage distribution provides the relative abundance of individuals in each stage. Reproductive value provides the contribution to future population growth of individuals in each stage. Sensitivity and elasticity combine these to tell us the relative importance of each transition in determining $\lambda$. Sensitivities of a population projection matrix are the direct contributions of each transition to determining $\lambda$. We would say, speaking in more mathematical terms, that the sensitivities for the elements $a_{ij}$ of a projection matrix are the changes in $\lambda$, given small changes in each element, or $\delta \lambda / \delta a_{ij}$. Not surprisingly, then, these are derived from the stable stage distribution and the reproductive values. Specifically, the sensitivities are calculated as \begin{equation} \label{eq:sens} \frac{\delta \lambda}{\delta a_{ij}}=\frac{v_{ij}w_{ij}}{\vec{v}\cdot\vec{w}} \end{equation} where $v_{i}w_{j}$ is the product of each pairwise combination of elements of the dominant left and right eigenvectors, \vec{v} and \vec{w}. The \emph{dot product}, $\vec{v} \cdot \vec{w}$, is the sum of the pairwise products of each vector element. Dividing by this causes the sensitivities to be relative to the magnitudes of $\vec{v}$ and $\vec{w}$. \medskip \noindent \begin{boxedminipage}{\linewidth} {\footnotesize \paragraph{Sensitivity of projection matrices} Let's calculate sensitivities now. First we calculate the numerator for eq. \ref{eq:sens}. <<>>= vw.s <- v %*% t(w) @ Now we sum these to get the denominator, and then divide to get the sensitivities. (The dot product $\vec{v} \cdot \vec{w}$ yields a $1 \times 1$ matrix; in order to divide by this quantity, the simplest thing is to cause the dot product to be a simple scalar rather than a matrix (using {\tt as.numeric}), and then \R~will multiply each element.) <<>>= (S <- vw.s/as.numeric(v%*%w)) @ We see from this that the most important transition exhibited by the plant is $s_{21}$, surviving from the seed stage to the second stage (the element $s_{31}$ is larger, but is not a transition that the plant undergoes). } \end{boxedminipage} \medskip \index{elasticity}Elasticities are sensitivities, weighted by the transition probabilities. Sensitivities are large when reproductive value and or the stable age distribution are high, and this makes sense biologically because these factors contribute a lot to $\lambda$. We may, however, be interested in how a \emph{proportional} change in a transition element influences lambda---how does a 10\% increase in seed production, or a 25\% decline in juvenile survival influence $\lambda$? For these answers, we need to adjust sensitivities to account for the relative magnitudes of the transition elements, and this provides the elasticities, $e_{ij}$, where \begin{equation} e_{ij}=\frac{a_{ij}}{\lambda}\frac{\delta \lambda}{\delta a_{ij}}. \end{equation} \medskip \noindent \begin{boxedminipage}{\linewidth} {\footnotesize \paragraph{Elasticity of projection matrices} In \R, this is also easy. <<>>= elas <- (A/L1) * S round(elas, 3) @ } \end{boxedminipage} \medskip Note that all the elasticities except the seed production by small adults appear equally important. Specifically, the same proportional change in any of these elements will result in approximately the same change in $\lambda$. There are two nice features of elasticities. First, \emph{impossible transitions have elasticities equal to zero}, because we multiply by the projection matrix itself. Second, the \emph{elasticities sum to zero}, and so it is easier to compare elasticities among differ matrices and different organisms. Once we have the sensitivities and elasticities, we can really begin to see what is controlling the growth rate of a stage (or age) structured population. Although these values do not tell us which stages and transitions will, in reality, be influenced by natural phenomona or management practices, they provide us with the predicted \emph{effects} on $\lambda$ of a proportional change in a demographic rate, $P$ or $F$. This is particularly important in the management of invasive (or endangered) species where we seek to have the maximum impact for the minimum amount of effort and resources \cite{Caswell:2001gu, Ellner:2006qe}. \subsection{More demographic model details} \paragraph{Births} For demographic models, a ``birth'' \index{birth}is merely the appearance in the first stage. If we census birds, a ``birth'' might be a fledging, if this is the youngest age class we sampled. If we census plants, we might choose to count seeds as the first age class, or we might use seedling, or some size threshold as the first stage. Regardless of the first stage- or age-class we use, a birth is the first appearance of an individual in the first stage. \paragraph{Pre- $vs$. post-breeding census} Note that you are sampling the population of mallwort at a particular time of year. This sampling happens to be a \emph{\index{postbreeding census}postbreeding census} because you captured everything right after breeding, when progeny were observed directly. The projection matrix would look different, and the interpretation of the matrix elements would differ, if we had used a \emph{\index{prebreeding census}prebreeding census}, sampling the population before breeding. In particular, the projection matrix would have only two stages (small and large adults), because no seeds would be present at the time of sampling. The contribution of adults to the youngest stage, therefore, would represent both fertility and survival to the juvenile stage in late spring. Nonetheless, both models would be equivalent, generating the same $\lambda$. \paragraph{Birth pulse $vs$. birth flow model} Another assumption we are making is that individuals set seed, or give birth, all at once. We refer to the relevant model as a \emph{\index{birth-pulse}birth-pulse model}. On the other hand, if we assume that we have continuous reproduction, we do things quite differently, and would refer to this as a \emph{\index{birth-flow}birth-flow model}. Whether a population is breeding continuously over a year, or whether reproduction is seasonal, will influence how we estimate fecundities. Even for synchronously breeding populations, many models pool years into a single age class or stage. As result, we need to be careful about how we approximate probabilities that will differ among individuals within the age- or stage-class. These details can get very confusing, and smart people don't always get it right. Therefore, consult an expert \cite{Caswell:2001gu,Ellner:2006qe}, and remember that the stages of life cycle graph and matrix are the stages that you collect at one point in time. \section{Confronting Demographic Models with Data} This section uses \R~extensively throughout. It is common to create a demographic matrix model with real data, and then use that model for an applied purpose (e.g.,\cite{Crouse1987, Endress:2004bs}). A central question, however, is just how confident we can be in our model, and the values we derive from it. It turns out that we can use our data to derive confidence intervals on important parameters. In Chapter 1, we used resampling to draw observed annual changes in bird counts at random to generate growth trajectories and confidence intervals on population size. Here we resample raw data to find confidence limits on $\lambda$. The method used here, \emph{\index{bootstrapping}bootstrapping}, and related data-based inference techniques have a large literature \cite{Manly1997}. Davison and Hinkley \cite{Davison1997} have an comprehensive \R-based text. Such randomization methods are very useful for a wide range of models in ecology, where the data do not conform clearly to parametric distributions or to situations like demographic models \cite{McPeek:1993rm} or null models \cite{Gotelli1996} for which analytical approximations are difficult or not possible. The basic idea of bootstrapping is to \begin{enumerate} \item calculate the \emph{observed parameter(s)} of interest (e.g., a mean, or $\lambda$) with your model and data, \item resample your data \emph{with replacement} to create a large number of datasets and recalculate your parameter(s) for each resampled dataset to generate a \emph{distribution} of the \emph{bootstrapped}% \footnote{``Bootstrapped'' estimates are thus named because you are picking yourself up by your own bootstraps -- a seemingly impossible task.}% parameter(s), \item Calculate a confidence interval for the bootstrapped parameter values --- this will provide an estimate of the confidence you have in your observed parameter. This will provide an {\em empirical confidence interval}. \end{enumerate} \subsection{An Example: \emph{Chamaedorea} palm demography} \emph{\index{Chamaedorea@\emph{Chamaedorea}}Chamaedorea radicalis} Mart. (Arecaceae) is an forest understory palm of northern Mexico, and it is one of approximately 100 \emph{Chamaedorea} species, many of which are economically valuable as either small, shade-tolerant potting plants or as harvested leaves in floral arrangements. Its demography is interesting for a number of reasons, including both management and as an example of a population that appears to be maintained through source-sink dynamics \cite{Berry:2008}. Berry {\em et al}. modeled \emph{Chamaedorea} demography with five stages (Fig. \ref{fig:Cham}). Demography is also influenced by substrate type, by livestock browsing, and harvesting \cite{Berry:2008}. Here we use a subset of the data to illustrate the generation of demographic parameters and confidence intervals. \index{harvesting, palm} \begin{figure}[ht] \centering \includegraphics[width=.75\textwidth]{ChamaedoreaLCG} \caption{Life cycle graph for \emph{Chamaedorea radicalis}. Classification criteria are based on the number of leaflets on the youngest fully-expanded leaf. Life-history stage transitions are indicated by arrows with solid lines and reproduction is indicated by dashed lines. Abbreviations: S-seed, Ss-seedling (bifid leaves), J-Juvenile (3--9 leaflets), A1-small adult (10--24 leaflets), A2-large adult (> 24 leaflets). Source: \cite{Endress:2004bs,Berry:2008}} \label{fig:Cham} \end{figure} This study was conducted in the montane mesophyll forests of Sierra de Guatemala mountain range, near the communities of San Jos\'{e} and Alta Cimas within the \index{El Cielo Biosphere Reserve}El Cielo Biosphere Reserve, Tamaulipas, Mexico (22$^\circ$55'--23$^\circ$30'N and 99$^\circ$02'--99$^\circ$30'W). Villagers within El Cielo (palmilleros) harvest adult \emph{C. radicalis} leaves for sale to international cut-foliage markets. Harvested leaves are usually $>=$ 40 cm in length, and have minimal damage from insects or pathogens \cite{Endress:2004bs}. These palm leaves are the only natural resource that these villagers are authorized to harvest, and provide the main source of income for most families. Although \emph{C. radicalis} is dioecious (more complications!), Berry et al. \cite{Berry:2008} used a one sex model, because its simplifying assumptions were well supported with data. Data collected allowed a postbreeding census model with a birth-pulse dynamic. \subsection{Strategy} There are an infinite number of ways to do anything in \R, and I am certain that my approach to this bootstrapping is not the very best way, but it is \emph{useful}. It gives valid answers in a reasonable amount of time, and that is what we want from a model. This is how we proceed in this instance: \begin{enumerate} \item We import the data and look at it. The appearance of the data, how the data are entered for instance, will influence subsequent decisions about how to proceed. \item We extract the relevant data and calculate the projection matrix elements (fecundities and transition probabilities). We first do it all piecewise, to figure out what we are doing. Then we can wrap it up inside a function putting \texttt{funcname <- function(data1, data2, data3)} at the beginning and collecting and returning all relevant parameters at the end (see sec. \ref{sec:writing-your-own} for writing functions). \item We also create a function to generate all the demographic parameters that we will eventually want ($\lambda$, elasticities, etc.). \item Last, we combine these two functions into one that also resamples the original data (with replacement), and then calls the data extraction and calculation functions to generate the new parameters for the bootstrapped data. \item The bootstrapping is repeated $B$ times. \item Having generated $B$ bootstrapped estimates of all the parameters, we can then calculate confidence intervals for any parameter that we like. \end{enumerate} \subsection{Preliminary data management} Let's import the data and have a look at it. For these purposes, we will assume that the data are clean and correct. Obviously, if I were doing this for the first time, data-checking and clean-up would be an important first step. Here we simply load them from the \texttt{primer} package. <<>>= data(stagedat) data(fruitdat) data(seeddat) @ Now I look at the structure of the data to make sure it is at least approximately what I think it is. <<>>= str(stagedat) @ The stage data provide the stage of each individual in the study. Each row is an individual, and its ID number is in column 1. Data in columns 2--4 identify its stage in years 2003--2005. We can count, or tabulate, the number of individuals in each stage in 2004. <<>>= table(stagedat[["Y2004"]]) @ We see, for instance, that in 2004 there were 165 individuals in stage 5. We also see that 17 individuals were dead in 2004 (stage = 0); these were alive in either 2003 or 2005. The fruit data have a different structure. Each row simply identifies the stage of each individual (col 1) and its fertility (number of seeds) for 2004. <<>>= str(fruitdat); @ We can tabulate the numbers of seeds (columns) of each stage (rows). <<>>= table(fruitdat[["Stage"]], fruitdat[["Y2004"]]) @ For instance, of the individuals in stage 4 (row 1), 28 individuals had no seeds, and one inidvidual had 6 seeds. Note also that only stage 4 and 5 had plants with \emph{any} seeds. The seed data are the fates of each seed in a sample of 400 seeds, in a data frame with only one column. <<>>= table(seeddat) @ Seeds may have germinated (2), remained viable (1), or died (0). @ \subsection{Estimating projection matrix} Now we work through the steps to create the projection matrix from individuals tagged in year 2003 and re-censused in 2004. If we convert the life cycle graph (Fig. \ref{fig:Cham}) into a transition matrix. \begin{equation} \label{eq:pm} \left( \begin{array}{ccccc} P_{11} & 0 &0&F_4&F_5\\ P_{21} & P_{22} & P_{23} & 0 & 0\\ 0 & P_{32} & P_{33} & P_{34} & 0\\ 0 & 0 & P_{43} & P_{44} & P_{45}\\ 0 &0&0&P_{54}&P_{55} \end{array} \right) \end{equation} Along the major diagonal (where $i=j$) the $P_{ij}$ represent the probability that a palm stays in the same stage. In the lower off-diagonal ($i>j$) the $P_{ij}$ represent the probability of growth, that an individual grows from stage $j$ into stage $i$. In the upper off-diagonal ($i>= mat1 <- matrix(0, nrow=5, ncol=5) @ \subsubsection{Fertilities} For each stage, we get mean fertility by applying \texttt{mean} to each stage of the 2004 \index{fertility}fertility data. Here {\tt Stage} is a factor and \texttt{tapply} will caculate a mean for each level of the factor. We will assume that half of the seeds are male. Therefore, we divide fertility by 2 to calculate the fertility associated with just the female seeds. <<>>= ferts <- tapply(fruitdat$Y2004,fruitdat$Stage, mean ) / 2 ferts @ These fertilities, $F_4$ and $F_5$, are the transitions from stages 4 and 5 (adults) to stage 1 (seeds). Next we insert the fertilities ({\tt ferts}) into the matrix we established above. <<>>= mat1[1,4] <- ferts[1] mat1[1,5] <- ferts[2] @ \subsubsection{Seed transitions} Now we get the frequencies of each seed fate (die, remain viable but dormant, or germinate), and then divide these by the number of seeds tested (the length of the seed vector); this results in proportions and probabilities. <<>>= seed.freqs <- table(seeddat[,1]) seedfates <- seed.freqs / length(seeddat[,1]); seedfates @ The last of these values is $P_{21}$, the transition from the first stage (seeds) to the stage 2 (seedlings). The second value is the transition of seed dormancy ($P_{1,1}$), that is, the probability that a seed remains a viable seed rather than dying or becoming a seedling. Next we insert the seed transitions into our projection matrix. <<>>= mat1[1,1] <- seedfates[2] mat1[2,1] <- seedfates[3] @ \subsubsection{Vegetative stage transitions} Here we calculate the transition probabilities for the vegetative stages. The pair of for-loops will calculate these transitions and put them into stages 2--5.The functions inside the for-loops (a) subset the data for each stage in 2003, (b) count the total number of individuals in each stage in 2003 (year $j$), (c) sum the number of individuals in each stage in 2004, given each stage for 2003, and then (d) calculate the proportion of each stage in 2003 that shows up in each stage in 2004. <<>>= for(i in 2:5) { for(j in 2:5) mat1[i, j] <- { x <-subset(stagedat, stagedat$Y2003 == j) jT <- nrow(x) iT <- sum(x$Y2004 == i) iT/jT } } @ @ Here we can see the key parts of a real projection matrix. <<>>= round(mat1,2) @ Compare these probabilities and fertilities to the life cycle graph and its matrix (Fig. \ref{fig:Cham}, eq. \eqref{eq:pm}). The diagonal elements $P_{j,j}$ are stasis probabilities, that an individual remains in that stage. Growth, from one stage to the next, is the lower off-diagonal, $P_{j+1,j}$. Regression, moving back one stage, is the upper off diagonal, $P_{j-1,j}$. The fertilities are in the top row, in columns 4 and 5. Note that there is a transition element in our data that is not in eq. \eqref{eq:pm}: $P_{53}$. This corresponds to very rapid growth --- a real event, albeit highly unusual. \subsubsection{A function for all transitions} What a lot of work! The beauty, of course, is that we can put all of those lines of code into a single function, called, for instance, {\tt ProjMat}, and all we have to supply are the three data sets. You could examine this function by typing \texttt{ProjMat} on the command line, with no parentheses, to see the code and compare it to our code above. You code also try it with data. <>= #source("DemoFunc.R") <>= ProjMat(stagedat, fruitdat, seeddat) @ This provides the observed transition matrix (results not shown). \subsection{Eigenanalyses} Next we want to do all those eigenanalyses and manipulations that gave us $\lambda$, the stable age distribution,reproductive value, and the sensitivity and elasticity matrices. All of this code is wrapped up in the function {\tt DemoInfo}. Convince yourself it is the same code by typing {\tt DemoInfo} with no parentheses at the prompt. Here we try it out on the projection matrix we created above, and examine the components of the output. <<>>= str( DemoInfo(mat1) ) @ We find that \texttt{DemoInfo} returns a \emph{list} with six named \emph{components}. The first component is a scalar, the second two are numeric vectors, and the last three are numeric matrices. The last of these is the projection matrix itself; it is often useful to return that to prove to ourselves that we analyzed the matrix we intended to. \subsection{Bootstrapping a demographic matrix} All of the above was incredibly useful and provides the best estimates of most or all the parameters we might want. However, it does not provide any idea of the certainty of those parameters. By bootstrapping these estimates by resampling our data, we get an idea of the uncertainty. Here we work through the steps of resampling our data, as we build a function, step by step, inch by inch. The basic idea of resampling is that we assume that our sample data are the best available approximation of the entire population. Therefore, we draw, with replacement, new data sets from the original one. See the last section in Chapter 1 for ideas regarding simulations and bootstrapping. We will create new resampled (bootstrapped) data sets, where the rows of the original data sets are selected at random with replacement. We then apply {\tt ProjMat} and {\tt DemoInfo}. The first step is to get the number of observations in the original data. <<>>= nL <- nrow(stagedat) nF <- nrow(fruitdat) nS <- nrow(seeddat) @ With these numbers, we will be able to resample our original data sets getting the correct number of resampled observations. Now we are going to use {\tt lapply} to perform everything multiple times. By ``everything,'' I mean \begin{enumerate} \item resample the observations to get bootstrapped data sets for vegetative stages, seed fates, and fertilities, \item calculate the projection matrix based on the three bootstrapped data sets, \item perform eigenanalysis and calculate $\lambda$, stage structure, sensitivities, and elasticities. \end{enumerate} All of that is one replicate simulation, $n=1$. For now, let's say $n=5$ times as a trial. Eventually this step is the one we will ask \R~to do 1000 or more times. <<>>= n <- 5 @ Next we use {\tt lapply} to do \emph{everything}, that is, a replicate simulation, $n$ times. It will store the $n$ replicates in a {\em list}, $n$ {\em components} long. Each of the $n$ components will be the output of {\tt DemoInfo}, which is itself a list. <<>>= n <- 5 out <- lapply(1:n, function(i) { stageR <- stagedat[sample(1:nL, nL, replace=TRUE),] fruitR <- fruitdat[sample(1:nF, nF, replace=TRUE),] seedR <- as.data.frame(seeddat[sample(1:nS,nS, replace=TRUE),]) matR <- ProjMat(stagedat=stageR, fruitdat=fruitR, seeddat=seedR) DemoInfo( matR ) } ) @ This code above uses {\tt sample} to draw row numbers at random and with replacement to create random draws of data ({\tt stageR, fruitR, and seedR}). We then use {\tt ProjMat} to generate the projection matrix with the random data, and use {\tt DemoInfo} to perform all the eigenanalysis and demographic calculations. Let's look at a small subset of this output, just the five $\lambda$ generated from five different bootstrapped data sets. The object \texttt{out} is a list, so using \texttt{sapply} on it will do the same thing to each component of the list. In this case, that something is to merely extract the bootstrapped $\lambda$. <<>>= sapply(out, function(x) x$lambda) @ We see that we have five different estimates of $\lambda$, each the dominant eigenvalue of a projection matrix calculated from bootstrapped data. We now have all the functions we need to analyze these demographic data. I have put all these together in a function called \texttt{DemoBoot}, whose arguments (inputs) are the raw data, and $n$, the number of bootstrapped samples. <<>>= args(DemoBoot) @ \subsection{The demographic analysis} Now we are armed with everything we need, including estimates and means to evaluate uncertainty, and we can move on to the ecology. We first interpret point estimates of of demographic information, including $\lambda$ and elasticities. Then we ask whether $\lambda$ differs significantly from 1.0 using our \index{bootstrapped confidence interval}bootstrapped confidence interval. First, point estimates based on original data. <<>>= estims <- DemoInfo( ProjMat(stagedat,fruitdat,seeddat)) estims$lambda @ Our estimate of $\lambda$ is greater than one, so the population seems to be growing. Which transitions seem to be the important ones? <<>>= round( estims$Elasticities, 4) @ It appears that the most important transition is persistence in the largest adult stage ($a_{5,5}=0.3$). Specifically, proportional changes to the persistence in this stage, neither regressing nor dying, are predicted to have the largest postive effect on the lambda of this population. We stated above that the population appears to be growing. However, this was based on a sample of the population, and not the entire population. One way to make inferences about the population is to ask whether the confidence interval for $\lambda$ lies above 1.0. Let's use {\tt DemoBoot} to bootstrap our confidence interval for $\lambda$.\footnote{The number of replicates needed for a bootstrap depend in part on how close the interval is to critical points. If, for instance, your empirical $P$-value seems to be very close to your cutoff of $\alpha=0.05$, then you should increase the replicates to be sure of your result. These days $n=1000$ is considered a bare minimum.} First, we'll run the bootstrap, and plot the $\lambda$'s. <<>>= system.time( out.boot <- DemoBoot(stagedat,fruitdat,seeddat,n=1000) ) lambdas <- sapply(out.boot, function(out) out$lambda) <>= myQ(mr=c(5,4,3,2)) <<>>= hist(lambdas, prob=T) lines(density(lambdas)) <>= dev.print(pdf, "lambdaCI.pdf") @ \begin{figure}[ht] \centering \includegraphics[width=0.5\textwidth]{lambdaCI} \caption{The frequency distribution for our bootstrapped $\lambda$. Note that it is fairly symmetrical, and largely greater than 1.0.} \label{fig:LCI} \end{figure} From this it seems clear that the population is probably growing ($\lambda > 1.0$), because the lower limit of the histogram is relatively large (Fig. \ref{fig:LCI}). We need to get a real confidence interval, however. Here we decide on a conventional $\alpha$ and then calculate quantiles, which will provide the median (the 50th percentile), and the lower and upper limits to the 95\% confidence interval.\footnote{Quantiles are ordered points in a cumulative probability distribution function.} <<>>= alpha <- 0.05 quantile(lambdas, c(alpha/2, .5, 1-alpha/2)) @ From this we see that the 95\% confidence interval (i.e. the 0.025 and 0.975 quantiles) does not include 1.0. Therefore, we conclude that under the conditions experienced by this population in this year, this \emph{Chamaedorea} population, from which we drew a sample, could be characterized as having a long-term asymptotic growth rate, $\lambda$, that is greaater than 1.0, and therefore would be likely to increase in abundance, if the environment remains the same. \paragraph{A caveat and refinement} Bootstrapping as we have done above, known variously as the basic or percentile bootstrap, is not a cure-all, and it can give inappropriate estimation and inferrence under some circumstances. A number of refinements have been proposed that make bootstrapping a more precise and accurate procedure \cite{Davison1997}. The problems are worst when the data are such that the bootstrap replicates are highly skewed, so that the mean and median are quite different. When the data are relatively symmetric, as ours is (Fig. \ref{fig:LCI}), the inference is relatively reliable. Often skewness will cause the mean of the bootstrap samples to differ from our observed estimate, and we refer to this as \emph{bias}. We should adjust the bootstrapped samples for this bias \cite{McPeek:1993rm}. Here we calculate the bias. <<>>= bias <- mean(lambdas)-estims$lambda bias @ We find that the bias is very small; this gives us confidence the our confidence intervals are pretty good. Nonetheless, we can be thorough and correct our samples for this bias. We subtract the bias from the bootstrapped $\lambda$ to get our confidence interval. <<>>= quantile(lambdas-bias, c(alpha/2, .5, 1-alpha/2)) @ These \index{bias-corrected quantiles}bias-corrected quantiles also indicate that this population in this year can be characterized by a $\lambda > 1$. If we want to infer something about the future success of this population, we need to make additional assumptions. First, we must assume that our sample was representative of the population; we have every reason to expect it is. Second, we need to assume that this year was representative of other years. In particular, we need to assume that the weather, the harvest intensity, and the browsing intensity are all representative. Clearly, it would be nice to repeat this for other years, and to try to get other sources of information regarding these factors. \section{Summary} Demography is the study of structured populations. Structure may be described by age or stage, and is represented by life cycle graphs and a corresponding projection or transition matrix of transition probabilities and fertilities. The finite rate of increase, $\lambda$, and the stable stage/age distribution are key characteristics of a population, and are estimated using eigenalysis; populations will grow geometrically at the per capita rate of $\lambda$ only when the population has reached its stable stage/age distribution. We measure the importance of transition elements with sensitivities and elasticities, the absolute or relative contributions $\lambda$ of transition elements. Demographic information is frequently useful for endangered and invasive species. % Problems or Exercises should be sorted chapterwise \section*{Problems} \addcontentsline{toc}{section}{Problems} % % Use the following environment. % Don't forget to label each problem; % the label is needed for the solutions' environment \begin{prob} \label{prob1} \textbf{Demographic analysis of a plant population}\\ Goldenseal (\emph{Hydrastis canadensis}) is a wild plant with medicinal properties that is widely harvested in eastern North American. Its rhizome (the thick underground stem) is dug up, and so harvesting can and frequently does have serious negative impacts on populations. A particular population of goldenseal is tracked over several years and investigators find, tag, and monitor several sizes of individuals \cite{Gorchov:2002uq}. After several years of surveys, they identify six relevant stages: dormant seed, seedling, small 1-leaved plant, medium 1-leaved plant, large 1-leaved plant, fertile plant (flowering, with 2 leaves). They determine that the population project matrix is: \begin{equation} \label{prob1:A} \vec{A}= \left(\begin{array}{cccccc} 0 & 0 & 0 & 0 & 0 & 1.642\\ 0.098 & 0 & 0 & 0 & 0 & 0.437 \\ 0 & 0.342& 0.591& 0.050& 0.095& 0\\ 0 &0.026& 0.295& 0.774& 0.177& 0.194\\ 0& 0& 0& 0.145& 0.596& 0.362\\ 0& 0& 0& 0.016 &0.277 &0.489\\ \end{array} \right) \end{equation} (a) Draw a life cycle graph of this population of goldenseal. Include the matrix elements associated with each transition.\\ (b) Start with $\vec{N} = \left(0\,10\,10\,10\,10\,10\right)$ and graph population dynamics for all stages for 10 years.\\ (c) Determine the stable stage distribution.\\ (d) Determine $\lambda$. Explain what this tells us about the population, including any assumptions regarding the stable stage distribution.\\ (d) Determine the elasticities. Which transition(s) are most influential in determining growth rate? \\ (e) Discuss which stages might be most suitable for harvesting; consider this question from both a financial and ecological perspective.\\ \end{prob} \begin{prob} \label{prob3} \textbf{Demographic analysis of an animal population}\\ Crouse et al. \cite{Crouse1987} performed a demographic analysis of an endangered sea turtle species, the loggerhead (\emph{Caretta caretta}). Management of loggerhead populations seemed essential for their long term survival, and a popular management strategy had been and still is to protect nesting females, eggs, and hatchlings. The ground breaking work by Crouse\footnote{Crouse was a graduate student at the time --- graduate students are the life-blood of modern science, doing cutting edge work and pushing their fields forward.} and her colleagues compiled data to create a stage-based projection matrix to analyze quantitatively which stages are important and least important in influencing long-term growth rate. This work led to US Federal laws requiring that US shrimp fishermen use nets that include Turtle Excluder Devices (TEDs, http://www.nmfs.noaa.gov/pr/species/turtles/ teds.htm ). Crouse et al. determined the transition matrix for their loggerhead populations: \begin{equation} \label{prob2:A} \vec{A} = \left(\begin{array}{ccccccc} 0& 0& 0& 0& 127& 4& 80\\ 0.6747& 0.7370& 0& 0& 0& 0&0\\ 0& 0.0486& 0.6610& 0& 0& 0& 0\\ 0& 0& 0.0147& 0.6907& 0& 0& 0\\ 0& 0& 0& 0.0518& 0& 0& 0\\ 0& 0& 0& 0& 0.8091& 0& 0\\ 0& 0& 0& 0& 0& 0.8091& 0.8089\\ \end{array} \right) \end{equation} (a) Draw a life cycle graph of this loggerhead population. Include the matrix elements associated with each transition.\\ (b) Determine the stable stage distribution.\\ (c) Determine $\lambda$. Explain what this tells us about the population, including any assumptions regarding the stable stage distribution.\\ (d) Determine the elasticities. Which transition(s) are most influential in determining growth rate? \\ (e) What is the predicted long-term relative abundance of all stages? What do we call this?\\ (f) If your interest is to maximize long-term growth rate, in which stage(s) should you invest protection measures? Which stages are least likely to enhance long-term growth rate, regardless of protective measures?\\ (g) Start with $\vec{N} = \left(0\,10\,10\,10\,10\,10\right)$ and graph dynamics for all stages for 10 years.\\ \end{prob}