The equation as a predictive tool in demography.

The Gompertz demographic model describes rates of aging and age-independent mortality with the parameters alpha and A, respectively. Estimates of these parameters have traditionally been based on the assumption that mortality rates are constant over short to moderate time periods. This assumption is questionable even for very large samples assayed over short time intervals. In this article, we compare several methods for estimating the Gompertz parameters, including some that do not assume constant mortality rates. A maximum likelihood method that does not assume constant mortality rates is shown to be best, based on the bias and variance of the Gompertz parameter estimates. Moreover, we show how the Gompertz equation can then be used to predict mean longevity and the time of the nth percentile of mortality. Methods are also developed that assign confidence intervals to such estimates. In some cases, these statistics may be estimated accurately from only the early deaths of a large cohort, thus providing an opportunity to estimate longevity on long-lived organisms quickly.


INTRODUCTION
THERE ARE compelling reasons to believe that a single phenomenon, the decline in the force of natural selection, is the cause of aging in plants and animals (Rose, 1991). But prior to the application of evolutionary theory to the problem of aging, phenomenological models were developed that attempted to capture the observation of increasing rates of mortality with age. One such model is that due to Benjamin Gompertz (1825), the "Gompertz model," which assumes that mortality rates increase exponentially with age. Although evolutionary theory shows that rates of mortality will increase with age in most multicellular organisms (Charlesworth, 1980), there is no deductive demonstration for the exponential form of the Gompertz model. Nevertheless, the Gompertz model has proven popular as a general description of many observed patterns of mortality (e.g., Finch, 1990).
Recently, several large studies with medflies (Carey et al., 1992) and Drosophila (Curtsinger et al., 1992) have concluded that rates of mortality may not increase exponentially in the oldest age classes, contrary to the Gompertz model. Brooks, Lithgow, and Johnson (1994) collected data that suggest that this phenomenon may be due, in part, to the existence of genetic or phenotypic subpopulations with very different rates of aging. Despite this uncertainty, there are still potentially useful applications of the Gompertz equation. In this article, we explore some practical issues concerning the use of the Gompertz equation to infer demographic parameters. We test these methods using a large array of data we have collected from populations of Drosophila melanogaster. For many species, it may be quite difficult and time consuming to follow a cohort through their entire life span in order to estimate the mean longevity, maximum longevity, and the parameters of the Gompertz equation. We show how one may use information from deaths early in life to estimate the mean longevity and the time to which the last 5 or 1% of the population is expected to survive. In addition, we develop techniques for placing confidence intervals on these estimates, and we test their accuracy using computer situations.

Gompertz equation
Mortality data is typically collected by following an initial cohort of N o individuals over time and recording the number of survivors, Nt, at fixed intervals. The Gompertz equation that describes the instantaneous rate of mortality in these types of experiments is, Consequently, the number of survivors at any time in the future can be described with eq. (1) as,

(3)
Pt, then, is the probability of any single individual surviving to age t. Let the observation vector of a single experiment be T = (tl, t2 ..... tN0), where t i is the time of death of the ith individual. GOMPERTZ STATISTICS 555

Methods of estimation
There are several approaches that have been taken to estimating the values of A and a from eq. (1). Here, we briefly review four of these methods.

