7
$\begingroup$

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.

$\endgroup$
  • $\begingroup$ 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. $\endgroup$ – probabilityislogic Apr 22 '14 at 4:51
  • $\begingroup$ 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? $\endgroup$ – Sameer Apr 22 '14 at 18:16
  • $\begingroup$ 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. $\endgroup$ – probabilityislogic Apr 24 '14 at 9:41
  • 1
    $\begingroup$ If it really only is an issue of sampling speed, then Stan might be an alternative. $\endgroup$ – coffeinjunky Apr 24 '14 at 22:04
  • 1
    $\begingroup$ 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... $\endgroup$ – coffeinjunky Apr 25 '14 at 1:08
2
+50
$\begingroup$

Have you looked into the stats literature on reduced-rank vector generalized linear models? It seems you might be trying to identify the model by introducing constraints that in effect make it a hybrid of the probit and ordinal regressions, plus with a multilevel structure on the right-hand side of the equation. This document (pdf alert) on the VGAM package in R offers a relatively straightforward introduction to reduced-rank generalized models.

$\endgroup$
  • $\begingroup$ Thanks for the recommendation. While it looks like rrvglm() is along the lines of what I need, VGAM only implements a bivariate logit/probit model and it's not clear to me how I could use constraints to coax a multinomial logit model into a multivariate logit. $\endgroup$ – Sameer Apr 21 '14 at 15:05

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.