Lasso：拉索中如何做统计推断

Stata连享会   主页 || 视频 || 推文 || 知乎

温馨提示： 定期 清理浏览器缓存，可以获得最佳浏览体验。

New！ lianxh 命令发布了：
随时搜索连享会推文、Stata 资源，安装命令如下：
  . ssc install lianxh
使用详情参见帮助文件 (有惊喜)：
  . help lianxh

✌ 课程详情： https://gitee.com/arlionn/Course   |   lianxh.cn

✨ 课程主页： https://gitee.com/arlionn/Course

作者：杨继超 (中山大学)
E-Mail: yangjch7@mail2.sysu.edu.cn

Source： Stata Blogs: Using the lasso for inference in high-dimensional models

目录

1. 引言

高维数据模型在应用研究当中越来越普遍。上一篇推文 Stata:拉索开心读懂-Lasso 入门 中讨论的 Lasso 方法便可用于高维数据模型中估计变量的系数。但是，上一篇推文始终没有解决一个问题就是统计推断。这其实也是 Lasso 方法面临的大难题，本篇推文则是延续上一篇推文，讨论 Stata16 中关于 Lasso 方法在高维数据模型中做统计推断的相关命令。

为了讨论这一问题，本文采用的例子是来自于 Sunyer et al. (2017) 这篇文章。具体而言，这篇文章是想估计空气污染对小学学生反应时间的影响，所采用的的模型为：

变量的具体含义如下表所示：

我们想要估计 no2_class 对于 htime 的影响系数及相应的置信区间。问题是向量 $x$ 中有 252 个控制变量，但我们只有 1084 个观测值，如果直接进行回归便会出现高维数据的过拟合问题。如果直接采用 Lasso 方法进行估计，因无法得到标准误，则无法获得估计系数的置信区间，即无法进行统计推断。

实际上，除了我们关心的变量 no2_class 以外，其余变量的统计推断问题甚至系数估计都是次要的。只要我们能找到控制变量 $x$ 的一个子集 $\tilde{x}$，使得能够估计出可靠的 $\gamma$ 即可。换言之，如果我们知道 $\tilde{x}$，我们便可以使用如下的模型进行估计：

但问题是我们并不知道 $x$ 中究竟哪些变量是属于 $\tilde{x}$ 的，因此要估计 $\gamma$，首先要进行变量筛选的工作。

对于变量筛选，可参考上篇推文 Stata:拉索开心读懂-Lasso 入门 中提到的一种方法 postselection 法。具体步骤分两步，第一，采用 Lasso 方法进行估计，筛选出核心的控制变量；第二，用被解释变量对自变量和筛选出来的控制变量进行回归。

可是，直接采用 postselection 法是无法进行可靠的统计推断的。Leeb and Pötscher (2008) 指出，postselection 法产生的估计量是不具备大样本情形下正态分布性质的，并且在有限样本的情形下采用一般的大样本理论也是无法进行可靠的统计推断的。

鉴于以上分析的问题，Belloni et al. (2012)、 Belloni, Chernozhukov, and Hansen (2014)、 Belloni, Chernozhukov, and Wei (2016) 以及 Chernozhukov et al. (2018) 推导了三种能够提供对于 $\gamma$ 可靠统计推断的方法。这三种方法分别为偏回归法（partialing-out，以下简称 PO 法）、双重选择回归法（double-selection，以下简称 DS 法）和交叉拟合偏回归法（cross-fit partialing-out，以下简称 XPO 法）。下表是 Stata16 中三种方法在不同模型中的命令。

关于命令的实操参见下一小节的分析。

2. Stata 命令及操作

本小节将延续上一节所讨论的空气污染对小学生反应时间的例子，截取 Sunyer et al. (2017) 所用数据的部分变量，采用 PO、DS 和 XPO 三种方法进行估计和统计推断，并给出方法背后的一些直觉。

2.1 PO 法

Step 1: 调用数据并将连续变量和因子变量及其交互项存入暂元

webuse breathe, clear //调用数据

***********将连续控制变量存入暂元`ccontrols'中********************

local ccontrols "sev_home sev_school age precip age0  siblings_old  "

local ccontrols "`ccontrols' siblings_young no2_home green_home noise_school "

***********将因子控制变量存入暂元`fcontrols'中********************

local fcontrols "grade sex lbweight breastfeed msmoke "

local fcontrols "`fcontrols' feducation meducation overweight "

