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 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 I . Proportional increases in fertility increase both the mean and variance of LRO, but reduce 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 individual1 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, 2011, 2014; van 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 2011, 2014; Hartemink 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 stage-specific 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 is a n × 1 vector of ones, I is the identity matrix of order n, and e 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 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 . 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 m n × 1 column vector. The vec-permutation matrix K satisfies vec X T = K 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 = 1,2,…. The (i,j) entry of R is the kth moment of the reproductive output associated with the transition j → i; that is, These matrices R can be obtained in several ways, discussed in detail in Caswell (2011).

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.

When only a single offspring is produced, mean offspring production equals the probability of reproducing. The matrices of second and third moments satisfy

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

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 stage-specific 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 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 . Our task is to derive the statistics of LRO from this information.

Let 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 corresponding to absorbing states are zero. Thus, we obtain the complete statistics of LRO if we can find the moment vectors , where contains the first τ entries of , corresponding to the living stages. This vector is given by where in terms of this matrix, . We also define , the τ × τ submatrix of R 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 are

and, in general,

where N =(I −U)−1 is the fundamental matrix of the Markov chain.

Proof

See Appendix A □

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 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 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 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, 2008, 2009, 2012; Caswell 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 sensitivity of the moment vector to the parameter vector is the derivative matrix That is, the (i,j) entry of this sensitivity matrix is the derivative of the ith entry of with respect to the jth entry of .

To write the sensitivity of the , let us define the following matrices. In Eqs. 27–29, K is the vec-permutation matrix of order (m,n); see Henderson and Searle (1981) and Magnus and Neudecker (1979).

Theorem 2

Let P, U, and M define the absorbing Markov chain in Eq. 3, with τ transient states, α absorbing states, and s = τ + α total states. Let R contain the ith moments of the reproductive rewards corresponding to each transition. Let be a vector of parameters. The vector contains the m th moments of remaining LRO for each of the τ transient stages. The sensitivity of to is

Proof

Proof is given in Appendix A. □

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 where 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).

Sensitivity of the statistics of LRO

The moments 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.

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 “Age-classified 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., , and G is a matrix of transition probabilities conditional on survival. Differentiating (43) gives Writing 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 must be compensated for by changes in other elements of column j. Caswell (2001, 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 denote column i of G; then where C = E −I.

Sensitivity to means and variances of fertility

Perturbations of fertility appear in Eqs. 27–29 as derivatives of the reproductive reward matrices R . 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 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 d v = 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 d f 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 Table 1. In the next section, we present age-classified and stage-classified examples of the analysis.

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.

2. 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)], or Poisson [ see Eqs. 6 and 7], or fixed [see Eqs. 8 and 9] reproduction.

4. Construct reward matrices from Eq. 11.

5. Compute desired moments of LRO from Eqs. 16–19.

6. Compute desired statistics of LRO from Eqs. 21–24.

7. Sensitivity analysis

(a) Specify parameter vector 𝜃 of interest

(b) Calculate derivatives of U, and R i to 𝜃. Take advantage of Eqs. 35 or 90–92 to compute derivatives of P to 𝜃.

(c) If the matrix is stage-classified,

i. Decompose U = GΣ.

ii. Use Eq. 46 to compute the derivative of U to mortality.

iii. Use Eq. 49 to compute the derivative of G to 𝜃, including compensation to preserve column sums of G

(d) Compute derivatives of the moment vectors ρ~_i for the moments of interest (i = 1, 2 suffice to analyze the variance, standard deviation,

CV, and I).

i. Compute V, W, and X using Eqs. 27–29.

ii. Compute derivatives of ρ~_i using Theorem 2.

(e) Compute sensitivity of desired statistics of LRO using Eqs. 38–41.

(f) If desired, compute elasticities of statistics of LRO using Eq. 42.

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 size-classified 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 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 , 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 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.

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

Population	Mean	V	V between (%)	V within (%)	SD	CV	I	Life exp.
NLD 1950	2.96	2.91	12.5	87.5	1.71	0.58	0.33	73.1
NLD 2009	1.78	1.61	1.4	98.6	1.27	0.72	0.51	83.1
SWE 1891	3.00	5.60	55.6	44.4	2.37	0.79	0.62	53.0
SWE 2010	1.97	1.79	1.4	98.6	1.34	0.68	0.46	84.0
JPN 1947	3.50	6.10	54.8	45.2	2.47	0.71	0.50	54.2
JPN 2009	1.35	1.26	0.9	99.1	1.12	0.83	0.69	86.9
Hadza	3.13	11.30	78.9	21.1	3.36	1.07	1.15	34.6
Ache	4.48	17.23	81.1	18.9	4.15	0.93	0.86	38.0
Hutterites	7.53	8.58	41.3	58.7	2.93	0.39	0.15	70.0

Not surprisingly, recent populations in developed countries have higher longevity and lower mean lifetime reproductive output. The hunter-gatherer and Hutterite populations show the highest mean LRO. In each of the developed countries (the Netherlands, Sweden, and Japan), the reductions in mean LRO and increases in life expectancy are accompanied by reductions in the variance in LRO. The Ache and Hadza show the highest variance in LRO, and the variance for the Hutterites is higher than any of the developed countries.

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 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 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 /d f 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).

The sensitivity of mean LRO, variance in LRO, and Crow’s 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.

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

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 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 combines both the mean and the variance. Increased mortality during the reproductive period increases 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 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 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 .

Fig. 2

Elasticity 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

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 population matrix A obtained from the Compadre Plant Matrix Database (Compadre Plant Matrix Database 2014). 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 . The variance in LRO exceeds that for any of the human populations by two orders of magnitude. The value of 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.

Table 3

The statistics of lifetime reproductive output for T. canadensis

Mean	V	V between (%)	V within (%)	SD	CV	I	Life exp.
1.42	1.41 × 103	99.9	0.1	37.54	26.39	696.23	12.16

As in Eq. 43, the transition matrix U is the product of the survival matrix and a growth matrix G. The derivations of the equations below are presented in Appendix A. The derivative of the transition matrix U with respect to mortality is The sensitivity of U to the growth matrix G, maintaining the column sums, is where dvec G/dvecT is given by Eq. 49.

The sensitivity of LRO to fertility is determined by the derivatives of the R 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 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.

Fig. 3

Elasticity of mean LRO and variance in LRO to changes in stage-specific mortality (left column) and stage-specific fertility (right column) for T. canadensis

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.

Fig. 4

Sensitivity of the mean LRO and variance in LRO of Tsuga canadensis to the growth matrix (left and right panels, respectively)

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 (2012). It calls for a comparison with models explicitly incorporating heterogeneity, either observed or unobserved, as suggested by Cam et al. (2016). For examples of this approach, see Caswell (2014); Hartemink et al. (2017); Jenouvrier et al. (2016).

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 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 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 to fertility is negative. The elasticity of mean and variance to fertility is positive, but the elasticity of Crow’s 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 hunter-gatherers 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, 2009, 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.