在这个大数据爆发的时代, 数据的收集成本不断降低, 高维数据在各个领域随处可见, 这样的数据不仅规模大, 而且关系复杂、相关性强, 传统的统计方法不再适用, 从而给统计学带来巨大的“压力”。那么怎样从冗余的数据中筛选出对我们有用的信息则变成了重中之重, 也是目前研究的热点方向。

处理高维数据的关键点在于进行降维。主成分分析是一种能够处理相关数据和进行数据降维的统计工具, 能够广泛应用于社会科学、医学、计量学等领域中^[1-3]。它将许多原始变量替换为少数主成分, 进而减少变量数量。但每个主成分都是原始变量的线性组合, 那么每一个主成分到底表达的是什么方面的特征就很难看出来, 使得分析结果通常难以解释。在大数据环境下, 该问题受到了很多关注, 比如针对原始变量的替代方案或修改, 以及与有偏回归进行联系。例如, 引入PCA^[4]的稀疏替代方案, Zou^[5]利用Lasso^[6]和elastic net^[7]的稀疏特性, 使得一些载荷能够缩小到0从而使变量稀疏化, 所得方法称为稀疏PCA(SPCA)。或是将主成分同回归进行有机结合, 主成分回归(PCR)^[8-9]是一种非常有吸引力的方法, 它能够解决数据共线性的问题, 并减少变量数量。在PCR中每个主成分都是原始变量的线性组合, 不再是原来的变量。Kawano等^[10]使用SPCA替代传统的主成分回归方法, 这被称为稀疏主成分回归(SPCR)。这不仅可以获得原始变量的稀疏主成分载荷同时还可以得到主成分的数量；张文明等^[11]在主成分分析的基础上, 利用elastic net对主成分的系数进行稀疏近似, 在保留主成分优点的基础上, 还能够提高模型的解释性；Tay等^[12]将主成分与Lasso回归相结合提出了一种新的监督学习方法, 即pcLasso。它将L1惩罚与二次惩罚相结合, 将系数向量缩小到特征值的主要主成分上, 通过模拟表明能够提高预测性能。

本文梳理了传统的PCA、SPCA、PCR的方法, 并与一些稀疏回归方法进行结合, 基于模拟来进行变量选择。利用敏感性与特异性等指标来分析比较不同方法的使用情况, 并扩展到高维数据情形下, 扩宽了模型的应用范围。

1 模型简介

先考虑一个最传统的线性回归模型：

$ y=X \xi+\varepsilon , $

(1)

其中y=(y₁, …, y_n)^T是具有np个观测值的解释变量; x=(x₁, …, x_p) 是n×p数据矩阵, 包含有n个解释变量的信息；ε为随机误差项, 假定相互独立且服从正态分布ε_t~N(0, σ²) 式(1)的最小二乘可以表示为：

$ \operatorname{argmin}\left[\sum_{i=1}^n\left(y_t-\sum_{j=1}^p \xi_j x_{i j}\right)^2\right], $

(2)

最小二乘估计是常用的一种系数估计方法, 在假设条件下, 具有无偏性。当自变量之间存在多重共线性时, 经典的OLS方法变得不再适用。为了解决这些问题, 主成分回归是一种常用的方法, 通过生成自变量的线性组合来解决多重共线性的问题, X分解为:

$ X=\boldsymbol{T P}^{\mathrm{T}}+E, $

(3)

其中E是误差项, P为载荷矩阵。将其带入原来的线性回归模型中, 结果如下：

$ y=\left(\boldsymbol{T P}^{\mathrm{T}}\right) \xi+\varepsilon_T=\boldsymbol{T} \boldsymbol{\theta}+\varepsilon_T, $

(4)

其中θ=P^Tβ, ε_T为误差项, 将OLS估计量应用以上的回归系数。

$ \hat{\theta}=\left(\boldsymbol{T}^{\mathrm{T}} \boldsymbol{T}\right)^{-1} \boldsymbol{T}^{\mathrm{T}} \boldsymbol{y}, $

得到的新的回归量, 也就是所谓的主成分, 不是原始的X的变量, 估计量$\hat{\theta}$对应的是PCA分数的空间, 那么原始数据中的最终回归系数由以下式子给出：

$ \dot{\xi}=p \hat{\theta} \text { 。} $

尽管, 主成分回归是分析多变量数据的有用工具, 但如果响应变量取决于具有小特征值的主成分, 那么该方法可能没有足够的预测精度。稀疏主成分(SPCA)能够获得与响应变量相关的稀疏主成分载荷, 并同时选择主成分的数量。

