ドキュメンテーション センター

  • 評価版
  • 製品アップデート

最新のリリースでは、このページがまだ翻訳されていません。 このページの最新版は英語でご覧になれます。

高次元のデータを分類する特徴の選択

この例では、高次元データを分類するための機能を選択する方法を示します。具体的には、最も一般的な特徴選択アルゴリズムのひとつである逐次特徴選択を実行する方法を示します。ホールドアウトと交差検定を使用して、選択した特徴の性能を評価する方法についても示します。

統計を学ぶうえで重要なのは、特徴の数 (次元) を減らすことです。特徴の数が多く観測の数が限定されたデータセット (バイオインフォマティクス データなど) が多く存在する場合、特徴の数が多いことは、目標とする学習結果を得るには適切ではありません。また、観測数が限定されていると、学習アルゴリズムがノイズに過剰適合する可能性があります。 特徴数を減らすと、ストレージと費やす時間が削減され、より正確に理解できるようになります。

特徴数を減らすには、特徴選択と特徴変換という 2 つの主な方法があります。特徴選択アルゴリズムでは、オリジナルの特徴セットから特徴のサブセットを選択します。特徴変換方法では、オリジナルの高次元特徴空間から次元の低い新しい空間にデータを変換します。

データの読み込み

血清プロテオミクス パターン診断を使用すると、疾患のある患者と疾患のない患者の観測を区別することができます。プロファイル パターンは、SELDI (表面増強レーザ脱着/イオン化) タンパク質質量分析を使用して生成されます。これらの特徴は、特定の質量/充填量値でのイオン強度レベルです。

この例のデータについては、「FDA-NCI 臨床プロテオミクス プログラム データ バンク」を参照してください。この例では、WCX2 タンパク質配列を使用して生成された高解像度の卵巣癌データセットを使用します。Bioinformatics Toolbox™ の Pre-processing Raw Mass Spectrometry Data の例で説明されている手順と同様の前処理手順を実行した後は、データセットには、2 つの変数 (obsgrp) が含まれます。変数 obs は、5000 の特徴と 216 の観測から構成されます。grp の各要素によって、obs の対応する行が属するグループが定義されます。

load ovariancancer;
whos
  Name        Size                Bytes  Class     Attributes

  grp       216x1                 26784  cell                
  obs       216x4000            3456000  single              

トレーニング セットとテスト セットへのデータの分割

この例で使用される関数の一部は、MATLAB の組み込みの乱数発生関数を呼び出します。この例で示された結果を再現するには、次のコマンドを実行して乱数発生器を既知の状態に設定します。コマンドを実行しないと、結果が異なる場合があります。

rng(8000,'twister');

訓練データの性能 (再置換の性能) は、独立した検定セットにおけるモデルの性能を推測するには適切ではありません。再置換性能は楽観的になりすぎる場合が多くあります。選択したモデルの性能を予測するには、モデルの構築に使用されなかった別のデータセットでの性能を評価する必要があります。この例題では、cvpartition を使用して、サイズ 160 のトレーニング セットとサイズ 56 のテスト セットにデータを分割します。トレーニング セットとテスト セットの比率は、grp のグループ比率とおおよそ同じです。訓練データを使用して特徴を選択し、学習データの選択された特徴の性能を検証します。このような方法は通常、ホールドアウト検定と呼ばれます。モデルを評価および選択するために簡単で広く使用されているもう 1 つの方法として交差検定があります。交差検定についてはこの例の後半で説明します。

holdoutCVP = cvpartition(grp,'holdout',56)
holdoutCVP = 

Hold-out cross validation partition
             N: 216
   NumTestSets: 1
     TrainSize: 160
      TestSize: 56
dataTrain = obs(holdoutCVP.training,:);
grpTrain = grp(holdoutCVP.training);

すべての特徴を使用してデータを分類する際の問題

最初に特徴数を減らさないと、特徴数が観測数を大幅に上回るため、この例で使用されるデータセットの分類アルゴリズムは失敗する場合があります。この例では、分類アルゴリズムとして QDA (2 次判別分析) を使用します。次のコマンドに示すように、すべての特徴を使用してデータに QDA を適用すると、各グループで共分散行列を推定するのに十分な標本を得られないため、誤差が発生します。

try
   yhat = classify(obs(holdoutCVP.test(),:), dataTrain, grpTrain,'quadratic');
