Plotfun.R
Code snipits and programs written in R, S or S-PLUS extractstats <- function(data, what="FDR", alpha=0.05){
- Options for what: U, V, T, S, R, UT, FDR, FNDR, TaddV, PrRgt0, Power, PrReq1, p0, p0NW, Alphap0
- 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
- if(what=="V" || what=="T" || what=="p0" || what=="p0NW") {
- return(get(what))
- }
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)) }
- extractstats2(F331, what="PrRgt0")
- 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, ...){
- 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)
}
}
- if(hexbinning){
require(hexbin)
- }
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)
- CIfudge <- c(seq(1,195,length=98),196:201)
xbin <- rep(p, each=length(data1)) ybin <- unlist(data)
- 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()
- 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] }
- 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] }
