








































Forensic Anthropology Vol. 2, No. 2: 72–86
DOI: 10.5744/fa.2019.1000

© 2019 University of Florida Press

Introduction

In osteological analysis, commingled assemblages present 
a situation in which discrete sets of remains are not readily 
apparent. Commingled assemblages, such as ossuaries, are 
a fairly common situation in bioarchaeology (Curtin 2008; 
Herrmann & Devlin 2008; Ubelaker & Rife 2008; Willey 
1990). With the increasing utilization of forensic anthropol-
ogists in arenas such as mass disaster (Hinkes 1989; Mun-
dorff 2008, 2012; Sledzik & Rodriguez 2001), cremation 
litigation (Steadman et al. 2008), and human rights investi-
gations (Primorac et al. 1996; Varas & Leiva 2012), resolu-
tion of commingling is becoming commonplace (Adams & 
Byrd 2008, 2014). Forensic analysis of commingled remains 
focuses on victim identification and reassociating remains 
into discrete individuals (Adams & Byrd 2006, 2008, 2014; 
Byrd & Adams 2003, 2009). This focus has led to an increase 
in research on methodology for resolving commingling 
(Adams & Byrd 2008, 2014).

Of the methods available for resolving commingling, 
osteometric reassociation, which uses statistical models to 
compare bone dimensions, is considered a reliable and rela-
tively objective technique (Adams & Byrd 2006; Buikstra 
et al. 1984; Byrd 2008; Byrd & Adams 2003; Byrd & LeGarde 
2014; Konigsberg & Frankenberg 2013; O’Brien & Storlie 
2011; Rosing & Pischtschan 1995; Snow & Folk 1970). Tra-
ditional osteometric sorting logic is a decision- making, error- 
mitigation approach (Byrd 2008; Byrd & Adams 2003; Byrd & 
LeGarde 2014). This approach does not seek to reassociate 
elements per se; rather the analyst tests the null hypothesis 
that the dimensions of two bones are similar enough to have 
derived from the same individual (Adams & Byrd 2006; 
Byrd 2008; Byrd & Adams 2003; Byrd & LeGarde 2014). 
Possible matches are eliminated if the calculated p- value 
exceeds an analyst- defined threshold, or alpha level. Bones 
are reassociated if all other possible matches can be elimi-
nated. This approach implies that, because of broad variation 
in intra- individual bone size, reassociation is achievable 
via osteometrics when the assemblage represents a closed 
population of a smaller number of different- sized individu-
als (Byrd 2008).

The logic of reassociation through elimination was first 
introduced by Byrd and Adams (2003). A regression model 
and associated 90% prediction interval, based on the natural 
logarithm of the summed measurements by element, was 

RESEARCH ARTICLE

A Novel Method for Osteometric Reassociation Using Hamiltonian 
Markov Chain Monte Carlo (MCMC) Simulation

Kyle McCormicka*

ABSTRACT: Traditional osteometric reassociation uses an error- mitigation approach, which seeks to eliminate possible matches, 
rather than a predictive approach, where possible matches are directly compared. This study examines the utility of a Bayesian approach 
for resolving commingling by using a probabilistic framework to predict correct matches.

Comparisons were grouped into three types: paired elements, articulating elements, and other elements. Ten individuals were ran-
domly removed from the total sample (N = 833), acting as a small- scale, closed- population commingled assemblage. One element was cho-
sen as the independent variable, with the ten possible matching elements representing the dependent variable. A Bayesian regression 
model was constructed using the remaining total sample, resulting in a distribution of possible values that were smoothed into a probabil-
ity density, and probabilities were calculated. The element with the highest posterior probability was considered the best match. This pro-
cess was repeated 500 times for each comparison. The correct match was identified 51.60% of the time. Paired elements performed the 
best, at 80.76%, followed by 42.10% for articulating and 33.63% for other comparisons.

These results suggest that metric analysis of commingled assemblages is complex and that both elimination- based and prediction- 
based approaches have a role in resolving commingling. In this regard, the strength of a Bayesian approach is versatility, allowing for 
prediction of the correct match and elimination of possible matches, as well as integration of independent lines of evidence within one 
cohesive model.

KEYWORDS: forensic anthropology commingling, Bayesian modeling, MCMC simulation

aDefense POW/MIA Accounting Agency, DoD, Hawai’i, USA
*Correspondence to: Kyle McCormick, Defense POW/MIA 

Accounting Agency, DoD, 590 Moffett St., Joint Base Pearl  
Harbor– Hickam, AFB, Hawai’i 96853, USA 

E- mail: kmccormick9@yahoo  .com

Received 19 June 2018; Revised 17 October 2018;  
Accepted 26 October 2018



McCormick 73

constructed. If the bone in question fell outside of the 90% 
prediction interval, the researcher concluded that the ele-
ments are too different in size to be from one individual. 
The form of decision making used by Byrd and Adams (2003) 
follows a Neyman- Pearson approach to hypothesis testing, 
where decisions concerning the null hypothesis are strictly 
based on whether a test statistic passes an a priori threshold 
value (Royall 2000). The researcher is making a dichotomous 
decision whether to reject or fail to reject the null hypothe-
sis. Under this paradigm, there is no degree of belief in the 
null hypothesis— it is either rejected or it is not (Royall 1997). 
The explicit decision- making rationale and ease of interpre-
tation of this approach to science has obvious strengths. The 
elements in question either derive from the same individual 
or they do not; there are only two possible outcomes (Byrd 
2008).

Byrd (2008) provides a more nuanced statistical frame-
work and presents specific osteometric reassociation models 
for paired, articulating, and other element comparisons. 
Again, possible matches are eliminated by comparing a 
 p- value to an alpha level (ranging from 0.05 to 0.10, depending 
on the comparison type). Byrd (2008) also provides a means 
for aggregating multiple test results in more complex com-
mingling situations and introduces the severity principle, 
which focuses on identifying and mitigating error in decision 
making (Mayo & Spanos 2010). Decisions concerning the 
null hypothesis are based on the output of a statistical test. A 
researcher feels confident in his or her decision concerning 
a hypothesis if the test has a high chance of detecting the fal-
sity of the hypothesis (Mayo & Spanos 2010). Severity is 
used to incorporate the strength of evidence into the decision- 
making process.

This interpretative shift blends two forms of testing sta-
tistical hypotheses: Neyman- Pearson hypothesis testing and 
Fisherian significance testing (Lew 2013; Royall 1997). These 
approaches have different purposes: the former sets an a pri-
ori criterion (alpha level) for deciding between two compet-
ing hypotheses, while the latter attempts to interpret the 
strength of evidence against the null hypothesis. Most con-
temporary frequentists blend these two forms of hypothesis 
testing into a third formulation, sometimes referred to as 
rejection trials (Royall 1997). Rejection trials use an a priori 
alpha level as a decision- making criterion, similar to the 
Neyman- Pearson approach, but the researcher subjectively 
interprets the p- value as a measure of the strength of evidence 
against the null hypothesis (Royall 1997).

While this shift toward including additional information 
into the decision- making process increases subjectivity, it 
also increases rationality. The decision to reassociate a set of 
remains should be based on multiple lines of evidence, of 
which osteometric reassociation is just one (Byrd 2008). 
Incorporating multiple lines of evidence into a decision is a 
subjective process, based in part on the experience of the 

