Constructing stagestructured matrix population models from life tables: comparison of methods
 Published
 Accepted
 Received
 Academic Editor
 Andrew Noymer
 Subject Areas
 Conservation Biology, Ecology, Statistics, Natural Resource Management
 Keywords
 Lefkovitch matrix, Leslie matrix, Life table analysis, Survivorship, Fecundity schedule, Life history, Matrix population models, Population matrix, EulerLotka equation
 Copyright
 © 2017 Fujiwara and DiazLopez
 Licence
 This is an open access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, reproduction and adaptation in any medium and for any purpose provided that it is properly attributed. For attribution, the original author(s), title, publication source (PeerJ) and either DOI or URL of the article must be cited.
 Cite this article
 2017) Constructing stagestructured matrix population models from life tables: comparison of methods. PeerJ 5:e3971 https://doi.org/10.7717/peerj.3971 (
Abstract
A matrix population model is a convenient tool for summarizing per capita survival and reproduction rates (collectively vital rates) of a population and can be used for calculating an asymptotic finite population growth rate (λ) and generation time. These two pieces of information can be used for determining the status of a threatened species. The use of stagestructured population models has increased in recent years, and the vital rates in such models are often estimated using a life table analysis. However, potential bias introduced when converting agestructured vital rates estimated from a life table into parameters for a stagestructured population model has not been assessed comprehensively. The objective of this study was to investigate the performance of methods for such conversions using simulated life histories of organisms. The underlying models incorporate various types of life history and true population growth rates of varying levels. The performance was measured by comparing differences in λ and the generation time calculated using the EulerLotka equation, agestructured population matrices, and several stagestructured population matrices that were obtained by applying different conversion methods. The results show that the discretization of age introduces only small bias in λ or generation time. Similarly, assuming a fixed age of maturation at the mean age of maturation does not introduce much bias. However, aggregating agespecific survival rates into a stagespecific survival rate and estimating a stagetransition rate can introduce substantial bias depending on the organism’s life history type and the true values of λ. In order to aggregate survival rates, the use of the weighted arithmetic mean was the most robust method for estimating λ. Here, the weights are given by survivorship curve after discounting with λ. To estimate a stagetransition rate, matching the proportion of individuals transitioning, with λ used for discounting the rate, was the best approach. However, stagestructured models performed poorly in estimating generation time, regardless of the methods used for constructing the models. Based on the results, we recommend using an agestructured matrix population model or the EulerLotka equation for calculating λ and generation time when life table data are available. Then, these agestructured vital rates can be converted into a stagestructured model for further analyses.
Introduction
A matrix population model is a convenient tool for summarizing per capita survival and reproduction rates (collectively vital rates) of a population, and is used for calculating an asymptotic finite population growth rate (commonly denoted by λ) and generation time. These two pieces of information can be used for determining the status of a threatened species (IUCN, 2012). Matrix population models are broadly categorized into agestructured (Leslie matrix; Leslie, 1945) and stagestructured (Lefkovitch matrix; Lefkovitch, 1965) models. Agestructured matrix models group individuals based on age whereas stagestructured matrix models group individuals based on other properties such as developmental stage and size. Stagestructured matrix models are often favored when a property of individuals besides age is a better indicator of survival and reproduction (e.g., Caswell, 2001; Cochran & Ellner, 1992). Although both types of models are common, the use of stagestructured population models has increased in recent years. Moreover, in many studies, vital rates in stagestructured population models may be estimated from agestructured vital rates through the grouping of ageclasses together to form a stage (e.g., Brault & Caswell, 1993; Caswell et al., 1998; Crouse, Crowder & Caswell, 1987; Crowder et al., 1994). Typically a life table analysis is one of the main methods for obtaining agespecific vital rates.
A life table is a list of agespecific population density and agespecific fecundity. From life table data, survivorship (the proportion of age 0 that is alive at age x) is estimated, and then agespecific survival rates (the proportion of individuals that survive from age x to x + 1) are estimated. Agespecific survival rates along with agespecific fecundity rates can be entered into an agestructured matrix population model almost directly (see the AgeStructured Matrix Population Models section), and λ and generation time can be calculated from the matrix. Matrix population models are also used for sensitivity and elasticity analyses (Caswell, 1978; De Kroon, Groenendael & Ehrlen, 2000). The inclusion of a large number of ageclasses for longlived organisms can make the interpretation of the sensitivity and elasticity analyses complicated because individuals in multiple age classes are often practically identical but separated in an agestructured model. Consequently, when longlived organisms are studied, it is common to convert agespecific vital rates into stagespecific vital rates, and to use stagestructured population matrices for calculating λ and generation time.
In order to convert agespecific into stagespecific vital rates, the former vital rates need to be aggregated for the various stages. Furthermore, a transition rate from one stage to another needs to be estimated. Several approaches, or conversion methods, to aggregate survival rates for calculating a transition rate exist. However, the performance of the conversion methods has not been investigated comprehensively. Intuitively, the performance should depend on how survivorship and reproductive schedule change with age (life history) and whether a population is growing or declining in its abundance. Therefore, it is critically important to investigate the performance of the conversion methods for different life history and population growth scenarios.
The objective of this study is to investigate the performance of the conversion methods to estimate vital rates for stagestructured matrix population models from life tables. The performance is measured by comparing λ’s and generation time calculated with the EulerLotka equation (see Kot, 2001), agestructured population matrices, and several stagestructured population matrices obtained using different conversion methods. The asymptotic finite population growth rate (λ) and the generation time calculated with the EulerLotka equation are considered the true values. Any discrepancies between the results obtained from the EulerLotka equation and those from an agestructured population matrix are considered to be due to bias introduced by discretizing age. In contrast, discrepancies between results from an age structured population matrix and those from stagestructured population matrices are considered to be due to bias introduced by conversion methods.
The comparisons were carried out with the life table data of organisms with different life history types. A wide range of life history types was incorporated into the analysis by artificially creating them using a competing risk model (Siler, 1979) and three types of fecundity functions. However, the analysis focuses on those with a prolonged duration in a stage of interest. When the duration is short, stage structure is often embedded within age structure, and it is not necessary to convert agestructured vital rates into stagestructured rates. The comparisons were also repeated with life table data from populations that are declining (λ < 1), maintaining the same population density (λ = 1), or increasing (λ > 1) in population density in order to determine whether a stagestructured population matrix can be used for determining a population growth rate correctly under different population growth conditions.
This paper is structured as follows. First, the life history types of organisms considered and the conversion methods to be compared are described. Then, methods for calculating λ’s and generation time with the EulerLotka equation and matrix population models are described followed by the description of the procedures specific to this study. Finally, results are presented and discussed.
Model  α_{1}  α_{2}  β_{2}  L_{∞}  κ  ${R}_{\lambda =0.9}^{1}$  ${R}_{\lambda =1.0}^{1}$  ${R}_{\lambda =1.1}^{1}$  x_{0}  x_{1}  x_{2}  x_{3} 

