###################################################
### chunk number 1: one
###################################################
library(gaga)
set.seed(10)
n <- 100; m <- c(6,6)
a0 <- 25.5; nu <- 0.109
balpha <- 1.183; nualpha <- 1683
probpat <- c(.95,.05)
xsim <- simGG(n,m,p.de=probpat[2],a0,nu,balpha,nualpha,equalcv=TRUE)


###################################################
### chunk number 2: two
###################################################
xsim
featureData(xsim)
phenoData(xsim)
a <- fData(xsim)[,c("alpha.1","alpha.2")]
l <- fData(xsim)[,c("mean.1","mean.2")]
x <- exprs(xsim)


###################################################
### chunk number 3: fig1aplot
###################################################
plot(density(x),xlab='Expression levels',main='')


###################################################
### chunk number 4: fig1bplot
###################################################
plot(l[,1],1/sqrt(a[,1]),xlab='Group 1 Mean',ylab='Group 1 CV')


###################################################
### chunk number 5: fig1a
###################################################
plot(density(x),xlab='Expression levels',main='')


###################################################
### chunk number 6: fig1b
###################################################
plot(l[,1],1/sqrt(a[,1]),xlab='Group 1 Mean',ylab='Group 1 CV')


###################################################
### chunk number 7: three
###################################################
groups <- pData(xsim)$group[c(-6,-12)]
groups
patterns <- matrix(c(0,0,0,1),2,2)
colnames(patterns) <- c('group 1','group 2')
patterns


###################################################
### chunk number 8: threebis
###################################################
patterns <- matrix(c(0,0,0,0,1,1,0,0,1,0,1,2),ncol=3,byrow=TRUE)
colnames(patterns) <- c('CONTROL','CANCER A','CANCER B')
patterns


###################################################
### chunk number 9: four
###################################################
patterns <- matrix(c(0,0,0,1),2,2)
colnames(patterns) <- c('group 1','group 2')
gg1 <- fitGG(x[,c(-6,-12)],groups,patterns=patterns,nclust=1,method='Gibbs',B=1000,trace=FALSE)
gg2 <- fitGG(x[,c(-6,-12)],groups,patterns=patterns,method='EM',trace=FALSE)  
gg3 <- fitGG(x[,c(-6,-12)],groups,patterns=patterns,method='quickEM',trace=FALSE)  


###################################################
### chunk number 10: five
###################################################
gg1 <- parest(gg1,x=x[,c(-6,-12)],groups,burnin=100,alpha=.05)
gg2 <- parest(gg2,x=x[,c(-6,-12)],groups,alpha=.05)
gg3 <- parest(gg3,x=x[,c(-6,-12)],groups,alpha=.05)
gg1
gg1$ci
gg2
gg3


###################################################
### chunk number 11: seven
###################################################
dim(gg1$pp)
gg1$pp[1,]
gg2$pp[1,]


###################################################
### chunk number 12: fig4aplot
###################################################
checkfit(gg1,x=x[,c(-6,-12)],groups,type='data',main='')


###################################################
### chunk number 13: fig4bplot
###################################################
checkfit(gg1,x=x[,c(-6,-12)],groups,type='shape',main='')


###################################################
### chunk number 14: fig4cplot
###################################################
checkfit(gg1,x=x[,c(-6,-12)],groups,type='mean',main='')


###################################################
### chunk number 15: fig4dplot
###################################################
checkfit(gg1,x=x[,c(-6,-12)],groups,type='shapemean',main='',xlab='Mean',ylab='1/sqrt(CV)')


###################################################
### chunk number 16: fig4a
###################################################
checkfit(gg1,x=x[,c(-6,-12)],groups,type='data',main='')


###################################################
### chunk number 17: fig4b
###################################################
checkfit(gg1,x=x[,c(-6,-12)],groups,type='shape',main='')


###################################################
### chunk number 18: fig4c
###################################################
checkfit(gg1,x=x[,c(-6,-12)],groups,type='mean',main='')


###################################################
### chunk number 19: fig4d
###################################################
checkfit(gg1,x=x[,c(-6,-12)],groups,type='shapemean',main='',xlab='Mean',ylab='1/sqrt(CV)')


###################################################
### chunk number 20: eight
###################################################
d1 <- findgenes(gg1,x[,c(-6,-12)],groups,fdrmax=.05,parametric=TRUE)
d1.nonpar <- findgenes(gg1,x[,c(-6,-12)],groups,fdrmax=.05,parametric=FALSE,B=1000)
dtrue <- (l[,1]!=l[,2])
table(d1$d,dtrue)
table(d1.nonpar$d,dtrue)


###################################################
### chunk number 21: fig5plot
###################################################
plot(d1.nonpar$fdrest,type='l',xlab='Bayesian FDR',ylab='Estimated frequentist FDR')


###################################################
### chunk number 22: fig1
###################################################
plot(d1.nonpar$fdrest,type='l',xlab='Bayesian FDR',ylab='Estimated frequentist FDR')


###################################################
### chunk number 23: eightbis
###################################################
d2 <- findgenes(gg2,x[,c(-6,-12)],groups,fdrmax=.05,parametric=TRUE)
d3 <- findgenes(gg3,x[,c(-6,-12)],groups,fdrmax=.05,parametric=TRUE)
table(d1$d,d2$d)
table(d1$d,d3$d)


###################################################
### chunk number 24: fc1
###################################################
mpos <- posmeansGG(gg1,x=x[,c(-6,-12)],groups=groups,underpattern=1)
fc <- mpos[,1]-mpos[,2]
fc[1:5]


###################################################
### chunk number 25: nine
###################################################
pred1 <- classpred(gg1,xnew=x[,6],x=x[,c(-6,-12)],groups,ngene=50,prgroups=c(.5,.5))
pred2 <- classpred(gg1,xnew=x[,12],x=x[,c(-6,-12)],groups,ngene=50,prgroups=c(.5,.5))
pred1
pred2