18.2 对角线性判别分析和最近收缩中心

第一章第十四章 中介绍过的基因表达序列,是生物学中一个重要的新技术。我们下一个例子中的数据来自一组微阵列实验,形成了一个 2308 个基因(列)和 63 个样本(行)的矩阵。每个表达水平取值是一个对数比率 $\log(R/G)$。$R$ 是目标样本中基因特异性(gene-specific)RNA 的数量,这些 RNA 混杂在微阵列一个特定的(基因特异性)位置上;$G$ 是参照样本中对应的 RNA 的数量。样本来自于在儿童中发现的小圆蓝细胞肿瘤(small round blue cell tumor,SRBCT),它们被分成了主要的四个类型:伯基特淋巴瘤(Burkitt lymphoma,BL)、尤文肉瘤(Ewing’s sarcoma,EWS)、神经母细胞瘤(neuroblastoma,NB)、和横纹肌肉瘤(rhabdomyosarcoma,RMS)。另外有 20 个观测样本的一个测试集。在这里不会再深入地介绍相关的学科背景知识。

因为 $p\gg N$,所以无法在数据上拟合一个完全的线性判别分析(LDA);这里需要进行某种程度的正则化。下面将介绍的方法与第 4.3.1 节中的方法类似,但对它有关键的改动从而达到特征选择的效果。最简单的正则化形式是假设在每个类别内特征之间是独立的,也就是说组内的协方差矩阵是对角矩阵。尽管实际上在类别内的特征之间不太可能是独立的,但是当 $p\gg N$ 时我们也没有足够的数据来估计它们之间的依赖关系。独立性的假设极大地减少了模型的参数个数,通常会得出一个有效的并且可解释的分类器。

所以我们使用 对角协方差(diagonal-covariance)LDA 规则来进行分类。类别 $k$ 的 判别得分(discriminant score)(见 110 页的式 4.12)为:

$$\delta_k(x^*) = - \sum_{j=1}^p \frac{(x_j^* - \bar{x}_{kj})^2}{x_j^2} + 2 \log \pi_k \tag{18.2}$$

其中的 $x^*=(x_1^*,x_2^*,\dots,x_p^*)$ 是一个测试样本点的表达水平取值向量,$s_j$ 是基因 $j$ 的混合组内标准差,$\bar{x}_{kj}=\sum_{i\in C_k}x_{ij}/N_k$ 是类别 $k$ 的 $N_k$ 的样本点中基因 $j$ 的平均值,$C_k$ 是类别 $k$ 的样本点下标集。将 $\tilde{x}=(\bar{x}_{k1},\bar{x}_{k2}\dots,\bar{x}_{kp})^T$ 称为类别 $k$ 的 中心点(centroid)。式 18.2 中的第一部分就是简单的 $x^*$ 与类别 $k$ 中心点之间(负)标准化平方距离。第二部分是基于类别的 先验概率(prior probability) $\pi_k$ 的一个修正,其中的 $\sum_{k=1}^K\pi_k=1$。那么分类的规则就是:

$$C(x^*) = \ell \text{ 如果 } \delta_\ell(x^*) = \max_k \delta_k(x^*) \tag{18.3}$$

可见对角 LDA 分类器等价于一个经过恰当标准化后的最近中心分类器(nearest centroid classifiler)。它也是第 6.6.3 节中介绍的朴素贝叶斯分类器的一个特例。它假设了在每个类别中的特征变量是独立的高斯分布,并且方差一样。

对角 LDA 分类器通常在高维度场景中比较有效。它也被 Bickel and Levina (2004) 称为“独立规则(independence rule)”,并且他们理论上演示了它在高维问题中通常回避标准的 LDA 表现更好。在这里的例子中,对角 LDA 分类器在 20 个测试样本里得到了 5 个误分类的误差。对角 LDA 分类器的一个缺点是它使用了所有的特征变量(基因),因此不方便对模型解释。通过更强的正则化我们可以在测试无差率和可解释性上都继续进行改进。

我们希望用正则化来自动地移除对类别预测没有贡献的特征变量。通过对每个特征分别地将类别上的均值向整体均值进行收缩,我们可以实现特征选择。这个结果是一个正则化版本的最近中心分类器,或者也等同于一个正则化版本的对角协方差 LDA。我们将这个计算过程称为 最近收缩中心(nearest shrunken centroid,NSC)

收缩的计算过程定义如下。令:

$$d_{kj} = \frac{\bar{x}_{kj} - \bar{x}_j}{m_k (s_j + s_0)} \tag{18.4}$$

其中的 $\bar{x}_j$ 是基因 $j$ 总体的均值,$m_k^2=1/N_k-1/N$,$s_0$ 是一个小的正常数,通常会选择为 $s_j$ 的中位数。这个常数是为了防止当 $s_j$ 比较接近于零从而导致非常大的 $d_{kj}$。如果组内方差是一个常数 $\sigma^2$,那么分子中对比的差值 $\bar{x}_{kj}-\bar{x}_j$ 的方差就是 $m_k^2\sigma^2$,所以分母中会采用这种标准化的表达式。 通过软阈值函数,将 $d_{kj}$ 向零收缩:

