0
$\begingroup$

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,]
dd <- datadist(training)
options(datadist = "dd")
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
$\endgroup$
4
  • $\begingroup$ lrm is for binary or ordinal $Y$. Your $Y$ is nominal/multinomial/polytomous. $\endgroup$ – Frank Harrell Apr 28 '18 at 11:46
  • $\begingroup$ 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. $\endgroup$ – josh myers Apr 28 '18 at 18:07
  • $\begingroup$ 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. $\endgroup$ – Frank Harrell Apr 28 '18 at 22:04
  • $\begingroup$ 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. $\endgroup$ – josh myers Apr 30 '18 at 7:02
0
$\begingroup$
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)
$\endgroup$
3
  • $\begingroup$ 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. $\endgroup$ – Frank Harrell Apr 30 '18 at 12:16
  • $\begingroup$ 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? $\endgroup$ – josh myers Apr 30 '18 at 23:09
  • 1
    $\begingroup$ 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. $\endgroup$ – Frank Harrell May 1 '18 at 11:42

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

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