基于Bayesian Bootstrap抽样的高维线性回归模型

High Dimensional Linear Regression Model Based on Bayesian Bootstrap Sampling

摘要

研究小样本下高维线性回归模型中的变量选择问题和模型预测能力。当自变量维数p远大于样本量n时，提出基于Bayesian bootstrap抽样的SCAD（smoothly clipped absolute deviation）压缩方法。仿真和实证分析表明，与SCAD和LASSO（least absolute shrinkage and selection operator）两种传统回归压缩方法相比，本算法受随机干扰影响较小。当样本量较小时，本算法的变量压缩结果更好，变量选择能力更强，模型的标准均方误差值也最小，且模型预测能力提升明显。

Abstract

The variable selection is proposed and the prediction ability of high dimensional linear regression model under small sample size is studied. When the dimension of independent variable p is much larger than the sample size n, the SCAD (smoothly clipped absolute deviation) compression method based on Bayesian bootstrap sampling is proposed. Simulation and empirical analysis show that the algorithm is less affected by random interference than SCAD and LASSO (least absolute shrink and selection operator), two traditional regression compression methods. It shows that the smaller the sample size, the better the result of compression of variables, the stronger the ability of variable selection, and the smaller the normalized mean square error of the model, and the prediction ability of the model is improved significantly.

Graphical abstract

关键词

高维线性回归 / 变量选择 / 小样本 / Bayesian bootstrap / LASSO（least absolute shrinkage and selection operator） / SCAD（smoothly clipped absolute deviation）

Key words

high dimensional linear regression / variable selection / small sample / Bayesian bootstrap / LASSO（least absolute shrinkage and selection operator） / SCAD（smoothly clipped absolute deviation）

0　 引 言

在医学、金融学等诸多领域中，经常会遇到高维线性回归模型预测问题，特别是当样本量较小，变量的维数远远超过样本量时，需要重点考虑模型的变量选择问题。现有的回归变量选择方法包括LASSO（least absolute shrinkage and selection operator）、自适应LASSO和SCAD（smoothly clipped absolute deviation）等。LASSO^［1］方法适用于高维、强相关、小样本量的数据，缺点是对所有参数进行相同强度压缩，往往会产生较大偏差，得到有偏估计的结果。自适应LASSO^［2］为LASSO惩罚项添加了一个权重，解决了LASSO估计的有偏问题，但是它不能直接有效解决多重共线性的问题，而且自适应LASSO需要一个初始估计量，通常这个估计量要求是\sqrt[]{n}相合的，在高维问题中这个要求不易满足。SCAD^［3］方法可应用于广义线性模型和强健的回归模型，而且具有Oracle性质，缺点是稳定性较差^［4］。Bayesian bootstrap抽样多用于处理小样本问题^［5］，Li等^［6］将Bayesian bootstrap用于回归模型的自适应LASSO估计，有效地降低了方差，但是变量的维数较小。本文使用Bayesian bootstrap与SCAD相结合的方法处理小样本量的高维数据，解决SCAD方法稳定性差的问题，并提升回归预测效果。

1　 模型与方法

研究一个高维线性回归模型。设因变量\mathit{y}=(y_{\mathrm{1}},y_{\mathrm{2}},\cdots ,y_n{)}^{\mathrm{T}}，y_i是第i个观测值，自变量\mathit{X}=(x_{\mathrm{1}},x_{\mathrm{2}},\cdots ,x_n{)}^{\mathrm{T}}，{\mathit{x}}_i是第i个样本自变量的值， X 是一个n×p维数据矩阵，其中p\gg n。该线性回归模型可以表示为

\mathit{y}=\mathit{X}\mathit{\beta }+k\mathit{\epsilon }（1）

其中，\mathit{\beta }是p维回归系数，\mathit{\epsilon }=({\epsilon }_{\mathrm{1}},{\epsilon }_{\mathrm{2}},\cdots ,{\epsilon }_n{)}^{\mathrm{T}}，{\epsilon }_i是均值为0，方差为1的随机变量，k是随机项系数。

