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

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

Этот пример использует набор данных о раке яичников высокого разрешения, который был сгенерирован с использованием массива WCX2 белка. После некоторых этапов предварительной обработки, аналогичных тем, которые показаны в Bioinformatics Toolbox™ пример Предварительная обработка необработанных данных масс-спектрометрии, набор данных имеет две переменные obs и grp. The 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 в этом примере). Методы обертки используют эффективность выбранного алгоритма обучения, чтобы оценить каждый подмножество функций кандидата. Методы Wrapper ищут функции, лучше подходящие для выбранного алгоритма обучения, но они могут быть значительно медленнее, чем методы фильтрации, если алгоритм обучения запускается долго. Концепции «фильтры» и «обертки» описаны в 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 ограничено, в противном случае в каждой группе может не хватить выборок для оценки ковариационной матрицы. На самом деле, для данных, используемых в этом примере, секция holdout и размеры двух групп диктуют, что самое большое допустимое количество функций для применения 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));

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