第 18.3 节中的方法使用了 $L_2$ 的惩罚项来对参数正则化,这与在岭回归中相同。所有的估计系数都非零,也因此没有进行特征选择。本节将介绍使用 $L_1$ 惩罚项的方法,也因此提供了自动的特征选择。
回溯第 3.4.2 节中套索回归(lasso):
$$\min_\beta \frac{1}{2} \sum_{i=1}^N (y_i - \beta_0 - \sum_{j=1}^p x_{ij} \beta_j)^2 + \lambda \sum_{j=1}^p |\beta_j| \tag{18.18}$$这里写成了它的拉格朗日形式(同式 3.52)。如同之前的讨论,对一个足够大的调节参数 $\lambda$ 取值,$L_1$ 的惩罚项导致估计系数 $\hat{\beta}_j$ 中的一个子集将等于零。
第 3.4.4 节介绍了最小角回归(LARS)算法,一个计算所有 $\lambda$ 的套索解的高效方法。当 $p>N$ 时(与本章一样),随着 $\lambda$ 趋近于零,套索回归会完全地拟合训练集数据。实际上,根据凸对偶性(convex duality)可以证明当 $p>N$ 时,在所有的 $\lambda$ 取值下非零系数的个数最大为 $N$(例如 Rosset and Zhu, 2007)。因此,套索方法提供了特征选择的一个(严格的)方式。
通过将输出变量编码为 $\pm1$,可以将套索回归应用在一个二类别的分类问题中,并使用一个临界值(通常为 0)来做出类别预测。在多于两个类别的问题中,已有很多方法,包括了第 18.3.3 节中描述的 OVA 和 OVO 方法。我们在第 18.3 节中的癌症数据上尝试了 OVA 方法。表 18.1 的第 4 行展示了结果。它是表现最好的几个模型之一。
分类问题一个更自然的方法是用套索惩罚项来正则化对数几率回归。在相关文献中已经提出了一些实现方法了,包括了与 LARS 类似的路径算法(Park and Hastie, 2007)。由于路径是分段平滑但非线性的,所以精确的方法会比 LARS 算法更慢,并且当 $p$ 比较大时可行性更低。
Friedman et al. (2010) 给出了一个拟合 $L_1$ 惩罚项的对数几率和多项分布回归模型的非常快的算法。他们使用了如第 18.3.2 节中式 18.10 的对称多项分布对数几率回归模型,并对类比式 18.11 的带惩罚对数似然函数最大化。
$$\max_{\{\beta_{0k},\beta_k\in\mathbb{R}^p\}_1^K} \left[ \sum_{i=1}^N \log \operatorname{Pr}(g_i|x_i) - \lambda \sum_{k=1}^K \sum_{j=1}^p |\beta_{kj}| \right]\tag{18.19}$$他们的算法在预先选定的一系列 $\lambda$ 的取值处,通过(第 3.8.6 节)循环坐标下降方法来计算精确的解。并且他们利用了两个性质:当 $p\gg N$ 时解是稀疏的;相邻 $\lambda$ 取值下的解往往会非常相似。表 18.1 的第 7 行使用了这个方法,其中整体的调节参数 $\lambda$ 是由交叉验证选取的。它的表现与最佳的几个方法类似,不过它的自动特征选择一共选中了 269 个基因。Genkin et al. (2007) 使用了一个类似的方法;虽然他们从一个贝叶斯的角度来描述他们的模型,不过他们实际上计算的后验众数就是带惩罚最大似然问题的解。
在基因学的应用中,变量之间通常存在强相关性;基因往往在分子通道(molecular pathway)中表现。套索惩罚项对一组强显著但彼此相关的变量的选择不太敏感(练习 3.28)。而另一方面,岭式惩罚项倾向于将彼此相关变量的系数向彼此的方向收缩(99 页的练习 3.29)。弹性网(elastic net) 惩罚项(Zou and Hastie, 2005)是两者的折中,其形式为:
$$\sum_{j=1}^p (\alpha|\beta_j| + (1-\alpha) \beta_j^2) \tag{18.20}$$第二项促进了高相关性的特征变量的平均化,而第一项倾向于得出这些被平均化的特征变量一个稀疏的解。弹性网惩罚项可以被用在任意的线性模型中,特别是回归和分类问题。
因此上述的多项分布问题加上了弹性网惩罚项后变成:
$$\max_{\{\beta_{0k},\beta_k\in\mathbb{R}^p\}_1^K} \left[ \sum_{i=1}^N \log \operatorname{Pr}(g_i|x_i) - \lambda \sum_{k=1}^K \sum_{j=1}^p (\alpha|\beta_{kj}| + (1-\alpha)\beta_{kj}^2) \right]$$ $$\tag{18.21}$$参数 $\alpha$ 决定了惩罚项的混合结构,通常会以定性的方式预先确定出来。当 $p>N$ 时,弹性网可以得到多于 $N$ 个非零参数,这是它对套索惩罚项的一个潜在优势。表 18.1 的第 8 行使用了这个模型,其中的 $\alpha$ 和 $\lambda$ 通过交叉验证选择。我们使用了从 0.05 到 1.0 之间的 20 个 $\alpha$ 取值,和整体取值域上以对数标尺均匀分布的 100 个 $\lambda$ 的取值组成的参数序列(网格)。在所有两个参数的组合中,$\alpha\in[0.75,0.80]$ 和 $\lambda<0.001$ 得到了最小的交叉验证误差。尽管这在所有的方法中是最小的测试误差,但差别较小而且并不显著。另外有意思的是,当在每个 $\alpha$ 取值上分别地进行交叉验证时,在 $\alpha=0.10$ 时得到了最小的测试误差 8.8,但是这与二维的交叉验证中选择的取值不同。
图 18.5 展示了在两个类别的白血病数据集(Golub et al., 1999)上的套索和弹性网的系数路径。数据中有 38 个样本的 7129 个基因表达测量,其中 27 个的类别为 ALL(急性淋巴细胞白血病),另外 11 个的类别为 AML(急性髓细胞白血病)。也有 34 个样本的一个测试集,类别数量分别为 20 和 14。由于数据是线性可分的,所以在 $\lambda=0$ 处的解是不确定的(练习 18.11),并且当 $\lambda$ 增加很小的值时模型表现就开始变差。所以所示的系数路径是随着拟合的概率接近 0 和 1 时就被截断了。左图中有 19 个非零系数,右图中有 39 个非零系数。
图 18.6 的左图展示了套索对数几率回归在训练集和测试集上的误分类误差,以及在训练集上的 10 折交叉验证误差。右图中使用了二项偏差(损失函数)1来衡量误差,曲线更平滑。尽管每条曲线看起来相对平滑,但小样本量会导致在这些曲线上有相当大的样本方差(sampling variance,参考 220 页的图 7.1)。这两个图都表明了 $\lambda\downarrow0$ 的极限解是足够好的,在测试集上的误分类率是 $3/34$。弹性网对应的图与之在性质上相似,就不再展示了。
当 $p\gg N$ 时,在所有的正则化对数几率回归模型中,极限下的系数是发散的,所以拟合软件在实践中会或显式地或隐含地设置一个 $\lambda>0$ 的最小值。然而,重新标准化后的系数是收敛的,这些极限下的解可被视为是对线性最优分离超平面(SVM)的一个很好的替代。当 $\alpha=0$ 时,极限解与 SVM 完全一致(参考第 18.3.2 节的末尾),但是所有 7129 个基因都被选择了。当 $\alpha=1$ 时,极限解与 $L_1$ 分离超平面(Rosset et al., 2004a)完全一致,它最多包含 38 个基因。随着 $\alpha$ 从 1 开始降低,弹性网的解会在分离超平面包含更多的基因。
18.4.1 套索在蛋白质质谱中的应用
蛋白质质谱(protein mass spectrometry)已经成为分析血液中蛋白质的一个普遍的技术,可以被用来诊断疾病或者理解疾病底层的发展过程。
在每个血清样本 $i$ 中,我们在很多个飞行时间(time of flight)$t_j$ 取值下观测强度 $x_{ij}$。粒子在经过机器的一套操作后从发射器移动到探测器,这个强度(intensity)就与所观测到在这个过程中耗费大约 $t_j$ 时间的粒子数量相关。飞行时间与血液中的成分蛋白质的质荷比(mass over charge ratio,m/z)存在已知的关系。所以在某个 $t_j$ 识别出的质频高峰就代表了存在着相对应质荷比的蛋白质。然后可以通过其他方法来确定这个蛋白质的种类。
图 18.7 展示了来自 Adam et al. (2003) 的一个例子。图中展示了健康患者和前列腺癌症患者的平均质谱图。共有 16,898 个 m/z 的取值位置,取值范围从 2000 到 40,000。全部的数据集中包括了 157 个健康患者和 167 个癌症患者,目的是找出可辨别出两组患者的 m/z 位置。这是一个 函数型数据(functional data) 的例子;自变量可被看成是 m/z 的一个函数。过去的几年中这个问题得到了很多关注;例如可参考 Petricoin et al. (2002)。
首先将数据进行标准化(减去基线以及标准化),并且只关注 m/z 在 2000 到 40,000 区间之内的数据(这个范围外的质谱和问题不相关)。然后在数据上使用最近收缩中心和套索回归,表 18.2 中展示了两个方法的结果。
方法 | 测试误差(/108) | m/z 位置的数量 |
---|---|---|
1. 最近收缩中心 | 34 | 459 |
2. 套索 | 22 | 113 |
3. Lasso on peaks | 28 | 35 |
表 18.2: 前列腺癌症数据例子的结果。测试误差的标准差大约为 4.5。
套索方法通过对数据的强拟合得到了一个相当低的测试误差率。但是它可能并不是一个在科学上有用的解决方案。理想的情况下,蛋白质质谱从一个生物样本中解析出它的蛋白质成分,而这些应该在质谱中对应了高峰的位置。套索方法对高峰位置没有任何的特殊处理方式,所以不出意外地在非零的套索权重中只有一部分在质谱中是位于高峰附近的。此外,同一个蛋白质在不同的质谱中可能是在稍微不同的 m/z 取值处产生高峰。所以为了识别出通用的高峰,需要逐样本地进行某种 m/z 的扭转变换(warping)。
为解决这个问题,我们对每个质谱使用了一个标准的谱峰提取算法,在 217 个质谱训练集中共得到了 5178 个谱峰。我们的思路是将来自所有患者的谱峰汇总到一起,从而构建一组公共的谱峰。为了实现这个目的,我们在这些谱峰沿着 $\log m/z$ 坐标轴的位置上使用了层次聚类(hierarchical clustering)。我们水平地在 $\log(0.005)$ 的高度剪切了所得到的树状图(dendrogram)2,并在每个得出的簇中计算平均的谱峰位置。这个计算过程得出了 728 个通用的簇以及它们对应的谱峰中心。
给定这些 728 个通用的谱峰,我们要判断在每个单独的质谱中出现了哪些谱峰,并且这些谱峰的高度。如果一个谱峰没有出现,那么会赋予为零的谱峰高度。这样就得出了一个作为特征变量的 $217\times728$ 谱峰高度矩阵,可以被用到一个套索回归中。在质谱测试集上也用同样的 728 个谱峰进行处理。
表 18.2 的最后一行展示了这个在谱峰特征上应用套索方法的预测结果;它的表现比较好,但要差于在原始质谱上的套索。但是这个拟合模型可能对生物学家来说是更有用的,因为它得出了 35 个谱峰的位置可供进一步的研究。而另一方面,这个结果意味着质谱中谱峰之间可能仍有一些对类别判断有效的信息,表格第 2 行套索方法中的位置信息可能值的进一步地进行研究。
18.4.2 函数型数据的融合套索
在上一个例子中,特征变量有一个天然的排序,即它们的质荷比 $m/z$。更一般地,可能存在函数式特征变量 $x_i(t)$,它可根据某个索引变量 $t$ 来排序。我们已经介绍过一些针对这种结构的方法。
可以将 $x_i(t)$ 用一组 $t$ 的基函数的系数来代表,比如样条函数、小波函数、或傅立叶基函数,然后将这些系数作为自变量来进行回归。或者等价地,可以用这些基函数来代表原始特征变量的这些系数。第 5.3 节介绍了这种方法。
在分类问题的场景中,第 12.6 节介绍了与带惩罚的判别分析相类比的方法。它使用了一个惩罚项来直接地控制了系数向量的平滑程度。
上述的方法一般是均匀地对系数进行平滑。而这里介绍一个更加自适应的策略,可通过调整套索惩罚项从而把特征变量的次序考虑进去。融合套索(fused lasso) 方法(Tibshirani et al., 2005) 对下式求解:
$$\max_{\beta\in\mathbb{R}^p} \left\{ \sum_{i=1}^N (y_i - \beta_0 - \sum_{j=1}^p x_{ij}\beta_j)^2 + \lambda_1 \sum_{j=1}^p |\beta_j| + \lambda_2 \sum_{j=1}^{p-1} |\beta_{j+1}-\beta_j| \right\}$$ $$\tag{18.22}$$这个准则对 $\beta$ 是严格凸性的,所以存在一个唯一解。第一个惩罚项促使解偏向于稀疏,而第二个惩罚项促使解对索引 $j$ 更平滑。
式 18.22 中的差分惩罚项假设了索引 $j$ 之间的间隔是均匀的。如果其背后索引变量 $t$ 的取值是不均匀的 $t_j$,式 18.22 一个自然的推广可以基于差商:
$$\lambda_2 \sum_{j=1}^{p-1} \frac{|\beta_{j+1} - \beta_j|}{|t_{j+1} - t_j|} \tag{18.23}$$也就是对差分系列中的每一项添加一个惩罚项修改因子。
当自变量矩阵 $\mathbf{X}=\mathbf{I}_N$,$N\times N$ 的单位矩阵时,这就称为了一个很有用的特例。这是融合套索的一个特例,被用来近似一个 $\{y_i\}_1^N$ 序列。融合套索信号逼近器(fused lasso signal approximator) 对下式求解:
$$\max_{\beta\in\mathbb{R}^p} \left\{ \sum_{i=1}^N (y_i - \beta_0 - \beta_i)^2 + \lambda_1 \sum_{i=1}^N |\beta_i| + \lambda_2 \sum_{i=1}^{N-1} |\beta_{i+1}-\beta_i| \right\}$$ $$\tag{18.24}$$图 18.8 展示了来自 Tibshirani and Wang (2007) 的一个例子。图中的数据是来自一个比较基因组杂交(comparative genomic hybridization,CGH)序列,在一个肿瘤样本中测量每个基因复制数量与正常样本相比的近似(以二为底的)对数比率。水平坐标轴代表了每个基因的染色体位置。其想法是,在癌症细胞中基因通常会被放大(复制)或被删除,所以去探测这样的事件是有意义的。此外,这样的事件往往会在一个连续的区域上出现。图中展示的深红色曲线是从融合套索信号逼近器得出的平滑信号估计(其中选用了适当的 $\lambda_1$ 和 $\lambda_2$ 取值)。显著非零的区域可以被用来探测在肿瘤中基因增强或损失的位置。
融合套索也有二维的版本,其中的参数排列在一个像素点网格上,对目标像素点的左侧、右侧、上侧、和下侧四个方向的一阶差分进行惩罚。这个方法可用于对图片的降噪(denoising)和分类。Friedman et al. (2007) 为一维和二维的融合套索开发出了快速的广义坐标下降算法。
本节练习
练习 18.11
假设有 $N$ 对 $(x_i,y_i)$ 的样本,其中 $y_i$ 是二元变量,$x_i\in\mathbb{R}^1$。 同时假设两个类别是可分的;例如,存在某个 $C>0$,使得对任意的 $i$ 和 $i’$ 并且 $y_i=0$ 和 $y_{i’}=1$,则有 $x_{i’}-x_i\geq C$。现在想要通过最大似然方法拟合一个线性的对数几率回归模型
$$\operatorname{logit} \operatorname{Pr}(Y=1 | X) = \alpha + \beta X$$试证明 $\hat{\beta}$ 无法确定。