********将连续型、因子型控制变量及其交乘项存入暂元`ctrls'中*********

local ctrls "i.(`fcontrols') c.(`ccontrols') "

local ctrls "`ctrls' i.(`fcontrols')#c.(`ccontrols') "

接着对变量进行一下简单的描述性分析，代码为

. describe react no2_class `fcontrols' `ccontrols'

variable name   type    format     label      variable label
--------------------------------------------------------------------
react           double  %10.0g              * Reaction time (ms)
no2_class       float   %9.0g                 Classroom NO2 levels
                                                (ug/m3)
grade           byte    %9.0g      grade      Grade in school
sex             byte    %9.0g      sex        Sex
lbweight        byte    %18.0g     lowbw    * Low birthweight
breastfeed      byte    %19.0f     bfeed      Duration of
                                                breastfeeding
msmoke          byte    %10.0f     smoke    * Mother smoked during
                                                pregnancy
feducation      byte    %17.0g     edu        Father's education
                                                level
meducation      byte    %17.0g     edu        Mother's education
                                                level
overweight      byte    %32.0g     overwt   * Overweight by WHO/CDC
                                                definition
sev_home        float   %9.0g                 Home socio-economic
                                                vulnerability index
sev_school      float   %9.0g                 School socio-economic
                                                vulnerability index
age             float   %9.0g                 Age (years)
precip          double  %10.0g                Daily total
                                                precipitation
age0            double  %4.1f                 Age started school
siblings_old    byte    %1.0f                 Number of older
                                                siblings in house
siblings_young  byte    %1.0f                 Number of younger
                                                siblings in house
no2_home        float   %9.0g                 Home NO2 levels
                                                (ug/m3)
green_home      double  %10.0g                Home greenness (NDVI),
                                                300m buffer
noise_school    float   %9.0g                 School noise levels
                                                (dB)

Step 2: 进行 PO 回归

采用 poregress 命令进行回归，controls() 设定潜在的控制变量，本例中加入了 Step1 当中设定的连续型、因子型变量及其交乘项作为控制变量。

