高次元のデータを分類する特徴量の選択
この例では、高次元データを分類するための特徴量を選択する方法を示します。具体的には、最も一般的な特徴選択アルゴリズムのひとつである逐次特徴選択を実行する方法を示します。ホールドアウトと交差検証を使用して、選択した特徴量の性能を評価する方法についても示します。
統計を学ぶうえで重要なのは、特徴量の数 (次元) を減らすことです。特徴量の数が多く観測の数が限定されたデータ セット (バイオインフォマティクス データなど) が多く存在する場合、特徴量の数が多いことは、目標とする学習結果を得るには適切ではありません。また、観測数が限定されていると、学習アルゴリズムがノイズに過適合する可能性があります。特徴量数を減らすと、ストレージと費やす時間が削減され、より正確に理解できるようになります。
特徴量数を減らすには、特徴選択と特徴変換という 2 つの主な方法があります。特徴選択アルゴリズムでは、オリジナルの特徴セットから特徴量のサブセットを選択します。特徴変換方法では、オリジナルの高次元特徴空間から次元の低い新しい空間にデータを変換します。
データの読み込み
血清プロテオミクス パターン診断を使用すると、疾患のある患者と疾患のない患者の観測を区別することができます。プロファイル パターンは、SELDI (表面増強レーザ脱着/イオン化) タンパク質質量分析を使用して生成されます。これらの特徴量は、特定の質量/充填量値でのイオン強度レベルです。
この例では、WCX2 タンパク質配列を使用して生成された高解像度の卵巣癌データ セットを使用します。Bioinformatics Toolbox™ の Preprocessing Raw Mass Spectrometry Data (Bioinformatics Toolbox) の例で説明されている手順と同様の前処理手順を実行すると、2 つの変数 obs
および grp
がデータ セットに含まれます。変数 obs
は、4000 個の特徴量と 216 個の観測値から構成されています。grp
の各要素によって、obs
の対応する行が属するグループが定義されます。
load ovariancancer;
whos
Name Size Bytes Class Attributes grp 216x1 25056 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 NumObservations: 216 NumTestSets: 1 TrainSize: 160 TestSize: 56 IsCustom: 0
dataTrain = obs(holdoutCVP.training,:); grpTrain = grp(holdoutCVP.training);
すべての特徴量を使用してデータを分類する際の問題
最初に特徴量数を減らさないと、特徴量数が観測数を大幅に上回るため、この例で使用されるデータ セットの分類アルゴリズムは失敗する場合があります。この例では、分類アルゴリズムとして QDA (2 次判別分析) を使用します。次のコマンドに示すように、すべての特徴量を使用してデータに QDA を適用すると、各グループで共分散行列を推定するのに十分な標本を得られないため、誤差が発生します。
try yhat = classify(obs(test(holdoutCVP),:), 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,'Vartype','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 のプロットを示します。選択したモデルの性能を合理的に推定するために重要なのは、160 の学習標本を使用して QDA モデルを当てはめて、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')
resubCVP = Resubstitution (no partition of data) NumObservations: 216 NumTestSets: 1 TrainSize: 216 TestSize: 216 IsCustom: 0
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');
便宜上、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 = 1×15
2814 2813 2721 2720 2452 2645 2644 2642 2650 2643 2731 2638 2730 2637 2398
逐次特徴選択の適用
上記の特徴選択アルゴリズムでは、特徴量間の交互作用が考慮されていません。また、各ランクに基づいてリストから選択された特徴量には、冗長な情報が含まれる可能性があるため、一部の特徴量は必要ありません。たとえば、最初に選択された特徴 (2814 列) と 2 番目に選択された特徴 (2813 列) の間の線形相関係数はほぼ 0.95 となります。
corr(dataTrain(:,featureIdxSortbyP(1)),dataTrain(:,featureIdxSortbyP(2)))
ans = single
0.9447
このような簡単な特徴選択手順は通常、高速であるため、前処理手順として使用されます。より高度な特徴選択アルゴリズムでは、性能が向上します。最も広く使用されている方法の 1 つとして、逐次特徴選択があります。この方法では、特定の停止条件に達するまで連続的に追加 (前方検索) または削除 (後方検索) を行うことで、特徴量のサブセットを選択します。
この例では、ラッパー形式で前方の逐次特徴選択を使用して重要な特徴量を検出します。より具体的には、分類における一般的な目的は、MCE を最小化することであるため、特徴選択手順では、候補特徴サブセットの性能を示すインジケーターとして、各候補特徴サブセットの学習アルゴリズム QDA の MCE を使用して逐次検索を行います。学習セットは特徴量の選択と QDA モデルの当てはめに使用し、テスト セットは最終的に選択された特徴の性能を評価するのに使用します。特徴選択手順で各候補特徴サブセットの性能を評価および比較するには、階層化された 10 分割交差検証を学習セットに適用します。交差検証を学習セットに適用することが重要である理由については、後述します。
まず、学習セットの階層化された 10 分割を生成します。
tenfoldCVP = cvpartition(grpTrain,'kfold',10)
tenfoldCVP = K-fold cross validation partition NumObservations: 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 IsCustom: 0
前のセクションで前処理手順として取得したフィルター結果を使用して特徴量を選択します。たとえば、次のように 150 の特徴量を選択します。
fs1 = featureIdxSortbyP(1:150);
これら 150 の特徴量に前方逐次特徴選択を適用します。関数 sequentialfs
は、必要な特徴量数を決定するための簡単な方法 (既定のオプション) を提供します。交差検証 MCE の最初の局所的最小値が検出されると、関数 sequentialfs は停止します。
fsLocal = sequentialfs(classf,dataTrain(:,fs1),grpTrain,'cv',tenfoldCVP);
選択された特徴量は次のようになります。
fs1(fsLocal)
ans = 1×3
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 = 1×10
2814 2721 2720 2452 2650 2731 2337 2658 864 3288
前方逐次特徴選択の手順で選択された順序でこれら 10 の特徴量を表示するには、historyCV
出力でこれらの特徴量がはじめて true になる行を検索します。
[orderlist,ignore] = find( [historyCV.In(1,:); diff(historyCV.In(1:10,:) )]' ); fs1(orderlist)
ans = 1×10
2337 864 3288 2721 2814 2658 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 倍となっています。ここでも、特徴量を評価および選択する際に性能の推定として再代入誤差を使用するのは適切ではありません。最終的な評価手順の場合だけでなく、特徴選択手順の場合にも、再代入誤差の使用を避けることをお勧めします。