unrestricted cumulative probit model with vglm()---extend to a multivariate setting?

Hi, all. I'm trying to model the transition of multiple, ordinal response variables on a continuous predictor (age). In this specific case, I'm investigating the ages at which a skeletal trait passes through a series of morphological stages. In anthropological parlance, this has become known as "transition analysis." Imagine that you have two datasets. The first, or "reference" dataset, comes from a museum collection of skeletons with known age-at-death and consists of three variables: age, pube3, and auric4 (if interested, feel free to access the reference dataset). The second, or "target" dataset, comes from a sample of skeletons whose age distribution you'd like to estimate and has two variables: tpube3 and tauric4 (the observed stages for the two skeletal traits; if interested, access the target dataset).

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 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
vglm(formula = pube3 ~ age, family = cumulative(link = "probit",
parallel = FALSE, reverse = TRUE), data = refdata)

(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
vglm(formula = auric4 ~ age, family = cumulative(link = "probit",
parallel = FALSE, reverse = TRUE), data = refdata)

(Intercept):1 (Intercept):2 (Intercept):3 age:1 age:2
-2.07719235 -2.43422370 -2.99123098 0.07319632 0.05133132

Degrees of Freedom: 2238 Total; 2232 Residual
Residual Deviance: 1583.47
Log-likelihood: -791.7348