高维回归Lasso之Stata操作案例

当计量经济学遭遇机器学习：揭开机器学习的神秘面纱（一）

当计量经济学遭遇机器学习（二）

当计量经济学遭遇机器学习（三）：高维回归之岭回归

当计量经济学遭遇机器学习（四）：高维回归之LASSO

套索估计量（Lasso）之简要回顾

所谓 “高维回归”，指的是回归的样本容量 n 超过变量个数 p 的情形；比如，样本为 100 个病人，而每个病人有 2 万个基因。此时由于存在严格多重共线性，无法使用 OLS，而套索估计量（Least Absolute Shrinkage and Selection Operator，简记 LASSO）则是非常流行的高维回归方法。

事实上，经济学家经常将 Lasso 应用于传统的低维数据，因为 Lasso 具有 “变量选择”（variable selection）的功能。

考虑以下线性回归模型：

Lasso 估计量（Tibshirani,1996）求解以下最小化问题：

 

其中， 为微调参数（tuning parameter），控制惩罚的力度（惩罚过大的参数）；而  为参数向量  的 1-范数（L1 norm），即

 

为各回归系数的绝对值之和。不难看出，由于惩罚项的存在，故 Lasso 为收缩估计量，即相较 OLS 估计量更为向原点收缩（如果没有此惩罚项，则就是 OLS）。

上述 Lasso 最小化问题可等价地写为如下约束极值问题：

 

其中， 为某常数（可视为微调参数）。不难看出，此约束极值问题的约束集  为菱形或高维的菱状体。以 p=2 为例，可将 Lasso 的约束极值问题图示如下。

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

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

微调参数（tuning parameter）的选择

显然，使用不同的微调参数 ，可得到不同的 Lasso 估计值。因此，如果变化 ，则可得到整个解的路径（solution path）。那么究竟应该怎样选择微调参数  呢？

一般的方法是，选择  ，使得此模型的预测能力最强（模型的预测能力显然是  的函数）。具体来说，常进行 k 折交叉验证（k-fold cross validation），其中一般取 k = 5 或 10。以 k =10为例，参见下图。

首先，将样本数据随机分为10个等分。将第1个子样本作为 “验证集”（validation set）而保留不用，而使用其余9个子样本作为 “训练集”（training set）来估计此模型，再以此预测第1个子样本，并计算第1个子样本的 “均方预测误差”（Mean Squared Prediction Error）。

其次，将第2个子样本作为验证集，而使用其余9个子样本作为训练集来预测第2个子样本，并计算第2个子样本的 MSPE。以此类推，将所有子样本的 MSPE 加总，即可得整个样本的 MSPE。

最后，选择微调参数 ，使得整个样本的 MSPE 最小，故具有最佳的预测能力。

Lasso 的 Stata 操作案例

在 Stata 中进行 Lasso 估计，可使用非官方命令 lassopack，其安装方法为

ssc install lassopack

Lassopack 包含三个与 Lasso 相关的子命令。其中，lasso2 可进行 Lasso 估计，以及 Lasso 的多个变种，包括 square-root lasso, elastic net, ridge regression, adaptive lasso 等。另一子命令 cvlasso 可进行 k 折交叉验证（k-fold cross validation）。

下面以 Tibshirani (1996) 所用的前列腺癌数据集作为案例。此数据集其实为低维数据（n = 97, p = 8），但它是创始人 Tibshirani 所用的案例，故常作为 Lasso 的演示数据。

该数据集包含了 97 位男病人的数据，被解释变量为“前列腺特异性抗原的对数”（log of prostate specific antigen，简记 lpsa），而解释变量则包括年龄（age），癌体积对数（log of cancer volume，简记 lcavol）等，参见 Tibshirani (1996) 的论文截图。

首先，载入数据，并察看数据特征。

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

. sum

其次，进行 Lasso 估计，并画出整个解的路径（solution path）。

. lasso2 lpsa lcavol lweight age lbph svi lcp gleason pgg45, plotpath(lambda)

其中，选择项 “plotpath(lambda)”，表示根据微调参数  的大小，画出整个解的路径。

上表显示，当微调参数  由大到小变化时，Lasso 估计量也筛选出越来越多的变量进入模型。比如，当  >163.62492 时，由于惩罚力度过大，所有变量的系数均为 0。当  =163.62492 时，常数项首先进入模型；而当  降为149.08894 时，变量 lcavol 进入模型，以此类推。

下图为整个解的路径（作为  的函数），画出了不同变量回归系数的变化过程。其中，当  = 0 时（下图最左边），不存在惩罚项，故此时 Lasso 等价于 OLS。而当  很大时（下图最右边），由于惩罚力度过大，所有变量系数均归于 0。

那么，究竟应如何选择微调参数  的具体取值呢？下面进行 k 折交叉验证 (k-fold cross validation)

cvlasso lpsa lcavol lweight age lbph svi lcp gleason pgg45, lopt seed(123)

其中，选择项 “lopt” （为 lambda optimal 之缩写）表示选择使 MSPE 最小的  ；选择项 “seed(123)” 表示将随机数种子设为 123 （可自行设定），以便结果具有可重复性；默认 k =10（即 10 折交叉验证）。

在上表中，打星号处的  = 14.566138，即为可使 MSPE 最小化的  值。与此最优  对应的 Lasso 估计结果为：

上表右边第1列即为 Lasso 所估计的变量系数。其中，除常数项外，只有 5 个变量的系数为非零，而其余 3 个变量（未出现在表中）的系数则为 0。

考虑到作为收缩估计量的 Lasso 存在偏差（bias），上表右边第2列汇报了 “Post Lasso” 估计量的结果，即仅使用 Lasso 进行变量筛选，然后扔掉 Lasso 的回归系数，再对筛选出来的变量进行 OLS 回归。

展望

无论在统计学、计量经济学，还是经济学领域的应用，以 Lasso 为代表的高维回归方法均方兴未艾，影响深远。更详细的介绍将在 “高级计量经济学与Stata现场班”（北京，国庆节）进一步展开，包括岭回归（ridge regression）、弹性网（elastic net）、适应性套索（adaptive Lasso）、post Lasso、post double Lasso、IV Lasso 等。