5

Bioconductor パッケージCMAを使用して、マイクロ アレイ データセットの SVM 分類器で内部のモンテカルロ交差検証 (MCCV) を実行しています。CMA は、SVM 作業に e1071 R パッケージを内部的に使用します。

データセットには、2 つのクラス (ラベル 0 または 1、約 1:1 の比率) のいずれかに属する 45 のサンプル (観察) に対する 387 の変数 (属性) があります。すべてのデータは数値であり、NA はありません。差次的遺伝子発現解析のリマ統計を使用して、SVM 用に選択された 15 の変数で 1000 反復 MCCV を試みています。MCCV では、45 サンプル セットの一部が SVM 分類器のトレーニングに使用され、残りの部分をテストするために使用されます。トレーニング セットの部分にさまざまな値を試しています。CMA は、内部ループ検証 (デフォルトでは、トレーニング セット内の 3 分割交差検証) も実行して、テスト セットに対する交差検証に使用される分類子を微調整します。これはすべて、CMA パッケージ内から実行されます。

トレーニング セットのサイズが小さい場合、CMA はコンソールにエラーを表示し、残りの分類コードの実行を停止します。

[中略]チューニング反復 575
チューニング反復 576
チューニング反復 577
predict.svm(ret, xhold, decision.values = TRUE) のエラー: モデルが空です!

変数の選択に limma 以外のテストを使用したり、分類器の生成に 15 個ではなく 2 個の変数を使用したりした場合でも発生します。私が使用する 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 で使用する R コードです。データ ファイル hy3ExpHsaMir.txt は、 http: //rapidshare.com/files/447062901/hy3ExpHsaMir.txt からダウンロードできます。

for(g in 0:10)ループでgが 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("モデルが空です!")
19: predict.svm(ret, xhold, decision.values = TRUE)
18: 予測 (ret、xhold、decision.values = TRUE)
17: na.action(predict(ret, xhold, decision.values = TRUE))
16: svm.default(コスト = 0.1、カーネル = "線形"、タイプ = "C 分類"、...
15: svm(cost = 0.1, kernel = "linear", type = "C-classification", ...
14: do.call("svm", args = ll)
13: 関数 (X、y、f、学習、確率、モデル = FALSE、...) ...
12: 関数 (X、y、f、学習、確率、モデル = FALSE、...) ...
11: do.call(classifier, args = c(list(X = X, y = y, learnind = 学習行列[i, ...
10: 分類(X = c(83.5832768669369, 83.146333099001, 94.253534443549, ...
9: 分類(X = c(83.5832768669369, 83.146333099001, 94.253534443549, ...
8: do.call("分類", args = c(list(X = Xi, y = yi, learningsets = lsi, ...
7: 調整(グリッド = リスト(コスト = c(0.1, 0.2, 0.5, 1, 2, 5, 10, 20, 50...
6: 調整(グリッド = リスト(コスト = c(0.1, 0.2, 0.5, 1, 2, 5, 10, 20, 50...
5: do.call("tune", args = c(tuninglist, ll))
4: 分類(X, y = as.numeric(y) - 1, learningsets = learningsets, ...
3: 分類(X, y = as.numeric(y) - 1, learningsets = learningsets, ...
2: classification(t(exp), factor(classes), learningsets = lsets, ...
1: 分類(t(exp), 因子(クラス), learningsets = lsets, ...
4

1 に答える 1

3

CMA パッケージのメンテナーは、私がこの問題について送ったメッセージにすぐに返信しました。

CMA は、トレーニング セット内の k 倍 CV ステップ (デフォルト k=3) でさまざまなパラメーター値をテストすることにより、トレーニング セットから生成された分類子を調整します。トレーニング セットのサイズが小さいと、1 つのクラスのみの観測がサブセット化されると、この内側のループが失敗する可能性があります。これが発生する可能性を減らす 2 つの方法は、2 倍の内部 CV を実行することと、層別サンプリングを指定することです。どちらの方法でも、CMA のtune()を介してチューニング ステップを個別に呼び出し、その出力をclassification()で使用する必要があります。

私が投稿したコードでは、チューニングはclassification()内から呼び出されます。これは、層別サンプリングまたは 2 倍の CV を許可しません。ただし、階層化サンプリングと 2 倍 CV のためにtune()を個別に呼び出しても、私の場合は役に立ちませんでした。小さなトレーニング セットでは、CMA は 1 つのクラスのみのメンバーのセットを持つケースに遭遇するため、これは驚くべきことではありません。

1 つの学習セットでこのような問題が発生した場合、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)
于 2011-02-10T20:51:08.547 に答える