R/calroc.R

In IRBLleida/UdBRpackage: Miscellaneous R Functions

#' Funció que avalua la calibració i la dicriminació d'un model logístic
#' 
#' @param model a logistic model (an object of class "glm", usually, a result of a call to glm).
#' @param nOutcome names of the outcome (should be the same than in the original dataset).
#' @param calibration logical. Should be the Hosmer-Lemeshow test for calibration performed? Default is TRUE.
#' @param roc logical. Should be discrimination assessed? Default is TRUE.
#' @param g numerical. Number of groups for the Hosmer-Lemeshow test. Default value is g = 10. Only needed if calibration is TRUE.
#' @param cal.plot logical. Should be the calibration plot obtained? Default is TRUE.
#' @param cal.title string. Title to be included in the calibration plot. Default is NULL (no title is included).
#' @param roc.plot logical. Should be the ROC curve plotted? Default is TRUE.
#' @param roc.title string. Title to be included in the ROC plot. Default is NULL (no title is included).
#' @param cal.name string. Name for the file to save the calibration plot. Should be specified if calibration and cal.plot are TRUE.
#' @param roc.name string. Name for the file to save the ROC plot. Should be specified if roc and roc.plot are TRUE.
#' @param save logical. Should be the plots saved at disk? Default is FALSE.
#' @param path Local directory to save the plots. Default is getwd().
#' @return A list with all results. 
#' @export
calroc <- function(model, nOutcome = NULL, calibration = T, roc = T, g = 10,
                 cal.plot = T, cal.title = NULL, roc.plot = T, roc.title = NULL, cal.name = NULL, roc.name = NULL,
                 add.text.ROC  =  FALSE,
                 p.width = 7, p.height = 4,
                 save = FALSE,
                 path = getwd()){

  if(isTRUE(cal.plot) & is.null(cal.name))cal.name = "CalibrationPlot"
  if(isTRUE(roc.plot) & is.null(roc.name))roc.name = "ROCPlot"

  if(is.null(path))path<-getwd()

  #calibration
  if(isTRUE(calibration)){
    ifelse(is.factor(model$model[,nOutcome]),resp<-as.numeric(model$model[,nOutcome])-1,resp<-model$model[,nOutcome])
    xcal<-hosmerlem(y = resp,yhat = fitted(model),g = g)
    cat("\nCalibration:\n");print(xcal)
    if(!isTRUE(cal.plot))pCal<-NULL
    if(isTRUE(cal.plot)){
      datag<-data.frame(obsrisk = as.numeric(xcal$observed[,2])/(as.numeric(xcal$observed[,1])+as.numeric(xcal$observed[,2])),
                        predrisk = as.numeric(xcal$expected[,2])/(as.numeric(xcal$expected[,1])+as.numeric(xcal$expected[,2])))
      pCal<-ggplot2::ggplot(datag,aes(x = predrisk,y = obsrisk)) + geom_point(cex = 1.5) + xlab("Predicted risk") + ylab("Observed risk") + geom_segment(aes(x = 0,xend = 1,y = 0,yend = 1), lty = 2) +
        geom_text(parse = TRUE, aes(x = 0.6, y = 0.1, label = paste("list(H-L~test:~chi[", xcal$df, "]^2==", deparse(format(roundp(xcal$chisq, 2), nsmall = 2)), ", p-value==", deparse(format(roundp(xcal$p.value, digits = 3), nsmall = 3)), ")"))) +
        ggtitle(cal.title) + theme(plot.title = element_text(hjust  =  0.5))
      #print(pCal)
      # tiff(file = paste(path, "/",cal.name, ".tiff",sep = ""), width = p.width, height = p.height, units  =  'in', res =  300)
      if(isTRUE(save)){
        pdf(file = paste0(path,cal.name, ".pdf"), width = p.width, height = p.height)
        print(pCal)
        dev.off()
      }
    }
  }

  #discrimination
  if(isTRUE(roc)) {
    area <- Hmisc::rcorr.cens(fitted(model), model$model[,nOutcome], outx  =  FALSE)
    cat("\nDiscrimination:\n"); print(area)
    if(!isTRUE(roc.plot))pROC<-NULL
    if(isTRUE(roc.plot)) {
      xroc <- ROCR::prediction(fitted(model), factor(model$model[,nOutcome]))
      xroc <- ROCR::performance(xroc, measure  =  "sens", x.measure  =  "spec")
      datag <- data.frame(sens  =  xroc@y.values[[1]], spe  =  xroc@x.values[[1]])
      pROC <- ggplot2::ggplot(datag, aes(x  =  1-spe, y  =  sens)) + geom_line(lwd  =  1) + xlab("1 - Specificity") + ylab("Sensitivity") + geom_segment(aes(x = 0,xend = 1,y = 0,yend = 1), lty = 2) +
        geom_text(aes(x = 0.6, y = 0.1,label = paste("AUC (95% CI) = ",format(round(area[[1]],3), nsmall  =  3), " (",format(round(area[[1]]-qnorm(1-0.05/2)*area[[3]]/2,3), nsmall  =  3), " - ",format(round(min(1,area[[1]]+qnorm(1-0.05/2)*area[[3]]/2),3), nsmall  =  3), ")",sep = ""))) +
        ggtitle(roc.title) + theme(plot.title  =  element_text(hjust  =  0.5)) +
        scale_x_continuous(breaks = seq(0,1,0.2)) + scale_y_continuous(breaks = seq(0,1,0.2))
      if(isTRUE(add.text.ROC)){
        add <- Epi::ROC(form = formula(model), data = model$model)
        add<-add$res[which.max(rowSums(add$res[,1:2])),]
        pROC <- pROC + geom_text(aes(x = 1-add$spec, y = add$sens, label = paste("Sensitivity  =  ",format(round(add$sens,3), nsmall  =  3), "\nSpecificity  =  ",format(round(add$spec,3), nsmall  =  3),sep = "")),
                                 nudge_x  =  0.0, nudge_y  =  -0.1, hjust  =  0) +
          geom_point(aes(x = 1-add$spec, y = add$sens), color = "red", cex  =  2)
      }
      #print(pROC)
      # tiff(file = paste(path, "/",roc.name, ".tiff",sep = ""), width = p.width, height = p.height, units  =  'in', res =  300)
      if(isTRUE(save)){
        pdf(file = paste0(path,roc.name, ".pdf"), width = p.width, height = p.height)
        print(pROC)
        dev.off()
      }
    }
  }

  if(isTRUE(calibration) & isTRUE(roc)) return(list(calibration = xcal, discrimination = area, plotCalibration = pCal, plotDiscrimination = pROC))
  if(isTRUE(calibration) & !isTRUE(roc)) return(list(calibration = xcal, plotCalibration = pCal))
  if(!isTRUE(calibration) & isTRUE(roc)) return(list(discrimination = area, plotDiscrimination = pROC))
}