The basic idea is to use the age-to-transition information from the reference dataset, the observed distribution of skeletal traits in the target dataset, and a parametric mortality model (e.g., Gompertz, Gompertz-Makeham...) to estimate the age distribution of the target group. So, one uses an MLE approach to find the mortality model parameters that most likely produced the observed distribution of skeletal traits in the target sample.

An example of the procedure from a workshop in 2007 is provided by Lyle Konigsberg. In the BARFAA workshop, Konigsberg used the polr() function from the MASS package, but more recently provided some examples of how to run the unrestricted cumulative probit model using the vglm() function in the VGAM package.

Using refdata as the reference dataset, I've been able to model the ages-to-transition for each skeletal trait separately, but I'd like to run the model with both traits (pube3 and auric4) at the same time, since they are not independent. Herrmann and Konigsberg (2002) model the multivariate case (see attachment), but I’m not sure how to actually accomplish this in R. Please feel free to work through the analysis with the datasets provided. I’ll probably be posting additional questions that apply to estimating the age-at-death distribution of the target sample (targdata). Any suggestions are welcome. Thanks.

--Trey

*******************

Trey Batey¬¬—Anthropology Instructor

Mt. Hood Community College

26000 SE Stark St.

Gresham, OR 97030

**********************************************

> mod.pube3<-vglm(pube3~age,data=refdata,cumulative(link="probit",parallel=FALSE,reverse=TRUE))

> mod.pube3

Call:

vglm(formula = pube3 ~ age, family = cumulative(link = "probit",

parallel = FALSE, reverse = TRUE), data = refdata)

Coefficients:

(Intercept):1 (Intercept):2 age:1 age:2

-1.65895567 -2.14755951 0.06688242 0.04055919

Degrees of Freedom: 1492 Total; 1488 Residual

Residual Deviance: 1188.909

Log-likelihood: -594.4543

> mod.auric4<-vglm(auric4~age,data=refdata,cumulative(link="probit",parallel=FALSE,reverse=TRUE))

> mod.auric4

Call:

vglm(formula = auric4 ~ age, family = cumulative(link = "probit",

parallel = FALSE, reverse = TRUE), data = refdata)

Coefficients:

(Intercept):1 (Intercept):2 (Intercept):3 age:1 age:2

-2.07719235 -2.43422370 -2.99123098 0.07319632 0.05133132

age:3

0.03797696

Degrees of Freedom: 2238 Total; 2232 Residual

Residual Deviance: 1583.47

Log-likelihood: -791.7348

***********************************************