researcher. It matters if a p- value is 0.049 or 0.000001— the 
latter can be regarded as stronger evidence against the null 
hypothesis than the former.

This frequentist logic has obvious strengths. With a 
focus on the hypothetical frequency of a correct rejection, the 
results are highly reliable and easy to interpret, if any inter-
pretation is needed. There is, however, an obvious downside 
to this approach; it does not directly address the primary 
question of interest, namely, which bones are from the same 
individual? The sole reliance on eliminating possible matches 
is peculiar compared to the predictive nature of most other 
forms of osteological analysis (e.g., age, sex, ancestry, stat-
ure estimation).

To a Bayesian, probability is the numeric representation 
of the “degree of belief” in a proposition or set of proposi-
tions (Stark & Freedman 2003). This usage is more in line 
with a layperson’s understanding of probability than the fre-
quentist view of probability as long- run frequencies of an 
event. A Bayesian understanding of probability has shown 
promise for resolving commingling (Konigsberg & Franken-
berg 2013; McCormick 2016) and other aspects of forensic 
investigation (Brennaman et al. 2017; Jantz & Ousley 2005; 
Konigsberg & Frankenberg 2013).

One way to operationalize a Bayesian approach is to 
assign prior probabilities to each possible match, either 
through prior information or uninformed (equal) probabili-
ties. Prior probabilities are multiplied by the likelihood of the 
data to obtain a posterior probability, which is interpreted as 
the relative probability of a correct match after incorporating 
model information (McCormick 2016; Byrd & LeGarde 2014; 
Konigsberg & Frankenberg 2013).

Prior probability distributions can be assigned to the 
parameters used in estimating the model, such as the slope, 
y- intercept, and error term in linear regression. These prior 
distributions are used along with the likelihood function of 
the data to explore parameter space (possible values of the 
parameter) and to arrive at a posterior distribution for that 
parameter (Kéry 2010; S. M. Lynch 2007). Model parame-
ters are explicitly treated as distributions, instead of point 
estimates with uncertainty around that estimate, typically 
associated with frequentist modeling. The consequence of 
these different views of parameters is obvious in predictive 
modeling, such as linear regression. A frequentist regression 
model results in a single value for model parameters, includ-
ing the dependent (y) variable. Some form of interval esti-
mation (typically confidence and prediction intervals) is 
required to better understand the uncertainty in parameter 
point estimates. These intervals are not direct properties of 
the parameter and are not probabilistic statements that a 
parameter’s true value lies within a specified boundary 
(Hoekstra et al. 2014; Mayo 1982). Rather, prior to observing 
the data, a 95% confidence interval means there is a 95% 
chance that the interval will contain the true parameter 



74 Osteometric Reassociation Using MCMC Simulation

value (Hoekstra et al. 2014; Mayo 1982). After the data are 
observed, the true value is either within the interval or it is 
not. The interpretation of these intervals is based in a fre-
quentist understanding of probability, leading to pervasive 
misunderstanding. The osteometric sorting model of Byrd 
and Adams (2003), where possible matches were rejected if 
they fell outside of the prescribed prediction interval, is an 
example of such a misunderstanding. Byrd and Adams (2003) 
is best viewed as a shortcoming of a frequentist approach to 
problems of prediction rather than statistical acumen. While 
there are valid criticisms of Bayesian modeling, such as sub-
jectivity of prior information and, by extension, posterior dis-
tributions, as well as directed sampling strategies (Gelman 
2008), Bayesian modeling does not contain the interpretative 
pitfalls of a frequentist design. The ease of interpretation, 
handling of model parameters, and flexibility in model con-
struction are major differences between frequentist and 
Bayesian modeling and are perhaps the main benefits of a 
Bayesian approach. Bayesian interpretation and modeling 
has yet to be applied to resolving commingling. The current 
study examines the utility of such an approach to osteomet-
ric reassociation.

Materials

The data consist of 24 standard limb measurements from a 
total of 833 individuals curated at the William M. Bass 
Donated Skeletal Collection at the University of Tennessee, 
Knoxville. Individuals in the current study are predominantly 
European American adults, ranging in age from 18 to 
70 years at death (Fig. 1), a majority of which are male 
(males = 583, females = 250). The number of individuals var-
ies by comparison, as only those with complete measure-
ments for the compared elements were used.

FIG. 1—Age- at- death distribution of the sample (n = 883).

Methods

Measurements

The measurements used in this study are from the Forensic 
Anthropology Data Bank (FDB). Some interobserver vari-
ability is expected, given the multiple contributors to the 
FDB. Bivariate plots comparing left-  and right- side homolo-
gous measurements were used to identify and remove obvi-
ous outliers. The number of measurements varies by element 
(see Table 1). The number and quality of measurements 
should have an influence on reassociation. Elements with a 
large number of highly correlated variables should show 
the highest accuracy rates.

Reassociation Model

This study tests the accuracy of a Bayesian approach to osteo-
metric reassociation by simulating small- scale (n = 10) 
closed- population commingled assemblages and predicting 
the best match using standard osteological measurements and 
Bayesian regression. This process is repeated 500 times for 
each comparison. Accuracy is defined as the correct classi-
fication rate, or the number of times the best match is the cor-
rect match divided by 500.

Following Byrd (2008), limb element comparisons are 
grouped into three types: paired, articulating, and other ele-
ment comparisons (Table 2). By virtue of being antimeres, 
measurements are directly comparable between paired 
elements. For articulating and other comparison types, 

TABLE 1—Forensic Data Bank Measurements by Element.

Element Measurement Name FDB #

Humerus Maximum length v40
n = 5 Epicondylar breadth v41

Maximum vertical head diameter v42
Maximum diameter at midshaft v43
Minimum diameter at midshaft v44

Radius Maximum length v45
n = 3 A/P diameter at midshaft v46

Transverse diameter at midshaft v47

Ulna Maximum length v48
n = 4 Dorso- Volar diameter v49

Transverse diameter v50
Physiological length v51

Femur Maximum length v60
n = 8 Bicondylar length v61

Epicondylar breadth v62
Maximum diameter of head v63
A/P subtrochanteric diameter v64
Transverse subtrochanteric diameter v65
A/P diameter at midshaft v66
Transverse diameter at midshaft v67

Tibia Condylar- malar length v69
n = 4 Maximum proximal epiphyseal breadth v70

Distal epiphyseal breadth v71
Maximum diameter at nutrient foramen v72



McCormick 75

TABLE 2—Osteometric Comparisons by Type.

Comparison Type

Femur/Femur Paired
Humerus/Humerus Paired
Radius/Radius Paired
Tibia/Tibia Paired
Ulna/Ulna Paired
Femur/Tibia Articulating
Humerus/Ulna Articulating
Humerus/Radius Articulating
Ulna/Radius Articulating
Femur/Humerus Other
Femur/Ulna Other
Femur/Radius Other
Tibia/Humerus Other
Tibia/Ulna Other
Tibia/Radius Other

transformation of raw measurements is required to compare 
elements.

Partial Least Squares

