Plotfun.R

Code snipits and programs written in R, S or S-PLUS extractstats <- function(data, what="FDR", alpha=0.05){

1. Options for what: U, V, T, S, R, UT, FDR, FNDR, TaddV, PrRgt0, Power, PrReq1, p0, p0NW, Alphap0
2. UT equals "m-R"

 U <- S <- R <- UT <- FDR <- FNDR <- TaddV<- PrRgt0 <- p0 <- p0NW <- Alphap0 <- list()
 Power <- PrFDReq1 <- Sens <- Spec <- list()
 Pchar <- names(data$V)
 p <- as.numeric(Pchar)
 m0 <- round(p * data$m)
 names(m0) <- Pchar
 m1 <- round((1-p)*data$m)
 names(m1) <- Pchar

1. if(what=="V" || what=="T" || what=="p0" || what=="p0NW") {
2. return(get(what))
3. }

 if(what=="V"){
   return(data$V)
 }
 if(what=="T"){
   return(data$T)
 }
 if(what=="p0"){
     return(data$p0)
 }	
 if(what=="p0NW"){
     return(data$p0NW)
 }	

 for(i in Pchar){
   if(what=="Alphap0"){
     Alphap0deparse(as.numeric(i)) <-  (alpha * as.numeric(i)) / data$p0i
   }

   if(what=="TaddV"){
     TaddVdeparse(as.numeric(i)) <- (data$Ti+data$Vi)
   }
   
   if(what=="U" || what=="Spec" || what=="UT" || what=="FNDR"){
     Udeparse(as.numeric(i)) <- m0[i]-data$Vi
     if(what=="Spec"){
           tmp <- Ui/m0[i]
           tmp[is.nan(tmp)] <- 0
           Specdeparse(as.numeric(i)) <- tmp
     }
     if(what=="UT" || what=="FNDR"){
       UTdeparse(as.numeric(i)) <- (Ui+data$Ti)
       if(what=="FNDR"){
         tmp <-  data$Ti/(UTi)
         tmp[is.nan(tmp)] <- 0
         FNDRdeparse(as.numeric(i)) <- tmp
       }
     }
   }
   if(what=="S" || what=="Sens" || what=="R" || what=="FDR"  || what=="PrFDReq1" || what=="PrRgt0" || what=="Power"){
   Sdeparse(as.numeric(i)) <- m1[i]-data$Ti
     if(what=="R" || what =="Sens" || what=="FDR" || what=="PrFDReq1" || what=="PrRgt0" || what=="Power"){
       Rdeparse(as.numeric(i)) <- (Si+data$Vi)
       if(what=="Sens"){

tmp <- Si/m1[i] tmp[is.nan(tmp)] <- 0 Sensdeparse(as.numeric(i)) <- tmp

       } 
       
       if(what=="FDR" || what=="PrFDReq1" || what=="Power"){

         tmp <-  data$Vi/(Ri)
         tmp[is.nan(tmp)] <- 0
         FDRdeparse(as.numeric(i)) <- tmp
       }
      
       if(what=="PrFDReq1"){
         tmp <- FDRi==1
         PrFDReq1deparse(as.numeric(i)) <- tmp
       }
       
       if(what=="PrRgt0" || what=="Power"){
         tmp <- Ri>0
         PrRgt0deparse(as.numeric(i)) <- tmp
       }
       if(what == "Power"){
         tmp <- PrRgt0i - FDRi
         Powerdeparse(as.numeric(i)) <- tmp
       }
     }
   }
 }

return(get(what)) }

1. extractstats2(F331, what="PrRgt0")
2. extractstats(F331, what="PrFDReq1")

