Stata 16 新功能之Lasso系列（一）：Lasso Basics

Stata 16 系列介绍之开篇

谨以此文献给山东大学

“我的山大，我的家”

Stata 16于2019年6月26日隆重推出，开启了在大数据时代的华丽转身。本公众号将陆续介绍Stata 16的一系列强大功能。

其中，Stata 16最引人注目的新功能为Lasso系列的命令。为此，Stata 16专门提供了一个全新的Lasso模块，相应的Lasso Manual达到355页。新增的Lasso模块包括与Lasso相关的丰富命令与最新方法：lasso，square-root lasso，adaptive lasso，ridge，elastic net；非线性lasso（logit lasso, probit lasso, Poisson lasso）；以及lasso模型的统计推断方法（double selection，partialing-out, cross-fit partialing out）。本公众号计划分三次推送介绍。

惩罚回归的基本原理

本质上，lasso 系列的方法均为针对“高维数据”（high-dimensional data）的“惩罚回归”（penalized regression），只是具体的“惩罚函数”（penalty function）有所不同。所谓高维数据，就是解释变量个数 p 超过样本容量 n 的数据。此时，由于存在严格多重共线性，故 OLS 不存在唯一解，虽然样本内可以完美地拟合数据，但模型的样本外预测能力可能很差，存在“过拟合”（overfit）。

为了解决过拟合问题，并降低估计量的方差（若存在严格多重共线性，则 OLS 方差无穷大），常使用“惩罚回归”（penalized regression），即在 OLS 的目标函数（残差平方和）之外，再加上一个惩罚项，惩罚太大的回归系数。

岭回归（Ridge Regression）

对于岭回归（ridge regression），使用 L2范数（即欧氏距离）作为惩罚函数： 

其中，为“调节参数”（tuning parameter），控制惩罚力度，通常使用“交叉验证”（cross-validation）来确定（即选择使得模型的预测误差最小，详见下文）。 为参数向量 的2-范数（L2 norm），即该向量的长度 。

求解此最小化问题，从其一阶条件，即可得到上述岭回归估计量的解析表达式。显然， 依赖于调节参数 ，可记为 ，此函数被称为“系数路径”（coefficient paths）。

为了解释岭回归，可将其目标函数等价地写为一个有约束的极值问题：

其中， 为某常数。对于此约束极值问题，可引入其拉格朗日乘子函数，并以  作为其乘子，即可得到前述的岭回归目标函数。由于约束集 为 p 维参数空间中的圆球，故可将此约束极值问题图示如下（假设 p =2）。

在上图中， 为OLS估计量，围绕  的椭圆为残差平方和的“等高线”，而灰色的圆球则为约束集（即可行的参数取值范围）。岭回归估计量即为椭圆形等高线与圆球形约束集相切的位置。

从上图还可看出，由于约束集为圆球，故等高线与约束集相切的位置一般不会碰巧在坐标轴上，故通常只是将所有的回归系数都收缩，而不会让某些回归系数严格等于0。因此，岭回归一般得不到稀疏解，不便于解释回归系数。

套索估计量（Lasso）

为此，套索估计量（Least Absolute Shrinkage and Selection Operator，简记 LASSO）应运而生。自 Tibshirani (1996) 提出 lasso 之后，很快成为大数据时代炙手可热的新宠。

事实上，套索估计量只是将岭回归的惩罚项（也称为“正则项”，regularization）作了“小小”的技术调整，即将2-范数改为1-范数：

其中， 依然为调节参数，控制惩罚的力度；而   为参数向量  的1-范数（L1 norm），即 ，为回归系数的绝对值之和。类似于岭回归，上述的 lasso 最小化问题可等价地写为如下约束极值问题：

其中， 为某常数。显然，此约束极值问题的约束集 不再是圆球，而一般是菱形或高维的菱状体。仍以 p =2为例，可将 lasso 的约束极值问题图示如下。

从上图可直观看出，由于 lasso 的约束集为菱形（而菱形的顶点恰好在坐标轴上），故椭圆等高线较容易与此约束集相交于坐标轴的位置，导致 lasso 估计量的某些回归系数严格等于0，从而得到一个稀疏模型（sparse model）。Lasso 的这种独特性质，使得它具备了“变量筛选”（variable selection）的功能，故也称为“筛选算子”（selection operator）。

弹性网估计量（Elastic Net）

Zou and Hastie (2005)将 lasso 与 ridge 相结合，同时使用 L1 与 L2 惩罚函数，得到“弹性网估计量”（Elastic Net）：