Early maturation type  
Dimension  τ^{−1}  τ^{−1}  τ^{−1}  l  τ^{−1}  μ^{−1}τl^{3}  μ^{−1}τ  μ^{−1}τ  τ  τ  τ  τ 
CHIF  .1535  0  –  10  .10  6,832  948  280  1.5  1.5  1.5  40.5 
CHCF  .1535  0  –  10  –  16,371  5,163  2,768  1.5  1.5  1.5  40.5 
CHDF  .1535  0  –  10  .10  6,262  3,134  1,975  1.5  1.5  1.5  40.5 
IHIF  0  .002  .205  10  .10  39,459  5,985  1,435  1.5  1.5  1.5  40.5 
IHCF  0  .002  .205  10  –  80,276  18,485  7,244  1.5  1.5  1.5  40.5 
IHDF  0  .002  .205  10  .10  22,872  8,100  4,176  1.5  1.5  1.5  40.5 
Model  α_{1}  α_{3}  β_{3}  m_{1}  m_{2}  a  ${R}_{\lambda =0.9}^{1}$  ${R}_{\lambda =1.0}^{1}$  ${R}_{\lambda =1.1}^{1}$  x_{0}  x_{1}  x_{2}  x_{3} 

Delayed maturation type  
Dimension  τ^{−1}  τ^{−1}  τ^{−1}  gτ^{−1}  τ^{−1}  τ^{−1}  μ^{−1}τl^{3}  μ^{−1}τ  μ^{−1}τ  τ  τ  τ  τ 
DHFM  .1000  .400  .300  –  –  .20  3.000  0.497  0.124  1.5  10.5  10.5  40.5 
CHFM  .2193  0  –  –  –  .20  3.005  0.498  0.124  1.5  10.5  10.5  40.5 
DHCM  .1000  .400  .300  .50  0  .20  2.958  0.502  0.131  1.5  8.5  10.5  40.5 
CHCM  .2193  0  –  .50  0  .20  3.082  0.538  0.143  1.5  8.5  10.5  40.5 
DHIM  .1000  .400  .300  .10  1.01  .20  2.942  0.491  0.124  1.5  8.5  10.5  40.5 
CHIM  .2193  0  –  .10  1.01  .20  3.030  0.510  0.130  1.5  8.5  10.5  40.5 
Notes:
 CH

Constant hazard (risk) model (Eq. (1))
 IH

Increasing hazard (risk) model (Eq. (2))
 DH

Decreasing hazard (risk) model (Eq. (6))
 IF

Increasing fecundity model (Eq. (3))
 CF

Constant fecundity model (Eq. (4))
 DF

Declining fecundity model (Eq. (5))
 FM

Fixed age of maturation at age x_{2}
 CM

Constant rate of maturation
 IM

Exponentially increasing rate of maturation
 τ

age
 l

length
 μ

number of births per adult
 g

