###################################################
### chunk number 1: loadData
###################################################
library(IBCLab4)
data(GolubData) 


###################################################
### chunk number 2: hclustCor
###################################################
library(stats) 
clust.cor <- hclust(as.dist(1-cor(golub)), method="complete") 
plot(clust.cor, cex=0.6) 


###################################################
### chunk number 3: hclustEuclid
###################################################
clust.euclid <- hclust(dist(t(golub)), method="average") 
plot(clust.euclid, cex=0.6) 


###################################################
### chunk number 4: heatmap
###################################################
library(sma)
golubvar <- apply(golub, 1, var, na.rm=TRUE) 
top100 <- stat.gnames(golubvar, 1:length(golubvar), crit=100)$gnames 
heatmap(golub[top100,]) 


###################################################
### chunk number 5: kmeans
###################################################
clust.kmeans <- kmeans(as.dist(1 - cor(golub)), 3) 
names(clust.kmeans$cluster) <- colnames(golub) 
clust.kmeans$cluster[1:10] 


###################################################
### chunk number 6: pam
###################################################
library(cluster) 
clust.pam <- pam(as.dist(1 - cor(golub)), 3, diss=TRUE) 
clusplot(clust.pam, labels=3, col.p=clust.pam$clustering) 


###################################################
### chunk number 7: GeneSelection
###################################################
golubSub <- golub[top100,] 
par(mfrow=c(2,2)) 
plot(hclust(as.dist(1-cor(golubSub)), method="complete"), cex=0.6) 
plot(hclust(dist(t(golubSub)), method="average"), cex=0.6) 
clust.pam <- pam(as.dist(1 - cor(golubSub)), 3, diss=TRUE) 
clusplot(clust.pam, labels=3, col.p=clust.pam$clustering) 
par(mfrow=c(1,1)) 


###################################################
### chunk number 8: knn
###################################################
library(class) 
LS.bw <- stat.bwss(t(LS), as.integer(LS.resp))$bw 
LStop100 <- stat.gnames(LS.bw, 1:length(LS.bw), crit=100)$gnames 
for(k in 4:7) 
{ 
  disc.knn<- as.vector(knn(LS[,LStop100], TS[, LStop100], cl=LS.resp, k)) 
  print(paste("K = ", k)) 
  print(table(disc.knn, TS.resp)) 
} 


###################################################
### chunk number 9: classif
###################################################
newLS <- LS[,LStop100[2:3]] 
newTS <- geneTS(newLS) 
res.knn <- knn(newLS, newTS, cl=as.numeric(LS.resp), k=3) 
plot.class2(newLS, newTS, res.knn, cl=LS.resp) 


###################################################
### chunk number 10: dQda
###################################################
disc.dqda <- stat.diag.da(LS[,LStop100], TS[,LStop100], cl=as.integer(LS.resp), pool=0)[[1]] 
table(disc.dqda, TS.resp) 
 
disc.dlda <- stat.diag.da(LS[,LStop100], TS[,LStop100], cl=as.integer(LS.resp), pool=1)[[1]] 
table(disc.dlda, TS.resp) 

res.dlda <- stat.diag.da(newLS, newTS, cl=as.integer(LS.resp), pool=1)$pred
plot.class2(newLS, newTS, res.dlda, cl=LS.resp)


###################################################
### chunk number 11: tree
###################################################
library(rpart) 
disc.rpart <- dorpart(LS=LS[,LStop100], TS=TS[,LStop100], cl=as.numeric(LS.resp)) 
table(disc.rpart, TS.resp) 


###################################################
### chunk number 12: baggtree
###################################################
library(ipred)
fit.bagtree <-  ipredbagg.factor(y=as.factor(LS.resp), X=as.data.frame(LS[,LStop100]), nbagg=100, coob=TRUE)
pred.bagtree <- predict(fit.bagtree, newdata= as.data.frame(TS[,LStop100]))
table(pred.bagtree, TS.resp)


###################################################
### chunk number 13: er1
###################################################
pred.LS <- as.vector(knn(LS, LS, cl=LS.resp, k=3)) 
table(pred.LS, LS.resp) 
sum(pred.LS != LS.resp) 


###################################################
### chunk number 14: er2
###################################################
pred.TS <- as.vector(knn(LS, TS, cl=LS.resp, k=3)) 
table(pred.TS, TS.resp) 
sum(pred.TS != TS.resp) 


###################################################
### chunk number 15: er3
###################################################
LStop5<- stat.gnames(LS.bw, 1:length(LS.bw), crit=5)$gnames 
pred.LS <- as.vector(knn(LS[,LStop5], LS[, LStop5], cl=LS.resp, k=2)) 
table(pred.LS, LS.resp) 
sum(pred.LS != LS.resp) 
pred.TS <- as.vector(knn(LS[,LStop5], TS[, LStop5], cl=LS.resp, k=2)) 
table(pred.TS, TS.resp) 
sum(pred.TS != TS.resp) 


###################################################
### chunk number 16: kcv
###################################################
knn.cvk<-function(LS, cl, k=1:5, ...) 
{ 
  cv.err <- cl.pred <- c() 
  for(i in k) 
  { 
    newpre <- as.vector(knn.cv(LS, cl, i)) 
    cl.pred <- cbind(cl.pred, newpre) 
    cv.err <- c(cv.err, sum(cl!=newpre)) 
  } 
  k0<-k[which.min(cv.err)] 
  return(list(k=k0,pred=cl.pred[,which.min(cv.err)])) 
} 
knn.cvk(LS, cl=LS.resp, k=3:7) 


###################################################
### chunk number 17: Rosetta1
###################################################
data(RosettaBreastData)
rindex.train <- 1:78
rindex.test <- 79:97