Cross Validated is a question and answer site for people interested in statistics, machine learning, data analysis, data mining, and data visualization. Join them; it only takes a minute:

Sign up
Here's how it works:
  1. Anybody can ask a question
  2. Anybody can answer
  3. The best answers are voted up and rise to the top

I want to predict a health problem. I have 3 outcome categories that are ordered: 'normal', 'mild', and 'severe'. I wish to predict this from two predictor variables, a test result (a continuous, interval covariate) and family history with this problem (yes or no). In my sample, the probabilities are 55% (normal), 35% (mild), and 10% (severe). In this sense, I could always just predict 'normal' and be right 55% of the time, although this would give me no information about individual patients. I fit the following model:

\begin{align} \text{the cut point for }\widehat{(y \ge 1)} &= -2.18 \\ \text{the cut point for }\widehat{(y \ge 2)} &= -4.27 \\ \hat\beta_{\rm test} &= 0.60 \\ \hat\beta_{\rm family\ history} &= 1.05 \end{align}

Assume there is no interaction and everything is fine with the model. The concordance, c, is 60.5%, which I understand to be the maximum predictive accuracy the model affords.

I come across two new patients with the following data: 1. test = 3.26, family = 0; 2. test = 2.85, family = 1. I want to predict their prognosis. Using the formula: $$ \frac{\exp(-X\beta - {\rm cutPoint})}{(1+\exp(-X\beta - {\rm cutPoint}))} $$ (and then taking the differences amongst the cumulative probabilities), I can calculate the probability distribution over the response categories conditional on the model. R code (n.b., due to rounding issues, the output doesn't match perfectly):

cut1 <- -2.18
cut2 <- -4.27
beta <- c(0.6, 1.05)
X    <- rbind(c(3.26, 0), c(2.85, 1))

pred_cat1      <- exp(-1*(X%*%beta)-cut1)/(1+exp(-1*(X%*%beta)-cut1))
pred_cat2.temp <- exp(-1*(X%*%beta)-cut2)/(1+exp(-1*(X%*%beta)-cut2))
pred_cat3      <- 1-pred_cat2.temp
pred_cat2      <- pred_cat2.temp-pred_cat1

predicted_distribution <- cbind(pred_cat1, pred_cat2, pred_cat3)

Namely: 1. 0 = 55.1%, 1 = 35.8%, 2 = 9.1%; and 2. 0 = 35.6%, 1 = 46.2%, 2 = 18.2%. My question is, how do I go from the probability distribution to a predicted response category?

I have tried several possibilities using the sample data, where the outcome is known. If I just pick max(probabilities), accuracy is 57%, a slight improvement over the null, but below the concordance. Moreover, in the sample, this approach never picks 'severe', which is what I really want to know. I tried a Bayesian approach by converting null and model probabilities into odds and then picking the max(odds ratio). This does pick 'severe' occasionally, but yields worse accuracy 49.5%. I also tried a sum of the categories weighted by the probabilities and rounding. This, again, never picks 'severe', and has low accuracy 51.5%.

What is the equation that takes the information above and yields optimal accuracy (60.5%)?

share|improve this question
The description of your "models" make no sense. Please edit to make this more clear. – Andy W Nov 9 '11 at 13:52
@AndyW I have only 1 model here, but would like it to be intelligible. If you can tell me what is confusing or missing, I will edit it. This is an ordinal logistic regression model, not an OLS reg model. The only thing I can think of is that I left out the exponentiation. Which I am now fixing. – gung Nov 9 '11 at 15:22
After the line I fit the following model in the first paragraph, what is y>=1-2.18 referring to? What is test .60 referring to (a value of .60 for the attribute test for one observation)? I also have no idea what the formula exp(-Xb -cutPoint) means. Is this the syntax for some software routine? Regardless, -Xb and -cutPoint are undefined so it is difficult to know what the operation does. – Andy W Nov 9 '11 at 15:35
@AndyW the former are model parameter estimates, the latter is a mathematical formula. I think part of the problem was CV was eliminating my whitespaces. I have edited it to make it as explicit as possible, and have added some R code that performs the computation. – gung Nov 9 '11 at 17:11
up vote 8 down vote accepted

You are making a leap that you need to classify predicted values. The fact that your method never picks the "severe" category is a consequence of the discrete nature of the problem and that "severe" is infrequent. With ordinal response models you can just use exceedance probabilities on their own (for all but one category) or just quote the individual probabilities. If $Y$ is roughly interval scaled you can also use the predicted mean. These are all available in the R rms package lrm and associated function predict.lrm. Many people assume that classification is the goal when in fact risk prediction is the underlying goal.

share|improve this answer
Thanks for your help. I suspected the low freq of severe was part of the issue. I think my raw Y, 0 1 2, is insufficiently equal interval. I gather my goal is mistaken. Unfortunately, I think I want to know what category a new patient will fall into / don't fully understand what my goal should be. Is it possible to provide a little more insight? (Actually, I suspect CV isn't the forum for a complete lesson; alternatively, do you know where I could learn about this issue? I have read the sections from Agresti's Intro & Hosmer & Lemeshow's Logistic, but to no avail.) – gung Nov 9 '11 at 15:16
The goal is driven by the needed decision or by the subject matter. If you state the ultimate goal I might be able to comment. – Frank Harrell Nov 9 '11 at 23:11
Sorry about my lack of clarity, it seems to have been a problem above also. At present, I want to be able to predict the outcome for new cases. In the long run, I want to understand ord log reg better, e.g. how do you get residuals if you don't have predicted categories? I gather slightly better accuracy is possible, but I don't know how to get it. I'm sure you don't have time to explain everything, but neither Agresti, nor H&L say anything about predictions or residuals, etc. and I couldn't find anything w/ google. Thus, I asked on CV. I do appreciate your continued assistance. – gung Nov 10 '11 at 2:50
Please explain your need for residuals. These are not in intrinsic part of the model. In terms of how to state predictions, providing the predicted probability that $Y\geq j$ for various $j$ is a good way, and requires no arbitrary decisions. Looking at the histogram of predicted probabilities you can judge how useful the model is, i.e., how often it provides more definitive probabilities that are closer to 0 or 1 than to 0.5. – Frank Harrell Nov 12 '11 at 14:28

Concordance is defined by checking the mean score, not the maximum score.

So for your examples, mean scores for 1 is 0*55.1% + 1*35.8% + 2*9.1% = 0.54, and 2 is (through similar calculations) 0.826.

It is this value that you should compare to get the concordance or any other association statistics.

Ref -

share|improve this answer
No, concordance is computed using the raw outcome variable $Y$ and the linear predictor $X\beta$ or any of the predicted probabilities (as they are all monotonically related to each other, i.e., just shifted by in intercept before computing the expit). Somers' $D_{xy}$ rank correlation coefficient uses this concordance measure. – Frank Harrell Nov 9 '11 at 14:02
P.S. Note that the SAS documentation, which is a re-write of my original documentation for the precursor SAS PROC LOGIST I wrote many years ago, is now incorrect, and its formula for the mean is incorrect unless $Y$ consists of consecutive integers. – Frank Harrell Nov 9 '11 at 14:13
Thanks for the info, guys. – gung Nov 9 '11 at 17:13

Your Answer


By posting your answer, you agree to the privacy policy and terms of service.

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