number of individuals reaching maturity per juvenile
Models and Methods
Life history types
The types of life histories considered are broadly categorized into two groups: one with a short juvenile stage (early maturation) and the other with a prolonged juvenile stage (delayed maturation). The early maturation type was used for comparing methods for converting agespecific survival rates and fecundity rates into a stagespecific adult survival rate and a stagespecific fertility rate. The delayed maturation type was used for comparing methods for converting agespecific survival rates into a stagespecific juvenile survival rate and a stage transition rate. Throughout this study, the unit for time and age was set to a year for convenience, but it can be scaled differently if both are changed in the same way. Age is denoted by x, and organisms experience five major life events: birth at age 0, beginning of the juvenile stage at age x_{0}, age at which the first individual mature x_{1}, the mean age of maturation x_{2}, and the last age of reproduction x_{3}. Hereafter, the subscripts for age x denotes life events. The parameters for all models to be described are listed in Table 1.
Early maturation types
For the early maturation types, two types of survivorship and three types of fecundity schedules were incorporated. For the first type of survivorship, individuals experience a constant risk of mortality. This leads to an exponentially declining survivorship curve ($l\left(x\right)$) with age: (1)$l\left(x\right)={\mathrm{e}}^{{\alpha}_{1}x},$ where α_{1} is the constant risk (i.e., hazard rate). For the second type of survivorship, individuals experience an exponentially increasing risk of mortality with age: (2)$l\left(x\right)={\mathrm{e}}^{\frac{{\alpha}_{2}}{{\beta}_{2}}\left(1{e}^{{\beta}_{2}x}\right)},$ where α_{2} and β_{2} are the parameters for the increasing risk. Under this model, the hazard function is given by $h\left(x\right)={\alpha}_{2}{e}^{{\beta}_{2}x}$ (Siler, 1979). This type of risk may be due to aging. These two types of survivorship curves are shown in Fig. 1A.
Figure 1B shows the three types of fecundity schedules. The first type assumes the fecundity $b\left(x\right)$ is proportional to the cube of the body length of the individuals, which in turn increases according to von Bertalanffy growth equation: (3)$b\left(x\right)={R}_{1}{\left[{L}_{\infty}\left(1{e}^{\kappa x}\right)\right]}^{3},$ where L_{∞} and κ are the parameters in the growth model (Von Bertalanffy, 1968). This assumes reproduction is approximately proportional to the mass of individuals. The second type assumes a constant fecundity with age (4)$b\left(x\right)={R}_{2}.$ Finally, the third type assumes exponentially declining fecundity: (5)$b\left(x\right)={R}_{3}{e}^{\gamma x},$ where γ is the parameter determining how fast fecundity declines with age. Parameter R_{m} (where m = 1, 2, or 3) are constants with different units depending on the type of fecundity, and these constants were used for adjusting λ’s in this study (see ‘Analytical Procedure’). For all of the early maturation types in this study, individuals are assumed to mature at age 1.5 (i.e., x_{1} = x_{2} = 1.5).
Delayed maturation types
For the delayed maturation types, two kinds of survivorship curves are incorporated into a juvenile stage. One incorporates a constant risk of mortality (Eq. (1)), which results in exponentially declining survivorship (dashed curve in Fig. 1C). The other assumes an exponentially declining risk of mortality in addition to a constant risk: (6)$l\left(x\right)={\mathrm{e}}^{{\alpha}_{1}x}{e}^{\left[\frac{{\alpha}_{3}}{{\beta}_{3}}\left(1{e}^{{\beta}_{3}x}\right)\right]},$ where α_{3} and β_{3} are the parameters for the declining risk. In this model, the exponentially declining hazard is given by $h\left(x\right)={\alpha}_{3}{e}^{{\beta}_{3}x}$ (Siler, 1979). This causes initial fast descent of survivorship (solid curve in Fig. 1C). This type of survivorship can be observed when organisms grow out of a mortality risk such as predation as they age. Once individuals mature, they reach a constant risk of mortality (Eq. (1)) and constant fecundity (Eq. (4)) under the delayed maturation types.
In addition, three types of maturation schedule are incorporated. The first type assumes a fixed age of maturation at age 10.5 years (x_{1} = x_{2} = 10.5). The second type assumes a constant rate of maturation, and the third type assumes an exponentially increasing rate of maturation with age. In the latter two cases, maturation begins at age 8.5 years (x_{1} = 8.5), and the parameters are chosen so that the mean age of maturity is approximately 10.5 years (x_{2} = 10.5); however, the mean age varies slightly depending on the level of mortality rate during the maturation period. Under the latter two types of maturation schedules, the densities of individuals in immature (all stages prior to maturation) and adult stages are described conveniently with a system of ordinary differential equations (ODEs): (7)$\begin{array}{c}\frac{d{n}_{1}\left(x\right)}{dx}=m\left(x\right){n}_{1}\left(x\right)h\left(x\right){n}_{1}\left(x\right)\hfill \\ \frac{d{n}_{2}\left(x\right)}{dx}=m\left(x\right){n}_{1}\left(x\right)a{n}_{2}\left(x\right)\hfill \end{array},$ where ${n}_{1}\left(x\right)$ and ${n}_{2}\left(x\right)$ are the densities of immature individuals and adults at age x, respectively, $m\left(x\right)$ is a percapita instantaneous maturation rate at age x, and $h\left(x\right)$ is the agedependent hazard function of mortality. For the second type of maturation schedule, $m\left(x\right)$ is set at m_{1} = 0.5 for x ≥ x_{1} and 0 otherwise. For the third type, $m\left(x\right)={m}_{1}{e}^{{m}_{2}\left(x{x}_{1}\right)}$ for x ≥ x_{1}, and $m\left(x\right)=0$ otherwise. The proportion of individuals that are mature at a given age for the three types of maturation schedule is shown in Fig. 1D. When the system of ODEs is solved with the initial condition $\left[\begin{array}{cc}{n}_{1}\left(0\right)\hfill & {n}_{2}\left(0\right)\hfill \end{array}\right]=\left[\begin{array}{cc}1\hfill & 0\hfill \end{array}\right]$ using a numerical ODE solver, the survivorship is given by $l\left(x\right)={n}_{1}\left(x\right)+{n}_{2}\left(x\right)$. Note Eq. (7) does not include reproduction, and the initial condition is 1 in the immature stage; therefore, $l\left(x\right)\le 1$ for all x ≥ 0.
Conversion methods
In order to convert agestructured vital rates into stagestructured vital rates, three types of vital rates need to be calculated: a stagespecific survival rate for stage i (S_{i}), a transition rate from stage i to stage j conditional on their survival (P_{j,i}), and a fertility rate for stage i (F_{i}). The methods for converting these vital rates are described in the following sections.
Stagespecific survival rate
All the methods for calculating a stagespecific survival rate from agespecific survival rates investigated in this study assume that the beginning and ending ages of a stage are given by the mean ages of transition into and from the stage, respectively, and that these ages are known. The former assumption can potentially produce bias when calculating λ and generation time because stage transitions often do not occur at a fixed age (i.e., individuals gradually transition from one stage to another over a range of ages). Therefore, the performance of conversion methods was also investigated when the age of stagetransitions was not fixed (Eq. (7)). It is also possible that stagetransitions are completely ageindependent (e.g., as in some plant species). This situation is beyond the scope of the current analysis, but it will be briefly discussed in the ‘Discussion’ section.
An agespecific survival rate from age x to x + 1 is given by (8)${s}_{x}=\frac{l\left(x+1\right)}{l\left(x\right)}.$ Hereafter, lower case s is used for denoting an agespecific survival rate whereas the upper case S is used for denoting a stagespecific survival rate. Suppose the stage begins at age ${x}^{\left(i\right)}$ and ends at ${x}^{\left(j\right)}1$ (i.e., age ${x}^{\left(j\right)}$ is the next stage), then the survival rate of individuals in the stage is given by the mean survival rate over age classes within the stage. Note that an index in the superscript within parentheses denote a stage. Here, nine different ways of calculating the mean are compared. First, the mean is given by geometric mean, arithmetic mean, or harmonic mean (Table 2). These are three common methods for obtaining a stagespecific survival rate. However, the number of individuals in age classes are different, and it is more appropriate to put weight based on the proportion of individuals in a given age class, which is given by a survivorship curve. This leads to a weighted geometric mean, weighted arithmetic mean, and weighted harmonic mean (Table 2). Although weighting with survivorship is fine when a population is neither growing nor declining (i.e., λ = 1), when a population is growing (or declining), there are more (or less) individuals in younger age classes than predicted by a survivorship curve. Therefore, it is necessary to discount the weight using λ. This leads to weighted means where the weight is given by both survivorship and λ (Table 2). Hereafter, the latter weight is termed “a discounted weight.”
Conditional transition rate
Transition rate from stage i to stage j conditional on their survival (p_{j,i}) is calculated in three ways (Table 3). The first method (T1) matches the expected duration in the stage assuming exponentially declining time to transition into the following stage. The second method (T2) matches the expected proportion of individuals making the transition into the following stage assuming the survivorship curve gives the distribution of individuals among ageclasses. The third method (T3) is the same as the second method except that the distribution is discounted by a population growth rate. These methods are described in more detail in Caswell (2001).
Type of mean  Formula 

