###################################################
### chunk number 1: options
###################################################
#line 49 "vignettes/SSPA/inst/doc/SSPA.Rnw"
options(continue=" ")


###################################################
### chunk number 2: 
###################################################
#line 57 "vignettes/SSPA/inst/doc/SSPA.Rnw"
library(multtest)
data(golub)  
table(golub.cl)    


###################################################
### chunk number 3: 
###################################################
#line 65 "vignettes/SSPA/inst/doc/SSPA.Rnw"
tst <- mt.teststat(golub, golub.cl) 


###################################################
### chunk number 4: 
###################################################
#line 71 "vignettes/SSPA/inst/doc/SSPA.Rnw"
library(SSPA)
pd <- pilotData(name="ALL/AML",testStatistics=tst,sampleSizeA=11,sampleSizeB=27)


###################################################
### chunk number 5: 
###################################################
#line 78 "vignettes/SSPA/inst/doc/SSPA.Rnw"
pd


###################################################
### chunk number 6: 
###################################################
#line 86 "vignettes/SSPA/inst/doc/SSPA.Rnw"
par(mar=c(4.1, 4.1, 2.1, 2.1))
layout(matrix(c(1,2),nrow=2))
hist(pd, cex.main=1)
plot(pd, cex.main=1)


###################################################
### chunk number 7: 
###################################################
#line 99 "vignettes/SSPA/inst/doc/SSPA.Rnw"
ss <- sampleSize(pd)
ss


###################################################
### chunk number 8: 
###################################################
#line 108 "vignettes/SSPA/inst/doc/SSPA.Rnw"
plotEffectSize(ss, type='l')


###################################################
### chunk number 9: 
###################################################
#line 121 "vignettes/SSPA/inst/doc/SSPA.Rnw"
layout(matrix(c(1:2), nrow=2))
par(mar=c(4.1, 4.1, 2.1, 2.1))
pwr <- Power(ss, plot = FALSE, samplesizes = c(5, 10, 15, 20), fdr=0.01) 
plot(c(5, 10, 15, 20), pwr, ylim = c(0, 1), type = "b", ylab = "Power", xlab = "Sample size per group") 
legend("bottomright", colnames(pwr), col=c(1:ncol(pwr)), pch=1, lty=1)
pwr <- Power(ss, plot = FALSE, samplesizes = c(5, 10, 15, 20), fdr=c(0.01, 0.05)) 
matplot(c(5, 10, 15, 20), pwr, ylim = c(0, 1), type = "b", pch=1, ylab = "Power", xlab = "Sample size per group")
legend("bottomright", colnames(pwr), col=c(1:ncol(pwr)), pch=1, lty=1)


###################################################
### chunk number 10: 
###################################################
#line 176 "vignettes/SSPA/inst/doc/SSPA.Rnw"
library(marray)
library(convert)
library(limma)
data(swirl)
swirl <- as(swirl, "RGList")
MA <- normalizeWithinArrays(swirl) #print tip loess
design <- c(-1,1,-1,1)
fit <- lmFit(MA, design)
ordinary.t <- fit$coef/fit$stdev.unscaled/fit$sigma
fitModerated <- eBayes(fit)


###################################################
### chunk number 11: 
###################################################
#line 192 "vignettes/SSPA/inst/doc/SSPA.Rnw"
par(mfcol=c(1,2))
qqt(fitModerated$t, df=fitModerated$df.prior + fitModerated$df.residual, pch=16, cex=0.5, main="Moderated t")
abline(0,1)
qqt(ordinary.t, df=fit$df.residual, pch=16, cex=0.5, main="Ordinary t")
abline(0,1)


###################################################
### chunk number 12: 
###################################################
#line 205 "vignettes/SSPA/inst/doc/SSPA.Rnw"
library(SSPA)

nu <- fit$df.residual[1] 
nu0 <- fitModerated$df.prior

pd <- pilotData(name = "Swirl", testStatistics = ordinary.t[,1], sampleSizeA = 2, sampleSizeB = 2, dof=nu, nullDist = "student")
pd

pdMod <- pilotData(name = "Swirl", testStatistics = fitModerated$t[,1], sampleSizeA = 2 + nu0/2, sampleSizeB = 2 + nu0/2, dof= nu + nu0, nullDist = "student")
pdMod


###################################################
### chunk number 13: 
###################################################
#line 225 "vignettes/SSPA/inst/doc/SSPA.Rnw"
par(mfcol=c(1,2))
plot(abs(pd@testStatistics), abs(pdMod@testStatistics), 
xlab="|ordinary t test statistics|", 
ylab="|moderated t test statistics|", pch=16, cex=0.5)
abline(a=0, b=1, col="red", lwd=2, lty=2)
plot(-log10(pd@pValues), -log10(pdMod@pValues), xaxt="n", yaxt='n', xlab=expression(paste(-log[10],"(ordinary p-values)")), 
ylab=expression(paste(-log[10],"(moderated p-values)")), pch=16, cex=0.5)
at <- axTicks(2)
axis(2, at=at, tcl=-1,	labels=parse(text=c("10^0", paste("10^-", at[-1], sep=""))))			
at <- axTicks(1)
axis(1, at=at, tcl=-1,	labels=parse(text=c("10^0", paste("10^-", at[-1], sep=""))))			
abline(a=0, b=1, col="red", lwd=2, lty=2)


###################################################
### chunk number 14: 
###################################################
#line 244 "vignettes/SSPA/inst/doc/SSPA.Rnw"
ss <- sampleSize(pd)
ss
ssMod <- sampleSize(pdMod)
ssMod


###################################################
### chunk number 15: 
###################################################
#line 256 "vignettes/SSPA/inst/doc/SSPA.Rnw"
plotEffectSize(ssMod, type='l' , lwd=2, sub=NULL)
lines(ss@theta, ss@lambda, col="red", lwd=2)
pi0Mod <- paste("pi0: ", signif(ssMod@pi0[[2]], 3))
pi0 <- paste("pi0: ", signif(ss@pi0[[2]], 3))
labels <- c(paste("moderated ", pi0Mod), paste("ordinary ", pi0))
legend("topleft", labels, col=c("black", "red"), pch=1)


###################################################
### chunk number 16:  eval=FALSE
###################################################
## #line 275 "vignettes/SSPA/inst/doc/SSPA.Rnw"
## samplesizes <- c(2, 4, 6, 8, 10)
## pwr <- Power(ss, plot = FALSE, samplesizes = samplesizes, fdr=c(0.05, 0.1))
## pwrMod <- Power(ssMod, plot = FALSE, samplesizes = samplesizes + nu0/2, fdr=c(0.05, 0.1))
## matplot(samplesizes, pwr, ylim = c(0, 1), type = "b", col=1, pch=1, ylab = "Power", xlab = "Sample size per group", main="Power Curves")
## matlines(samplesizes, pwrMod, col=2, type="b", pch=1)
## legend("topleft", colnames(pwr), col=1, lty=c(1, 2))
## legend("bottomright", c("ordinary", "moderated"), col=1:2, lty=1)


###################################################
### chunk number 17: 
###################################################
#line 319 "vignettes/SSPA/inst/doc/SSPA.Rnw"
toLatex(sessionInfo())