第一节
LASSO

一、LASSO降维思想和原理

（一）LASSO降维思想

在样本量小、训练数据不足时，经常会出现过拟合（over-fitting）现象，从而导致模型对训练集外的数据拟合效果较差。正则化是机器学习中最常见的过拟合解决方法，在损失函数中加入正则项来惩罚模型的参数，以此降低模型的复杂度。正则化，又称权重衰减，是防止过拟合的重要方法之一。LASSO降维本质就是约束或限制要优化的参数，实际上就是在不可约平面代数曲线的奇点处，把不同切线的曲线分支区分开，进而消除奇异性。采用正则化的方法可以减弱不重要的变量，提取重要的特征变量，从而达到减小特征变量数量级的作用，降低模型的复杂度，保留解释性较强的特征子集。

（二）LASSO方法原理

LASSO方法是1996年由Robert Tibshirani首次提出的回归方法。该方法是估计稀疏系数的线性模型，通过构造一个惩罚函数得到一个较为精炼的模型，压缩一些回归系数，即强制回归系数绝对值之和小于某个固定值，同时设定一些回归系数为零，从而保留子集收缩的优点，可以解决多重共线性的问题。

考虑多元线性模型（式5-1）：

（式5-1）

其中， X _{n × p} 为自变量矩阵， β _{p ×1} 为回归系数向量， ε _{n ×1} 为误差向量， Y _{n ×1} 为响应向量， β ₀ 为截距。

LASSO方法选用 L 1正则化，其估计定义为（式5-2）：

（式5-2）

其中， ， ， λ ≥0是一个惩罚项。

如果选择 L 2正则化，就是岭回归（ridge regression）分析。 L 1、 L 2正则化都能防止过拟合，提升模型的泛化能力， L 1正则化得到的解更加稀疏。LASSO回归估计解决了加惩罚项 λ || β || ₁ 的最小二乘的最小化。当 λ =0时，即为普通线性回归的最小二乘算法。LASSO的求解方法包括坐标下降法、近端梯度下降法和最小角回归等方法。

LASSO是一种变量筛选方法，这种方法适用于自变量个数远远大于样本量个数的数据。例如基因位点数据，每个人测序后位点可能成千上万个，但测序的人数可能只有几十例，因此传统回归的前进法、后退法、逐步回归法和Wald法等均不再适用。此时可以通过LASSO回归筛选变量，再用LASSO筛选后的变量构建模型。LASSO回归的特点是在拟合广义线性模型的同时进行变量筛选和复杂度调整。因此，无论因变量是连续变量还是分类变量，都可以使用LASSO回归建立模型并进行预测。LASSO回归的变量筛选是有选择地把变量放入模型从而得到更好的性能参数。复杂度调整是指通过一系列参数控制模型的复杂度，从而避免过度拟合。复杂度调整的程度由参数 λ 控制， λ 越大对变量较多的线性模型的惩罚力度就越大，从而获得一个变量较少的模型。LASSO回归模型最终能确定具有最小规模但其解释性达到最大化的预测特征子集。

目前较为好用的拟合广义线性模型的R软件包是glmnet软件包，由LASSO回归的发明人，斯坦福统计学家Trevor Hastie领衔开发。其特点是对一系列不同 λ 值进行拟合，每次拟合都用到上一个 λ 值拟合的结果，可同时进行系数估计和特征选择，从而大大提高了运算效率。此外它还包括了并行计算的功能，可以调动一台计算机的多个核或多个计算机的运算网络，进一步缩短运算时间。glmnet可以拟合线性回归模型、二分类Logistic回归模型、多分类Logistic回归模型、Poisson回归模型和Cox比例风险回归模型等。

以下将通过案例讲解如何使用R软件的glmnet包拟合基于LASSO方法的回归模型。

二、实例应用

例5-1　基于LASSO方法的线性回归模型

某种疾病患者有12名，每名患者测量了8 000个蛋白的表达水平。通过蛋白的差异表达分析，从中提取了65个差异蛋白，如表5-1所示。其中 x 1～ x 65为65个差异蛋白的表达值， y 为一个临床指标（响应变量），数据命名为linearlasso.csv，保存于d盘中。请用LASSO方法从65个差异蛋白中提取与 y 相关的特征蛋白。

在R窗口中输入语句：