Partial least squares (PLS) analysis is a class of techniques 
for data reduction and latent variable analysis (Boulestiex & 
Strimmer 2006; Chen & Hoo 2011; Haenlein & Kaplan 2004; 
Rosipal & Krämer 2006; Wegelin 2000). These techniques 
share a common method of extracting components— via 
ordinary least squares regression. PLS analysis is similar to 
principal component analysis (PCA) and canonical correla-
tion analysis (CCA), which extract orthogonal (uncorrelated) 
score vectors that are weighted composites of the original 
data set (Rosipal & Krämer 2006). Typically, the goal with 
any type of predictive data reduction analysis is twofold: (1) 
to find linear combinations that well represent the original 
variables and (2) to find highly correlated linear combina-
tions. Because PCA captures a maximum amount of varia-
tion from the original variables, it is an optimal solution to 
the first goal. In a predictive framework, where one block of 
variables is used to predict another block, PCA fails to 
achieve the second goal, because components between blocks 
of variables have no relationship. On the other hand, CCA 
optimally achieves the second goal by creating linear com-
binations of each block that are maximally correlated with 
one another. However, CCA fails at the first goal because 
these components are not designed to capture information or 
variance within a block and are based on the correlation 
matrix of raw variables, obscuring the biological meaning of 
components and making the interpretation of components 
difficult (Bookstein 1991; Wegelin 2000). Furthermore, CCA 
components are unstable in instances of multicollinearity, 
and solutions are not uniquely defined when the number of 
variables is large compared to the sample size (Wegelin 
2000). Simply, PCA explains variation within a block of vari-
ables and CCA explains variation between two blocks of 
variables. While not optimal, PLS achieves both goals by 

finding linear combinations of variables through the covari-
ance of raw variables that both capture variability and are 
highly correlated (Bookstein 1991; Wegelin 2000). Compo-
nents of the X- block (independent variables) are orthogonal, 
are good representations of X, and are good at explaining Y 
(dependent variables). Components of the Y- block are orthog-
onal, are good representations of Y, and are highly correlated 
with the X- block components. Stated another way, PLS mod-
els create components that predict a set of dependent vari-
ables from a set of independent variables that have the best 
predictive power on the dependent variables (Chen & Hoo 
2011). The package “plsdepot” (Sanchez 2016) was used in R 
(R Core Team 2015) to extract relevant PLS components.

Simulated Commingling

Ten individuals were randomly removed from the total data 
set. These 10 individuals act as a simulated commingled 
assemblage. One element is chosen as the independent (x) 
variable, with the 10 possible matching elements acting as the 
dependent (y) variable. For example, if we are interested in 
reassociating a left femur with 10 possible right femora, then 
the left femur is predicting the right femur. In this situation, 
the left femur is the independent variable and the right femur 
is the dependent variable. A left femur is selected from the 
commingled assemblage and compared to the 10 possible 
right femur matches. These comparisons are made using the 
model described below, with the remaining sample (total 
sample excluding the commingled individuals) acting as 
training data.

Bayesian Regression

The model used for assessing each variable is a simple lin-
ear regression, which takes the form of:

 yi =α + βxi + ε i  (1)

where yi and xi are the ith case of the dependent and indepen-
dent variables, respectively. The y- intercept is represented 
by α (alpha), and β (beta) represents the slope, or coefficient 
by which the independent variable changes in relation to the 
dependent variable. The error term is εi (sigma) and rep-
resents the stochastic part of the model that accounts for all 
other factors that influence the value of the dependent vari-
able. The y- intercept and slope are the deterministic portions 
of the model.

Typically, the regression line is fit by finding the line that 
minimizes the squared vertical distance between all data 
points. Although point estimates for the y- intercept and 
slope are calculated, uncertainty is not incorporated into 
those estimates. Confidence and prediction intervals attempt 
to deal with this limitation but are often misinterpreted and 



76 Osteometric Reassociation Using MCMC Simulation

misapplied. Linear regression of this type is associated with 
frequentist inference and does not provide an intuitive or 
easily interpretable way for comparing multiple possible 
values of yi. Bayesians specify regression models in terms of 
probability distributions, eliminating these inferential lim-
itations. Bayes’ theorem is used to specify probability dis-
tributions, taking the form of:

 ρ(θ | y,x)∝ ρ( y |θ ,x)ρ(θ ,x)  (2)

In this un- normalized form, the posterior probability ρ (θ |y,x) 
of parameter, θ, given data, y, and constant, x, is propor-
tional (for fixed y and x) to the product of the likelihood 
function ρ (y|θ,x) and prior ρ (θ,x) (Stan Development Team 
2016).

The Bayesian regression model used in this study assigns 
a normal distribution to the y- variable, with improper (uni-
form) prior distributions for regression parameters. 
Unbounded (– ∞ to +∞) uniform priors are assigned to the 
alpha and beta regression parameters, with a positive uniform 
(0< to +∞) assigned to sigma. These uniform priors are essen-
tially non- informative, leading the posterior distribution of 
the regression parameters to be driven by the likelihood of 
the training data. While on its face this model may seem sub-
optimal by assigning non- informative prior distributions to 
the regression model, on a practical level this model is needed 
because of its flexibility. Variable values will change based 
on the type of comparison and the variable values of the indi-
viduals in the training set. Thus, an abstracted regression 
model is needed to help ensure that predictions are realistic 
for all variables.

Markov Chain Monte Carlo

Bayesians view parameters as observed realizations of ran-
dom variables drawn from a probability distribution. As 
such, parameters are modeled as distributions. Modeling 
parameters as distributions requires calculus, and calculus is 
difficult, even for computers. This difficulty and the associ-
ated computational modeling time is reduced through Mar-
kov chain Monte Carlo (MCMC) simulation. MCMC 
methods provide a means for exploring the parameter space 
utilizing equation 2. Given a model, a likelihood, and data, 
MCMC simulates draws from the posterior distribution using 
quasi- dependent sequences of random variables (Kéry 2010; 
S. M. Lynch 2007). This process is repeated a large number 
of times to approximate the posterior distribution of the 
parameter, or parameter space.

Many algorithms are available for searching this param-
eter space. All of them require an initial burn- in or warmup 
period (Kéry 2010; S. M. Lynch 2007; Stan Development 
Team 2016). This period is the initial sequence of random 
draws that are strongly influenced by initial starting values 

and are not representative of the posterior distribution of the 
parameter (S. M. Lynch 2007). The Markov chain is consid-
ered representative of the posterior parameter space once the 
chain has converged to equilibrium, or entered a high prob-
ability area of the stationary distribution of the parameter 
(Stan Development Team 2016).

The effectiveness of a MCMC algorithm is measured by 
its ability to quickly reach convergence and exhaustively 
explore the parameter space. Many algorithms are inefficient 
in these respects because they can rely heavily on initial start-
ing values and incoherently search parameter space (Carpen-
ter et al. 2017). Hamiltonian Monte Carlo sampling, however, 
is both coherent and efficient (Carpenter et al. 2017). This 
method is based on modeling the behavior of particles using 
the properties of physical system (Hamiltonian) dynamics 
(Carpenter et al. 2017; Neal 2011). This system state consists of 
the position of the particle, q, and the momentum of the parti-
cle, p (Neal 2011). The position and momentum of the particle 
are described by its potential and kinetic energy, respectively 
(Neal 2011). These energy forms are inversely related. As 
this particle moves across a surface, its potential and kinetic 
energy change with the slope of the surface.

