################################################### ### chunk number 1: loadlibs ################################################### library(Biobase) library(annotate) library(golubEsets) library(genefilter) library(nnet) library(mva) library(e1071) library(class) library(bioclabs) library(MASS)

################################################### ### chunk number 2: setup ###################################################

mmfilt <- function(r=5, d=500, na.rm=TRUE) { function(x) { minval <- min(2^x, na.rm=na.rm) maxval <- max(2^x, na.rm=na.rm) (maxval/minval > r) && (maxval-minval > d) } }

GolubTrans <-function(eSet) { X<-exprs(eSet) X[X<100]<-100 X[X>16000]<-16000 X <- log2(X) eSet@exprs <- X eSet }

data(golubTrain) data(golubMerge) data(golubTest) gTrn <- GolubTrans(golubTrain) gTest <- GolubTrans(golubTest) gMerge <- GolubTrans(golubMerge)

mmfun <- mmfilt()

ffun <- filterfun(mmfun) sub <- genefilter(gTrn, ffun ) ##three are right on r=5, d=500 -- they are dropped if we ##don't fool around with logs -- so we drop them here sub[c(2401,3398,4168)] <- FALSE sum(sub) ## Should get 3051

## we will subset all data sets according to this same criterion

gTrnS <- gTrn[sub,] gTestS <- gTest[sub,] gMergeS <- gMerge[sub,]

Ytr <-paste(golubTrain$ALL.AML,golubTrain$T.B.cell) Ytr<-sub("NA","",Ytr)

Ytest <- paste(golubTest$ALL.AML, golubTest$T.B.cell) Ytest <- sub("NA","",Ytest)

Ymerge <- paste(golubMerge$ALL.AML, golubMerge$T.B.cell) Ymerge <- sub("NA","",Ymerge)

################################################### ### chunk number 3: Anovagenes ###################################################

data(gfaF)

gTrA <- gTrnS[gfaF,] gTeA <- gTestS[gfaF,] gMeA <- gMergeS[gfaF,]

whBad1 <- (apply(gTrA@exprs, 1, mad) == 0)

whBad2 <- (apply(gTeA@exprs, 1, mad) == 0)

whBad <- whBad1 | whBad2

sum(whBad)

gTrA <- gTrA[!whBad,] gTeA <- gTeA[!whBad,] gMeA <- gMeA[!whBad,]

################################################### ### chunk number 4: removeBad ###################################################

star <- function(x) (x-median(x))/mad(x)

TrExprs <- t(apply(exprs(gTrA), 1, star)) TeExprs <- t(apply(exprs(gTeA), 1, star))

MeExprs <- t(apply(exprs(gMeA), 1, star))

################################################### ### chunk number 5: lda ###################################################

gTr.lda <- lda(t(TrExprs), Ytr)

plot(gTr.lda)

preds.lda <- predict(gTr.lda, t(TeExprs))

table(preds.lda$class, Ytest)

################################################### ### chunk number 6: logisticda ###################################################

library(nnet)

gTr.mult <- multinom(factor(Ytr) ~ ., data=data.frame(t(TrExprs)), maxit=250)

tEdf <- data.frame(t(TeExprs))

log.preds <- predict(gTr.mult, data.frame(t(TeExprs)))

table(log.preds, Ytest)

################################################### ### chunk number 7: knn1 ###################################################

knn1 <- knn(t(TrExprs), t(TeExprs), factor(Ytr), k=1)

table(knn1, Ytest)

################################################### ### chunk number 8: knn3 ###################################################

knn3 <- knn(t(TrExprs), t(TeExprs), factor(Ytr), k=3) table(knn3, Ytest)

################################################### ### chunk number 9: knn5 ###################################################

knn5 <- knn(t(TrExprs), t(TeExprs), factor(Ytr), k=5) table(knn5, Ytest)

################################################### ### chunk number 10: knnCV ################################################### knn1.cvpreds <- knn.cv(t(MeExprs), factor(Ymerge), k=1)

table(knn1.cvpreds, Ymerge)

################################################### ### chunk number 11: knnCV3 ################################################### knn3.cvpreds <- knn.cv(t(MeExprs), factor(Ymerge), k=3)

table(knn3.cvpreds, Ymerge)

################################################### ### chunk number 12: knnCV5 ################################################### knn5.cvpreds <- knn.cv(t(MeExprs), factor(Ymerge), k=5)

table(knn5.cvpreds, Ymerge)