注：family规定了模型的类型，“gaussian”适用于响应变量为一维连续变量；“mgaussian”适用于响应变量为多维连续变量；“poisson”适用于响应变量为非负计数型变量；“binomial”适用于响应变量为二元离散变量；“multinomial”适用于响应变量为多元离散变量。输出结果（部分）如图5-1所示：

在图5-1的输出结果中，每一行代表了一个模型。 Df 表示自由度，代表了非零线性模型拟合系数的个数。%Dev表示由模型解释的残差比例，对线性模型来说相当于决定系数 R ² 值，该值范围为0～1，越接近1说明模型的表现越好。Lambda表示每个模型对应的 λ 值。从结果可以看出，随着 λ 值逐渐变小，越来越多的自变量被纳入模型中，%Dev也越来越大。其中第64行中（图5-2），模型包含了8个自变量，%Dev也达到了93.86%，说明获得的包含8个自变量（特征蛋白）的模型已经能很好地描述这组数据。

下面输出回归系数，在R窗口中输入语句：

coefficients＜ -coef（fit，s= fit$lambda[64]）

coefficients

注：这里s表示取第64行的 λ 值。输出结果中提取的8个特征蛋白的回归系数见表5-2。

还可以通过绘图观察模型的回归系数是如何变化的，在R窗口中输入语句：

plot（fit，xvar="lambda"，label=TRUE）

输出的图形如图5-3所示：

图5-3中每一条曲线代表每一个指标变量回归系数的变化轨迹，纵坐标表示回归系数值，底部横坐标表示log（ λ ），顶部横坐标表示此时模型中非零系数的个数。例如，指标变量 x 23， λ 值在0.1～0.4时有非零的系数，随着 λ 值的变小系数不断增大。指定好 λ 值，就可以对新数据进行预测。此处仍然指定第64行的 λ 值，在R窗口中输入语句：

nx=matrix（rnorm（325，mean=0.75，sd=0.7），5，65）#通过产生随机数构建5个新样本的自变量矩阵

predict（fit，newx=nx，s= fit$lambda[64]） #预测响应变量 y 值

此时，就预测出了这5个样本的 y 值，输出结果如图5-4所示。

当响应变量是二分类变量时，经常使用二分类Logistic回归。此时只需要将构建模型语句中的family="gaussian"修改为family="binomial"就可以实现基于LASSO方法的Logistic回归。

例5-2　基于LASSO方法的Logistic回归模型

某医生收集了某地区成年人样本10 000份，记录了age（年龄）、triglyceride（甘油三酯）、cholesterol（胆固醇）和LDL（低密度脂蛋白），并将其作为指标变量，outcome（高脂血症）作为二分类的响应变量。为了便于说明，此处随机抽取157份样本，将数据命名为

logisticlasso.csv，如表5-3所示，并保存于d盘。请用LASSO方法提取与高脂血症相关的指标。

在R窗口中输入语句：

在上述基于LASSO方法的线性回归模型中，通过先拟合模型，而后选取最优的 λ 值。但在这种方法下，所有数据都被做了一次拟合，有可能造成过拟合的现象。因此，此处采用交叉验证（cross validation，CV）的方法拟合和选取模型，同时这种方法对模型的性能会有更准确的估计。下面应用交叉验证法构建模型，在R窗口中输入语句：

cv.fit＜-cv.glmnet（x，y，family="binomial"，type.measure="class"） #应用交叉验证法构建模型

plot（cv. fit） #绘制交叉验证曲线图

注：type.measure是用来指定交叉验证选取模型时希望最小化的目标参量；“class”表示使用模型错分误差（misclassification error，ME）；“deviance”为默认值，表示似然函数值自然对数的-2倍（-2log-likelihood），常用来反映模型的拟合程度，其值越小，表示模型拟合程度越好；“mse”表示使用拟合响应变量与实际响应变量的均方误差（mean squared error，MSE）；“mae”表示使用拟合响应变量与实际响应变量的平均绝对误差（mean absolute error，MAE）；“auc”表示使用综合考量模型性能的受试者操作特征（receiver operating characteristic，ROC）曲线下面积（area under curve，AUC）。

绘制的交叉验证曲线图如图5-5所示。