$ \underset{A, B}{\operatorname{argmin}}\left\{\sum\limits_{i=1}^n\left\|\boldsymbol{x}_t^{\mathrm{T}}-\boldsymbol{A} \boldsymbol{B}^{\mathrm{T}} x_t\right\|^2+\lambda \sum\limits_{j=1}^k\left\|\xi_j\right\|^2+\sum\limits_{j=1}^k \lambda_{i, j}\left\|\xi_j\right\|_1\right\}, \boldsymbol{A}^{\mathrm{T}} \boldsymbol{A}=I_k, $

(5)

其中B=(β₁, …, β_k) 是p×k的载荷矩阵, k表示主成分的数量, A=(α₁, …, α_k) 是一个p×k的矩阵, λ_{i, j}(j=1, …, k) 是正则化参数, ‖·‖₁是ξ的L₁范数。

对于回归模型来说, 另外一个具有解决问题的方案是Kawano(2015)提出的稀疏PCR(SPCR), 他同样能够改善当响应变量对应小特征值而无法预测的情况。它通过下式定义：

$ \underset{A, \xi, \xi, \xi_0}{\operatorname{argmin}}\left\{\begin{array}{l} (1-\omega) \sum\limits_{i=1}^n\left(y_1-\xi_0-\boldsymbol{\xi}^{\mathrm{T}} \boldsymbol{B}^{\mathrm{T}} x_i\right)^2+\omega \sum\limits_{i=1}^n\left\|\boldsymbol{x}_i^{\mathrm{T}}-\boldsymbol{A} \boldsymbol{B}^{\mathrm{T}} x_i\right\|^2 \\ +\lambda_\beta(1-\xi) \sum\limits_{j=1}^k\left\|\xi_j\right\|_1+\lambda_\beta \xi \sum\limits_{j=1}^k\left\|\beta_j\right\|^2+\lambda_\gamma\|\xi\|_1 \end{array}\right\}, \boldsymbol{A}^{\mathrm{T}} \boldsymbol{A}=I_k, $

(6)

其中ξ₀是截距, λ_β和λ_ξ是非负的正则化参数, w和ξ是介于0和1之间的调整参数。公式中的第一项表示相应变量与主成分B^Tx之间的最小二乘损失, 第二项代表数据X的主成分丢失, 调整参数w是用来权衡第一项和第二项的权数, 可以根据目的来选择w的值, 要是想要获得更高的预测精度时, 可以选择较小的w值, 如果想要获得主成分载荷的精确形式时, 可以使用较大的w值。第三项和第五项则是用来表示B和γ的稀疏程度。B上的稀疏性能够轻松进行主成分载荷的选择, γ的稀疏使得主成分个数选择更加方便。调整参数ξ用来权衡B的L₁和L₂范数之间权重。

由此来看, 可以利用一些稀疏的估计量来代替OLS来改进预测效果, 如使用Lasso回归来构建回归模型, 其中一些大的主成分对应大的特征值, 这可以看作是PCA与Lasso的组合, 同样, 可以将其推广到PCA与ALasso, 这样, 所得到的的估计量成为pcLasso、pcALasso。同样的, 使用稀疏主成分来来代替传统的主成分, 再将其推广至稀疏回归。所得到的估计称为spcLasso、spcALasso。

2 随机模拟

这部分利用模拟研究来展示所用方法的性能, 情形1和情形2适用于低维情况, 情形3和情形4则用于高维情况。模拟的方法和参数选择, 参照Kawano所提供的方法。

2.1 低维数据

通过以下几种按情况响应变量由下式计算：y=χξ+ε, 其中误差项ε来自于均值为0, 方差为σ的正态分布。数据矩阵X由多元正态分布N(0, ∑)生成, ξ是真实的系数参数。

情形1:(a)令p=10, 另外∑=I₁₀, 设回归系数真实值为β=(2, 1, 0, …, 0);(b)∑=diag{1, 3², 1, …, 1}, 设回归系数的真实值为β=(8, 1, 0, …, 0)取样本量为50和100。

情形2：在这个例子中, 令β=(-1, 0, 1, 1, 0, -1, -1, 0, 1, 1), x_i和x_j之间的相关性ρ=0.5, 0.9, 分别用来表示信息变量之间的中度相关和高度相关。有∑=ρ^|j-k|, 1≤j, k≤p, 取样本量为50和100。

为了评估以上不同的估计量, 这些评估是通过100次重复试验所计算出来的, 选择以下几种方法作为评价的性能指标。

$ \operatorname{MSE}(\hat{\xi})=\frac{1}{n} \sum\limits_{i=1}^n\left(y_i-\hat{y}_i\right)^2 \text { 。} $

