Lifetime reproductive output: individual stochasticity, variance, and sensitivity analysis

Lifetime reproductive output (LRO) determines per-generation growth rates, establishes criteria for population growth or decline, and is an important component of fitness. Empirical measurements of LRO reveal high variance among individuals. This variance may result from genuine heterogeneity in individual properties, or from individual stochasticity, the outcome of probabilistic demographic events during the life cycle. To evaluate the extent of individual stochasticity requires the calculation of the statistics of LRO from a demographic model. Mean LRO is routinely calculated (as the net reproductive rate), but the calculation of variances has only recently received attention. Here, we present a complete, exact, analytical, closed-form solution for all the moments of LRO, for age- and stage-classified populations. Previous studies have relied on simulation, iterative solutions, or closed-form analytical solutions that capture only part of the sources of variance. We also present the sensitivity and elasticity of all of the statistics of LRO to parameters defining survival, stage transitions, and (st)age-specific fertility. Selection can operate on variance in LRO only if the variance results from genetic heterogeneity. The potential opportunity for selection is quantified by Crow’s index \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\mathcal {I}$\end{document}I, the ratio of the variance to the square of the mean. But variance due to individual stochasticity is only an apparent opportunity for selection. In a comparison of a range of age-classified models for human populations, we find that proportional increases in mortality have very small effects on the mean and variance of LRO, but large positive effects on \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\mathcal {I}$\end{document}I. Proportional increases in fertility increase both the mean and variance of LRO, but reduce \documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$\mathcal {I}$\end{document}I. For a size-classified tree population, the elasticity of both mean and variance of LRO to stage-specific mortality are negative; the elasticities to stage-specific fertility are positive.


Introduction
Like all men in Babylon, I have been proconsul; like all, I have been a slave. I have known omnipotence, ignominy, imprisonment . . . I owe this almost atrocious variety to an institution which other republics know nothing about, and which operates among them imperfectly and in secret: the lottery.