对每一个 λ 值，图5-5中的圆点及上下误差线表示交叉验证获得的目标参量错分误差的平均值及95%置信区间。图中的两条虚线表示两个特殊的 λ 值，分别为lambda.min和lambda.1se。lambda.min是平均最小错分误差对应的 λ 值，而lambda.1se是平均最小错分误差1倍 SE 对应的 λ 值，其对应的模型更为简洁。 λ 值达到一定大小后，继续增加模型的指标变量个数，也就是缩小 λ 值，并不能很显著地提高模型性能，因此lambda.1se给出的是一个具备优良性能但指标变量个数最少的模型。也就是说，lambda.1se构建的模型最简洁，使用的指标变量数量最少，而lambda.min则准确率更高，使用的指标变量数量更多。

指定 λ 值取lambda.min，获得回归系数，在R窗口中输入语句：

coefficients＜-coef（cv. fit，s=cv. fit$lambda.min）

coefficients

输出结果如图5-6所示。

图5-6的输出结果中保留了不为0的回归系数，其中只有age（年龄）是筛选出的与高脂血症相关的影响因素，则预测概率模型为 。假设有一个新的样本，该样本的各个指标分别为：age=50，triglyceride=0.97，cholesterol=8.45，LDL=2.16，则该样本发生高脂血症的概率为 。下面应用R软件对新的5个样本进行预测，在R窗口中输入语句：

x1 ＜ -c（50，0.97，8.45，2.16）

x2 ＜ -c（30，1.12，5.51，1.87）

x3 ＜ -c（70，2.56，8.78，4.20）

x4 ＜ -c（38，1.06，4.14，3.02）

x5 ＜ -c（65，2.33，5.20，1.98）

newx=matrix（cbind（x1，x2，x3，x4，x5），5，4，byrow=TRUE） #构建5个新样本的指标变量矩阵

predict（cv.fit，newx，type="response"，s=cv.fit$lambda.min）#对样本结局事件进行预测

注：type=link给出的是线性预测值，即进行Logit变换之前的值；type=response给出的是概率预测值，即进行Logit变换之后的值；type=class给出0和1预测值。输出结果展示了5个样本的概率预测值，如图5-7所示。

当数据含有时间变量和结局事件时，经常采用Cox回归模型。此时输入family="cox"就可以实现基于LASSO方法的Cox回归。

例5-3　基于LASSO方法的Cox回归模型

假设下载了GEO数据库中某肿瘤样本的基因表达数据，数据中含有52个肿瘤样本和20 152个基因表达值。每个肿瘤样本都记录了随访时间（time）和结局状态（status）。为了便于说明，从数据中提取了20个差异表达基因（表5-4）。将数据命名为coxlasso.csv，保存于d盘。请用LASSO方法提取与生存相关的差异表达基因。

在R窗口中输入语句：

read.table（"d：\\coxlasso.csv"，header=TRUE，sep="，"）-＞a #读入数据

x＜-a[，3：22] #提取指标变量

x＜-as.matrix（x） #转化为指标变量矩阵

y=cbind（time=a[，1]，status=a[，2]）#提取时间变量和结局变量

cvfit=cv.glmnet（x，y，family="cox"） #应用交叉验证法构建模型

plot（cvfit） #绘制交叉验证曲线图

绘制的交叉验证曲线图如图5-8所示。

同上述说明，图中的两条黑色虚线分别表示了lambda.min和lambda.1se。此处仍然指定 λ 值为lambda.min，获得回归系数，在R窗口中输入语句：

coefficients＜ -coef（cvfit，s=cv.fit$lambda.min）

coefficients

输出结果如图5-9所示。

从图5-9的结果可以看出，差异表达基因 DNAL4 、 ECM1 、 EGFR 、 FAM3B 和 GFPT2 是与生存相关的基因。下面应用构建的模型，通过原始基因表达矩阵x预测患者的结局并进行Harrel’s C 统计量的计算，在R窗口中输入语句：

pred=predict（cvfit，newx=x） #通过基因数据矩阵x预测患者结局

Cindex（pred，y）#计算Harrel’s C统计量

输出结果如图5-10所示。

从图5-10的输出结果可以看出，模型将患者结局均预测为0，Harrel’s C 统计量的值也仅为0.5，该结果说明模型的预测效果一般。在这种情况下，如果样本量较少但指标变量成千上万，可以先通过LASSO回归筛选出变量，再将筛选出的变量重新进行Cox回归分析从而获得相应的分析结果。

适应性LASSO（adaptive LASSO，ALASSO）是在LASSO回归的基础上改进而成。对于LASSO的有偏性，适应性LASSO对系数采用不同权重进行二次惩罚，其作用是对越重要的变量惩罚越小，这样就可以使重要的变量更容易被挑选出来，而不重要的变量更容易被剔除，很好地弥补了LASSO的缺陷。