Hamiltonian dynamics are extended to searching param-
eter space by interpreting the parameter, θ, as the position of 
a fictional particle at a point in time, with a potential energy 
defined by the negative log of the probability density of θ and 
a stochastic momentum variable (Neal 2011; Stan Develop-
ment Team 2016). Stated simply, Hamiltonian MCMC is an 
efficient and effective way of exploring parameter space, 
allowing for the explicit modeling of uncertainty in parame-
ter estimates, including the dependent variable. Thus, instead 
of a point estimate for an expected bone value, Hamiltonian 
MCMC provides a distribution of values. These values are 
weighted by their relative simulated frequency. Convergence 
of the MCMC simulations is required for the simulated y- val-
ues to be a good predictive representation (Stan Develop-
ment Team 2016). Visual inspection of autocorrelation and 
chain mixture plots as well as metrics, including r- hat and 
effective sample size values, are methods for assessing model 
convergence used in this study.

The Hamiltonian MCMC sampler STAN implemented 
with the package “rstan” (Stan Development Team 2016) in 
R (R Core Team 2015) was used to simulate y- values. Spe-
cifically, each variable was modeled using 1,000 iterations 
across four chains with three simulated y- values per iteration. 
Four chains of 1,000 iterations was chosen over one chain of 
4,000 iterations for several reasons. Chains can be ran in par-
allel, or simultaneously, increasing computational efficiency 
and reducing run time. Additionally, chains have random 
starting values. The convergence and proper mixing of each 
chain provide another check of correct model behavior. The 
package “shinystan” (Stan Development Team 2016) was 
used in R (R Core Team 2015) to periodically assess model 



FIG. 2—Relationship between predictive ability of a variable and the 
distribution of simulated y- values, or the standard error around the 
mean estimate. Each sample is 100,000 random draws from a normal 
distribution with a mean of 0 and different standard deviations. The 
lighter sample (low) has a standard deviation of 1, and the darker 
sample (high) has a standard deviation of 3.

FIG. 3—The density distributions of the samples in Figure 2. A high 
standard error in the estimation of y results in low density estimates, 
especially for the mean predicted y- value. Conversely, a low standard 
error results in high density estimates for values around the mean.

McCormick 77

diagnostics to confirm proper mixing and Markov chain con-
vergence. The default in STAN is to treat the first half of 
iterations as the burn- in period (Stan Development Team 
2016). Thus, for each variable, 6,000 y- values were simu-
lated. Further treatment is required to normalize these val-
ues into a probability density function to assess the relative 
probabilities of each possible match.

Kernel Density Estimation

Kernel density estimation is a means of estimating a proba-
bility density function based on the frequency of sample 
values (Duong 2007). This family of techniques fits a con-
tinuous line to the shape of the data with a kernel and 
bandwidth. The kernel is a non- negative function centered 
on zero that integrates to one (Duong 2007). The bandwidth 
is a free parameter that determines the width of the data range 
on which the kernel function is fit. A small bandwidth for the 
data results in an under- smoothed density estimate, contain-
ing spurious data artifacts, and is essentially “connecting the 
dots” between data points. An overly wide bandwidth results 
in an over- smoothed density and obscures the underlying 
structure of these data. The bandwidth used in this study 
approaches an optimal solution for the density estimate by 
selecting a bandwidth that is the standard deviation of the 
kernel function (R Core Team 2015). The function density( ) 
in the package “stats” (R Core Team 2015) was used in R to 
fit a kernel density to the simulated y- values.

Estimating Best Matches

The result of this analysis is a probability density function 
of y- values for a given x- value for each variable on which the 
values for the 10 possible matches can be evaluated. The 
function approx( ) in the package “stats” (R Core Team 2015) 
was used in R to evaluate densities for each possible match. 
These densities are used in two ways to estimate the best 
match: density weight and equal weight. In the first best 
match estimate, each possibility is weighted by its density 
estimate for each variable. This calculation takes the form of:

 Pri =
∑(di1…n )
∑(dtot1…n )

 (3)

where Pri is the match probability for the ith possible match, 
din is the density estimate of the ith possible match for the nth 
predictive variable, and dtotn is the density estimate of all 
possible matches for the nth predictive variable. Calculating 
match probability in this way does not weigh each predictive 
variable equally. Predictive variables that have high correla-
tions between x- values and y- values will result in tightly dis-
persed simulated y- values, because uncertainty in its prediction 
is low (Fig. 2). Conversely, predictive variables that have low 

correlations also have high uncertainty in y- value predic-
tions, leading to widely dispersed y- values (Fig. 2). This 
relationship between predictive ability of a variable and the 
standard error of simulated y- values affects the resulting 
density estimates (Fig. 3). With this calculation of match 
probability, predictive variables with higher correlations will 
lead to higher density estimates and larger relative contribu-
tions to the overall match probability. However, these larger 
relative contributions may swamp the contribution of other, 
lower correlated variables, leading to spurious classifications 
if the best match from predictive variables with high cor-
relations is not the correct match.



78 Osteometric Reassociation Using MCMC Simulation

The second calculation of overall match probability 
weighs all predictive variables equally and takes the form of:

 Pri =
di1
dtot1

+
di2
dtot2

+!
din
dtotn

⎛

⎝⎜
⎞

⎠⎟ / n  (4)

where the notation is the same as formula 3. Here, densities 
are normalized into probabilities for each variable. The over-
all match probability is the sum of these probabilities divided 
by the number of variables. This way of calculating the best 
match artificially increases the relative importance of vari-
ables with low predictive ability. Each method for assess-
ing the best match has possible strengths and weaknesses. 
Thus, each type is employed to empirically address which 
performs best.

Quantile Tests

Best match probabilities are a poor metric for recognizing 
model error. Similar to other methods that classify using 
Bayesian probability, one of the possible matches will be 
classified as the best match even when the actual match is 
not among the possible choices. Thus, it is useful to have 
another metric by which to assess possible matches. To this 
end, the 5% and 95% quantiles of the simulated range of 
y- values were identified. A possible match failed this test if 
it fell outside of these boundaries. Quantile tests can be 
interpreted as a two- tailed significance test with an alpha 
level of 0.10. These tests may be used to reject possible 
matches, similar to the traditional logic, or to aid in identi-
fying model error. There is a major difference between tra-
ditional rejection- based logic, which arrives at a single 
p- value, and the quantile tests of this study (Adams & Byrd 
2006; Byrd 2008; Byrd & Adams 2003; Byrd & LeGarde 
2014; J. J. Lynch 2018; Warnke- Sommer et al. 2019). Quan-
tile tests for possible matches were conducted for each vari-
able, with the number of variables ranging from eight for the 
femur to three for many other comparisons (see Table 1). 
Quantile tests allow for the examination of this metric as a 
means of assessing model error and as a rejection criterion. 
Conducting a quantile test on each variable relates Type 1 
error rates to the number of variables rather than directly 
to the possible match. Thus, comparisons with more vari-
ables increase the number of chances for Type 1 error (reject-
ing a possible match if any variable failed the quantile test). 
The equation for the expected Type 1 error rate for correct 
matches is:

 Expected Error = 1− (1− pfail )n  (5)

where pfail is the expected chance of failing a quantile test 
and n is the number of variables.

Results

Classification Accuracy

Among all comparison types, the correct match was identi-
fied in 51.60% of the simulations (3,870/7,500). Correct clas-
sification varied by prediction and comparison types 
(Fig. 4). In all but two instances, equal- weight comparisons 
provided the best classification (Table 3 and Fig. 5). With an 
average difference of 12.12%, paired elements exhibit the 
largest difference between prediction types. Interestingly, 

