本节介绍一个在 $p\gg N$ 时特别有效的用于回归和广义回归的简单方法。我们用零一个微阵列数据的例子来演示这个方法。数据来自 Rosewald et al. (2002)，包含了 240 个弥漫大 B 细胞淋巴瘤（DLBCL）患者的样本，样本中有 7399 个基因表达的测量。输出变量是观测到的或者右删失1的生存时间。随机地将淋巴瘤样本分成大小为 160 的训练集和 80 的测试集。
尽管监督主成分可用于线性回归，但它最有吸引力的应用可能在于生存分析，也是这个例子所关注的。
本书中尚未讨论过对删失（censored）生存数据的回归；它是一个回归的一个推广形式，其中一部分个体的输出变量（生存时间）只被部分地观测到。例如假设在一个持续 365 天的医学试验中，简单起见假设所有的被试都是在第一天加入的。可能会观测到某个体在试验开始的 200 天之后死亡，而另一个体可能在 365 天后试验结束时仍然存活。后者就被称为在 365 天“右删失”。我们只知道该个体生存了至少 365 天。虽然并不知道在 365 天之后该个体实际又生存了多长时间，但删失的观测仍然提供了很多信息。图 18.11 展示了删失数据的例子。
图 18.12 展示了在 80 位患者的测试集上通过 Kaplan-Meier 方法估计出的生存曲线。更多关于 Kaplan-Meier 方法的介绍，可参考 Kalbfleisch and Prentice (1980)。
在这个例子中，我们的目标是找出能够预测一组独立的患者集生存时间的特征变量（基因）集合。这可以用作一个预后指示变量来协助治疗方案的选取，或者有助于对疾病的生物学基础的理解。
图 18.13 展示了监督主成分的底层概念模型。假设有两种细胞类型，平均来说携带好细胞类型的患者会生存更长时间。不过在这两组生存时间分布中有很大程度的重叠区域。可以将生存时间看成是对细胞类型的一个“嘈杂替代变量（noisy surrogate）”。一个完全监督方法可以给那些与生存时间关系最强的基因赋予最高的权重。这些基因与细胞类型部分相关，但并不完全相关。患者隐含的细胞类型通常可以被在通路共同起作用的基因的一个较大的特征所反映出来，而如果我们可以发掘出这个细胞类型，那么就可以更好地预测患者的生存时间。
虽然图 18.13 中的细胞类型是离散的（好和坏），但它也可以是由特征变量的某种线性组合所定义出的连续细胞类型。我们会将细胞类型作为一个连续的数量来进行估计，然后在展示和解释中将它离散化。
那么如何才能找到定义这个重要的隐含细胞类型的线性组合呢？主成分分析（第 14.5 节）是在特征变量的线性组合中寻找数据集中变动较大的线性组合的一个很有效的方法。但在此例中我们要寻找的线性组合需要同时具有高方差以及与输出变量有显著的相关性。图 18.14 的右下图展示了在此例中使用标准的主成分的结果；第一主成分与生存时间并没有很强的相关性（具体细节参考图下的介绍）。
所以我们想要促使主成分分析能够找到与输出变量高度相关的特征变量线性组合。为实现这个目的，我们仅考虑那些本身就与输出变量存在相当的相关性的特征变量。算法 18.1，监督主成分（supervised principal components），概括了这个方法；图 18.14 做了一个演示。
算法 18.1：监督主成分
- 分别地用每个特征变量对输出变量做回归，计算它们的标准化单变量回归系数。
- 对一个阈值列表 $0<\theta_1<\theta_2<\dots<\theta_K$ 中的每一个取值 $\theta$：
- 构造一个缩小的数据矩阵，其中只包含单变量系数绝对值超过 $\theta$ 的那些特征变量，然后计算这个矩阵的前 $m$ 个主成分。
- 通过一个使用这些主成分的回归模型来预测输出变量。
- 通过交叉验证来选择 $\theta$（以及 $m$）。
步骤（1）和（2.2）中的操作细节会依赖于输出变量的类型。对于常规的回归问题，在步骤（1）中可使用单变量线性最小二乘法的系数，在步骤（2.2）中可使用一个线性最小二乘模型。在生存时间问题中，Cox 比例风险回归（Cox’s proportional hazards regression）模型被广泛使用；因此在步骤（1）中可使用来自这个模型的得分检验（score test），在步骤（2.2）中可使用多元 Cox 模型。具体的计算细节不影响对这个方法基础的理解；可以参考 Bair et al. (2006)。
图 18.14 展示了此例中监督主成分的结果。模型使用了 3.53 的 Cox 得分阈值点，得到了 27 个基因，这个 3.53 的取值是通过十折交叉验证选取的。然后在只基于这个子集数据来计算第一主成分（$m=1$），以及它在每个测试观测点上的取值。将这个主成分变量作为 Cox 回归模型中的一个数值自变量，它的似然比显著性为 $p=0.005$。利用模型结果对测试样本进行二分（dichotomized，用训练集数据的得分均值作为划分的阈值），可明显看到测试集中的患者被分成了低风险和高风险的两个组（图 18.14 的右中图，$p=0.006$）。
图 18.14 的右上图只用最高得分的基因（二分后）来预测生存时间。结果在测试集上并不显著。类似地，右下图展示了使用所有训练集数据的主成分二分，结果也不显著。在算法过程的步骤（2.1）中可以使用 $m>1$ 个主成分。但步骤（1）中的监督回归促使了主成分与输出变量的方向靠近，所以在大部分场景中往往只有第一个或前几个主成分对预测是有效的2。在下文的数学推导中将只考虑第一主成分，不过到多个主成分的推广可从类似的方法推导出来。
18.6.1 与潜变量模型的关联
通过数据的一个潜变量模型，可以看出监督主成分与隐含细胞类型模型（图 18.13）之间的联系。假设一个输出变量 $Y$ 与一个隐含潜变量 $U$ 呈线性模型：
$$Y = \beta_0 + \beta_1 U + \varepsilon \tag{18.32}$$此外还有一组特征变量 $X_j$ 的测量值，索引 $j\in\mathcal{P}$（P 来自通路 pathway），使得：
$$X_j = \alpha_{0j} + \alpha_{1j} U + \varepsilon_j, j \in \mathcal{P} \tag{18.33}$$假设误差项 $\varepsilon$ 和 $\epsilon_j$ 都为均值零并且独立于它们各自模型中所有其他的随机变量。
还有很多额外的特征变量 $X_k$，$k\notin\mathcal{P}$，它们与 $U$ 独立。我们想要识别出 $\mathcal{P}$、估计出 $U$、并且拟合预测模型（式 18.32）。这是隐结构（latent structure）模型，或者单成分因子分析（single component factor analysis）模型的一个特例（参考 Mardia et al., 1979 以及第 14.7 节）。潜因子 $U$ 在概念上就是图 18.13 中细胞类型变量的连续版本。
监督主成分算法可被理解为是拟合这个模型的一个方法：
- 步骤（1）的变量筛选估计出集合 $\mathcal{P}$。
- 给定了 $\hat{\mathcal{P}}$，步骤（2.1）中的最大主成分估计出潜因子 $U$。
- 最后，步骤（2.2）中的回归拟合估计出模型 18.32 中的系数。
步骤（1）的做法很自然，因为平均来说只有当 $\alpha_{1j}$ 非零时回归模型的系数才会非零。所以这个步骤会选择出特征变量 $j\in\mathcal{P}$。如果假设了误差项 $\epsilon_j$ 是同方差的高斯分布，那么步骤（2.1）也很自然。这样主成分分析就是单因子模型的最大似然估计（Mardia et al., 1979）。步骤（2.2）中的回归则是很容易理解的最后一步。
假设一共有 $p$ 个特征变量，其中 $p_1$ 个属于有意义的特征集合 $\mathcal{P}$。如果 $p$ 和 $p_1$ 都增大但 $p_1$ 相对于 $p$ 较小，那么可以证明（在一些合理的条件下）第一监督主成分对隐含的潜因子是一致的。而常规的第一主成分可能不是一致的，因为它可能被非常多的“噪声”特征变量所污染。
最后，如果监督主成分步骤（1）中所使用阈值得出了非常多的用于计算主成分的特征变量。那么为了便于理解和解释以及在实践中的使用，我们希望可以找出能够近似这个模型的一个特征变量缩小集合。预处理（第 18.6.3 节）就是完成这项工作的一个方法。
18.6.2 与偏最小二乘的关系 😱
监督主成分与偏最小二乘回归（PLS，第 3.5.2 节）有非常紧密的关联。Bair et al. (2006) 发现了监督主成分效果很好的关键在于在步骤（2.1）中过滤掉了干扰的特征变量。偏最小二乘方法给干扰特征变量降低了权重，但并没有把它们过滤掉；结果导致非常多的干扰特征变量可能会影响预测。不过已经出现了一些对偏最小二乘方法的改进，使其与监督主成分比较相似（例如 Brown et al. (1991)、Nadler and Coifman (2005)）。简单来说，先和监督主成分的步骤（1）和（2.1）一样选取特征变量，但是之后在这些特征变量上使用 PLS（而不是主成分）。在下文中将这个方法称为“阈值 PLS”。
阈值 PLS 可被看作是监督主成分的一个嘈杂（noisy）的版本。所以预计在实践中预计它可能不会达到一样的效果。假设所有的变量已经被标准化。第一个 PLS 变量的形式为：
$$\mathbf{z} = \sum_{j\in\mathcal{P}} \langle \mathbf{y}, \mathbf{x}_j \rangle \mathbf{x}_j \tag{18.34}$$它可被理解为是对模型 18.33 中的潜因子 $U$ 的一个估计。与之相比，监督主成分的方向 $\hat{\mathbf{u}}$ 满足：
$$\hat{\mathbf{u}} = \frac{1}{d^2} \sum_{j\in\mathcal{P}} \langle \hat{\mathbf{u}}, \mathbf{x}_j \rangle \mathbf{x}_j \tag{18.35}$$其中的 $d$ 是 $\mathbf{X}_\mathcal{P}$ 的最大奇异值。这是从第一主成分的定义得来的结论。所以阈值 PLS 使用的权重是 $\mathbf{y}$ 与每个特征变量的内积，而监督主成分则是使用这些特征变量来推导出一个“自洽（self consistent）”的估计 $\hat{\mathbf{u}}$。由于在 $\hat{\mathbf{u}}$ 的估计中使用到了很多的特征变量，而不只是单个的输出变量 $\mathbf{y}$，可以预见 $\hat{\mathbf{u}}$ 没有 $\mathbf{z}$ 那么嘈杂。实际上，如果在集合 $\mathcal{P}$ 中有 $p_1$ 个特征变量、$p$ 和 $p_1$ 趋向于无穷大、并且 $p_1/N\to0$，那么利用 Bair et al. (2006) 中的方法可以证明：
$$\begin{align} \mathbf{z} &= \mathbf{u} + O_p(1) \\ \hat{\mathbf{u}} &= \mathbf{u} + O_p(\sqrt{p_1/N}) \tag{18.36}\end{align}$$其中的 $\mathbf{u}$ 是模型 18.32 和 18.33 中真实的（未观测到的）潜变量。
下面我们通过一个模拟例子来从数值上比较这些方法。模拟中有 $N=100$ 个样本和 $p=5000$ 个基因。以下是数据生成的方法：
$$\begin{align} x_{ij} &= \begin{cases} 3 + \epsilon_{ij} & \text{if } i \leq 50 \\ 4 + \epsilon_{ij} & \text{if } i > 50 \end{cases} && j = 1, \dots, 50 \\ x_{ij} &= \begin{cases} 1.5 + \epsilon_{ij} & \text{if } 1 \leq i \leq 25 \text{ or } 51 \leq i \leq 75 \\ 5.5 + \epsilon_{ij} & \text{if } 26 \leq i \leq 50 \text{ or } 76 \leq i \leq 100 \end{cases} && j = 51, \dots, 250 \\ x_{ij} &= \epsilon_{ij} && j = 251, \dots, 5000 \\ y_i &= 2 \cdot \frac{1}{50} \sum_{j=1}^{50} x_{ij} + \varepsilon_i && \end{align}$$ $$\tag{18.37}$$其中的 $\epsilon_{ij}$ 和 $\varepsilon_i$ 是独立的正态分布随机变量，均值为 0，标准差分别为 1 和 1.5。因此在前 50 个基因中，样本 1~50 和 51~100 之间平均相差 1 个单位，而且两组样本之间的这个差别与输出变量 $y$ 是相关的。在接下来的 200 个基因中，样本 1~25 和 51~75 与样本 26~50 和 76~100 之间平均相差 4 个单位，但这个差别与输出变量不想管。其他的基因是噪声特征变量。图 18.15 展示了一个典型的模拟实现的热力图，其中左边是输出变量，右边是前 500 个基因。
按照这个模型进行 100 次模拟，图 18.16 概括了测试误差结果。图中最右侧的结果展示了（常规的）主成分和 PLS 的测试误差；两者都被数据中的噪声特征变量严重地干扰了。监督主成分和阈值 PLS 则在选中特征数量的一个很宽泛的区间上都表现良好，并且前者总比后者的测试误差更低。
尽管此例看起来是专门按监督主成分的优势而设计的，但在其他的模拟和真实数据中它也仍然保持着良好的表现。（Bair et al., 2006）
18.6.3 特征选择的预处理
如图 18.16 中所示，监督主成分可以达到比其他方法更低的测试误差。不过它并不总会得出一个只包含了少数特征变量（基因）的稀疏模型。即使在算法的步骤（1）中调整阈值只保留相对少数的特征变量，被排除的特征变量仍可能与监督主成分有相当大的内积（因此仍可作为一个不错的替代变量）。另外，彼此高度相关的特征变量往往会被同时选中，因此在选中的特征变量集合中可能存在很大程度的冗余。
另一方面，已知套索方法（第 18.4 节和第 3.4.2 节）可以从数据中得出一个稀疏的模型。那么在上文中的模拟例子中，这两个方法的测试误差比较起来会是怎么的呢？
图 18.17 展示了在模型 18.37 的一次模拟实现上，套索回归、监督主成分、和预处理套索回归的测试误差。
图中可见监督主成分（橙色曲线）在模型中包括了大约 50 个特征变量时达到了它的最低误差点，这也是模拟模型背后真实的相关特征个数。虽然（从真实模型来看）基于前 50 个特征变量上的一个线性模型应该是最优的模型，但套索回归（绿色曲线）严重地被太多的噪声特征变量干扰，并且在模型中包含（比监督主成分）更少数量的特征时就开始出现了过拟合。
那么是否可以同时兼得监督主成分的低测试误差以及套索回归的稀疏性呢？这就是 预处理（pre-conditioning） 方法（Paul et al., 2008）想要要达成的目标。在这个方法中，首先为每个训练集的观测样本计算出它的监督主成分预测结果 $\hat{y}_i$（通过交叉验证来选择其中的阈值）。然后在用 $\hat{y}_i$，而不是通常的 $y_i$，作为输出变量进行一个套索回归。在这个套索回归拟合中使用了所有的特征变量，而不只是那些被监督主成分的阈值过滤步骤保留下来的特征变量。方法的思想是先对输出变量进行降噪（denoise），那么套索回归就不会像之前那样被太多的噪声变量严重干扰。图 18.17 中可见预处理方法（紫色曲线）在此例中比较奏效，达到了比常规的套索回归低得多的测试误差，也是（在此例中）与监督主成分一样的误差水平。而且同时它使用了更少的特征变量就达到了这个低误差。应用在原始输出变量上的常规套索回归比预处理的套索回归更快地开始过拟合。（在预处理套索中）由于输出变量已经被降噪处理，就不会有过拟合的问题了。对预处理套索的条件参数的选择通常会更多出于主观的考虑，比如基于模型的简约性（parsimony）。
预处理方法可以被用在很多不同的场景中，所使用的初始化估计方法可以不局限于监督主成分，而所使用的后程处理也不局限于套索回归。更多细节可以参考 Paul et al. (2008)。