multivariate ordinal probit regression with vglm()

#1
Hello, all.

I'm investigating the rate at which skeletal joint surfaces pass through a series of ordered stages (changes in morphology). Current statistical methods in this type of research use various logit or probit regression techniques (e.g., proportional odds logit/probit, forward/backward continuation ratio, or restricted/unrestricted
cumulative probit). Data typically include the predictor (age) and one or more response variables (the stages of skeletal morphology). For example, the pubic symphysis and auriclar surface (two joint surfaces of the pelvis) may be observed in three and four stages, respectively (see sample dataframe "refdata" below).

age pube3 auric4
1 32 3 2
2 42 3 2
3 27 2 1
4 39 2 1
5 85 3 4

I've had some success in fitting the ordinal probit model using both polr(method="probit") in the MASS library and vglm() in the VGAM library, but I've hit a wall when it comes to fitting a model that includes both response variables ("pube3" and "auric4" in the sample dataframe above). I've included the two univariate models below, but I'd like to know how to model the two response variables on age---returning the coefficients for each response AND the correlation parameter, since the two responses are not independent. If it would help, please feel free to access the dataframe (https://docs.google.com/open?id=0B5zZGW2utJN0TEctcW1oblFRcTJrNDVLOVBmRWRaQQ).
Thanks in advance.

--Trey

***************************
Trey Batey---Anthropology Instructor
Mt. Hood Community College
26000 SE Stark St.
Gresham, OR 97030
trey.batey@mhcc.edu or ekt.batey@gmail.com


## unrestricted cumulative probit model for pubic scores
> mod.pube3
Call:
vglm(formula = refdata$pube3 ~ refdata$age, family = cumulative(link =
"probit",
parallel = FALSE, reverse = TRUE))

Coefficients:
(Intercept):1 (Intercept):2 refdata$age:1 refdata$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

## unrestricted cumulative probit model for auricular scores
> mod.auric4
Call:
vglm(formula = refdata$auric4 ~ refdata$age, family = cumulative(link
= "probit",
parallel = FALSE, reverse = TRUE))

Coefficients:
(Intercept):1 (Intercept):2 (Intercept):3 refdata$age:1 refdata$age:2
-2.07719235 -2.43422370 -2.99123098 0.07319632 0.05133132
refdata$age:3
0.03797696

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