plotfun <- function(FDlist, what="FDR", probs=c(0.025,0.975), xlim, ylim, add=F, col="black", lty=1, xpos, ypos=0.1, psfile=F,closefile=F,method="", main, CIs = T, hexbinning=T, binNo=100, ...){

1. Problems fixed... main was missing, method has no default if options not found...

if(is.character(psfile)){

 postscript(psfile,horizontal=F, height=297/25.4, width=210/25.4)
 parcalls <- list(...)
 par(pty="s", mai=rep(1,4))
 for(i in names(parcalls)){
   par(i=parcallsi)
 }

}

1. if(hexbinning){

 require(hexbin)

1. }

sim <- deparse(substitute(FDlist))

if(charmatch("EBarrays",sim, nomatch=0)){

 method <- "Empirical Bayes"

} if(charmatch("Multtest",sim, nomatch=0)){

 method <- "Within Gene t-tests"

} if(missing(main)){

 main <- "Model Unknown"

} if(length(grep("GG", sim))){

 main <- "GG-B model"

} if(length(grep("LNN",sim))){

 main <- "LNN-B model"

} if(length(grep("GU", sim))){

 main <- "GU-B model"

} data <- extractstats(FDlist, what)

ylab <- switch(what,

             "FDR"  = "False Discovery Rate",
             "FNDR" = "False Non-Discovery Rate",
             "TaddV"= "# Type I and II errors",
             "U"    = "# not rejected given EE",
             "V"    = "# rejected given EE",
             "T"    = "# not rejected given DE",
             "S"    = "# rejected given DE",
             "UT"   = "# not rejected",
             "R"    = "# rejected",
             "PrRgt0" = "P(R>0)",
             "Power" = "P(R>0) - False Discovery Rate",
             "PrFDReq1" = "P(FDR=1)",
             "Spec" = "Specificity, # rejected over m1",
             "Sens" = "Sensitivity, # not rejected over m0",
             "p0" = "Proportion of genes not changing",               
             "Alphap0" = "alpha* . pi[0]"                 	 
             )

p <- as.numeric(names(data)) alpha <- FDlist$alpha # Check this bit

if(missing(ylim)){

 ylim <- c(0,0.4)

} if(missing(xlim)){

 xlim <- c(0,1)

} if(missing(xpos)){

 xpos <- 0.1

}

lprob <- probs[1] uprob <- probs[2]

xbar <- sapply(data, mean) lCI <- sapply(data, quantile, probs=lprob) uCI <- sapply(data, quantile, probs=uprob)

1. CIfudge <- c(seq(1,195,length=98),196:201)

xbin <- rep(p, each=length(data1)) ybin <- unlist(data)

1. browser()

bins <- hexbin(xbin, ybin, xbins=binNo)

P <- plot(bins, type="n", main=main, ylab = ylab, xlab=expression(paste("Proportion of true null hypotheses ", pi[0], sep=""))) pushHexport(P$plot.vp)

if(add){

 grid.lines(p, xbar)#, col=col, lty=lty)
 if(CIs){ #### WRONG!!!
   grid.lines(p, lCI)#, col="darkgrey",lty=lty)
   grid.lines(p, uCI)#, col="darkgrey",lty=lty)
 }

}else{

 if(hexbinning){
   grid.hexagons(bins)
 } 
 if(what=="FDR") {
   grid.lines(p, rep(alpha, length(p)), default.units="native", gp=gpar(lty=lty, col="red"))
   grid.lines(p, p*rep(alpha, length(p)), default.units="native", gp=gpar(lty=lty, col="red"))
 }

 if(CIs){
   grid.lines(x=p, y=lCI, default.units="native", gp=gpar(col=col,lty=4))
   grid.lines(x=p, y=xbar, default.units="native", gp=gpar(col=col,lty=lty))
   grid.lines(x=p, y=uCI, default.units="native", gp=gpar(col=col,lty=4))
 }

} if(closefile){

   dev.off()
 }

return() }

 save.image()
 
 

1. Little function to find the maximum quickly for any statistic of interest

maxquantile<- function(dset, what, prob, subsetpi0=NULL){ if(!is.null(subsetpi0)){ tmp <- sapply(extractstats(dset, what)[subsetpi0], quantile, prob) }else{ tmp <- sapply(extractstats(dset, what), quantile, prob) } maxquantile <- max(tmp) tmp[tmp==maxquantile] }

minquantile<- function(dset, what, prob, subsetpi0=NULL){ if(!is.null(subsetpi0)){ tmp <- sapply(extractstats(dset, what)[subsetpi0], quantile, prob) }else{ tmp <- sapply(extractstats(dset, what), quantile, prob) } maxquantile <- min(tmp) tmp[tmp==maxquantile] }

1. Little function to find the maximum quickly for any statistic of interest

maxmean<- function(dset, what, subsetpi0=NULL){ if(!is.null(subsetpi0)){ tmp <- sapply(extractstats(dset, what)[subsetpi0], mean) }else{ tmp <- sapply(extractstats(dset, what), mean) }

maxmean <- max(tmp) tmp[tmp==maxmean] }