Geometric  ${S}_{i}={\left({\prod}_{x={x}_{i}}^{{x}_{j}1}{s}_{x}\right)}^{\frac{1}{\left({x}_{j}{x}_{i}\right)}}$ 
Arithmetic  ${S}_{i}=\frac{1}{{x}_{j}{x}_{i}}{\sum}_{x={x}_{i}}^{{x}_{j}1}{s}_{x}$ 
Harmonic  ${S}_{i}={\left(\frac{{\sum}_{x={x}_{i}}^{{x}_{j}1}{s}_{x}^{1}}{{x}_{j}{x}_{i}}\right)}^{1}$ 
Weighted geometric  ${S}_{i}={\left({\prod}_{x={x}_{i}}^{{x}_{j}1}{s}_{x}^{{\omega}_{1,x}}\right)}^{\frac{1}{{\sum}_{x={x}_{i}}^{{x}_{j}1}{\omega}_{1,x}}}$ where ${\omega}_{1,x}=l\left(x\right)$ 
Weighted arithmetic  ${S}_{i}=\frac{{\sum}_{x={x}_{i}}^{{x}_{j}1}{\omega}_{1,x}{s}_{x}}{{\sum}_{x={x}_{i}}^{{x}_{j}1}{\omega}_{1,x}}$ where ${w}_{1,x}=l\left(x\right)$ 
Weighted harmonic  ${S}_{i}={\left(\frac{{\sum}_{x={x}_{i}}^{{x}_{j}1}{\omega}_{1,x}{s}_{x}^{1}}{{\sum}_{x={x}_{i}}^{{x}_{j}1}{\omega}_{1,x}}\right)}^{1}$ where ${\omega}_{1,x}=l\left(x\right)$ 
Discountedweight geometric  ${S}_{i}={\left({\prod}_{x={x}_{i}}^{{x}_{j}1}{s}_{x}^{{\omega}_{2,x}}\right)}^{\frac{1}{{\sum}_{x={x}_{i}}^{{x}_{j}1}{\omega}_{2,x}}}$ where ${\omega}_{2,x}=\frac{1}{{\lambda}^{x{x}_{i}}}l\left(x\right)$ 
Discountedweight arithmetic  ${S}_{i}=\frac{{\sum}_{x={x}_{i}}^{{x}_{j}1}{\omega}_{2,x}{s}_{x}}{{\sum}_{x={x}_{i}}^{{x}_{j}1}{\omega}_{2,x}}$ where ${\omega}_{2,x}=\frac{1}{{\lambda}^{x{x}_{i}}}l\left(x\right)$ 
Discountedweight harmonic  ${S}_{i}={\left(\frac{{\sum}_{x={x}_{i}}^{{x}_{j}1}{\omega}_{2,x}{s}_{x}^{1}}{{\sum}_{x={x}_{i}}^{{x}_{j}1}{\omega}_{2,x}}\right)}^{1}$ where ${\omega}_{2,x}=\frac{1}{{\lambda}^{x{x}_{i}}}l\left(x\right)$ 
MODEL  Type of transition rate  Formula 

