Resolviendo un error 'modelo vacío' en la validación cruzada para la clasificación SVM cuando utilizo el paquete CMA Bioconductor para R

Question

Estoy utilizando el paquete Bioconductor CMA para realizar validación cruzada interna de Monte Carlo (MCCV) en clasificadores SVM en un microarray conjunto de datos. CMA usa internamente el paquete e1071 R para el trabajo SVM.

El conjunto de datos tiene 387 variables (atributos) para 45 muestras (observaciones) que pertenecen a una de dos clases (etiquetas 0 o 1, a una proporción de 1: 1). Todos los datos son numéricos sin NA. Estoy probando un MCCV de 1000 iteraciones con 15 variables seleccionadas para SVM usando el limma statistics para el análisis de expresión genética diferencial. Durante MCCV, una fracción del conjunto de 45 muestras se usa para entrenar un clasificador SVM, que luego se usa para probar la fracción restante, y estoy probando diferentes valores para la fracción del conjunto de entrenamiento. CMA también realiza validaciones de bucle interno (validación cruzada de 3 veces dentro de los conjuntos de entrenamiento, por defecto) para ajustar con precisión los clasificadores que se utilizarán para la validación cruzada contra los conjuntos de prueba. Todo esto se hace desde dentro del paquete CMA.

A veces, para los tamaños de conjunto de entrenamiento bajo, CMA muestra un error en la consola y detiene el resto del código para que se ejecute la clasificación.

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

Se produce incluso cuando se utiliza una prueba de que no sea de limma para la selección de variables, o el uso de dos en lugar de 15 las variables para la generación de clasificador. El código R que uso debe garantizar que los conjuntos de entrenamiento siempre tengan miembros de ambas clases. Agradecería cualquier idea sobre esto.

A continuación se muestra el código R que utilizo, con Mac OS X 10.6.6, R 2.12.1, Biobase 2.10.0, CMA 1.8.1, limma 3.6.9 y WilcoxCV 1.0.2. El archivo de datos hy3ExpHsaMir.txt se puede descargar desde http://rapidshare.com/files/447062901/hy3ExpHsaMir.txt.

todo va bien hasta g es 9 en el para (g en doce y diez) bucle (para variar los tamaños de formación/prueba-set).

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

de salida para el rastreo():

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

El responsable del paquete CMA respondió rápidamente a un mensaje que había enviado sobre este tema.

CMA sintoniza un clasificador generado a partir de un conjunto de entrenamiento mediante la prueba de diferentes valores de parámetros en un paso dentro del entrenamiento, k-fold CV (valor predeterminado k = 3). Con pequeños tamaños de conjuntos de entrenamiento, este ciclo interno puede fallar si las observaciones de una sola clase se sub-configuran. Dos formas de reducir la posibilidad de que esto ocurra es hacer un doble CV interior, y especificar el muestreo estratificado, los cuales requieren que el paso de sintonización se invoque por separado a través de la canción de CMA() y usar su salida en la clasificación ().

En el código que publiqué, la sintonización se invoca desde clasificación(), que no permite el muestreo estratificado ni CV doble. Sin embargo, invocar tune() por separado, para muestreo estratificado y CV doble, no ayudó en mi caso. Esto no es sorprendente ya que, con pequeños grupos de entrenamiento, CMA encuentra casos con grupos de miembros de una sola clase.

Deseo que, en lugar de terminar abruptamente todo, CMA continúe trabajando en los conjuntos de aprendizaje restantes cuando se encuentra con un conjunto de aprendizaje.También sería bueno si, al tener este problema, CMA probaría diferentes valores de k para el paso interno de k-fold CV.

[Editado el 14 de febrero] La generación de juegos de aprendizaje de CMA para CV no verifica si existen suficientes miembros de ambas clases en un conjunto de entrenamiento. El siguiente reemplazo para una parte del código en la publicación original debería por lo tanto hacer que las cosas funcionen correctamente:

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)