########################################################## # # FpAUC index: FITTED PARTIAL AREA UNDER THE ROC CURVE # FOR HIGH SENSITIVITY INTERVAL (TPR0, 1) # The FpaucHS function computes the FpAUC index # # Notation: # xsample = Gold standard # ysample = test scores # lower.sen = tpr0 (the minimum sensitivity required)# ########################################################## ########################################################## # CHECK THE GOLD STANDARD # The 1st is the negative group # The 2nd is the positive group CheckGoldStandard <- function(xsample, positive=NULL) { xs=NULL; xs.levels=NULL; xs <- as.factor(xsample) xs.levels <- levels(xs) if (length(xs.levels) !=2 ) { stop("Error in the Gold Standard: it must be binary") } if (!is.null(positive)) { if (!any(xs.levels == positive)) { stop("Error in Gold Standard: positive case must be one of the levels") } } else { positive <- xs.levels[2] } if (xs.levels[2] != positive) { xs <- relevel(xs, ref=positive) } xs <- (xs == positive)*1 return(xs)} ########################################################## # POINTS OF THE ROC CURVE # ROCpoints function calculates the coordinates (FPR, TPR) ROCpoints <- function(xsample, ysample) { fpr.p=NULL; sen.p=NULL; xy.roc=NULL; x.p=NULL; y.p=NULL; pts=NULL; x.p <- xsample[which(is.na(xsample)==FALSE & is.na(ysample)==FALSE)] y.p <- ysample[which(is.na(xsample)==FALSE & is.na(ysample)==FALSE)] x.p <- CheckGoldStandard(x.p) pts <- sort(y.p) pts <- append(pts[-length(pts)]+diff(pts)/2, min(y.p)-1, 0) pts <- append(pts, max(y.p)+1, length(pts)) for (i.pt in 1:length(pts)) { pre.p <- NULL pt=NULL pt <- pts[i.pt] pre.p <- (y.p > pt)*1 fpr.p[i.pt] <- sum((pre.p == 1)*(x.p == 0))/sum(x.p == 0) sen.p[i.pt] <- sum((pre.p == 1)*(x.p == 1))/sum(x.p == 1)} if (is.unsorted(sen.p)) { sen.p <- rev(sen.p) fpr.p <- rev(fpr.p)} xy.roc <- cbind(fpr=fpr.p, tpr=sen.p) return(xy.roc)} ########################################################## # POINTS OF THE PARTIAL ROC CURVE FOR HIGH SENSITIVITY # pHSpoints function calculates the points in (tpr0, 1) pHSpoints <- function(xsample, ysample, lower.sen) { pts.roc=NULL; fpr.roc=NULL; sen.roc=NULL; fpr.p=NULL; sen.p=NULL; ppoints=NULL; i.low=NULL; j.low=NULL; lscale=NULL; if ((lower.sen >= 1) || (lower.sen<0)) {stop("Error in the prefixed TPR range")} pts.roc <- ROCpoints(xsample, ysample) fpr.roc <- pts.roc[,1] sen.roc <- pts.roc[,2] i.low <- min(which(sen.roc >= lower.sen)) j.low <- max(i.low -1, 1) fpr.p <- fpr.roc[i.low:length(fpr.roc)] sen.p <- sen.roc[i.low:length(sen.roc)] if ((sen.roc[i.low] > lower.sen) && (i.low>1)) { sen.p <- append(sen.p, lower.sen, 0) lscale <- (sen.p[1]-sen.roc[j.low])/(sen.roc[i.low]-sen.roc[j.low]) fpr.p <- append(fpr.p, fpr.roc[j.low]+(fpr.roc[i.low]- fpr.roc[j.low])*lscale, 0)} ppoints=cbind(fprp=fpr.p, tprp=sen.p) return(ppoints)} ########################################################## # PARTIAL AREA UNDER THE ROC CURVE FOR HIGH SENSITIVITY # paucHS function calculates the pAUC between (tpr0, 1) paucHS <- function(xsample, ysample, lower.sen) { pts.proc=NULL; fpr.p=NULL; sen.p=NULL; esp.p=NULL; auc.p=NULL; pts.proc <- pHSpoints(xsample, ysample, lower.sen) fpr.p <- pts.proc[,1] sen.p <- pts.proc[,2] esp.p <- 1-fpr.p auc.p <- sum(diff(sen.p)*apply(cbind(esp.p[-1],esp.p[-length(esp.p)]), 1, mean)) return(auc.p)} ########################################################## # PARTIAL NEGATIVE LIKELIHOOD RATIO FOR HIGH SENSITIVITY # pNLR_HS function calculates the NLR between (tpr0, 1) pNLR_HS <- function(xsample, ysample, lower.sen) { pts.proc=NULL; fpr.p=NULL; sen.p=NULL; pnlr=NULL; pts.proc <- pHSpoints(xsample, ysample, lower.sen) fpr.p <- pts.proc[,1] sen.p <- pts.proc[,2] fpr.p <- fpr.p[-length(fpr.p)] sen.p <- sen.p[-length(sen.p)] pnlr <- (1-sen.p)/(1-fpr.p) return(pnlr)} ########################################################## # PARTIAL ROC CURVE SHAPE FOR HIGH SENSITIVITY # shapepROC function classifies the ROC curve by its pNLR shapepROC <- function(xsample, ysample, lower.sen) { nlr.p=NULL; nlr0=NULL; sproc=NULL; nlr.p <- pNLR_HS(xsample, ysample, lower.sen) nlr0 <- nlr.p[1] if ((all(nlr.p <= nlr0)) & (is.finite(nlr0))) { sproc <- "BpNLR" } else { if (all(nlr.p <= 1)) { sproc <- "pProp" } else { sproc <- "other" } } return(sproc)} ########################################################## # FITTED PARTIAL AREA INDEX FOR HIGH SENSITIVITY # FpaucHS function calculates the FpAUC index in (tpr0, 1) FpaucHS <- function(xsample, ysample, lower.sen) { fpr.p=NULL; sen.p=NULL; pauc.p=NULL; sproc.p=NULL; fpauc.p=NULL; fpr.p <- pHSpoints(xsample, ysample, lower.sen)[,1] sen.p <- pHSpoints(xsample, ysample, lower.sen)[,2] pauc.p <- paucHS(xsample, ysample, lower.sen) sproc.p <- shapepROC(xsample, ysample, lower.sen) fpauc.p <- switch(sproc.p, BpNLR = pauc.p/ ((1-fpr.p[1])*(1-sen.p[1])), pProp = (pauc.p+(1-sen.p[1])*(sen.p[1]-fpr.p[1]))/ ((1+sen.p[1]-2*fpr.p[1])*(1-sen.p[1])), other = (pauc.p+(1-sen.p[1])*(1-fpr.p[1]))/ (2*(1-fpr.p[1])*(1-sen.p[1])) ) return(fpauc.p)} ##########################################################