Bayesian bootstrap^［7］是bootstrap的一种贝叶斯模拟。给定权重向量\mathit{\omega }=({\omega }_{\mathrm{1}},{\omega }_{\mathrm{2}},\cdots ,{\omega }_n{)}^{\mathrm{T}}。假设先验分布n^{-\mathrm{1}}\mathit{\omega }\sim \mathrm{D}\mathrm{i}\mathrm{r}\mathrm{i}(\mathit{\alpha })，其中\mathit{\alpha }=({\alpha }_{\mathrm{1}},{\alpha }_{\mathrm{2}},\cdots ,{\alpha }_n{)}^{\mathrm{T}}，每个观测都是独立且随机的，则后验分布n^{-\mathrm{1}}\mathit{\omega }\sim \mathrm{D}\mathrm{i}\mathrm{r}\mathrm{i}(\mathit{\alpha }+\mathrm{1})，\mathrm{1}={(\mathrm{1,1},\cdots ,\mathrm{1})}^{\mathrm{T}}。通常取\mathit{\alpha }=\mathrm{0}，则n^{-\mathrm{1}}\mathit{\omega }的后验分布为\mathrm{D}\mathrm{i}\mathrm{r}\mathrm{i}(\mathrm{1})，渐近后验分布通过蒙特卡洛随机模拟得到。该方法将随机模拟的权重向量与bootstrap的随机抽样相结合，即先对样本进行M次bootstrap随机抽样，与随机权重向量的M次蒙特卡罗抽样相结合，给每一个bootstrap随机抽样的样本分配一组分布为n\mathrm{D}\mathrm{i}\mathrm{r}\mathrm{i}(\mathrm{1})的随机权重，使得样本量从固定样本量变成一个更大的后验样本量。Newton 等^［8］和Lyddon 等^［9］得到渐近分布性质。

假设F_{\mathrm{0}}(\mathit{x})是总体的理论分布函数，理论参数{\mathit{\beta }}_{\mathrm{0}}可由下式估计

{\widehat{\mathit{\beta }}}_{\mathrm{0}}=\underset{\mathit{\beta }}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}\int L(\mathit{\beta };\mathit{x})\mathrm{d}{F}_{\mathrm{0}}(\mathit{x})（2）

其中L(\mathit{\beta };\mathit{x})是损失函数。因为F_{\mathrm{0}}(\mathit{x})未知，所以通过经验分布函数F_n^{(m)}(\mathit{x})来做估计，其中m=\mathrm{1,2},\cdots ,M。经验分布函数可以表示为

（3）

其中，{\omega }_i^{(m)}>\mathrm{0}且{\displaystyle \sum _{i=\mathrm{1}}^n}{\omega }_i^{(m)}=n，I(\cdot )为示性函数。

第m次重抽样后，\mathit{\beta }的Bayesian bootstrap的参数可以通过下式进行估计

{\widehat{\mathit{\beta }}}^{(m)}=\underset{\mathit{\beta }}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}\int L(\mathit{\beta };\mathit{x})\mathrm{d}{F}_n^{(m)}(\mathit{x})=

\underset{\mathit{\beta }}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}\int {‖\mathit{y}-\mathit{X}\mathit{\beta }‖}^{\mathrm{2}}\mathrm{d}{F}_n^{(m)}(\mathit{x})（4）

由于传统的最小二乘不能进行系数压缩，所以无法筛选出非零变量，就容易导致过拟合和高方差的问题，因此考虑引入带有权重的惩罚项p_{\lambda }({\beta }_j)，参数估计形式变为

{\widehat{\mathit{\beta }}}^{(m)}=\underset{\mathit{\beta }}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}\left\{\int {‖\mathit{y}-\mathit{X}\mathit{\beta }‖}^{\mathrm{2}}\mathrm{d}{F}_n^{(m)}(\mathit{x})+{\displaystyle \sum _{j=\mathrm{1}}^p}{p}_{\lambda }\left({\beta }_j\right)\right\}（5）

若p_{\lambda }({\beta }_j)=\lambda \left|{\beta }_j\right|，LASSO方法通过（5）式中的一阶惩罚项将部分回归系数压缩为零，从而实现变量选择。

若p_{\lambda }({\beta }_j)满足

（6）

