This part, a vanilla probit, works perfectly -- # Simulate from probit model -- x1 <- 2*runif(5000) x2 <- 5*runif(5000) ystar <- 7 + 3*x1 - 4*x2 + rnorm(5000) y <- as.numeric(ystar>0) table(y) # Estimation using micEcon::probit() library(micEcon) summary(probit(y ~ x1 + x2)) # Estimation using glm() -- summary(glm(y~x1+x2, family=binomial(probit)) ) But this part, where I move on to ordered probit, gives me wonky values -- # Simulate from ordered probit model y <- cut(ystar, breaks=c(-100, -5, 0, 5, 100)) table(y) # Estimate ordered probit model -- library(MASS) summary(polr(y ~ x1 + x2, method="probit")) I get coefficients which are quite far from c(7,3,-4) which were used in simulation above and I get breaks which are quite far from c(-5,0,5) which were used in simulation above. I wonder what I'm doing wrong. -- Ajay Shah http://www.mayin.org/ajayshah ajayshah at mayin.org http://ajayshahblog.blogspot.com <*(:-? - wizard who doesn't know the answer.