###################################################
### chunk number 1: 
###################################################
library(siggenes)


###################################################
### chunk number 2: 
###################################################
library(siggenes)


###################################################
### chunk number 3: 
###################################################
data(golub)


###################################################
### chunk number 4: 
###################################################
args(sam)


###################################################
### chunk number 5: 
###################################################
args(d.stat)


###################################################
### chunk number 6: 
###################################################
args.sam(summary)


###################################################
### chunk number 7: 
###################################################
n <- 10
rep(1, n)


###################################################
### chunk number 8: 
###################################################
n1 <- n2 <- 5
rep(c(0, 1), c(n1, n2))


###################################################
### chunk number 9: 
###################################################
K <- 5
c((-1:-5), 1:5)


###################################################
### chunk number 10: 
###################################################
K <- 5
rep(1:K, e = 2) * rep(c(-1 ,1), K)


###################################################
### chunk number 11: 
###################################################
K <- 5
cbind(rep(c(-1, 1), 5), rep(1:5, e = 2))


###################################################
### chunk number 12: 
###################################################
sam.out <- sam(golub, golub.cl, rand = 123, gene.names = golub.gnames[,3])
sam.out


###################################################
### chunk number 13: 
###################################################
sam.out2 <- sam(golub, golub.cl, method = wilc.stat, rand = 123)


###################################################
### chunk number 14: 
###################################################
summary(sam.out)


###################################################
### chunk number 15: 
###################################################
print(sam.out, seq(1.5, 2.4, 0.1))


###################################################
### chunk number 16: 
###################################################
sum.sam.out <- summary(sam.out, 3.3)
sum.sam.out


###################################################
### chunk number 17: 
###################################################
print(sum.sam.out, varNames = "Proteins")


###################################################
### chunk number 18: 
###################################################
sum.sam.out@row.sig.genes


###################################################
### chunk number 19: 
###################################################
sum.sam.out@mat.fdr


###################################################
### chunk number 20: 
###################################################
sum.sam.out@mat.sig


###################################################
### chunk number 21: 
###################################################
list.siggenes(sam.out, 3.3)


###################################################
### chunk number 22: 
###################################################
findDelta(sam.out, fdr = 0.05)


###################################################
### chunk number 23: 
###################################################
findDelta(sam.out, genes = 200)


###################################################
### chunk number 24: 
###################################################
find.out <- find.a0(golub, golub.cl, rand = 123)


###################################################
### chunk number 25: 
###################################################
find.out


###################################################
### chunk number 26: 
###################################################
print(find.out, 0.95)


###################################################
### chunk number 27: 
###################################################
ebam(find.out)


###################################################
### chunk number 28: 
###################################################
ebam(find.out, which.a0 = 2)


###################################################
### chunk number 29: 
###################################################
ebam(golub, golub.cl, a0 = 0, fast = TRUE, rand = 123)


###################################################
### chunk number 30: 
###################################################
ebam.out <- ebam(golub, golub.cl, a0 = 0, rand = 123)


###################################################
### chunk number 31: 
###################################################
print(ebam.out, seq(0.91, 0.99, 0.01))


###################################################
### chunk number 32: 
###################################################
summary(ebam.out, 0.99997)


###################################################
### chunk number 33: 
###################################################
ebam(golub, golub.cl, a0 = 0, var.equal = TRUE, rand = 123)


###################################################
### chunk number 34: 
###################################################
ebam(golub, golub.cl, quan.a0 = 0.5, rand = 123)


###################################################
### chunk number 35: 
###################################################
ebam(golub, golub.cl, method = wilc.ebam, rand =123)


###################################################
### chunk number 36: 
###################################################
t.stat <- function(data, cl){
    require(genefilter) ||
        stop("genefilter required.")
    cl <- as.factor(cl)
    row.out <- rowttests(data, cl)
    d <- row.out$statistic
    m <- length(na.exclude(d))
    d.bar <- qt(((1:m) - 0.5)/m, length(cl) - 2)
    p.value <- row.out$p.value
    vec.false <- m * p.value/2
    s <- row.out$dm/d
    msg <- paste("SAM Two-Class Analysis",
         "Assuming Normality\n\n")
    list(d = -d, d.bar = d.bar, p.value = p.value,
        vec.false = vec.false, s = s, s0 = 0,
        mat.samp = matrix(numeric(0)),
        msg = msg, fold = numeric(0))
}


###################################################
### chunk number 37: 
###################################################
sam(golub, golub.cl, method = t.stat)


###################################################
### chunk number 38: 
###################################################
t.find <- function(data, cl, B = 50){
    require(genefilter)
    z.fun <- function(data, cl){
        cl <- as.factor(cl)
        out <- rowttests(data, cl)
        r<- out$dm
        s<- r / out$statistic
        return(list(r = -r, s = s))
    }
    mat.samp <- matrix(0, B, length(cl))
    for(i in 1:B)
        mat.samp[i, ] <- sample(cl)
    z.out <- z.fun(data, cl)
    msg <- paste("EBAM Analysis with a Moderated t-Statistic\n\n")
    list(r = z.out$r, s = z.out$s,
        mat.samp = mat.samp, z.fun = z.fun, msg = msg)
}


###################################################
### chunk number 39: 
###################################################
t.out <- find.a0(golub, golub.cl, method = t.find, B = 100, rand =123)
t.out


###################################################
### chunk number 40: 
###################################################
find.a0(golub, golub.cl, var.equal = TRUE, rand =123)


###################################################
### chunk number 41: 
###################################################
ebam(t.out)


###################################################
### chunk number 42: 
###################################################
t.ebam<-function(data, cl){
    require(genefilter)
    cl <- as.factor(cl)
    out <- rowttests(data, cl)
    z <- -out$statistic
    z.dens <- denspr(z)$y
    m <- length(z)
    vec.pos <- m * out$p.value / 2
    z.null <- dt(z, length(cl) - 2)
    msg<-paste("EBAM Analysis with t-Statistic Assuming Normality.\n\n")
    list(z = z, ratio = z.null/z.dens, vec.pos = vec.pos,
        vec.neg = vec.pos, msg = msg)
}


###################################################
### chunk number 43: 
###################################################
ebam(golub, golub.cl, method = t.ebam)