poregress react no2_class, controls(`ctrls') //PO 法

结果如下：

Partialing-out linear model          Number of obs                =      1,084
                                     Number of controls           =          0
                                     Number of selected controls  =          0
                                     Wald chi2(1)                 =      11.39
                                     Prob > chi2                  =     0.0007

------------------------------------------------------------------------------
             |               Robust
       react |      Coef.   Std. Err.      z    P>|z|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
   no2_class |   1.472821   .4363526     3.38   0.001      .617586    2.328057
------------------------------------------------------------------------------
Note: Chi-squared test is a Wald test of the coefficients of the variables
      of interest jointly equal to zero. Lassos select controls for model
      estimation. Type lassoinfo to see number of selected variables in each
      lasso.

将结果存储为 poplug

estimates store poplug

以上结果表明，每立方米再增加一微克 NO2 ，平均反应时间增加 1.47 毫秒。可以注意到，结果中只有核心解释变量 no2_class 的估计结果，其余控制变量并没有结果，这也是采用此类模型的特点。

PO 法采用的是多次 Lasso 方法来筛选控制变量。为说明 PO 方法的机理，我们考虑一个简单的线性模型：

其中 $y$ 为被解释变量，$d$ 为核心自变量，$X$ 为控制变量，要采用 PO 法估计 $\gamma$ ，一共有 5 步：

• 1.采用 Lasso 方法，用 $y$ 对 $X$ 进行回归，选出能够最佳预测 $y$ 的控制变量 ${\tilde{X}}_y$；
• 2.用 $y$ 对 ${\tilde{X}}_y$ 回归，得到残差 $\tilde{y}$;
• 3.采用 Lasso 方法，用 $d$ 对 $X$ 进行回归，选出能够最佳预测 $d$ 的控制变量 ${\tilde{X}}_d$；
• 4.用 $d$ 对 ${\tilde{X}}_d$ 回归，得到残差 $\tilde{d}$;
• 5.用 $\tilde{y}$ 对 $\tilde{d}$ 回归，得到 $\gamma$ 的估计值和标准误。

关于以上 5 步若想进一步了解，可参考 Chernozhukov, Hansen, and Spindler (2015a, b) 更为详细的分析及推导。

注： poregress 命令默认采用的是基于 plugin 方法的 Lasso 回归，也可以通过在选项中添加 selection() 来使用其他检验方法，这部分也可以可参考上篇推文 Stata:拉索开心读懂-Lasso 入门。

Step 3: 进一步分析

我们可以采用 lassocoef 命令看下每个 Lasso 回归哪些变量被选择。

lassocoef ( ., for(react)) ( ., for(no2_class))//变量筛选结果对比

结果为：

-------------------------------------------
                     |   react   no2_class
---------------------+---------------------
                 age |     x
                     |
  grade#c.green_home |
                4th  |     x
                     |
grade#c.noise_school |
                2nd  |     x
                     |
           sex#c.age |
                  0  |     x
                     |
    feducation#c.age |
                  4  |     x
                     |
          sev_school |               x
              precip |               x
            no2_home |               x
          green_home |               x
        noise_school |               x
                     |
  grade#c.sev_school |
                2nd  |               x
                     |
               _cons |     x         x
-------------------------------------------
Legend:
  b - base level
  e - empty cell
  o - omitted
  x - estimated

从上面的结果可以看出，对于 react 的回归中， age 和四个交乘项被选择。对于对于 no2_class 的回归中， sev_school、precip、no2_home、green_home、noise_school 和一个交乘项被选择。可见，两次 Lasso 回归中除了交乘项中的部分变量是相同的，其余被选择的控制变量差别是很大的。

当然，以上变量筛选的结果也可以用 lassoknots 命令，具体代码为：

lassoknots , for(react)//变量 react 的变量筛选
lassoknots , for(no2_class)//变量 no2_class 的变量筛选

结果为：

 lassoknots , for(react) //变量 react 的变量筛选

----------------------------------------------------------------------------------------------------
       |              No. of           |
       |             nonzero In-sample |                       Variables (A)dded, (R)emoved,
    ID |   lambda      coef. R-squared |                            or left (U)nchanged
-------+-------------------------------+------------------------------------------------------------
 (*) 1 | .1375306          5    0.1619 | A age               0.sex#c.age       3.grade#c.green_home
       |                               |   1.grade#c.noise_school
       |                               |   4.feducation#c.age
----------------------------------------------------------------------------------------------------
* lambda selected by plugin assuming heteroskedastic.

. lassoknots , for(no2_class) //变量 no2_class 的变量筛选

----------------------------------------------------------------------------------------------------
       |              No. of           |
       |             nonzero In-sample |                       Variables (A)dded, (R)emoved,
    ID |   lambda      coef. R-squared |                            or left (U)nchanged
-------+-------------------------------+------------------------------------------------------------
 (*) 1 | .1375306          6    0.3411 | A precip            no2_home          sev_school
       |                               |   green_home        noise_school      1.grade#c.sev_school
----------------------------------------------------------------------------------------------------
* lambda selected by plugin assuming heteroskedastic.

2.2 DS 法

DS 法是 PO 法的拓展。简单来讲， DS 法是通过引入额外的控制变量使得估计结果更稳健。和 PO 法的语法结构类似，采用命令 dsregress 即可实现，最后将估计结构存储为 dsplug。

. dsregress react no2_class, controls(`ctrls') //DS 法

. estimates store dsplug

Double-selection linear model  Number of obs               =  1,084
                               Number of controls          =      0
                               Number of selected controls =      0
                               Wald chi2(1)                =  11.39
                               Prob > chi2                 = 0.0007

-------------------------------------------------------------------
          |            Robust
    react |    Coef.  Std. Err.    z    P>|z|  [95% Conf. Interval]
----------+--------------------------------------------------------
no2_class | 1.472821  .4363526   3.38   0.001   .617586    2.328057
-------------------------------------------------------------------
Note: Chi-squared test is a Wald test of the coefficients of the variables
      of interest jointly equal to zero. Lassos select controls for model
      estimation. Type lassoinfo to see number of selected variables in each
      lasso.

DS 法估计的结果和 PO 法相似，结果的解读也是相似的，不再此处赘述。

对于 DS 法的机理，一共有 4 步，具体为：

• 采用 Lasso 方法，用 $y$ 对 $X$ 进行回归，选出能够最佳预测 $y$ 的控制变量 ${\tilde{X}}_y$；
• 采用 Lasso 方法，用 $d$ 对 $X$ 进行回归，选出能够最佳预测 $d$ 的控制变量 ${\tilde{X}}_d$；
• 令 ${\tilde{X}}_u$ 为 ${\tilde{X}}_y$ 和 ${\tilde{X}}_d$ 的合集； 用 $y$ 对 $d$ 和 ${\tilde{X}}_u$ 回归，得到 $\gamma$ 的估计值和标准误。

Belloni, Chernozhukov, and Wei (2016) 指出，尽管在大样本情况下 DS 法和 PO 法拥有同样的性质，但在他们的模拟中，前者表现更好。在有限样本情形下，DS 法表现更好可能是因为这种方法引入的控制变量是在一次回归当中而非两次单独的回归。

同理， DS 法也可以采用 lassocoef 和 lassoknots 命令进一步分析，因这部分分析和 PO 方法类似，故不再赘述。

2.3 XPO 法

XPO 法也被称为双重机器学习法 (DML)， Chernozhukov et al. (2018) 经过推导发现，这种方法相比于 PO 法而言在理论上具有更好的性质，并且在有限样本情形下也表现更好。其中，二者最为重要的差别是 XPO 法要求的稀疏条件更弱。在实际操作中，这意味着 XPO 法可以提供更可靠的统计推断，因为这种方法采用了样本拆分技术（split-sample techniques），从而可以引入更多的控制变量。

具体实操和 PO 法也是类似的。采用命令 xporegress 即可实现，最后将估计结果存储为 xpoplug。

xporegress react no2_class, controls(`ctrls') //XPO 法

estimates store xpoplug

结果为：

Cross-fit partialing-out    Number of obs                =  1,084
linear model                Number of controls           =      0
                            Number of selected controls  =      0
                            Number of folds in cross-fit =     10
                            Number of resamples          =      1
                            Wald chi2(1)                 =  11.39
                            Prob > chi2                  = 0.0007

------------------------------------------------------------------
          |            Robust
    react |    Coef.  Std. Err.   z    P>|z|  [95% Conf. Interval]
----------+-------------------------------------------------------
no2_class | 1.472821  .4363526  3.38   0.001   .617586    2.328057
------------------------------------------------------------------
Note: Chi-squared test is a Wald test of the coefficients of the variables
      of interest jointly equal to zero. Lassos select controls for model
      estimation. Type lassoinfo to see number of selected variables in each
      lasso.

估计和结果及相应的解释与 PO 法也是类似的，故也不再赘述。

关于 XPO 法的机理，首先是将样本进行拆分，拆分的子样份数叫做折（folds），以 2 折为例，一共包含 6 步：

1. 将样本拆分为 2 折，分别为 A 和 B；
2. 使用子样本 A 进行变量筛选：

   • 采用 Lasso 方法，用 $y$ 对 $X$ 进行回归，选出能够最佳预测 $y$ 的控制变量 ${\tilde{X}}_y$；
   • 用 $y$ 对 ${\tilde{X}}_y$ 回归，${\tilde{\beta }}_A$ 为估计系数；
   • 采用 Lasso 方法，用 $d$ 对 $X$ 进行回归，选出能够最佳预测 $d$ 的控制变量 ${\tilde{X}}_d$；
   • 用 $d$ 对 ${\tilde{X}}_d$ 回归，${\tilde{\delta }}_A$ 为估计系数。
3. 用 A 子样中估计系数，计算 B 子样中的 $y$ 和 $d$ 的残差，具体公式分别为：

   • $\tilde{y}=y-{\tilde{X}}_y{\tilde{\beta }}_A$
   • $\tilde{d}=y-{\tilde{X}}_d{\tilde{\delta }}_A$
4. 将上述操作 2-3 反过来再做一遍，即先使用子样本 B 进行变量筛选：

   • 采用 Lasso 方法，用 $y$ 对 $X$ 进行回归，选出能够最佳预测 $y$ 的控制变量 ${\tilde{X}}_y$；
   • 用 $y$ 对 ${\tilde{X}}_y$ 回归，${\tilde{\beta }}_B$ 为估计系数；
   • 采用 Lasso 方法，用 $d$ 对 $X$ 进行回归，选出能够最佳预测 $d$ 的控制变量 ${\tilde{X}}_d$；
   • 用 $d$ 对 ${\tilde{X}}_d$ 回归，${\tilde{\delta }}_B$ 为估计系数。
5. 用 B 子样中估计系数，计算 A 子样中的 $y$ 和 $d$ 的残差，具体公式分别为：

   • $\tilde{y}=y-{\tilde{X}}_y{\tilde{\beta }}_B$
   • $\tilde{d}=y-{\tilde{X}}_d{\tilde{\delta }}_B$
6. 用 $\tilde{y}$ 对 $\tilde{d}$ 回归，便可得到 $\gamma$ 的估计值及标准误。

以上是将样本拆分为 2 折的算法，10 折甚至 k 折的算法其实都是类似的。

3. 结论

综上，对比三种方法，最为推荐使用的是 XPO 法，因为这种方法在大样本和有限样本情形下都具有更好的性质，其缺点是计算会比较耗时。

4. 参考文献

• Belloni, A., D. Chen, V. Chernozhukov, and C. Hansen. 2012. Sparse models and methods for optimal instruments with an application to eminent domain. Econometrica 80: 2369–2429.-PDF-
• Belloni, A., V. Chernozhukov, and C. Hansen. 2014. Inference on treatment effects after selection among high-dimensional controls. Review of Economic Studies 81: 608–650.-PDF-
• Belloni, A., V. Chernozhukov, and Y. Wei. 2016. Post-selection inference for generalized linear models with many controls. Journal of Business & Economic Statistics 34: 606–619.-PDF-
• Chernozhukov, V., D. Chetverikov, M. Demirer, E. Duflo, C. Hansen, W. Newey, and J. Robins. 2018. Double/debiased machine learning for treatment and structural parameters. Econometrics Journal 21: C1–C68.-PDF-
• Chernozhukov, V., C. Hansen, and M. Spindler. 2015a. Post-selection and post-regularization inference in linear models with many controls and instruments. American Economic Review 105: 486–90.-PDF-
• ——. 2015b. Valid post-selection and post-regularization inference: An elementary, general approach. Annual Review of Economics 7: 649–688.-PDF-
• Leeb, H., and B. M. Pötscher. 2008. Sparse estimators and the oracle property, or the return of Hodges estimator. Journal of Econometrics 142: 201–211.-PDF-
• Sunyer, J., et al. 2017. "Traffic-related Air Pollution and Attention in Primary School Children." 28(2): 181-189. -PDF-
• Wooldridge, J. M. 2020. Introductory Econometrics: A Modern Approach. 7th ed. Boston, MA: Cengage-Learning.
• 杨继超， 连享会推文 Stata:拉索开心读懂-Lasso 入门

5. 相关推文

Note：产生如下推文列表的命令为：
  lianxh lasso 高维 随机推断
安装最新版 lianxh 命令：
  ssc install lianxh, replace

• 专题：Stata命令

  • Stata新命令-pdslasso：众多控制变量和工具变量如何挑选？
• 专题：Stata绘图

  • Stata绘图：随机推断中的系数可视化
• 专题：回归分析

  • Stata：拉索回归和岭回归-(Ridge,-Lasso)-简介
  • Stata Blogs - An introduction to the lasso in Stata (拉索回归简介)
• 专题：内生性-因果推断

  • Stata：ritest-随机推断（Randomization Inference）
• 专题：机器学习

  • Stata：拉索开心读懂-Lasso入门

相关课程

连享会-直播课 上线了！
http://lianxh.duanshu.com

免费公开课：

• 直击面板数据模型 - 连玉君，时长：1小时40分钟，课程主页
• Stata 33 讲 - 连玉君, 每讲 15 分钟.
• Stata 小白的取经之路 - 龙志能，时长：2 小时，课程主页
• 部分直播课 课程资料下载 (PPT，dofiles等)

课程一览

支持回看

Note: 部分课程的资料，PPT 等可以前往 连享会-直播课 主页查看，下载。

关于我们

• Stata连享会 由中山大学连玉君老师团队创办，定期分享实证分析经验。直播间 有很多视频课程，可以随时观看。
• 连享会-主页 和 知乎专栏，400+ 推文，实证分析不再抓狂。
• 公众号关键词搜索/回复 功能已经上线。大家可以在公众号左下角点击键盘图标，输入简要关键词，以便快速呈现历史推文，获取工具软件和数据下载。常见关键词：课程, 直播, 视频, 客服, 模型设定, 研究设计, stata, plus, 绘图, 编程, 面板, 论文重现, 可视化, RDD, DID, PSM, 合成控制法 等

连享会小程序：扫一扫，看推文，看视频……

扫码加入连享会微信群，提问交流更方便

✏ 连享会学习群-常见问题解答汇总：
✨ https://gitee.com/arlionn/WD

New！ lianxh 命令发布了：
随时搜索连享会推文、Stata 资源，安装命令如下：
  . ssc install lianxh
使用详情参见帮助文件 (有惊喜)：
  . help lianxh