[R] multivariate ordinal probit regression vglm()

[Trey Batey]

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 at mhcc.edu or ekt.batey at 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     ref$age:1     ref$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     ref$age:1     ref$age:2
  -2.07719235   -2.43422370   -2.99123098    0.07319632    0.05133132
    ref$age:3
   0.03797696

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