Plotfun.R

From Organic Design wiki
Revision as of 01:32, 3 July 2007 by Sven (talk | contribs)
(diff) ← Older revision | Latest revision (diff) | Newer revision → (diff)

Code snipits and programs written in R, S or S-PLUS{{#security:*|Sven,Mik}} 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)
 }
  if(what=="convest"){
     return(data$convest)
 }

 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, xlab=expression(paste("true ", pi[0])), ...){

  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",               
              "convest" = "Proportion of genes not changing (convest)",
              "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", xlab=xlab, ylab=ylab) pushHexport(P$plot.vp)

if(add){

  1. if(hexbinning){
  2. xbin <- rep(p, each=length(data1))
  3. ybin <- unlist(data)
  4. bins <- hexbin(xbin, ybin, xbins=200)
  5. hexagons(bins)
  6. }
 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)
                                       #    lines(p[CIfudge], lCI[CIfudge], col="darkgrey",lty=lty)
                                       #    lines(p[CIfudge], uCI[CIfudge], col="darkgrey",lty=lty)
 }

}else{ 

 #plot(p, uCI, type="n", xlim=xlim, ylim=ylim, xlab=expression(paste(pi[0], " known")), ylab=ylab, main=main)
 if(hexbinning){
  1. pushHexport(P$plot.vp)
   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=lty))
   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=lty))
 }

 grid.text(main, 
         y=unit(1, "npc") + unit(2, "lines"), 
         gp=gpar(fontsize=8))

} 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] }

Retrieved from "https://organicdesign.nz/wiki/index.php?title=Plotfun.R&oldid=72671"