T1  Matching duration  ${P}_{j,i}=1\frac{1}{{x}_{j}{x}_{i}}$ 
T2  Matching proportion transitioning  ${P}_{j,i}=\frac{l\left({x}_{j}1\right)}{{\sum}_{x={x}_{i}}^{{x}_{j}1}l\left(x\right)}$ 
T3  Matching proportion transitioning with discount  ${P}_{j,i}=\frac{{\lambda}^{\left({x}_{j}{x}_{i}1\right)}l\left({x}_{j}1\right)}{{\sum}_{x={x}_{i}}^{{x}_{j}1}{\lambda}^{\left(x{x}_{i}\right)}l\left(x\right)}$ 
Fertility rate
Fertility rates in matrix population models are different from fecundity in a life table. A fertility rate gives the per capita rate of contribution from one stage to the next over one year so that adults will have to survive to reproduce, and/or offspring must survive to appear in the next stage. On the other hand, fecundity is the number of offspring produced. The latter does not include any survival of adult or offspring. Consequently, there are two steps in converting life table data into a stagestructured fertility rate: calculating agespecific fertility rates and converting them into a stagespecific fertility rate. There are many approaches for the first step depending on the reproductive schedules of organisms within a year (see Caswell, 2001) or among years (e.g., Crowder et al., 1994). Evaluating them is beyond the scope of this study. Here, reproduction was assumed to occur at any time of year (i.e., a birth flow model), and an agespecific fertility rate (f_{x}) was estimated as follows: (9)${f}_{x}=l\left(0.5\right){\int}_{z=x0.5}^{x+0.5}b\left(z\right)\frac{l\left(z\right)}{l\left(x0.5\right)}dz,$ where $b\left(z\right)$ is the instantaneous percapita fecundity rate at age x. In the equation, $l\left(z\right)\u2215l\left(x0.5\right)$ gives the survival rate of adults from age x − 0.5 to age z (where z > x − 0.5), and the integral calculates the total number of offspring produced between ages x − 0.5 and x + 0.5 per adult that was alive at age x − 0.5. Then, all births are attributed to the half way point (i.e., age x) so that offspring will have to survive half a year on average to appear in the first stage; this survival rate is given by $l\left(0.5\right)$, which is the survival rate from age 0 to 0.5.
Once agespecific fertility rates are obtained, a stagespecific fertility rate F_{i} is calculated in three different ways. First, an arithmetic mean is taken with equal weights on all age classes, (10)${F}_{i}=\frac{1}{{x}_{3}{x}_{2}}\sum _{x={x}_{2}+0.5}^{{x}_{3}0.5}{f}_{x}.$ The second approach is to take a weighted arithmetic mean, (11)${F}_{i}=\frac{1}{\sum _{x={x}_{2}+0.5}^{{x}_{3}0.5}{\omega}_{1,x}}\sum _{x={x}_{2}+0.5}^{{x}_{3}0.5}{\omega}_{1,x}{f}_{x},$ where ${\omega}_{1,x}=l\left(x\right)$. The third method is to use a discounted weight by replacing ω_{1,x} in Eq. (11) with ${\omega}_{2,x}=\frac{1}{{\lambda}^{x{x}_{0}}}l\left(x\right)$. In all calculations, the same weight (or no weight) is used for calculating a stagespecific survival rate and a stagespecific fertility rate.
Asymptotic population growth rate and generation time
An asymptotic population growth rate (a finite per capita population growth rate λ) and generation time are calculated using the EulerLotka equation, agestructured (Leslie) matrices, and stagestructured (Lefkovitch) matrices. These models and the calculations of λ and generation time are described in this section.
EulerLotka equation
The EulerLotka equation (Kot, 2001) is given by (12)$1={\int}_{0}^{{x}_{3}}{\lambda}^{x}b\left(x\right)l\left(x\right)dx.$ The fecundity $b\left(x\right)$ and survivorship $l\left(x\right)$ are defined for a specific life history type as described in the previous section (Fig. 1). Provided R_{m} is fixed, the only unknown is λ, which can be found numerically by searching the value that satisfies Eq. (12) using a rootfinding algorithm.
Generation time in this study is defined as the mean age of parents where the mean is calculated over all offspring born at a given time. More specifically, generation time G_{1} is given by (13)${G}_{1}={\int}_{0}^{{x}_{3}}x{\lambda}^{x}b\left(x\right)l\left(x\right)dx.$ (Keyfitz & Caswell, 2005). In the case where individuals mature over a range of ages, some individuals do not reproduce in certain age ranges. Therefore, generation time G_{2} is instead given by (14)${G}_{2}={\int}_{0}^{{x}_{3}}x{\lambda}^{x}b\left(x\right)q\left(x\right)l\left(x\right)dx,$ where $q\left(x\right)$ is the proportion of individuals that are mature at age x.
Agestructured matrix population models (Leslie matrix)
An agestructured population matrix is given as (15)$\mathbf{A}=\left[\begin{array}{ccccc}\hfill 0\hfill & \hfill {f}_{2}\hfill & \hfill \dots \hfill & \hfill {f}_{{x}_{3}1.5}\hfill & \hfill {f}_{{x}_{3}0.5}\hfill \\ \hfill {s}_{0.5}\hfill & \hfill 0\hfill & \hfill \dots \hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill {s}_{1.5}\hfill & \hfill \dots \hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill \dots \hfill & \hfill \dots \hfill & \hfill \dots \hfill & \hfill \dots \hfill & \hfill \dots \hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill \dots \hfill & \hfill {s}_{{x}_{3}2}\hfill & \hfill 0\hfill \end{array}\right].$ In this study, the first age class is assumed to begin at age 0.5, and subsequent age class is incremented by one year. It should be noted that survival rate s_{x} gives the proportion of individuals that survive from time x to x + 1 and fertility f_{x} gives the mean number of offspring of age 0.5 produced by a parent alive at age x − 0.5 by reproducing between age x − 0.5 and x + 0.5 (see Eq. (9)). It is assumed that the reproduction begins at the mean age of maturation x_{2}.
In the case where individuals mature over a range of age, the population matrix is given instead as (16)$\mathbf{A}=\left[\begin{array}{ccccc}\hfill 0\hfill & \hfill {q}_{2}{f}_{2}\hfill & \hfill \dots \hfill & \hfill {q}_{{x}_{3}1.5}{f}_{{x}_{3}1.5}\hfill & \hfill {q}_{{x}_{3}0.5}{f}_{{x}_{3}0.5}\hfill \\ \hfill {s}_{0.5}\hfill & \hfill 0\hfill & \hfill \dots \hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill {s}_{1.5}\hfill & \hfill \dots \hfill & \hfill 0\hfill & \hfill 0\hfill \\ \hfill \dots \hfill & \hfill \dots \hfill & \hfill \dots \hfill & \hfill \dots \hfill & \hfill \dots \hfill \\ \hfill 0\hfill & \hfill 0\hfill & \hfill \dots \hfill & \hfill {s}_{{x}_{3}2}\hfill & \hfill 0\hfill \end{array}\right].$
Once the population matrix is constructed, λ is obtained by taking its dominant eigenvalue. Furthermore, generation time is given by (17)$G}_{3}=\frac{\lambda \mathbf{vw}}{\mathbf{vFw}$ where v and w are the left and right eigenvectors (row and column vectors), respectively, of the population matrix, and F is a modified population matrix only consisting of fertility terms (Bienvenu & Legendre, 2015).
Stagestructured matrix population models (Lefkovitch matrix)
Stagestructured matrix population models are given in two different forms in this study. For early maturation types, a twostage population matrix is used (18)$\mathbf{A}=\left[\begin{array}{cc}0\hfill & {F}_{2}\hfill \\ {S}_{1}\hfill & {S}_{2}\hfill \end{array}\right]$ where the duration of stage 1 is from age x_{0} = 0.5 to age x_{1} = x_{2} = 1.5 and the duration of stage 2 is from age 1.5 to age x_{3} = 40.5.
For organisms with delayed maturation, threestage matrix population models are used (19)$\mathbf{A}=\left[\begin{array}{ccc}\hfill 0\hfill & \hfill 0\hfill & \hfill {F}_{3}\hfill \\ \hfill {S}_{1}\hfill & \hfill {S}_{2}\left(1{P}_{3,2}\right)\hfill & \hfill 0\hfill \\ \hfill 0\hfill & \hfill {S}_{2}{P}_{3,2}\hfill & \hfill {S}_{3}\hfill \end{array}\right].$ The duration in stages 1, 2, and 3 are assumed to be from age x_{0} = 0.5 to 1.5, from age 1.5 to x_{1} = x_{2} = 10.5, and from age x_{2} = 10.5 to x_{3} = 40.5. Using the stagestructured matrix, λ and generation time can be calculated in the same way as with the agestructured matrices.
Analytical procedure
Thirtysix different scenarios of life history strategies and a population growth rate were investigated in this study. For the early maturation type, there were six different life history strategies (two types of survivorships and three types of fecundity schedules; Figs. 1A and 1B). With each of the six life history strategies, the true finite asymptotic population growth rate (λ) was adjusted to 0.900, 1.000, and 1.100 by adjusting coefficient R_{m} in the fecundity function (Eqs. (3)–(5)). This adjustment was done with the EulerLotka equation by searching the value R_{m} that gives the corresponding value of λ. Consequently, there were 18 early maturation type scenarios to investigate. For each scenario, one agestructured population model and nine twostage population models were constructed; the nine models differ in the ways that adult survival rate and a fertility rate were calculated (see Table 2). Then, with each model, λ and generation time were calculated.
Similarly, for the delayed maturation type, there were six different life history strategies (two types of survivorships and three types of maturation schedules; Figs. 1C and 1D). For each type, finite asymptotic population growth rate (λ) was also adjusted to 0.900, 1.000, and 1.100. Then, two agestructured population models (Eqs. (15) and (16)) and nine threestage population models were constructed; these nine models differ in the ways that a juvenile survival rate was aggregated and a conditional transition rate was obtained (see Table 2). Then, with each model, λ and generation time were calculated.
In order to obtain the true generation time, Eq. (13) was used when organisms had a fixed age of maturation, while Eq. (14) was used when organisms matured over a range of age. All of the calculations were done with MathWorks (2012). For solving a system of ODE’s, a builtin ODE solver “ode45” was used in most cases. When a maturation rate was exponentially increasing with age, there was a numerical problem with the function (equations became stiff). Therefore, “ode15s” was used, instead. For the methods that require λ in aggregating a survival rate, a transition rate, and/or a fertility rate, the true λ was used.
In actual analyses conducted with real data, estimating λ is one of the main purposes of constructing a population matrix. In other words, it is not known a priori. To overcome this issue, Caswell (2001) describes an iterative method. In this method, an initial λ is arbitrarily selected and used for obtaining factors needed to calculate the conditional transition rate. Then, λ is estimated from the obtained population matrix. Subsequently, the new λ is used to obtain a new population matrix and λ is calculated again. This process is repeated until λ converges. For example, this method was used for developing a stagestructured population matrix of loggerhead sea turtles (Crowder et al., 1994). In the current study, this procedure is applied to delayed maturation life history strategy. Discountedweight arithmetic mean was used for aggregating juvenile and adult survival rates and a transition rate was obtained by matching proportion transitioning also using λ as a discounting factor.
Results
Finite asymptotic population growth rate
Figures 2–7 show estimated finite per capita population growth rates (λ). In each figure, each panel shows the estimated λ of the same life history type calculated using different methods. These methods are listed in Table 4. The first bar (black) is the true value, and the rest of the bars are obtained with age or stagestructured population matrices.
Models  Description  Generation time calculation 

a  EulerLotka equation  G_{1} or G_{2} 
b  Age structured model with proportion mature  G_{3} 
c  Age structured model with fixed age of maturity  G_{3} 
d  Twostage model with a geometric mean survival rate  G_{3} 
e  Twostage model with an arithmetic mean survival rate  G_{3} 
f  Twostage model with a harmonic mean survival rate  G_{3} 
g  Twostage model with a weighted geometric mean survival rate  G_{3} 
h  Twostage model with a weighted geometric mean survival rate  G_{3} 
i  Twostage model with a weighted harmonic mean survival rate  G_{3} 
j  Twostage model with a weighted geometric mean survival rate with λ used as discounting factor  G_{3} 
k  Twostage model with a weighted geometric mean survival rate with λ used as discounting factor  G_{3} 
l  Twostage model with a weighted harmonic mean survival rate with λ used as discounting factor  G_{3} 
m  Threestage model with an arithmetic mean survival rate and matching duration (T1)  G_{3} 
n  Threestage model with an arithmetic mean survival rate and matching proportion transitioning (T2)  G_{3} 
o  Threestage model with an arithmetic mean survival rate and matching proportion transitioning with discount (T3)  G_{3} 
p  Threestage model with a weighted arithmetic mean survival rate and matching duration (T1)  G_{3} 
q  Threestage model with a weighted arithmetic mean survival rate and matching proportion transitioning (T2)  G_{3} 
r  Threestage model with a weighted arithmetic mean survival rate and matching proportion transitioning with discount (T3)  G_{3} 
s  Threestage model with a discountedweight arithmetic mean survival rate and matching duration (T1)  G_{3} 
t  Threestage model with a discountweight arithmetic mean survival rate and matching proportion transitioning (T2)  G_{3} 
u  Threestage model with a discountweight arithmetic mean survival rate and matching proportion transitioning with discount (T3)  G_{3} 
True λ  Life history types (see Table 1)  

D.HF.M.  C.H.F.M.  D.H.C.M.  C.H.C.M.  D.H.I.M.  C.H.I.M.  
λ = 0.90  0.900  0.903  0.899  0.902  0.900  0.902 
λ = 1.00  0.994  1.000  0.992  0.995  0.994  0.999 
λ = 1.10  1.090  1.100  1.085  1.089  1.090  1.096 
For the early maturation types (Figs. 2–4), one agestructured model and nine stagestructured models that differ in the methods for calculating an adult survival rate were used. The three figures are different in the true λ’s. For all life history types and all true λ’s (Figs. 2–4), the EulerLotka equation (Model a) and agestructured model (model c) show the similar values (black and gray bars in Figs. 2–4). This suggests there is no bias introduced by discretizing age.
When adults experience a constant risk of mortality (Eq. (1)) and constant fecundity (Eq. (4)), all conversion methods performed well (within 2% of the true value; see of Figs. 2B–4B). In all other cases, some of the methods performed poorly (see Figs. 2A, 2C–2F, 3A, 3C–3F and 4A, 4C–4F).
Comparing the use of the geometric mean, arithmetic mean, and harmonic mean of adult survival rate (Models d–f in Figs. 2A–2C, 3A–3C and 4A–4C), all gave the similar results for adults experiencing a constant risk. However, for adults experiencing an exponentially increasing risk (Figs. 2D–2F, 3D–3F and 4D 4F, the three means (Models d–f) gave different results. In particular, the harmonic mean (Model f) grossly underestimated λ, as did the geometric mean (Model d) although to a lesser extent. Among the three, the arithmetic mean (Model e) performed better, but it still underestimated λ.
Weighted means (Models g–i) improved the estimation of λ, compared with unweighted means (Models d–e). For example, when the weighted arithmetic mean was used, declining populations (λ = 0.9) were always identified as declining (λ < 1), and increasing populations (λ = 1.1) were always identified as increasing (λ > 1). The use of discounted weight (Models j–k) further improved the estimation of λ, especially, when λ was not at 1.0 (Figs. 3 and 4).
For the delayed maturation types, two agestructured models and nine stagestructured models were used (Figs. 5–7). The two agestructured models differed in whether the age of maturation was assumed to be fixed at the mean age of maturity (Model b; Eq. (15)) or the proportion of mature individuals at a given age is incorporated (Model c; Eq. (16)). The nine stagestructured models differed in the method for calculating the conditional transition rate (Table 3) and in how juvenile survival was aggregated (Table 2). For the latter aggregation, arithmetic mean, weighted arithmetic mean, and discountedweight arithmetic mean were used. First, agestructured models gave λ consistent with the true value (cf. Model a and Models b–c of all panels of Figs. 5–7). This suggests discretization of age did not introduce bias in λ. Furthermore, the incorporation of the proportion of mature individuals did not improve λ (cf. Model b and c of all panels of Figs. 5–7), suggesting the assumption of the fixed age of maturation was appropriate.
For stagestructured models, bias introduced in the estimation of λ was not as large as that with the early maturation types. Among the three conditional transition rate estimation methods (Table 3), the third method (T3) to match the proportion of individuals making the transition after discounting with λ (Models o, r, and u) performed best, and the first method to match the average duration (Models m, p, and s) performed the worst. The estimation was less sensitive to how the juvenile survival rate was aggregated.
Generation time
Generation time is shown in Figs. 8–13. In all cases, agestructured models (Model b in Figs. 8–10 and Models b and c in Figs. 11–13) gave a generation time that was consistent with the true values (Model a). This suggests discretization of age does not introduce bias. However, stagestructured models gave substantially different values for most scenarios. They gave similar results to the true values only for populations experiencing constant hazard and constant fecundity for early maturation types as long as the population was steady (λ = 1, Fig. 8B) or growing (λ > 1, Fig. 10B). These results suggest generation time obtained from stagestructured models are generally not accurate.
Iterative methods
All of the iterative methods in this study converged very quickly. The estimated λ value obtained using the iterative method is shown in Table 5. Overall, the iterative method performed well in estimating λ. It was accurate to the first digit of λ at least. However, generation time was estimated almost equally poorly with that obtained assuming λ was known a priori (i.e., Model u of all panels of Figs. 11–13).
Discussion
Population growth rate (λ) and generation time are two of the most important pieces of information in determining the status of threatened species (see IUCN, 2012). The former tells us how quickly a population is expected to be declining or growing over time, and the latter tells us an appropriate timescale for a population. Accurate determination of these parameters is essential. For example, the overestimation of generation time or underestimation of λ will result in species being erroneously placed in categories of higher threat than they should be. In this study, we have evaluated the performance of stagestructured population models in calculating λ and generation time, when life table data are used to parameterize the models.
Overall, discretization of age does not introduce much bias in population growth rate (λ) or generation time (comparing the first and second bars of all panels of Figs. 2–13). Similarly, assuming a fixed age of maturation at the mean age of maturity, does not introduce much bias (comparing the second and third bars of all panels of Figs. 5–7 and 11–13). Although only three types of maturation schedule have been investigated in this study, the latter conclusion is expected to remain true as long as the agerange of maturation is not substantially wide. On the other hand, the aggregations of agespecific vital rates into a stagestructured vital rate introduces bias.
In order to aggregate survival rates for stagestructured population matrices, the use of the discountedweight arithmetic mean is the most robust method to estimate λ (Figs. 2–4). However, the method requires the use of λ to discount the weight, which defeats the purpose of constructing a stagestructured population model to obtain λ. Therefore, the use of the weighted arithmetic mean is a better option although it can still introduce some bias when the population is either growing or declining. For the calculation of the conditional transition rate, the best approach is to match the proportion transitioning, with λ being used as a discounting factor; however, the use of the method without the discounting factor also performs well (Figs. 5–7).
The study also suggests that aggregating adult stage needs to be done more carefully than aggregating juvenile stage (compare Figs. 2–4 and 5–7). Adults experience survival and reproduction whereas juveniles experience survival and maturation. Because they experience different population processes, it is not surprising to see the difference. However, the results may also be because the adult stage is much longer than the juvenile stage in this study. For example, some organisms exhibit a shortlived adult stage with delayed maturation. For such organisms, the aggregation for the juvenile stage would need to be done more carefully than for the adult stage.
One of the problems is the heterogeneity among individuals within a stage. In our model, the heterogeneity exists because either the risk of mortality or the fecundity rate changes with age. However, a matrix population model assumes that all individuals are the same within a stage. When both mortality risk and fecundity rate are constant over age, λ and generation time are estimated accurately with stagestructured population models (see second columns of Figs. 2 and 7). However, variation in mortality with age is commonly experienced. For example, young individuals may grow out of a predation risk, and older individuals may experience senescence. The former is common among any animal populations, and the latter is common among longlived mammals. Fecundity also changes with age. For example, it is often a function of the size of adults, which increases with age, a common change seen in fish. One way to overcome the issue of a nonconstant mortality risk or nonconstant fecundity may be to aggregate smaller number of age classes into one stage by including a large number of stages in a model. This could be done after fitting a competing risk model (Siler, 1979) and a flexible fecundity function to life table data and examining the age classes with a similar mortality risk and fecundity.
Our results suggest estimating generation time is more problematic than estimating λ. For example, when fecundity was declining with age under the early maturation model (Fig. 7C), generation time was always overestimated with stagestructured models. This is because although in reality younger individuals contribute more to reproduction than older individuals do, stagestructured models assume homogeneity among adults. In other words, older adults contribute more in the models than the true population. Consequently, generation time is overestimated. When fecundity was declining with age (Fig. 7A), the opposite effect was observed. When survival rate was changing with age, it can lead to over or underestimation of generation time, and it is difficult to predict. For these reasons, the construction of an age structured model or the use of the EulerLotka equation is recommended for calculating λ and generation time.
Lebreton (2005) discusses potential problems with using stagestructured population models to calculate generation time. Lebreton recommends the use of stagestructured models in which stages are embedded within age classes. This effectively increases the number of stages in a model. Such a population matrix can be constructed if one collects information on the stage of individuals in addition to age in order to construct a life table in which individuals are categorized into different age and stage classes. This approach is also useful when life history strategies are complex or when they do not exhibit agedependent changes in vital rates. For such organisms, a simple life table alone may not be informative because each age class can include multiple stages, but a stagestructured model embedded in age structure should provide accurate estimations of both λ and generation time. SalgueroGomez & Plotkin (2010) also discusses the relationship between the number of stages and various statistics obtained from stagestructured population matrices.
An additional problem is the change in stage/age distribution due to growing or declining population abundance. A population growth rate needs to be discounted when a survivorship curve is estimated from a life table. Discounting is required for the same reason that it is included in the EulerLotka equation. There are two types of life table, dynamic or static life table. A dynamic life table is also called a cohort life table because it is constructed using data obtained by following the same cohort over time/age. On the other hand, a static life table is obtained by examining the age distribution of samples collected at one sampling occasion. Throughout this paper, a dynamic life table has been assumed. When a static life table is collected instead, a survivorship curve needs to be discounted by a population growth rate to obtain the true survivorship curve. The discounting is necessary because individuals in different age classes in the data were born at different time and the population is growing or declining so that increasing or decreasing numbers of individuals are born from one year to the next. Interestingly, the age distribution obtained directly from a static life table gives the discounted weight used in this study.
The iterative method allows the use of λ in constructing a stagestructured population matrix without knowing λ a priori. This method performed well in identifying whether λ is greater or smaller than 1. In this sense, the iterative method is a viable option for constructing a stagestructured population model. However, when a life table is available, an agestructured matrix population model can be used for calculating λ more accurately. Furthermore, an agestructured matrix allows accurate estimation of generation time, which can be substantially biased in stagestructured population models.
In this study, bias was introduced when constructing matrix population models using information from life table data. As an alternative to life table data, Cormack–Jolly–Seber (CJS) type capturerecapture data (Lebreton et al., 1992) can be obtained. With additional information to categorize captured individuals into different stage classes, a multitype/multistage capture recapture method can be used for estimating parameters in stagestructured matrix population models directly (Fujiwara & Caswell, 2002; Nichols et al., 1992). Provided there is no heterogeneity in a capture rate within a given stage, the method accounts for an underlying agedistribution although an actual age distribution may not be observable. When age is not known, the use of stagestructured population models is the only option. However, generation time obtained from such models may not be accurate. Future investigations of the performance of generation time calculations with stagestructured population models when vital rates are estimated from individual capturerecapture data are needed.
Conclusions
When life table data are collected, we recommend fitting a competing risk model and a flexible fecundity function to the data and estimating population growth rate λ and generation time using an agestructured population matrix or the EulerLotka equation. Calculating generation time using stagestructured population models should be avoided. If a researcher is interested in constructing stagestructured population models (e.g., for the purpose of sensitivity and elasticity analyses), the conversion from agestructured vital rates to stagestructured vital rates should be done by aggregating age classes with similar mortality risk and fecundity into the same stage. When aggregating survival rates for constructing stagestructured population models, discountedweight arithmetic mean should be used. When calculating the conditional transition rate, one should use the method for matching the proportion making transitions with λ as a discounting factor.
Supplemental Information
MATLAB code and parameters
Complete MATLAB code used for the analysis and Excel file containing parameters. All of the files should be extracted and saved under the same folder. MATLAB code will read the Excel file and run the analysis. Please see “Supplemental Material.docx” for more information.