Changing the colour of a calibration plot

Question

I've been generating calibration plots for my cph models of survival data. However, the default setting puts the "ideal" line in grey, which makes it difficult to discriminate. I've tried to specify the colour parameters in plot(), but this obviously only changes the line for "observed". What can I pass in plot() to change the line of the "ideal" line in a calibration plot generated in rms?

Answer 1

Here is one option:

Let's say you have code to create a cph model of survival data and use calibrate from the rms package:

library(rms)

set.seed(1)
n <- 200
d.time <- rexp(n)
x1 <- runif(n)
x2 <- factor(sample(c('a', 'b', 'c'), n, TRUE))
f <- cph(Surv(d.time) ~ pol(x1,2) * x2, x=TRUE, y=TRUE, surv=TRUE,time.inc=1.5)
cal <- calibrate(f, u=1.5, cmethod='KM', m=50, B=20)

This will generate a calibrate object:

R> class(cal)
[1] "calibrate"

If you are using plot on this object, you can discover the function being called in rms:

R> getAnywhere("plot.calibrate.default")
A single object matching ‘plot.calibrate.default’ was found
It was found in the following places
  registered S3 method for plot from namespace rms
  namespace:rms
with value

function (x, xlab, ylab, xlim, ylim, legend = TRUE, subtitles = TRUE, 
    cex.subtitles = 0.75, riskdist = TRUE, scat1d.opts = list(nhistSpike = 200), 
    ...) 

You can create your own function based on this function, and alter the color of the ideal line. In this case, we make the ideal line green (and revise the text labels to match):

myplot <- function (x, xlab, ylab, subtitles = TRUE, conf.int = TRUE, cex.subtitles = 0.75, 
          riskdist = TRUE, add = FALSE, scat1d.opts = list(nhistSpike = 200), 
          par.corrected = NULL, ...) 
{
  at <- attributes(x)
  u <- at$u
  units <- at$units
  if (length(par.corrected) && !is.list(par.corrected)) 
    stop("par.corrected must be a list")
  z <- list(col = "blue", lty = 1, lwd = 1, pch = 4)
  if (!length(par.corrected)) 
    par.corrected <- z
  else for (n in setdiff(names(z), names(par.corrected))) par.corrected[[n]] <- z[[n]]
  predicted <- at$predicted
  if ("KM" %in% colnames(x)) {
    type <- "stratified"
    pred <- x[, "mean.predicted"]
    cal <- x[, "KM"]
    cal.corrected <- x[, "KM.corrected"]
    se <- x[, "std.err"]
  }
  else {
    type <- "smooth"
    pred <- x[, "pred"]
    cal <- x[, "calibrated"]
    cal.corrected <- x[, "calibrated.corrected"]
    se <- NULL
  }
  un <- if (u == 1) 
    paste(units, "s", sep = "")
  else units
  if (missing(xlab)) 
    xlab <- paste("Predicted ", format(u), units, "Survival")
  if (missing(ylab)) 
    ylab <- paste("Fraction Surviving ", format(u), " ", 
                  un, sep = "")
  if (length(se) && conf.int) {
    ciupper <- function(surv, d) ifelse(surv == 0, 0, pmin(1, 
                                                           surv * exp(d)))
    cilower <- function(surv, d) ifelse(surv == 0, 0, surv * 
                                          exp(-d))
    errbar(pred, cal, cilower(cal, 1.959964 * se), ciupper(cal, 
                                                           1.959964 * se), xlab = xlab, ylab = ylab, type = "b", 
           add = add, ...)
  }
  else if (add) 
    lines(pred, cal, type = if (type == "smooth") 
      "l"
      else "b")
  else plot(pred, cal, xlab = xlab, ylab = ylab, type = if (type == 
                                                            "smooth") 
    "l"
    else "b", ...)
  err <- NULL
  if (riskdist && length(predicted)) {
    do.call("scat1d", c(list(x = predicted), scat1d.opts))
    if (type == "smooth") {
      s <- !is.na(pred + cal.corrected)
      err <- predicted - approxExtrap(pred[s], cal.corrected[s], 
                                      xout = predicted, ties = mean)$y
    }
  }
  if (subtitles && !add) {
    if (type == "smooth") {
      Col <- par.corrected$col
      substring(Col, 1, 1) <- toupper(substring(Col, 1, 
                                                1))
      title(sub = sprintf("Black: observed  Green: ideal\n%s : optimism corrected", 
                          Col), adj = 0, cex.sub = cex.subtitles)
      w <- if (length(err)) 
        paste("B=", at$B, " based on ", at$what, "\nMean |error|=", 
              round(mean(abs(err)), 3), "  0.9 Quantile=", 
              round(quantile(abs(err), 0.9, na.rm = TRUE), 
                    3), sep = "")
      else paste("B=", at$B, "\nBased on ", at$what, sep = "")
      title(sub = w, adj = 1, cex.sub = cex.subtitles)
    }
    else {
      title(sub = paste("n=", at$n, " d=", at$d, " p=", 
                        at$p, ", ", at$m, " subjects per group\nGreen: ideal", 
                        sep = ""), adj = 0, cex.sub = cex.subtitles)
      title(sub = paste("X - resampling optimism added, B=", 
                        at$B, "\nBased on ", at$what, sep = ""), adj = 1, 
            cex.sub = cex.subtitles)
    }
  }
  abline(0, 1, col = "green")
  if (type == "stratified") 
    points(pred, cal.corrected, pch = par.corrected$pch, 
           col = par.corrected$col)
  else lines(pred, cal.corrected, col = par.corrected$col, 
             lty = par.corrected$lty, lwd = par.corrected$lwd)
  invisible()
}

Then you can use your custom function with your calibrate object:

myplot(cal)