A few minutes ago I had asked why this didn't seem to work: # Simulate from probit model -- x1 <- 2*runif(5000) x2 <- 5*runif(5000) ystar <- 7 + 3*x1 - 4*x2 + rnorm(5000) y <- cut(ystar, breaks=c(-100, -5, 0, 5, 100)) table(y) library(MASS) summary(polr(y ~ x1 + x2, method="probit")) A little thinking revealed why. In the ordered probit model, the intercept "7" is not identified. So the estimation recovers -7 + c(-5,0,5) = c(-12,-7,-2) For the rest, it works fine. The slopes c(3,-4) are recovered correctly. :-) -- Ajay Shah http://www.mayin.org/ajayshah ajayshah at mayin.org http://ajayshahblog.blogspot.com <*(:-? - wizard who doesn't know the answer.