Загрузка данных

Сывороточная протеомная диагностика может быть использована для дифференциации наблюдений от пациентов с заболеванием и без него. Профили формируются с использованием масс-спектрометрии белка с улучшенной лазерной десорбцией и ионизацией (SELDI). Эти признаки представляют собой уровни интенсивности ионов при определенных значениях массы/заряда.

В этом примере используется набор данных о раке яичников с высоким разрешением, который был создан с использованием массива белков WCX2. После некоторых этапов предварительной обработки, аналогичных тем, которые показаны в биоинформатике Toolbox™ пример предварительной обработки необработанных масс-спектрометрических данных, набор данных имеет две переменные obs и grp. obs переменная состоит из 216 наблюдений с 4000 элементами. Каждый элемент в 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. Мы выбираем функции с помощью обучающих данных и оцениваем эффективность выбранных функций на тестовых данных. Это часто называется проверкой удержания. Другим простым и широко используемым методом оценки и выбора модели является перекрестная проверка, которая будет проиллюстрирована ниже в этом примере.

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

holdoutCVP = 
Hold-out cross validation partition
   NumObservations: 216
       NumTestSets: 1
         TrainSize: 160
          TestSize: 56

dataTrain = obs(holdoutCVP.training,:);
grpTrain = grp(holdoutCVP.training);

Проблема классификации данных с использованием всех функций

Без предварительного уменьшения количества признаков некоторые алгоритмы классификации могут привести к сбою в наборе данных, используемом в этом примере, поскольку количество признаков намного больше, чем количество наблюдений. В этом примере в качестве алгоритма классификации используется квадратичный дискриминантный анализ (QDA). Если мы применим 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.

Выбор элементов с использованием простого подхода фильтра

Наша цель - уменьшить размерность данных, найдя небольшой набор важных характеристик, которые могут обеспечить хорошую эффективность классификации. Алгоритмы выбора элементов можно грубо сгруппировать в две категории: методы фильтрации и методы обертки. Методы фильтрации основаны на общих характеристиках данных для оценки и выбора подмножеств признаков без участия выбранного алгоритма обучения (QDA в этом примере). Методы обертки используют производительность выбранного алгоритма обучения для оценки каждого подмножества признаков-кандидатов. Методы обертки ищут функции, лучше подходящие для выбранного алгоритма обучения, но они могут быть значительно медленнее, чем методы фильтрации, если выполнение алгоритма обучения занимает много времени. Понятия «фильтры» и «обертки» описаны в John G. Kohavi R. (1997) «Обертки для выбора подмножества признаков», Artificial Intelligence, Vol.97, No.1-2, pp.272-324. В этом примере показан один экземпляр метода фильтра и один экземпляр метода обертки.

Фильтры обычно используются в качестве этапа предварительной обработки, поскольку они просты и быстры. Широко используемый метод фильтрации для данных биоинформатики заключается в применении одномерного критерия отдельно к каждому элементу, предполагая, что между элементами нет взаимодействия.

Например, мы можем применить t-тест к каждой функции и сравнить p-значение (или абсолютные значения t-статистики) для каждой функции как меру того, насколько она эффективна при разделении групп.

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

Чтобы получить общее представление о том, насколько хорошо разделены две группы по каждой особенности, мы построим график эмпирической кумулятивной функции распределения (CDF) p-значений:

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

Существует около 35% признаков, имеющих p-значения, близкие к нулю, и более 50% признаков, имеющих p-значения, меньшие 0,05, что означает, что среди исходных 5000 признаков имеется более 2500 признаков, которые обладают сильной способностью различения. Можно отсортировать эти функции по их p-значениям (или абсолютным значениям t-статистики) и выбрать некоторые функции из отсортированного списка. Однако, как правило, трудно решить, сколько функций необходимо, если только не известно о какой-либо области или максимальное количество функций, которые можно рассмотреть, было заранее продиктовано на основе внешних ограничений.

Одним из быстрых способов определения количества необходимых элементов является построение графика MCE (ошибка неправильной классификации, т.е. количество неправильно классифицированных наблюдений, деленное на количество наблюдений) на тестовом наборе в зависимости от количества элементов. Поскольку в обучающем наборе всего 160 наблюдений, наибольшее количество признаков для применения QDA ограничено, в противном случае в каждой группе может не хватить выборок для оценки ковариационной матрицы. Фактически, для данных, используемых в этом примере, раздел удержания и размеры двух групп диктуют, что наибольшее допустимое количество признаков для применения QDA составляет около 70. Теперь мы вычисляем MCE для различных чисел признаков между 5 и 70 и показываем график MCE как функцию от числа признаков. Чтобы разумно оценить производительность выбранной модели, важно использовать 160 обучающих образцов для соответствия модели QDA и вычислить MCE на 56 тестовых наблюдениях (синие круглые метки на следующем графике). Чтобы проиллюстрировать, почему ошибка повторного замещения не является хорошей оценкой ошибки теста, мы также показываем 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

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 определяется как анонимная функция. Оно соответствует QDA для данного обучающего набора и возвращает количество неправильно классифицированных образцов для данного тестового набора. При разработке собственного алгоритма классификации его можно поместить в отдельный файл следующим образом:

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