$$d'_{kj} = \operatorname{sign}(d_{kj}) (|d_{kj}| - \Delta)_+ \tag{18.5}$$
**图 18.02**:橙色线是软阈值函数 $\operatorname{sign}(d_{kj})(|d_{kj}|-\Delta)\_+$,红色线是 45 度直线。
图 18.02:橙色线是软阈值函数 $\operatorname{sign}(d_{kj})(|d_{kj}|-\Delta)_+$,红色线是 45 度直线。

参考图 18.2。其中的 $\Delta$ 是一个待决定的参数;在这个例子中是通过十折交叉验证来估计的(参考图 18.4 中的上图)。每个 $d_{kj}$ 都在绝对值上被缩减了 $\Delta$ 数量,并且当它的绝对值小于 $\Delta$ 时,它被收缩到零。图 18.2 展示了软阈值函数;第 5.9 节在小波系数中也应用了同样的软阈值函数。或者另一种收缩方法是通过硬阈值函数:

$$d'_{kj} = d_{kj} \cdot I(|d_{kj}| \geq \Delta) \tag{18.6}$$

我们更倾向于软阈值,因为它是一个平滑的计算操作并且通常效果更好。然后通过对式 18.4 的转化进行逆运算就可以得出收缩版本的 $\bar{x}_{kj}$:

$$\bar{x}'_{kj} = \bar{x}_j + m_k(s_j + s_0) d'_{kj} \tag{18.7}$$

然后用收缩中心 $\bar{x}’_{kj}$ 替代判别得分函数(式 18.2)中原始的 $\bar{x}_{kj}$。式 18.7 中的估计量可被看成是一个类别均值的套索方式的估计量(练习 18.2)。

需要注意只要至少在一个类别中有非零 $d’_{kj}$ 的那些基因才会存在于分类规则中,所以绝大部分的基因通常会被舍弃掉。在这个例子中,除了 43 个基因外所有都被舍弃掉了,只留下用于描述每个类别的少量易于解释的基因。图 18.3 用一个热力图表现了这些基因。

**图 18.03**:所选的 43 个基因的热力图。在每个水平的分区(行)内,我们通过层次聚类对基因进行了排序;在每个垂直的分区(列)内,也对样本做了类似的处理。黄色代表了过表达,蓝色代表了欠表达。
图 18.03:所选的 43 个基因的热力图。在每个水平的分区(行)内,我们通过层次聚类对基因进行了排序;在每个垂直的分区(列)内,也对样本做了类似的处理。黄色代表了过表达,蓝色代表了欠表达。
**图 18.04**:
图 18.04

图 18.4 的上图演示了收缩的有效性。在没有收缩时,我们在测试数据上有 5/20 个误分类,并且在训练集和验证集上也有一些误分类。在 $\Delta$ 相当宽的一个取值范围内,最近收缩中心都在测试集上没有误分类。图 18.4 的下图展示了 SRBCT 数据中相对于总体中心点的四个(类别的)中心点(灰色)。利用 $\Delta=4.3$ 的软阈值对灰色柱状图收缩,就得到了这些中心点收缩版本的蓝色柱状图。可以使用判别得分函数(式 18.2)来构建类别概率估计:

$$\hat{p}_k(x^*) = \frac{e^{\frac{1}{2} \delta_k(x^*)}} {\sum_{\ell=1}^K e^{\frac{1}{2} \delta_\ell(x^*)}} \tag{18.8}$$

这些概率可以用来评价分类结果,或者决定是否根本不应对某个特定样本进行分类。

在这个场景中也可以使用其他形式的特征选择,其中包含了硬阈值。Fan and Fan (2008) 在理论上展示了在高维问题的对角 LDA 中,进行某种特征选择的重要性。


本节练习

练习 18.2

Nearest shrunken centroids and the lasso

Consider a (naive Bayes) Gaussian model for classification in which the features $j=1,2,\dots,p$ are assumed to be independent within each class $k=1,2,\dots,K$. With observations $i=1,2,\dots,N$ and $C_k$ equal to the set of indices of the $N_k$ observations in class k, we observe $x_{ij}\sim\mathcal{N}(\mu_j+\mu_{jk},\sigma_j^2)$ for $i\in C_k$ with $\sum_{k=1}^K\mu_{jk}=0$. Set $\hat{\sigma}_j^2=s_j^2$, the pooled within-class variance for feature j, and consider the lasso-style minimization problem

$$\min_{\{\mu_j, \mu_{jk}\}} \left\{ \frac{1}{2} \sum_{j=1}^p \sum_{k=1}^K \sum_{i \in C_k} \frac{(x_{ij} - \mu_j - \mu_{jk})^2}{s_j^2} + \lambda \sqrt{N_k} \sum_{j=1}^p \sum_{k=1}^K \frac{|\mu_{jk}|}{s_j} \right\}$$ $$\tag{18.55}$$

Show that the solution is equivalent to the nearest shrunken centroid estimator (18.7), with s0 set to zero, and m2k equal to 1/Nk instead of 1/Nk − 1/N as before.

下一页