R help - Apr 2012 - multivariate ordered probit regression---use standard bivariate normal distribution?

Hello.

I have yet to receive a response to my previous post, so I may have
done a poor job asking the question. So, here is the general question:
how can I run a run a multivariate (more than one non-independent,
response variables) ordered probit regression model? I've had success
doing this in the univariate case using the vglm() function in the
VGAM package. For example:

unvar.prob<-vglm(y~x,cumulative(link="probit",parallel=FALSE,r
everse=TRUE))

This would fit a unrestricted cumulative probit model (or a
"thresholds of change" model) with the ordered response variable
treated as a latent trait modeled with a standard normal distribution.
What I'd like to do is include an additional, ordinal response
variable, so that the two (non-independent) response variables are
treated as latent traits and modeled with a standard bivariate normal
distribution.

The VGAM package includes the function binormal(), which seems like it
should be of some use here, but if so, I don't know where/how it fits
into the model. I've included the results of running the model with
the two response variables separately. Any suggestions are welcome.
Again, if it helps in your response, I've included a link to the data
below.  Thanks.

--Trey


Get the data (a *.csv file) here:
https://docs.google.com/open?id=0B5zZGW2utJN0TEctcW1oblFRcTJrNDVLOVBmRWRaQQ


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

##########################################################

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


*********************************************
Trey Batey--Anthropology Instructor
Division of Social Sciences
Mt. Hood Community College
Gresham, OR 97030