catch ME
   display(ME.message);
end
The covariance matrix of each group in TRAINING must be positive definite.

簡単なフィルター方法を使用した特徴の選択

ここでの目的は、優れた分類性能を得ることができる重要な特徴の小規模なセットを見つけ出して、データの次元を減らすことです。特徴選択アルゴリズムは大きく分けて 2 種類あります。フィルター方法とラッパー方法です。フィルター方法は、データの一般的特性に依存して、選択された学習アルゴリズム (この例では QDA) を含まない特徴のサブセットを評価および選択します。ラッパー方法では、選択された学習アルゴリズムの性能を使用して、各候補特徴のサブセットを評価します。ラッパー方法では、選択された学習アルゴリズムに最適な特徴を検索しますが、学習アルゴリズムの実行時間が長い場合、フィルター方法より速度が大幅に低下する可能性があります。「フィルター」と「ラッパー」の概念は、文献 (John G. Kohavi R. (1997) "Wrappers for feature subset selection", Artificial Intelligence, Vol.97, No.1-2, pp.272-324) を参照してください。この例では、フィルター方法とラッパー方法の例を 1 つずつ示します。

フィルターは通常、簡単かつ高速であるため、前処理手順として使用されます。バイオインフォマティクス データで広く使用されるフィルター方法では、特徴間に交互作用がないことを前提に、各特徴に一変量基準を個別に適用します。

たとえば、各特徴に t 検定を適用して、グループの分割における有効性の測定として各特徴の p 値(または t 統計量の絶対値) を比較するとします。

dataTrainG1 = dataTrain(grp2idx(grpTrain)==1,:);
dataTrainG2 = dataTrain(grp2idx(grpTrain)==2,:);
[h,p,ci,stat] = ttest2(dataTrainG1,dataTrainG2,[],[],'unequal');

各特徴によって 2 つのグループがどの程度適切に分割されているかについて一般的な概念を得るには、p 値の実測の累積分布関数 (CDF) をプロットします。

ecdf(p);
xlabel('P value');
ylabel('CDF value')

0 値に近い p 値をもつ特徴は約 35%、0.05 より小さな p 値をもつ特徴は 50%以上 存在します。つまり、識別力の大きいオリジナルの 5000 の特徴のうち 2500 以上の特徴ということになります。特徴の p 値 (または t 統計量の絶対値) に応じてこれらの特徴を並べ替えて、並べ替えられたリストから特徴を選択することが可能ですが、ドメインに熟知している、または検討可能な特徴の最大数が外部の制約に基づいてあらかじめ指定されている場合を除き、通常、必要な特徴数を決定することは困難です。

必要な特徴数を決定する簡単な方法の 1 つには、テスト セットの MCE (誤判別エラー、観測数によって除算した誤判別観測数など) を特徴数の関数としてプロットします。トレーニング セットには 160 の観測しか存在しないため、QDA に適用できる特徴の最大数は限られます。それ以外の場合、共分散行列を推定するのに十分な標本が各グループに存在しないことがあります。実際、この例で使用されるデータの場合、2 つのグループのホールドアウト分割とサイズによって、QDA を適用できる特徴の最大許容数は約 70 であると指定されます。ここで、5 から 70 までの特徴数の MCE を計算し、特徴数の関数として MCE のプロットを示します。選択したモデルの性能を合理的に推定するために重要なのは、QDA モデルに適合する 160 の訓練標本を使用して、56 のテスト観測の MCE (次のプロットで青の円マーク) を計算することです。また、誤判別率がテスト誤差の優れた推定誤差ではない理由を説明するため、誤判別 MCE を赤の三角形マークで示しています。

[~,featureIdxSortbyP]= sort(p,2); % sort the features
testMCE =zeros(1,14);
resubMCE = zeros(1,14);
nfs = 5:5:70;
classf = @(xtrain,ytrain,xtest,ytest) ...
             sum(~strcmp(ytest,classify(xtest,xtrain,ytrain,'quadratic')));
resubCVP = cvpartition(length(grp),'resubstitution')
for i=1:14
   fs = featureIdxSortbyP(1:nfs(i));
   testMCE(i) = crossval(classf,obs(:,fs),grp,'partition',holdoutCVP)...
       /holdoutCVP.TestSize;
   resubMCE(i) = crossval(classf,obs(:,fs),grp,'partition',resubCVP)/...
       resubCVP.TestSize;