其中， 与  都是调节参数。显然，如果 ，则弹性网就是 lasso；而如果 ，则弹性网就是 ridge。

平方根套索（Square-root Lasso）

平方根套索（square-root lasso）由Belloni, Chernozhukov and Wang (2011) 所提出，其目标函数为：

虽然平方根套索在理论上有其方便之处（更易从理论上估计最优调节参数，即所谓 plugin estimator），但由于在实践中一般使用交叉验证选择调节参数 ，故平方根套索较少使用。

交叉验证（Cross-validation）

对于以上的 lasso 系列惩罚估计量，在实践中都需要确定最优的调节参数，比如  或 。Stata 16 提供了三种选择最优调节参数的方式，即交叉验证（cross-validation）、adaptive lasso 与 plugin estimator。由于实践中一般使用交叉验证，故这里仅介绍交叉验证法。

Stata 16的lasso系列命令默认进行“10折交叉验证”（10-fold cross-validation），如下图所示。

首先，将样本数据随机地10等分，得到大致相等的10个子样本（即所谓的10折，每折对应于一个子样本）。然后使用其中的9折（即9个子样本）估计此模型，并将所得模型预测其余1折的被解释变量，得到相应的均方误差（Mean Squared Error，简记 MSE）。如此依次进行，可得到10个均方误差，其平均值即为整个样本数据的交叉验证MSE，也称“交叉验证误差”（CV Error）。

显然，CV Error 为调节参数 的函数，可记为 ，将此函数画图即为“交叉验证图”（CV plot）。在实践中，可在  的可能取值范围 （如果  的惩罚力度足够大，所有系数将为0，故存在最大的 ）的网格上（比如 100 等分）评估此  函数 ，找到其最低点，即最优调节参数 。

Stata 实例

下面使用 Tibshirani (1996) 所用的前列腺癌数据（prostate cancer data）进行演示。该数据集包含 97 位即将进行前列腺手术男子的数据。被解释变量为前列腺特异性抗原的对数（log of prostate specific antigen，记为lpsa），解释变量共有8个，参见下面的论文截图。虽然此例并非高维数据，但由于为 lasso 创始人 Tibshirani (1996) 所用，故经常作为案例。

首先，从 Stanford 大学网站下载此数据集，并查看数据的统计特征。

insheet using https://web.stanford.edu/~hastie/ElemStatLearn/datasets/prostate.data, clear tab

sum

其次，由于这是低维数据，不妨进行 OLS 回归，作为参照系：

reg lpsa lcavol lweight age lbph svi lcp gleason pgg45,r

从上表可知，样本内的为 0.6634，但这很可能高估了此模型在样本外的预测能力。

下面，使用命令 lasso linear 进行 lasso 回归，其基本格式与 reg 类似：

lasso linear lpsa lcavol lweight age lbph svi lcp gleason pgg45, selection(cv, alllambdas) stop(0) rseed(12345) nolog

其中，选择项“selection(cv, alllambdas)”表示使用交叉验证（默认为 10折）选择调节参数 ，而“alllambdas”与“stop(0)”表示搜索默认的全部100 个  值；否则，Stata 将在交叉验证误差不再下降之后即停止搜索。选择项“rseed(12345)”表示在将样本随机分为10 等分时，使用随机种子12345（可自行指定），便于重复结果。选择项“nolog”表示不显示交叉验证的过程。

上表显示，根据 10 折交叉验证，调节参数  的最优取值为 0.0034851，共有 8 个非零回归系数（在此例中所有变量的系数均非零）。相应的样本外（Out-of-sample R-squared）为 0.5982，比上述 OLS 的样本内0.6634 更低；而交叉验证误差（CV mean prediction error）为 0.5299。

下面，使用命令 coefpath 来画 lasso 的系数路径（coefficient paths）：

coefpath,legend(on position(12) cols(4)) 

其中，选择项“legend(on position(12) cols(4))”表示将图例（legend）放在 12 点钟的位置（即图像正上方），并以 4 列（columns）来表示。 

命令 coefpath 默认将 L1 范数（L1 norm）作为横轴来画系数路径。显然，当 L1 范数为 0 时，所有系数均为 0；而当 L1 范数足够大时，则为 OLS 估计（对于低维数据而言）。

如果想以调节参数 （log scale，对数尺度）作为横轴来画系数路径，并在最优值  = 0.0034851 处画一条垂直线，可使用如下命令：

coefpath, legend(on position(12) cols(4)) xunits(lnlambda)  xline(.0034851)

