Variance in longevity among individuals may arise as an effect of heterogeneity (differences in mortality rates experienced at the same age or stage) or as an effect of individual stochasticity (the outcome of random demographic events during the life cycle). Decomposing the variance into components due to heterogeneity and stochasticity is crucial for evolutionary analyses.In this study, we analyze longevity from ten studies of invertebrates in the laboratory, and use the results to parti- tion the variance in longevity into its components. To do so, we fit finite mixtures of Weibull survival functions to each data set by maximum likelihood, using the EM algorithm. We used the Bayesian Information Criterion to select the most well supported model. The results of the mixture analysis were used to construct an age × stage-classified matrix model, with heterogeneity groups as stages, from which we calculated the variance in longevity and its components. Almost all data sets revealed evidence of some degree of heterogeneity. The median contribution of unobserved heterogeneity to the total vari- ance was 35%, with the remaining 65% due to stochasticity. The differences among groups in mean longevity were typically on the order of 30% of the overall life expectancy. There was considerable variation among data sets in both the magnitude of heterogeneity and the proportion of variance due to heterogeneity, but no clear patterns were apparent in relation to sex, taxon, or environmental conditions. Keywords Age-stage classified model · Heterogeneity · Individual stochasticity · Mixture models · Variance in longevity · Weibull distribution Introduction longevity may arise as a result of two different underlying causes: stochastic processes and heterogeneity between indi- Individual variance in fitness components is central to evo- viduals. That is, even in a population without heterogeneity, lutionary demography and ecology, since variation between in which all individuals experience identical age-specific individuals in their traits and the resulting consequences mortality rates, death would be a probabilistic event, lead- for fitness are the basis for natural selection. Longevity, or ing to variance in longevity among individuals. This source age at death, is such a fitness component that varies among of variance is called individual stochasticity (Caswell 2009). individuals within a cohort or population. This variance in On top of that, genuine heterogeneity in age-specific mortal- ity risk among individuals can be a cause of variance. Such differences between individuals with respect to their mortal- ity risk, especially unobserved differences, are often referred Electronic supplementary material The online version of this article (https ://doi.org/10.1007/s1014 4-018-0616-7) contains to as heterogeneity in individual frailty (Vaupel et al. 1979), supplementary material, which is available to authorized users. where frailty is defined as proneness to mortality. The impact of heterogeneity on demographic outcomes, * Hal Caswell eco-evolutionary processes, and population dynamics has email@example.com been the topic of several studies (e.g., Kendall and Fox Institute for Biodiversity and Ecosystem Dynamics, 2002; Robert et al. 2003; Kendall et al. 2011; Vindenes and University of Amsterdam, P.O. Box 94248, Langangen 2015; Cam et al. 2016). However, the extent to 1090 GE Amsterdam, The Netherlands Vol.:(0123456789) 1 3 90 Population Ecology (2018) 60:89–99 which variance in fitness components can be accounted for by We begin by describing the statistical estimation procedures individual stochasticity is still open (Cam et al. 2016). This and then develop the age × stage-classified matrix population question is fundamental to evolutionary demography because model. Subsequent sections calculate and decompose the vari- the two sources of variance have very different implications. ance in longevity. We show the results in a series of tables and Although it can arise from many other causes, variance due figures; detailed results for each data set are found in Elec- to heterogeneity may have a genetic basis, and hence play tronic Supplementary Materials [ESM-1 and 2]. a role in selection. Because variance due to individual sto- chasticity arises from individuals experiencing identical vital Materials and methods rates, by definition it cannot have a genetic basis. It may even slow down selection by obscuring genetic variance that does Finite mixture models for survival exist (Steiner and Tuljapurkar 2012). Automatically attribut- ing observed variance in fitness components to heterogeneity We described heterogeneity by a finite mixture model, in overestimates the potential for selection. which a discrete number g of groups are defined, each group The key to quantifying the relative contributions of het- having its own mortality parameters, and with group i com- erogeneity and stochasticity is to construct a demographic prising a proportion of the cohort at the initial age. The model in which both factors operate, and partition the result- mortality parameters of each group and the mixing distribu- ing variance into components due to each process. Such a tion are estimated by maximum likelihood from data on the method has been developed by using age×stage-classified observed distribution of age at death. matrix population models (Caswell 2014; Hartemink et al. Such finite mixture models are widely used in statis- 2017). The primary state variable is age, and some aspect tics (e.g., McLachlan and Peel 2004; Frühwirth-Schnatter of heterogeneity (e.g., frailty) is included as a stage. When 2006). They have long been applied in survival analysis as variance in, for example, longevity is calculated from such an alternative to continuous frailty models (Farewell 1982; a matrix, it can be decomposed into a variance component Heckman and Singer 1982; McLachlan and McGiffin 1994; between frailty classes, which is due to heterogeneity and Erişoğlu et al. 2012) such as the Gamma–Gompertz model a variance component within frailty classes, which is due (Vaupel and Carey 1993). Bijwaard (2014) and Putter and to stochasticity. A first, crude attempt to explore this ques- van Houwelingen (2015) have explored finite mixtures in tion in laboratory animal cohort data was made in Caswell multistate models. (2014), but that analysis relied on previously published We used the Weibull distribution to model survival. This estimates of uncertain methodology, based on a restrictive distribution is more flexible than the Gompertz model, because mortality model, and with only an approximate variance it permits increasing, decreasing, or constant hazard rates, thus decomposition. Here we present a more rigorous analysis. incorporating the type I, II, and III survivorship curves familiar We quantify the relative contribution of heterogeneity in ecology (Pinder et al. 1978). The Weibull distribution has and stochasticity to the variance in longevity in a range of appealing biological interpretations as the time to failure of a invertebrate laboratory animal studies, comprising 25 data system that relies on the continued operation of a large num- sets on 9 species of nematodes and insects, totaling about 3.2 ber of processes and fails when any one of them does so (e.g., million individuals. Heterogeneity in mortality was captured Horvath 1968), or as the result of accumulation of damage by fitting finite mixtures of Weibull functions to data on beyond a certain threshold (Rinne 2008). The Weibull hazard individual ages at death, using the expectation–maximiza- function is: tion (EM) algorithm. We used model selection criteria to k−1 k x choose the mixture model most well supported by the data, (x, k)= (1) and constructed the corresponding matrix model to estimate the components of variance. where x is age, (x) is the mortality hazard, is a scale In this paper we address several questions: (1) is there evi- parameter, and k is a shape parameter. In medical statis- dence for unobserved heterogeneity in mortality of inverte- tics, alternative parameterizations are sometimes used (e.g., brates under controlled conditions, (2) if so, how much, and Mills 2011), in which the shape parameter is the same as (−k) how are individuals distributed among heterogeneity groups, above, but the scale parameter is in our parametriza- and (3) how much of the variance in longevity is due to het- tion. In Mtlaba , the model is specified using a for the scale erogeneity and how much to individual stochasticity. Because parameter and b for the shape parameter. If k ≤ 1 , the hazard we are using data on a variety of species, on both sexes, and increases with time. If k > 1 , the hazard decreases over time. sometimes under different conditions, we will look for patterns If k = 1 , the hazard is constant and the model reduces to an related to these variables. But because our results are based on exponential model. The probability density function of age an arbitrary and non-random selection of data, constrained by at death for the Weibull distribution is: sample size, rigorous comparative analyses are impossible. 1 3 Population Ecology (2018) 60:89–99 91 X is X . The matrix K is the vec-permutation matrix (Hen- k−1 k k x x f (x, k)= exp − . (2) derson and Searle 1981). The estimated number of groups (which we will denote by g), the proportion of individuals in each group (described by Maximum likelihood estimation the mixing distribution vector ) and the Weibull parameters and k for each group serve as input for our age×stage-clas- The mixture models were fit to the data using maximum sified matrix model. The stages are groups, each with its own likelihood. The estimation of mixture models is, in general, mortality function, where the age-specific hazard is specified difficult, but the expectation–maximization (EM) algorithm by the estimated Weibull parameters for that particular group. (Dempster et al. 1977) has made it widely possible. The EM The state of an individual is given by its age and its het- algorithm is an iterative procedure that alternates between an erogeneity group. To include both these variables, we use expectation (E) step and a maximization (M) step, until the an age×stage-classified matrix model, in which individu- estimates converge (for details, see McLachlan and Krishnan als are jointly classified by age and stage (Caswell 2012). (2007)). Conceptually, it treats group membership as missing In this case, stages are heterogeneity groups. Each group data. In the E step, the expected value of the unknown group has an age-dependent mortality schedule specified by its membership is calculated for each individual, given the sur- Weibull parameters. In general, an age × stage-classified vival parameters in each group. The M step then finds param- model describes both progression through age classes and eters that maximize the likelihood, given the expected group transitions among stages (Caswell 2012). However, the het- memberships of each individual. Then the expectation step erogeneity groups here are fixed, so we need not include is repeated with the new parameters, and so on. See McLa- transitions among them (but see Hartemink et al. (2017), chlan and McGin ffi ( 1994) for a general discussion of the EM Caswell (2014) and Caswell et al. (2018) for more details algorithm in relation to survival analysis. We programmed on how to include such transitions). our analysis following the application of the EM algorithm Let be the number of age classes and g be the number to survival data by Mohammed et al. (2013). The approach of heterogeneity groups. The population vector ̃ is has been shown by simulation studies to be capable of distin- guishing mixtures of Weibull (and other) distributions (e.g., ⎛ ⎞ Erişoğlu et al. 2011, 2012; Mohammed et al. 2015). ⎜ ⎟ ⎜ ⎟ To help ensure convergence to global rather than local ⎜ ⎟ 1g maxima of the likelihood function, we sampled at least ten ⎜ ⎟ = ⋮ (3) ⎜ ⎟ initial values for the parameters. We selected the best esti- ⎜ 𝜔1 ⎟ mate of the number of groups; following the suggestion of ⋮ ⎟ Frühwirth-Schnatter (2006) we used the minimum Bayesian ⎜ ⎟ 𝜔g ⎝ ⎠ Information Criterion (BIC) as our criterion. This lessens the risk of overfitting the number of heterogeneity groups. Runs that resulted in values of k greater than 10 were excluded, where the jth block of entries in ̃ is a sub-vector describing because values of k > 10 produce extremely narrow distribu- the abundance of the g stage classes within age class j. tions of age at death, indicating that the model is trying to fit For each group i, define a survival matrix of dimension a few data points instead of a general distribution. i × that contains age-specific survival probabilities on the first subdiagonal and zeros elsewhere, A matrix model including heterogeneity 00 ⋯ 0 ⎛ ⎞ − (0) Notation Matrices are denoted by upper-case boldface letters ⎜ e 0 ⋯ 0 ⎟ (4) i ⎜ ⎟ (e.g., U), and vectors by lower-case boldface letters (e.g., ⋮⋱ ⋮ ⎜ ⎟ − (−1) n). Block-diagonal matrices are denoted by blackboard font ⎝ 0 ⋯ e 0 ⎠ (e.g., ). A tilde is used to distinguish matrices and vectors associated with the full age×stage-classified model, e.g., by where (x) is the mortality rate at age x, given by Eq. 1, , ̃ ; these matrices are block-structured and contain entries for group i. Create a block-diagonal matrix (of dimension for all combinations of age classes and heterogeneity groups. g × g ) by placing the U on the diagonal, The identity matrix of order s is denoted I , and 1 is a s × 1 s s vector of ones. The unit vector e is a vector with a 1 in the i ⎛ ⋯ 0 ⎞ ith entry and zeros elsewhere. The symbol ◦ denotes the ⎜ ⎟ = ⋮ ⋱ ⋮ . (5) ⎜ ⎟ Hadamard, or element-by-element product; the symbol ⊗ 0 ⋯ ⎝ ⎠ denotes the Kronecker product. The transpose of the matrix 1 3 92 Population Ecology (2018) 60:89–99 The joint age × stage composition of the cohort at time x is Variance decomposition: heterogeneity projected as and stochasticity ̃(x + 1)= ̃(x) (6) The first age class is a mixture of individuals with a mix- where the projection matrix is ing distribution (which is a vector giving the fractions of the population in each group); is estimated by the (7) EM algorithm. The variance in longevity of age class 1, with K = K the vec-permutation matrix (Henderson and g, considered as a mixture of groups, is Searle 1981; Hunter and Caswell 2005; Caswell 2012), which rearranges the population vector to permit multipli- V()= E V + V E (14) groups groups cation by the block diagonal matrix. = V + V 1 × 1. within between (15) Calculating longevity The first term is the within-group variance; it is the weighted The matrix is the transient matrix of an absorbing Markov mean of the group variances in the vector V( ), as given groups chain, with death as an absorbing state (e.g., Caswell 2001, by Eq. 13, 2009, 2014). The fundamental matrix of this chain (of dimen- sion g × g ) is V = V( ) (16) within groups −1 ̃ ̃ = − (8) 𝜔g = ⊗ V(̃) 1 × 1. (17) ̃ The second term is the between-group variance; it is where is an identity matrix. The (x, y) entry of is the the weighted variance of the group means in the vector expected number of visits to state y by an individual in state E( ), as given by Eq. 12, x, where state refers to the specific combination of age and groups stage. V = E( )◦E( ) ̃ between groups groups The statistics of longevity are calculated from (e.g., Cas- (18) well 2009). The vectors of first and second moments of lon- − E( ) groups gevity, are given by = ⊗ ̃ ◦ ⊗ ̃ g 1 g 1 ̃ = 1 g𝜔 × 1 (9) 1 1 (19) − ⊗ ̃ 1 × 1. ̃ = ̃ 2 − g𝜔 × 1 (10) 2 𝜔g 1 The within-group variance component measures the vari- ance due to individual stochasticity among individuals expe- These vectors contain the moments of the longevity of all g riencing the same group-specific mortality schedule. The age×stage combinations. The vector of mean life expectan- between-group component measures the variance due to the cies (mean longevities) of each age×stage combination is ̃ . differences in the mortality schedules among the groups. The vector of variances in longevity is In the absence of heterogeneity, the variance among group V(̃)=̃ − ̃ ◦̃ g𝜔 × 1 (11) means would be zero and all variance would be due to sto- 2 1 1 chasticity. In the absence of stochasticity, all the group vari- We are interested in the remaining longevity from the start of ances would be zero and all variance would be due to het- the cohort (age class 1), so we extract the mean and variance erogeneity. Thus the between-group variance, as a fraction of longevity at age 1 from the full vectors. Define a vector of the total, is a measure of the contribution of heterogeneity , of dimension g × 1 , that contains the longevity, at age groups to variance in longevity. 1, of individuals in each of the heterogeneity groups. The mean and variance of are groups The magnitude of heterogeneity E( )= ⊗ ̃ g × 1 (12) groups g 1 We measured the amount or magnitude of heterogeneity V( )= ⊗ V(̃) g × 1 (13) groups g in two ways. First, we consider the concentration of indi- viduals within groups. Heterogeneity is less if individuals where e is a vector of length with a 1 in the first entry and are concentrated in one or a few groups than if they are zeros elsewhere and I is an identity matrix of size g. spread out among groups in relatively equal proportions. 1 3 Population Ecology (2018) 60:89–99 93 We measure this concentration by the entropy of the mix- Table 1 Characteristics of the data sets analyzed in the paper, show- ing sample size (N), life expectancy (LE) in days, the observed vari- ing distribution ance in age at death, and the maximum observed life span in days Species strain/sex Raw mortality data H =− log (20) i i N LE Variance Max. life span i=1 C. elegans which has its maximum value H = log g when all the N2 1000 14.32 26.20 33 are equal. Because the entropy is affected by the number CLK-1 800 18.24 106.48 55 of groups as well as the distribution of individuals among DAF-2 800 30.19 224.64 62 the groups, we scale it relative to its maximum to obtain the All 2600 20.40 156.93 62 evenness, Human louse Females 400 17.08 84.31 46 J = (21) Males 400 17.17 63.59 44 log g Both sexes 800 17.12 73.36 46 which ranges from 0 (in the limit as all individuals are con- Housefly centrated in one group) to 1 (individuals equally distributed Females 3875 28.74 146.14 65 among groups). Males 4627 16.93 45.60 58 Second, we consider the magnitude of differences among Both sexes 8502 22.22 123.01 65 Anastrepha ludens the groups. Heterogeneity is less when the differences Females 363,971 28.54 198.86 163 among groups are small than when they are large. We meas- Males 487,128 30.31 196.53 155 ured the magnitude of the differences among groups by the Both sexes 851,099 29.56 198.29 163 between-group standard deviation; i.e., the square root of the A. obliqua between-group variance (Eq. 18). In order to compare spe- Females 134,807 17.91 147.07 84 cies with different life expectancies, we scaled the standard Males 162,280 15.43 100.67 67 deviation by the overall life expectancy. Both sexes 297,087 16.55 123.24 84 A. serpentina Data Females 169,031 17.96 96.72 90 Males 172,283 18.57 113.74 77 We obtained individual survival data from the literature Both sexes 341,314 18.27 105.40 90 or from the DATLife database (DATLife 2017), choosing D. longicaudata studies with large sample sizes to permit rigorous statisti- Females 14,184 8.60 38.25 70 cal analysis. All data were obtained from laboratory studies Males 13,358 7.97 30.79 64 under constant (to the best efforts of the original investiga- Both sexes 27,542 8.29 34.73 70 Drosophila tors) conditions. In the end, we analyzed 25 data sets on Long-winged females 5426 38.00 406.93 97 nine species of invertebrates (one nematode and 8 insects). Long-winged males 4568 40.40 433.58 95 Some of the data were additionally broken down by sex or Short-winged females 906 15.26 112.16 46 genetic strains. The sample sizes and characteristics of the Short-winged males 854 13.67 79.91 51 data are shown in Table 1. More detailed information on spe- All 11,754 35.42 449.17 97 cies, sources of data, and experimental conditions is given Medfly diet study in Appendix 1. Females sugar only 101,362 12.89 58.40 97 Males sugar only 106,662 14.23 56.18 79 Females sugar+protein 99,312 14.35 64.51 133 Results Males sugar+protein 108,953 14.18 51.98 64 All 416,289 13.92 57.95 133 From each data set, we obtained the estimated number g of Million medfly study groups (selected by minimizing BIC), the Weibull param- Females 605,528 19.58 87.56 171 Males 598,118 22.13 80.45 164 eters for each group, and the proportions of each group in Both sexes 1,203,646 20.84 85.65 171 the initial cohort. From these, we constructed the age×stage- classified matrix model and partitioned the variance in lon- gevity into components due to heterogeneity and individual In Electronic Supplementary Material [ESM-1], we provide stochasticity, and obtained the proportion of the variance due the complete set of all estimates, not just those for the model to heterogeneity. The resulting values are shown in Table 2. selected by minimizing BIC. 1 3 94 Population Ecology (2018) 60:89–99 Fig. 1 Age-at-death of C. elegans strain CLK-1: raw data (red) and Fig. 2 C. elegans CLK-1: fitted Weibull functions for age-at-death for modelled by a single Weibull function (black) the weighted mixture and for the two groups Detailed results for C. elegans 4.5−1 4.5 x (x)= (22) 11.9 11.9 To clarify the analyses and as an example of the procedure, we present here the detailed results for a cohort of 800 indi- 2.5−1 viduals of the CLK-1 strain of the nematode C. elegans from 2.5 x (x)= (23) an experiment described in Chen et al. (2007). Fitting a sin- 27.3 27.3 gle Weibull function to the age-at-death data of this strain These hazard functions determine the age-specific sur- yields an estimate of = 20.7 and k = 1.9. From Fig. 1, it is vival probabilities in the matrices in Eq. 4. From this, clear that a single Weibull function does not provide a good the block-diagonal matrix and the projection matrix are fit to the data. derived using Eqs. 5 and 7. The mean longevity in the het- The results of fitting mixtures of two, three, four, or five erogeneous cohort is 19.2 d. The variance is 106 d ; of this Weibull functions are shown in Table 3 (models with mix- variance, 41.4% is due to heterogeneity between the groups, tures of six, seven, or eight Weibull functions either did not and the remaining 58.6% is due to individual stochasticity. converge or produced results with k > 10). The among-group standard deviation is 35% of the mean Based on the BIC values, we conclude that a mixture longevity. of two Weibull functions is the model most well sup- The Mtalab scripts for estimating the parameters and ported by these data. The first group, comprising 45% of BIC values for each model using the EM algorithm and the individuals, is characterized by Weibull parameters for calculating longevity statistics and decomposing the = 11.9 and k = 4.5 . The other group, comprising the variance, can be found in the Electronic Supplementary remaining 55% of the individuals, has Weibull parameters Material [ESM-3]. = 27.3 and k = 2.5 . These two functions, scaled by their mixing proportions, and their mixture are shown in Fig. 2. Results: species comparison The first group is characterized by a shorter and less vari- able longevity (modal age at death 11 days), the second Table 2 shows the results for the number of heterogeneity group by a longer and more variable life span (modal age groups identified, the evenness of the distribution of individ- at death 22 days). When comparing the raw data, a single uals among groups, the magnitude of the differences among Weibull, and a mixture of two Weibull functions (Fig. 3), groups, and the fraction of variance due to heterogeneity. In it is clear that the mixture model provides the better fit. only four cases (both sexes of the human louse, the N2 strain These estimated parameters are used as input for our of C. elegans, and males of the short-winged strain of Dros- age×stage-classified matrix model. The number of groups ophila) did we fail to find evidence of heterogeneity. In the is g = 2 in this case. We used 200 age classes ( = 200 ) . other cases, populations were quite evenly distributed among The mixing distribution = 0.45 0.55 ; this nearly groups, with a median evenness (J) of 0.75 (interquartile equal division into two groups yields an evenness of 0.99. range 0.41–0.89). The age-specific mortality hazards for the two groups are 1 3 Population Ecology (2018) 60:89–99 95 Partially because of its evolutionary implications, much of the interest in unobserved heterogeneity in fitness com - ponents focuses on accounting for variance (e.g., Caswell 2011, 2014; Steiner and Tuljapurkar 2012; Vindenes and Langangen 2015; Cam et al. 2016; Hartemink et al. 2017; van Daalen and Caswell 2017; Jenouvrier et al. 2018). In this case, we found that heterogeneity could typically account for less than half of the variance in longevity (35%, with interquartile range 23–44%). We found substantial differences among species in the number of groups distinguished and in the proportion of the variance attributable to heterogeneity. However, we found no clear patterns involving differences between sexes, treatments, or strains. Application of this approach to other experimental studies, in which large numbers of individuals are exposed to different treatments, would be valuable. Fig. 3 Age-at-death of C. elegans strain CLK-1: raw data (red), mod- The very large datasets (the Anastrepha species and the elled as single Weibull function (black) and modelled as a mixture of Million Medfly experiment) seem to reveal higher numbers two Weibull functions (blue) of groups; this may reflect an increased ability to detect het- erogeneity with large sample sizes. There is no correlation The magnitude of the heterogeneity (the among-group between the number of groups and the fraction of variance standard deviation) had a median value of 28% of life expec- due to heterogeneity; both high and low numbers of groups tancy (interquartile range 24–34%). Heterogeneity accounts can result in high fraction of variance due to heterogeneity. for a substantial but not overwhelming fraction of the vari- For example, the fraction of variance due to heterogene- ance in longevity. The median contribution of heterogeneity ity was high with only two groups (e.g., 65% in the DAF-2 is 35% (interquartile range 23–44%). The highest contribu- strain of C. elegans) and with six groups (e.g., 75% in female tions are 75% in Anastrepha obliqua females and 65% in the A. obliqua). DAF-2 strain of C. elegans. Note that the value of the estimated longevity is in all cases approximately one unit higher than the longevity in the raw data, this is caused by the matrix model assumption Discussion of one remaining unit of life expectancy, even at the time of death. If we subtract this unit here, the estimated longevities We set out to address three questions: is there evidence for match the raw data very well. The estimated variance also heterogeneity, if so how much, and what fraction of the vari- closely matches the variance as calculated from the raw data. ance in longevity is due to heterogeneity, and what fraction There exist a few studies to which this one can be com- to stochasticity. We found statistical support for heterogene- pared. Heterogeneity makes a larger contribution to vari- ity in 31 out of 35 cases. In only four data sets (the N2 strain ance in longevity in this study than was previously found for of C. elegans, both sexes of the human louse, and males of humans in an analysis of cohort and period mortality pat- the short-winged strain of Drosophila) was a homogeneous terns, over many years, for populations of Sweden, France, model, with only a single group, the best supported. These and Italy (Hartemink et al. 2017). In that analysis, heteroge- were among the smallest data sets in our studies (1,000 indi- neity accounted for less than 10%, and usually less than 5%, viduals for C. elegans, 400 of each sex for the louse, and 854 of the variance in longevity. It was based on a continuous for the short-winged Drosophila). It would not be surprising heterogeneity model, in which a Gamma-distributed frailty if heterogeneity is more difficult to detect in small samples. term, acting as a proportional hazard on mortality, was Had these experiments been performed with more individu- applied to a Gompertz–Makeham mortality model (Missov als, multiple groups might have been identified. 2013; Missov and Lenart 2013). The Gompertz–Makeham For all other data sets, our analysis identified from 2 to model is applicable to human populations only after the age 6 heterogeneity groups. Generalizing from these results, of 30–40 years, so the human data corresponded to a later we can say that individuals are relatively evenly spread out “adult” age than is the case for the invertebrate species stud- among the groups, with an evenness of about 75% of its ied here. maximum. The differences among groups in life expectancy Caswell (2014) made a brief exploration of laboratory are about 28% of overall life expectancy. data on six species of invertebrates using Gamma–Gompertz parameters reported by Horiuchi (2003). A crude, 1 3 96 Population Ecology (2018) 60:89–99 Table 2 Results of mixture model analysis Species strain/sex Best fitting mixture of Weibull functions g LE Variance Within Between Entropy Evenness Ratio % of variance C. elegans N2 1 15.33 26.82 26.82 0.00 0.00 – 0.0 0.0 CLK-1 2 19.22 106.45 62.38 44.07 0.69 0.99 0.35 41.4 DAF-2 2 31.21 223.75 78.90 144.85 0.61 0.88 0.39 64.7 All 4 21.42 156.78 46.83 109.95 1.11 0.80 0.39 70.1 Human louse Females 1 18.06 85.28 85.28 0.00 0.00 – 0.0 0.0 Males 1 18.14 62.76 62.76 0.00 0.00 – 0.0 0.0 Both sexes 1 18.10 74.27 74.27 0.00 0.00 – 0.0 0.0 Housefly Females 2 29.73 146.83 91.01 55.83 0.44 0.64 0.25 38.0 Males 4 17.94 45.84 32.95 12.89 1.10 0.79 0.20 28.1 Both sexes 4 23.30 125.81 70.11 55.70 0.97 0.70 0.32 44.3 Anastrepha ludens Females 6 29.55 198.62 126.76 71.86 1.44 0.80 0.29 36.2 Males 5 31.32 196.35 158.19 38.16 0.81 0.50 0.20 19.4 Both sexes 5 30.56 198.08 160.74 37.34 1.04 0.65 0.20 18.9 A. obliqua Females 6 18.93 146.65 37.29 109.35 1.72 0.96 0.55 74.6 Males 6 16.45 100.38 32.12 68.26 1.67 0.93 0.50 68.0 Both sexes 6 17.59 122.77 69.77 53.01 1.28 0.71 0.41 43.2 A. serpentina Females 5 18.97 96.47 78.43 18.04 0.80 0.50 0.22 18.7 Males 5 19.59 113.40 82.87 30.53 1.04 0.65 0.28 26.9 Both sexes 5 19.29 105.07 82.09 22.98 1.14 0.71 0.25 21.9 D. longicaudata Females 4 9.64 38.07 27.70 10.37 1.31 0.95 0.33 27.2 Males 4 9.05 30.45 25.77 4.69 0.71 0.51 0.24 15.4 Both sexes 4 9.36 34.36 28.95 5.41 0.85 0.61 0.25 15.7 Drosophila Long-winged females 2 39.06 403.70 240.11 163.59 0.57 0.82 0.33 40.5 Long-winged males 2 41.43 431.64 239.23 192.41 0.52 0.76 0.33 44.6 Short-winged females 2 16.28 111.82 50.83 60.99 0.65 0.93 0.48 54.5 Short-winged males 1 14.76 77.23 77.23 0.00 0.00 – 0.00 0.0 All 2 36.46 446.91 213.69 233.22 0.64 0.92 0.42 52.2 Medfly diet study Females sugar only 4 13.94 58.10 44.85 13.25 1.24 0.89 0.26 22.8 Males sugar only 4 15.27 55.88 50.47 5.41 1.05 0.75 0.15 9.7 Females sugar + protein 4 15.40 64.04 44.80 19.23 1.31 0.94 0.28 30.0 Males sugar + protein 4 15.23 51.60 38.72 12.88 1.30 0.94 0.24 25.0 All 4 14.97 57.55 44.06 13.49 1.38 1.00 0.25 23.4 Million medfly study Females 6 20.58 87.58 53.13 34.44 1.37 0.76 0.29 39.3 Males 8 23.14 80.47 53.00 27.47 1.71 0.82 0.23 34.1 Both sexes 8 21.85 85.66 50.10 35.56 1.56 0.75 0.27 41.5 The number of groups (g) that results in the lowest BIC, the life expectancy (LE), the total variance in longevity, the within- and between-class variance, the entropy (H) and evenness (J) indices, the ratio between the square root of the between variance and the life expectancy, and the percentage of variance due to heterogeneity 1 3 Population Ecology (2018) 60:89–99 97 Table 3 Best fits for models g Group 1 Group 2 Group 3 Group 4 Group 5 ΔBIC using mixtures of up to five Weibull functions for the CLK-1 k % k % k % k % k % strain of C. elegans 1 20.7 1.9 100 167.23 2 11.9 4.5 45.0 27.3 2.5 55.0 0 3 12.0 4.3 49.2 0.8 1.9 0.3 28.6 2.8 50.5 8.51 4 19.6 4.6 15.7 1.2 1.5 0.4 32.8 3.7 34.0 11.8 4.4 0.50 20.80 5 37.1 4.5 14.8 20.7 4.1 25.2 11.8 4.4 49.5 33.5 7.3 10.0 1.2 1.5 0.4 40.12 For each model the number of groups (g), the estimated Weibull parameters and k, the percentage of the population made up by each group are shown. Also shown is ΔBIC, the difference from the minimum value of BIC, which corresponds to the best model approximate, variance decomposition found about 60% of 2014). In our approach, the number of groups, the distri- the variance in longevity to be due to heterogeneity. How- bution of individuals among the groups, and the scale and ever, because it is now known that neglecting the Makeham shape parameters of each of the Weibull functions are esti- mortality term can bias estimates of the Gompertz param- mated without restriction. eters (Missov and Németh 2016), and because of the ad hoc An important direction for future research is the incorpo- variance decomposition used in Caswell (2014), we view ration of dynamic heterogeneity, in which group member- these results as suggestive but not reliable. ship is not fixed over the life of the individual. It will not be In a recent field study of the Southern fulmar (Fulmarus easy to estimate unobserved heterogeneity in these models glacialoides), Jenouvrier et al. (2018) used multievent cap- (e.g., Putter and van Houwelingen 2015). However, when ture-recapture analysis to identify three unobserved hetero- the heterogeneity can be observed or measured, dynamic geneity groups, where the groups were allowed to differ in transitions may be incorporated following the methods of any of the transition probabilities in a stage-classified matrix multistate event history analysis (Willekens 2014). population model. They decomposed variance in longevity, Demographic components of fitness (longevity, lifetime age at first breeding, and lifetime reproductive output into reproductive output, age at first breeding, etc.) are important contributions from heterogeneity and stochasticity. Hetero- components of evolutionary demography. The variance in geneity accounted for only 5.9% of the variance in longevity, these components, if due to genetic heterogeneity, would 3.7% of the variance in age at first breeding, and 22% of the provide material for natural selection. The results presented variance in lifetime reproductive output. here, and those for recent human populations (Hartemink Insects undergo metamorphosis before reaching the adult et al. 2017) and one long-lived seabird (Jenouvrier et al. stage, and in some laboratory conditions (e.g., Drosophila 2018), suggest that for longevity, most of the variance is culture bottles), crowded larval conditions may create het- due to individual stochasticity. An understanding of the fac- erogeneity through competition, with some individuals com- tors that influence this proportion, and the patterns shown pleting the larval stage, but less well equipped for the adult by other taxa, are important research questions. stage, leading to ‘early failure’ . Such early deaths would Acknowledgements We thank Mijke Rhemtulla for suggestions about probably form a group with very low mean longevity, and statistical methods, Roberto Salguero-Gomez and an anonymous this would contribute substantially to the variance, and also reviewer for helpful comments, and we particularly thank Alexander increase the contribution of heterogeneity to the variance Scheuerlein for assistance with the DATLife database. This research in longevity. was supported by the European Research Council under the European Union’s Seventh Framework Program, ERC Advanced Grant 322989. The finite mixture approach to estimating heterogeneity has advantages over other frailty analyses. It does not require Open Access This article is distributed under the terms of the Crea- an assumption of a parametric mixing distribution (e.g., the tive Commons Attribution 4.0 International License (http://creat iveco Gamma distribution), nor does it require an assumption of mmons.or g/licenses/b y/4.0/), which permits unrestricted use, distribu- how the heterogeneity acts. In the Gamma–Gompertz–Make- tion, and reproduction in any medium, provided you give appropriate credit to the original author(s) and the source, provide a link to the ham model, for instance, heterogeneity acts as a proportional Creative Commons license, and indicate if changes were made. factor multiplying a baseline hazard (Vaupel and Missov Early failure is a concept coined in reliability engineering (Barlow and Proschan1996) for machines or parts that break down very soon after the onset of the use, due to some mechanical failure. 1 3 98 Population Ecology (2018) 60:89–99 below. The data were obtained from the DATLife data- Appendix: Species and data base (DATLife 2017). 6. Anastrepha serpentina The sapote or serpentine fruit This appendix gives details of the species and experimen- fly (Insecta:Diptera:Tephritoidea) is a major pest in tal sources for the data analyzed. Detailed comparisons of Mexico. We used daily mortality data from a longevity the raw data and the estimated Weibull functions are given experiment on 172,283 male and 169,031 female flies. for all data sets in ESM-2. The study conditions were identical to Million medfly study described below. The data were obtained from 1. Caenorhabditis elegans This is a free-living nematode the DATLife database (DATLife 2017). species and an often-used model organism. Mortal- 7. Diachasmimorpha longicaudata This is a solitary parasi- ity data on C. elegans was obtained from a experi- toid wasp (Insecta:Hymenoptera:Braconidae). It is a para- ment by Chen et al. (2007). In this experiment, the sitoid of Caribbean fruit fly larvae. Daily mortality data longevity of a standard wild-type strain (N2) and two were obtained from a longevity experiment on 13,358 long-lived mutant strains (CLK-1 and DAF-2) were male and 14,184 female wasps. The study conditions reported. Daily survival data on 1000 individuals of were essentially identical to the Million medfly study the N2 strain and 800 individuals of each of the mutant conducted by J. Carey and described below. Data were strains were available. obtained from the DATLife database (DATLife 2017). 2. Pediculus humanus (Human louse) The human louse 8. Drosophila melanogaster The common fruit fly (Pediculus humanus L.) is a small, wingless insect. (Insecta:Diptera:Drosophilidae) is probably the most Human lice are obligate parasites of humans, that is, widely studied laboratory organism. Life tables for they normally feed exclusively on human blood, but they this species were obtained from an early publication can be reared successfully in the lab on rabbits. The life by Pearl and Parker (1921) on longevity experiments cycle consists of the egg, three larval instars and the on long-winged and short-winged (Quintuple stock) adult stage. Data on adult survival of head lice were strains. Longevity was measured for 4586 long-winged obtained from a study by Evans and Smith (1952), in males, 5426 long-winged females, 854 short-winged which 800 freshly emerged adults (400 males and 400 males and 906 short-winged females. females) were kept in mixed colonies. The mean dura- 9. Medfly caloric restriction experiment, Cerati- tion of life was 17.6 for both males and females; there tis capitata. The Mediterranean fruit fly (often was no significant difference between the two sexes. called Medfly for short), is a species of fruit fly 3. Musca domestica (Common House fly) The common (Insecta:Diptera:Tephritidae) and an important fruit house fly (Musca domestica L.) is a well-known fly pest. It is native to the Mediterranean area, but has species (Insecta:Diptera:Muscidae). We use data on spread invasively to many parts of the world. We used adult life span from life tables for 4,627 males and data from an experiment on the effect of caloric restric- 3,875 females published in Rockstein and Lieberman tion (Müller et al. 1997). Daily mortality of 200,674 (1959). The data originally came from 2 separate stud- males and 215,615 females, maintained in grouped ies, but no difference in mean length of life was found cages, was observed for two caloric restriction groups: between these two datasets, so we treated it as one sugar and sugar plus protein. The data were obtained dataset. The house flies belonged to the strain NAIDM from the DATLife database (DATLife 2017). and the flies were inbred for about 200 generations. 10. The Million medfly experiment, Ceratitis capitata. The 4. Anastrepha ludens Also known as the Mexican fruit fly Million Medfly dataset consists of longevity data for (Insecta:Diptera:Tephritoidea) this is a major pest of very large cohorts of Ceratitis capitata. The experiment agriculture across the Americas. We used daily mor- was performed in 1991 at the Moscamed medi fl es mass- tality data from a longevity experiment conducted on rearing facility. The purpose of this study was to exam- 487,128 male flies and 363,971 female flies (Vaupel ine mortality at the extreme ages. Approximately 7200 et al. 1998). The study conditions were identical to the medflies (both sexes) were maintained in each of the Million medfly study described below. The data were 167 cages. Adults were given a diet of sugar and water, obtained from the DATLife database (DATLife 2017). ad libitum. Parts of these data were originally published 5. Anastrepha obliqua The West Indian or Antillian fruit in Carey et al. (1992) and Carey (1993). The dataset fruit fly (Insecta:Diptera:Tephritoidea) is a major used here contains information on the daily numbers of pest of mangoes. We used daily mortality data from age-cage-and-sex-specific deaths of the total 1,203,646 a longevity experiment on 162,280 male and 134,807 medflies (598,118 males and 605,528 females). The data female flies (Vaupel et al. 1998). The study conditions were obtained from the DATLife database (DATLife were identical to the Million medfly study described 2017). 1 3 Population Ecology (2018) 60:89–99 99 Kendall BE, Fox GA (2002) Variation among individuals and reduced References demographic stochasticity. Conserv Biol 16:109–116 Kendall BE, Fox GA, Fujiwara M, Nogeire TM (2011) Demographic Barlow RE, Proschan F (1996) Mathematical theory of reliability. SIAM, heterogeneity, cohort selection, and population growth. Ecology Philadelphia 92:1985–1993 Bijwaard GE (2014) Multistate event history analysis with frailty. Demogr McLachlan G, Krishnan T (2007) The EM algorithm and extensions, vol Res 30:1591–1620 382, 2nd edn. Wiley, New York Cam E, Aubry LM, Authier M (2016) The conundrum of heterogeneities McLachlan GJ, McGiffn DC (1994) On the role of finite mixture models in life history studies. Trends Ecol Evol 31:872–886 in survival analysis. Stat Methods Med Res 3:211–226 Carey JR (1993) Applied demography for biologists: with special empha- McLachlan G, Peel D (2004) Finite mixture models. Wiley, New York sis on insects. Oxford University Press, Oxford Mills M (2011) Introducing survival and event history analysis. Sage, Carey JR, Liedo P, Orozco D, Vaupel JW (1992) Slowing of mortality London rates at older ages in large medfly cohorts. Science 258:457–491 Missov TI (2013) Gamma–Gompertz life expectancy at birth. Demogr Caswell H (2001) Matrix population models: construction, analysis, and Res 28:259–270 interpretation, 2nd edn. Sinauer Associates, Sunderland Missov TI, Lenart A (2013) Gompertz–Makeham life expectancies: Caswell H (2009) Stage, age and individual stochasticity in demography. expressions and applications. Theor Popul Biol 90:29–35 Oikos 118:1763–1782 Missov TI, Nemeth L (2016) Sensitivity of model-based human mortal- Caswell H (2011) Beyond R : demographic models for variability of life- ity measures to exclusion of the Makeham or the frailty parameter. time reproductive output. PLoS One 6(6):e20809 Genus 71:113–135 Caswell H (2012) Matrix models and sensitivity analysis of populations Mohammed YA, Yatim B, Ismail S (2013) A simulation study of a para- classified by age and stage: a vecpermutation matrix approach. metric mixture model of three die ff rent distributions to analyze het - Theor Ecol 5:403–417 erogeneous survival data. Mod Appl Sci 7:1–9 Caswell H (2014) A matrix approach to the statistics of longevity in het- Mohammed YA, Yatim B, Ismail S (2015) Mixture model of the expo- erogeneous frailty models. Demogr Res 31:553–592 nential, gamma and Weibull distributions to analyse heterogeneous Chen J, Senturk D, Wang JL, Müller HG, Carey JR, Caswell H, Cas- survival data. J Sci Res Rep 5:132–139 well-Chen EP (2007) A demographic analysis of the fitness cost of Müller HG, Wang JL, Capra WB, Liedo P, Carey JR (1997) Early mor- extended longevity in Caenorhabditis elegans. J Gerontol Ser A Biol tality surge in protein-deprived females causes reversal of sex dif- Sci Med Sci 62:126–135 ferential of life expectancy in Mediterranean fruit flies. Proc Natl Caswell H, de Vries C, Hartemink N, Roth G, van Daalen S (2018) Age Acad Sci USA 94:2762–2765 × stage-classified demographic analysis: a comprehensive approach. Pearl R, Parker SL (1921) Experimental studies on the duration of life. Ecol Monogr. https ://doi.org/10.1002/ecm.1306 I. Introductory discussion of the duration of life in Drosophila. Am DATLife (2017) Database: demography of aging across the tree of life. Nat 55:481–509 Max Planck Institute for Demographic Research. https ://datli fe.org Pinder JE, Wiener JG, Smith MH (1978) The Weibull distribution: a new Dempster AP, Laird NM, Rubin DB (1977) Maximum likelihood from method of summarizing survivorship data. Ecology 59:175–179 incomplete data via the EM algorithm. J Roy Stat Soc Ser B (Meth- Putter H, van Houwelingen HC (2015) Frailties in multi-state models: odol) 39:1–38 are they identifiable? Do we need them? Stat Methods Med Res Erişoğlu U, Erişoğlu M, Erol H (2011) A mixture model of two different 24:675–692 distributions approach to the analysis of heterogeneous survival data. Rinne H (2008) The Weibull distribution: a handbook. CRC, Boca Raton Int J Comput Math Sci 5:544–548 Robert A, Sarrazin F, Couvet D (2003) Variation among individuals, Erişoğlu U, Erişoğlu M, Erol H (2012) Mixture model approach to the demographic stochasticity, and extinction: response to Kendall and analysis of heterogeneous survival data. Pak J Stat 28:115–130 Fox. Conserv Biol 17:1166–1169 Evans FC, Smith FE (1952) The intrinsic rate of natural increase for the Rockstein M, Lieberman HM (1959) A life table for the common house- human louse, Pediculus humanus L. Am Nat 86:299–310 fly, Musca domestica. Gerontology 3:23–36 Farewell VT (1982) The use of mixture models for the analysis of survival Steiner UK, Tuljapurkar S (2012) Neutral theory for life histories and data with long-term survivors. Biometrics 38:1041–1046 individual variability in fitness components. Proc Natl Acad Sci Frühwirth-Schnatter S (2006) Finite mixture and Markov switching mod- USA 109:4684–4689 els. Springer, New York van Daalen SF, Caswell H (2017) Lifetime reproductive output: indi- Hartemink N, Missov TI, Caswell H (2017) Stochasticity, heterogeneity, vidual stochasticity, variance, and sensitivity analysis. Theor Ecol and variance in longevity in human populations. Theor Pop Biol 10:355–374 114:107–117 Vaupel JW, Carey JR (1993) Compositional interpretations of medfly Heckman JJ, Singer B (1982) Population heterogeneity in demographic mortality. Science 260:1666 models. In: Land KC, Rogers A (eds) Multidimensional mathemati- Vaupel JW, Missov TI (2014) Unobserved population heterogeneity: a cal demography. Academic, New York, pp 567–599 review of formal relationships. Demogr Res 31:659–686 Henderson HV, Searle SR (1981) The vec-permutation matrix, the vec Vaupel JW, Manton KG, Stallard E (1979) The impact of heterogene- operator and Kronecker products: a review. Linear Multilinear Alge- ity in individual frailty on the dynamics of mortality. Demography bra 9:271–288 16:439–454 Horiuchi S (2003) Interspecies differences in the life span distribution: Vaupel JW, Carey JR, Christensen K, Johnson TE, Yashin AI, Holm humans versus invertebrates. Popul Dev Rev 29:127–151 NV, Iachine IA, Kannisto V, Khazaeli AA, Liedo P, Longo VD, Horvath WJ (1968) A statistical model for the duration of wars and Zeng Y, Manton KG, Curtsinger JW (1998) Biodemographic tra- strikes. Syst Res Behav Sci 13:18–28 jectories of longevity. Science 280:855–860 Hunter CM, Caswell H (2005) The use of the vecpermutation matrix in Vindenes Y, Langangen Ø (2015) Individual heterogeneity in life histo- spatial matrix population models. Ecol Model 188:15–21 ries and eco-evolutionary dynamics. Ecol Lett 18:417–432 Jenouvrier S, Aubry LM, Barbraud C, Weimerskirch H, Caswell H (2018) Willekens F (2014) Multistate analysis of life histories with R. Interacting effects of unobserved heterogeneity and individual sto- Springer, New York chasticity in the life cycle of the Southern fulmar. J Anim Ecol 87:212–222 1 3
Population Ecology – Springer Journals
Published: May 28, 2018
It’s your single place to instantly
discover and read the research
that matters to you.
Enjoy affordable access to
over 18 million articles from more than
15,000 peer-reviewed journals.
All for just $49/month
Query the DeepDyve database, plus search all of PubMed and Google Scholar seamlessly
Save any article or search result from DeepDyve, PubMed, and Google Scholar... all in one place.
All the latest content is available, no embargo periods.
“Whoa! It’s like Spotify but for academic articles.”@Phil_Robichaud