end
 plot(nfs, testMCE,'o',nfs,resubMCE,'r^');
 xlabel('Number of Features');
 ylabel('MCE');
 legend({'MCE on the test set' 'Resubstitution MCE'},'location','NW');
 title('Simple Filter Feature Selection Method');
resubCVP = 

Resubstitution (no partition of data)
             N: 216
   NumTestSets: 1
     TrainSize: 216
      TestSize: 216

便宜上、classf は無名関数として定義されています。また、classf は、指定されたトレーニング セットの QDA と適合し、指定されたテスト セットの誤判別の標本数を返します。独自の分類アルゴリズムを作成する場合、次に示すように、classf を別のファイルに配置することをお勧めします。

%  function err = classf(xtrain,ytrain,xtest,ytest)
%       yfit = classify(xtest,xtrain,ytrain,'quadratic');
%        err = sum(~strcmp(ytest,yfit));

再置換 MCE はかなり楽観的になっています。MCE では、より多くの特徴が使用されると常に減少し、60 を超える特徴が使用されると 0 になります。ただし、誤判別率が減少しているときにテスト誤差が増加すると、過剰適合が発生する場合があります。この簡単なフィルター特徴選択方法では、15 の特徴が使用されると、テスト セットで最小の MCE が取得されます。プロットで 20 以上の特徴が使用されると、過剰適合を示します。テスト セットの最小の MCE は 12.5% です。

testMCE(3)
ans =

    0.1250

これらは、最小 MCE を達成する最初の 15 の特徴です。

featureIdxSortbyP(1:15)
ans =

  Columns 1 through 6

        2814        2813        2721        2720        2452        2645

  Columns 7 through 12

        2644        2642        2650        2643        2731        2638

  Columns 13 through 15

        2730        2637        2398

逐次特徴選択の適用

上記の特徴選択アルゴリズムでは、特徴間の相互作用が考慮されていません。また、各ランクに基づいてリストから選択された特徴には、冗長な情報が含まれる可能性があるため、一部の特徴は必要ありません。たとえば、最初に選択された特徴 (2814 列) と 2 番目に選択された特徴 (2813 列) の間の線形相関係数はほぼ 0.95 となります。

corr(dataTrain(:,featureIdxSortbyP(1)),dataTrain(:,featureIdxSortbyP(2)))
ans =

    0.9447

このような簡単な特徴選択手順は通常、高速であるため、前処理手順として使用されます。より高度な特徴選択アルゴリズムでは、性能が向上します。最も広く使用されている方法の 1 つとして、逐次特徴選択があります。この方法では、特定の停止条件に達するまで連続的に追加 (前方検索) または削除 (後方検索) を行うことで、特徴のサブセットを選択します。

この例では、ラッパー形式で前方の逐次特徴選択を使用して重要な特徴を検出します。より具体的には、分類における一般的な目的は、MCE を最小化することであるため、特徴選択手順では、候補特徴サブセットの性能を示すインジケーターとして、各候補特徴サブセットの学習アルゴリズム QDA の MCE を使用して逐次検索を行います。特徴を選択して QDA モデルに適合するにはトレーニング セットを使用し、最終的に選択された特徴の性能を評価するにはテスト セットを使用します。特徴選択手順で各候補特徴サブセットの性能を評価および比較するには、階層化された 10 分割交差検定をトレーニング セットに適用します。交差検定をトレーニング セットに適用することが重要である理由については、後述します。

まず、トレーニング セットの階層化された 10 分割を生成します。

tenfoldCVP = cvpartition(grpTrain,'kfold',10)
tenfoldCVP = 

K-fold cross validation partition
             N: 160
   NumTestSets: 10
     TrainSize: 144  144  144  144  144  144  144  144  144  144
      TestSize: 16  16  16  16  16  16  16  16  16  16

前のセクションで前処理手順として取得したフィルター結果を使用して特徴を選択します。たとえば、次のように 150 の特徴を選択します。

fs1 = featureIdxSortbyP(1:150);

これら 150 の特徴に前方逐次特徴選択を適用します。関数 sequentialfs は、必要な特徴数を決定するための簡単な方法 (既定のオプション) を提供します。交差検定 MCE の最初の局所的最小値が検出されると、関数 sequentialfs は停止します。

 fsLocal = sequentialfs(classf,dataTrain(:,fs1),grpTrain,'cv',tenfoldCVP);

