# Inferring parameters for a regression with features of both multivariate probit and ordinal regression?

Based on my data generated by a complex process and the problem below detailed, I have tried various approaches, to no avail.

I am trying to answer one or more of the following questions:

a) Has someone solved my problem already? Are there other clever approaches I could use instead (perhaps something like PLS or CCA)?

b) What is the right literature to look for previous work addressing my problem?

c) If this problem is novel (which I personally doubt), what is the right approach to solve it?

d) In my proposed approach below described, is there something obvious I am missing?

Here is my problem:

For $$N$$ individuals, with $$M$$ associated covariates, I observe outcomes in $$J$$ possible classes; thus my design matrix $$X$$ is of dimension $$M$$ x $$N$$ and my response matrix $$Y$$ is of dimension $$N \times J$$. For a given class $$j$$ there is a latent variable $$U_{i,j}$$ s.t. when $$U_{i,j} < 0$$, $$Y_{i,j}=0$$ and $$Y_{i,j}=1$$, otherwise . Up to this point, if $$U \sim MVN(\beta X, \Sigma)$$, I have described a standard multivariate probit regression scenario. However, if a given $$U_{i,j}$$ is greater than all other classes $$U_{i,-j}$$ and greater than some constant $$c_1$$, then $$Y_{i,j} =t$$ if $$c_{t-1} < U_{i,j} < c_t$$ where $$t\in \{1,2,...,T-1\}$$; $$Y_{i,j}=T$$ if $$U_{i,j} > c_{T-1}$$. Here the cutpoints $$c_t$$ are unknown. On its own, the latter half corresponds to an ordinal regression problem. However, together this problem seems to be more complex than either probit or ordinal regression.

Below is the R code to generate data from this process.

The values of $$U$$ are generated directly without specifying a linear relationship with $$X$$:

library(MASS) # for mvrnorm

N <- 1e3 # num. of observations
J <- 12  # num. of classes
cuts <- c(1.5, 2, 2.25, 2.5, 2.7, 3) # T=6

S <- rWishart(n=1, Sigma=diag(J), df=20) # sample cov. mat.
U <- mvrnorm(n=N, mu=rep(0, J), Sigma=S[,,1]) # generate U
Y <- apply(U, 1, function(u) { # generate Y from U
u <- as.numeric(scale(u)) # s.d. is 1 for identifiability
y <- ifelse(u < 0, 0, 1)  # U's extreme enough are observed
cond <- u > cuts[1] & u == max(u) # most extreme U's are ranked
if(any(cond))
y[which(cond)] <- findInterval(u[which(cond)], cuts)
return(y)
}
)


In my case, $$N$$ is approximately $$3500$$, $$J=12$$ and $$M=30$$, and $$T=6$$. The outcomes of the classes are correlated and the individuals are hierarchically related.

Here are the solutions I have tried so far:

Initially, I thought I could use JAGS to solve this problem. Before implementing the full model, I tried a hierarchical multivariate probit regression. However, JAGS seemed to take an inordinately long time to sample from the posterior.

Thus, I have since been searching through the simulated likelihood literature e.g. GHK, but found no work reproducing the model described above. Specifically, I am exploring approaches such as approximate Bayesian computation (ABC), indirect inference, and simulated method of moments.

Any help would therefore be much appreciated.

• Do you really need the notion of a latent class? It is very general and may be confusing you. Also GLMMs are notorious for having bad mixing - it really depends on how you set up the sampler. BayesX and INLA may be better tools than jags. – probabilityislogic Apr 22 '14 at 4:51
• What is the alternative to a latent class when dealing with ordinal data? Wouldn't I lose power by not modeling the hidden process underlying the observed data? – Sameer Apr 22 '14 at 18:16
• Only example is to model $Pr (Y> y)$. and make this an increasing/decreasing function of $y$. I think this is what is done in hazard/survival modelling. – probabilityislogic Apr 24 '14 at 9:41
• If it really only is an issue of sampling speed, then Stan might be an alternative. – coffeinjunky Apr 24 '14 at 22:04
• By the way, just to make sure I understand your problem correctly, let me phrase it in my own words: you want to estimate an ordered probit, and the complication is that you have several outcomes, so that the ordered probit becomes a multivariate ordered probit. Does this sum up your problem? I am asking because it sounds like that to me, but it seems that there are papers out there dealing with multivariate ordered probits, so I conclude I must misinterpret sth... – coffeinjunky Apr 25 '14 at 1:08