FIG. 4—Accuracy by comparison and prediction type.

TABLE 3—Accuracy of Osteometric Comparison Types.

Comparison Type Equal Weight Density Weight

Femur/Femur Paired 93.00%* 91.60%
Humerus/Humerus Paired 86.20% 74.00%
Tibia/Tibia Paired 85.60% 75.80%
Radius/Radius Paired 71.60% 54.40%
Ulna/Ulna Paired 67.40% 47.40%

Paired Overall 80.76% 68.64%

Femur/Tibia Articulating 55.40% 51.80%
Ulna/Radius Articulating 42.40% 43.20%
Humerus/Radius Articulating 35.80% 34.00%
Humerus/Ulna Articulating 34.80% 32.20%

Articulating Overall 42.10% 40.30%

Femur/Humerus Other 35.60% 34.40%
Femur/Ulna Other 35.60% 32.60%
Tibia/Radius Other 34.60% 31.40%
Tibia/Ulna Other 32.80% 32.40%
Femur/Radius Other 32.60% 29.00%
Tibia/Humerus Other 30.60% 30.80%

Other Overall 33.63% 31.77%
Overall 51.60% 46.33%

*Bold indicates prediction types with the highest correct classification rates.



McCormick 79

femur/femur equal- weight comparisons only showed a 1.4% 
improvement. Articulating elements and other comparisons, 
at an average increase of 1.80% and 1.86%, respectively, 
showed a minimal difference between prediction types.

Differences in prediction type should be identifiable 
through the distribution of posterior probabilities. The dis-
tributions of density and equal- weight posterior probabilities 
are quite similar (Fig. 6). Moreover, the differences between 
median values for equal- weight and density predictions are 
negligible (0.08% for other elements, 0.13% for paired ele-
ments, and 0.23% for articulating elements). The strong 
similarities in prediction type distributions for paired ele-
ments are unexpected given the difference in accuracy. 
Given the better or similar classification and similar distribu-
tional properties, correct classification refers to equal- weight 
predictions unless otherwise specified.

Paired elements performed markedly better than the 
two other comparison types, almost doubling the accuracy 
of articulating elements. In contrast, the difference between 
articulating and other comparisons was less than 10% (see 
Table 2). Unlike prediction type, the distribution of poste-
rior probabilities by comparison type shows distinct differ-
ences (Fig. 7). The relatively low accuracies of articulating 
and other element comparisons result in slightly positively 

FIG. 5—Accuracy of all osteometric comparisons by prediction type.

skewed normal distribution. Paired elements, on the other 
hand, show a bimodal distribution with a strong positive 
skew. The shape of these distributions is in line with expec-
tations based on comparison type accuracies. The low 
accuracy of articulating and other comparisons is due to 
uncertainty. This randomness leads to posterior probabilities 
approximating a normal distribution over a large number of 
trials. The high accuracy of paired elements results in less 
uncertainty. This structure results in a model that not only 
predicts the best match well but is also good at identifying 
bad and not- so- bad matches, leading to a high density of val-
ues near zero, another peak near the median, and a long pos-
itive tail.

Besides identifying the best match, posterior probabili-
ties rank all possible matches. This aspect is most useful 
when the analyst is trying to cull down possible matches, in 
open- population situations, or in non- paired comparisons, 
where correct classification rates are relatively low. Tables 4– 6 
provide the best match rank of the correct match. For paired 
element comparisons, the correct match is among the top 
three best matches in over 97% of the simulations. For artic-
ulating and other comparisons, the correct match is among 
the top five best matches for 85.85% and 81.74% of simula-
tions, respectively.



80 Osteometric Reassociation Using MCMC Simulation

Quantile Tests

For each variable, possible match values were compared 
against the 5% and 95% quantiles of the predicted match 
distribution, resulting in 300,000 quantile tests. There are 
interesting trends in the behavior of the quantile tests (Table 7). 
Correct best match variables failed less often than the 
expected 10%, with articulating and other comparison vari-
ables failing roughly an order of magnitude less. Incorrect 
best match variables failed more often than correct match 
variables— 10.5% of the time for paired elements, but rarely 
for articulating and other comparison variables. As expected, 
variables for other possible matches failed quantile tests 
more often than best matches. Correct matches failed at 
least one quantile test (Type 1 error) less than the expected 
39.69% for paired, 30.84% for articulating, and 28.37% 
for other comparisons. Again, Type 1 error for articulating 
and other comparisons was well below expected rates (see 
Table 7).

Model Diagnostics

Figure 8 shows an example of typical model diagnostics 
plot results. These plots show that the MCMC model is 
working quite well and the parameter estimates are reli-
able. Density plots should approximate a normal distri-
bution; autocorrelation plots should look like an inverse 
exponential curve in histogram form, where autocorrela-
tion is initially high and quickly drops off. Chain mixture 
plots should show no discernible pattern, where each chain 
moves around parameter space without getting “stuck” in a 
particular area.

Metric model diagnostics were also periodically checked, 
including r- hat values and effective sample sizes. An r- hat 

FIG. 7—Distribution of posterior probabilities for other elements 
(n = 30,000 for each prediction type), articulating elements (n = 20,000 
for each prediction type), and paired elements (n = 25,000 for each 
prediction type).

FIG. 6—Distribution of posterior probabilities for (A) other elements 
(n = 30,000 for each prediction type), (B) articulating elements 
(n = 20,000 for each prediction type), and (C) paired elements 
(n = 25,000 for each prediction type).



McCormick 81

value is an estimate of convergence based on the mean and 
standard deviation estimated from each chain (Stan Devel-
opment Team 2016). Chains have properly converged with 
r- hat values between 1.0 and 1.2; the closer to 1.0, the better 

the convergence. Rarely were r- hat values above 1.0, and in 
no case was an r- hat value above 1.2. Effective sample size 
is an estimate of the information available from each simu-
lation; the closer the effective sample size is to the number 

TABLE 6—Other Element Comparison Correct Match Rank (n = 3,000 Simulations). 

Rank F/H F/U F/R T/H T/U T/R Total % Correct Cumulative %

1 178 178 163 153 164 173 1,009 33.63 33.63
2 109 87 71 96 94 99 556 18.53 52.16
3 73 53 63 84 71 48 392 13.07 65.23
4 35 42 72 69 48 54 320 10.67 75.90
5 22 36 40 36 39 27 200 6.67 82.57
6 19 21 10 23 18 19 110 3.37 85.94
7 23 24 17 21 23 21 129 4.30 90.24
8 20 16 22 7 20 21 106 3.53 93.77
9 11 13 17 7 13 19 80 2.67 96.44
10 10 30 25 4 10 19 98 3.27 100.00

F = femur, H = humerus, T = tibia, R = radius, U = ulna.

TABLE 7—Quantile Test Results for All Measurements and by Individual.

Type Position Tests Fails % Individuals Ind. Fails %

Paired Best/Correct 14,047 640 4.56 2,019 514 25.46
Best/Incorrect 2,953 310 10.50 481 48 41.16
Other Possibilities 108,000 66,349 61.43 22,500 21,686 96.38

Articulating Best/Correct 4,666 56 1.20 842 198 5.70
Best/Incorrect 6,334 107 1.69 1,158 97 8.38
Other Possibilities 63,000 19,105 30.33 18,000 13,158 73.10