Jorge Luis Borges, The Lottery in Babylon
Lifetime reproductive output (LRO) is, as the name implies, the total production of offspring over the lifetime of an individual 1 and is one of the most important characteristics of an individual life history. The expectation of LRO, calculated in terms of female offspring per female, is the net reproductive rate R 0 . In ecology, the critical value R 0 = 1 defines the boundary separating population growth and persistence from population decline and extinction. In epidemiology, R 0 for a pathogen determines whether a disease will or will not cause an outbreak. In evolutionary biology, R 0 is a critical component of fitness (sometimes considered to be fitness, although that is sometimes problematic).
Genetic variance in LRO is considered to be the raw material on which natural selection operates. Crow (1958) introduced an index, which now bears his name, of the "opportunity for selection." If X denotes some measure of fitness, then the opportunity for selection is also known as the standardized variance. It gives the maximum rate of evolutionary change that could be produced by selection if all the variance in X were genetic. To rigorously investigate opportunities for selection would require more explicit population genetic models. Even so, the opportunity for selection is widely used in studies of both animal and human populations (e.g., Jones 2009;Brown et al. 2009;Robbins et al. 2011;Moorad et al. 2011;Courtiol et al. 2012). Empirical studies of individual LRO routinely find large variance and usually a positive skew. Most individuals produce few, or no, offspring, while a few rare individuals produce many offspring (Clutton-Brock 1988;Newton 1989). These differences in LRO have two possible sources. One is heterogeneity-differences in the properties of individuals-including genetic heterogeneity, physiological differences, phenotypic plasticity, and environmental heterogeneity.
However, the differences may also be due to individual stochasticity (Caswell 2009(Caswell , 2011(Caswell , 2014van Daalen and Caswell 2015). Individual stochasticity refers to differences among individuals due to the accumulation of random outcomes of the stochastic processes of mortality, growth, development, breeding, etc. Individual stochasticity would lead to variance among individuals even if they were totally identical and experienced exactly the same demographic rates. 2 Depending on the outcome under consideration, variance due to individual stochasticity can be as great as or even exceed that caused by unobserved heterogeneity (Tuljapurkar et al. 2009;Steiner et al. 2010;Caswell 2011Caswell , 2014Hartemink et al. 2017). Thus, before invoking heterogeneity as the explanation for differences among individuals, it is important to calculate the variation due to stochasticity as a kind of "neutral model" for variation (Steiner and Tuljapurkar 2012).
Individual stochasticity itself has two components. Consider an individual at some stage in its life cycle (say, at birth). The growth, development, and eventual death of this individual define a path through the stages of the life cycle. The pathways of two or more identical individuals, subject to the same rates at every stage of the life cycle, will differ randomly among themselves, and this variation among pathways is one source of individual stochasticity. At each step along an individual's path, it may or may not reproduce. If it reproduces, the number of its offspring will be drawn from some probability distribution. This within-pathway stochastic variation in (stage-specific) fertility is the second source of individual stochasticity.
Our goal is to calculate the variance (and other statistics) of LRO, due to individual stochasticity, from basic demographic information, just as R 0 can be calculated from age-classified or stage-classified demographic models (Rhodes 1940;Cushing and Zhou 1994;De Camino-Beck and Lewis 2007;Cushing and Ackleh 2012). To calculate the variance in LRO, we need to account for the generally infinite number of pathways through the life cycle, calculate the probabilities of each path, calculate the distribution of reproductive output at each stage on each path, and then integrate those probabilities to calculate the mean, variance, etc. of LRO.
The calculation of variance in LRO has been approached in several ways. Early studies used simulation to generate random trajectories through stages, including stages defined by reproductive output, to create a sample of lives from which variance in LRO could be calculated (e.g., Tuljapurkar et al. 2009;Steiner et al. 2010). An analytical solution for all the moments of LRO was provided by Caswell (2011) and provides the framework for our results here. That result took the form of an iterative calculation rather than a closed form expression. Steiner and Tuljapurkar (2012) presented a closed form analytical solution for all the moments of one component of LRO, using a moment generating function approach. They also reported simulations of the complete distribution of LRO and explored effects of heterogeneity. However, their solution included only part of the variance in LRO. They assumed the fertility of each age or stage to be a fixed deterministic quantity, neglecting the within-pathway component of variance. For example, the age-specific fertility of a 26-year-old Japanese woman in 1950 was 0.25. The analysis of Steiner and Tuljapurkar (2012) assumes that every 26-year-old woman, without exception, produces one fourth of a baby. In our framework, every 26-year-old woman, without exception, produces one baby with a probability of 0.25 and zero babies with a probability of 0.75 (Caswell 2011). This component of variance is biologically realistic and quantitively important. For example, analysis of 40 developed countries during the second demographic transition found that, as life expectancy increased, the fraction of variance in LRO due to within-trajectory randomness increased from ∼ 50% to ∼ 99% (Van Daalen and Caswell 2015). We show further examples and present the methodology for decomposing variance into these components, below.
All the studies published so far agree in finding that the variance in LRO due to individual stochasticity can comprise a high proportion of the observed variance in LRO, and that to ignore it is to miss a major source of variation (e.g., Tuljapurkar et al. 2009;Steiner et al. 2010;Steiner and Tuljapurkar 2012;Caswell 2011;van Daalen and Caswell 2015).
In this paper, we present exact, closed form, analytical formulae for all the moments of LRO, extending and replacing the iterative formulae of Caswell (2011). We include both the between-trajectory and within-trajectory components of variance and incorporate variance in stagespecific fertility either empirically or by a statistical model. We also present the sensitivity and elasticity analysis of means, variances, and all moments of LRO due to individual stochasticity. Our results can be applied to any age-or stage-classified matrix population model and to constant, periodic, or stochastic environments.
Our results rely on a mathematical model called a Markov chain with rewards, which we describe in "Markov chains with rewards as a model for LRO." The construction of a Markov chain with rewards requires demographic information on survival, stage transitions, and fertility. Section "The statistics of LRO" presents the calculation of all the moments and other descriptive statistics of LRO. Section "Sensitivity analysis of LRO" presents the sensitivity and elasticity analysis of LRO. Section "Examples" presents examples of age-classified and stage-classified populations, and "Discussion" concludes with a discussion of results and possible extensions. Proofs and derivations appear in Appendix A.

Markov chains with rewards as a model for LRO
Notation Matrices are denoted by uppercase boldface letters (e.g., P), and vectors by lowercase boldface letters (e.g., ρ). Vectors are column vectors by default; X T is the transpose of X. The vector 1 n is a n × 1 vector of ones, I n is the identity matrix of order n, and e i is the ith unit vector, with a 1 in the ith entry and zeros elsewhere. The matrix E is a matrix of ones, and E ij is a matrix with a 1 in the (i, j ) entry and zeros elsewhere. The diagonal matrix with the vector x on the diagonal and zeros elsewhere is denoted D(x). The symbol • denotes the Hadamard, or element-by-element product, and ⊗ denotes the Kronecker product. The vec operator vec X stacks the columns of an m × n matrix X into an mn × 1 column vector. The vecpermutation matrix K m,n satisfies vec X T = K m,n vec X. In cases where it will help understanding, we indicate the dimension of displayed matrix expressions.

The life cycle as a Markov chain
Our approach uses a mathematical model called a Markov chain with rewards. These models have a long history in stochastic process theory (e.g., Howard 1960;Puterman 1994;Sheskin 2010) but have only recently been applied to study lifetime reproduction (Caswell 2011;Van Daalen and Caswell 2015).
The individual life cycle is described by a finite-state, discrete-time, absorbing Markov chain; absorbing states represent death. The population projection matrix is written where U contains the transition probabilities for extant individuals and F contains stage-specific fertilities. Let τ be the number of transient (living) states, α the number of absorbing states, and s = τ + α the total number of states. Then the Markov chain transition matrix, including both transient and absorbing states, can be written where the dimensions of the submatrices are noted. The matrix P is column stochastic. We assume that the spectral radius of U is strictly less than 1; thus, any individual eventually dies with probability 1. The probability of dying while in each of the transient states is contained in the mortality matrix M.
Reproductive rewards An individual experiences a sequence of transitions according to the probabilities in P. Associated with each transition (including the transition of remaining in a stage) is a "reward" representing, in our case, reproductive output. The reproductive reward is a random variable with a specified set of moments. Rewards accumulate over the lifetime of the individual; the total accumulation at the time of death is the LRO. We make the reasonable assumption that individuals in the absorbing state stop accumulating rewards. 3 The reproductive rewards associated with each transition are given by a set of matrices R k , k = 1, 2, . . .. The (i, j ) entry of R k is the kth moment of the reproductive output associated with the transition j → i; that is, (4) These matrices R m can be obtained in several ways, discussed in detail in Caswell (2011).
Empirical distribution. The moments can be calculated empirically from data on individual reproduction. Such data are frequently obtained, but typically only mean reproductive output is reported. In the absence of this information, the reward matrices can be modelled by any probability distribution that is determined by its mean, including the following. Bernoulli distribution. When only a single offspring is produced, mean offspring production equals the probability of reproducing. The matrices of second and third moments satisfy Poisson distribution. When multiple offspring are produced, the Poisson distribution describes a situation where every individual has the same chances of producing those offspring. The moments satisfy Fixed rewards. It is possible to eliminate variance in fertility by creating reward matrices where every individual in a given stage produces exactly the mean number of offspring, in which case In most demographic models, mean fertility is stagespecific rather than transition-specific. That is, mean fertility is described by a vector f that typically appears in the first row of the matrix A. An individual in stage i produces, on average, f i offspring per time step, regardless of what transition the individual may make. In this case If parental survival is required for reproduction, then the last row[s] of R 1 should be set to zero.

The statistics of LRO
In this model, every individual is subject to the same rates at any given stage, so there is no heterogeneity. Even so, each individual may experience a different life course. The resulting variation among individuals is due to the individual stochasticity implied by the vital rates in P and the fertility process described by the R i . Our task is to derive the statistics of LRO from this information. Let ρ k be a vector (dimension s×1) whose ith entry is the kth moment of the remaining lifetime reproductive output of an individual starting in state i, i.e., Because we assume that the dead do not reproduce, we know that the entries of ρ k corresponding to absorbing states are zero. Thus, we obtain the complete statistics of LRO if we can find the moment vectorsρ k , whereρ contains the first τ entries of ρ, corresponding to the living stages. This vector is given bỹ where in terms of this matrix, ρ = Z Tρ . We also defineR k , the τ × τ submatrix of R k corresponding to transitions among the transient states: In terms of these quantities, the moments of LRO are given by the following theorem.

Theorem 1 The moment vectorsρ of lifetime accumulated reproductive output arẽ
and, in general, where N = (I τ − U) −1 is the fundamental matrix of the Markov chain.
Proof See Appendix A.1 In terms of the moment vectors provided by Theorem 1, the mean, variance, standard deviation, and coefficient of variation of LRO, and Crow's index I of opportunity for selection, are given by We focus in this paper on statistics obtained from the first and second moments; these statistics, such as the variance, describe variability of the distribution. Calculations of the skewness, as a measure of the shape of the distribution, are given in Caswell (2011) and Van Daalen and Caswell (2015). All of the methods we present here can be applied to these statistics.

Partitioning variance: within and between pathways
Variance in LRO is partly due to individuals following different pathways through the life cycle and partly due to variance in stage-specific fertility along those pathways. The overall variance V (ρ) can be decomposed into these components using the law of conditional variance (Rényi 1970, Theorem 1, p. 275). To calculate the variance between trajectories, we eliminate the within-trajectory variance by calculating V (ρ) using the fixed reward model, in which the second moment matrix R 2 is given by Eq. 8. The within-trajectory variance is obtained by subtraction as

Sensitivity analysis of LRO
The statistics of LRO depend on the life history parameters that determine the transition matrices U and M, and the moment matrices R 1 , R 2 , . . . . Sensitivity analysis reveals how these parameters affect LRO; our goal is to derive the sensitivity and elasticity of LRO to changes in any parameters. Although sensitivity analysis of the net reproductive rate R 0 has been presented before (Matser et al. 2009;Caswell 2009), there has been no such analysis for the variance or other measures of variation in LRO. We derive sensitivity formulae using matrix calculus formalism (Caswell 2007(Caswell , 2008(Caswell , 2009(Caswell , 2012Caswell and Salguero-Gómez 2013). Let θ be a vector (p × 1) of parameters of interest; these could be mortalities, transition probabilities, means or variances of fertility, etc. The sensi-tivity of the moment vectorρ m to the parameter vector θ is the derivative matrix dρ m That is, the (i, j ) entry of this sensitivity matrix is the derivative of the ith entry ofρ m with respect to the j th entry of θ .
To write the sensitivity of theρ m , let us define the following matrices.
Proof Proof is given in Appendix A.2.
The sensitivities of the first three moments of LRO are of particular interest; they are If, as is often the case, the model includes only a single absorbing state (death), then M is a (1 × τ ) matrix, given by In this case, the derivative of P that appears in Eq. 27 is Caswell and Van Daalen (2016). This permits expressing all the effects of θ on the transition matrix P in terms of effects on U (see Appendix A.2).

Sensitivity of the statistics of LRO
The momentsρ i provide the statistics (21)-(24) describing the inter-individual variability of LRO, including the variance, standard deviation, coefficient of variation, and the scaled variance (Crow's index) (Caswell 2011;Van Daalen and Caswell 2015). The sensitivities of these quantities are For derivations, see Appendix A.3.

Elasticity
The derivatives in Theorem 2 and Eqs. 38-41 measure the effect of an additive perturbation of the parameter vector θ . Elasticities, which measure the proportional change resulting from a proportional change in θ , are easily calculated. Let ξ be any quantity calculated from the demographic model. The elasticity of ξ with respect to θ is the matrix As usual, elasticity calculations can be applied only if ξ > 0 and θ ≥ 0.

Some special perturbations
Here we consider some perturbations that are of special interest: the survival and transitions in stage-classified models and sensitivity to means and variances of stage-specific fertility.

Mortality and transitions in stage-classified models
In an age-classified model, the transient matrix U is completely determined by the mortality schedule (see "Ageclassified human populations" for age-classified examples).
In a stage-classified model, however, U depends on both survival and probabilities of transitions among stages. To calculate sensitivities of LRO to survival and transition probabilities, we write where has the survival probabilities σ on the diagonal; i.e., = D(σ ), and G is a matrix of transition probabilities conditional on survival. Differentiating (43) gives Writing D(dσ ) = I • 1 τ dσ T and simplifying lead to and since σ = exp(−μ), From (44), the sensitivity of U to the growth matrix G is However, G is column-stochastic; i.e., its columns must sum to 1. An arbitrary perturbation of G would result in the loss of the column-stochastic property. The only relevant perturbations are those that maintain the column sums. Thus, perturbations of the entry g ij must be compensated for by changes in other elements of column j . Caswell (2001Caswell ( , 2013) presents a method for compensating for perturbation of each entry in a column-stochastic matrix in a way that maintains the proportional structure of the column. In matrix notation, we write where the matrix is a matrix of perturbations that include the compensation for the column sums. The matrix of derivatives of G with respect to is given by Caswell (2013). Let g i denote column i of G; then

Sensitivity to means and variances of fertility
Perturbations of fertility appear in Eqs. 27-29 as derivatives of the reproductive reward matrices R i . When the distributions of stage-specific fertilities are specified by a parametric distribution, the moments may be linked, so that changes in mean fertility also affect the variance (e.g., in the Poisson distribution, the variance is equal to the mean). Sometimes, however, it is of interest to treat the mean and variance of fertility as independent traits and calculate the sensitivity of LRO to the mean, holding the variance fixed, and to the variance, holding the mean fixed. This subtle but important distinction was emphasized by Tuljapurkar et al. (2003) and Haridas and Tuljapurkar (2005) in the context of the elasticity of the stochastic growth rate to the entries of a stochastically varying matrix. One might compute the effect of changing one of the matrix entries (perhaps by a change in energy allocation strategy), recognizing that this would change both the mean and the variance. Or, one might be interested in the effects of variance per se and manipulate the moments to calculate elasticities with respect to the mean and variance independently. As in Eq. 11, suppose that fertility is defined by a fertility vector f with first and second moments f 1 and f 2 , and a variance vector (50) The first and second moment matrices are given by Sensitivity to mean fertility, variance fixed From Eq. 50, it follows that To hold the variance fixed, we require dv = 0, which implies that To evaluate sensitivity to the mean, we set the parameter vector θ = f 1 , subject to Eq. 53, and obtain Substituting these expressions into Eqs. 31 and 32, and then into Eq. 38, we obtain the sensitivity of the mean and variance of LRO to changes in mean fertility, with variance in fertility held constant.
Sensitivity to variance in fertility, mean fixed To hold the mean fixed, we require that df 1 = 0, in which case (52) implies that Now we set the parameter vector θ = v to obtain Substituting these expressions into Eqs. 31 and 32, and then into Eq. 38, gives the sensitivity of the mean and variance of LRO to changes in the variance in fertility, with mean fertility held constant.

A protocol for the analysis of lifetime reproductive output
The results presented to this point provide a protocol for analysis of lifetime reproductive output, applicable to any matrix population model. A stepwise version of this protocol is given in

Examples
This section presents examples of the calculation of the statistical properties of lifetime reproductive output and its subsequent sensitivity analysis, for both age-classified and stage-classified population. In "Age-classified human populations," we analyze a set of age-classified human populations that span a wide range of demographic conditions. In "A stage-classified tree population," we analyze a sizeclassified model for Canadian hemlock (Tsuga canadensis), a coniferous tree.

Age-classified human populations
The transition matrix U for an age-classified model contains survival probabilities on the subdiagonal and zeros Table 1 A step-by-step protocol for analysis of lifetime reproductive output and its sensitivity, from any stage-or age-classified matrix population model 1. Obtain a transition matrix U, perhaps from decomposing a population projection matrix as A = U + F.

Locate reproductive transitions.
(a) If fertility is transition specific, identify the transitions (e.g., to reproductive states). (b) If fertility is stage-specific, extract the vector f from F. 3. Obtain statistical moments of fertility: (a) From empirical measurements of the moments of stage-specific fertility, or (b) From an assumption of Bernoulli [see equation (5) elsewhere. One absorbing state, death, is included, and the mortality matrix M is calculated according to Eq. 34.
We ignore multiple births and treat the entries of the fertility vector f as the probability of producing a single offspring. The offspring production is given by the reward r ij , following a Bernoulli distribution, The matrix R 1 containing the first moment of the reward matrix is built from the fertility vector using Eq. 11. The higher moments of the reward matrix follow from the Bernoulli model of reproduction, in which the higher moments are all equal to the first, as in Eq. 5, so that R 1 = R 2 = R 3 . We present results for nine populations: the Netherlands (1950 and 2009), Sweden (1891 and 2010), Japan (1947 and 2009), two hunter-gather populations (the Ache of subtropical Paraguay (Gurven and Kaplan 2007;Hill and Hurtado 1996) and the Hadza of the Tanzanian savanna Blurton Jones 2011), and the Hutterites of North America. The Netherlands, Sweden, and Japan are typical of developed countries progressing through the demographic transition. The hunter-gatherer populations have higher mortality, lower life expectancy, and higher fertility than the developed countries. The Hutterites, an Anabaptist religious community in the United States and Canada, are known for having the highest total fertility for any known human population (Eaton and Mayer 1953), but are assumed to have experienced mortality rates similar to that of the USA around 1946-1950. Data for the Netherlands, Sweden, and Japan were obtained from the Human Mortality Database (2013) and the Human Fertility Database (2013). Data for the Ache were obtained from Gurven and Kaplan (2007) and Hill and Hurtado (1996), and for the Hadza from Blurton Jones (2011). The fertility and mortality schedules for the Hutterites were taken from Eaton and Mayer (1953).
Variance in LRO Using Theorem 1 and Eqs. 21-24, and following the protocol in Table 1, we computed the mean, variance (both within and between pathways), standard deviation, coefficient of variation, and Crow's I from U, R 1 , and R 2 . We also calculated life expectancy as the column sums of the fundamental matrix. The results are given in Table 2.
Not surprisingly, recent populations in developed countries have higher longevity and lower mean lifetime repro- Table 2 The statistics of lifetime reproductive output for the Netherlands (NLD), Sweden (SWE), and Japan (JPN), with two points in time for each country, two hunter-gatherer populations, the Hadza and the Ache, and a population of high-fertility Hutterites Most of the variance in the hunter-gatherer populations is due to variance among pathways, which in an age-classified model is determined by survival from birth through the reproductive ages. The most recent years in developed countries show extremely low between-pathway variance because almost all women survive through their reproductive years (Van Daalen and Caswell 2015). The Ache and Hadza, with the lowest life expectancy, show very high V between . The Hutterites have a long life expectancy, but their high fertility amplifies the effect of differences in longevity, so V between is similar to that of Japan in 1947 and Sweden in 1891.
The opportunity for selection I varies less (about sevenfold) among these populations than does the variance in LRO (about 13-fold). The Hutterites showed the lowest opportunity for selection of all the populations we included, yet the other high-fertility populations, the Ache and the Hadza, show the highest values for the opportunity for selection. Taking into account the fact that the hunter-gatherer populations have the highest variances in LRO, mostly due to variation in the pathways individuals take through life, we posit that the high opportunity for selection reflects room for improvement in survival rates from birth to reproductive ages.
The variance in LRO documented in Table 2 is calculated on the assumption that every individual experiences the same vital rates at every age and is thus due to individual stochasticity. Crow's I is a measure of the potential relative increase in fitness per generation, but the variance here is stochastic, not genetic, so the opportunity for selection is only apparent, not real.

Sensitivity analysis
The age-classified models for human populations are parameterized by the mortality rate vector μ and the mean fertility vector f. The sensitivity of the statistics of LRO to these parameters, obtained from Theorem 2 and Eqs. 38-41, requires dvec U/dμ T and dvec R i /df T for i = 1, 2. The derivative of U is where Y is an indicator matrix defining the subdiagonal structure of U, and P = exp(−μ) is the vector of survival probabilities. The Bernoulli distribution assumption implies that the derivatives of the first and second (and all other) moments are equal, with (for derivations, see Appendix A.4). The sensitivity of mean LRO, variance in LRO, and Crow's I to mortality and fertility are shown in Fig. 1. Not surprisingly, increased mortality at any age up to the end of the reproductive period reduces mean LRO, while increased fertility increases mean LRO. The mortality effect is greatest for the Hutterites, because their high fertility magnifies the impact of mortality changes, and least for the recent developed countries. The sensitivity of mean LRO to fertility is given by the survivorship function and thus is smallest for the Ache and Hadza. It is highest for recent developed countries in which low mortality means that almost everyone would survive to benefit from an increase in fertility.
The variance in LRO increases with an increase in mortality rate for most populations; the effect is greatest for the Hutterites (Fig. 1). For the developed countries, the sensitivity is positive over most of the reproductive ages. For the Crow's Index NLD 1950NLD 2009SWE 1891SWE 2010JPN 1947JPN 2009 Hadza Ache Hutterite Fig. 1 Sensitivity of mean LRO, variance in LRO, and Crow's index to changes in age-specific mortality (left column) and age-specific fertility (right column) for nine human populations Hadza and the Ache, variance in LRO decreases with higher age-specific mortality rates. Variance in LRO increases with increasing fertility rates for all countries. The Hutterites, Sweden, the Netherlands, and Japan show a reduced sensitivity of variance in LRO to fertility around the reproductive ages. Crow's opportunity for selection I combines both the mean and the variance. Increased mortality during the reproductive period increases I in all the populations. It is most sensitive to mortality in the Ache and Hadza populations and least sensitive in the Hutterites. An increase in fertility reduces I in all the populations. Thus, the net result of environmental changes that affect both mortality and fertility cannot be predicted a priori.
Both mortality and fertility vary widely across ages in these populations, so it may be useful to standardize the responses by calculating elasticities (Fig. 2). The elasticities of the mean and variance of LRO with respect to mortality are generally low, except for effects of infant mortality, especially for the hunter-gatherer populations, in which infant mortality is high. However, the elasticity of Crow's I is large and positive, in fact, the largest of any of the elasticities obtained. The elasticities with respect to fertility are naturally confined to the reproductive ages. Proportional increases in fertility increase the mean and variance of LRO, but reduce the value of Crow's I.

A stage-classified tree population
As an example of a typical stage-classified population, we analyze a model for the Canadian hemlock (T. canadensis L.). Lamar and McGraw (2005) developed a model based on six size classes (from < 5.0 cm dbh (diameter at breast height) to > 42.5 cm dbh). They reported population projection matrices for trees in a low-disturbance plot in Shenandoah National Park in the eastern USA, between 1997 and 1999. We obtained U and F from the mean popu- Reproduction was measured as new recruits (rather than seeds or seedlings) per individual, per year, as a function of size (Lamar and McGraw 2005). In the absence of information on the empirical distribution of size-specific fertility, we use the Poisson distribution to define the reward moment matrices R 1 and R 2 following Eqs. 11 and 6.

Results
The mean, variance, and other statistics of lifetime reproductive output for T. canadensis, obtained using Theorem 1 and Eqs. 21-24, are shown in Table 3. Because of the high mortality of small trees, the mean LRO is small and life expectancy is only 12 years. However, the variance in LRO is very high, as is the apparent opportunity for selection I. The variance in LRO exceeds that for any of the human populations by two orders of magnitude. The value of I is 600 times higher than the highest for human populations in Table 2. Almost all of the variance in LRO is due to variance among pathways. The variance due to stochastic rewards along those pathways is small and approximately equal to the mean, which reflects the assumption of Poisson distributed fertility. It is swamped by the huge differences among pathways of trees that die small and those that become large and reproductive.  The sensitivity of U to the growth matrix G, maintaining the column sums, is

Sensitivity analysis
where dvec G/dvec T is given by Eq. 49. The sensitivity of LRO to fertility is determined by the derivatives of the R i to fertility. The derivative of the matrix R 1 of mean fertility is given by Eq. 63. The derivatives of R 2 and R 3 with respect to f are calculated from the Poisson distribution using Eqs. 6 and 7: Incorporating (67) and (68) into Eqs. 27-29 and applying Theorem 2 provides the sensitivity of lifetime reproductive output to mortality, fertility, and the growth matrix (see Appendix A.4 for the derivations). To make comparisons across these variables measured on different scales, we calculate the elasticity or proportional sensitivity of LRO using Eq. 42. An increase in mortality in the first size-class or the last size-class reduces both mean lifetime reproduction of Hemlock individuals and variance in LRO (Fig. 3). Increasing fertility rates increases both mean and variance in LRO for trees in the last size-class. These are the trees with the highest survival and fertility rates.
The elasticity of the mean and variance to the growth matrix G is shown in Fig. 4. The results are dominated by the extremely large negative elasticities to the probability of remaining in stage 1 (thereby reducing growth) and the large positive elasticity to the probability of remaining in the largest size class. Increasing the probability of staying in any size-class (again, reducing growth) reduces both mean and variance in LRO.
The mean and variance, the contribution to the variance of different processes, and the sensitivity of these indices to parameters differ between these examples. This reflects the different life history strategies of trees and humans, the difference between between age-classified populations with low fertility, and a size-classified population with strongly size-dependent fertility, and the difference between assumptions of Bernoulli or Poisson distributed rewards. Vastly different life histories can be incorporated into the Markov chain with reward framework, allowing for the investigation of life history in many species, and from different perspectives.

Discussion
Lifetime reproductive output is an outcome of the life cycle. Any demographic model implies a distribution for LRO, just as it implies more familiar measures such as R 0 and life expectancy. Even in a deterministic environment, the LRO is a random variable; the stochasticity arises from two sources: the random pathway that the individual follows through its life and the random fertility it exhibits at each stage on that pathway.
Our results in "Markov chains with rewards as a model for LRO" provide the analytical machinery needed to calculate all the statistical properties of LRO that follow from a specified set of stages, a transition matrix, and the moments of stage-specific fertility. These statistical properties include measures of variability among individuals (the variance, CV, skewness, opportunity for selection, etc.). It is important to recognize that this variability does not reflect heterogeneity, genetic or otherwise. Every individual is subject to the same set of vital rates at any stage in the life cycle. Only the outcomes of applying those rates vary; the variation is thus due to individual stochasticity.
The variance, or standardized variance, calculated from demographic models for a variety of species, is large (see also Caswell 2011;van Daalen and Caswell 2015;Steiner et al. 2010). Individual stochasticity creates a large apparent opportunity for selection that is not, in fact, a true opportunity. In this sense, as emphasized by Steiner and Tuljapurkar (2012), the calculations of individual stochasticity provide a neutral model for LRO. A number of comparisons of calculated variance (due to stochasticity) and observed variance (due to some mix of stochasticity and heterogeneity) have shown that stochasticity may explain a significant amount, or all, of the observed variance (Caswell 2011;Van Daalen and Caswell 2015;Steiner et al. 2010).
It is important to remember what neutral model results do, and do not, imply (Caswell 1976). The results show that a certain amount of variance can be accounted for by stochasticity, and hence, that the mere observation of such variation is no evidence for heterogeneity. It does not prove that the observed variance is stochastic, or that there is no heterogeneity, as pointed out by Steiner and Tuljapurkar The results of Theorem 1 provide an exact solution to the calculation of the statistics of LRO, including both components (within and between pathways; see "Partitioning variance: within and between pathways"). We also provide a complete sensitivity analysis for LRO. Theorem 2 makes it possible to calculate the sensitivity and elasticity of all the moments of LRO and all the statistics calculated from those moments, with respect to changes in mortality, transition probabilities, and the moments of stage-specific fertility.
The results include the sensitivity to changes in mean fertility (holding variance constant) and variance in fertility (holding the mean constant).
The formulae for the moments in Theorem 1 and the sensitivities in Theorem 2 are complicated and opaque, because the relationships between LRO and the life cycle structure, the moments of reproduction, the probabilities of survival, and the infinite diversity of pathways through the life cycle, are complicated. Simplifications that permit qualitative generalities are always welcome, and more work on this will be valuable.
We applied sensitivity analysis to several populations of humans and a population of trees. The patterns of sensitivity and elasticity of LRO that we report for these populations Elasticity to growth have not been described before. Some suggestive patterns appear; they warrant further investigation.
In long-lived age-classified populations with low reproductive output, as diverse as the nineteenth century Swedes, mid-twentieth century Hutterites, the twenty-first century Dutch, and Hadza and Ache hunter-gatherers, the sensitivity of mean LRO to mortality is negative. Most populations show a positive sensitivity of variance in LRO to mortality for the first 40 years of life, only showing a small negative sensitivity between ages 20 and 40. The Ache and Hadza have more pronounced negative sensitivity around these ages, with the Hadza showing negative sensitivity across the first 40 years of life. The sensitivity of Crow's I to mortality in the first 40 years of life is positive, showing that an increase in mortality would increase the apparent opportunity for selection of lifetime reproduction. Patterns for elasticities are similar, though smaller in magnitude.
All populations show broadly similar patterns for the sensitivity of LRO to fertility. The sensitivities of mean LRO and variance in LRO to fertility are positive. The sensitivity of Crow's I to fertility is negative. The elasticity of mean and variance to fertility is positive, but the elasticity of Crow's I is negative. According to these results, populations in modern countries would reduce the apparent opportunity for selection of LRO with increasing fertilities, and high-fertility populations such as Hutterites and huntergatherers would only slightly reduce the opportunity for selection should fertilities increase.
On average, the variance in LRO is 59% between pathways and 41% within pathways. However, the individual populations differ significantly in these contributions by sources of variance. In the twenty-first century in Japan, Sweden, and the Netherlands, only 1.2% of variance in LRO is between pathways and 98.8% is within pathways. Variance in LRO of the high-fertility populations, the Ache, Hadza, and Hutterites, is 71% between pathways and 29% within pathways.
In contrast, in T. canadensis, with size-dependent demography, high fertility, and strongly increasing fertility with size, we find that the elasticities of both mean and variance in LRO to mortality are negative across all size classes. The elasticities of these statistical properties to changes in fertility in the reproductive classes are positive. In case of the elasticity to growth transition rates, mean and variance in LRO once again show similar patterns. Elasticity to stasis, i.e., not growing, is negative for the first five size classes, and positive for the last size class.
Lifetime reproductive output interests population ecologists and epidemiologists (for whom R 0 is a measure of population growth), evolutionary biologists (for whom variance in LRO is a measure of potential selection), and human demographers (for whom declines in LRO following the demographic transition pose serious social policy challenges). It is thus important that our analysis is not restricted to any one class of population models. It applies equally to age-structured, stage-structured, and multistate models, and to any reproductive strategy. It also applies to periodic and stochastic time-varying models, by applying Theorems 1 and 2 to the models in Caswell (2011). The method can also be applied to rewards other than reproductive output, including health status (Caswell and Zarulli 2015) and economic transfers (Caswell and Kluge 2015).
The transition matrix U is part of any population projection matrix; the COMPADRE and COMADRE matrix databases provide many examples (Salguero-Gómez et al. 2016;Salguero-Gómez et al. 2015). The mean reproductive reward matrix R 1 can be obtained from the projection matrix, but the higher moments cannot and require assumptions of a parametric distribution for fertility. We encourage researchers with the appropriate reproduction data to report not only mean fertility but also the higher moments, or even the entire distribution.
Individual stochasticity arises in both reproductive output and in survival or longevity. Our results here complement the analysis of variation in longevity using Markov chain methods, which are widely used in ecology (e.g., Cochran and Ellner 1992, Caswell 2001, 2006, Horvitz and Tuljapurkar 2008, Tuljapurkar and Horvitz 2006 and have well-developed sensitivity analyses. The Markov chain with reward model now permits a similarly deep analysis of lifetime reproduction and its sensitivity analysis. Open Access This article is distributed under the terms of the Creative Commons Attribution 4.0 International License (http:// creativecommons.org/licenses/by/4.0/), which permits unrestricted use, distribution, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons license, and indicate if changes were made.

A.1 Proof of Theorem 1: moments of LRO
We assume a finite state absorbing Markov chain with transition matrix P given by Eq. 3, with τ transient states, α absorbing states, and s = τ + α total states. The (i, j ) entry of the s × s matrix R k is the kth moment of the reproductive reward associated with the transition from state j to state i. We assume that no rewards accrue to individuals in absorbing states, so columns τ + 1 to s of R k are zero.
The ith entry of the s × 1 vector ρ m is the mth moment or remaining lifetime reproductive output for an individual in stage i. Caswell (2011, Proposition 1) derived an iterative approximation for the moment vectors, 2, . . . (69) which converges to the lifetime accumulation, for each stage, as t → ∞. To obtain an analytical solution for the moment vectors, we solve (69) for its equilibrium. The equilibrium moment vector ρ m satisfies Because ρ 0 = 1 s and R 0 is a s × s matrix of ones, this simplifies to Solving (71) directly for ρ m is impossible because the matrix I − P T is singular. However, since no rewards accumulate in the absorbing states, so we need solve only for the vector of remaining LRO from the transient states, which we denote byρ m . Multiplying both sides of Eq. 70 by the matrix Z in Eq. 14 gives Noting that we see thatρ m must satisfỹ Solving (75) forρ m gives (19) and completes the proof of Theorem 1.

A.2 Proof of Theorem 2: sensitivity analysis of moments of LRO
To derive the sensitivity result of Theorem 2, we will break the result (19) into three types of terms and differentiate each of these in turn. We rewrite (19) as Equation 76 contains three types of terms: We differentiate each of these in turn. Differentiating (77) gives Substituting gives and applying the vec operator yields Differentiating term B gives Applying the vec operator to both sides of Eq. 84 yields Differentiating term C gives Applying the vec operator yields dvec C = ρ T ⊗ I K τ τ D(vecR)dvec U + D(vec U)dvecR Next, differentiate (76) and substitute (83)-(87) where appropriate; this leads to Using the terms V m , W k,m−k , and X m identified in Eq. 88, multiplying both sides by N T , and replacing differentials with derivatives with respect to a vector θ of parameters give which completes the derivation of Theorem 2. The term V i in Eq. 89 contains the derivative dvec P/dθ T . Because P is a block structured matrix and 0 and I are constants, the derivative of P can be written as a linear combination of the derivatives of U and M, where (Caswell and van Daalen 2016). In many applications, only a single absorbing state (death) exists, in which case α = 1 and Eqs. 91 and 92 simplify to Eq. 35.

A.3 Sensitivity analysis of the statistics of LRO
The sensitivity analysis of the statistical properties of LRO is obtained by differentiating and applying the vec operator to Eqs. 21-24. Several matrix calculus results are used frequently in these sensitivity calculations. Roth's theorem, i.e., is applied often, as is the rule for taking the vec of a For notational clarity, we redefine some of the statistical properties of LRO as follows in this section only: Variance Differentiating the variance in lifetime reproductive output as given by Eq. 21 gives Applying the vec operator results in which, using the rule for taking the vec of a Hadamard product, yields Taking the derivatives with respect to a vector of parameters θ and rearranging give as found in Eq. 38.

Standard deviation
To obtain the sensitivity of the standard deviation of LRO, we rewrite (22) as and differentiate this giving us Using the vec operator, we obtain where, after replacing the differentials with derivatives with respect to θ and applying the Hadamard-rule on the lefthand side, we find 2D (s) ds Solving for ds/dθ T yields Eq. 39.

Coefficient of variation
The sensitivity of the coefficient of variation is obtained by differentiating from Eq. 23, resulting in However, for any nonsingular matrix Y, Thus, Eq. 109 can be written as Applying the vec operator and Roth's theorem yields Noting that and that D can be rewritten as we obtain Using the rule for taking the vec of a Hadamard product and applying Roth's theorem once more yields Replacing the differentials with derivatives with respect to θ gives Crow's index I The opportunity for selection, also known as Crow's index, can be calculated from the coefficient of variation, as shown in Eq. 24. Differentiating this equation results in dI = dvec (C • C) .
Applying the vec operator to the Hadamard product yields Replacing the differentials with derivatives of Crow's I with respect to θ yields dI dθ T = 2D(C) dC dθ T (120)

A.4 Deriving sensitivities for the examples
To obtain the sensitivities for mean LRO and other statistics, the following pieces are required: In the case of both humans and trees, a single stage of death is incorporated, so that the derivative of P with respect to θ can be obtained from the derivative of U with respect to θ , as shown in Eq. 35. As the transition matrix depends on mortality rates and the reward matrix depends on fertility rates, the required pieces for sensitivity analysis become dvec U dμ T and The way these matrices depend on the underlying vectors of parameters differs between species, as shown below for humans and trees.
Humans In humans, U is a function of the mortality rates μ through its survival probabilities σ . U depends on σ as where Y is a τ × τ matrix with ones on the subdiagonal and zeros elsewhere. Given that the derivative of σ becomes dσ = −D(σ )dμ.
As such, the complete analytical expression for the sensitivity of the transition matrix to the mortality schedule is The sensitivity of the moments of R to a vector of fertility rates depends on the kind of reproductive strategy used. In the case of humans, fertility rate is the chance of producing a single offspring. With the subsequent assumption that the moments of the reward matrix follow a Bernoulli distribution, Therefore, we only need to derive the sensitivity of R 1 to f to know all the moments of the reward matrix. Equation 11 can be rewritten as Differentiating this entails taking the derivative on both sides and applying Roth's theorem. Then, the equation for the sensitivity is simply Trees In trees, the transition matrix consists not only of survival rates but also growth/stasis/regression rates (see Equation 43). By differentiating (43) with respect to , we obtain dvec U dvec T = (I T ⊗ G).
As U now depends on mortality rate through , the sensitivity of the transition matrix to μ is Given that can be written as = I • (1σ T ), differentiating, applying the vec operator, using the rule for taking the vec of a Hadamard product, and applying Roth's theorem give dvec dσ T = D(vec I)(I T ⊗ 1).
The sensitivity of the transition matrix to the growth matrix is found by differentiating (43) with respect to G dvec U dvec T G = ( T ⊗ I).
However, a complication arises from the fact that G is column-stochastic, that is, the probabilities of staying or moving between stages always sums to 1. Sensitivity analysis works on the premise of a small change in the underlying parameters effecting a chance in the model outcome. In the case of G, a small additive change would result in the loss of the column-stochastic property. Fortunately, we can instead use a proportional change. Caswell (2001Caswell ( , 2013 presents a method of proportionally compensating for a perturbation of a single entry in a column-stochastic matrix by proportionally subtracting the perturbation from the other entries in that column. In matrix notation, this becomes In trees, reproduction is assumed to adhere to a Poisson distribution as opposed to a Bernoulli distribution. Higher moments of the reward matrix therefore depend on R 1 as follows: By differentiating the expression for R 2 we obtain Applying the vec operator yields dvec R 2 = dvec R 1 + 2D(vec R 1 )dvec R 1 .
Replacing the differentials with the derivative of R 2 with respect to R 1 gives dvec R 2 dvec R 1 = I s 2 + 2D(vec R 1 ).
Then, the derivative of R 2 with respect to f becomes which, using the expression above and recalling (129), becomes dvec R 2 df T = I s 2 + 2D(vec R 1 ) (Z T ⊗ 1 s ).
Similarly, we can obtain the sensitivity of R 3 to f by differentiating (138), resulting in Applying the vec operator then gives where the differentials can be replaced with derivatives of R 3 with respect to R 1 , yielding dvec R 3 dvec R 1 = I s 2 + 6D(vec R 1 ) + 3D (vec(R 1 • R 1 )) . (146) Given that