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

#1
Hello. I'm trying to fit a multivariate, unrestricted cumulative probit model. I've had success in fitting a proportional odds probit model using polr(method="probit") in the MASS package and the unrestricted cumulative probit model using vglm(cumulative(link="probit)) in the VGAM package. I've only been able to do so with one response variable, though, but I'd really like to include two response variables in one model because they are not independent. Below, I've included the results of the two univariate cases using vglm(). If interested, you can access the data as a *.csv file. Also, the attached file describes the mathematical approach that I'm following as it relates to the data. Any suggestions are welcome. Thanks in advance.

--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="probi t",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="prob it",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