Other Best/Correct 5,193 40 0.77 1,009 34 3.39
Best/Incorrect 10,307 116 1.28 1,991 104 5.21
Other Possibilities 85,500 23,639 27.65 27,000 17,452 64.64

Total 300,000 110,362 36.79 75,000 53,291 71.05

TABLE 4—Paired Element Comparison Correct Match Rank (n = 2,500 Simulations).

Rank Femur Humerus Tibia Radius Ulna Total % Correct Cumulative %

1 465 431 428 358 337 2,015 80.76 80.76
2 30 50 61 101 107 349 13.96 94.72
3 4 9 6 23 27 69 2.76 97.48
4 0 6 3 5 18 32 1.28 98.76
5 0 1 1 8 8 18 0.72 99.48
6 1 2 0 4 1 8 0.32 99.80
7 0 1 0 0 0 1 0.04 99.84
8 0 0 0 0 1 1 0.04 99.88
9 0 0 1 0 1 2 0.08 99.96
10 0 0 0 1 0 1 0.04 100.00

TABLE 5—Articulating Element Comparison Correct Match Rank (n = 2,000 Simulations).

Rank F/T U/R H/R H/U Total % Correct Cumulative %

1 277 212 179 174 842 42.10 42.10
2 102 98 98 89 387 19.35 61.45
3 50 65 80 74 269 13.45 75.90
4 17 36 36 36 125 6.25 81.15
5 15 19 30 30 94 4.70 85.85
6 15 16 14 25 70 3.50 89.35
7 7 14 16 21 58 2.90 92.25
8 6 12 18 13 49 2.45 94.70
9 5 12 15 17 49 2.45 97.15
10 6 16 14 21 57 2.85 100.00

F = femur, H = humerus, T = tibia, R = radius, U = ulna.



FIG. 8—An example of typical model diagnostics. The density plot visualizes the posterior distribution of a parameter. The autocorrelation bar  
graph represents the correlation or dependency of MCMC draws. The chain mixture or trace plot measures how well the sampler is  
exploring parameter space.

82 Osteometric Reassociation Using MCMC Simulation

of simulations, the better the chain convergence. Rarely was 
the effective sample size below 75% of the total number of 
drawn samples. Most effective sample sizes were between 
80 and 90% of the total number of draws, yet another confir-
mation that model chains are properly converging.

Discussion

The strength of a Bayesian approach to resolving com-
mingling is its versatility. The posterior distribution of 
y- values allows for the prediction of the correct match and 



McCormick 83

TABLE 8—Descriptive Statistics of the Correlation between Variables.

Comparison Type
# of 

Vars. Avg. r Max. r Min. r

Femur/Femur Paired 8 0.923 0.988 0.852
Tibia/Tibia Paired 4 0.863 0.980 0.755
Humerus/Humerus Paired 5 0.840 0.964 0.722
Radius/Radius Paired 3 0.792 0.964 0.671
Ulna/Ulna Paired 4 0.765 0.958 0.485

All Paired 24 0.837

Femur/Tibia Articulating 4 0.646 0.910 0.306
Humerus/Ulna Articulating 3 0.593 0.778 0.344
Ulna/Radius Articulating 3 0.532 0.898 0.096
Humerus/Radius Articulating 4 0.505 0.838 0.079

All Articulating 14 0.569

Tibia/Humerus Other 3 0.597 0.823 0.346
Femur/Humerus Other 4 0.574 0.866 0.307
Tibia/Ulna Other 3 0.551 0.775 0.175
Tibia/Radius Other 3 0.540 0.779 0.133
Femur/Ulna Other 3 0.508 0.778 0.155
Femur/Radius Other 3 0.501 0.809 0.114

All Others 19 0.545

rejecting possible matches. Like most practical applica-
tions in forensic anthropology, the analyst must have a 
clear question to address and a strong understanding of 
the strengths and weaknesses of the method employed. 
This study represents a start to understanding those meth-
odological aspects of a Bayesian approach to resolving 
commingling.

Equal- weight variable predictions perform better than 
density- weighted predictions. A more nuanced look at trends 
between prediction types, however, suggests underlying fac-
tors that may be affecting classification accuracy by type: 
the number, type, and predictive ability of measurements 
(Table 8). As expected, the more highly correlated variables 
used in the model, the better the accuracy. This trend may 
explain the almost nonexistent difference between prediction 
types for the femur and the large difference for other paired 
comparisons, like the radius. The correlations between left 
and right length measurements are the strongest for all ele-
ments, and are likely driving density- weighted predictions. 
For the radii, besides maximum length, there are only two 
moderately correlated measurements of the midshaft. Mid-
shaft measurements are swamped by length in density- 
weighted comparisons, but are able to adjust the best match 
to the correct match often in equal- weight predictions. In 
femur comparisons, the other strongly correlated variables 
are able to adjust predictions when length- based predictions 
are wrong, leading to comparable accuracy between the two 
types. This trend also suggests measurements that quantify 
different aspects of a bone increase model performance. A 
likely reason for the high correct classification rates of the 
femur is the novel information provided by femur measure-
ments. Stated another way, if maximum length is in the model, 

the addition of bicondylar length is unlikely to appreciably 
improve performance, as these measurements are, at least 
statistically, essentially the same (r2 = 0.995). Adding infor-
mation on epipcondylar breadth, femoral head size, and sub-
trochanteric dimensions is likely to show a marked increase 
in model performance at each step. Directly testing this 
assertion is an avenue for future research.

At an overall performance of just over 50%, it would 
appear that a Bayesian approach to osteometric reassociation 
is impractical in many situations. It is important to consider 
the difficulty of the question that these models are attempt-
ing to answer: What is the best match from among these ten 
possibilities? This question is an order of magnitude more 
difficult than the question typically asked in osteometric 
reassociation: Is this one possible match different enough that 
it can be reliably rejected as a possibility? In this light, an over-
all performance of just over 50% does not seem so bad. The 
overall correct classification rate, however, is a misleading 
metric that obscures some important aspects of osteometric 
reassociation identified in this study.

Paired element comparisons are superior to articulating 
and other comparison types. Femora, for example, were cor-
rectly matched in 93% of simulations. Paired elements are 
developmentally and (to varying degrees) functionally inte-
grated elements with directly comparable measurements. 
Composite variables are required to directly compare non- 
paired elements. While composite variables are orthogonal 
and there is a good degree of redundancy by treating each 
paired element measurement as independent, composite 
variables are likely obscuring important size and shape rela-
tionships that paired element models are able to exploit. 
This assertion is supported by the lower percentage of quan-
tile rejections for articulating and other type comparisons, 
where the composite variables may artificially make ele-
ments more homogeneous. The use of composite variables 
should have the largest negative effect in osteometric models 
based on a rejection criterion. Issues identified in rejection- 
based osteometric reassociation models (LeGarde 2012; 
McCormick 2016; Vickers et al. 2015) have been mitigated 
in recent improvements to and expansions on the paired 
element model described by Byrd (2008) (J. J. Lynch et al. 
2018; Warnke- Sommer et al. 2019). While these changes 
improve model performance, the underlying logic of the 
approach has remained the same since its description in 
Byrd and Adams (2003).

