Risoluzione di un errore "modello vuoto" nella convalida incrociata per la classificazione SVM quando si utilizza il pacchetto CMA Bioconductor per R

Question

Sto utilizzando il pacchetto Bioconductor CMA per eseguire la convalida incrociata Monte Carlo interna (MCCV) sui classificatori SVM in un microarray set di dati. CMA utilizza internamente il pacchetto e1071 R per il lavoro SVM.

Il set di dati ha 387 variabili (attributi) per 45 campioni (osservazioni) che appartengono a una delle due classi (etichette 0 o 1, con proporzione di circa 1: 1). Tutti i dati sono numerici senza NA. Sto provando un MCCV a 1000 iterazioni con 15 variabili selezionate per SVM usando lo limma statistics per l'analisi di espressione genica differenziale. Durante MCCV, una frazione del set di 45 campioni viene utilizzata per addestrare un classificatore SVM, che viene quindi utilizzato per testare la frazione rimanente e sto provando diversi valori per la frazione dell'insieme di allenamento. CMA esegue anche convalide del ciclo interno (3 volte la convalida incrociata all'interno dei set di allenamento, per impostazione predefinita) per mettere a punto i classificatori da utilizzare per la convalida incrociata con i set di test. Tutto questo viene fatto all'interno del pacchetto CMA.

A volte, per le dimensioni del set di allenamento basso, CMA mostra un errore nella console e interrompe il resto del codice per la classificazione dall'esecuzione.

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

Si verifica anche quando si utilizza un test diverso limma di per la selezione variabile, o utilizzare due invece di 15 variabili per la generazione classificatore. Il codice R che uso dovrebbe garantire che i set di formazione abbiano sempre membri di entrambe le classi. Gradirei qualsiasi idea su questo.

Di seguito è riportato il codice R che uso, con Mac OS X 10.6.6, R 2.12.1, Biobase 2.10.0, CMA 1.8.1, limma 3.6.9 e WilcoxCV 1.0.2. Il file di dati hy3ExpHsaMir.txt può essere scaricato da http://rapidshare.com/files/447062901/hy3ExpHsaMir.txt.

Tutto va bene fino a quando g è 9 nel per (g in 00:10) loop (per variare le dimensioni/prova-set di formazione).

# 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='')

Uscita per 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, ...

Answer 1

Il manutentore del pacchetto CMA prontamente risposto a un messaggio che avevo inviato su questo problema.

La CMA sintonizza un classificatore generato da un set di allenamento testando diversi valori di parametro in un passo di formazione in k set, k-fold (predefinito k = 3). Con piccole dimensioni del set di allenamento, questo ciclo interno può fallire se le osservazioni di una sola classe vengono sub-settate. Due modi per ridurre la possibilità che questo accada sono di fare un CV interno doppio e di specificare il campionamento stratificato, entrambi richiedono che il passo di sintonia sia invocato separatamente tramite tune() di CMA e utilizzi il suo output nella classificazione ().

Nel codice che ho postato, l'ottimizzazione viene invocata dalla classificazione (), che non consente il campionamento stratificato o il CV a 2 pieghe. Tuttavia, invocare il tune() separatamente, per il campionamento stratificato e il CV 2 volte, non ha aiutato nel mio caso. Questo non è sorprendente in quanto, con piccoli insiemi di formazione, il CMA incontra i casi con gruppi di membri di una sola classe.

Desidero che invece di interrompere bruscamente tutto, CMA continuerà a lavorare sui restanti set di apprendimento quando incontra un problema di questo tipo con un set di apprendimento.Sarebbe anche bello se, avendo questo problema, il CMA provasse diversi valori di k per il passo CV k-fold interno.

[Modificato il 14 febbraio] La generazione di set di apprendimento per CV di CMA non verifica se esistono membri sufficienti di entrambe le classi in un set di allenamento. La seguente sostituzione di una parte del codice nel post originale dovrebbe quindi far funzionare correttamente le cose:

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)