SCAD方法通过（6）式中的惩罚项将部分回归系数压缩为零，其中a>\mathrm{2}为常数，是预先给定的参数。相比于LASSO方法的，SCAD方法具有无偏性，但计算量比LASSO大。

徐国盛等^［4］的数值模拟发现LASSO和SCAD方法的变量选择性能相当，但SCAD的稳定性要比LASSO方法差。Li等^［6］的实验结果表明，Bayesian bootstrap方法可以有效降低方差。因此为增强SCAD方法的稳定性，本文使用Bayesian bootstrap和SCAD结合的方法降低方差，进而提高SCAD的预测性能。Bayesian bootstrap SCAD（BBSCAD）的参数估计为^［6］

{\widehat{\mathit{\beta }}}^{\mathrm{B}\mathrm{B}\mathrm{S}\mathrm{C}\mathrm{A}\mathrm{D}}=

\underset{\mathit{\beta }}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}\left\{\int {‖\mathit{y}-\mathit{X}\mathit{\beta }‖}^{\mathrm{2}}\mathrm{d}{F}_n+n{\displaystyle \sum _{j=\mathrm{1}}^p}{p}_{\lambda }\left(\left|{\beta }_j\right|\right)\right\}=

\underset{\mathit{\beta }}{\mathrm{a}\mathrm{r}\mathrm{g}\mathrm{m}\mathrm{i}\mathrm{n}}\left\{{\displaystyle \sum _{i=\mathrm{1}}^n}{\omega }_i{\left(y_i-{\mathit{X}}_i^{\mathrm{T}}\mathit{\beta }\right)}^{\mathrm{2}}+n{\displaystyle \sum _{j=\mathrm{1}}^p}{p}_{\lambda }\left(\left|{\beta }_j\right|\right)\right\}

（7）

算法1给出了一次BBSCAD重复的详细算法：

通过BBSCAD算法可得到M个采样参数估计：\left({\widehat{\mathit{\beta }}}^{\mathrm{B}\mathrm{B}\mathrm{S}\mathrm{C}\mathrm{A}\mathrm{D}(\mathrm{1})},{\widehat{\mathit{\beta }}}^{\mathrm{B}\mathrm{B}\mathrm{S}\mathrm{C}\mathrm{A}\mathrm{D}(\mathrm{2})},\cdots ,{\widehat{\mathit{\beta }}}^{\mathrm{B}\mathrm{B}\mathrm{S}\mathrm{C}\mathrm{A}\mathrm{D}(M)}\right)，其中{\widehat{\mathit{\beta }}}^{\mathrm{B}\mathrm{B}\mathrm{S}\mathrm{C}\mathrm{A}\mathrm{D}(m)}=\left({\widehat{\beta }}_{\mathrm{1}}^{\mathrm{B}\mathrm{B}\mathrm{S}\mathrm{C}\mathrm{A}\mathrm{D}(m)},{\widehat{\beta }}_{\mathrm{2}}^{\mathrm{B}\mathrm{B}\mathrm{S}\mathrm{C}\mathrm{A}\mathrm{D}(m)},\cdots ,{\widehat{\beta }}_p^{\mathrm{B}\mathrm{B}\mathrm{S}\mathrm{C}\mathrm{A}\mathrm{D}(m)}\right)，m=\mathrm{1,2},\cdots ,M。这里有两种方法可以得到最终的回归系数估计：方法一是对M个参数估计取中位数；方法二是对M个参数估计取平均值。取平均值的方法如下

{\widehat{\beta }}_i=\frac{\mathrm{1}}{M}{\displaystyle \sum _{j=\mathrm{1}}^M}{\widehat{\beta }}_i^{\mathrm{B}\mathrm{B}\mathrm{S}\mathrm{C}\mathrm{A}\mathrm{D}(j)},i=\mathrm{1,2},\cdots ,p（8）

可得最终的参数估计为

\widehat{\mathit{\beta }}=\left({\widehat{\beta }}_{\mathrm{1}},{\widehat{\beta }}_{\mathrm{2}},\cdots ,{\widehat{\beta }}_p\right)（9）