(4)
To estimate the dependent variables in (4), observations of mortality are divided into discrete intervals. The number of deaths in each interval is used as an estimate of u(t). For instance, to estimate u(t) from the number of survivors the following approximation is generally used, The length of the time intervals used must be large enough to include at least one observed death. This technique implicitly assumes that u(t) does not change in the interval t to t + 1 which, of course, is contrary to the form of eq. (1). Such an approximation can only be expected to work if the time interval is very small.
Nonlinear regression. This method will find values of A and a, as the parameters that minimize the residual sum of squares below (Gavrilov and Gavrilova, 1991 ;Hirsh et al., 1994) which can be found using standard numerical search algorithms (Marquardt, 1963). The last death is not utilized in the estimation process, because eq. (6) is undefined in that single case. Using natural logarithms reduces the chance of obtaining numerical overflows during the search of the A -a space for the least-squares solutions because the logarithm of Pt is invariably a smaller number (either in absolute value or in its exponent) than Pt.
Maximum likelihood with approximation. To describe the likelihood function (Fukui et al., 1993) we first note that the last death among a cohort of N O individuals is on day tNo. Thus, we can also record the number of deaths on each day as, dl, d2 ..... dtNo.
Let q(t) be the probability that an individual that lived to day t, dies by time t + 1, then the likelihood function can then be defined as, Equation (7) is then evaluated assuming that -In(1 -q(t)) = Aexp(e~t), which is formally the same as eq. (5). The likelihood function follows from the assumption that the number of deaths (di) on day i is a binomial random variable with probability q(i).
Maximum likelihood. The most accurate application of the maximum likelihood technique would be to use eq. (7), estimating q(t) as 1 -Pt+l (Pt)-l from eq. (3).

Important results derived from the Gompertz equation
Because Pt is the chance of surviving to time t, the chance of an individual dying before time t is simply, Prob(t i~< t) = 1 -pt

