Simulating ordinal variables using fitted probit models I have fitted a probit model for an ordinal response and a number of predictors, using polr function in R. Now I want to use this fitted model in order to sample from the conditional distribution of the ordinal variable, given a set of specific values for the predictors. I was thinking of using multinomial distribution, with the probabilities computed based on the coefficient estimates from the fitted probit model. But then I am thinking that this approach ignores any lack of fit for the model. I am not sure how this simulation could be performed. Perhaps I am missing or confusing some key concepts. Any help would be appreciated.
 A: If the name of your model object is fit, predict(fit, type = 'class') will return the predicted class for the data you fit your model to. If you simulate data, (however you like), you can pass this to the predict function via the newdata argument, and get predicted categories or probabilities, whichever you prefer. expand.grid is great for doing that also. Something like predict(fit, type = 'class', newdata = sim.df will work.
What this is calculating is $Pr(Y_i = j | X_i \hat{\beta})$
We are essentially cutting up a latent continuous variable $Y^*$ with $J-1$ cutpoints (denoted $\tau$ with $\tau_1 = -\infty$ and $\tau_{j-1} = \infty$. Thus
$$Pr(\tau_{j-1} \leq Y^*_i \leq \tau_j | X)$$
$$ = Pr(\tau_{j-1} \leq X_i\beta + u_i \leq \tau_j)$$
$$ = Pr(\tau_{j-1} - X_i\beta \leq u_i \leq \tau_j - X_i \beta)$$
$$ = \int^{\tau_{j} - X_i \beta}_{-\infty} f(u_i) du - \int^{\tau_{j-1} - X_i \beta}_{-\infty} f(u_i)du$$
$$ = F(\tau_j - X_i \beta) - F(\tau_{j-1} - X_i \beta)$$
Where $f$ is the density for $u$ and $F$ is the cdf.
For the probit this just means that (for 3 categories)
$$Pr(Y_i = 1) = \phi(\tau_1 - X_i \beta) $$
$$Pr(Y_i = 2) = \phi(\tau_2 - X_i \beta) - \phi(\tau_1 - X_i \beta)$$
$$Pr(Y_i = 3) = 1 - \phi(\tau_3 - X_i \beta)$$
When you call predict with the class argument it returns the category with the highest predicted probability.
A: For simulation you use the predicted probabilities of belonging to each category. You draw a single column vector from a standard continuous uniform distribution, and assign a an observation to the first category if the random number is less than the first predicted probability, to the second catogry if it is larger than the first predicted probability but less than the sum of the first and second predicted probability, etc. I am not an R user, so it makes no sense for me to provide you R code, but in Stata that would look like this:
webuse fullauto, clear
oprobit rep77 foreign length mpg

// get predicted probabilities
predict double pr1 pr2 pr3 pr4 pr5

// turn them into cumulative probabilities
gen pr0 = 0 
replace pr2 = pr1 + pr2
replace pr3 = pr2 + pr3
replace pr4 = pr3 + pr4
replace pr5 = 1

// draw a random number from a standard continuous uniform distribution
gen u = runiform()

// create a new simulated variable
gen byte sim = cond(u > pr0 & u < pr1, 1, ///
               cond(u > pr1 & u < pr2, 2, ///
               cond(u > pr2 & u < pr3, 3, ///
               cond(u > pr3 & u < pr4, 4, 5 ))))

// check it out
oprobit sim foreign length mpg