但由于高维数据的稀疏性、抽样的随机性和添加权重的随机性，导致每次获得的变量特征呈现多样化，因此通过添加阈值的方法将系数估计平均值较小的剔除。为了使模型最优，通过C_p统计量的方法确定阈值大小，实现变量选择。由于中位数和平均数都是描述数据集中趋势的统计量，反映数据的一般水平，又因为每个样本随机权重的分布为n\mathrm{D}\mathrm{i}\mathrm{r}\mathrm{i}(\mathrm{1})，所以回归系数的中位数和平均数相差不大，不会对最终结果产生过大影响。

选择NMSE^［10］（normalized mean square error）作为模型评价的标准，定义

\begin{array}{l}\mathrm{N}\mathrm{M}\mathrm{S}\mathrm{E}=\overline{{(y-\widehat{y})}^{\mathrm{2}}}/\overline{{(y-\overline{y})}^{\mathrm{2}}}=\\ \sum {(y-\widehat{y})}^{\mathrm{2}}/\sum {(y-\overline{y})}^{\mathrm{2}}\end{array}（10）

其中，\overline{y}为因变量的均值，\widehat{y}为训练集得到的模型对测试集数据的预测值。

NMSE不仅可以作为多个模型的横向比较标准（同一个测试集NMSE越小的模型预测性能越好），而且还可判断该模型对测试集的预测效果（即如果不用任何的模型只用均值来预测，令\widehat{y}=\overline{y}，那么NMSE等于1，因此如果用训练模型预测得到的NMSE比1大，那么就可以判断该模型的预测效果比较差）。

2　 数值模拟

模拟两个数据集，分别用LASSO、SCAD和BBSCAD方法处理数据， X 和 ε 是标准正态的。为保证各个回归系数之间的相关性，令相关系数为{\rho }^{\left|i-j\right|}，其中\rho =\mathrm{0.5}，对于BBSCAD方法的M个回归系数估计，这里采用取中位数的方法得到最终的回归系数估计。

模拟数据集1：令前三维的变量系数非0，其他系数都为0，即\mathit{\beta }=(-\mathrm{9,7},\mathrm{11,0},\mathrm{0},\dots ,\mathrm{0})，维数p=\mathrm{100}，bootstrap抽样和随机权重抽样次数M=\mathrm{100}，变量选择个数为d。令样本量n=\mathrm{30}，随机项系数k=\mathrm{0.5}，然后保持样本量不变，将随机干扰系数k增加到1，最后保持k不变，把样本量增加到50。运用LASSO、SCAD、BBSCAD分别进行变量选择，判断它们正确选择变量的能力。

由图1可见，在变量选择过程中，LASSO方法以均方误差最小为标准，SCAD方法和BBSCAD方法则选择交叉验证误差最小。图1中阴影部分表示对应\lambda估计值约68%的置信区间。SCAD和BBSCAD方法倾向于选择3个变量，LASSO方法倾向于选择7个变量。相比于实际上的3个变量，LASSO方法多选择了4个变量。具体数据结果见表1。由表1可见：在样本量为30，随机干扰系数k为0.5时，BBSCAD方法的NMSE为0.001 032，小于SCAD方法的0.002 283和LASSO方法的0.001 333，三种方法都正确地选择了3个变量；当样本量不变，增大k为1.0时，BBSCAD方法的NMSE为0.003 982，小于SCAD方法的0.005 583，LASSO方法的NMSE为0.003 974，虽然略小于BBSCAD方法的NMSE，但是BBSCAD和SCAD方法正确地选择了3个变量，LASSO方法则多选择了4个变量；当样本量n=\mathrm{50}，k为1.0时，BBSCAD方法的NMSE为0.004 888，小于SCAD方法的0.009 991，LASSO方法的NMSE为0.004 676，虽然略小于BBSCAD方法的NMSE，但是BBSCAD和SCAD方法正确的选择了3个变量，LASSO方法则多选择了5个变量。

