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

library(ALL) library(e1071) library(Biobase) library(genefilter)

################################################### ### chunk number 2: subset ################################################### data("ALL") ## subset of interest: 37+42 samples sel1 <- grep("^B", as.character(ALL$BT)) sel2 <- which(as.character(ALL$mol) %in% c("BCR/ABL","NEG")) eset <- ALL[, intersect(sel1, sel2)]

################################################### ### chunk number 3: genefilter ################################################### library(genefilter) ## intensities above 200 in at least 25% of the samples ##just so we have fewer genes and it runs in finite time f1 <- pOverA(.25, log2(200)) f2 <- function(x)(IQR(x)>0.5) ff <- filterfun(f1, f2) selected <- genefilter(eset, ff) sum(selected) esetSub <- eset[selected,]

################################################### ### chunk number 4: ################################################### Y <- factor(esetSub$mol.biol) #$ table(Y) X <- t(exprs(esetSub))

################################################### ### chunk number 5: tree1 ################################################### library(rpart) df <- data.frame(Y=Y, X=X) tr <- rpart(Y ~ ., data=df)

################################################### ### chunk number 6: ################################################### plotcp(tr)

################################################### ### chunk number 7: treesum ################################################### print(tr) print(summary(tr))

################################################### ### chunk number 8: testtrain ################################################### a1=which(as.character(Y)=="NEG") a2 = which(as.character(Y)=="BCR/ABL") set.seed(100) sub1N = sample(a1, 21) sub1B = sample(a2, 19) sub2N = a1[!(a1 %in% sub1N)] sub2B = a2[!(a2 %in% sub1B)] testSet = esetSub[, c(sub1N, sub1B)] testY = Y[c(sub1N, sub1B)] trainSet = esetSub[, c(sub2N, sub2B)] trainY = Y[c(sub2N, sub2B)]

################################################### ### chunk number 9: nn1 ################################################### library(nnet) df <- data.frame(Y=Y, X=X[, sample(ncol(X), 150)]) nn1 <- nnet(Y~., data=df, size=5, decay=.01, MaxNWts=5000)

################################################### ### chunk number 10: nn.print ################################################### print(nn1) print(table(predict(nn1,type="class"), Y))

################################################### ### chunk number 11: ################################################### library(class) knn1 <- knn(t(exprs(trainSet)), t(exprs(testSet)), trainY, k=1) table(knn1, testY)

################################################### ### chunk number 12: fitsvm ################################################### Xm <- t(exprs(trainSet)) svm1 <- svm(Xm, trainY, type="C-classification", kernel="linear")

################################################### ### chunk number 13: training-error ###################################################

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

################################################### ### chunk number 14: trcv ###################################################

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

summary(trcv)

################################################### ### chunk number 15: comptr ###################################################

Xmtr <- t(exprs(testSet)) tepred <- predict(svm1, Xmtr) sum(tepred != testY) table(tepred, testY)

################################################### ### chunk number 16: exp1 ###################################################

newlabs <- sample(c(rep("B", 18), rep("N",21)), 39) table(newlabs, trainY) funnysvm <- svm(Xm, newlabs, type="C-classification", kernel="linear")

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

################################################### ### chunk number 17: cv2 ###################################################

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

summary(trfcv)

################################################### ### chunk number 18: loadlibs ###################################################

library(randomForest)

################################################### ### chunk number 19: rf1 ###################################################

set.seed(123) rf1 <- randomForest(Xm, trainY, ntree=2000, mtry=55, importance=TRUE) rf1

rf2 <- randomForest(Xm, trainY, ntree=2000, mtry=35, importance=TRUE)

rf2

################################################### ### chunk number 20: rfpred ###################################################

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

################################################### ### chunk number 21: ################################################### var.imp.plot(rf1, n.var=15)

################################################### ### chunk number 22: ################################################### var.imp.plot(rf2, n.var=15)

################################################### ### chunk number 23: varimp ################################################### 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("hgu95av2") library(annotate) isyms <- getSYMBOL(names(iv.rf1),data="hgu95av2")