################################################### ### chunk number 1: loadlibs ###################################################

library(golubEsets) library(e1071) library(Biobase) library(genefilter)

################################################### ### chunk number 2: windsorize ################################################### X <- exprs(golubTrain) Wlow <- 100 Whigh <- 16000 X[XWhigh] <- Whigh

################################################### ### chunk number 3: filter ################################################### mmfilt <- function(r=5, d=500, na.rm=TRUE) { function(x) { minval <- min(x, na.rm=na.rm) maxval <- max(x, na.rm=na.rm) (maxval/minval > r) && (maxval-minval > d) } }

mmfun <- mmfilt()

ffun <- filterfun(mmfun) sub <- genefilter(X, ffun ) sum(sub) ## Should get 3051

################################################### ### chunk number 4: subset ################################################### X <- X[sub,] dim(X) golubTrainSub<-golubTrain[sub,] golubTrainSub@exprs <- X

Y <- golubTrainSub$ALL.AML

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

################################################### ### chunk number 5: testset ################################################### Xt <- exprs(golubTest) Xt[XtWhigh] <- Whigh

golubTestSub <- golubTest[sub,]

golubTestSub@exprs <- Xt[sub,]

Ytest <- golubTestSub$ALL.AML

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

################################################### ### chunk number 6: fitsvm ################################################### Xm <- t(exprs(golubTrainSub)) resp <- golubTrainSub$ALL svm1 <- svm(Xm, resp, type="C-classification", kernel="linear")

################################################### ### chunk number 7: training-error ###################################################

trpred <- predict(svm1, Xm) sum( trpred != resp) table(trpred, resp)

################################################### ### chunk number 8: trcv ###################################################

trcv <- svm(Xm, resp, type="C-classification", kernel="linear", cross=10)

summary(trcv)

################################################### ### chunk number 9: comptr ###################################################

Xmtr <- t(exprs(golubTestSub)) tepred <- predict(svm1, Xmtr) sum(tepred != golubTestSub$ALL) table(tepred, golubTestSub$ALL)

################################################### ### chunk number 10: exp1 ###################################################

newlabs <- sample(c(rep("A",27), rep("B",11)), 38) table(newlabs, golubTrain$ALL) funnysvm <- svm(Xm, newlabs, type="C-classification", kernel="linear")

fpred <- predict(funnysvm, Xm) sum( fpred != newlabs) table(fpred, newlabs)

################################################### ### chunk number 11: cv2 ###################################################

trfcv <- svm(Xm, newlabs, type="C-classification", kernel="linear", cross=10)

summary(trfcv)

################################################### ### chunk number 12: loadlibs ###################################################

library(randomForest)

################################################### ### chunk number 13: rf1 ###################################################

set.seed(123) rf1 <- randomForest(t(X), golubTrainSub$ALL, ntree=2000, mtry=55, importance=TRUE) rf1

rf2 <- randomForest(t(X), golubTrainSub$ALL, ntree=2000, mtry=35, importance=TRUE)

rf2

################################################### ### chunk number 14: rfpred ###################################################

p1 <- predict(rf1, t(Xt), prox=TRUE) table(p1$pred) p2 <- predict(rf2, t(Xt), prox=TRUE) table(p2$pred)

################################################### ### chunk number 15: evenup ###################################################

gT2 <- golubTestSub[,7:34] rfx <- randomForest(t(exprs(gT2)), gT2$ALL, ntree=2000, mtry=55, importance=TRUE) rfx

rfx2 <- randomForest(t(exprs(gT2)), gT2$ALL, ntree=2000, mtry=35, importance=TRUE) rfx2

px2 <- predict(rfx2, t(X))

table(px2, golubTrainSub$ALL)

################################################### ### chunk number 16: varimp ###################################################

var.imp.plot(rf1, n.var=15) var.imp.plot(rf2, n.var=15)

impvars <- function(x, which=2, k=10) { v1<-order(x$importance[,which]) l1 <- length(v1) x$importance[v1[(l1-k+1):l1], which] } iv.rf1 <- impvars(rf1, k=25) library(hu6800) library(annotate) isyms <- getSYMBOL(names(iv.rf1),data="hu6800")