模拟数据集2：增加变量维数，保持前三维变量系数非0，其他系数都为0，即令\mathit{\beta }=(-\mathrm{9,7},\mathrm{11,0},\mathrm{0},\dots ,\mathrm{0})，维数p=\mathrm{1}\mathrm{000}，bootstrap抽样和随机权重抽样次数M=\mathrm{100}，变量选择个数为d，令样本量n=\mathrm{40}，k=\mathrm{1}.0，然后将样本量增加到70、100。分别运用LASSO、SCAD和BBSCAD进行回归预测，对比结果如表2所示。

分析表2，在p\gg n时，当n=\mathrm{40},k=\mathrm{1.0}时，BBSCAD方法的NMSE为0.006 009，小于SCAD方法的0.014 807和LASSO方法的0.022 161；变量选择个数方面，BBSCAD方法和SCAD方法都能正确地选择3个变量，LASSO方法则多选择了11个变量。当n=\mathrm{70},k=\mathrm{1.0}时，BBSCAD方法的NMSE为0.005 452，小于SCAD方法的0.008 766和LASSO方法的0.010 234；变量选择个数方面，BBSCAD方法和SCAD方法都能正确的选择3个变量，LASSO方法多选择了2个变量。当n=\mathrm{100},k=\mathrm{1.0}时，BBSCAD方法的NMSE为0.007 413，小于SCAD方法的0.016 994，LASSO方法的NMSE为0.007 190；变量选择个数方面，BBSCAD方法和SCAD方法都能正确的选择3个变量，LASSO方法多选择了8个变量。

综合表1和表2，在高维小样本线性回归模型假设下，BBSCAD方法继承了SCAD方法变量选择方面的优势，在模型预测方面相比SCAD方法提高明显。当样本量很小时，无论是变量选择还是模型预测，BBSCAD相比其他两种方法均表现出较好的性质；当随机干扰系数增大时，BBSCAD方法变量选择个数不变，而且相比SCAD方法，NMSE提高依然明显，LASSO方法相比BBSCAD方法虽然在NMSE方面略小，但变量选择个数影响明显，代价比较大；随着样本量的增大，LASSO方法逐渐表现出预测能力的优势，但相比于选择变量个数的代价，预测能力的优势并不明显。

3　 实例分析

通过实例数据，比较BBSCAD、SCAD和LASSO方法的变量选择性能和预测性能，数据集来自UCI数据库（http：//archive.ics.uci.edu/ml/datasets/Relative+location+of+CT+slices+on+axial+axis）。因变量为人体的CT切片在轴向上的相对位置，自变量是从CT图像中提取的384个特征。

鉴于研究p\gg n的情形，分别选择样本量为100、150、200作为训练集，剩下的样本作为验证集，bootstrap抽样和随机权重抽样次数M=\mathrm{500}。在计算BBSCAD方法的最终回归系数估计时，由于实例数据的复杂性和小样本问题，会使得很多特征在重抽样过程中表现不显著，从而导致过多的回归系数为0，因此采用对500个系数估计取平均值的方法，再利用C_p统计量确定阈值，对均值较小的回归系数估计进行筛选。

图2显示的是当n=\mathrm{200}时，三种方法的变量选择情况，图中阴影部分表示对应于\lambda估计值为68%的置信区间。由图2可知，BBSCAD方法表现出变量选择方面的优势，把384个变量压缩到13个，SCAD方法和LASSO方法的变量选择个数都介于27~74之间。

从表3可知，当n=\mathrm{100}，BBSCAD方法的NMSE值为0.550 023 2，比SCAD方法的0.781 555 4和LASSO方法的0.556 694 4都小，SCAD方法选择的变量个数最少。为了方便比较，调整BBSCAD方法的阈值，当变量选择个数为11时，NMSE变为0.556 690 3，依然是最小的。当n=\mathrm{150}时，BBSCAD方法的NMSE值为0.432 989 4，比SCAD方法的0.589 560 2小，LASSO方法的NMSE为0.286 142 7，虽然LASSO方法的NMSE最小，但它选择了39个变量，相比于SCAD和BBSCAD方法的11个变量，代价比较大。当n=\mathrm{200}时，BBSCAD方法的NMSE值为0.391 053 2，变量选择的个数为15。为了方便比较，调整阈值，令BBSCAD阈值减小到1.6，NMSE变为0.334 898 8，变量选择个数为31，相比SCAD方法的NMSE值0.362 713 8、变量选择个数32，从整体水平来看BBSCAD方法更优，虽然LASSO方法的NMSE为0.307 484 9，但该方法选择了36个变量。

