选择用于高维数据分类的特征 -凯发k8网页登录

加载数据

血清蛋白质组图谱诊断可用于区分患者和非患者的观测值。图谱是使用表面增强激光解吸电离 (seldi) 蛋白质质谱法产生的。这些特征是特定质荷比下的离子强度水平。

此示例用到了使用 wcx2 蛋白质阵列生成的高分辨率卵巢癌数据集。经过一些预处理步骤（类似于 bioinformatics toolbox™ 示例 (bioinformatics 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) "wrappers for feature subset selection", 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');

为了大致了解每个特征区分两个组的效果，我们可以绘制 p 值的经验累积分布函数 (cdf)：

ecdf(p);
xlabel('p value'); 
ylabel('cdf value')

大约 35% 的特征的 p 值接近零，超过 50% 的特征的 p 值小于 0.05，这意味着 5000 个原始特征中有超过 2500 个特征具有强大的区分能力。您可以根据 p 值（或 t 统计量的绝对值）对这些特征进行排序，然后从排序的列表中选择一些特征。但是，除非您具备一定的专业知识，或者已经基于外部约束提前确定可以考虑的最大特征数，否则通常很难确定需要多少个特征。

要确定需要多少个特征，一种快速的方法是绘制测试集上的 mce 对特征数的函数，mce 是误分类误差，即误分类的观测值数除以观测值总数。由于训练集中只有 160 个观测值，因此应用 qda 的最大特征数是有限的，否则，每个组中可能没有足够的样本来估计协方差矩阵。实际上，对于本示例中使用的数据，留出法分区和两个组的大小决定应用 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

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 时，它会减少到零。但是，如果测试误差升高而再代入误差继续下降，则可能发生过拟合。当使用 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 列）和第二个特征（2813 列）之间的线性相关系数几乎为 0.95。

corr(datatrain(:,featureidxsortbyp(1)),datatrain(:,featureidxsortbyp(2)))

ans = single
    0.9447

这种简单的特征选择过程速度很快，因而通常用作预处理步骤。更高级的特征选择算法则用于提高性能。序列特征选择是使用最广泛的方法之一。它通过按顺序添加特征（前向搜索）或删除特征（后向搜索）来选择特征子集，直到满足某些停止条件为止。

在此示例中，我们以包装器方法使用前向序列特征选择来查找重要特征。具体而言，由于分类的目的通常是为了使 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

然后，我们将上一节中的过滤作为预处理步骤，使用其结果来选择特征。例如，我们在此选择 150 个特征：

fs1 = featureidxsortbyp(1:150);

对这 150 个特征应用前向序列特征选择。函数 sequentialfs 提供一种简单的方法（默认选项）来决定需要多少个特征。它在找到交叉验证 mce 的第一个局部最小值时停止。

 fslocal = sequentialfs(classf,datatrain(:,fs1),grptrain,'cv',tenfoldcvp);

所选特征如下：

fs1(fslocal)

ans = 1×3
        2337         864        3288

为了评估具有这三个特征的所选模型的性能，我们基于 56 个测试样本计算 mce。

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

当使用 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 值，且使用 16 个特征时，再代入 mce 值开始变为零。我们可以计算这 16 个特征针对测试集的 mce 值，以查看它们的真实性能：

fsresubfor16 = fs1(historyresub.in(16,:));
testmceresubfor16 = crossval(classf,obs(:,fsresubfor16),grp,'partition',...
    holdoutcvp)/holdoutcvp.testsize

testmceresubfor16 = 0.0714

这 16 个特征（在特征选择过程中通过再代入选出）针对测试集的性能 testmceresubfor16 大约是 10 个特征（在特征选择过程中通过 10 折交叉验证选出）针对测试集的性能 testmcecvfor10 的两倍。这再次表明，在评估和选择特征时，再代入误差往往并不能很好地反映性能。不管是在最终评估期间，还是在特征选择过程中，我们都要避免使用再代入误差。