МСЕ по реституции чрезмерно оптимистичен. Оно последовательно уменьшается при использовании большего количества функций и падает до нуля при использовании более 60 функций. Однако если ошибка теста увеличивается, в то время как ошибка повторного замещения все еще уменьшается, то может произойти переоборудование. Этот простой метод выбора функции фильтра получает наименьший MCE в тестовом наборе при использовании 15 функций. График показывает, что подгонка начинает происходить, когда используются 20 или более элементов. Наименьшее значение MCE в тестовом наборе составляет 12,5%:

testMCE(3)

ans = 0.1250

Это первые 15 функций, которые достигают минимального MCE:

featureIdxSortbyP(1:15)

ans = 1×15

        2814        2813        2721        2720        2452        2645        2644        2642        2650        2643        2731        2638        2730        2637        2398

Применение последовательного выбора элементов

Вышеупомянутый алгоритм выбора признаков не учитывает взаимодействие между признаками; кроме того, признаки, выбранные из списка на основе их индивидуального ранжирования, могут также содержать избыточную информацию, так что не все признаки необходимы. Например, коэффициент линейной корреляции между первым выбранным признаком (столбец 2814) и вторым выбранным признаком (столбец 2813) составляет почти 0,95.

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

ans = single
    0.9447

Этот вид простой процедуры выбора признаков обычно используется в качестве этапа предварительной обработки, поскольку он является быстрым. Более совершенные алгоритмы выбора функций повышают производительность. Последовательный выбор признаков является одним из наиболее широко используемых методов. Он выбирает подмножество элементов, последовательно добавляя (прямой поиск) или удаляя (обратный поиск) до тех пор, пока не будут выполнены определенные условия остановки.

В этом примере для поиска важных элементов используется прямой последовательный выбор элементов оберткой. Более конкретно, поскольку типичной целью классификации является минимизация MCE, процедура выбора признака выполняет последовательный поиск с использованием MCE алгоритма QDA обучения по каждому подмножеству признаков-кандидатов в качестве индикатора эффективности для этого подмножества. Обучающий набор используется для выбора элементов и соответствия модели 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

Затем мы используем результаты фильтра из предыдущего раздела в качестве шага предварительной обработки для выбора элементов. Например, здесь мы выбираем 150 функций:

fs1 = featureIdxSortbyP(1:150);

Мы применяем прямой последовательный выбор функций для этих 150 функций. Функция sequentialfs предоставляет простой способ (параметр по умолчанию) определения количества необходимых элементов. Он останавливается при обнаружении первого локального минимума MCE перекрестной проверки.

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

Выбранными элементами являются:

fs1(fsLocal)

ans = 1×3

        2337         864        3288

Чтобы оценить производительность выбранной модели с этими тремя функциями, мы вычисляем MCE на 56 тестовых образцах.

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

testMCELocal = 0.0714

При выборе только трех элементов MCE составляет лишь немногим более половины наименьшего MCE с использованием простого метода выбора элементов фильтра.

Возможно, алгоритм преждевременно остановлен. Иногда меньший MCE достижим путем поиска минимума MCE перекрестной проверки в разумном диапазоне характеристик. Например, мы рисуем график MCE перекрестной проверки как функцию количества элементов для 50 элементов.

[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');

MCE перекрестной проверки достигает минимального значения при использовании 10 элементов, и эта кривая остается плоской в диапазоне от 10 элементов до 35 элементов. Кроме того, кривая поднимается, когда используется более 35 элементов, что означает, что там происходит переоборудование.

Обычно предпочтительнее иметь меньше элементов, поэтому здесь мы выбираем 10 элементов:

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

fsCVfor10 = 1×10

        2814        2721        2720        2452        2650        2731        2337        2658         864        3288

Чтобы показать эти 10 функций в том порядке, в котором они выбраны в процедуре последовательного переноса, мы найдем строку, в которой они впервые становятся истинными в historyCV выходные данные:

[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 функций, мы вычисляем их MCE для QDA на тестовом наборе. Мы получаем наименьшее значение 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 этих 16 функций на тестовом наборе, чтобы увидеть их реальную производительность:

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

testMCEResubfor16 = 0.0714

testMCEResubfor16, производительность этих 16 функций (выбранных путем повторного замещения во время процедуры выбора функции) на тестовом наборе примерно вдвое выше, чем для testMCECVfor10, производительность 10 функций (выбранных 10-кратной перекрестной проверкой во время процедуры выбора функций) на тестовом наборе. Это снова указывает на то, что ошибка повторного замещения, как правило, не является хорошей оценкой производительности для оценки и выбора признаков. Возможно, потребуется избежать использования ошибки повторного предоставления не только на заключительном этапе оценки, но и во время процедуры выбора функций.