使用敏感性和特异性两个指标来评价回归模型中参数选择的好坏, 敏感性(Sensitivity)和特异性(Specificity)的定义如下：

$ \text { Sensitivity }=(\hat{\xi}, \xi)=\frac{\#\left\{(i, j): \hat{\xi}_{i j} \neq 0, \xi_{i j} \neq 0\right\}}{\#\left\{(i, j): \xi_{i j} \neq 0\right\}}, $

$ \text { Specificity }=(\hat{\xi}, \xi)=\frac{\#\left\{(i, j): \hat{\xi}_{i j}=0, \xi_{i j}=0\right\}}{\#\left\{(i, j): \xi_{i j}=0\right\}}, $

其中#代表计数, Sensitivity∈[0, 1], Specificity∈[0, 1], 值越接近于1, 变量选择效果越好。

表 1展示了在低维情形下的模拟结果, 从中可以看出: (1)当样本量不断增大时, 模型的估计效率就越好; (2)在情形1下, 提取的主成分为2时, SPCR的表现最佳, 不论是在n=50或是在n=200的情况下都有不错的估计精度, 在提取的主成分越多时, spcLasso则具有更高的精确度。(3)在情形2下, 当变量间的共线程度不断增大时, spcALasso结合了SPCR和自适应Lasso回归的优点, 在变量间高度相关下能够提高模型的精度, 且比SPCR方法具有优势。(4)随着样本量的增大, pcLasso和pcALasso也表现出了不俗的估计效率。由于低维情形下各模型的敏感性和特异性差异不大, 故不作展示。

2.2 高维数据

当自变量的个数大于样本量的个数(即p>n)时, 即为高维数据, 在上部分讨论了低维情形下的几种模型的变量选择情况, 想要提高预测精度, 只需不断增大样本量, 就可以提高预测的精度。而在实际问题中, 经常出现p→∞的情况, 如医学问题, 维数通常会远大于样本量, 这时消除无关基因, 寻找致病因子, 确定某种癌变细胞的基因组合, 从而寻找并制定合适的医疗方案就变得尤为重要, 这就是接下来将要讨论的高维数据问题。

假设变量来自传统的的线性回归模型y=χξ+ε, 其中误差项ε来自于均值为0, 方差为1的正态分布。用ρ^|i-j|表示任意两个解释变量之间的相关系数, 分别取ρ=0.5和ρ=0.9, 样本量n=30和70, 重复进行100次试验。考虑以下几种情形：

情形3：样本量n=30, 回归系数为：$\xi=\underbrace{3, 3, 3, 3, 0, 0, 0, 0, 0, 0, 0, 0, 0, 0}_{\times 5}$此时p>n。

情形4：样本n=70, 回归系数与上情形一样, 此时p=n。

对于每一个p≥n组合, 分别计算以上两种情形下模型的敏感性、特异性以及MSE, 结果见表 2。

通过模拟, 表 2列出了在情形3和情形4中所描述的高维数据场景的三个指标结果。从结果中看：第一, 在样本量不变时, 提取的主成分越少时, SPCR能够准确提取变量的比例会高于其他几种模型。另外，增大样本量, 敏感性和特异性也会增大；第二, 当维数和相关系数ρ不变时, 增大样本量增加, 模型的均方误差也会减小；第三, 对于模型能够正确判别数据的比率, 在不同的情况下会有不同的结果。在样本量越大, 提取的主成分越多时, 数据间的相关性越高, 模型能够正确选择变量的比率就越高, 可以看到敏感性和特异性的值越接近于1；第四, 相关系数增大时, 模型的MSE值有所减少, 在这种高度相关的情况下, L₂惩罚会变得不稳定, 共线问题会导致Lasso预测能力的降低, 自适应Lasso具有同样的问题, 但由于加入了SPC, 这减少了变量间的共线性, 明显提高了预测精度。

3 结论

通过模拟结果, 能够发现在低维模型中, 当提取的主成分较少时, 且变量间的相关程度不高时, SPCR的预测精度会优于其他的模型；另外, 当提取的主成分较多时, spcLasso则展现出了它的优势。当维数不断增大, 到达p≥n时, 能够发现：PCR适用于提取主成分较少时, spcALasso适用于主成分较多时；增大相关系数时, 我们所选择的模型MSE值都有所减少, 这说明当变量间相关性很高时, 本文所提出的方法能够很好的提高模型的预测精度。需要强调的是, 本文并未表明某种模型具有绝对优势, 是为了说明不同模型适用于不同的数据类型。