|
Lasso其实是“least absolute shrinkage and selection operator”的缩写，即最小绝对值压缩和选择算子，但现在往往直接翻译成套索。从本质上来讲，lasso回归并不和线性回归、logistic回归、Poisson回归、Cox回归等一样是某一类回归，而是和逐步回归一样，是一种算法、一种策略，诸多回归都可以使用这种lasso策略，用于结局预测、模型选择和因果推断。 在传统的回归模型中，模型系数估计常用最小二乘法(Ordinary Least Square, OLS)，OLS是通过最小化残差平方和（residual sum of squares，RSS，结局变量的预测值与实际观测值之间差值的平方和）来估计回归模型中的参数，即估计参数使得模型 (Y-βX)2最小，可以针对当前样本提供最准确的线性无偏估计。OLS的问题在于容易导致过拟合（Overfit），缺乏概化能力（Generalizability），即基于当前样本获得的回归模型在预测来自同一总体的其他观测数据时表现不佳。模型的预测误差可以被分解为偏差（Bias）和方差（Variance）两部分， 其中偏差指预测值和真实值之间的差异, 方差指预测值的离散情况。所谓模型偏差小，指的是模型拟合效果好，模型曲线尽可能切合样本的所有数据点；而方差小则反应的是对于同一个总体的数据，每个抽样的样本拟合的模型变异不大，而不是换一个样本模型相差十万八千里，在曲线上表现为比较平滑。一般来说，偏差很小方差会很大，所以在实际中需要能很好地处理这种偏差−方差权衡(Bias-Variance Tradeoff)问题。OLS是通过控制估计偏差来降低模型的预测误差, 更关注对当前数据集的精确估计，但却导致样本间方差增大, 使得当前的参数估计结果可能仅适用于当前数据集。另外，OLS非常容易受到预测变量间的多重共线（Multicollinearity）的影响。当遇到包含了许多变量的大数据集时，多重共线的问题其实是很难避免的。最后，如果想得到相对稳定的模型，OLS对样本量是有要求的，这里大家可能都知道10倍EPV原则（Events Per Variable，每个变量对应10个阳性事件。当结果阳性事件更多时则是阴性事件满足10EPV），我们常常会遇到样本量不足的问题，甚至变量数大于样本量。Lasso回归是引入正则化方法通过在模型估计中增加惩罚项的方式将过小的回归系数压缩到0, 以一定的估计偏差为代价获得更高的模型预测准确度和模型概化能力。简单说，就是利用压缩（shrinkage）思想，压缩系数估计值，减小方差。由于lasso可以将预测变量的估计系数压缩为0, 因此可以在压缩系数的同时起到变量筛选的作用。不明白为什么把Regularization翻译成“正则化”这样一个稀奇古怪不容易理解的名词，其实本意就是规范化，既然要规范，就需要限制一些东西以不让你太膨胀，限制什么呢？就是在标准的OLS损失函数上加上一个惩罚项。在OLS中模型系数估计值使得(Y-βX)2值最小，(Y-βX)2就是OLS的损失函数。而在lasso中，模型系数估计值使得(Y-βX)2+λΣ|β|值最小化，这里面Σ|β|就是惩罚函数，也称L1范数（L1-Norm），是所有回归系数绝对值之和。λ(≥0)是调整参数(Tuning Parameter), 用于控制模型系数的压缩程度, λ值越大则惩罚力度越强，不难理解当λ=0时, lasso即OLS。Lasso很好地解决了模型过拟合、变量共线、样本量不足带来的问题，在满足稀疏性的条件下，甚至可以允许变量个数大于样本量。即使对于相关性不强或不相关的小数据集，利用Lasso和弹性网络也可以获得更好的预测模型，或选择用于模型简化或未来研究的变量。 Lasso采用的是L1范数正则化方法，即惩罚项是所有回归系数绝对值之和。另外一种采用L2正则化方法的压缩估计方法是领回归（Ridge Regression），L2正则化对应的惩罚函数是Σ|β|2，即所有系数的平方和，与L1正则化唯一不同之处就是将惩罚函数中的绝对值变成了平方，但领回归无法解决变量选择的问题。领回归和lasso回归可以看做是弹性网络（Elastic Net）的特例，弹性网络可以保持两者的一些特性，其有个调整参数α，α=0就是领回归，α=1就是lasso回归。Lasso回归和领回归最关键的步骤在于寻找一个恰当的λ值。参数λ值决定了回归系数被压缩的程度, λ值越大模型系数越小，不同的λ值可能产生不同的结果。过大的λ值虽然可以降低方差，但却增大了偏差，我们需要的是偏差−方差均衡。怎么选择λ值呢？一般有两种方法，一种是采用来获得不同λ值对应的模型的残差均方（MSE），再以MSE大小决定参数λ的取值。理论上应该选择MSE最小时的λ值，但是过小的λ值意味着回归系数压缩幅度较小，不一定能完全解决模型过拟合、共线性的问题，因此实际应用中常常选择大于最小MSE一个标准误时对应的λ值。有时候可能是其他指标，如-对数似然值等。第二种方法选择λ值的方法是信息准则(Information Criteria)，即先获取不同λ值下lasso回归模型的信息准则的值(如, Akaike Information Criterion, AIC、Bayesian Information Criterion, BIC)，然后选取信息准则值最小时对应的λ值。心理科学进展,2020,28(10): 1777-1788.除了λ值的计算，Lasso回归模型另外一个非常重要的参数是P值。Lasso可以选择变量并估计其系数，但却没有提供用于这些参数的统计推断所需的标准误，因此计算P值需要专门的方法。你可能觉得这并不是一个问题。因为我们可以先用套索回归来筛选变量，然后再利用筛选的变量拟合回归模型（线性、logistic、probit、poisson等）不就可以了，有大量的高分文献都是这么做的。但实际上这种做法在理论上存在着不合适的地方：①我们建立的一个统计模型，系数标准误代表的是该系数在重复抽样中的分布，尽管我们仅有一个总体的样本，但95%CI代表的是从总体的抽样中，有95%的重复抽样包含了系数的真实值。但当我们采用lasso进行变量选择时（其他变量选择也是如此），每个样本选择的变量都可能不一样，一般的回归模型系数的标准误不能对此作出解释；②lasso回归倾向于省略具有小系数的变量，而不是忽略不在真实模型中的变量，这可能会使模型发生偏倚。③非常重要的一点就是lasso可能会选择与真实模型中变量相关的一些变量，而不一定会选择在真实模型中的变量，也就是说lasso选择的变量与真实模型中的变量可能是不同的。如果lasso选中了变量X，则意味着X可能是真实模型中的一个变量，或者X与真实模型中的某个/些变量相关；如果lasso省略了变量Z，则意味着Z不属于真实模型，或者Z虽属于模型，但是与已经被lasso选择的某个/些变量相关。基于以上这些原因，我们需要特殊的lasso方法来进行统计推断。在stata中提供的用于因果推断的算法有双重选择（double selection）、偏出（partialing out）、交叉验证偏出（cross-fifit partialing out），而用于模型选择和预测的办法有lasso、弹性网络（elasticnet）和平方根lasso。通过“完全模型”与嵌套其内的限制模型之间的偏差进行卡方检验也是获得嵌套模型间差异的显著性的一种方法。在lasso回归的基础上，研究者采用不同形式的惩罚函数，建立和发展出了多种正则化模型, 例如Relaxed Lasso、Adaptive Lasso、Bayesian Lasso等。
计算机学报,2015,38(7):1307-1325. 当观测变量数p远大于观测样本量N时，Relaxed Lasso先通过普通的Lasso回归筛选出合适的预测变量，然后再对筛选出的变量进行系数估计，具有更快的收敛速度。但需要注意Relaxed Lasso中的参数估计并非OLS，而是调整参数也参与其中的一种算法。Adaptive Lasso不同于普通的lasso回归的地方在于惩罚项中每个系数都被赋予了一个权重，使得不同变量的回归系数受到的惩罚程度不同（普通Lasso方法对所有变量施加相同程度的惩罚），从而重要的变量更易进入模型, 而不重要的变量更易被剔除。Adaptive Lasso也更适用于变量数和样本量的比值非常大的情况。示例来自R语言VIM包的diabetes数据集（Indian Prime Diabetes Data）。该数据集最初来自美国国家糖尿病、消化和肾脏疾病研究所，用来根据数据集中的某些诊断测量值来诊断性的预测患者是否患有糖尿病，所有患者都是至少21岁的皮马印度血统的女性。指标释义如下：
我们先通过R命令将数据导出： data(diabetes,package="VIM") head(diabetes) write.csv(diabete,file = "D:/Temp/IPDiabetes.csv") #数据导出成csv文件，文件名称为IPDiabetes，csv文件可以直接用excel打开。 案例存在大量的缺失数据，是VIM包中用于缺失值分析的的一个数据集。本例用来演示lasso回归。考虑到演示的目的，并未对缺失值进行填补。分析软件采用JMP。分析>>拟合模型，将Outcome选入[Y]，通过构造模型效应中的[添加]将其他解释变量选入，[特质]中选入广义回归，[分布]会自动默认二项，[目标水平]为yes（Outcome中yes为糖尿病患者），点击运行。 
运行后的报表显示两大部分内容：模型启动控制面板和最大似然报表。不同的因变量分布对应的这两部分内容会有所不同，本例选择的是满足二项分布的二分类结局，可供选择的模型估计方法和验证方法比较多，最大似然报表对应的不进行惩罚的Logistic回归，常作为模型比较的基线。最大似然结果（本例为logistic回归）结果见下图右侧，相当于直接进行二分类logistic回归，如果想获得变量的OR值，可在[logistic回归]前面的红色倒三角中选中“优势比”即可。二分类logistic回归的笔者曾在《因变量二分类资料的logistic回归》、《R笔记：二分类资料的logistic回归》做过最为详尽的解读，大家可以参考。本次笔记示例的是lasso的JMP操作，因此不再对结果做过多的解释。模型启动控制面板中的[估计方法]默认lasso，本例去掉[自适应]选项；[验证方法]选择K重，[折数]修改为10，[执行]。[高级控件]中的[最初显示的解]默认的是[最佳拟合]，也可以选择经常使用的其他几个比如[绿色区域最大]等，绿色区和红色区分在选择模型时很重要，具体定义可参见后面的结果部分。[强制]选项可以选入保留在模型中不进行筛选的变量。
报表中会显示使用“K重”验证的lasso结果。主要结果包括解路径（Solution Path）、原始测量变量的参数估计值及效应检验。回归报表还可以通过相应部分标题前的红色倒三角增加“中心化和统一尺度的预测变量的参数估计值”和“活跃参数估计值”。解路径报表中有两张图，一张是解路径图，另外一张是验证图。两图横坐标都是“统一尺度的参数估计值量值（Magnitude of Scaled Parameter Estimates）”，名字啰里啰嗦的，实际就是l1范数，即所有回归系数绝对值之和。不知道为何，JMP里面没有给出λ值，但是较大l1范数对应较小的λ值，此时参数估计值接近最大似然估计的结果；而较小的l1范数对应较大的λ值，此时参数估计值会被给予较重的惩罚。解路径图的纵坐标是参数估计值，图中每条线代表每个变量的参数随横坐标的变化而出现的变化，随着统一尺度的参数估计值量值的减小（横坐标从右向左），λ值增大，惩罚力度增大，各变量的系数逐渐被压缩至0。验证图的纵坐标是选定验证方法的验证统计量的值，有-对数似然的相关指标和信息准则两大类，图中有两个区域：绿色区和黄色区，若模型的验证统计量的值落在某区域内，则该模型会落在该区域内。
我们在选择模型时，考虑到压缩的需要，未必会选择最佳模型拟合处的λ值，比如可以选择最佳拟合标准一个标准误时对应的λ值以保证适当地参数压缩幅度，此时你可以选择绿色区域的最大值。 图中的红色竖线是当前模型指示符，指示符默位置是验证准则选择的最佳拟合模型处，当然这个可以在前面提到的高级控件中的[最初选择的解]处进行修改。
结果显示，随着统一尺度的参数估计值量值的减小（横坐标从右向左），λ值增大，惩罚力度增大，各变量的系数逐渐被压缩至0，选择进入模型的变量也会逐渐变少。本例选择的是K重验证，当-对数似然相关值最小时，模型最佳，结果默认显示的也是最佳模型的参数估计值和效应检验，具体结果大家可以具体看相应的报表。如果想获得变量的OR值，可在[使用“K重”检验的“lasso”]前面的红色倒三角中选中“优势比”即可。jmp最值令人称道的就是它的交互性，比如点击解路径中的某条路径，对应的变量会在下方参数估计表中用深色背景显示，你可以拖动指示符顶部的黑色箭头来更改罚值的量值，从而获取相应的模型（最佳拟合模型处仍有一条垂直虚线），下方的参数估计值也会发生相应的变化。
在R语言中，Lasso回归常用glmnet函数来完成，对这个函数的作图可以获得MSE最小时对应的λ值和最小MSE一个标准误时对应的λ值。笔者尝试了几次，发现JMP的lasso参数估计结果和glmnet函数的不太一致，原因可能是模型的验证标准不一样；另外采用cv.glmnet函数进行K折交叉验证法选择合适的λ值时，为确保每次运行K折交叉验证的结果一样，是需要设定随机种子的。JMP中提供了多种模型验证的标准，比如-对数似然相关值和信息准则。一般情况下这些验证标准值最小时对应的模型就是最佳模型。但我们在文献里常常不一定采用这个模型，原因是这种模型对应的λ值可能会过小，这就意味着回归系数压缩幅度较小, 不一定能完全解决过拟合的问题，因此常常采用最小标准一个标准误时对应的参数 λ 值。JMP里面没有给出具体的λ值（也可能是我没找到），不过这不影响我们寻找这几个关键标准的位置，一种是直接对指示符进行拖拽，二是可以在[模型启动]控制面板中通过[高级控件]中的[最初显示的解]进行修改，比如可以选择绿色区域中的最大值对应的就是V最佳+4或者L最佳+LSE。除了菜单式的交互操作，采用JMP进行lasso回归比用R的glmnet函数还有一个很大的优势就是JMP可以直接对多分类变量进行筛选，这个在glmnet函数中是做不到的，目前glmnet函数尚不能处理分类变量，在遇到多分类变量的时候还需设置哑变量变通一下数据。下次笔记的时候我们会介绍一下。在进行lasso回归时，我们可能会更愿意采取一些改进模型，比如近似无偏的自适应lasso回归，这个在JMP中是默认的选项。最后，再看一眼数据改用自适应的lasso回归、AICc进行验证的结果：
如果你细心可能会发现采用（自适应）lasso回归和直接进行logistic回归，在模型中变量一样的情况下，参数估计值也是不一样的，原因我们在笔记的开始部分就说过了。我们常常仅把lasso回归当作变量筛选的工具，将lasso筛得的变量再进行普通的logistic回归是常见的分析套路。不是说不可以这么做，也不讲理论上存在的问题，单就我们花了这么大力气用lasso回归筛到了变量并获得了变量的系数，为什么不直接用lasso的参数结果直接进行下一步的预测或推断呢？！
|
