AMERICAN JOURNAL OF PHYSICAL ANTHROPOLOGY 117:297–309 (2002) Deconstructing Death in Paleodemography Lyle W. Konigsberg* and Susan R. Frankenberg Department of Anthropology, University of Tennessee, Knoxville, Tennessee 37996-0720 KEY WORDS life tables; hazards analysis; human skeletal biology ABSTRACT In 1992 in this Journal (Konigsberg and Frankenberg  Am. J. Phys. Anthropol. 89:235–256), we wrote about the use of maximum likelihood methods for the “estimation of age structure in anthropological demography.” More specifically, we presented a particular method (the “iterated age-length key”) from the fisheries literature and suggested that the method could be used in human and primate demography and paleodemography as well. In our paper (section titled “Some Future Directions”), we spelled out two broad areas that we expected to see develop over the ensuing years. First, we felt that the use of explicit likelihood methods would open up interest in basic estimation issues, such as the calculation of stan- dard errors for demographic estimates and the formulation of tests for whether samples differed in their demographic structure. Second, we felt that the time was ripe for hazards analyses that would incorporate the uncertainty in estimation that follows from using age “indicators” rather than known ages. While some of these developments have occurred during the last decade, few have been reported in the American Journal of Physical Anthropology. In this paper we resolve some issues from our 1992 paper, and attempt to redress this deficit in the literature by reviewing some recent developments in paleodemography over the past decade. Am J Phys Anthropol 117:297–309, 2002. © 2002 Wiley-Liss, Inc. In a 1992 article “Estimation of Age Structure in Anthropological Demography,” published in this Journal, we noted, “The stumbling-blocks faced in many anthropological demography studies, and virtually all paleodemographic studies are that: 1) population growth rate may not be known, 2) samples many not be representative, and 3) ages are usually estimated rather than known” (Konigsberg and Frankenberg, 1992, p. 235). Now, a decade later, it seems reasonable to reevaluate these problems and see what progress has been made in the last 10 years. Our primary focus in this paper will be on the third problem of age estimation, though we also touch tangentially on the second issue of how to assess statistically whether “mortality” samples are biased by nonbiological factors (for a brief discussion of such “natural and cultural filters,” see Konigsberg and Frankenberg, 1994). We will not deal with the first problem of population growth, as others have recently proposed promising solutions for the simultaneous estimation of growth rate and mortality (described in Wood et al., 2002). Our 1992 discussion of developments in the estimation of age structure in paleodemography was roundly criticized by Bocquet-Appel and Masset (1996) in an article entitled “Paleodemography: Expectancy and False Hope.” At the time this critique appeared, we felt that the methods that we and others were working on were sufficiently different and advanced to overcome their objections. These new developments in paleodemographic methods, plus many others that we did not anticipate, were recently published in a Cambridge University Press volume coedited by Hoppa and Vaupel (2002). We focus here on an erroneous application in our 1992 article that Bocquet-Appel and Masset (1996) did not criticize, and that has been repeated by BocquetAppel and Bacro (1997) and by Jackes (2000). Specifically, we are culpable for introducing in this Journal some bad thinking in model estimation by providing in our 1992 paper (our Table 2) an example application with negative 5 degrees of freedom. More recently, Bocquet-Appel and Bacro (1997, their Table 2) framed most of their analysis around model estimation with negative one degree of freedom, and Jackes (2000, her Fig. 15.7) presented an analysis with negative 11 degrees of freedom. In part then, this paper makes amends for some misapplications from the past. This paper also argues that it is high time to move past life table approaches in paleodemography. We first briefly describe “contingency table” approaches to paleodemography, reanalyzing the data presented in Bocquet-Appel and Bacro (1997) to motivate the discussion. These are essentially nonparametric methods in which one or more “age indica- © 2002 WILEY-LISS, INC. Grant sponsor: National Science Foundation; Grant number: SBR9727386. *Correspondence to: Lyle W. Konigsberg, Department of Anthropology, University of Tennessee, 250 South Stadium Drive, Knoxville, TN 37996-0720. E-mail: firstname.lastname@example.org Received 2 February 2001; accepted 11 October 2001. DOI 10.1002/ajpa.10039 Published online in Wiley InterScience (www.interscience.wiley. com). 298 L.W. KONIGSBERG AND S.R. FRANKENBERG tors” are cross-tabulated against age categories in a reference (known age-at-death) collection in order to estimate the proportions of deaths in age categories for a target (unknown age-at-death) sample. Contingency table models are limited by the fact that there must be at least as many indicator states as the number of age categories in order for the model to be identified (i.e., estimable). Next, we will consider hazard models as an alternative to contingency table analyses. The hazard model is used with reference sample information in order to predict the proportion of individuals that should be in each indicator state in the target sample, and the hazard parameters are then adjusted until there is the best possible fit between the expected proportions in indicator states and the observed proportions. Because hazard models have a small number of parameters, usually considerably less than the number of age classes in the contingency table approach, hazard models often can be estimated where contingency table methods would fail. Within our discussion of hazard models, we will consider a number of specific issues, including: 1) estimation of confidence limits for hazard parameters when ages are estimated rather than known, 2) calculation of mean age-at-death and its standard error from hazard models, again incorporating the uncertainty in ages, and 3) calculation of entropy and its standard error for a target age-at-death distribution based on a hazard model. Then, again using hazard models, we will discuss the problem of how to identify whether a target sample is a reasonable representation of deaths that occurred in a prehistoric population. Finally, we will discuss the overall problem of how age estimation should be done in paleodemography. In addition to reviewing some previous methods, we note the utility of fully parametric models where the target sample mortality is fit using a hazard model and the development of the “age indicator” is also fit with some form of parametric model. CONTINGENCY TABLE PALEODEMOGRAPHY In Konigsberg and Frankenberg (1992), we suggested using a technique originally developed in fisheries management (Kimura and Chikuni, 1987) referred to as the “iterated age length key” (IALK). IALK was shown by Kimura and Chikuni (1987) to be an expectation-maximization algorithm (Dempster et al., 1977) that leads to maximum likelihood estimates for the age-at-death structure for a target sample. Indeed, all of the methods we will consider in this paper are based on the method of maximum likelihood. IALK uses a contingency table of an age indicator against known age classes in a reference sample in order to infer the age-at-death structure in a target sample. Bocquet-Appel and Masset (1996) suggested that they had historical priority on the development of the method, which they had developed out of an “iterative proportional fitting procedure” (IPFP) commonly used to fit log-linear mod- els. Bocquet-Appel and Masset (1996) were trying to fit an unobserved marginal distribution (age-atdeath), using the same type of information that we were using in 1992. The relationship between IPFP and IALK was unclear in Bocquet-Appel and Masset (1996), because they obtained slightly different results using the two techniques. A year later, Bocquet-Appel and Bacro (1997, p. 571) wrote in this Journal, “One can demonstrate analytically that both algorithms, although written differently, have the same step two and lead in fact to the same solution.” Yet, when they applied the two techniques to the same data, they also obtained slightly different answers. We programmed both procedures (see http://konig.la. utk.edu/paleod.html) and, provided that they are started from the same starting values for the estimated age-at-death distribution, the algorithms take identical steps through the parameter space. Though they are written differently, IPFP and the IALK are the same, and following historical priority, we will refer to the method as IPFP. IPFP does differ, however, from an algorithm developed by Hoenig and Heisey (1987) that accounts for finite reference sample size. Various authors (Konigsberg and Frankenberg, 1992; Bocquet-Appel and Bacro, 1997; Jackes, 2000) could have avoided erroneous application of IPFP resulting in overfitting (i.e., negative degrees of freedom) by reading more closely the following from Clark (1981, p. 299): If I ⬍ J (i.e., the number of length intervals is less than the number of age-groups with distinct length distributions), there will usually be a multiplicity of algebraic solutions and therefore no useful estimate. The recent analyses by Bocquet-Appel and Bacro (1997) and Jackes (2000), which estimate the age structure for a target sample as a single event, should have detected this transgression in their evaluation statistics. Using the goodness-of-fit test given in Kimura and Chikuni (1987, p. 28) and in Hoenig and Heisey (1987, p. 235), they would have found that the test gives negative degrees of freedom if there are more age classes than indicator stages. Bocquet-Appel and Bacro (1997) and Jackes (2000) also should have recovered excessively large calculated standard errors when they fit more age classes than indicator states, or if they had looked at the likelihood surface, they would have noted remarkably flat surfaces, a clear sign that the model is overparameterized. (Goodness-of-fit and standard errors were not calculated in the example in Konigsberg and Frankenberg (1992), because it was a simulation study.) In her recent analysis, Jackes (2000) attempted to estimate a model with 17 age classes (i.e., a model with 16 parameters, as the sum of the proportions in age classes must equal 1.0), using the SucheyBrooks six-phase system. Her results yielded a number of age classes with zero frequencies, a clear sign DECONSTRUCTING DEATH IN PALEODEMOGRAPHY 299 TABLE 1. Distribution of classifiable adults by age and stage of trabecular involution of the femoral head for the Coimbra reference sample and Loisy-en-Brie target sample1 Stage Age I II III IV V VI 23–29 30–39 40–49 50–59 60–69 70–79 80–89 Loisy-en-Brie 8 2 0 0 0 0 0 2 19 18 6 2 0 0 0 8 30 43 29 26 9 7 2 31.5 7 25 27 37 28 28 10 40.5 1 1 5 13 9 10 10 12 0 0 1 0 1 3 5 2 1 From Bocquet-Appel and Bacro (1997: their Table 2 and p. 573). We use the values from within their Table 2, which sum to N ⫽ 422, and not the N ⫽ 421 given in their footnote and text. that the model was not identified. Bocquet-Appel and Bacro (1997) presented a less extreme estimation problem in trying to fit a seven age-class model using six stages. In the remainder of this paper we will deal with their application, discuss some related issues from Bocquet-Appel and Masset (1996), and sketch some of the directions in which we see paleodemography now heading. Bocquet-Appel and Bacro (1997) presented tables of age and stage data for their reference sample and target sample that we reproduce here in Table 1. Their analyses for adults were based on stage data for femoral head trabecular involution, using 422 known-age individuals from an anatomical collection (Coimbra) and 96 individuals from a “target” archaeological sample (Loisy-en-Brie). Rather than use IPFP to estimate the age structure of the target sample here, we explicitly maximized the log-likelihood using the Nelder-Mead simplex algorithm (Nelder and Mead, 1965) available in the “optim” function of the statistical/graphics package “R.” “R” is a major package distributed under the GNU Public License, and as a consequence it can be freely downloaded (see cran.r-project.org for binaries and source code, Ripley (2001) for an overview, and Cribari-Neto and Zarkos (1999) for a recent review). Using “optim” to maximize the log-likelihood leads to estimated proportions in each of the seven age classes of 0.1645, 0.0269, 0.2781, 0.2061, 0.1197, 0.1935, and 0.0113. “Optim” can also be used to find the standard errors of these estimates by: 1) requesting calculation of the Hessian matrix (i.e., the matrix of second partial derivatives of the log-likelihood with respect to the parameters), 2) changing the sign of all elements of the matrix, 3) inverting the matrix, and then 4) taking the square root of diagonal elements. Rounding to the nearest whole number the standard errors on these proportions are 576, 3,156, 9,496, 14,355, 25,786, 22,084, and 3,737, respectively. The analytical standard errors for the estimated proportions are even larger. Clearly, these are not usable results. As indicated above, the reason things went awry here is that the model, with negative 1 degree Fig. 1. Normed likelihood for proportion of deaths in the 40 –50-year age interval (curved line). Horizontal line gives asymptotic 95% confidence set at the line’s intersections with the normed likelihood. of freedom, is overparameterized. Another way to see this is to plot the normed likelihood (Lindsey, 1996, p. 82) for one of the parameters, while fixing all other parameters save one at their maximum likelihood estimates. In Figure 1, we plotted the normed likelihood for the proportion of deaths in the 40 – 49-year age interval, allowing only the proportion of deaths in the 50 –59-year age interval to vary so that the age distribution sums to 1. The asymptotic 95% confidence set for this parameter consists of all values of the death proportion with a normed likelihood greater than or equal to 0.1465 (see equation 7.20 in Shao, 1999, p. 447), which we have plotted as a horizontal line. From Figure 1 it is clear that the normed likelihood never reaches this horizontal line within the parameter space, so that the confidence set is undefined. Further, because Figure 1 is drawn as a profile likelihood, the normed likelihood appears even less flat than it actually should be. As Lindsey (1996, p. 111) notes, “A profile likelihood gives a narrower range of likely values of the parameter of interest than is necessarily appropriate” because it does not allow for uncertainty in the other parameters (death proportions in the case used here). One possible solution to this problem of overparameterization is to decrease the number of age intervals by combining some adjacent age categories. For example, if we combine the last two age categories (70 –79 and 80 – 89 years), then there are six classes for the indicator state, and six age classes. This leads to six equations with six unknowns, and can be solved using any of a number of linear algebraic approaches (e.g., QR decomposi- 300 L.W. KONIGSBERG AND S.R. FRANKENBERG tion). In the more complicated situation where there are more equations than unknowns, as discussed by Clark (1981), some age classes could include negative estimated counts when fit by solving the system of linear equations. When there are more unknowns than equations, as in the original example in Bocquet-Appel and Bacro (1997), there is no unique solution. For the six age-class case (with the 70 –79-year and 80 – 89-year age classes combined), we solved the linear equations to estimate the proportions in age classes as 0.1678, 0.0079, 0.3360, 0.1182, 0.2761, and 0.0939. Unfortunately, the associated standard errors that we obtained from the Hessian matrix both analytically and numerically are still unacceptably large (0.2305, 0.7879, 1.0613, 0.7593, 0.5534, and 0.2237, respectively). In addition to the fact that the 95% confidence intervals in this case include negative values for the proportion of deaths in every age class, the estimated age distribution is very bumpy in nature. For example, about 17% of deaths are estimated as occurring between 23–29 years, about 34% between 40 – 49 years, but less than 1% of deaths occur between these two age classes at 30 –39 years. This “bumpiness” occurs because in the contingency table approach to paleodemography, there is no parametric model describing mortality. Next, we address this problem using hazards analysis. HAZARD MODELS IN PALEODEMOGRAPHY A solution to the estimation problems we noted above could have been developed from our 1992 article. There we suggested (Konigsberg and Frankenberg, 1992, p. 252), “A future direction that we expect to see in anthropological demography and paleodemography is the incorporation of uncertainty of age estimates into reduced parameterizations of life table functions.” This was not a particularly novel concept, as Gage (1988) had already shown the utility of hazard models in paleodemography. We further suggested that one could maximize the loglikelihood that accounted for age uncertainty in a hazard model, an application inaccessible to Gage (1988) working from life table tabulations that lacked the osteological information used to estimate ages. In the following, we show how hazard models can be fit to the Loisy-en-Brie data while accounting for the fact that ages are estimated rather than known. The first step in this process is counterintuitive, in that we will fit a model that perfectly recovers the observed proportions of Loisy-en-Brie adult skeletons in the six femoral head stages, rather than fitting a hazard model. In statistical parlance, such a model is usually called a “saturated” model. This model is easy to find, as we simply use the observed proportions in the femoral stages for Loisy-en-Brie, which are 0.02083, 0.08333, 0.32813, 0.42188, 0.12500, and 0.02083. The model has five parameters, because given the constraint that the propor- tions must sum to 1.0, if we know any five parameters then we also know the sixth. The log-likelihood for this saturated model is defined (up to an additive constant) as the sum of the products of observed counts with their natural log proportions. Consequently, we have for the log-likelihood of the saturated model: 2.0 ⫻ ln(0.02083) ⫹ 8.0 ⫻ ln(0.08333) ⫹ . . . ⫹ 2.0 ⫻ ln(0.02083), which equals ⫺130.3731. Lindsey (1996, p. 77–78) refers to this quantity as the “empirical” log-likelihood. Now we can fit a hazard model, which for simplicity’s sake will be a single-parameter exponential model. We assume a survivorship of 1.0 at age 23 (i.e., life begins at age 23, so we subtract 23 years from all ages). The exponential model gives survivorship to exact age t as: S共t兲 ⫽ exp共⫺a 2 ⫻ t兲, (1) where a2 is the single parameter to be estimated. (We use the subscript to differentiate this parameter from later ones to be considered, following the numbering system in Gage, 1988.) If we assume a starting value for a2, say 0.1, then we can fit survivorship at each of the age divisions given in the life table of Bocquet-Appel and Bacro (1997), and then difference the series to get the expected proportion of deaths in each age interval, which we place in a 7 ⫻ 1 column vector d. We then define a row stochastic matrix P ila from Table 1 that gives the probabilities from the reference sample of being in the ith stage of the indicator conditional on being in the ath age interval. If ni is a column vector of counts of the number of skeletons from the target sample in each indicator state, then the log-likelihood is (up to an additive constant) ln(d⬘Pila)ni, which with a2 equal to 0.1 gives a log-likelihood of ⫺156.1552. If we numerically maximize this model (again using “optim” to find the value of a2 at which the log-likelihood is highest), then the maximum likelihood estimate for a2 is 0.02924. The log-likelihood at this point is ⫺132.1952. Negative two times the difference between the log-likelihood for this model and the log-likelihood from the saturated model is a number called the “deviance” or a “likelihood ratio” (see Agresti, 1996, p. 96, or Kennedy, 1992, p. 59 – 62). The deviance is asymptotically distributed as a chi-square random variable, with degrees of freedom equal to the difference in the number of parameters between the hazard model and the saturated model. In this case, the deviance is equal to 3.6442 with 4 degrees of freedom, which gives an associated probability value of 0.4563. This relatively high probability value suggests that little is lost in “moving” from the saturated model with five parameters “down” to the more parsimonious exponential model with only one parameter. The one-parameter exponential model is generally viewed as too simplistic to capture adult mortality, even though it appears to do an adequate job of recovering the observed proportions of individuals in indicator states for Loisy-en-Brie. We can fit a DECONSTRUCTING DEATH IN PALEODEMOGRAPHY 301 Fig. 2. Normed likelihood surface for the two hazard parameters found by using a Gompertz model of survivorship and the likelihood of target skeletons in particular indicator states being within particular age intervals. slightly less parsimonious model, which is the twoparameter Gompertz model. The Gompertz model gives survivorship to exact age t as: S共t兲 ⫽ exp 冉 冊 a3 共1 ⫺ exp共b 3 ⫻ t兲兲 , b3 (2) again following the notation of Gage (1988), where a and b are subscripted to indicate that they represent the third component of mortality in a Siler model. Numerically maximizing the log-likelihood in the same way as for the exponential model, we find the hazard parameters equal to 0.0129 and 0.0416, with standard errors of 0.0065 and 0.0221, respectively. Figure 2 contains a plot of the normed likelihood surface for the two parameters, which shows that the surface is quite convex. Again, the confidence set in any one direction (i.e., for one of the parameters, holding the other constant) is given where the normed likelihood falls above 0.1465. The Gompertz model gives a log-likelihood value of ⫺130.4562, yielding a deviance of 0.1662 with 3 degrees of freedom, and an associated probability value of 0.9829. This very high probability indicates that the Gompertz model does an excellent job of recovering the observed proportions in indicator stages for Loisyen-Brie. We also estimated a three-parameter Makeham model, but the fit is trivially different from the Gompertz model. The deviance for the Makeham model is 0.0414 with 2 degrees of freedom, resulting in a probability value of 0.9795 for the goodnessof-fit test. Because the Gompertz model is “nested” within the Makeham, we can subtract the deviances and degrees of freedom for one model from the other to calculate an improvement 2. The resulting value of 0.1248 with 1 degree of freedom gives a probabil- 302 L.W. KONIGSBERG AND S.R. FRANKENBERG Fig. 3. Kernel density plots of the two hazard parameters for asymptotic normal distributions from the original analysis (solid line), bootstrapped target (open circles), and bootstrapped target and reference (filled circles) samples for 500 samplings. ity value of 0.7239, indicating that there is little to be gained by adding the additional parameter in the Makeham model. Simulation and confidence limits Bocquet-Appel and Masset (1996, p. 574) raised the question, “Besides the expectancy of an estimate, what is its variance of error?” They attempted to answer this question using simulation studies. We present a more modest simulation study here, which takes the estimated Gompertz model for the 96 adults from Loisy-en-Brie and simulates in order to establish the statistical properties of the estimates. This is largely an unnecessary exercise, as we have presented a maximum likelihood estimator whose statistical properties are fairly well understood. In order to simulate counts of the number of individuals in each stage for the femoral head, we drew 96 individuals from the multinomial age-atdeath vector d and then applied the appropriate row from Pila to draw on the multinomial distribution of indicator conditional on age. We repeated the simulation 500 times, and estimated the Gompertz parameters just as we did for the actual observations. This simulation consequently is a parametric bootstrap in the sense that we take the estimated Gompertz parameters from Loisy-en-Brie to use as the generating values in the simulation runs. Figure 3 shows kernel density plots for the a3 and b3 parameters from the 500 simulations, as well as the asymptotic normal distributions from our original analysis. The normal distributions are drawn using the estimated values as the means, and the standard error of estimates from “optim.” While the bootstraps show a bit more dispersion as well as some asymmetry relative to the asymptotic results, the results do not suggest that there is “no hope with the expectation” (Bocquet-Appel and Masset, 1996, p. 578). As Bocquet-Appel and Masset (1996) also were concerned about the effects of a finite reference sample, we bootstrapped the reference sample by drawing from the two-dimensional contingency table in Table 1, and by following the previously described simulation procedure. As seen in Figure 3, bootstrapping the reference has virtually no effect on the simulation results. Bocquet-Appel and Masset (1996, p. 576) argued, “Bootstrapping the reference sample in order to generate an expectation for the estimates would be a wrong approach, because the expectation would be in the expectancy domain of the reference sample.” If by this statement they meant that such a bootstrap inappropriately carries with it the “uniformitarian assumption” of equivalent aging rates in the reference and target sample (Howell, 1976), then this is an insurmountable stumbling block. It is not possible to test the assumption of uniform rates of aging in a univariate, or single-age indicator, setting. When there are multiple age indicators, in contrast, it is possible to test for “inconsistency” in information (Brown and Sundberg, 1987) in the target sample relative to the reference sample. Although one of us applied such an analysis to stature estimation (Konigsberg et al., 1998), we are unaware of similar applications made for age estimation in paleodemography. A statistical summary from hazard models: mean age-at-death Bocquet-Appel and Masset (1996, p. 578) felt on the basis of their simulation work that “it is not possible to estimate reasonably the shape of an age distribution,” but that “one can nevertheless obtain a fairly good estimation for the direction of the overall trend of an undefined distribution: its average.” They consequently suggested calculating the mean age-at-death for adults. While it is debatable whether mean age-at-death and other summary measures given in Bocquet-Appel and Masset (1996) and Bocquet-Appel and Bacro (1997) are all that can be calculated in paleodemographic studies, hazard models do not preclude us from extracting such summary measures. In fact, any desired summary measures can be calculated from the hazard model parameters, although this is less immediate and direct than working with the parameters themselves. Assuming zero population growth, life expectancy can be found from a hazard model by integrating the survivorship, which we do in “R,” using the add-on library “integrate.” For Loisy-en-Brie, the life expectancy at age 23 is estimated at 28.92, for a mean age-at-death for those who die after age 23 of 51.92 years, which is close to estimates in Bocquet-Appel and Bacro (1997) of 52.36 and 52.27. The standard DECONSTRUCTING DEATH IN PALEODEMOGRAPHY error of the mean that we calculated can be found using the delta method (Shao, 1999, p. 45), for which we will need numerical derivatives of the integral with respect to the two parameters in the Gompertz. We find these derivatives using “numericDeriv” in the “R” add-on library “nls.” The asymptotic standard error is 3.18 years, similar to the figure of ⫾3 years that Bocquet-Appel and Bacro (1997) reported from simulation work (see also Bocquet-Appel, 1994). The mean age-at-death we estimated from the IPFP is 52.40 years, with a standard error also from the delta method of 119,991 years, which is much greater than the 3 years given in BocquetAppel and Bacro (1997). Death and entropy in ancient France Now we are faced with a dilemma. We have an estimate of the mean age-at-death and its standard error for Loisy-en-Brie, but we have no mortality model. We could, for example, assume that everyone at Loisy-en-Brie who survived past age 23 years had the temerity to die at precisely 23.0 ⫹ 28.92 ⫽ 51.92 years of age, which would produce a rectangular survival curve. Or we could assume that starting at age 23 years, the instantaneous hazard is 1/28.92 for the rest of life, which would produce a very different negative exponential adult survival curve. By arguing that the shape of mortality cannot be estimated, Bocquet-Appel and Masset (1996) deprive us of not only the distribution, but also the dispersion of deaths. Some idea of the dispersion can be recovered by estimating the standard deviation for age-atdeath for those who died past age 23 years. From the Gompertz model, we again need to integrate in order to find the standard deviation for age-at-death. The variance for age-at-death for those who die above age 23 is: 2 冕 2 t ⫻ S共t ⫺ 23兲 dt ⫺ e 23 , (3) 23 2 where e23 is the square of life expectancy at age 23. For Loisy-en-Brie, the estimated standard deviation for age-at-death in those individuals surviving to age 23 years is 16.41 years, with an asymptotic standard error of 7.85 years. But the standard deviation of age-at-death is not a particularly useful demographic parameter, and it is difficult to compare across populations or samples. In place of the standard deviation, we suggest using the entropy (Demetrius, 1978, 1979; Goldman and Lord, 1986; Vaupel, 1986) of deaths at and over age 23, which is defined as: ⫺ H⫽ 冕 S共t ⫺ 23兲 ⫻ ln共S共t ⫺ 23兲兲 dt 23 e 23 . (4) Entropy can range from a minimum of 0.0, in which case all deaths occur at a single age, to a maximum 303 of 1.0, in which the hazard rate is a constant across age. For Loisy-en-Brie the estimated entropy for those who die past age 23 is 0.5218, with a standard error of 0.1340. The 95% confidence interval for entropy easily includes most of the entropy values for age-at-death for deaths past age 20 in the model life tables of Coale and Demeny (1966). ARE PALEODEMOGRAPHIC SAMPLES REPRESENTATIVE? One of the most niggling problems in paleodemography is the question of whether archaeological skeletal samples are truly representative of deaths that occurred in a prehistoric population. We touched very briefly on this issue before (Konigsberg and Frankenberg, 1994), and Hoppa (1996) wrote his dissertation on the topic. Bocquet-Appel and Masset (1996, p. 574) ran their simulations of mortality patterns “representing various models of age distributions which could possibly be found in cemeteries,” though their “bimodal” age distribution is not well-fit by any reasonable model of human adult mortality. If the point of their simulations was to show that paleodemographic reconstruction falls apart when samples are nonrepresentative, demographically unusual, or simply out of the range of expectation, then it is difficult to agree with their intent. There is a misperception in the field that the analyst has no recourse (and no interest) when an age structure estimated for a target sample is driven by cultural and archaeological sampling. To the contrary, even though we will get back a Gompertz (or any other) model when we fit a Gompertz (or any other) hazard model, this does not mean that we lack internal checks on reasonableness. We have already seen that the Gompertz model for Loisy-en-Brie yields a deviance value of 0.1662, which with 3 degrees of freedom gives a probability of about 0.9828. Clearly, the Gompertz does a quite reasonable job of reproducing the observed distribution of the indicator state. As an example of the failure of the Gompertz to reproduce an observed distribution, we can take a made-up vector with counts in the indicator states of 10, 20, 50, 100, 50, and 10. Using the reference information in BocquetAppel and Bacro (1997), we would estimate a Gompertz model with a3 ⫽ 0.01855 and b3 ⫽ 0.0019. While these parameters look reasonable, they give a likelihood ratio chi-square of 13.86, which with 3 degrees of freedom produces a probability of about 0.003. The Gompertz consequently is not doing a good job of reproducing the distribution of indicator states in this fictitious example. Loisy-en-Brie also contained the remains of individuals who died at ages less than 23 years, as well as 12 individuals who died at ages greater than 23 years but for whom the femoral head was unobservable. Bocquet-Appel and Bacro (1997) chose not to include these individuals in their analysis, presumably because of the problem of possible underenumeration of subadults. This approach is logically 304 L.W. KONIGSBERG AND S.R. FRANKENBERG TABLE 2. Loisy-en-Brie age distribution for those not aged by the femoral head1 1 Age Count 0 1–4 5–9 10–14 15–19 20–22 23–89 4 15 15 7 9 6 12 From Bocquet-Appel and Bacro (1997: their Table 1 and p. 573). Fig. 5. Hazard rate for Loisy-en-Brie when four 0 –1-year-olds are added to the sample, as compared with the hazard rate from model West 3 of Coale and Demeny (1966) for females. Fig. 4. Hazard rate for Loisy-en-Brie, with a 95% confidence interval (solid line) as compared with hazard rate from model West 3 of Coale and Demeny (1966) for females. inconsistent, because it excludes the subadults out of hand, without ever considering whether or not these individuals are properly represented in the assemblage. We can incorporate the subadults (listed in Table 2) and the 12 adults in whom the age indicator was unavailable, treating these individuals as having interval-censored ages-at-death. In the case of the 12 adults, we treat their ages as being equal to or greater than 23 years. With this additional information, we can fit a fiveparameter Siler model (Gage, 1988; Wood et al., 2002). The survivorship from this model is: 冉 S共t兲 ⫽ exp ⫺ a1 共1 ⫺ exp共⫺b 1 䡠 t兲兲 b1 ⫻ exp共⫺a 2 䡠 t兲exp 冉 冊 冊 a3 共1 ⫺ exp共b 3 䡠 t兲兲 . b3 (5) The estimated parameters are 0.0220, 0.1144, 0.0054, 0.0026, and 0.0507, with standard errors of 0.0110, 0.1491, 0.0108, 0.0011, and 0.0098, respectively. Figure 4 shows the 95% confidence interval around the hazard rate for Loisy-en-Brie, as well as the hazard rate from model Coale and Demeny’s (1966) West 3 for females. The only departure from the model life table appears to be the lower hazard rate for infants at Loisy-en-Brie. In Figure 5, we show the hazard rate redrawn by adding four 0 –1year-olds to the Loisy-en-Brie sample, again compared to model West 3. The hazard rates are remarkably similar after the addition of these four infants. HOW SHOULD AGE ESTIMATION BE DONE? There have been a number of suggestions in this Journal as to how age estimation should proceed in paleodemography (Bedford et al., 1993; Jackes, 1985; Lovejoy et al., 1985; Meindl et al., 1985, 1990). Any discussion of paleodemographic methods, whether life table- or hazard analysis-based, cannot ignore the central issue of determination of age-atdeath. Both the method proposed by Jackes (1985) and that developed by Meindl et al. (1985) in this Journal suffer from the problem that the reference sample age structure will influence age-at-death estimation. Bocquet-Appel and Masset (1982) first discussed the magnitude of this problem in paleodemography. We illustrate this problem for the methods of Jackes (1985) and Meindl et al. (1985, 1990), and suggest ways to conduct age-at-death estimation that both avoid such “mimicry” (Mensforth, 1990) and allow incorporation of uncertainty of age estimates. Jackes (1985) suggested that the mean and standard deviation for age-at-death within pubic symphyseal stages could be used to probabilistically assign ages-at-death in a target sample. To implement this, she assigned individual probabilities to integer ages, using the cumulative normal truncated to the middle 95% of the distribution. This method is dif- DECONSTRUCTING DEATH IN PALEODEMOGRAPHY 305 ing information on the reference sample age distribution. We can see this by writing the mean age and variance of age within indicator state i as: a 兩i ⫽ 冕 冕 a Pr共i兩a兲f共a兲 da 0 V a 兩i ⫽ (6) 共a ⫺ 共a 兩i兲兲 2 Pr共i兩a兲f共a兲 da, 0 Fig. 6. Log-normal distributions of age for stages 0 –2 (“stage 0”) and 3–5 (“stage 1”) of third symphyseal component in the Korean War dead of McKern and Stewart (1957) and Terry Anatomical Collection samples. ficult to justify, as it is highly unlikely that age is normally distributed within stages (discussed further below). We explore a slightly different approach here, which is to use the full normal probability density function for age within indicator states. As an abbreviated example, we use a scoring of the third component of McKern and Stewart (1957) of the symphyseal surface, where we have simply dichotomized between stages 0 –2 and stages 3–5. Essentially, this divides pubic symphyses into those with incomplete ventral symphyseal rims and those with complete rims. We have scores for this trait on casts of 351 male pubic bones from the Korean War dead study of McKern and Stewart (1957), and 413 male pubic bones from the Terry Anatomical Collection. Rather than using the distribution of raw age within these two stages, we use the natural log of age. For the Korean War dead, the mean and standard deviation of log age are 3.07 and 0.156 within stage 0 (original 0 –2 of McKern and Stewart, 1957), and 3.44 and 0.173 for stage 1 (original 3–5 of McKern and Stewart, 1957). The comparable figures from the Terry Collection are 3.47 and 0.361 for stage 0 and 3.90 and 0.324 for stage 1. Figure 6 shows a plot of these log-normal distributions. Now we are presented with quite a dilemma. If we plan to use this aging information for a target sample, should we use the log normal from the Korean War dead, or from the Terry Collection? More importantly, why are these distributions so different? The reason that the information from these two collections is so different is that conditioning on indicator, instead of on age, has the effect of includ- where Pr(i兩a) is the probability of being in stage i conditional on being exact age a, and f(a) is the probability density function for age-at-death. To model Pr(i兩a) we use a probit regression of the indicator on log age in the Korean War dead and in the Terry Collection. For the Korean War dead, the intercept is 18.88 and the slope is 5.61, whereas for the Terry Collection, the intercept is 6.65 and the slope is 1.99. Viewed as a transition analysis (Boldsen et al., 2002), the mean age-to-transitions for the Korean War dead and Terry Collection samples are exp(18.88/5.61) and exp(6.65/1.99), or about 28.9 and 28.3 years, respectively. The only essential difference between these two collections in the distribution of age-to-transition is in the dispersion, with the log scale standard deviation being 1/5.61 and 1/1.99, or 0.178 and 0.502, for the Korean War dead and Terry Collection samples, respectively. These two samples have markedly different ageat-death structures, and it is this factor that contributes to the different distributions for age-at-death within pubic stages. Fitting a Gompertz model to mortality at and past age 16, we obtain a3 ⫽ 0.082 and b3 ⫽ 0.053 for the Korean War dead sample, and a3 ⫽ 0.011 and b3 ⫽ 0.039 for the Terry Collection sample. Applying Equation 6 with the appropriate Pr(i兩a) for each sample, we can recover the mean and variance of the age-at-death distribution within stage for each collection. On the log scale, these values for the Korean War dead are 3.06 and 0.186 within stage 0, and 3.43 and 0.175 within stage 1. For the Terry Collection these values are 3.41 and 0.376 within stage 0, and 3.91 and 0.341 within stage 1. All of these figures agree well with the direct calculations of means and standard deviations given above. As a consequence of the age determination method of Jackes (1985, p. 294) depending on the reference sample age distribution, we must disagree with her statement that her method would “alleviate the problems of adult age assessment.” In fact, her method only serves to formalize the critique that Bocquet-Appel and Masset (1982) leveled against age estimation in paleodemography. The suggestion (Lovejoy et al., 1985; Meindl et al., 1985, 1990) to seriate target samples on separate age indicators, apply point estimates to each skeleton, and then weight the estimates via a principal components analysis, suffers from the same problem that we have noted for the method of Jackes (1985). We do not, in principle, disagree with seriating tar- 306 L.W. KONIGSBERG AND S.R. FRANKENBERG get sample skeletal material on the basis of morphological indicators, although later this may call into question the usual statistical assumption of individual variate values being independently distributed. However, a major problem arises when ages are assigned across the seriated sample. These ages must be assigned in regard to a reference sample, and as such they again depend in part on the reference sample age distribution. For example, suppose we seriated the bones (or casts) when we scored Todd phases for the Korean War dead and Terry Collection samples. As a consequence, within each phase for each sample we would have a smear with, e.g., phase III represented at one end by pubes that look more phase two-ish and at the other end by pubes that look more four-ish. Now we have the nontrivial problem of deciding how to establish breaks in the seriation so that we know where one integer age ends and a new age begins. Meindl et al. (1985, 1990) gloss over this part of the analysis by referring to seriating “age,” whereas in fact we have seriated the indicator. They also provide no suggestion for mapping the morphological series into discrete integer ages. Should ages be assigned uniformly within indicator class on the basis of the reference information? Further, as we have chosen to seriate, do indicator classes even exist anymore? Finally, how should we treat reversals where older individuals are seriated morphologically into the younger segment? In desperation we could use the mean and standard deviation of age within indicator classes from the reference sample to establish ages in the target sample, assigning integer ages so that they are approximately normally distributed within the target indicator classes. This brings us full circle to the original problem with the method of Jackes (1985), which we summarized in Equation 6. It is becoming increasingly clear that the only logical way to proceed in estimating age in paleodemography is via Bayes’ theorem (Milner et al., 2000; Hoppa and Vaupel, 2002). Specifically, we have: f共a兩i兲 ⫽ 冕 Pr共i兩a兲f共a兲 , (7) Pr共i兩a兲f共a兲 da 0 which includes the unknown age-at-death distribution f(a). If Pr(i兩a) is a proper probability function, then the integral across age for Equation 7 shown in the denominator gives P̂r(i) the estimated probability that an individual will be in the ith class of the indicator. As a consequence, the sum of the products of counts within indicator classes with log(P̂r(i)) is the log-likelihood for the age-at-death distribution. In this article, we assume that the age-at-death distribution can be represented by a hazard model, so that the task now is to maximize the log-likelihood across the hazard parameters. This is precisely what we did in the example from Loisy-en-Brie, and had we wished to estimate individual ages, we could then apply Equation 7. Although we make use of Bayes’ theorem, it is important to understand that the age estimates we obtain are maximum-likelihood estimates, not Bayesian ones. This is true because Equation 7 uses an estimated hazard, rather than a prior distribution. Figure 7 shows the full posterior densities for age-at-death in Loisy-en-Brie conditional on being in each of the six indicator states. The sawtooth appearance of these densities is due to the fact that we used the multinomial distribution (following Bocquet-Appel and Bacro, 1997) for being in a particular indicator state conditional on being in a particular age class. Aside from producing a very ragged-looking distribution of age-at-death, such an assumption does not make very good use of the reference sample data, because we should have exact ages from that sample. A number of authors (Boldsen et al., 2002; Holman et al., 2002; Konigsberg and Herrmann, 2002; Konigsberg and Holman, 1999) are currently using parametric models that relate the indicator to exact age, while others (Aykroyd et al., 1999; Love and Müeller et al., 2002) have used nonparametric models for the same purpose. Equation 7 also shows why the method of Jackes (1985, 2000, shown as “95% CLP” in some of her figures) for age estimation is difficult to justify. We can simplify the equation a bit by removing the normalization (division by the integral) and rewriting Pr(i兩a) as a likelihood L(a兩i) to get: f 共a兩i兲 ⬀ L共a兩i兲f共a兲, (8) where f (a兩i) is the distribution for age within an indicator stage (i.e., the age distribution posterior to our knowledge about the indicator) and f (a) is the prior distribution for age. The method of Jackes (1985) relies on the posterior following a normal distribution, but makes no stipulations about the distributional forms for the likelihood or the prior. We will assume the simplest case, where the likelihood follows a normal distribution. This is a reasonable assumption for all but the first and last stages, provided that we have modeled the age indicator using a cumulative probit model, and that the mean ages-to-transition (see Boldsen et al., 2002) are not too widely separated across stages. In this case, the likelihood for the first stage follows a normal distribution function (a cumulative normal), while the likelihood for the last stage is a normal survivor function (the complement of a cumulative normal). If the likelihood is normal, then the prior must also be normal in order for the posterior to be normally distributed. However, there is no justification for assuming that mortality, which is what the prior represents, is normally distributed. The argument we give here is based on the Bayesian concept of “conjugacy” (Lee, 1989; Gelman et al., 1995), where once we have established the form for the likelihood, we pick a form for the prior such that the posterior has the DECONSTRUCTING DEATH IN PALEODEMOGRAPHY Fig. 7. 307 Full posterior densities for age-at-death in the Loisy-en-Brie sample conditional on indicator state. same distributional form as the prior (hence the prior is referred to as a “conjugate prior”). For a normal likelihood, the conjugate prior is also a normal distribution. It is possible that we could, in theory, combine some as yet unspecified likelihood with some other prior to get a normal posterior distribution. However, our practical experience from examining distributions of known ages within osteological stages is that ages are not normally distributed. There are two cautionary notes that should be made in regard to age-at-death estimation. The first concerns the nature of estimated individual ages-atdeath. Although Equation 7 shows how to provide age estimates for individual skeletons once a hazard model has been estimated for the target, it is difficult to think of examples where one would actually want individual age estimates in paleodemography. Unlike forensic anthropology, paleodemography has no particular currency in the stock and trade of individual age estimation. On the other hand, a researcher might want to look at the relationship between age-at-death and some other variable, or set of variables, from the bioarchaeological context. If this is the goal, then a simple substitution of a point estimate for age-at-death for each skeleton is an affront to the data. Substituting point estimates neglects the fact that ages are estimated rather than known, and gives the analyst false power in tests. Konigsberg and Holman (1999) previously made this 308 L.W. KONIGSBERG AND S.R. FRANKENBERG point for studies of prehistoric somatic growth, while Milner et al. (2000) extensively considered the statistical basis for individual age estimates. The second cautionary note concerns the fact that our reevaluation of the analysis by Bocquet-Appel and Bacro (1997) of Loisy-en-Brie is based on the summary data they provided. Far more could be learned from reevaluating the raw data. For example, models that replace the F matrix (which we have written as Pila) of Bocquet-Appel and Masset (1996) with smooth functions will usually require access to the raw data from the reference sample. It is for this reason that we make available on the Web our raw reference sample data (http://konig.la. utk.edu/paleod.html). CONCLUSIONS Bocquet-Appel (1994, p. 201) closed a chapter he wrote on “Estimating the Average Age for an Unknown Age Distribution in Anthropology” by suggesting, “One can again take an interest in paleodemography.” We agree with the general tenor of this statement, although it would appear from the voluminous response to the original critique by BocquetAppel and Masset (1982) that paleodemography is a subject in which physical anthropologists have never particularly lost interest. We would strengthen the statement of Bocquet-Appel (1994) to, “We should take an interest in paleodemography,” and we take the remainder of this conclusion to “deconstruct” this statement. First, the emphasis in paleodemography should be on “we,” not “me,” “I,” “one,” or “us” (as vs. “them”). There have been too many fractious statements made concerning paleodemography (Bocquet-Appel and Masset, 1982, 1996; Petersen, 1975), and rarely do these comments provide suggestions on how to improve the field. In contrast, the past few years have seen a number of collaborative efforts in paleodemography (Hoppa and Vaupel, 2002; Paine, 1997), and these have had and should continue to have a very positive effect in advancing the field beyond the rote calculation of life tables. Second, we should take an interest in paleodemography because of the inextricable link between demography and genetics in evolution, whether that evolution is in humans, primates, or mammals (e.g., the treatment of the “evolution of the human life cycle” in Bogin and Smith, 2000). Additionally, the relationships between demography and ecology remain to be explored, and at least some of that exploration will likely include human prehistory as a component (e.g., the comparative discussion of primate demography in Gage, 1998). There is a vast amount to be learned about past demography, and we are just beginning to enter a period when many questions stand a reasonable chance of leading to reasonable answers. ACKNOWLEDGMENTS The Korean War dead and Terry Collection sample data were collected under funding from the Na- tional Science Foundation (SBR-9727386). We are especially grateful to Nick Herrmann and Danny Wescott for helping the first author collect the data, and to Dave Hunt at the Smithsonian for his tireless assistance with the Terry Anatomical Collection. We benefited from the useful and detailed comments provided by Jim Wood, Emőke Szathmáry, Rob Hoppa, Darryl Holman, and the anonymous reviewers. LITERATURE CITED Agresti A. 1996. An introduction to categorical data analysis. New York: John Wiley and Sons. Aykroyd RG, Lucy D, Pollard AM, Roberts CA. 1999. Nasty, brutish, but not necessarily short: a reconsideration of the statistical methods used to calculate age at death from adult human skeletal and dental age indicators. Am Antiq 64:55–70. Bedford ME, Russell KF, Lovejoy CO, Meindl RS, Simpson SW, Stuart-Macadam PL. 1993. Test of the multifactorial aging method using skeletons with known ages-at-death from the Grant collection. Am J Phys Anthropol 91:287–297. Bocquet-Appel J-P. 1994. Estimating the average for an unknown age distribution in anthropology. In: Borgognini-Tarli S, Di Bacco M, Pacciani E, editors. Statistical tools in human biology. Singapore: World Scientific. p 197–202. Bocquet-Appel J-P, Bacro JN. 1997. Brief communication: estimates of some demographic parameters in a Neolithic rock-cut chamber (approximately 2000 BC) using iterative techniques for aging and demographic estimators. Am J Phys Anthropol 102:569 –575. Bocquet-Appel J-P, Masset C. 1982. Farewell to paleodemography. J Hum Evol 11:321–333. Bocquet-Appel J-P, Masset C. 1996. Paleodemography: expectancy and false hope. Am J Phys Anthropol 99:571–583. Bogin B, Smith BH. 2000. Evolution of the human life cycle. In: Stinson S, Bogin B, Huss-Ashmore R, O’Rourke D, editors. Human biology: an evolutionary and biocultural perspective. New York: Wiley-Liss. p 377– 424. Boldsen JL, Milner GR, Konigsberg LW, Wood JW. 2002. Transition analysis: a new means of estimating age from skeletons. In: Hoppa R, Vaupel J, editors. Paleodemography: age distributions from skeletal samples. Cambridge Studies in Biological and Evolutionary Anthropology 31. New York: Cambridge University Press. p 73–106. Brown PJ, Sundberg R. 1987. Confidence and conflict in multivariate calibration. J R Stat Soc [B] 49:46 –57. Clark WG. 1981. Restricted least-squares estimates of age composition from length composition. Can J Fish Aquat Sci 38:297– 307. Coale AJ, Demeny P. 1966. Regional model life tables and stable populations. Princeton, NJ: Princeton University Press. Cribari-Neto F, Zarkos SG. 1999. R: Yet another econometric programming environment. J Appl Econometrics 14:319 –329. Demetrius L. 1978. Adaptive value, entropy, and survivorship curves. Nature 257:213–214. Demetrius L. 1979. Relations between demographic parameters. Demography 16:329 –338. Dempster AP, Laird NM, Rubin DB. 1977. Maximum likelihood from incomplete data via the EM algorithm. J R Stat Soc [B] 39:1–38. Gage TB. 1988. Mathematical hazard models of mortality: an alternative to model life tables. Am J Phys Anthropol 76:429 – 441. Gage TB. 1998. The comparative demography of primates: with some comments on the evolution of life histories. Annu Rev Anthropol 27:197–221. Gelman A, Carlin JB, Stern HS, Rubin DB. 1995. Bayesian data analysis. New York: Chapman & Hall. Goldman N, Lord G. 1986. A new look at entropy and the life table. Demography 23:275–282. Hoenig JM, Heisey DM. 1987. Use of a log-linear model with the EM algorithm to correct estimates of stock composition and convert length to age. Trans Am Fish Soc 116:232–243. DECONSTRUCTING DEATH IN PALEODEMOGRAPHY Holman DJ, Wood JW, O’Connor KA. 2002. Estimating age-atdeath distributions from skeletal samples: a multivariate latent trait approach. In: Hoppa R, Vaupel J, editors. Paleodemography: age distributions from skeletal samples. Cambridge Studies in Biological and Evolutionary Anthropology 31. New York: Cambridge University Press. p 193–221. Hoppa RD. 1996. Representativeness and bias in cemetery samples: implications for paleodemographic reconstruction of past populations. Ph.D. dissertation, McMaster University, Hamilton, Ontario, Canada. Hoppa RD, Vaupel JW, editors. 2002. Paleodemography: age distributions from skeletal samples. Cambridge Studies in Biological and Evolutionary Anthropology 31. New York: Cambridge University Press. Howell N. 1976. Toward a uniformitarian theory of human paleodemography. In: Ward R, Weiss K, editors. The demographic evolution of human populations. New York: Academic Press. p 25– 40. Jackes MK. 1985. Pubic symphysis age distributions. Am J Phys Anthropol 68:281–299. Jackes M. 2000. Building the bases for paleodemographic analysis: adult age estimation. In: Katzenberg M, Saunders S, editors. Biological anthropology of the human skeleton. New York: Wiley-Liss. p 417– 466. Kennedy JJ. 1992. Analyzing qualitative data: log-linear analysis for behavioral research, 2nd ed. New York: Praeger Publishers. Kimura DK, Chikuni S. 1987. Mixtures of empirical distributions: an iterative application of the age-length key. Biometrics 43: 23–35. Konigsberg LW, Frankenberg SR. 1992. Estimation of age structure in anthropological demography. Am J Phys Anthropol 89:235–256. Konigsberg LW, Frankenberg SR. 1994. Paleodemography: “not quite dead.” Evol Anthropol 3:92–105. Konigsberg LW, Herrmann NP. 2002. Markov chain Monte Carlo estimation of hazard model parameters in paleodemography. In: Hoppa R, Vaupel J, editors. Paleodemography: age distributions from skeletal samples. Cambridge Studies in Biological and Evolutionary Anthropology 31. New York: Cambridge University Press. p 222–242. Konigsberg L, Holman D. 1999. Estimation of age at death from dental emergence and implications for studies of prehistoric somatic growth. In: Hoppa R, FitzGerald C, editors. Human growth in the past: studies from bones and teeth. New York: Cambridge University Press. p 264 –289. Konigsberg LW, Hens SM, Jantz LM, Jungers WL. 1998. Stature estimation and calibration: Bayesian and maximum likelihood perspectives in physical anthropology. Yrbk Phys Anthropol 41:65–92. Lee PM. 1989. Bayesian statistics: an introduction. New York: Edward Arnold. 309 Lindsey JK. 1996. Parametric statistical inference. New York: Oxford University Press. Love B, Müeller H-G. 2002. A solution to the problem of obtaining a mortality schedule for paleodemographic data. In: Hoppa R, Vaupel J, editors. Paleodemography: age distributions from skeletal samples. Cambridge Studies in Biological and Evolutionary Anthropology 31. New York: Cambridge University Press. p 181–192. Lovejoy CO, Meindl RS, Mensforth RP, Barton TJ. 1985. Multifactorial determination of skeletal age at death: a method and blind tests of its accuracy. Am J Phys Anthropol 68:1–14. McKern TW, Stewart TD. 1957. Skeletal age changes in young American males. Quartermaster Research and Development Command Technical Report EP-45. Meindl RS, Lovejoy CO, Mensforth RP, Walker RA. 1985. A revised method of age determination using the os pubis, with a review and tests of accuracy of other current methods of pubic symphyseal aging. Am J Phys Anthropol 68:29 – 45. Meindl RS, Russell KF, Lovejoy CO. 1990. Reliability of age at death in the Hamann-Todd collection: validity of subselection procedures used in blind tests of the summary age technique. Am J Phys Anthropol 83:349 –357. Mensforth RP. 1990. Paleodemography of the Carlston Annis (Bt-5) Late Archaic skeletal population. Am J Phys Anthropol 82:81–99. Milner GR, Wood JW, Boldsen JL. 2000. Paleodemography. In: Katzenberg M, Saunders S, editors. Biological anthropology of the human skeleton. New York: Wiley-Liss. p 467– 497. Müller H-G, Love B, Hoppa RD. 2002. Semiparametric method for estimating paleodemographic profiles from age indicator data. Am J Phys Anthropol 117:1–14. Nelder JA, Mead R. 1965. A simplex algorithm for function minimization. Comput J 7:308 –313. Paine RR, editor. 1997. Integrating archaeological demography: multidisciplinary approaches to prehistoric population. Carbondale, IL: Center for Archaeological Investigations, SIU. Petersen W. 1975. A demographer’s view of prehistoric demography. Curr Anthropol 16:227–245. Ripley BD. 2001. The R project in statistical computing. MSOR Connections 1:23–25 (available from mathstore.ac.uk/ newsletter/). Shao J. 1999. Mathematical statistics. New York: Springer-Verlag. Vaupel JW. 1986. How change in age-specific mortality affects life expectancy. Popul Stud 40:147–157. Wood JW, Holman DJ, O’Connor KA, Ferrell RJ. 2002. Mortality models for paleodemography. In: Hoppa R, Vaupel J, editors. Paleodemography: age distributions from skeletal samples. Cambridge Studies in Biological and Evolutionary Anthropology 31. New York: Cambridge University Press. p 129 –168.