The behavior of quantile tests in this study suggests that 
deriving one composite p- value as a rejection criterion would 
be of little value outside of paired element comparisons. 
Indeed, recent research has focused exclusively on paired ele-
ment comparisons (J. J. Lynch 2018; J. J. Lynch et al. 2018; 
Warnke- Sommer et al. 2019). Rejecting a possible match, and 
by extension, identifying model error for best matches, if any 
variable failed a quantile test appears to be a viable approach 



84 Osteometric Reassociation Using MCMC Simulation

for articulating and other comparisons. This assertion does 
not apply to paired element comparisons. Over 41% of 
incorrect best matches and 96% of non- best possibilities 
fail at least one quantile test, which seems excellent for 
identifying model error and rejecting possible matches. 
However, over 25% of the correct best matches also fail at 
least one quantile test. Despite being below the expected 
Type 1 error rate (39.69%), it is an unacceptably high rate 
compared to an aggregate p- value. These results suggest the 
traditional rejection- based approach is optimal for paired 
elements.

Conclusions

This study simulated 7,500 closed- population commingled 
assemblages and assessed the accuracy of predicting the cor-
rect match using a complete set of limb measurements. The 
correct match was identified in 3,870 of the simulations, for 
an overall correct classification of 51.60%. There are several 
factors to consider when interpreting the results of this 
study. The sample used to construct the Bayesian regres-
sion model should be near- identical to the average simulated 
commingled assemblage. While the simulated commingled 
assemblage was removed from the overall sample prior to the 
construction of the model, the ten random individuals were 
drawn from the same population of predominantly European 
American males. Drawing the commingled assemblage from 
the same population as the reference sample should have two 
main influences on these results. First, the reference sample 
used to create the regression model is very appropriate for 
the commingled assemblage and represents a “best case” for 
predicting the best match. Second, the simulated commin-
gled assemblages are, on average, quite homogeneous, mak-
ing discriminating among possible matches difficult. The 
relative influence of these factors is beyond the scope of this 
study. Reference sample composition and homogeneity of 
the commingled assemblage are additional areas of future 
research. Other areas of future research include examining 
the effect of assemblage size, missing measurements, and 
different methods of quantifying skeletal elements on classi-
fication rates.

The reassociation model described above is firmly 
placed within a Bayesian framework, in both model con-
struction and inference. A Bayesian understanding of prob-
ability is easily interpreted and is in line with practical 
applications of forensic anthropology, where deductive rea-
soning is required to make statements about a particular case 
based on a larger theory of knowledge. This approach is not 
to say the frequentist paradigm is not without merit. In fact, 
this study has a major aspect that most frequentists would 
laud— the simulation of commingled assemblages to directly 
test model performance over the “long run.” The frequency 
of correct matches over an extended series of trials is an 

inductive way to build the theoretical foundation on which 
deductive statements are made. Bayesian modeling is flexi-
ble, can be tailored to various types of data, and assumptions 
can be explicitly built into the model. Modeling parameters 
as distributions provides an intuitive way to directly compare 
possible matches. The posterior distribution of y- values can 
be interpreted in different ways, depending on the goal of the 
analysis. Although rejecting possible matches has been the 
purview of the traditional, frequentist approach, there is no 
reason to limit a Bayesian approach to just prediction. The 
beauty of the model presented here is the analyst can have 
the “best of both worlds” through the ability to predict the 
best match and reject possible matches to create a short list 
of possibilities. Furthermore, Bayesian inference allows for 
incorporating additional lines of evidence into the calcula-
tion of posterior probabilities. Thus, in theory, other meth-
ods or information, such as the spatial relationship between 
elements recovered in the field, can be incorporated into an 
overall match probability.

Acknowledgments

The author would like to thank Dr. John Byrd for his men-
torship and encouragement to learn more about Bayesian 
statistics. Thanks are certainly due to Dr. Dawnie Stead-
man, Dr. Amy Mundorff, Dr. Benjamin Auerbach, Dr. Rich-
ard Jantz, and Dr. James Fordyce for their guidance and 
advice, and to the anonymous reviewers for their thoughtful 
comments.

References Cited

Adams BJ, Byrd JE, eds. Commingled Human Remains: Methods 
in Recovery, Analysis, and Identification. San Diego: Academic 
Press; 2014.

Adams BJ, Byrd JE, eds. Recovery, Analysis, and Identification of 
Commingled Human Remains. Totowa, NJ: Humana Press; 
2008.

Adams BJ, Byrd JE. Resolution of small- scale commingling: A 
case report from the Vietnam War. Forensic Science Interna-
tional 2006;156(1):63– 69.

Bookstein FL. Morphometric Tools for Landmark Data: Geome-
try and Biology. Cambridge: Cambridge University Press; 
1991.

Boulesteix A- L, Strimmer K. Partial least squares: A versatile tool 
for the analysis of high- dimensional genomic data. Briefings 
in Bioinformatics 2007;8(1):32– 44.

Brennaman AL, Love KR, Bethard JD, Pokines JT. A Bayesian 
approach to age-at-death estimation from osteoarthritis of the 
shoulder in modern North Americans. Journal of Forensic 
Sciences 2017;62(3):573– 584.

Buikstra JE, Gordon CC, St. Hoyme L. The case of the severed 
skull: Individuation in forensic anthropology. In: Rathburn 
TA, Buikstra JE, eds. Human Identification: Case Studies in 
Forensic Anthropology. Springfield: Charles C. Thomas; 
1984:121– 135.



McCormick 85

Byrd JE. Models and methods for osteometric sorting. In: Byrd 
JE, Adams BJ, eds. Recovery, Analysis, and Identification of 
Commingled Human Remains. Totowa, NJ: Humana Press; 
2008:199– 220.

Byrd JE, Adams BJ. Analysis of commingled human remains. In: 
Blau S, Ubelaker DH, eds. Handbook of Forensic Anthro-
pology and Archaeology. Walnut Creek: Left Coast Press; 
2009:174– 185

Byrd JE, Adams BJ. Osteometric sorting of commingled human 
remains. Journal of Forensic Sciences 2003;48(4):717– 724.

Byrd JE, LeGarde CB. Osteometric sorting. In: Adams BJ, Byrd JE, 
eds. Commingled Human Remains: Methods in Recovery, 
Analysis, and Identification. San Diego: Academic Press; 2014; 
167– 191.

Carpenter B, Gelman A, Hoffman MD, Lee D, Goodrich B, Betan-
court M, et al. Stan: A probabilistic programming language. 
Journal of Statistical Software 2017;76(1):1– 32. doi: 10.18637 
/jss.v076.i01

Chen Y, Hoo KA. Application of partial least square regression in 
uncertainty study area. American Control Conference 2011; 
1958– 1962.

Curtin AJ. Putting together the pieces: Reconstructing mortuary 
practices from commingled ossuary cremains. In: Schmidt 
CW, Symes SA, eds. The Analysis of Burned Human Remains. 
San Diego: Academic Press; 2008;219– 227.

Duong T. ks: Kernel density estimation and kernel discriminant 
analysis for multivariate data in R. Journal of Statistical Soft-
ware 2007;21(7):1– 16.

Gelman A. Objections to Bayesian statistics. Bayesian Analysis 
2008;3:445– 449.

Haenlein M, Kaplan AM. A beginner’s guide to partial least squares 
analysis. Understanding Statistics 2004;3(4):283– 297.