从上图可知，当惩罚力度 = 0 时，即为 OLS 估计；而当  足够大时，惩罚过于严厉，以至于所有回归系数均为 0。

进一步，可以画交叉验证图，其命令为：

cvplot

上图给出了调节参数  的最优值 ，即使得函数 最小的  值。从上图可知，在最优值  附近，函数  非常平坦，这意味着在最优值附近变化 ，对于模型的预测能力影响很小。

如想显示 lasso 回归的系数，可使用命令

lassocoef, display(coef) sort(coef)

其中，选择项“display(coef)”表示显示回归系数（默认仅显示选中的变量，即系数非零的变量），而选择项“sort(coef)”表示将变量按照系数大小排列。

以上显示的回归系数事实上为“标准化”（standardized）之后的系数。对于惩罚回归而言，变量的单位或取值范围对回归结果有实质影响（因为回归系数为惩罚对象），故一般在进行 lasso 系列的回归之前，须先将所有解释变量进行标准化，即减去其样本均值，再除以其样本标准差，使得每个变量的均值都为0，而标准差均为 1 。

如果想显示标准化之前的回归系数，可使用命令

lassocoef, display(coef,penalized) sort(coef,penalized)

注意到以上两表中，均未显示回归系数的标准误或 p 值。这是因为 lasso 系列的估计量，仍没有公认的标准误。关于 lasso 系列命令的统计推断，我们将在本系列推文中专门介绍。

下面，我们演示平方根套索的操作，其命令也类似：

sqrtlasso lpsa lcavol lweight age lbph svi lcp gleason pgg45, selection(cv, alllambdas) stop(0) rseed(12345) nolog

上表显示，最优  值为 0.3545764。注意平方根 lasso 的最优  值与 lasso 相应的最优  值并没有可比性（因为目标函数不同）。

接下来看平方根 lasso 的系数路径：

coefpath,legend(on position(12) cols(4)) xunits(lnlambda) xline(.3545764)

从上图看到，平方根 lasso 的系数路径有时并不光滑，这或许是由于其目标函数中的平方根所致。

下面进行弹性网估计，其命令依然类似：

elasticnet linear lpsa lcavol lweight age lbph svi lcp gleason pgg45,alpha(0(0.1)1) selection(cv, alllambdas) stop(0) rseed(12345) nolog

其中，选择项“alpha(0(0.1)1)”表示对于调节参数  ，从0（对应于 ridge）到 1（对应于 lasso），以 0.1 为间隔，进行网格搜索。

上表显示，交叉验证选择了  =0 而  = 0.0356706（上表中打星号 * 处），故对于此数据集而言，岭回归的预测效果最佳。

下期公众号将继续介绍非线性 lasso（logit lasso, probit lasso, Poisson lasso 等），敬请期待。

备注：本公众号不提供 Stata 16。如需正版 Stata 16，请咨询 Stata 软件官方授权经销商及合作伙伴：北京友万信息科技有限公司（www.uone-tech.cn），希望能为 Stata 中国用户提供更多服务与支持，联系人：徐老师，Tel/Wechat：18610597626。

参考文献

陈强，《计量经济学及Stata应用》，高等教育出版社，2015年（配套教学视频，可在网易云课堂学习，详见https://study.163.com/course/introduction/1006076251.htm）

陈强，《高级计量经济学及Stata应用》，第2版，高等教育出版社，2014年（配套高级计量六天现场班，北京，2019年10月1-6日，详见https://bbs.pinggu.org/thread-3156565-1-1.html）

陈强，《机器学习及R应用》，高等教育出版社，2020年（即将出版）

震撼来袭  >>  机器学习及Stata、R三天现场班

上海，2019年8月17日-19日

主办：第三届Stata中国用户大会、友万科技

主讲：陈强教授 （山东大学）

授课方式：思想原理 + 数学精髓 + Stata、R案例

陈强老师将首次推出全新的“机器学习及Stata、R应用”三天现场培训班。结合Stata与R的实操案例，深入浅出地介绍最为流行的机器学习方法，包括KNN，判别分析、朴素贝叶斯、决策树、随机森林、提升法、支持向量机、神经网络等。

跟着陈老师，三天入门机器学习，赶上时代步伐！

更多详情，请点击页底“阅读原文”，或识别下图中二维码：

(c) 2019， 陈强，山东大学经济学院

www.econometrics-stata.com

转载请注明作者与出处

Our mission is to make econometrics easy, and facilitate convincing empirical works.