例5-4　适应性LASSO回归

此处仍然采用基于LASSO方法的线性回归模型的案例数据linearassoc.csv，采用适应性LASSO回归进行分析并将结果与LASSO回归进行比较。

适应性LASSO回归采用R软件的msgps包实现。在R窗口中输入语句：

主要输出结果见表5-5。

表5-5中，复杂度参数（complexity parameter，CP）准则筛选出了3个指标变量；校正的赤池信息准则（Akaike information criterion corrected，AICC）筛选出了13个指标变量；广义交叉验证（generalized cross validation，GCV）准则筛选出了16个指标变量；贝叶斯信息准则（Bayesian information criterion，BIC）筛选出了2个指标变量。由于BIC准则的选取结果优于其他三种准则，因此按照BIC准则，适应性LASSO方法比LASSO方法筛选出的指标变量个数要少，也更加精确（等于0或接近0的指标变量系数更多），而且筛选出的指标变量 x 23和 x 56在LASSO回归中的系数也较高，说明这两个指标变量可能是影响响应变量的潜在重要因素。适应性LASSO方法和LASSO方法都可以相对准确地筛选出变量，并且适应性LASSO方法比LASSO方法选出的变量数相对更少，结果也更为准确。

在数据分析时，考虑到协变量之间存在一些组结构，在进行变量筛选时应该同时选入模型。在分类数据分析中，对组结构的处理方式通常是将其哑变量化。例如肺肿瘤的切除术式有三种：楔形、肺段和肺叶。假设以楔形术式为参照，编码为（0 0），则肺段术式编码为（1 0），肺叶术式编码为（0 1）。在变量筛选时，将其归为一个组。Yuan等在2006年将LASSO方法推广到组（group）上面，诞生了组LASSO（group LASSO）。该方法针对含有组的回归方程，加上组惩罚项。在运算前，需要对定义的组变量进行正交标准化。由于组LASSO的求解过程较慢，因此需要利用统计优化算法对求解过程进行优化，进而加快组LASSO的求解。组LASSO是对LASSO方法的推广，对组特征的选择进行研究，改进了LASSO在组结构数据上的缺陷。

例5-5　组LASSO回归

采用R软件的grpreg软件包进行组LASSO回归。该软件包中自带数据集Birthwt，该数据集为1986年马萨诸塞州斯普林菲尔德贝州医疗中心收集的婴儿出生体重数据，包括189个样本以及16个指标变量和1个结局变量。结局变量分为连续型（婴儿出生体重数值）和二分类（婴儿出生体重＜2.5kg和≥2.5kg）。请采用该数据进行组LASSO回归。

在R窗口中输入语句：

注：x表示提取数据集中的指标变量矩阵，y表示提取数据集中结局变量的二分类数据类型。group是一个向量，用来描述指标变量的分组。其中，age1、age2和age3分别表示婴儿母亲年龄的1次、2次和3次正交多项式；lwt1、lwt2和lwt3分别表示婴儿母亲末次月经时体重的1次、2次和3次正交多项式；white和black表示婴儿母亲的种族（以其他种族作为参照）；smoke表示婴儿母亲孕期吸烟状况；ptl1和ptl2m表示婴儿母亲有1个或2个及以上婴儿早产经历（以没有早产经历作为参照）；ht表示婴儿母亲的高血压病史；ui表示婴儿母亲出现子宫易激；ftv1、ftv2和ftv3m分别表示婴儿母亲在妊娠期前3个月就诊一次、两次或三次及以上（以未就诊作为参照）（表5-6）。在R窗口中输入语句：

cvfit＜-cv.grpreg（x，y，group） #应用交叉验证法计算模型交叉验证误差

summary（cvfit）

plot（cvfit） #绘制交叉验证曲线图

输出结果如图5-11所示。

图5-11输出结果显示，当 λ =0.0261时，模型的交叉验证误差最小（0.20），此时模型筛选出16个非0指标和8个非0组变量指标。绘制的交叉验证曲线图如图5-12所示。

下面通过组LASSO回归获取指标变量的参数估计系数。在R窗口中输入语句：

fit＜-grpreg（x，y，group，penalty="grLasso"，family="binomial"） #组LASSO回归

fit

按照 λ =0.0261，模型最终纳入的指标及参数估计系数见表5-6。