(8)
If we view the time of death as a random process, then 1 -Pt is the distribution function (which by definition is Prob(ti ~< t)) of this process. Further, the density function can be derived as, The Gompertz equation then yields as an estimate of the mean longevity, fo x tf(t)dt.
In addition to the mean longevity, the time (t~) by which some critical fraction ([3) of the population is expected to be dead may be calculated from eqs.

Bootstrap
The bootstrap is a numerically intensive technique for making statistical inferences (for a recent review, see Efron and Tibshirani, 1993). It is usually reserved for problems that are not readily solved by traditional statistical theory. The success of the bootstrap depends on the recognition that the actual sample may be used as an empirical estimate of the true distribution of the underlying random variables. From this empirical distribution, multiple samples may be drawn and the behavior of the statistics of interest observed in these samples. A crucial, but nontrivial, component of the bootstrap is the use of a proper sampling scheme from the original observations. For the mortality data used in this article, there are at least three different ways of sampling: sampling individuals, sampling residuals, and sampling the stochastic process. Because preliminary work suggested that only the last of these methods worked well, we discuss it in more detail below.
Sampling the stochastic process. This process involves using the estimated values of A and a from the original data to estimate the probability of surviving to each census point. For the data analyzed here, populations are checked at daily intervals for deaths. The chance of an individual dying between day x and day x + 1 is given by, and the "hat" indicates that estimates of A and ot have been substituted in eq. (3). The original No individuals are then followed through time. On day 1, a random number with a uniform distribution on (0,1) is chosen for each of the No individuals. If the number is less than Ul, that individual dies; otherwise, it lives. The process is repeated for the survivors at day 2, 3 .... until No deaths are recorded. The whole process of generating N O deaths is repeated m times. The behavior of these m replicate data sets can then be used to make statistical inferences such as confidence intervals.
Confidence intervals. The "BCa" method (Efron and Tibshirani, 1993) is used for generating confidence intervals. With this method, the m values of the statistic of interest are ordered from smallest to largest. The lower 100oqth percentile of these order statistics is used for the lower confidence limit and the 100a2th percentile for the upper limit. Where, cq = ~ zo + 1 -gt(~o + z('~))j ' °L2 = dp zo + 1 --~(~ +-~-'~)) ~(.) is the standard-normal cumulative distribution function, and z (~) is the lO0oLth percentile of a standard normal distribution. If × is the fraction of all bootstrap estimates of the statistic that are less than the original estimate, then, ~o = q~-I(X).
Efron and Tibshirani call c~ the acceleration parameter, which is estimated from the jackknife. Under the assumption that the statistic of interest has a normal distribution, the variance of that statistic will be the same no matter what the true value is. In practice, the variance may, in fact, change for different population values of the statistic and require correction. The acceleration parameter helps to achieve this correction.
The acceleration parameter is computed from the following relationship, 558 L.D. MUELLER et al.
where s; is the ith pseudovalue computed by deleting the time of death for the ith individual from T and finding the maximum likelihood estimates of A and o~ from eq. (7). Because there are No points in T, there are N o pseudovalues whose mean is ~.

Simulated data sets
To test the accuracy of the estimates provided by eqs. (10) and (11), as well as the bootstrap confidence intervals, computer simulations were conducted. The first step in these simulations was to generate samples from a simulated population with Gompertz mortality rates. These samples consisted of the recorded day of death for each of N individuals (where N = 50, 100, or 150). We assumed that the population was checked at daily intervals for new deaths. For a fixed combination of A and o~, we computed the probability of surviving to time t as Pt, from (3). Suppose there are N t surviving adults at this time. The chance of surviving to time t + 1 is Pt+~. Thus, the probability of dying, u t, between day t and day t + 1 is 1 -(Pt+ 1/Pt). For each of the N t survivors, a uniform random number from (0,1) was chosen. If the number was less than ut, then the individual died on day t + 1; otherwise, the individual lived. This process was repeated until all N adults had died. For each of the three sample sizes, 1000 artificial data sets were created by following the sampling scheme just described.
With each of the 1000 data sets, the techniques described in the previous section were used to estimate A and ~. From these estimates of A and o~, the mean longevity, to. 5, t0.7, to.9, and t0.95 were calculated.
In general, we would like to minimize the variance and bias of estimators. We are also interested in investigating the ability of the Gompertz model to predict these statistical quantities, when only a fraction of the data set has been used. As a means of studying this problem, only portions of the artificial data sets were used to carry out the estimation process. In particular, we utilized the first 20, 35, and 50% of each of the data sets to estimate A and ~.

Empirical data sets
To illustrate the methods described here and determine if there are differences in the ability of each demographic model to describe real data, we have applied these methods to 20 cohorts of Drosophila melanogaster. These populations have been subjected to age-specific selection (Rose, 1984) and are of two treatment types: O, selected for late life reproduction; and B, selected for early life reproduction. There are five replicates of each type (e.g., B 1, B 2 ..... B 5) and data has been collected separately for males and females. Previously we have shown that the Bs and Os differ in demographic parameters, including their "rates of aging," cx (Nusbaum et al., submitted).

Numerical methods
Numerical integration of eq. (10) was performed by a modified Romberg integration procedure, using the Pascal procedure QROMO described in Press et al. (1986). The upper limit to this integral must be set to reasonable limits; otherwise, calls to the Gompertz density functions may result in numerical overflows. One convenient way to estimate a reasonable upper limit is to use eq. (11) and set t3 to 0.9999. In general, whatever limit is chosen must also be used in any bootstrap procedure to avoid systematic bias. The maximum of the likelihood function was found numerically with a downhill simplex method (Press et al., 1986). This method does not require finding the derivatives of eq. (7), which would be an onerous task for the full Gompertz equation.
Random numbers were calculated with the RAN3 function (Press et al., 1986). This function uses two calls to the random number generator for each number produced. The first call picks a random cell in an array of random numbers, while the second call replaces the random number chosen by the first call. Although slower than a procedure that simply generates a single number per call, RAN3 effectively increases its cycle time by the square of the cycle time of the source generator. Of course, the longer it takes the random number generator to cycle, the less likely it is that cycles will induce correlations in the observations. The normal distribution function was approximated by the polynomial expansion in Johnson and Kotz (1970).

The instantaneous death rate approximation
As mentioned earlier, a number of studies (Tatar et al., 1993;Brooks et al., 1994;Tatar and Carey 1994) utilize the approximation in eq. (5). Using the notation developed in this article, it follows that, Pt-NoPt +l Nop____~t -exp{ A [(1-exp(a(t + 1))-(1-exp(ott))]} After some algebra, one can show that, Ot The right side of eq. (12) will equal the Gompertz equation (1) only if we use the approximation that e ~ ~-1 + et.
More generally the probability of surviving over a single time interval is sometimes estimated from the geometric average of observations over many time intervals, e.g., ff = (etPt+ 1 • • • Pt+k-1) l/k. In this case, it can be shown that the approximation u(t) = -In(if) depends on the accuracy of the approximation e k~ ~ I + ktx. Clearly, the accuracy of these approximations will depend on the time interval used and the value of or.
For instance, many of the Drosophila populations examined in this article have an close to 0.1. In this case, when eq. (5) is used on daily observations of mortality, the difference between the approximation and the true value is only 0.5% of the true value. However, if eq. (5) had been used on weekly data, the difference between the true value and the approximation would be 16% of the true value. The reason this problem arises in the first place is that the Gompertz model assumes that mortality rates change 560 L.D. MUELLER et al. continuously with time, while methods for estimating mortality from survival data can only do so over finite time intervals. The larger the time interval, the worse the approximation. Because there are essentially no problems that can't be solved with the full Gompertz model without this approximation from eq. (5), we see little value in using it given the potentially large errors that can result.

Methods of estimation
Four different methods were described above for estimating the two parameters of the Gompertz equation. Here, we examine the ability of each method to provide accurate estimates of A and a. To do this, we created 1000 computer-generated data sets from the Gompertz equation with A = 0.00113 and et = 0.0928. These values were chosen to be close to values we have observed in our Drosophila populations. For each of these data sets, the values of A and et were estimated from each method. For each technique we then determined the bias, variance, and mean square error (= variance + bias 2) (Bickel and Doksum, 1977).
The results (Tables 1-2) show that in general the maximum likelihood method was best, and the linear regression technique the worst. Thus, the maximum likelihood method without approximation had smaller variance and bias (and, thus, smaller mean squared error) relative to all other methods. Using eq. (5) with maximum likelihood, even when the time intervals were only one day, increased both the variance and bias of the estimator relative to the maximum likelihood without this approximation.

Bootstrap confidence intervals
Utilizing the data sets described in the previous section we have tested the accuracy of 95% confidence intervals produced from the bootstrap and maximum likelihood techniques. For each sample, a bootstrap confidence interval was created for mean longevity, tso, t7o, t90, and t95. We then noted if the interval included the true value or not. The results of these simulations (Table 3) show that the confidence levels of these intervals are very close to the nominal value of 95%.

Simulation results
We also simulated a population with A and et values similar to those in B populations. With the simulated data sets, we estimated the statistics of interest with only the first 20, 35, and 50% of the data. The effects of using different fractions of these data sets and different sample sizes of individuals are shown in Fig. 1 for the variance in mean longevity, t0. 5, /0.7, to.9, and t0.95. The variance of these estimators clearly decreases with increasing sample size and increasing fractions of the data (Fig. 1). But the shape of these curves is, in fact, very similar for each statistic, suggesting a common numerical response to the size of the data sets used. Empirically, we can estimate from Fig. 1 that the variance of these estimators is decreasing in proportion to N, if N is the total sample size.
The effects of the fraction of the data used are more complicated. For instance,

Fraction of Data Used
FIG. 1. The variance of several demographic statistics as a function of (i) the initial size of the cohort and (ii) the fraction of all deaths that were observed prior to estimation of the Gompertz parameters.
inspection of Fig. 1 shows that the drop in variance in mean longevity from 20 to 35% is about the same for each sample size (50,100,150). In addition, the decrease in variance from 35 to 50% of the data is about the same as the decrease seen from 20 to 35%. In general, the variance depends on both the fraction of data used and the sample size.

Application to B and 0 data
The Gompertz equation was used with mortality data (Figs. 2-3) from the B and O populations to predict the mean longevity and to.95 (Tables 4-5). In addition, bootstrap confidence intervals were calculated utilizing the stochastic sampling technique. For each population in Tables 4-5, these estimates are made with various fractions of the initial data set: I00, 50, 35, or 20%.
Examination of the size of the confidence intervals in Tables 4-5 shows that often there are only small increases in interval size with decreasing fractions of the data used. This is especially so when moving from 100 to 50% or 50 to 35%. Although we do not know the true mean longevity of these populations, we can use the observations to estimate the mean and a confidence interval. The confidence interval on these observed longevities almost always overlaps the confidence interval of the predicted longevities.
We have also enumerated the number of deaths that exceed t0.95 (Table 6). Because we have different estimates of A and a for each fraction of the data used, we also get different values of to.95 in each case. Accordingly, we have counted the number of deaths that exceed each of the different t0.95. Because the sample sizes were small, we pooled all counts from like sex and population. In 15 out of 20 comparisons the observed number of deaths in this tail of the distribution is within the expected 95% confidence interval (these are exact confidence intervals computed from the binomial distribution). To evaluate these results, it is important to remember that the value of to.95 is, itself, an estimate. Had we looked at the performance of the Gompertz using the range provided by the 95% confidence intervals (Tables 4-5) on t0.95, the level of agreement would have been much better.
These results suggest that considerable savings in time and energy may be achieved by computing the expected size of these confidence intervals ahead of time for various combinations of N and fraction of data used. This can be done with the techniques outlined here if one can make a reasonable guess at A or or. If this cannot be done, a range of A and e~ values might be examined. Estimates of A and eL can be obtained if one has estimates of two values of t~. If the values of 13 used are not very close (e.g., 50 and 90%), then estimates may be obtained by numerically solving the systems of nonlinear equations that result from eq. (11), using techniques like the method of bisection or Newton's method (Philips and Taylor, 1973). This technique is also useful for estimating the starting values of A and et when doing maximum likelihood.

DISCUSSION
Our first conclusion is that the linear regression methods typically used to estimate parameters of the Gompertz do not work well and should probably never be used. The estimates obtained by linear regression depend on unstable features of the data set, which make them useless for comparison across populations. These problems can be easily avoided by using the maximum likelihood techniques described here. The max- imum likelihood method, however, assumes sampling is binomial. This is a reasonable assumption for the laboratory populations studied here, but may not be justified in other circumstances. Clearly, appropriateness of this assumption must be evaluated on a case-by-case basis. For populations in which the Gompertz equation appears to work well, estimates of mean longevity, t0.95, t0.99, etc., can be made by observing the first 20% or so of all deaths. Confidence intervals on these estimates may be made as small as desired by simply increasing the initial cohort size. The design of experiments can be aided by first  I  I  I  I  I  1   I  I  I  I  I  I   x~,l" 02   ~  I  I  I  I  I   0 3   I  I  I  I  I  I   I  I  I  I  I  I   I  I  T  making a guess at A and (x for the population of interest. Following the bootstrap process outlined here, estimates of the size of confidence intervals for any of the important statistics that follow from the Gompertz equation can be made. The effects of initial cohort size and fraction of the data examined can be easily evaluated with these techniques to ascertain the most efficient manner to collect mortality data. For populations that do not obey Gompertz mortality kinetics, the techniques in this article can be easily modified. For instance, Carey et al. (1993) suggest that medflies violate crucial assumption of the Gompertz equation. Hirsh et al. (1944) describe sev- eral two parameter models that appear to provide much better descriptions of medfly mortality kinetics. All the techniques described here can be used easily with these other models. In fact, one simply has to replace two lines of computer code in all the programs used in this article to compare equivalent results for the JMA and gamma models described by Hirsh et al. (these programs are available from the authors upon request). Not all populations will be expected to obey the Gompertz equation. Nor do we expect the Gompertz to work well in the very old age classes, in which natural selection has not been an effective force at molding rates of mortality (Nusbaum, Mueller, and Rose, submitted). Nevertheless, for populations that are relatively long-lived and follow Gompertz mortality curves, or some other simple mortality kinetic model, sub- Although all populations have been pooled the particular to.95 from each population was used in these calculations.
stantial savings in estimating longevity could be achieved by applying the methods presented here.