use*_*694 5 r classification svm bioconductor
我正在使用Bioconductor软件包CMA在微阵列数据集中对SVM分类器执行内部蒙特卡罗交叉验证(MCCV).CMA内部使用e1071 R软件包进行SVM工作.
数据集有45个样本(观察)的387个变量(属性),它们属于两个类别之一(标签0或1;比例约为1:1).所有数据都是数字的,没有NA.我正在尝试使用limma统计数据进行差异基因表达分析,为SVM选择15个变量的1000次迭代MCCV .在MCCV期间,45个样本集的一部分用于训练SVM分类器,然后用于测试剩余分数,并且我正在尝试训练集分数的不同值.CMA还执行内循环验证(默认情况下在训练集内进行3次交叉验证),以微调用于针对测试集进行交叉验证的分类器.所有这些都是在CMA包中完成的.
有时,对于低训练集大小,CMA在控制台中显示错误并停止执行分类的其余代码.
[snip]tuning iteration 575 tuning iteration 576 tuning iteration 577 Error in predict.svm(ret, xhold, decision.values = TRUE) : Model is empty!
即使我使用除limma之外的测试进行变量选择,或者使用两个而不是15个变量进行分类器生成,它也会发生.我使用的R代码应确保训练集始终具有两个类的成员.我很感激任何见解.
下面是我使用的R代码,Mac OS X 10.6.6,R 2.12.1,Biobase 2.10.0,CMA 1.8.1,limma 3.6.9和WilcoxCV 1.0.2.数据文件hy3ExpHsaMir.txt可以从http://rapidshare.com/files/447062901/hy3ExpHsaMir.txt下载.
一切顺利,直到g在for(g在0:10)循环中为9 (用于改变训练/测试集大小).
# exp is the expression table, a matrix; 'classes' is list of known classes
exp <- as.matrix(read.table(file='hy3ExpHsaMir.txt', sep='\t', row.names=1, header=T, check.names=F))
#best is to use 0 and 1 as class labels (instead of 'p', 'g', etc.) with 1 for 'positive' items (positive for recurrence, or for disease, etc.)
classes <- c(1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 1, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0)
yesPredVal = 1 # class label for 'positive' items in 'classes'
library(CMA)
library(WilcoxCV)
myNumFun <- function(x, y){round(y(as.numeric(x), na.rm=T), 4)}
set.seed(631)
out = ''
out2 = '\nEffect of varying the training-set size:\nTraining-set size\tSuccessful iterations\tMean acc.\tSD acc.\tMean sens.\tSD sens.\tMean spec.\tSD spec.\tMean PPV\tSD PPV\tMean NPV\tSD NPV\tTotal genes in the classifiers\n'
niter = 1000
diffTest = 'limma'
diffGeneNum = 15
svmCost <- c(0.1, 0.2, 0.5, 1, 2, 5, 10, 20, 50)
for(g in 0:10){ # varying the training/test-set sizes
ntest = 3+g*3 # test-set size
result <- matrix(nrow=0, ncol=7)
colnames(result) <- c('trainSetSize', 'iteration', 'acc', 'sens', 'spec', 'ppv', 'npv')
diffGenes <- numeric()
# generate training and test sets
lsets <- GenerateLearningsets(n=ncol(exp), y=classes, method=c('MCCV'), niter=niter, ntrain=ncol(exp)-ntest)
# actual prediction work
svm <- classification(t(exp), factor(classes), learningsets=lsets, genesellist= list(method=diffTest), classifier=svmCMA, nbgene= diffGeneNum, tuninglist=list(grids=list(cost=svmCost)), probability=TRUE)
svm <- join(svm)
# genes in classifiers
svmGenes <- GeneSelection(t(exp), classes, learningsets=lsets, method=diffTest)
actualIters=0
for(h in 1:niter){
m <- ntest*(h-1)
# valid SVM classification requires min. 2 classes
if(1 < length(unique(classes[-lsets@learnmatrix[h,]]))){
actualIters = actualIters+1
tp <- tn <- fp <- fn <- 0
for(i in 1:ntest){
pred <- svm@yhat[m+i]
known <- svm@y[m+i]
if(pred == known){
if(pred == yesPredVal){tp <- tp+1}
else{tn <- tn+1}
}else{
if(pred == yesPredVal){fp <- fp+1}
else{fn <- fn+1}
}
}
result <- rbind(result, c(ncol(exp)-ntest, h, (tp+tn)/(tp+tn+fp+fn), tp/(tp+fn), tn/(tn+fp), tp/(tp+fp), tn/(tn+fn)))
diffGenes <- c(diffGenes, toplist(svmGenes, k=diffGeneNum, iter=h, show=F)$index)
} # end if valid SVM
} # end for h
# output accuracy, etc.
out = paste(out, 'SVM MCCV using ', niter, ' attempted iterations and ', actualIters, ' successful iterations, with ', ncol(exp)-ntest, ' of ', ncol(exp), ' total samples used for training:\nThe means (ranges; SDs) of prediction accuracy, sensitivity, specificity, PPV and NPV in fractions are ',
myNumFun(result[, 'acc'],mean), ' (', myNumFun(result[, 'acc'], min), '-', myNumFun(result[, 'acc'], max), '; ', myNumFun(result[, 'acc'], sd), '), ',
myNumFun(result[, 'sens'], mean), ' (', myNumFun(result[, 'sens'], min), '-', myNumFun(result[, 'sens'], max), '; ', myNumFun(result[, 'sens'], sd), '), ',
myNumFun(result[, 'spec'], mean), ' (', myNumFun(result[, 'spec'], min), '-', myNumFun(result[, 'spec'], max), '; ', myNumFun(result[, 'spec'], sd), '), ',
myNumFun(result[, 'ppv'], mean), ' (', myNumFun(result[, 'ppv'], min), '-', myNumFun(result[, 'ppv'], max), '; ', myNumFun(result[, 'ppv'], sd), '), and ',
myNumFun(result[, 'npv'], mean), ' (', myNumFun(result[, 'npv'], min), '-', myNumFun(result[, 'npv'], max), '; ', myNumFun(result[, 'npv'], sd), '), respectively.\n', sep='')
# output classifier genes
diffGenesUnq <- unique(diffGenes)
out = paste(out, 'A total of ', length(diffGenesUnq), ' genes occur in the ', actualIters, ' classifiers, with occurrence frequencies in fractions:\n', sep='')
for(i in 1:length(diffGenesUnq)){
out = paste(out, rownames(exp)[diffGenesUnq[i]], '\t', round(sum(diffGenes == diffGenesUnq[i])/actualIters, 3), '\n', sep='')
}
# output split-size effect
out2 = paste(out2, ncol(exp)-ntest, '\t', actualIters, '\t', myNumFun(result[, 'acc'], mean), '\t', myNumFun(result[, 'acc'], sd), '\t', myNumFun(result[, 'sens'], mean), '\t', myNumFun(result[, 'sens'], sd), '\t', myNumFun(result[, 'spec'], mean), '\t', myNumFun(result[, 'spec'], sd), '\t', myNumFun(result[, 'ppv'], mean), '\t', myNumFun(result[, 'ppv'], sd),
'\t', myNumFun(result[, 'npv'], mean), '\t', myNumFun(result[, 'npv'], sd), '\t', length(diffGenesUnq), '\n', sep='')
} # end for g
cat(out, out2, sep='')traceback()的输出:
20: stop("Model is empty!")
19: predict.svm(ret, xhold, decision.values = TRUE)
18: predict(ret, xhold, decision.values = TRUE)
17: na.action(predict(ret, xhold, decision.values = TRUE))
16: svm.default(cost = 0.1, kernel = "linear", type = "C-classification", ...
15: svm(cost = 0.1, kernel = "linear", type = "C-classification", ...
14: do.call("svm", args = ll)
13: function (X, y, f, learnind, probability, models = FALSE, ...) ...
12: function (X, y, f, learnind, probability, models = FALSE, ...) ...
11: do.call(classifier, args = c(list(X = X, y = y, learnind = learnmatrix[i, ...
10: classification(X = c(83.5832768669369, 83.146333099001, 94.253534443549, ...
9: classification(X = c(83.5832768669369, 83.146333099001, 94.253534443549, ...
8: do.call("classification", args = c(list(X = Xi, y = yi, learningsets = lsi, ...
7: tune(grids = list(cost = c(0.1, 0.2, 0.5, 1, 2, 5, 10, 20, 50...
6: tune(grids = list(cost = c(0.1, 0.2, 0.5, 1, 2, 5, 10, 20, 50...
5: do.call("tune", args = c(tuninglist, ll))
4: classification(X, y = as.numeric(y) - 1, learningsets = learningsets, ...
3: classification(X, y = as.numeric(y) - 1, learningsets = learningsets, ...
2: classification(t(exp), factor(classes), learningsets = lsets, ...
1: classification(t(exp), factor(classes), learningsets = lsets, ...
CMA 包的维护者立即回复了我发送的有关此问题的消息。
CMA 通过在训练集内的 k 倍 CV 步骤（默认 k=3）中测试不同的参数值来调整从训练集生成的分类器。对于较小的训练集大小，如果仅对一个类的观察进行子集化，则此内部循环可能会失败。减少发生这种情况的机会的两种方法是执行 2 倍内部 CV，并指定分层采样，这两种方法都要求通过 CMA 的une()单独调用调整步骤并在classification()中使用其输出。
在我发布的代码中，从classification()内部调用调整，它不允许分层采样或2倍CV。然而，对于分层采样和 2 倍 CV 单独调用tunnel()对我的情况没有帮助。这并不奇怪，因为对于小型训练集，CMA 会遇到只有一个类别的一组成员的情况。
我希望 CMA 在一个学习集遇到这样的问题时，不要突然结束一切，而是继续处理剩余的学习集。如果在遇到此问题时，CMA 能够为内部 k 重 CV 步骤尝试不同的 k 值，那就太好了。
[2 月 14 日编辑] CMA 的 CV 学习集生成不会检查训练集中是否存在两个类别的足够成员。因此，以下替换原始帖子中的部分代码应该可以使事情正常工作：
numInnerFold = 3 # k for the k-fold inner validation called through tune()
# generate learning-sets with 2*niter members in case some have to be removed
lsets <- GenerateLearningsets(n=ncol(exp), y=classes, method=c('MCCV'), niter=2*niter, ntrain=ncol(exp)-ntest)
temp <- lsets@learnmatrix
for(i in 1:(2*niter)){
unq <- unique(classes[lsets@learnmatrix[i, ]])
if((2 > length(unique(classes[lsets@learnmatrix[i, ]])))
| (numInnerFold > sum(classes[lsets@learnmatrix[i, ]] == yesPredVal))
| (numInnerFold > sum(classes[lsets@learnmatrix[i, ]] != yesPredVal))){
# cat('removed set', i,'\n')
temp <- lsets@learnmatrix[-i, ]
}
}
lsets@learnmatrix <- temp[1:niter, ]
lsets@iter <- niter
# genes in classifiers
svmGenes <- GeneSelection(t(exp), classes, learningsets=lsets, method=diffTest)
svmTune <- tune(t(exp), factor(classes), learningsets=lsets, genesel=svmGenes, classifier=svmCMA, nbgene=diffGeneNum, grids=list(cost=svmCost), strat=T, fold=numInnerFold)
# actual prediction work
svm <- classification(t(exp), factor(classes), learningsets=lsets, genesel=svmGenes, classifier=svmCMA, nbgene=diffGeneNum, tuneres=svmTune)