Herrmann NP, Devlin JB. Assessment of commingled human 
remains using a GIS- based approach. In: Adams BJ, Byrd JE, 
eds. Recovery, Analysis and Identification of Commingled 
Human Remains. Totowa, NJ: Humana Press; 2008;257– 269.

Hinkes MJ. The role of forensic anthropology in mass disaster res-
olution. Aviation, Space, and Environmental Medicine 1989;60: 
A60– 3.

Hoekstra R, Morey RD, Rouder JN, Wagenmakers EJ. Robust mis-
interpretation of confidence intervals. Psychonomic Bulletin & 
Review 2014;21(5):1157– 1164.

Jantz RL, Ousley SD. FORDISC 3.0: Personal computer forensic 
discriminant functions. University of Tennessee, Knoxville; 
2005.

Kéry M. Introduction to WinBUGS for Ecologists: Bayesian 
Approach to Regression, ANOVA, Mixed Models and Related 
Analyses. New York: Academic Press; 2010.

Konigsberg LW, Frankenberg SR. Bayes in biological anthropology. 
American Journal of Physical Anthropology 2013;152(S57): 
153– 184.

LeGarde CB. Asymmetry of the Humerus: The Influence of Hand-
edness on the Deltoid Tuberosity and Possible Implications for 
Osteometric Sorting [master’s thesis]. Missoula: University of 
Montana; 2012.

Lew MJ. To P or not to P: On the evidential nature of P- values 
and their place in scientific inference stat. arXiv:1311.0081; 
2013.

Lynch JJ. An analysis on the choice of alpha level in the osteomet-
ric pair- matching of the os coxa, scapula, and clavicle. Journal 
of Forensic Sciences 2018;63(3):793– 797.

Lynch JJ, Byrd J, LeGarde CB. The power of exclusion using auto-
mated osteometric sorting: Pair- matching. Journal of Forensic 
Sciences 2018;63(2):371– 380.

Lynch SM. Introduction to Applied Bayesian Statistics and Esti-
mation for Social Scientists. New York: Springer Science and 
Business Media; 2007.

Mayo DG. In defense of the Neyman- Pearson theory of confidence 
intervals. Philosophy of Science 1981;48(2):269– 280.

Mayo DG, Spanos A. Error and Inference: Recent Exchanges on 
Experimental Reasoning, Reliability, and the Objectivity and 
Rationality of Science. Cambridge: Cambridge University 
Press; 2010.

McCormick KA. A Biologically Informed Structure to Accuracy 
in Osteometric Reassociation [PhD dissertation]. Knoxville: 
University of Tennessee; 2016.

Mundorff AZ. Anthropologist- directed triage: Three distinct mass 
fatality events involving fragmentation of human remains. In: 
Adams BJ, Byrd JE, eds. Recovery, Analysis and Identifica-
tion of Commingled Human Remains. Totowa, NJ: Humana 
Press; 2008;123– 144.

Mundorff AZ. Integrating forensic anthropology into disaster vic-
tim identification. Forensic Science, Medicine, and Pathology 
2012;8:131– 139.

Neal RM. MCMC using Hamiltonian dynamics. In: Brooks S, Gel-
man A, Jones GL, Meng X- L, eds. Handbook of Markov Chain 
Monte Carlo. Boca Raton, FL: CRC Press: 2011: 113– 162.

O’Brien MJ, Storlie CB. An alternative bilateral refitting model for 
zooarchaeological assemblages. Journal of Taphonomy 2011;9: 
245– 268.

Primorac D, Andelinovic S, Definis- Gojanovic M, Drmic I, Rezic 
B, Baden MM, et al. Identification of war victims from mass 
graves in Croatia, Bosnia, and Herzegovina by the use of stan-
dard forensic methods and DNA typing. Journal of Forensic 
Sciences 1996;41(5):891– 894.

R Core Team. R: A language and environment for statistical comput-
ing. R Foundation for Statistical Computing, Vienna, Austria, 
2015. http:  //www  .R  - project  .org  /.

Rosing FW, Pischtschan E. Re- individualisation of commingled 
skeletal remains. In: Jacob B, Bonte W, eds. Advances in 
Forensic Sciences. Berlin: Verlag Dr. Köster; 1995:1– 9.

Rosipal R, Krämer N. Overview and recent advances in partial least 
squares. In: Saunders C, Grobelnik M, Gunn S, Shawe- Taylor 
J, eds. Subspace, Latent Structure and Feature Selection. Ber-
lin: Springer; 2006:34– 51.

Royall R. On the probability of observing misleading statistical evi-
dence. Journal of the American Statistical Association 2000; 
95(451):760– 768.

Royall R. Statistical Evidence: A Likelihood Paradigm. Boca 
Raton, FL:CRC Press; 1997. Monographs on Statistics and 
Applied Probability; vol. 71.

Sanchez G. plsdepot: Partial least squares (PLS) data analysis meth-
ods. R package version 0.1, 17; 2016.

Sledzik PS, Rodriguez WC. Damnum fatale: The taphonomic fate 
of human remains in mass disasters. In: Haglund WD, Sorg 
MH, eds. Advances in Forensic Taphonomy: Method, Theory, 
and Archaeological Perspectives. Boca Raton, FL: CRC Press; 
2001;322– 330.

Snow CC, Folk ED. Statistical assessment of commingled skeletal 
remains. American Journal of Physical Anthropology 1970;32: 
423– 427.

Stan Development Team. Stan modeling language: User’s guide and 
reference manual, Version 2.17.0; 2016.

Stark PB, Freedman DA. What is the chance of an earthquake? 
NATO Science Series IV: Earth and Environmental Sciences 
2003;32:201– 213.

Steadman DW, Sperry K, Snow F, Fulginiti L, Craig E. Anthropo-
logical investigations of the Tri- State Crematorium incident. 
In: Adams BJ, Byrd JE, eds. Recovery, Analysis and Identifi-
cation of Commingled Human Remains. Totowa, NJ: Humana 
Press; 2008:81– 96.

Ubelaker DH, Rife JL. Approaches to commingling issues in arche-
ological samples: A case study from Roman era tombs in 
Greece. In: Adams BJ, Byrd JE, eds. Recovery, Analysis and 

http://www.R-project.org/


86 Osteometric Reassociation Using MCMC Simulation

Identification of Commingled Human Remains. Totowa, NJ: 
Humana Press; 2008:97– 122.

Varas CG, Leiva MI. Managing commingled remains from mass 
graves: Considerations, implications and recommendations 
from a human rights case in Chile. Forensic Science Interna-
tional 2012;219:e19– e24.

Vickers S, Lubinski PM, DeLeon LH, Bowen JT. Proposed method 
for predicting pair matching of skeletal elements allows too 
many false rejections. Journal of Forensic Sciences 2015;60(1): 
102– 106.

Warnke- Sommer JD, Lynch JJ, Pawaskar SS, Damann FE. 
Z- Transform method for pairwise osteometric pair- matching. 
Journal of Forensic Sciences. 2019;64(1):23–33.

Wegelin JA. A survey of partia1 least squares (PLS) methods, 
with emphasis on the two- block case. Technical report, 
Department of Statistics, University of Washington, Seattle; 
2000.

Willey PS. Prehistoric Warfare on the Great Plains: Skeletal Anal-
ysis of the Crow Creek Massacre Victims. New York: Garland 
Publications; 1990.


