Estimation of the most likely number of individuals from commingled human skeletal remains.код для вставкиСкачать
AMERICAN JOURNAL OF PHYSICAL ANTHROPOLOGY 125:138 –151 (2004) Estimation of the Most Likely Number of Individuals From Commingled Human Skeletal Remains Bradley J. Adams1 and Lyle W. Konigsberg2* 1 2 Central Identiﬁcation Laboratory, Hickam Air Force Base, Oahu, Hawaii 96853 Department of Anthropology, University of Tennessee, Knoxville, Tennessee 37996-0720 KEY WORDS Lincoln index; MNI; forensic anthropology; skeletal biology; pair-matching; recovery probability ABSTRACT This study examines quantiﬁcation techniques applicable to human skeletal remains, and in particular the Lincoln index (LI), the minimum number of individuals (MNI), and what we refer to as the most likely number of individuals (MLNI), which is a modiﬁcation of the LI by Chapman ( Univ. Calif. Publ. Stat. 1:131–159). As part of the study, a test of pair-matching between commingled homologous elements, e.g., right and left femora, was performed based on gross morphology. The results show that pair-matching can be accurately performed, and that the MLNI is a useful technique for dealing with well-preserved commingled remains recovered from archaeological excavations and/or forensic investigations. Our results show that it is potentially misleading to draw population A wide variety of quantiﬁcation techniques are available to the zooarchaeologist. However, physical anthropologists tend to focus almost exclusively on calculation of the minimum number of individuals, or MNI, as the best way of quantifying the number of individuals represented by commingled human remains (e.g., Ubelaker, 1974; Willey, 1990). The Lincoln index (LI) has been effectively used to quantify zooarchaeological samples and should be directly applicable to the analysis of commingled human remains. The key difference between these two quantiﬁcation techniques is that the LI estimates the original number of individuals represented by the osteological assemblage, while the MNI only estimates the number of individuals represented by the recovered assemblage. In cases involving significant taphonomic loss, the MNI may provide misleading number estimates. The LI, on the other hand, will provide a more accurate estimate of the original population size useful for paleodemographic or forensic purposes. In the following, we discuss the MNI, variants of the LI, and methods that can be applied to multiple paired elements to estimate the most likely number of individuals (MLNI) to have generated a commingled human skeletal assemblage. © 2004 WILEY-LISS, INC. conclusions based on the MNI, except in instances where recovery is near 100%. The MLNI was found to be the best method to compensate for the potential underestimates of the MNI and potential bias in the original LI estimates resulting from small sample sizes. We demonstrate the use of MLNI in estimating the number of individuals from Lodge 21 at the Larson site, a late protohistoric structure at which the inhabitants were massacred and subsequently had their skeletal elements commingled by further taphonomic processes. We also show how to calculate estimates and standard errors for the recovery probabilities of skeletal elements. Am J Phys Anthropol 125:138 –151, 2004. © 2004 Wiley-Liss, Inc. METHODS FOR QUANTIFYING NUMBERS OF INDIVIDUALS Minimum number of individuals The MNI is arguably the most popular method of quantiﬁcation in any type of commingled osteological analysis. Many researchers credit White (1953) with its initial use for abundance studies in archaeology. For interpreting population size from a skeletal assemblage, the MNI (as the name suggests) presents the minimum estimate for the number of individuals that contributed to the sample. In order to deal with fragmentary remains, speciﬁc segments of an element (e.g., distal femur) can be used for the Grant sponsor: National Science Foundation; Grant numbers: GS837, GS1635; Grant sponsor: National Geographic Society; Grant number: 699,912. *Correspondence to: Lyle Konigsberg, Department of Anthropology, University of Tennessee, 252 South Stadium Hall, Knoxville, TN 37996-0720. E-mail: firstname.lastname@example.org Received 8 November 2002; accepted 4 July 2003. DOI 10.1002/ajpa.10381 Published online 12 January 2004 in Wiley InterScience (www. interscience.wiley.com). QUANTIFICATION OF HUMAN REMAINS calculation of the MNI. Every fragment must share a unique landmark to ensure that fragments of the same element of the same individual are not counted as two distinct individuals. The basic principle of an MNI estimate is to avoid counting the same individual twice. There has been considerable debate in the zooarchaeological literature regarding the MNI and other quantiﬁcation techniques (e.g., Grayson, 1973, 1978, 1979, 1981, 1984). Many criticisms of the MNI for the analysis of zooarchaeological remains (and the recommended alternatives) are not germane to the study of human remains (e.g., differential transport and butchery of body parts). In most instances, human remains would have been initially deposited at a site in their entirety, rather than in parts. Regardless of the target population being studied, an important shortcoming of the MNI is that in many instances it will not provide an accurate estimate of the original number of individuals. In their discussion of faunal remains, Fieller and Turner (1982, p. 56) wrote that “the very presence of unmatched bones indicates that the MNI estimate is necessarily an underestimate of the number comprising the death assemblage.” The MNI simply states how many individuals would have been necessary to provide the recovered skeletal elements, but says nothing about the original death population. As we will show, the MNI varies depending on the recovery probability (i.e., the percentage of recovered elements), and for this reason it may be of limited value for the quantiﬁcation of osteological assemblages. Calculation of MNI There are three main variations for the calculation of the MNI, where L signiﬁes the number of left bones, R signiﬁes the number of right bones, and P signiﬁes the number of pairs. The MNI estimators are: Max 共L, R兲 (1) 共L ⫹ R兲/2 (2) L ⫹ R ⫺ P (3) The most common variant of the MNI used for the analysis of human remains is the Max(L, R) method. This value is obtained by sorting elements into lefts and rights, and then taking the greatest number as the estimate. It is equivalent to assuming that all of the less frequently observed bones are paired with the more frequently observed bones, and subsequently provides a rather unlikely estimate. The variant (L ⫹ R)/2 will usually provide the most minimal counts. It is an average for a paired element, and unless the lefts and rights are equal, the value will always be less than the maximum side. The L ⫹ R ⫺ P method usually provides a higher estimate than the other two options, because unpaired bones from different sides are assumed to come from dif- 139 ferent individuals. Horton (1984) referred to Equation (3) as the Grand minimum total. We will not deal extensively with MNI calculations in this paper. Under a simple model where each bone has a probability r of being recovered (which we refer to as the recovery probability) and there are N individuals, we would expect to recover Nr lefts, Nr rights, and Nr2 paired bones. Equations (1) and (2) then underestimate the actual N by (r ⫺ 1) ⫻ 100%, while Equation (3) underestimates the true N by (r(2 ⫺ r) ⫺ 1) ⫻ 100%. For example, when the recovery probability equals 0.7, Equations (1) and (2) underestimate the true N by 30%, while Equation (3) underestimates it by 9%. Equation (3) would appear preferable to Equations (1) and (2), but as we will see below, there are better methods than the MNI for estimating the actual number of individuals when considering paired elements. Further, because the recovery probability is unknown in most osteological assemblages, we have no real way of determining to what extent MNI methods underestimate the true N. This general trend in the MNI’s behavior indicates the need for caution when using MNI to quantify an abundance of human remains. Lincoln/Petersen index (LI) The Lincoln index (LI) has been almost exclusively applied to zooarchaeological remains and population studies of living animals (Allen and Guy, 1984; Borchers et al., 2002; Chapman, 1951; Chase and Hagaman, 1987; Fieller and Turner, 1982; Horton, 1984; Klein and Cruz-Uribe, 1984; Plug, 1984; Ringrose, 1993; Seber, 1973; Turner, 1980, 1983, 1984; Turner and Fieller, 1985; Wolter, 1990; and many others). The LI was ﬁrst developed for population studies of living animals based on capturerecapture techniques, and was later adapted for application to zooarchaeological faunal assemblages. Fishery biologists refer to it as “Petersen’s method,” which was used by C.G.J. Petersen in 1889, while ornithologists and mammalogists call it the “Lincoln index” due to its use by F.C. Lincoln in 1930 (LeCren, 1965). Although the distinction may be trivial, it is important to be aware that both names refer to the same technique. The same logic as used in the LI is also applied to human census data in order to estimate the number of missed individuals from a large census. In this setting, the method is referred to as “dual-system estimation,” and is used when there is both a census and a “post-enumeration survey” (Freedman, 1991). The key reason to use the LI for skeletal remains is that accurate estimates of the original population can be derived from samples in which taphonomic biasing has occurred. The theoretical basis of the formula is that one is studying populations in which all of the animals (living or dead) need not or cannot be observed. This is of particular utility with archaeological or paleoanthropological samples in which site sampling takes place, or in forensic situations 140 B.J. ADAMS AND L.W. KONIGSBERG when remains may have been exposed to taphonomic processes. The LI, as an inherent feature of its derivation, accounts for some degree of data loss, but it is important that the data loss occurs randomly. Ringrose (1993) stated that if both sides of an element were initially deposited on a site or transported to it, and neither was subsequently transported away, then the LI is an effective tool. Ringrose (1993, p. 129) concluded that this technique “is perhaps the best way we have of trying to get back to the Death Assemblage.” A similar conclusion was drawn by Fieller and Turner (1982, p. 59), who wrote that: “if matches can be made and if bias in selection has not operated, then the method is capable of reconstructing the original death assemblage size even though a number of selection processes may have operated. This is of particular importance to the archaeologist dealing with material which has been subjected to selection during initial deposition, preservation and eventual recovery.” Of course, estimates will be more precise with high levels of recovery, but it is extremely important to realize that this technique does account for random data loss. Calculation of the LI and MLNI from single elements When working with skeletal elements, the calculation of the Lincoln index requires a determination of element pair-matching. Pair-matching involves the comparison of right and left elements to decide if they are from a single individual. In theory, any paired element in the body could be used for the calculation of the LI. In practice, though, particular parts of the skeleton will be more useful than others. For example, teeth may not be appropriate for this technique due to the fact that asymmetric wear and antemortem loss resulting from decay or trauma could confound the process. Important criteria to consider when choosing appropriate skeletal elements are element size, the presence of distinct morphological traits, the potential for age and sex determination, and the likelihood of survival. Numerous researchers (e.g., Brain, 1976; Lyman, 1993, 1994; Waldron, 1987; Willey et al., 1997) found that the larger and denser bone elements have a better chance of survival than the delicate elements. With respect to human remains, some of the best bones to include would be the femur, tibia, humerus, and os coxa due to their morphology, survivability, and/or potential for providing signiﬁcant biological information. The bones from one side of the skeleton, e.g., left (L), are analogous to the initial capture stage in a capture-recapture study. The bones from the other side of the skeleton, e.g., right (R), are analogous to the second stage of the capture-recapture study. It is irrelevant which side of the skeleton is treated as the initial “catch.” The number of elements which can be matched from the right and left sides (P), indicating that they originated from the same indi- TABLE 1. Contingency table representations for sampling of paired elements1 a. Sampling layout assuming known population size Sampling of lefts ⫹ ⫺ Total Sampling of rights ⫹ ⫺ Total P R⫺P R L⫺P N⫺L⫺R⫹P N⫺R L N⫺L N b. Example discussed in text Sampling of rights Sampling of lefts ⫹ ⫺ Total ⫹ ⫺ Total 12 8 20 6 4 10 18 12 30 c. MNI approach based on observed numbers of lefts, rights, and pairs Sampling of lefts ⫹ ⫺ Total Sampling of rights ⫹ ⫺ Total P R⫺P R L⫺P ⱖ0 ⱖL ⫺ P L ⱖR ⫺ P ⱖL ⫹ R ⫺ P d. MNI approach based on observed lefts and rights (with L ⱖ R), ignoring the existence of pairs Sampling of rights Sampling of lefts ⫹ ⫺ Total ⫹ ⫺ Total R 0 R L⫺R ⱖ0 ⱖL ⫺ R L ⱖ0 ⱖL 1 (L,number of lefts; R, number of rights; P, number of pairs; N, population size. ⫹ and ⫺, observed (sampled) and unobserved (not sampled), respectively. vidual, is analogous to the recapture of initially tagged animals. An estimate of the original death assemblage N̂ represented by the skeletal elements is: N̂ ⫽ L ⫻ R , P (4) In the event that no pairs are discovered, Fieller and Turner (1982) and Turner (1983, 1984) suggested that the formula should be modiﬁed to: N̂ ⫽ 共L ⫹ 1兲共R ⫹ 1兲. (5) While the typical presentation of the LI is in terms of capture-recapture studies, we will use a contingency table presentation (Cormack, 1989) focused on the speciﬁc example of estimating the most likely number of individuals from paired elements. In a hypothetical death assemblage of N individuals, a number of right (R) and a number of left (L) elements of a particular bone will be recovered. Some of the elements can be matched and suggested to have derived from single individuals, where the number of such left with right pairings is P. Note that max(L, R) ⱖ P. If we know the original size of the population, a 2 ⫻ 2 contingency table can be constructed, as shown in Table 1a. As an example, if the recovered sample contained 18 left and 20 right QUANTIFICATION OF HUMAN REMAINS 141 bones, of which 12 are paired and the original population consisted of 30 individuals, the contingency table in Table 1b results. In Table 1b, the marginal totals present the numbers of left and right bones that were sampled or observed, and the numbers that were not sampled or observed. The fourth cell in the table gives the number of individuals from whom neither a right nor a left was sampled, and who were thus unobserved. If we assume independence for sampling of rights and lefts, then the expected number of paired bones from Table 1a is: E共P兲 ⫽ 冉 冊冉 冊 R N L N. N (6) Replacing the number of expected pairs with the observed number and solving for N gives the Lincoln index shown in Equation (4). Chapman (1951) referred to this estimator as N0, and noted that as it must take an integer value, the correct estimate is given by the integer that satisﬁes the following inequality: L ⫻ R L ⫻ R ⫺ 1 ⬍ N0 ⱕ , P P (7) where N0 is the integer value used in place of N̂. In an alternative scheme, where we only wish to ﬁnd the minimum number of individuals allowing for pairs, we assume that N ⫺ L ⫺ R ⫹ P is at least zero, as shown in Table 1c. The total (L ⫹ R ⫺ P) is then the MNI based on observed left, right, and paired bones (as in Equation (3)). Curiously, if we ignore the existence of pairs, then this is equivalent to saying that all of the bones from the less frequently occurring side must be paired with all of the bones from the more frequently occurring side. In this case, P ⫽ min(L, R), so that L ⫹ R ⫺ P becomes max(L, R) ⫹ min(L, R) ⫺ min(L, R) ⫽ max(L, R) (as in Equation (1)). The case where there are more lefts than rights is shown in Table 1d for this ﬁnal setting, where the MNI is simply L. The LI was shown to be the best asymptotically normal estimator for N, given L, R, and P (Chapman, 1951). However, under conditions of small sample size, the index is well-known to produce biased estimates (Bailey, 1951; Chapman, 1951; Robson and Regier, 1964). Chapman (1951) and Bailey (1951) both proposed modiﬁcations to the LI, but as Robson and Regier (1964) emphasized, these modiﬁed estimators have substantial negative biases when the number of recovered pairs is low. The variant by Chapman (1951) of the Lincoln index is: N* ⫽ 共L ⫹ 1兲共R ⫹ 1兲 ⫺ 1 , P ⫹ 1 (8) where the symbol means “ﬂoor” (i.e., truncation to the integer value). The estimator of Chapman (1951) is the value we will use to estimate the “most likely number of individuals” (MLNI). In order to demonstrate the relationship between the original LI (Equation (4)) and the estimator by Fig. 1. Comparison of percent bias for N̂ (original Lincoln index) and N* (modiﬁcation of LI by Chapman, 1951) against recovery probability when true N is equal to 50. Chapman (1951) (Equation (8)), Figure 1 shows a plot of the expected percent bias for each estimator. The plot shows the percent bias against recovery probability for a case where the true N is 50, and under the simplifying assumption that L ⫽ R. The graph was constructed assuming that the probability function for P conditional on L (and R, which is equal to L) is speciﬁed by the hypergeometric distribution. Figure 1 shows that at very low recovery probabilities, the original LI (N̂) ﬁrst underestimates the true N, and then overestimates the true N. In contrast, the estimator of Chapman (1951) (N*) initially underestimates the true N, but then rapidly becomes an unbiased estimator. Figure 1 demonstrates that the LI estimator of Chapman (1951) (Equation (8)) outperforms the traditional index, and given the simplicity of its calculation, it is preferable to both the traditional LI and any MNIbased approach. However, the estimator by Chapman (1951) is not immune to bias. Figure 2 shows the expectation for the estimator of Chapman (1951) across different population sizes (N ⫽ 15, 30, and 50) for various probabilities of bone recovery and for various expected numbers of pairs. The salient point in these graphs is that this estimator is unbiased, provided that the recovery probability is 0.5 or greater and/or there are ﬁve or more expected pairs. More correctly, the estimator is unbiased if L ⫹ R ⱖ N. Robson and Regier (1964) pointed out that the estimator of Chapman (1951) was unbiased if L ⫹ R ⱖ N, but as neither N nor the recovery probability is known in applications of the LI, it is necessary to rely on some other criterion. Robson and Regier (1964, p. 217) suggested that if the number of observed pairs is 7 or more, then “with 95 percent conﬁdence the bias is negligible.” This is born out in Figure 2, where there is no appreciable bias when there are 5 or more pairs (the lower number of 5 as opposed to 7 is again due to our assumption that L ⫽ R). The results of Robson and Regier (1964) indicate 142 B.J. ADAMS AND L.W. KONIGSBERG for N, which is used to specify exact intervals on estimates, and which allows us to extend the estimator of Chapman (1951) to multiple elements. To ﬁnd the probability function for N, we must write the probability of N conditional on L, R, and P. This probability can be found by combining the probability of getting P number of pairs conditional on L, R, and N with an uninformative prior for the recovery probability r. The probability of getting P number of pairs conditional on L, R, and N is found using the hypergeometric distribution (Burford, 1967; Evans et al., 1993; Feller, 1950; Freund, 1992; Harris and Stocker, 1998; Hays, 1988; Hogg and Tanis, 1983; Larson, 1969; Lee, 1997; Mendenhall et al., 1990; Mosteller et al., 1970; Neter et al., 1982; Nolan and Speed, 2000; Sokal and Rohlf, 1981; Wadsworth and Bryan, 1974; Weiss, 1961; Zehna, 1969), while the uninformative prior for the recovery probability is a uniform distribution between zero and one, speciﬁed by a beta(1,1) distribution (Evans et al., 1993; Freund, 1992; Harris and Stocker, 1998; Hogg and Tanis, 1983; Larson, 1969; Lee, 1997; Mendenhall et al., 1990; Mosteller and Tukey, 1977; Nolan and Speed, 2000; Wadsworth and Bryan, 1974). Using results from Roberts (1967, Equation 2.11; see also Lee, 1997, p. 200 –202), we write the probability of N conditional on L, R, and P as: Pr共N兩L,R,P兲 ⫽ Fig. 2. Comparison of percent bias for modiﬁcation by Chapman (1951) of LI graphed against recovery probability and expected number of pairs with N ⫽ 15, 30, and 50. that, provided there are at least 7 pairs present in an assemblage, the variant by Chapman (1951) of the LI (Equation (8), above) will have only a trivial negative sampling bias. In addition to estimates of the number of individuals, it would be useful to form conﬁdence intervals around these estimates. Seber (1973) discussed the calculation of the variance of N, which can be used to form approximate conﬁdence intervals. However, the distribution for estimates of N is generally asymmetric, and it is also a discrete distribution, making the conﬁdence intervals based on variance estimates misleading. The main advantage of the original LI N̂ and of the modiﬁed form by Chapman (1951) of the LI (N*) is that they can be easily and quickly calculated. If correct calculation of conﬁdence intervals is at issue, then it is far better to use a probabilistic model, as we show below. Estimation of the MLNI from the complete probability function Here, we derive the estimator of Chapman (1951) to use as the MLNI. The derivation is necessary in order to generate the complete probability function L!R!共N ⫺ R兲!共N ⫺ L兲! 共L ⫺ P兲!共R ⫺ P兲!共N ⫺ L ⫺ R ⫹ P兲! 共N ⫹ 1兲!共P ⫺ 1兲! ⫽ Pr共P兩L,R,N兲 ⫻ (9) P , 共N ⫹ 1兲 where Pr(P兩L,R,N) is the probability from the hypergeometric distribution of getting P pairs, upon drawing L lefts and R rights from N individuals. Equation (9) is easy to implement because many statistical packages as well as some spreadsheets include the hypergeometric probability function. We implemented Equation (9) in an Excel™ spreadsheet (available from http://konig.la.utk.edu/ MLNI.html). To ﬁnd the maximum likelihood estimator from (9), we write the ratio of adjacent probabilities as: 共N ⫺ R兲共N ⫺ L兲 Pr共N兩L,R,P兲 ⫽ . Pr共N ⫺ 1兩L,R,P兲 共N ⫺ L ⫺ R ⫹ P兲共N ⫹ 1兲 (10) For values of N less than the maximum likelihood estimate, this ratio will be greater than one, while for values of N greater than the maximum likelihood estimate, this ratio will be less than one. We can solve the inequality: 共N ⫺ R兲共N ⫺ L兲 ⱖ 共N ⫺ L ⫺ R ⫹ P兲共N ⫹ 1兲 (11) 143 QUANTIFICATION OF HUMAN REMAINS TABLE 2. Probability values from Equation (9) and cumulative probabilities, for L ⫽ 21, R ⫽ 15, and P ⫽ 11, tabled at various N values N P Cumulative 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 0.0457 0.0930 0.1196 0.1251 0.1167 0.1017 0.0847 0.0686 0.0545 0.0427 0.0332 0.0257 0.0198 0.0153 0.0118 0.0091 0.0071 0.0055 0.0043 0.0033 0.0026 0.0020 0.0016 0.0013 0.0010 0.0008 0.0457 0.1386 0.2582 0.3833 0.5000 0.6017 0.6864 0.7550 0.8094 0.8521 0.8853 0.9111 0.9309 0.9462 0.9580 0.9671 0.9742 0.9797 0.9839 0.9872 0.9898 0.9919 0.9935 0.9948 0.9958 0.9966 to ﬁnd that the maximum likelihood estimate occurs at: 共L ⫹ 1兲共R ⫹ 1兲 ⫺ 1 , P ⫹ 1 identical with Equation (8). We refer to Equation (8) as the MLNI because it gives the value of N with the highest posterior probability after starting from a uniform prior for N and a uniform prior for the recovery probability. Table 2 shows the probability values across N from Equation (9) for the case where L ⫽ 21, R ⫽ 15, and P ⫽ 11. These probability values were calculated in the Microsoft Excel™ spreadsheet at http:// konig.la.utk.edu/MLNI.html, and are shown from N ⫽ 25 (the MNI) up to N ⫽ 50. Table 2 also contains the cumulative sum of the probability function for N, from which we can establish intervals around the maximum likelihood estimate of N. Because N follows a discrete distribution, it is not usually possible to give customary conﬁdence intervals, such as 95% intervals. We can, however, determine the “highestdensity region” (HDR) (see Lee 1997, p. 49), which we prefer to classic conﬁdence intervals. The concept of an HDR is best treated graphically, which we do in Figure 3. At the top of Figure 3 we draw the frequency histogram for N from the probabilities in Table 2. We then move a horizontal line (dashed line in Fig. 3) down through the histogram until it comes to rest on a block of our choosing, and we hatch the histogram columns that terminate above this line. By doing so, we ﬁnd all values of N at which the probability is larger than the “baseline” value established by our horizontal line (hence, the “highest density region”). This region (Fig. 3, top) includes 70.9% of the distribution and runs from N ⫽ 26 to N ⫽ 32, inclusive. The percentage in the HDR is found as F(32) ⫺ F(25), where F() is the cumulative density up to and including the number in parentheses. A 70.9% HDR is rather small, so we drop the horizontal line even further, until we obtain Figure 3 (bottom). This second HDR runs from N ⫽ 25 (the MNI) to N ⫽ 38, inclusive, and contains 94.6% of the distribution (equal to F(38) ⫺ F(24), where F(24) ⫽ 0). From Figure 3, or the values presented in Table 2, we can see that the maximum likelihood estimate is at 28 individuals, since this is the highest part of the histogram. This is also the value we would obtain using Equations (7) or (8). As an alternative to direct calculation of the distribution for N (Equation (9)), George and Robert (1992) and Manly (1997) discussed Markov chain Monte Carlo (MCMC) estimators in mark-recapture settings. The edited volume by Gilks et al. (1996) gives a general account of MCMC methods. In addition to the spreadsheet calculations described above, we used a Metropolis algorithm (Gelman et al., 1995, p. 323–324) to sample from the distribution for N. We used Equation (10) to calculate the “acceptance probability” for proposed “jumps” to N ⫺ 1 or N ⫹ 1 in the Monte Carlo chain. This MCMC progresses rather slowly, so we cannot recommend it over the exact analytical calculations from Equation (9). Estimating the recovery probability In addition to estimating N, it may be desirable to estimate the “recovery probability,” or “r.” The recovery probability (r) is deﬁned here as the probability that a bone will make its way into the sample being analyzed. The recovery probability can be affected by myriad factors, including mortuary practice, taphonomy, and archaeological recovery methods. We use the term here simply to represent the ﬁnal probability of recovery in a sample. In the mark-recapture setting, the recovery probability is usually referred to as the “probability of capture” or “catchability.” Assuming that the probability of recovering a left bone is equal to that of recovering a right, the maximum likelihood estimate of r is: r̂ ⫽ 2P , L ⫹ R (12) and the asymptotic standard error of the estimate is: 冋 册 共r̂ ⫺ 1兲2 共r̂ ⫺ 2兲2 r̂2 s.e.共r̂兲 ⫽ 2 r̂ 共L ⫹ R兲共3 ⫺ 2r̂兲 ⫹ 2P共2 ⫺ 6r̂ ⫹ 3r̂2) 1 2 (13) As an example, with L ⫽ 21, R ⫽ 15, and P ⫽ 11, Equation (12) gives an estimated recovery probability of 0.6111, and Equation (13) gives a standard error of 0.0957. The calculations for the analytical estimate and standard error of the recovery proba- 144 B.J. ADAMS AND L.W. KONIGSBERG Fig. 3. Graphical representation of calculation of “highest-density regions” (HDR) for N when L ⫽ 21, R ⫽ 15, and P ⫽ 11. Top: 70.9% HDR. Bottom: 94.6% HDR. bility are available in the spreadsheet we provide at http://konig.la.utk.edu/MLNI.html, and the data in the spreadsheet can be edited to solve different problems. Extension of the MLNI to multiple elements The methods we described above can be easily extended to cover multiple paired elements. In the mark-recapture literature, such an application is most analogous to what is known as a mark-resight study. In mark-resight studies, animals are ﬁrst marked (generally with a radio transmitter), and then at later times, animals are tallied in a survey as to whether or not they have a radio transmitter (White and Garrott, 1990). The multiple resightings are analogous to having multiple elements available for estimating N. In the wildlife literature, N is estimated using a joint hypergeometric estimator (JHE) (Bartmann et al., 1987; Neal et al., 1993; White, 1993, White 1996; White and Garrott, 1990). In our setting, this means that we use the probability distribution for N (Equation (9)) across each element, multiply the probabilities across elements for each value of N, and then normalize the probabilities so that they sum to 1.0 (we give a speciﬁc example from the Larson site, below). The assumption in making this calculation is that the recovery of, say, femora is independent of the recovery of, say, tibiae. White (1996) wrote and made available the program NOREMARK (see http://www.cnr.colostate. edu/⬃gwhite/software.html) which ﬁts the JHE, and his program can also be used for the single-element estimation problem described above. His likelihood calculation is from the product rule applied to inde- QUANTIFICATION OF HUMAN REMAINS pendent hypergeometric distributions, and is consequently a generalization of the original Lincoln index. Our calculations for multiple elements are available in a spreadsheet from http://konig.la.utk.edu/MLNI.html, and serve as a generalization of the estimator of Chapman (1951) to the multiple element case. The use of multiple elements to estimate N introduces additional complications in the analysis, particularly because the assumption of independence in recovery cannot be tested. It is possible to use a likelihood ratio chi-square test to test for consistency of the estimator (see http://konig.la.utk.edu/ MLNI.html), but this test assumes independence in recovery of elements. Provided the likelihood ratio chi-square test is not signiﬁcant, the multiple element estimator stands a reasonable chance of being unbiased. However, if the assumption of independence is violated, then the HDR calculation will provide regions that are too narrow. APPLICATIONS OF MLNI AND RECOVERY PROBABILITY ESTIMATION Test for the accuracy of pair-matching As accurate pair-matching is critical in the estimation of the LI, MLNI, and recovery probabilities, it was necessary to carefully examine this facet of the technique. In order to examine the reliability of pair-matching bones of the same element type (e.g., right and left femora), an experimental test was performed using an archaeologically recovered assemblage of well-preserved human skeletal remains. The observation of gross morphological characteristics was the only criterion used in the pair-matching process. It should be noted that only pair-matching of the same element was examined. Morphological sorting between element types (e.g., associating a femur with a humerus) may be feasible when only a small number of individuals are commingled and there is signiﬁcant size variation between individuals, but in more complicated situations this technique can be highly subjective and problematic. It should be noted that morphological sorting between different element types is not necessary for the calculation of the MLNI or LI. In order to perform a veriﬁable test of gross pairmatching abilities, a random sample was drawn from an archaeological skeletal collection housed at the University of Tennessee, Knoxville. Larson Cemetery (39WW2) is an Arikara site associated with the Postcontact Variant of the Coalescent Tradition (Owsley, 1994; Owsley et al., 1977). The Arikara were a Caddoan-speaking tribe who constructed dome-shaped, earth-covered lodges (Bass and Rucker, 1976; Catlin, 1989; Denig, 1961). Initial occupancy at Larson is suspected to have begun about A.D. 1680 and to have ended by A.D. 1700 (Johnson, 1994). Larson Cemetery was excavated during the 1966 –1968 ﬁeld seasons under the direction of William M. Bass through the University of 145 Kansas. In total, 621 individuals were recovered, most as primary ﬂexed interments (Owsley et al., 1977). It is believed that this sample represents approximately 90% of the cemetery population (Bass and Rucker, 1976; Owsley and Bass, 1979). Based on the demographic analysis of the remains, the cemetery appears to have been the sole repository for the village dead (Bass and Rucker, 1976). Two randomly generated lists simulated a 60% recovery of tibiae, femora, and humeri from both adult and subadult primary burials, one derived from an original sample of 15 skeletons, and one from 30 skeletons. The samples drawn for the test consisted of both complete and fragmentary elements. The randomly selected elements were commingled on a table without any assistance by the observer (B.J.A.), and any identifying labels (e.g., catalog numbers) were covered. These procedures removed any chance of bias by the observer during the test. The observer also did not know the recovery probability or the original number of individuals represented by the remains. After an inventory of the test sample was completed, fragments from the same element were conjoined, and then pair-matching was performed. This task involved placing all specimens of a particular element (e.g., complete and fragmentary tibiae) on a table. The bones were then sorted by side, age, and robusticity, to assist in the conjoining of fragmentary remains and eventual pairing of elements. A pair-match can be determined based on general morphology and taphonomic indicators. Morphological indicators include robusticity, muscle markings, epiphyseal shape, bilateral expression of a periosteal reaction, and general symmetry between elements. Furthermore, biological information (e.g., age and sex) may also be useful for pair-matching. Taphonomic variables may include the state of preservation (e.g., degree of weathering and color), presence of burning, cut marks, and animal damage. Because of the variation that is possible from differential preservation, taphonomic indicators should not be weighted as heavily as gross morphological features for pair-matching. Results of the 15-individual group drawn from Larson Cemetery (Table 3) show that all femora, all tibiae, and all but one pair of humeri were matched correctly. The error with the humeri was due to the fact that one of the specimens exhibited signiﬁcant remodeling resulting from a probable dislocation. Table 3 also contains the realized recovery probabilities (L ⫹ R divided by two times the known N) and their standard errors from the binomial distribution. Table 3 also lists the estimated recovery probabilities and their standard errors from the observed counts of pairs and from the true counts. There is good agreement of the estimated recovery probabilities with the realized probabilities. The standard errors of the estimated recovery probabilities are slightly higher than for the realized recovery probabilities, but this is to be expected because of the 146 B.J. ADAMS AND L.W. KONIGSBERG TABLE 3. Results of 15-individual test for pair-matching TABLE 4. Results of 30-individual test for pair-matching a. Bone counts, realized recovery probabilities, and their standard errors a. Bone counts, realized recovery probabilities ,and their standard errors Element Lefts Rights Pairs Realized recovery probability Femur Tibia Humerus 9 7 7 9 7 9 6 4 3 (4)1 0.600 0.467 0.533 Realized standard error Element Lefts Rights 0.089 0.091 0.091 Femur Tibia Humerus 15 19 20 21 18 18 b. Estimates of recovery probabilities (r) and their standard errors, based on observed pair-matches and actual (correct) pair match for humerus Pairs Realized recovery probability Realized standard error 11 10 (11)1 11 (12)1 0.600 0.617 0.633 0.063 0.063 0.062 b. Estimates of recovery probabilities (r) and their standard errors, based on observed pair-matches and actual (correct) pair match for tibia and humerus Element r s.e. Element r s.e. Femur Tibia Humerus Humerus (correct) 0.667 0.571 0.375 0.500 0.128 0.158 0.154 0.153 Femur Tibia Tibia (correct) Humerus Humerus (correct) 0.611 0.541 0.595 0.579 0.632 0.096 0.099 0.096 0.095 0.092 1 Error was made, and number in parentheses is correct answer. Fig. 4. Probability function for N in 15-individual test case, where solid line is from observed number of pairs, and dashed line is from actual (true) number of pairs. additional uncertainty that arises from not knowing the value of N. Figure 4 shows the probability histogram for N for both the observed pair counts and the true pair counts. The maximum likelihood estimate (MLNI) for the number of individuals is 14 for the observed pair counts and 13 for the true pair counts (using the combined probability across bones from Equation (9)). In both cases, the true value of 15 individuals is within the (approximate) 95% HDR. For the observed pairs, the 95.7% HDR is from 13–19 individuals, while for the true pair count, the 96.3% HDR is from 12–17 individuals. Results of the blind pair-matching test involving 30 individuals (Table 4) revealed that one pair of humeri was missed, one pair of tibiae was missed, and all femora were correctly pair-matched. Again, there is good agreement of the estimated recovery probabilities with the realized probabilities, and again the standard errors of the estimated recovery 1 Error was made, and number in parentheses is correct answer. Fig. 5. Probability function for N in 30-individual test case, where solid line is from observed number of pairs, and dashed line is from actual (true) number of pairs. probabilities are slightly higher than for the realized recovery probabilities. The maximum likelihood estimate (MLNI) for the number of individuals is 31 (with a 94.8% HDR of 28 –36 individuals) for the observed pair counts, and 29 (with a 93.5% HDR of 27–33 individuals) for the true pair counts (see Fig. 5). Results of the pair-matching test derived from the Larson Cemetery sample indicate that pair-matching can be accurately performed for all of the elements examined. Furthermore, if errors in pairmatching are committed, they are more likely to occur from overlooking true pairs, as opposed to the pairing of unrelated elements. It is possible that this tendency may change if the sample size is substantially increased, since variation between individuals may not be as obvious. The main reason for difﬁculty in pair-matching is that fragmentation or weather- 147 QUANTIFICATION OF HUMAN REMAINS TABLE 5. Results of analysis of Larson Lodge 21 a. Element counts and estimates of N Element Lefts Rights Pairs Max(L, R) L⫹R⫺P Tibia Os Coxa Humerus Femur Overall 30 37 31 43 34 39 37 36 20 29 22 31 34 39 37 43 43 44 47 46 48 48 MLNI1 45, 47, 47, 48, 48, 50, 49, 51, 49, 50, 62 55 62 54 52 % HDR 94.8% 95.2% 94.6% 94.5% 95.2% b. Estimates of recovery probabilities (r) and their standard errors Element r s.e. Tibia Os coxa Humerus Femur 0.625 0.763 0.647 0.785 0.071 0.054 0.067 0.051 1 Central number in triplet is MLNI, while ﬁrst and last numbers give highest-density region (HDR). Column to immediate right gives percentage in HDR. ing, even minimal, can obliterate key areas used for identifying a match. Overall, the Larson Cemetery pair-matching test provided encouraging results. If pair-matching and conjoining are performed by an experienced osteologist, accurate results can be derived, especially when the original number of individuals is relatively small. Estimation of MNLI and recovery probabilities for Larson Village Lodge 21 A sample of commingled skeletal remains from Larson Village (39WW2) were utilized as a test for the validity of the quantiﬁcation methods described above. Larson Village is an Arikara habitation site that is associated with Larson Cemetery (described above, in Test for the Accuracy of Pair-Matching). Larson Village was excavated during the summers of 1963 and 1964 as part of an archaeological salvage program of the Smithsonian Institution River Basin Surveys. Although 29 circular lodge depressions were visible on the ground surface at the village site, only three (Lodges 1, 21, and 23) were completely excavated (Owsley et al., 1977). Discovery of scattered human remains on the lodge ﬂoors initially led to suspicion that villagers had been struck down by an epidemic disease (smallpox; see Bass and Rucker, 1976), but subsequent analysis revealed extensive perimortem trauma suggesting warfare (Owsley et al., 1977). The state of the skeletons encountered within the lodges implied that intentional burial did not occur. Evidence suggested that the individuals were either massacred and placed within the lodges, or they were left to decompose inside the lodges where they had been killed. The massacre at Larson Village is associated with the terminal occupation period of the site (A.D. 1700). Collapse of the lodge roofs buried the remains, but considerable commingling had occurred prior to this sealing of deposits. Owsley et al. (1977, p. 121) stated, “Due to disturbance and mixture of the skeletons . . . it was necessary to treat the bones as if they came from an ossuary.” The majority of the remains were discovered on the ﬂoor of Lodge 21, which had an MNI calculated at 44 individuals (Owsley et al., 1977). In order to estimate the MLNI and recovery probabilities for Lodge 21, four paired elements were selected for study. Sorting and pair-matching procedures were performed on the humeri, os coxae, femora, and tibiae. Counts of left, right, and paired elements are given in Table 5. The MNI results from the current study (Table 5) are reasonably consistent with MNI totals obtained by Owsley et al. (1977, p. 121), which were derived by “counting the major bones, and dividing them into rights and lefts.” This suggests that the Max(L, R) variant of the MNI was employed. Sorting of the femora as part of the current study revealed an MNI of 43. It is suspected that this difference may be due to the fact that a complete inventory of the entire skeleton was not performed for the present research; Owsley et al. (1977) did not indicate which element gave them their MNI of 44. Table 5 also lists the MNI values (from L ⫹ R ⫺ P) and the MLNI values for each element. The largest MNI value, 48 individuals, was derived from the femur. The MLNI values are 50 from the tibiae, 49 from the os coxae, 51 from the humerii, and 49 from the femora. Figure 6 graphs the probability distributions across N for each bone (from Equation (9)), as well as the combined probability across all four elements. Figure 6 shows the MLNI across all elements at 50 individuals, and the 95.2% HDR is from 48 –52 individuals (see spreadsheet at http://konig.la.utk.edu/MLNI.html for the calculations). The likelihood ratio chi-square comparing the MLNI estimates from each element (50, 49, 51, and 49 from the tibia, os coxa, humerus, and femur, respectively) to the overall estimate of 50 is 0.3634, which with three degrees of freedom (number of elements minus one) yields a probability value equal to 0.9477. Each element therefore provides a consistent estimate of the number of individuals. Table 5 also contains the estimated recovery probabilities and their standard errors. The estimated 148 B.J. ADAMS AND L.W. KONIGSBERG Fig. 6. Graph of probability function for N from tibiae, os coxae, humerii, and femora from Larson Lodge 21. Graph also shows combined probability function for all bones. recovery probabilities are 0.625, 0.763, 0.647, and 0.785 for the tibia, os coxa, humerus, and femur, respectively, which represent fairly high recovery probabilities. Indeed, from these probabilities the probability of obtaining no tibiae, os coxae, humerii, or femora from an individual would be about 4.55 ⫻ 10⫺5. Figure 7 shows a comparison of the distribution function for the recovery probability (the cumulative probability), assuming a normal distribution against 1,000 nonparametric bootstrap samples (Borchers et al., 2002, p. 113) for each of the elements. The agreement between our parametric approach (Equations (12) and (13)) and the nonparametric bootstrap is quite good, indicating that parametric estimation can be used in place of the more tedious bootstrap. Results from the Larson Village analysis suggest that in situations of good skeletal preservation, near complete recovery, and deposition resulting from a single event, the MNI estimates will accurately reﬂect the original population numbers. Similarly, in this type of situation, MLNI estimates should also be accurate. Agreement between the two techniques can be used to suggest that a high degree of recovery was achieved. However, it is considered advisable to compute the MLNI whenever feasible, since it will be much more accurate than the MNI when recovery probabilities drop. DISCUSSION Data loss from skeletal assemblages can occur as a result of many different taphonomic forces, including disarticulation, dispersal, bone weathering, and fragmentation (Lyman, 1987, 1994; Haglund and Sorg, 1997, 2002). Furthermore, recovery procedures, curational factors, and analytical issues need to be considered as possible taphonomic forces that can potentially bias both ancient archaeological sites, as well as more recent forensic contexts. Any study that attempts to draw conclusions related to numbers of individuals based on an analysis of commingled human remains needs to implement the most appropriate quantiﬁcation techniques available. This study compared two different quantiﬁcation techniques, the MNI and the LI, as well as a generalization of the LI that we refer to as the “most likely number of individuals,” or MLNI. With respect to the analysis of human skeletal remains, typically the only method considered for estimating the original population size has been the MNI. The results described here show that the LI and MLNI are better estimators of the original death assemblage. The strength of these techniques lies in the fact that they can account for taphonomic biases common to both archaeological samples and forensic contexts, with negligible effects upon the estimate. Rogers (2000) described another maximum likelihood estimator (“analysis of bone counts by maximum likelihood,” or “abcml”) that can be used to estimate the MLNI and taphonomic parameters of commingled zooarchaeological skeletal remains. His method could presumably be applied to samples of human remains. However, “abcml” ﬁrst requires external estimates of recovery probabilities in order to generate probabilities of what Rogers (2000) referred to as “conﬁgurations.” In fact, we could have applied “abcml” using our internally estimated recovery probabilities, but this would be a completely circular analysis. In the future, if recovery probabilities are shown to have some generality for human skeletal remains, then they could be used as ﬁxed parameters in applications of the “abcml” method to commingled human remains. Further, the “abcml” method has the distinct advantage over other quantiﬁcation methods in that it is calculated directly from bone counts, and does not require the timeintensive and often subjective demands of element pair-matching. We have not yet discussed variation in recovery probability across individuals, but it may be the case that there is variation in this parameter across an archaeological site. For example, the skeletal material may have been recovered from an ossuary pit where part of the pit was ﬁlled with a heavy, poorly drained clay, while the other part was ﬁlled with a well-drained loess. We would expect poor preservation within the clay with a concomitant low recovery probability, and good preservation within the loess with a concomitant high recovery probability. The effect that this variation has on the LI depends on how these different regimes correlate with individuals. For example, if there is complete correlation, such that all individuals were exposed to only one depositional environment, then the LI should be calculated for each context and summed. On the other hand, if the depositional environment is ran- QUANTIFICATION OF HUMAN REMAINS 149 Fig. 7. Distribution function for recovery probabilities of tibiae, os coxae, humerii, and femora from Larson Lodge 21. Smooth lines are from cumulative normal distributions with means and standard deviations from Equations (12) and (13), while step functions are empirical cumulative distributions from nonparametric bootstraps with 1,000 samples drawn for each graph. dom with respect to individuals, then the LI should be calculated for the entire site. If, for example, the fact that the left femur from an individual ended up in the clay does not affect whether the right femur ended up in the clay, then ignoring the heterogeneity in recovery/preservation has no effect on the LI calculated for the entire site. It is only when there is an intermediate correlation that estimation of the LI is subject to unknowable error. Due primarily to precedence and ease of calculation, the MNI will likely continue to serve as the quantiﬁer of choice with commingled human remains. However, as demonstrated here, superior methods of quantiﬁcation are available. The LI and MLNI provide a more realistic reconstruction of population counts from commingled samples, and supply researchers with more accurate estimates from which to present demographic results. It is suggested that the MNI should be reported alongside the LI or MLNI. Similar results between techniques substantiate the reliability of estimates and indicate a high recovery probability. Discrepancies between methods demonstrate the MNI’s correlation with the recovery probability and its inability to estimate the original population when data loss has occurred. The LI’s and MLNI’s behavior is not as susceptible to variable recovery probabilities, and their estimates provide a more accurate reconstruction of the original numbers of dead. 150 B.J. ADAMS AND L.W. KONIGSBERG Overall, positive results were found for the reliability of pair-matching from selected skeletal elements. It should be stressed that the pair-matching of homologous elements (e.g., right and left femora) was tested, and not the association of different element types with each other. If fragmentation is extensive or preservation is extremely poor, so that accurate pair-matches are impossible to determine, the LI and MLNI are prone to gross miscalculations due to the multiplicative nature of the procedures. In this type of situation, no method (save for “abcml,” with its own stringent data requirements) can provide accurate estimates of original populations. In situations of reasonable recovery and good preservation, the LI or MLNI provides an accurate assessment of the actual population that contributed to the recovered osteological assemblage. The total population can be divided into subgroups (e.g., age cohorts), using modiﬁcations to the LI or MLNI for more detailed demographic reconstructions. It is also possible to apply HDRs to these estimates, a feature not available with some other techniques. As recovery probabilities surpass 50%, the HDRs can become quite small. Overall, the LI and MLNI provide results that are of much more value for interpretation than any other quantiﬁcation technique currently employed. ACKNOWLEDGMENTS We thank William M. Bass, Richard L. Jantz, and Walter E. Klippel for their input on B.J.A.’s M.A. thesis, which provided the genesis for this paper. We also thank Gary C. White for providing L.W.K. with the source code for the program NOREMARK. Finally, we thank the anonymous reviews for their helpful comments and Susan R. Frankenberg for her comments on the ﬁnal draft. The Larson Cemetery was excavated under grants from the National Science Foundation (GS837 and GS1635) and a grant from the National Geographic Society (699,912), all to William M. Bass. LITERATURE CITED Allen J, Guy JBM. 1984. Optimal estimations of individuals in archaeological faunal assemblages: how minimal is the MNI? Archaeol Oceania 19:41– 47. Bailey NTJ. 1951. On estimating the size of mobile populations from recapture data. Biometrika 38:293–306. Bartmann RM, White GC, Carpenter LH, Garrott RA. 1987. Aerial mark-recapture estimates of conﬁned mule deer in pinyon-juniper woodlands. J Wildl Manage 51:41– 46. Bass WM, Rucker MD. 1976. Preliminary investigation of artifact association in an Arikara Cemetery (Larson Site), Walworth County, South Dakota. Washington, DC: National Geographic research reports, 1968 projects. p. 33– 48. Borchers DL, Buckland ST, Zucchini W. 2002. Estimating animal abundance: closed populations. New York: Springer. Brain CK. 1976. Some principles in the interpretation of bone accumulations associated with man. In: Isaac GL, McCown ER, editors. Human origins. Menlo Park: W.A. Benjamin. p 97–116. Burford RL. 1967. Introduction to ﬁnite probability. Columbus, OH: Charles E. Merrill. Catlin G. 1989. North American Indians. New York: Viking. Chapman DG. 1951. Some properties of the hypergeometric distribution with applications to zoological sample census. Univ Calif Publ Stat 1:131–159. Chase PG, Hagaman RM. 1987. Minimum number of individuals and its alternatives: a probability theory perspective. Ossa 13:75– 86. Cormack RM. 1989. Log-linear models for capture-recapture. Biometrics 45:395– 413. Denig ET. 1961. Five Indian tribes of the Upper Missouri. Norman: University of Oklahoma Press. Evans M, Hastings N, Peacock B. 1993. Statistical distributions, 2nd ed. New York: John Wiley and Sons. Feller W. 1950. An introduction to probability theory and its applications, volume 1. New York: John Wiley and Sons. Fieller NRJ, Turner A. 1982. Number estimation in vertebrate samples. J Archaeol Sci 9:49 – 62. Freedman DA. 1991. Adjusting the 1990 census. Science 252: 1233–1236. Freund JE. 1992. Mathematical statistics, 5th ed. Englewood Cliff, NJ: Prentice-Hall. Gelman A, Carlin JB, Stern HS, Rubin DB. 1995. Bayesian data analysis. New York: Chapman and Hall Press. George EI, Robert CP. 1992. Capture-recapture estimation via Gibbs sampling. Biometrika 79:677– 683. Gilks WR, Richardson S, Spiegelhalter DJ, editors. 1996. Markov chain Monte Carlo in practice. New York: Chapman and Hall Press. Grayson DK. 1973. On the methodology of faunal analysis. Am Antiq 39:432– 439. Grayson DK. 1978. Minimum numbers and sample size in vertebrate faunal analysis. Am Antiq 43:53– 65. Grayson DK. 1979. On the quantiﬁcation of vertebrate archaeofaunas. In: Schiffer M, editor. Advances in archaeological methods and theory, vol 2. p. 199 –237. Grayson DK. 1981. A critical view of the use of archaeological vertebrates in paleoenvironmental reconstruction. J Ethnobiol 1:28 –38. Grayson DK. 1984. Quantitative zooarchaeology. New York: Academic Press, Inc. Haglund WD, Sorg MH, editors. 1997. Forensic taphonomy: the postmortem fate of human remains. Boca Raton: CRC Press. Haglund WD, Sorg MH, editors. 2002. Advances in forensic taphonomy: method, theory, and archaeological perspectives. Boca Raton: CRC Press. Harris JW, Stocker H. 1998. Handbook of mathematics and computational science. New York: Springer-Verlag. Hays WL. 1988. Statistics, 4th ed. Chicago: Holt, Rinehart and Winston. Hogg RV, Tanis EA. 1983. Probability and statistical inference, 2nd ed. New York: Macmillan. Horton DR. 1984. Minimum numbers: a consideration. J Archaeol Sci 11:255–271. Johnson CM. 1994. A chronology of Middle Missouri Plains village sites. Draft report prepared for the Department of Anthropology. Washington, DC: National Museum of Natural History, Smithsonian Institution. Klein RG, Cruz-Uribe K. 1984. The analysis of animal bones from archeological sites. Chicago: University of Chicago Press. Larson HJ. 1969. Introduction to probability theory and statistical inference. New York: John Wiley and Sons. LeCren ED. 1965. A note on the history of mark-recapture population estimates. J Anim Ecol 34:453– 454. Lee PM. 1997. Bayesian statistics: an introduction, 2nd ed. New York: Oxford University Press. Lyman RL. 1987. Zooarchaeology and taphonomy: a general consideration. J Ethnobiol 7:93–117. Lyman RL. 1993. Density-mediated attrition of bone assemblages: new insights. In: Hudson J, editor. From bones to behavior. Carbondale: Center for Archaeological Investigations. p 324 – 341. Lyman RL. 1994. Vertebrate taphonomy. Cambridge: Cambridge University Press. Manly BFJ. 1997. Randomization, bootstrap and Monte Carlo methods in biology. New York: Chapman and Hall. QUANTIFICATION OF HUMAN REMAINS Mendenhall W, Wackerly DD, Scheaffer RL. 1990. Mathematical statistics with applications, 4th ed. Belmont, CA: Duxbury. Mosteller F, Tukey JW. 1977. Data analysis and regression: a second course in statistics. Reading, MA: Addison-Wesley. Mosteller F, Rourke REK, Thomas GB Jr. 1970. Probability with statistical applications, 2nd ed. Reading, MA: Addison-Wesley. Neal AK, White GC, Gill RB, Reed DF, Olterman JH. 1993. Evaluation of mark-resight model assumptions for estimating mountain sheep numbers. J Wildl Manage 57:436 – 450. Neter J, Wasserman W, Whitmore GA. 1982. Applied statistics, 2nd ed. Boston: Allyn and Bacon. Nolan D, Speed T. 2000. Stat labs: mathematical statistics through applications. New York: Springer. Owsley DW. 1994. Warfare in coalescent tradition populations of the Northern Plains. In: Owsley DW, Jantz RL, editors. Skeletal biology in the Great Plains. Washington, DC: Smithsonian Institution Press. p 333–344. Owsley DW, Bass WM. 1979. A demographic analysis of skeletons from the Larson Site (39WW2), Walworth County, South Dakota: vital statistics. Am J Phys Anthropol 51:145–154. Owsley DW, Berryman HE, Bass WM. 1977. Demographic and osteological evidence for warfare at the Larson Site, South Dakota. Plains Anthropol 13:119 –131. Plug I. 1984. MNI counts, pits and features. In: Hall M, Avery G, Avery DM, Wilson ML, Humphreys AJB, editors. Frontiers: southern African archaeology today. Oxford: Oxford University Press. p 357–362. Ringrose TJ. 1993. Bone counts and statistics: a critique. J Archaeol Sci 20:121–157. Roberts HV. 1967. Information stopping rules and inference about population size. J Amer Statist Assoc 62:763–775. Robson DS, Regier HA. 1964. Sample size in Petersen markrecapture experiments. Trans Am Fish Soc 93:215–226. Rogers AR. 2000. Analysis of bone counts by maximum likelihood. J Archaeol Sci 27:111–125. Seber GAF. 1973. The estimation of animal abundance and related parameters. London: Grifﬁn. Sokal RR, Rohlf JF. 1981. Biometry: the principles and practice of statistics in biological research, 2nd ed. San Francisco: W.H. Freeman. Turner A. 1980. Minimum number estimation offers minimal insight in faunal analysis. Ossa 7:199 –201. Turner A. 1983. The quantiﬁcation of relative abundances in fossil and subfossil bone assemblages. Ann Trans Mus 33:311–321. 151 Turner A. 1984. Behavioural inferences based on frequencies in bone assemblages from archaeological sites. In: Hall M, Avery G, Avery DM, Wilson ML, Humphreys AJB, editors. Frontiers: southern African archaeology today. Oxford: Oxford University Press. p 363–366. Turner A, Fieller NRJ. 1985. Considerations of minimum numbers: a response to Horton. J Archaeol Sci 12:477– 483. Ubelaker DH. 1974. Reconstruction of demographic proﬁles from ossuary skeletal samples: a case study from the Tidewater Potomoc. Smithsonian contributions to anthropology 18. Washington, DC: Smithsonian Institution Press. Wadsworth GP, Bryan JG. 1974. Applications of probability and random variables, 2nd ed. New York: McGraw-Hill. Waldron T. 1987. The relative survival of the human skeleton: implications for palaeopathology. In: Boddington A, Garland AN, Janaway RC, editors. Death, decay, and reconstruction. Manchester: Manchester University Press. p 55– 64. Weiss L. 1961. Statistical decision theory. New York: McGrawHill. White GC. 1993. Evaluation of radio tagging marking and sighting estimators of population size using Monte Carlo simulations. In: Lebreton JD, North PM, editors. Marked individuals in the study of bird populations. Basel, Switzerland: Birkhäuser Verlag. p 91–103. White GC. 1996. NOREMARK: population estimation from markresighting surveys. Wildl Soc Bull 24:50 –52 (see http://www. cnr.colostate.edu/⬃gwhite/software.html). White GC, Garrott RA. 1990. Analysis of wildlife radio-tracking data. New York: Academic Press. White TE. 1953. A method of calculating the dietary percentage of various food animals utilized by aboriginal peoples. Am Antiq 4:396 –398. Willey P. 1990. Prehistoric warfare on the Great Plains: skeletal analysis of the Crow Creek Massacre victims. New York: Garland. Willey P, Galloway A, Snyder L. 1997. Bone mineral density and survival of elements and element portions in the bones of the Crow Creek Massacre victims. Am J Phys Anthropol 104:513– 528. Wolter KM. 1990. Capture-recapture estimation in the presence of a known sex ratio. Biometrics 46:157–162. Zehna PW. 1969. Finite probability. Boston: Allyn and Bacon.