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
Improve this question
$\endgroup$
4
  • $\begingroup$ lrm is for binary or ordinal $Y$. Your $Y$ is nominal/multinomial/polytomous. $\endgroup$
    – Frank Harrell
    Commented Apr 28, 2018 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$
    – tauft
    Commented Apr 28, 2018 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
    Commented Apr 28, 2018 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$
    – tauft
    Commented Apr 30, 2018 at 7:02

1 Answer 1

Reset to default
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)
Improve this answer
$\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
    Commented Apr 30, 2018 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$
    – tauft
    Commented Apr 30, 2018 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
    Commented May 1, 2018 at 11:42

Your Answer

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

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