Calculating predicted probabilities for ordinal logistic regression I've run an ordered logistic regression model in R with Zelig and am looking to calculate predicted probabilities. Zelig has a series of simple one line commands to generate the information I want on first differences and so forth. Unfortunately, I keep getting an error when running the zelig function and was wondering if there was a quick alternative for generating predicted probabilities for a ordered logit in R.
For what it's worth, here's the error from my Zelig code.
> x.out <- setx(mod, credit=1)
Error in dta[complete.cases(mf), names(dta) %in% vars, drop = FALSE] : 
  incorrect number of dimensions

I just need an alternative solution that I can use to generate the probabilities.
 A: A similar issue was raised on Stack Overflow more than one year ago. I don't know if re-installing Zelig and its dependencies will solve your problem (especially because I would prefer to understand why this error message came up before reinstalling).
Anyway, you can use the lrm() function from the rms package, as it allows to fit several models for categorical outcomes including proportional odds model. There is a predict() (but also Predict()) function to get the desired predicted values. As an alternative, you may want to look at the ordinal package (see the clm() function).
A: There is also polr in MASS that can fit the proportional odds cumulative logit model, which I like because you can show the fitted model easily using the effects package (for lrm and clm this is not the case I believe, except maybe with Effect() in the development version), and also get 95% confidence limits on your predictions. Here is a little example :
data=read.table("https://onlinecourses.science.psu.edu/stat504/sites/onlinecourses.science.psu.edu.stat504/files/lesson07/cheese.dat",col.names=c("cheese","response","count"))
data$response=factor(data$response, ordered=T)
fit=polr(response ~ cheese, weights=count, data=data)
summary(fit)
z=summary(fit)$coefficients/summary(fit)$standard.errors
p=(1 - pnorm(abs(z), 0, 1)) * 2 # 2-tailed Wald z tests to test significance of coefficients
p
Anova(fit, type="III")  # likelihood ratio tests for effect of different factors
library(colorRamps)
plot(allEffects(fit),ylab="Response",rescale.axis=F,style="stacked",colors=colorRampPalette(c("white","red"))(9))