選択された特徴は次のようになります。

fs1(fsLocal)
ans =

        2337         864        3288

これらの 3 つの特徴をもつ選択されたモデルの性能を評価するには、56 のテスト標本について MCE を計算します。

testMCELocal = crossval(classf,obs(:,fs1(fsLocal)),grp,'partition',...
    holdoutCVP)/holdoutCVP.TestSize
testMCELocal =

    0.0714

3 つの特徴だけが選択されているため、MCE は、簡単なフィルター特徴選択方法を使用して最小の MCE の半分強となります。

アルゴリズムは途中で停止する場合があります。合理的な範囲の特徴数ではなく交差検定 MCE の最小値を検索することで、より小さな MCE を達成できる場合があります。たとえば、最大 50 の特徴の特徴数の関数として交差検定 MCE のプロットを次のように描画します。

[fsCVfor50,historyCV] = sequentialfs(classf,dataTrain(:,fs1),grpTrain,...
    'cv',tenfoldCVP,'Nf',50);
plot(historyCV.Crit,'o');
xlabel('Number of Features');
ylabel('CV MCE');
title('Forward Sequential Feature Selection with cross-validation');

10 の特徴が使用されると、交差検定 MCE は最小値に達し、10 から 35 までの特徴の範囲でこの曲線は一定に保たれます。35 を超える特徴が使用されると、曲線は右上がりになり、過剰適合が発生したことを意味します。

通常、より少ない特徴数が適切であるため、ここでは 10 の特徴を選択します。

fsCVfor10 = fs1(historyCV.In(10,:))
fsCVfor10 =

  Columns 1 through 6

        2814        2721        2720        2452        2650        2731

  Columns 7 through 10

        2337        2658         864        3288

前方逐次特徴選択の手順で選択された順序でこれら 10 の特徴を表示するには、historyCV 出力でこれらの特徴がはじめて true になる行を検索します。

[orderlist,ignore] = find( [historyCV.In(1,:); diff(historyCV.In(1:10,:) )]' );
fs1(orderlist)
ans =

  Columns 1 through 6

        2337         864        3288        2721        2814        2658

  Columns 7 through 10

        2452        2731        2650        2720

これら 10 の特徴を評価するには、検定セットで QDA の MCE を計算します。これまでで最小の MCE 値が取得されます。

testMCECVfor10 = crossval(classf,obs(:,fsCVfor10),grp,'partition',...
    holdoutCVP)/holdoutCVP.TestSize
testMCECVfor10 =

    0.0357

興味深いのは、特徴数の関数として (たとえば、特徴選択手順の実行時に交差検定を実行せずに) トレーニング セットの誤判別 MCE 値のプロットを調べた場合です。

[fsResubfor50,historyResub] = sequentialfs(classf,dataTrain(:,fs1),...
     grpTrain,'cv','resubstitution','Nf',50);
plot(1:50, historyCV.Crit,'bo',1:50, historyResub.Crit,'r^');
xlabel('Number of Features');
ylabel('MCE');
legend({'10-fold CV MCE' 'Resubstitution MCE'},'location','NE');

ここでも、再置換 MCE の値はかなり楽観的です。誤判別 MCE の値は交差検定 MCE より小さい値が多く、16 の特徴が使用されると、誤判別 MCE は 0 になります。テスト セットでこれら 16 の特徴の MCE 値を計算して、実際の性能を確認できます。

fsResubfor16 = fs1(historyResub.In(16,:));
testMCEResubfor16 = crossval(classf,obs(:,fsResubfor16),grp,'partition',...
    holdoutCVP)/holdoutCVP.TestSize
testMCEResubfor16 =

    0.0714

testMCEResubfor16 は、テスト セットの (特徴選択手順の誤判別によって選択された) 16 の特徴の性能で、テスト セットの (特徴選択手順の 10 分割交差検定によって選択された) 10 の特徴の性能である testMCECVfor10 の約 2 倍となっています。ここでも、特徴を評価および選択する際に性能の推定として誤判別率を使用するのは適切ではありません。最終的な評価手順の場合だけでなく、特徴選択手順の場合にも、誤判別率の使用を避けることをお勧めします。

この情報は役に立ちましたか?