### R code from vignette source 'vignettes/SSPA/inst/doc/SSPA.Rnw'
### Encoding: UTF-8

###################################################
### code chunk number 1: options
###################################################
options(continue=" ")


###################################################
### code chunk number 2: SSPA.Rnw:56-59
###################################################
library(multtest)
data(golub)  
table(golub.cl)    


###################################################
### code chunk number 3: SSPA.Rnw:64-65
###################################################
tst <- mt.teststat(golub, golub.cl) 


###################################################
### code chunk number 4: SSPA.Rnw:70-72
###################################################
library(SSPA)
pd <- pilotData(name="ALL/AML",testStatistics=tst,sampleSizeA=11,sampleSizeB=27)


###################################################
### code chunk number 5: SSPA.Rnw:77-78
###################################################
pd


###################################################
### code chunk number 6: SSPA.Rnw:85-89
###################################################
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)


###################################################
### code chunk number 7: SSPA.Rnw:98-100
###################################################
ss <- sampleSize(pd)
ss


###################################################
### code chunk number 8: SSPA.Rnw:107-108
###################################################
plotEffectSize(ss, type='l')


###################################################
### code chunk number 9: SSPA.Rnw:120-128
###################################################
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)


###################################################
### code chunk number 10: SSPA.Rnw:175-185
###################################################
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)


###################################################
### code chunk number 11: SSPA.Rnw:191-196
###################################################
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)


###################################################
### code chunk number 12: SSPA.Rnw:204-214
###################################################
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


###################################################
### code chunk number 13: SSPA.Rnw:224-236
###################################################
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)


###################################################
### code chunk number 14: SSPA.Rnw:243-247
###################################################
ss <- sampleSize(pd)
ss
ssMod <- sampleSize(pdMod)
ssMod


###################################################
### code chunk number 15: SSPA.Rnw:255-261
###################################################
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)


###################################################
### code chunk number 16: SSPA.Rnw:274-281 (eval = FALSE)
###################################################
## 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)


###################################################
### code chunk number 17: SSPA.Rnw:318-319
###################################################
toLatex(sessionInfo())