#!/usr/bin/Rscript library(Biobase) library(ALL) library(hgu95av2.db) library(pamr) data(ALL) samples_b <- which(substr(ALL$BT,1,1)=="B") samples_bcr.abl.neg <- which(ALL$mol.biol %in% c("BCR/ABL", "NEG")) subset <- intersect(samples_b, samples_bcr.abl.neg) subset.data <- ALL[,subset] subset.data$label <- factor(subset.data$mol.biol) test <- sample(1:length(subset),6) testset.data <- subset.data[,test] trainset.data <- subset.data[,-test] e <- exprs(trainset.data) genes <- featureNames(trainset.data) #genenames <- sapply(genes, get, hgu95av2SYMBOL) train.data <- list(x = e, y = trainset.data$label, genenames = genenames, geneid = genes) model <- pamr.train(train.data) print(model) readline() model.cv <- pamr.cv(model, train.data, nfold = 10) pamr.plotcv(model.cv) readline() par(mfrow=c(3,3)) for(Delta in seq(1,5,.5)) { pamr.plotcvprob(model, train.data, Delta) p <- pamr.predict(model, exprs(testset.data), Delta, type="posterior") d <- abs(p[,1]-p[,2]) print(names(sort(d))) print(mean(d)) } readline() for(i in 1:9) { data <- rnorm(10000*11, 0, 1) dim(data) <- c(10000, 11) data[,1] <- data[,1]>0 train.data <- list(x = data[,2:11], y = data[,1]) model <- pamr.train(train.data) model.cv <- pamr.cv(model, train.data, nfold = 10) print(model) print(model.cv) pamr.plotcv(model.cv) readline() }