# How to predict using ordered probit regression and calculate prediction accuracy?

I want to do an ordered probit regression, then cross-validate model prediction accuracy with 80% data for training and 20% for validation, and calculate RMSE for predictions.

Consider this dataset:

  X     Y
----------
2.3    1
3.1    2
3.5    2
10.0   5
6.8    4
5.0    3
5.4    2
3.2    1


I did this:

x=c(2.3,3.1,3.5,10.0,6.8,5.0,5.4,3.2)
y=c(1,2,2,5,4,3,2,1)
myData=data.frame(cbind(x,y))

library("MASS")
reg=polr(as.factor(myData$y)~myData$x,data=myData,method="probit")


I saw this question, but I couldn't fully understand. Suppose myValidationData contains 20% of data which I want to use for validation. So, I would do:

fit=predict(reg,type="probs")
x=c(5.6, 5.1)
y=c(3,3)
myValidationData=data.frame(cbind(x,y))


This is how I tried to predict, but is it correct, when I want to cross-validate?

fit=predict(reg,data=myValidationData,type="probs")


How should I measure RMSE? And, how can I plot the prediction?

The R rms package has many capabilities for validating ordinal regression models. Start with the orm function. Note that split-sample validation takes an extremely large sample size to work. You might be better off with bootstrap validate as implemented in the rms validate and calibrate functions.

Measures of predictive accuracy for ordinal $Y$ include

• Generalized $c$-index (generalized ROC area) from Somers' $D_{xy}$ rank correlation
• Spearman $\rho$
• Other rank correlation measures - these are all measures of pure predictive discrimination
• Generalized $R^2$ based on model likelihood ratio $\chi^2$ statistic
• Calibration accuracy for $Prob(Y \geq y | X)$ using a nonparametric smooth calibration curve
• Thanks for the answer. My actual data is about ~2000 elements for training and ~200 elements for validation, so I think it would be fine. – Ho1 Sep 26 '15 at 5:03
• There is a fundamental problem: How can I calculate prediction accuracy on ordered categorical data? stats.stackexchange.com/questions/174255 – Ho1 Sep 26 '15 at 10:51
• That's what the answer above was trying to accomplish. – Frank Harrell Sep 27 '15 at 15:54

That looks correct; note that after you make the myData data frame, you can use:

myData$y = factor(myData$y)
reg <- polr(y ~ x, data = myData, method = "probit")


Later, you can make the validation data with:

myValidationData <- data.frame(x = c(5.6, 5.1), y = c(3,3))


Your syntax works fine, I just thought this was a bit "cleaner".

Here's a great link on ordinal regression that includes syntax to plot the predicted values (it's actually about SEM with categorical variables, but there's a great section on ordinal regression in the middle):

http://www.personality-project.org/r/tutorials/summerschool.14/rosseel_sem_cat.pdf

I'm not sure how to measure the RMSE, but you should be able to use

(residuals(reg))


To get the model errors.

• I get NULL when I invoke (residuals(reg)). – Ho1 Sep 26 '15 at 5:09
• Hmmm, it looks like polr() doesn't supply a residuals object. In that case, try predict(reg) - myData\$y to get the residuals (without a newdata argument, predict() uses the dataset used to fit the model). – Nester Sep 26 '15 at 22:15