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

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

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

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 в этом примере). Методы обертки используют производительность выбранного алгоритма обучения, чтобы оценить каждое подмножество функции кандидата. Методы обертки ищут функции, лучше подходящий для выбранного алгоритма обучения, но они могут быть значительно медленнее, чем методы фильтра, если алгоритм обучения занимает много времени, чтобы запуститься. Концепции "фильтров" и "оберток" описаны в Джоне Г. Кохэви Р. (1997) "Обертки для выбора подмножества функции", Искусственный интеллект, 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, означая, что существует больше чем 2 500 функций среди исходных 5000 функций, которые имеют сильную силу дискриминации. Можно отсортировать эти функции согласно их p-значениям (или абсолютные значения t-статистической-величины) и выбрать некоторые функции из отсортированного списка. Однако это обычно затрудняет, чтобы решить, сколько функций необходимо, если у каждого нет некоторых знаний проблемной области или максимального количества функций, которые могут быть рассмотрены, был продиктован заранее на основе внешних ограничений.

Один быстрый способ решить количество необходимых функций состоит в том, чтобы построить MCE (misclassification ошибка, т.е. количество неправильно классифицированных наблюдений, разделенных на количество наблюдений) на наборе тестов как функция количества функций. С тех пор существует только 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')
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)
   NumObservations: 216
       NumTestSets: 1
         TrainSize: 216
          TestSize: 216

Для удобства, 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 =

  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) и второй выбранной функцией (столбец 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 =

        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 =

  Columns 1 through 6

        2814        2721        2720        2452        2650        2731

  Columns 7 through 10

        2337        2658         864        3288

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

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