attach("JSMLab6.RData")

index.train <- 1:78

index.test <- 79:97

trainData <- bcdata[,index.train]

testData <- bcdata[,index.test]



#############################################

## Q1.

#############################################

library(stats)

clust.cor <- hclust(as.dist(1 - cor(bcdata)), method = "complete")

plot(clust.cor, cex = 0.6)





#############################################

## Q2.

#############################################

clust.euclid <- hclust(dist(t(bcdata)), method = "average")

plot(clust.euclid, cex = 0.6)





#############################################

## Q3.

#############################################

library(sma)

bcdatavar <- apply(bcdata, 1, var, na.rm = TRUE)

top100 <- stat.gnames(bcdatavar, 1:length(bcdatavar), crit = 100)$gnames

heatmap(bcdata[top100, ])





#############################################

## Q4.

#############################################

clust.kmeans <- kmeans(as.dist(1-cor(bcdata)),2)

names(clust.kmeans$cluster) <- colnames(bcdata)

clust.kmeans$cluster[1:10]





#############################################

## Q5.

#############################################

library(cluster)

clust.pam <- pam(as.dist(1-cor(bcdata)), 2, diss =TRUE)

clusplot(clust.pam, labels = 2, col.p = clust.pam$clustering)



#############################################

## Q7.

#############################################



golubSub <- bcdata[top100, ]

par(mfrow = c(2, 2))

plot(hclust(as.dist(1 - cor(bcdata)), method = "complete"), cex = 0.6)

plot(hclust(dist(t(bcdata)), method = "average"), cex = 0.6)

clust.pam <- pam(as.dist(1 - cor(bcdata)), 2, diss = TRUE)

clusplot(clust.pam, labels = 2, col.p = clust.pam$clustering)

par(mfrow = c(1, 1))





#############################################

## Q8.

#############################################



library(sma)

rosetta.mad <- apply(bcdata[, index.train], 1, mad, na.rm = TRUE)

top100.mad <- stat.gnames(rosetta.mad, 1:length(rosetta.mad), crit = 100)$gnames





#############################################

## Q9.

#############################################



rosetta.tpool <- apply(bcdata[, index.train], 1, t.func, surv.resp[index.train], var.equal = TRUE)

top100.tpool <- stat.gnames(abs(rosetta.tpool), 1:length(rosetta.tpool), crit = 100)$gnames





#############################################

## Q10.

#############################################



labcols <- c(rep("green", length(which(surv.lab[index.train] == "good"))), 

             rep("red", length(which(surv.lab[index.train] == "bad"))))

par(pty = "s")



dt <- plotCorOrder()

plot.cor(cor.na(dt), new = F, nrgcols = 50, labels = colnames(trainData), labcols = labcols, 

         title = "All genes", zlim = c(-1,1))



dt <- plotCorOrder(top100.mad)

plot.cor(cor.na(dt), new = F, nrgcols = 50, labels = colnames(trainData), labcols = labcols, 

         title = "100 most variable genes", zlim = c(-1, 1))



dt <- plotCorOrder(top100.tpool)

plot.cor(cor.na(dt), new = F, nrgcols = 50, labels = colnames(trainData), labcols = labcols, 

         title = "Top 100 genes with the highest association",  zlim = c(-1, 1))





############################################# 

## 

############################################# 



LS <- t(bcdata[top100.tpool, index.train])

TS <- t(bcdata[top100.tpool, index.test])

LS.resp <- surv.resp[index.train] 

TS.resp <- surv.resp[index.test]





#############################################

## Q11.

#############################################



library(class)

for (k in 4:7) {

  disc.knn <- as.vector(knn(LS, TS, cl = LS.resp, k))

  print(paste("K = ", k))

  print(table(disc.knn, TS.resp)); cat("\n\n")

}





#############################################

## Q12.

#############################################



plotLS <- LS[, 1:2]

plotTS <- geneTS(plotLS)

res.knn <- knn(plotLS, plotTS, cl = LS.resp, k = 2)

plot.class2(plotLS, plotTS, as.vector(res.knn), cl = LS.resp)





#############################################

## Q13.

#############################################



disc.dqda <- stat.diag.da(LS, TS, cl = LS.resp, pool = 0)[[1]]

table(disc.dqda, TS.resp)



disc.dlda <- stat.diag.da(LS, TS, cl = LS.resp, pool = 1)[[1]]

table(disc.dlda, TS.resp)



res.dlda <- stat.diag.da(plotLS, plotTS, cl=LS.resp, pool=1)[[1]]

plot.class2(plotLS, plotTS, res.dlda, cl=LS.resp)





#############################################

## Q14.

#############################################



fit.logit <- glm(LS.resp ~ ., data = as.data.frame(LS), family = "binomial")

pred.logit <- predict(fit.logit, newdata = as.data.frame(TS), type = "response")



## converting output probabilities into classes.

pred.logit.type <- rep(0, length(pred.logit))

pred.logit.type[pred.logit >= 0.5] <- 1

pred.logit.type[pred.logit < 0.5] <- 0

table(pred.logit.type, TS.resp)



#############################################

## Q15.

#############################################



library(rpart)

disc.rpart <- dorpart(LS, TS, cl = as.numeric(LS.resp))

table(disc.rpart, TS.resp)





#############################################

## Q16.

#############################################



library(ipred)

fit.bagging <- bagging(as.factor(LS.resp) ~ ., data = as.data.frame(LS), nbagg = 100, coob = TRUE)

pred.bagging <- predict(fit.bagging, newdata = as.data.frame(TS))

table(pred.bagging, TS.resp)



## fit.bagging <- bagging(as.factor(LS.resp) ~ ., data = as.data.frame(plotLS), nbagg = 100, coob = TRUE)

## res.bagging <- predict(fit.bagging, plotTS, cl = as.integer(LS.resp))

## plot.class2(plotLS, plotTS, res.bagging, cl = LS.resp)





#############################################

## Q17.

#############################################



pred.LS <- as.vector(knn(LS, LS, cl = LS.resp, k = 3))

table(pred.LS, LS.resp)

sum(pred.LS != LS.resp)



pred.TS <- as.vector(knn(LS, TS, cl = LS.resp, k = 3))

table(pred.TS, TS.resp)

sum(pred.TS != TS.resp)



LStop5 <- stat.gnames(top100.tpool, 1:length(top100.tpool), 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)



#############################################

## Q18.

#############################################



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)