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:



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:

x=c(5.6, 5.1)

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


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


2 Answers 2


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
  • $\begingroup$ 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. $\endgroup$
    – Ho1
    Sep 26, 2015 at 5:03
  • $\begingroup$ There is a fundamental problem: How can I calculate prediction accuracy on ordered categorical data? stats.stackexchange.com/questions/174255 $\endgroup$
    – Ho1
    Sep 26, 2015 at 10:51
  • 1
    $\begingroup$ That's what the answer above was trying to accomplish. $\endgroup$ Sep 27, 2015 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):


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


To get the model errors.

  • $\begingroup$ I get NULL when I invoke (residuals(reg)). $\endgroup$
    – Ho1
    Sep 26, 2015 at 5:09
  • $\begingroup$ 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). $\endgroup$
    – Nester
    Sep 26, 2015 at 22:15

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge that you have read and understand our privacy policy and code of conduct.

Not the answer you're looking for? Browse other questions tagged or ask your own question.