综合来看BBSCAD方法倾向于选择更简单的模型，特别是在p\gg n时，BBSCAD方法的NMSE的值表现最优，选择变量个数也相对较少，随着样本量的增大，LASSO和SCAD方法的NMSE逐渐提高，但付出了选择更复杂模型的代价，BBSCAD方法则依然保持着对变量的高压缩，预测能力也稳步提高。

4　 结 语

利用Bayesian bootstrap方法与SCAD方法相结合得到一个新的BBSCAD方法，通过仿真模拟和实例数据分析，在p\gg n时，BBSCAD方法利用M次重抽样估计模拟回归系数的分布，无论是在模型预测能力还是在变量选择方面，相比SCAD和LASSO方法都有提高。综合来看，BBSCAD适合处理高维小样本线性数据，且效果良好。

参考文献

[1] TIBSHIRANI R. Regression shrinkage and selection via the LASSO ［J］. Journal of the Royal Statistical Society： Series B Methodological， 1996， 58（1）：267-288. DOI： 10.1111/j.1467-9868.2011.00771.x .

[2] ZOU H. The adaptive LASSO and its oracle properties ［J］. Journal of the American Statistical Association， 2006， 101（476）：1418-1429. DOI： 10.1198/016214506000000735 .

[3] FAN J Q， LI R Z. Variable selection via nonconcave penalized likelihood and its oracle properties ［J］. Journal of the American Statistical Association， 2001， 456（96）：1348-1360. DOI： 10.1198/016214501753382273 .

[4] 徐国盛，赵晓兵.变量选择方法在医疗保险赔付评估中的应用［J］.统计与信息论坛，2014，29（11）：59-64. DOI： 1007-3116（2014）11-0059-06 .

[5] XU G S， ZHAO X B. Application of variable selection in estimation of compensation for medical insurance via penalized functions ［J］. Statistics & Information Forum， 2014， 29（11）：59-64. DOI： 1007-3116（2014）11-0059-06（Ch ）.

[6] 张海如，欧阳缮，王国富，等.基于Bayes Bootstrap统计降噪方法的磁共振测深信号检测［J］.中南大学学报（自然科学版），2014，45（9）：3144-3149.

[7] ZHANG H R， OU Y S， WANG G F， et al. Magnetic resonance sounding signal detection based on statistical noise reduction method of Bayes Bootstrap ［J］. Journal of Central South University （Science and Technology）， 2014， 45（9）：3144-3149 （Ch）.

[8] LI B H， WU J. Bayesian bootstrap adaptive lasso estimators of regression models ［J］. Journal of Statistical Computation and Simulation， 2021， 91（2）：1-30. DOI： 10.1080/00949655.2020.1865959

[9] RUBIN D B. The Bayesian bootstrap ［J］. The Annals of Statistics， 1981， 9（1）：130-134. DOI： 10.1214/aos/1176345338 .

[10] NEWTON M， RAFTERY A. Approximate Bayesian inference by the weighted likelihood bootstrap （with discussion） ［J］. Journal of the Royal Statistical Society Series B Methodological， 1994， 56（1）：3-48. DOI： 10.1111/j.2517-6161.1994.tb01956.x .

[11] LYDDON S， HOLMES C， WALKER S. General Bayesian updating and the loss-likelihood bootstrap ［J］. Biometrika， 2019， 106（2）：465-478. DOI： 10.1093/biomet/asz006 .

[12] 吴喜之.复杂数据统计方法——基于R的应用［M］. 第三版. 北京：中国人民大学出版社， 2015：34-42.

[13] WU X Z. Statistical methods of complex data—Application based on R ［M］. 3nd Ed. Beijing： China Renmin University Press， 2015：34-42（Ch）.

基金资助

国家自然科学基金(41972319)

中国高等教育学会理科教育专业委员会高等理科教育研究课题(20ZSLKJYYB32)

华中科技大学教学研究项目(2020100)