2011-02-08 7 views
5

मैं BioConductor पैकेज CMA उपयोग कर रहा हूँ आंतरिक प्रदर्शन करने के लिए सीएमए BioConductor पैकेज का उपयोग, मोंटे कार्लो पार सत्यापन (MCCV) एक माइक्रोएरे में SVM classifiers पर डाटासेट। सीएमए आंतरिक रूप से एसवीएम काम के लिए e1071 आर पैकेज का उपयोग करता है।SVM वर्गीकरण के लिए पार सत्यापन में एक 'मॉडल खाली' त्रुटि का समाधान जब आर

डेटासेट में 45 नमूने (अवलोकन) के लिए 387 चर (गुण) हैं जो दो वर्गों में से एक हैं (लेबल 0 या 1; लगभग 1: 1 अनुपात में)। सभी डेटा संख्यात्मक के साथ संख्यात्मक है। मैं अलग-अलग जीन अभिव्यक्ति विश्लेषण के लिए limma statistics का उपयोग कर एसवीएम के लिए चुने गए 15 चर के साथ एक 1000-पुनरावृत्ति एमसीसीवी की कोशिश कर रहा हूं। एमसीसीवी के दौरान, 45-नमूना सेट का एक अंश एसवीएम वर्गीकरण को प्रशिक्षित करने के लिए उपयोग किया जाता है, जिसका उपयोग शेष अंश का परीक्षण करने के लिए किया जाता है, और मैं प्रशिक्षण-सेट अंश के लिए अलग-अलग मानों की कोशिश कर रहा हूं। सीएमए परीक्षण-सेट के खिलाफ क्रॉस-सत्यापन के लिए क्लासिफायरों का उपयोग करने के लिए आंतरिक-लूप सत्यापन (प्रशिक्षण सेट के भीतर 3-गुना क्रॉस-सत्यापन) भी करता है। यह सब सीएमए पैकेज के भीतर से किया जाता है।

कभी-कभी, कम प्रशिक्षण-सेट आकारों के लिए, सीएमए कंसोल में एक त्रुटि दिखाता है और वर्गीकरण के लिए शेष कोड को निष्पादित करने से रोकता है।

[snip]tuning iteration 575 
tuning iteration 576 
tuning iteration 577 
Error in predict.svm(ret, xhold, decision.values = TRUE) : Model is empty! 

जब मैं limma चर चयन के लिए की तुलना में एक परीक्षण अन्य उपयोग करें, या वर्गीकारक पीढ़ी के लिए दो 15 के बजाय चर का उपयोग यह भी होता है। मेरे द्वारा उपयोग किए जाने वाले आर कोड को यह सुनिश्चित करना चाहिए कि प्रशिक्षण-सेट में हमेशा दोनों वर्गों के सदस्य हों। मैं इस पर किसी अंतर्दृष्टि की सराहना करता हूं।

मैक ओएस एक्स 10.6.6, आर 2.12.1, बायोबेज 2.10.0, सीएमए 1.8.1, लिमा 3.6.9, और विल्कोक्ससीवी 1.0.2 के साथ मैं आर कोड का उपयोग करता हूं। डेटा फ़ाइल hy3ExpHsaMir.txt को http://rapidshare.com/files/447062901/hy3ExpHsaMir.txt से डाउनलोड किया जा सकता है।

सब कुछ जी (00: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[[email protected][h,]]))){ 
    actualIters = actualIters+1 
    tp <- tn <- fp <- fn <- 0 
    for(i in 1:ntest){ 
    pred <- [email protected][m+i] 
    known <- [email protected][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='')

आउटपुट():

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, ...
+0

डेटा के बिना परीक्षण करना असंभव है। –

+0

यह कुछ ऐसा हो सकता है जिसे आपको पैकेज लेखक के साथ चर्चा करने का प्रयास करना चाहिए। –

+0

मैंने मूल पोस्ट में डेटा फ़ाइल के लिए एक लिंक जोड़ा है। – user594694

उत्तर

3

सीएमए पैकेज के मेंटेनर तुरंत संदेश मैं इस मुद्दे के बारे भेजा था करने के लिए प्रतिक्रिया व्यक्त की।

सीएमए प्रशिक्षण-सेट से उत्पन्न एक वर्गीकरण को एक प्रशिक्षण-सेट, के-फ़ोल्ड सीवी चरण (डिफ़ॉल्ट के = 3) में विभिन्न पैरामीटर मानों का परीक्षण करके ट्यून करता है। छोटे प्रशिक्षण-सेट आकारों के साथ, यह आंतरिक लूप विफल हो सकता है यदि केवल एक वर्ग के अवलोकन उप-सेट हो जाते हैं। इस घटना के मौके को कम करने के दो तरीके 2 गुना आंतरिक सीवी करना है, और स्तरीकृत नमूनाकरण निर्दिष्ट करना है, जिनमें से दोनों को सीएमए के ट्यून() के माध्यम से अलग-अलग ट्यूनिंग चरण लागू किया जाता है और वर्गीकरण में इसका आउटपुट उपयोग किया जाता है।()

मेरे द्वारा पोस्ट किए गए कोड में, ट्यूनिंग वर्गीकरण() के भीतर से लागू की जाती है, जो स्तरीकृत नमूनाकरण या 2-गुना सीवी के लिए अनुमति नहीं देती है। हालांकि, ट्यून() अलग से, स्तरीकृत नमूनाकरण और 2 गुना सीवी के लिए, मेरे मामले में मदद नहीं की। यह आश्चर्यजनक नहीं है, छोटे प्रशिक्षण-सेट के साथ, सीएमए केवल एक वर्ग के सदस्यों के सेट के साथ मामलों का सामना करता है।

मेरी इच्छा है कि अचानक सब कुछ खत्म करने की बजाय, सीएमए शेष सीखने के सेट पर काम करना जारी रखेगा जब यह एक सीखने के सेट के साथ इस तरह के मुद्दे का सामना करेगा।यह भी अच्छा होगा अगर, इस मुद्दे को होने पर, सीएमए आंतरिक के-गुना सीवी चरण के लिए के विभिन्न मानों को आजमाएगा।

[14 फरवरी को संपादित] सीएमए की सीवी के लिए सीखने के सेट की पीढ़ी यह जांच नहीं करती है कि प्रशिक्षण कक्षा में दोनों वर्गों के पर्याप्त सदस्य मौजूद हैं या नहीं। मूल पोस्ट में कोड के एक हिस्से के लिए निम्नलिखित प्रतिस्थापन इसलिए चीजों को सही तरीके से काम करना चाहिए:


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 <- [email protected] 
for(i in 1:(2*niter)){ 
unq <- unique(classes[[email protected][i, ]]) 
if((2 > length(unique(classes[[email protected][i, ]]))) 
    | (numInnerFold > sum(classes[[email protected][i, ]] == yesPredVal)) 
    | (numInnerFold > sum(classes[[email protected][i, ]] != yesPredVal))){ 
    # cat('removed set', i,'\n') 
    temp <- [email protected][-i, ] 
} 
} 
[email protected] <- temp[1:niter, ] 
[email protected] <- 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) 
संबंधित मुद्दे