# How to plot the calibration curve for an ordinal logistic regression model applied to a test sample?

I'm doing a validation study of an ordinal logistic regression model that was made with the lrm function of the rms package in R. How can I plot the calibration curve for the model when applied to new data? I want to create the cal1 and cal2 plots below (without bootstrapping), but using a new sample:

library(rms)
library(heplots)
data(Diabetes)
d = Diabetes
id = 1:145
d$id = id training = d[sample(1:nrow(d), 75, replace = F),] testing = d[-training$id,]
mod = lrm(group ~ glufast, data = training, x = T, y = T)
# normal is reference level, Chemical_Diabetic is middle level, and then Overt_Diabetic
cal1 = plot(calibrate(mod, kint = 1)) # plot for group >= Chemical_Diabetic
cal2 = plot(calibrate(mod, kint = 2)) # plot for group >= Overt_Diabetic


I have previously used the val.prob function to calibrate a new sample with a binary response, but this function doesn't take an ordinal outcome. I tried changing my response to be binary for each level (y >= Chemical_Diabetic and y >= Overt_Diabetic), but it seems incorrect:

pred.test = as.data.frame(predict(mod, testing, type = "fitted"))
chemical.pred = pred.test[,1] # select predictions for >= Chemical_Diabetic
overt.pred = pred.test[,2] # select predictions for >= Overt_Diabetic
y = as.numeric(testing$group) y.chemical <- replace(y, y==3, 2) # set Overt_Diabetic to Chemical_Diabetic y.overt <- replace(y, y==2, 1) # set Chemical_Diabetic to Normal val.prob(chemical.pred, y.chemical, logistic.cal = F, smooth = T) # trying to be cal1 val.prob(overt.pred, y.overt, logistic.cal = F, smooth = T) # trying to be cal2  • lrm is for binary or ordinal$Y$. Your$Y$is nominal/multinomial/polytomous. – Frank Harrell Apr 28 '18 at 11:46 • I have updated my example with a more appropriate dataset. It now uses the Diabetes data from the heplots package which has an ordinal outcome. – josh myers Apr 28 '18 at 18:07 • RMS course notes chapter on ordinal regression will get you much of the way there. You ask for predicted probabilities for a specific intercept with a special argument to predict, then use val.prob to validate that against actual Y >= cutoff binary variable. – Frank Harrell Apr 28 '18 at 22:04 • Thanks Professor Harrell. According to the predict.lrm documentation, you can only specify a specific intercept (the kint argument) for type = 'lp'. I believe I need the fitted values for the calibration plot, in which case the output is a list with an element for prob Y >= j for each intercept. – josh myers Apr 30 '18 at 7:02 ## 1 Answer library(rms) library(heplots) data(Diabetes) d = Diabetes id = 1:145 d$id = id
training = d[sample(1:nrow(d), 75, replace = F),]
testing = d[-training$id,] dd <- datadist(training) options(datadist = "dd") mod = lrm(group ~ glufast, data = training, x = T, y = T) # normal is reference level, Chemical_Diabetic is 2nd level, and then Overt_Diabetic cal.plot = plot(calibrate(mod, kint = 1)) # plot for Y >= Chemical_Diabetic # predict on new data pred.test = as.data.frame(predict(mod, testing, type = "fitted")) chemical.pred = pred.test[,1] # select predictions for >= Chemical_Diabetic y.chemical = as.numeric(testing$group)
y.chemical <- replace(y.chemical, y.chemical==1, 0) # make normal 0, diseased 1
y.chemical <- replace(y.chemical, y.chemical==2, 1)
y.chemical <- replace(y.chemical, y.chemical==3, 1) # set Overt_Diabetic to Chemical_Diabetic
cal.plot.test = val.prob(chemical.pred, y.chemical)

• Just keep in mind that it's very inefficient to use split-sample validation unless you have a huge dataset to begin with (say n > 20,000). Details are here. – Frank Harrell Apr 30 '18 at 12:16
• Yes, in my dataset I have 2 observations for each subject, 1 right, 1 left (no time difference). My approach has been to sample 1 observation for model development, and use the other for validation (rather than just discarding those data). It may be more efficient to use both observations develop the model, but I think that would require a hierarchical model. Can I fit that sort of model with rms? – josh myers Apr 30 '18 at 23:09
• Yes you can use a GEE-like approach with an after-the-fit intra-subject correlation adjustment using the robcov or bootcov functions. But a Bayesian hierarchical model would have better performance and would result in exact inference. – Frank Harrell May 1 '18 at 11:42