# TODO: Add comment # # Author: E.Korsching Dec 7, 2012 ############################################################################### gliding.proportions.table <- function(x, pos.level=c("max","max")){ # setting stepwise different thresholds for paired variables # (assumed that one var represents a hypothesized class and the other a true class) # calculate the corresponding 2x2 contigency table and the performance parameters # x: matrix/data.frame with two variables (col1: hypothesized class: 2 levels, col2: true class: 2 to many levels) # --table()-- is included in this function # pos.level: define as true positive : col1: max/min, col2: max/min : e.g.: c("max","max") #ini nr <- nrow(x) nc <- ncol(x) if(nr<2){ cat("\nMinimum of two rows"); return() } if(nc>2){ cat("\nExact two columns: col1: hypothesized class: 2 levels, col2: true class: 2 to many levels"); return() } if( length(unique(x[,1]))>2 ){ cat("\ncol1: hypothesized class: 2 levels, col2: true class: 2 to many levels"); return() } if( length(unique(x[,2]))<2 ){ cat("\ncol1: hypothesized class: 2 levels, col2: true class: 2 to many levels"); return() } #ini2 x <- table(x[,1],x[,2]) nr <- nrow(x) nc <- ncol(x) #tp pair must be upper left if((pos.level[1]!="min" & pos.level[1]!="max") | (pos.level[1]!="min" & pos.level[1]!="max")){ cat("\nerror in pos.level (min/max)"); return() } if(pos.level[1]=="min"){ if(pos.level[2]=="max"){ x <- x[,nc:1,drop=F] } } if(pos.level[1]=="max"){ if(pos.level[2]=="min"){ x <- x[nr:1,,drop=F] }else{ x <- x[,nc:1,drop=F] x <- x[nr:1,,drop=F] } } # y <- matrix(0,2,2) #second roc matrix if x is greater than 2x2 erg <- data.frame(matrix(0,nc-1,16)) names(erg) <- c("tp","fp","fn","tn","pp","np","p","n","fpr","tpr","spc","fdr","ppv","F","F1","MCC") #roc statistics f1 <- function(x){ # true class, actual value # hypotesized T F # class Y TP FP P' # predicted N FN TN N' # outcome P N erg <- matrix(0,1,16) erg[1,1] <- x[1,1] # tp: true positive erg[1,2] <- x[1,2] # fp: false positive erg[1,3] <- x[2,1] # fn: false negative erg[1,4] <- x[2,2] # tn: true negative erg[1,5] <- sum(x[1, ]) # pp: true positive + false positive erg[1,6] <- sum(x[2, ]) # np: false negative + true negative erg[1,7] <- sum(x[ ,1]) # p: true positive + false negative erg[1,8] <- sum(x[ ,2]) # n: false positive + true negative erg[1,9] <- erg[1,2]/erg[1,8] # fpr: false positive rate (x) erg[1,10] <- erg[1,1]/erg[1,7] # tpr: true positive rate (y) sensitivity hit rate recall erg[1,11] <- 1-erg[1,9] # spc: specificity also tn/n erg[1,12] <- erg[1,2]/(erg[1,2]+erg[1,1]) # fdr: false discovery rate erg[1,13] <- erg[1,1]/(erg[1,2]+erg[1,1]) # ppv: precision positive predictive value erg[1,14] <- 2/(1/erg[1,13]+1/erg[1,10]) # F: measure erg[1,15] <- 2*erg[1,1]/(erg[1,5]+erg[1,7]) # F1: measure erg[1,16] <- (erg[1,1]*erg[1,4]-erg[1,2]*erg[1,3])/sqrt(erg[1,5]*erg[1,6]*erg[1,7]*erg[1,8]) # MCC: Matthews correlation coefficient return(erg) } #fill roc statistics if(nc==2){ erg[1,] <- f1(x) } if(nc>2){ for(i in 2:nc){ y[1,1] <- sum(x[1,(1:(i-1))]) y[2,1] <- sum(x[2,(1:(i-1))]) y[1,2] <- sum(x[1,(i:nc)]) y[2,2] <- sum(x[2,(i:nc)]) erg[(i-1),] <- f1(y) } } #plot par(mfrow=c(2,2)) # tp fp fn tn pp np p n y.max <- max(erg[ ,1:8]) plot(0,0, xlim=c(0,nc), ylim=c(0,y.max), xlab="sequence of thresholds", ylab="counts", type="n", axes=T) leg.str <- c("TP","FP","FN","TN","pp","pn","p","n") lp.col <- palette()[1:8] # lp.col <- colorRampPalette(c("green", "grey", "orange"))( xx ) lp.pch <- c(rep(20,4),rep(4,4)) for(i in 1:(nc-1)){ points(x=i, y=erg[i,1], type="p", col=lp.col[1], cex=0.5, pch=lp.pch[1]) #y: TP points(x=i, y=erg[i,2], type="p", col=lp.col[2], cex=0.5, pch=lp.pch[2]) #y: FP points(x=i, y=erg[i,3], type="p", col=lp.col[3], cex=0.5, pch=lp.pch[3]) #y: FN points(x=i, y=erg[i,4], type="p", col=lp.col[4], cex=0.5, pch=lp.pch[4]) #y: TN points(x=i, y=erg[i,5], type="p", col=lp.col[5], cex=0.5, pch=lp.pch[5]) #y: pp points(x=i, y=erg[i,6], type="p", col=lp.col[6], cex=0.5, pch=lp.pch[6]) #y: pn points(x=i, y=erg[i,7], type="p", col=lp.col[7], cex=0.5, pch=lp.pch[7]) #y: p points(x=i, y=erg[i,8], type="p", col=lp.col[8], cex=0.5, pch=lp.pch[8]) #y: n } legend(x=nc*0.25, y=y.max*0.65, legend=leg.str, col=lp.col, border="black", pch=lp.pch, bty="o", title=NULL, cex=0.5) # FPR TPR plot(0,0, xlim=c(0,1), ylim=c(0,1), xlab="FPR", ylab="TPR", type="n", axes=T) segments(x0=0,y0=0,x1=1,y1=1, col="blue", lty=2, lwd=1.2) #diagonal line leg.str <- paste("FPR, TPR points: ",nc-1,sep="") lp.col <- palette()[2] lp.pch <- 3 for(i in 1:(nc-1)){ points(x=erg[i,9], y=erg[i,10], type="p", col=lp.col, cex=0.5, pch=lp.pch) #x/y: FPR/TPR } legend(x=0.6, y=0.4, legend=leg.str, col=lp.col, border="black", pch=lp.pch, bty="o", title=NULL, cex=0.5) # SPC FDR PPV plot(0,0, xlim=c(0,nc), ylim=c(0,1), xlab="sequence of thresholds", ylab="SPC FDR PPV", type="n", axes=T) leg.str <- c("SPC","FDR","PPV") lp.col <- palette()[1:3] lp.pch <- c(3,3,3) for(i in 1:(nc-1)){ points(x=i, y=erg[i,11], type="p", col=lp.col[1], cex=0.5, pch=lp.pch[1]) #y: SPC points(x=i, y=erg[i,12], type="p", col=lp.col[2], cex=0.5, pch=lp.pch[2]) #y: FDR points(x=i, y=erg[i,13], type="p", col=lp.col[3], cex=0.5, pch=lp.pch[3]) #y: PPV } legend(x=1, y=0.6, legend=leg.str, col=lp.col, border="black", pch=lp.pch, bty="o", title=NULL, cex=0.5) # F F1 MCC plot(0,0, xlim=c(0,nc), ylim=c(-1,1), xlab="sequence of thresholds", ylab="F F1 MCC", type="n", axes=T) leg.str <- c("F","F1","MCC") lp.col <- palette()[1:3] lp.pch <- c(3,3,3) for(i in 1:(nc-1)){ points(x=i, y=erg[i,14], type="p", col=lp.col[1], cex=0.5, pch=lp.pch[1]) #y: F points(x=i, y=erg[i,15], type="p", col=lp.col[2], cex=0.5, pch=lp.pch[2]) #y: F1 points(x=i, y=erg[i,16], type="p", col=lp.col[3], cex=0.5, pch=lp.pch[3]) #y: MCC } legend(x=1, y=1, legend=leg.str, col=lp.col, border="black", pch=lp.pch, bty="o", title=NULL, cex=0.5) # return(erg) } # ac <- gliding.proportions.table(x=aa[[1]], pos.level=c("max","max"))