Stata：标准误！标准误！

发布时间：2022-01-28 阅读 5687

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

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

New！ lianxh 命令发布了：
随时搜索推文、Stata 资源。安装：
. ssc install lianxh
详情参见帮助文件 (有惊喜)：
. help lianxh
连享会新命令：cnssc, ihelp, rdbalance, gitee, installpkg

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

⛳ 课程主页： https://gitee.com/lianxh/Course

⛳ Stata 系列推文：

全部 | Stata入门 | Stata教程 | Stata资源 | Stata命令
计量专题 | 论文写作 | 数据分享 | 专题课程
结果输出 | Stata绘图 | 数据处理 | Stata程序
回归分析 | 面板数据 | 交乘项-调节 | IV-GMM
内生性-因果推断 | 倍分法DID | 断点回归RDD | PSM-Matching | 合成控制法
Probit-Logit | 时间序列 | 空间计量 | 分位数回归 | 生存分析 | SFA-DEA
文本分析-爬虫 | Python-R-Matlab | 机器学习
Markdown | 工具软件 | 其它

☝ PDF下载 - 推文合集

作者：苏湘晗 (中山大学)
邮箱：945294986@qq.com

编者按：本文主要摘译自下文，特此致谢！
Source：Michler J D, Josephson A. Recent Developments in Inference: Practicalities for Applied Economics[J]. arXiv preprint arXiv:2107.09736, 2021. -PDF-

1. 标准误校正新发展

近年来，应用经济学领域在校正标准误方面取得了许多新进展。本文将综合这些进展，并结合相应的 Stata 命令来解决传统 OLS 回归中出现的问题，如异方差、聚类与序列相关性，以及检验多重假设。除此之外，还涉及数值方法的改进，如自举法以及随机化推理。

关于计算正确的标准误和检验统计量，我们提出三项具体建议：

首先，应用经济学家应该明确知道数据在假设检验的过程中会出现的问题，以及问题的根源；
其次，现代统计软件和计算功能十分强大，各类解决方法都有对应的 Stata 命令，应加强实际操作；
最后，由于复杂的抽样策略和研究设计方法，标准误有较大误差。在这些情况下，运用自抽样 (bootstrap) 可以计算出估计量或检验统计量分布的近似值，这个值接近于基于渐近理论计算出的值甚至更加准确。所以，通过自抽样进行渐近细化应该成为一种标准做法。

2. 传统假设检验存在的问题

2.1 异方差问题

传统假设检验的一个基本假设是，误差是同方差的，即：

E = [ε_{i}^{2} | X_{i}] = Ω

这个默认的假设被嵌入 Stata 计算的标准误中。然而，在大多数实际应用中，同方差性的假设不太可能成立，残差将随着 $X_{i}$ 值的大小而变化。

在大样本中，我们可以依赖于 Eicker-Huber-White (EHW) 方差估计，它是渐近一致的。为了得到 EHW 稳健标准误，我们只需将 OLS 方差估计为：

V_{E H W} (\hat{β}) = (X_{i}^{'} X_{i})^{- 1} (\sum_{i} Ω_{i} X_{i} X_{i}^{'}) (X_{i}^{'} X_{i})^{- 1}

其中， $Ω$ 不再是一个常数，而是估计的残差的方差：

{HC}_{0} : {\hat{Ω}}_{i} = {\hat{ε}}_{i}^{2}

在这里，HC 代表异方差是一致的 (heteroskedasticity-consistent)。但是，如果残差确实是同方差的，那么 EHW 标准误实际上将比传统标准误有更大的偏差。即使使用了适当的标准误差，传统的标准误差和 EHW 标准误差都可能会渐近偏置，导致 ${\hat{Ω}}_{i}$ 偏小。

MacKinnon 和 White (1985) 为 EHW 标准误中的偏差提供了三种可能的修正：

2.1.1 vce(robust)

第一种修正是根据样本量 (N) 和自由度 (k) 简单调整自由度：

{HC}_{1} : {\hat{Ω}}_{i} = \frac{N}{N - k} {\hat{ε}}_{i}^{2}

以 nlsw88.dta 数据集为例，我们用 vce(robust) 修正，具体命令如下：

. sysuse nlsw88.dta, clear
. global x "hours ttl_exp tenure south collgrad married" 
. reg wage $x, vce(robust)

Linear regression                        Number of obs =   2,227
                                         F(6, 2220)    =   75.34
                                         Prob > F      =  0.0000
                                         R-squared     =  0.1517
                                         Root MSE      =  5.3183
-----------------------------------------------------------------
         |               Robust
    wage | Coefficient  std. err.      t    P>|t|     [95% CI]
---------+-------------------------------------------------------
   hours |      0.053      0.011     4.78   0.000   0.031   0.075
 ttl_exp |      0.252      0.030     8.48   0.000   0.193   0.310
  tenure |      0.022      0.025     0.89   0.372  -0.026   0.070
   south |     -1.528      0.226    -6.75   0.000  -1.973  -1.084
collgrad |      3.167      0.300    10.56   0.000   2.579   3.755
 married |     -0.200      0.248    -0.80   0.422  -0.687   0.288
   _cons |      2.546      0.455     5.59   0.000   1.653   3.439
-----------------------------------------------------------------

. est store m1

2.1.2 vce(hc2)

第二种修正是根据 $X_{i}$ 的杠杆率来调整方差。我们将杠杆率定义为：

h_{i i} = h_{i}^{'} h_{i} = X_{i}^{'} (X^{'} X)^{- 1} X_{i}

如果数据中的所有 $X_{i}$ 与平均值的距离相等，那么 ${HC}_{1}$ 的修正是正确的。但是，如果 $X_{i}$ 有不同的杠杆率 (绝大多数情况)，那么 ${HC}_{1}$ 就并不能完全纠正偏差。所以，调整杠杆率为：

{HC}_{2} : {\hat{Ω}}_{i} = \frac{N}{1 - h_{i i}} {\hat{ε}}_{i}^{2}

我们用 vce(hc2) 修正，具体命令如下：

. reg wage $x, vce(hc2)

Linear regression                    Number of obs =   2,227
                                     F(6, 2220)    =   75.27
                                     Prob > F      =  0.0000
                                     R-squared     =  0.1517
                                     Root MSE      =  5.3183
-------------------------------------------------------------
         |         Robust HC2
    wage |   beta   std. err.      t    P>|t|     [95% CI]
---------+---------------------------------------------------
   hours |  0.053      0.011     4.77   0.000    0.031  0.075
 ttl_exp |  0.252      0.030     8.48   0.000    0.193  0.310
  tenure |  0.022      0.025     0.89   0.372   -0.026  0.070
   south | -1.528      0.226    -6.75   0.000   -1.973 -1.084
collgrad |  3.167      0.300    10.56   0.000    2.579  3.755
 married | -0.200      0.248    -0.80   0.422   -0.687  0.288
   _cons |  2.546      0.456     5.59   0.000    1.653  3.440
-------------------------------------------------------------

. est store m2

2.1.3 vce(hc3)

第三种修正是在第二种修正的基础上将分母平方：

{HC}_{3} : {\hat{Ω}}_{i} = \frac{N}{(1 - h_{i i})^{2}} {\hat{ε}}_{i}^{2}

当残差确实是异方差时，这个修正比其他修正表现得更好；但如果残差是同方差时，这个估计将会有偏差。Angrist 和 Pischke (2009) 表明， ${HC}_{3}$ 校正倾向于产生大于真实方差的估计值。

我们用 vce(hc3) 修正，具体命令如下：

. reg wage $x, vce(hc3)

Linear regression                     Number of obs =   2,227
                                      F(6, 2220)    =   74.97
                                      Prob > F      =  0.0000
                                      R-squared     =  0.1517
                                      Root MSE      =  5.3183
-------------------------------------------------------------
         |         Robust HC3                               
    wage |   beta   std. err.     t    P>|t|      [95% CI]   
---------+---------------------------------------------------
   hours |  0.053      0.011    4.76   0.000   0.031    0.075
 ttl_exp |  0.252      0.030    8.46   0.000   0.193    0.310
  tenure |  0.022      0.025    0.89   0.373  -0.026    0.070
   south | -1.528      0.227   -6.74   0.000  -1.973   -1.084
collgrad |  3.167      0.301   10.54   0.000   2.578    3.757
 married | -0.200      0.249   -0.80   0.422  -0.688    0.288
   _cons |  2.546      0.457    5.58   0.000   1.651    3.442
-------------------------------------------------------------

. est store m3

三种修正的回归结果如下所示。根据对比结果可以发现，回归的系数大小不变，但 SE 小幅度变大。

. esttab m1 m2 m3, mtitle(vce(robust) vce(hc2) vce(hc3)) nogap se(%6.4f)  ///
>     star(* 0.1 ** 0.05 *** 0.01)

------------------------------------------------------------
                      (1)             (2)             (3)   
              vce(robust)        vce(hc2)        vce(hc3)   
------------------------------------------------------------
hours              0.0531***       0.0531***       0.0531***
                 (0.0111)        (0.0111)        (0.0112)   
ttl_exp             0.252***        0.252***        0.252***
                 (0.0297)        (0.0297)        (0.0297)   
tenure             0.0219          0.0219          0.0219   
                 (0.0246)        (0.0246)        (0.0246)   
south              -1.528***       -1.528***       -1.528***
                 (0.2265)        (0.2265)        (0.2269)   
collgrad            3.167***        3.167***        3.167***
                 (0.2999)        (0.3000)        (0.3006)   
married            -0.200          -0.200          -0.200   
                 (0.2485)        (0.2485)        (0.2489)   
_cons               2.546***        2.546***        2.546***
                 (0.4553)        (0.4556)        (0.4566)   
------------------------------------------------------------
N                    2227            2227            2227   
------------------------------------------------------------
Standard errors in parentheses
* p<0.1, ** p<0.05, *** p<0.01

2.1.4 reg_sandwich

第四种对 EHW 标准误的修正是 Bell 和 McCaffffrey (2002) 提出的 BM 校正。BM 校正不使用正态分布或自由度调整后的 $t$ 分布，而是将方差估计量的矩与 $χ^{2}$ 分布相匹配。Stata 命令是 reg_sandwich，可以对小样本数据进行聚类稳健标准误纠偏处理的线性回归命令，它提供了聚类稳健标准误估计回归模型。

* 命令安装
ssc install reg_sandwich, replace
* 命令语法
reg_sandwich depvar indepvars, cluster(varname)

其中，depvar 为因变量，indepvars 为自变量和控制变量，cluster(varname) 是指用于聚类纠偏的变量。首先，使用官方命令对行业 industry 进行聚类标准误纠偏，但未作小样本纠偏处理，具体命令如下：

. reg wage $x, cluster(industry)  

Linear regression                               Number of obs     =      2,213
                                                F(6, 11)          =      81.70
                                                Prob > F          =     0.0000
                                                R-squared         =     0.1499
                                                Root MSE          =     5.3317

                              (Std. err. adjusted for 12 clusters in industry)
------------------------------------------------------------------------------
             |               Robust
        wage | Coefficient  std. err.      t    P>|t|     [95% conf. interval]
-------------+----------------------------------------------------------------
       hours |      0.053      0.018     2.98   0.012        0.014       0.092
     ttl_exp |      0.252      0.040     6.30   0.000        0.164       0.339
      tenure |      0.020      0.027     0.75   0.471       -0.040       0.081
       south |     -1.522      0.239    -6.37   0.000       -2.047      -0.996
    collgrad |      3.175      0.390     8.15   0.000        2.317       4.032
     married |     -0.202      0.268    -0.75   0.466       -0.792       0.387
       _cons |      2.574      0.924     2.78   0.018        0.540       4.609
------------------------------------------------------------------------------

. est store m4

其次，对小样本进行纠偏处理，回归命令如下：

. reg_sandwich wage $x, cluster(industry) 

Robust Small Sample Corrected standard error estimation using OLS
                                                      Number of obs =     2227
                                                      R-squared     =   0.1499
                                                      Adj R-squared =   0.1476
                                                      Root MSE      =   5.3317

                                  (Std. Err. adjusted for 12 clusters in industry)
------------------------------------------------------------------------------
             |              Robust
        wage |      Coef.   Std. Err.  dfs       p-value   [95%Conf. Interval]
-------------+----------------------------------------------------------------
       hours |      0.053    0.02       3.54     0.0611    -0.0042   0.1097
     ttl_exp |      0.252    0.04       4.70     0.0020    0.1445    0.3586
      tenure |      0.020    0.03       4.28     0.5206    -0.0585   0.0993
       south |      -1.522   0.24       4.22     0.0028    -2.1828   -0.8604
    collgrad |      3.175    0.50       2.69     0.0110    1.4636    4.8860
     married |      -0.202   0.30       4.58     0.5363    -1.0026   0.5980
       _cons |      2.574    1.04       3.92     0.0705    -0.3475   5.4964
------------------------------------------------------------------------------

. est store m5

可以使用以下命令对比分析纠偏前后的回归结果：

. esttab m4 m5, mtitle(reg reg_sandwich) nogap se(%6.4f)  ///
>     star(* 0.1 ** 0.05 *** 0.01)

--------------------------------------------
                      (1)             (2)   
                      reg    reg_sandwich   
--------------------------------------------
hours              0.0528**        0.0528***
                 (0.0177)        (0.0195)   
ttl_exp             0.252***        0.252***
                 (0.0399)        (0.0409)   
tenure             0.0204          0.0204   
                 (0.0273)        (0.0292)   
south              -1.522***       -1.522***
                 (0.2387)        (0.2431)   
collgrad            3.175***        3.175***
                 (0.3895)        (0.5037)   
married            -0.202          -0.202   
                 (0.2679)        (0.3027)   
_cons               2.574**         2.574** 
                 (0.9245)        (1.0443)   
--------------------------------------------
N                    2213            2227   
--------------------------------------------
Standard errors in parentheses
* p<0.1, ** p<0.05, *** p<0.01

根据纠偏前后的结果对比可以发现，纠偏前的 SE 通常偏小，导致 $t$ 值偏大，回归系数的显著性较低。通过 ·reg_sandwich· 命令对回归进行纠偏处理后，回归的系数大小不变，而显著性将提升，降低了假设被错误拒绝的可能性。

2.1.5 edfreg

第五种修正是由 Young (2016b) 提出的有效自由度 (effective degrees of freedom，简称 EDOF) 校正，以消除偏差并调整 $t$ 分布，保证检验统计量的方差不会膨胀。

* 命令安装
ssc install edfreg, replace
* 命令语法
edfreg depvar indepvars, cluster(varname)

. edfreg wage $x, cluster(industry) 

Linear regression                               Number of obs     =      2,213
                                                F(6, 11)          =      81.70
                                                Prob > F          =     0.0000
                                                R-squared         =     0.1499
                                                Root MSE          =     5.3317
                              (Std. err. adjusted for 12 clusters in industry)
------------------------------------------------------------------------------
             |               Robust
        wage | Coefficient  std. err.      t    P>|t|     [95% conf. interval]
-------------+----------------------------------------------------------------
       hours |      0.053      0.018     2.98   0.012        0.014       0.092
     ttl_exp |      0.252      0.040     6.30   0.000        0.164       0.339
      tenure |      0.020      0.027     0.75   0.471       -0.040       0.081
       south |     -1.522      0.239    -6.37   0.000       -2.047      -0.996
    collgrad |      3.175      0.390     8.15   0.000        2.317       4.032
     married |     -0.202      0.268    -0.75   0.466       -0.792       0.387
       _cons |      2.574      0.924     2.78   0.018        0.540       4.609
------------------------------------------------------------------------------
 
Nominal degrees of freedom == 11

          |                 Bias     Effective
          |               Adjusted    Degrees
 Variable |     Coef.    Std. Err.  of Freedom  t   P>|t|  [95% Conf. Interval]
-------------+-----------------------------------------------------------------
    hours |  .05277815  .01952374     4.5     2.70  0.048  .00079738  .10475892
  ttl_exp |  .25153637  .04236519     5.3     5.94  0.002  .14441144  .35866131
   tenure |  .02038819  .02928611     5.0     0.70  0.517  -.0548942  .09567057
    south | -1.5216007  .25593342     5.0    -5.95  0.002 -2.1804233 -.86277802
 collgrad |  3.1747799  .47055587     4.3     6.75  0.002  1.9085245  4.4410352
  married | -.20226289  .28530694     5.2    -0.71  0.509  -.9272895  .52276373
    _cons |  2.5744341  1.0081355     4.7     2.55  0.054 -.06791532  5.2167836

. est store m6

根据纠偏前后的结果对比可以发现，纠偏后的 $t$ 值小于纠偏前的 $t$ 值，显著性提升，降低了假设被错误拒绝的可能性。

2.2 聚类与序列相关性

除了同方差的假设之外，在应用经济学研究中，传统的假设是每个观察结果都是从某些群体中随机抽取的。这种独立性假设意味着回归中的残差将彼此不相关。然而，无论是从抽样方法还是实验的研究设计，大多数微观数据都是聚类、分组或分层的。随着时间的推移，从同一单位进行的重复采样会产生序列相关性，这意味着单位内的残差将是相关的。

2.2.1 vce(cluster clustvar)

组内相关的标准校正是由 Liang 和 Zeger (1986) 提出的，随后 LZ 将 White (1980) 的稳健协方差矩阵推广到分组数据中。这种修正让标准误对于未指定形式的聚类内部相关性是稳健的。为了得到 LZ 类稳健标准误差，只需估计 OLS 方差为：

V_{L Z} (\hat{β}) = (X_{i}^{'} X_{i})^{- 1} (\sum_{c} X_{c}^{'} Ω_{c} X_{c}) (X_{i}^{'} X_{i})^{- 1}

其中， $Ω_{c}$ 是集群 $c$ 中对应单位的子矩阵。为了得到方差的一致估计，我们使用估计的残差：

{\hat{Ω}}_{c} = a {\hat{ε}}_{c} {\hat{ε}}_{c}^{'}

用 vce(cluster industry)修正，对行业 industry 进行聚类标准误纠偏，具体命令如下：

. reg wage $x, vce(cluster industry)

Linear regression                               Number of obs     =      2,213
                                                F(6, 11)          =      81.70
                                                Prob > F          =     0.0000
                                                R-squared         =     0.1499
                                                Root MSE          =     5.3317

                              (Std. err. adjusted for 12 clusters in industry)
------------------------------------------------------------------------------
             |               Robust
        wage | Coefficient  std. err.      t    P>|t|     [95% conf. interval]
-------------+----------------------------------------------------------------
       hours |      0.053      0.018     2.98   0.012        0.014       0.092
     ttl_exp |      0.252      0.040     6.30   0.000        0.164       0.339
      tenure |      0.020      0.027     0.75   0.471       -0.040       0.081
       south |     -1.522      0.239    -6.37   0.000       -2.047      -0.996
    collgrad |      3.175      0.390     8.15   0.000        2.317       4.032
     married |     -0.202      0.268    -0.75   0.466       -0.792       0.387
       _cons |      2.574      0.924     2.78   0.018        0.540       4.609
------------------------------------------------------------------------------

2.2.2 clusteff

在大样本、横截面数据中，聚类校正是直接的。然而，如果我们的分组很少，或者随着时间的推移对数据进行分组，这种方法就会失效。“太少” 组的问题类似于异方差的小样本问题。从根本上说，LZ 修正依赖于渐近理论，要求簇的数量趋于无穷大的一致性。

clusteff 使用自变量向量、聚类变量检查集群大小和杠杆的异质性程度，然后计算有效的集群数量。

* 命令安装
net install st0531.pkg, replace

* 命令语法
clusteff varlist, cluster(varname)

计算 industry 集群的有效数量，使用命令及结果如下：

. clusteff $x, cluster(industry)
Warning: no variable specified for the test option.
Calculating G* under the null hypothesis: b1 = 0.
Number of clusters: 12
Estimated effective number of clusters: 2.0841147
Warning: G* estimated to be below 50.

2.2.3 ivreg2

Cameron 和 Miller (2015) 提出了一个更好的解决方案，即实现多路聚类，有效地考虑了回归 (分组数据) 中的聚类内相关性，以及残差中的聚类内相关性 (串行相关性)。多路聚类可以使用 ivreg2 实现。

ivreg2 先在第一级 (组) 计算聚类稳健方差-协方差矩阵；接着在第二级 (时间) 计算矩阵，然后在两个级别的交点 (组×时间) 计算矩阵；最后，将前两个矩阵相加后减去第三个矩阵。虽然这种方法允许沿着两个维度进行集群，但它需要在两个级别中拥有 “足够多” 的集群。

* 命令安装
ssc install ivreg2, replace

* 命令语法
ivreg2 depvar [varlist1] (varlist2=varlist_iv) [if] [in] [weight]
       [, gmm gmm2s cluster(varname) partial(varlist)]

在这里我们使用 griliches76.dta 数据集，并以 year 为集群，进行面板 2-step GMM 估计。

. use http://fmwww.bc.edu/ec-p/data/hayashi/griliches76.dta, clear
. xi i.year
. ivreg2 lw s expr tenure rns _I* (iq=kww age), cluster(year) partial(_I*) gmm2s

2-Step GMM estimation
---------------------
Number of clusters (year) =          7                Number of obs =      758
                                                      F(  5,     6) =   115.34
                                                      Prob > F      =   0.0000
Total (centered) SS     =  104.5637468                Centered R2   =   0.2126
Total (uncentered) SS   =  104.5637468                Uncentered R2 =   0.2126
Residual SS             =  82.33325095                Root MSE      =    .3296
------------------------------------------------------------------------------
             |               Robust
          lw | Coefficient  std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
          iq |      0.003      0.004     0.67   0.502       -0.006       0.012
           s |      0.072      0.017     4.17   0.000        0.038       0.106
        expr |      0.038      0.011     3.29   0.001        0.015       0.060
      tenure |      0.047      0.005     8.74   0.000        0.036       0.057
         rns |     -0.107      0.041    -2.63   0.009       -0.187      -0.027
------------------------------------------------------------------------------
Underidentification test (Kleibergen-Paap rk LM statistic):              5.653
                                                   Chi-sq(2) P-val =    0.0592
------------------------------------------------------------------------------
Weak identification test (Cragg-Donald Wald F statistic):               26.253
                         (Kleibergen-Paap rk Wald F statistic):         43.005
Stock-Yogo weak ID test critical values: 10% maximal IV size             19.93
                                         15% maximal IV size             11.59
                                         20% maximal IV size              8.75
                                         25% maximal IV size              7.25
Source: Stock-Yogo (2005).  Reproduced by permission.
NB: Critical values are for Cragg-Donald F statistic and i.i.d. errors.
------------------------------------------------------------------------------
Hansen J statistic (overidentification test of all instruments):         5.251
                                                   Chi-sq(1) P-val =    0.0219
------------------------------------------------------------------------------
Instrumented:         iq
Included instruments: s expr tenure rns
Excluded instruments: kww age
Partialled-out:       _Iyear_67 _Iyear_68 _Iyear_69 _Iyear_70 _Iyear_71
                      _Iyear_73 _cons
                      nb: total SS, model F and R2s are after partialling-out;
                          any small-sample adjustments include partialled-out
                          variables in regressor count K
------------------------------------------------------------------------------

2.3 多重假设

假设人们想要测试一项政策在减少食物损失和浪费方面是否有统计学上的显著效果，有两种不同的方式来衡量。其中一个是执行两个独立的测试，当 $α = 0.5$ 时，不犯第一类错误的概率是 $0.95 \times 0.95 = 0.9025$ 。当实际上是不显著时，得出有统计上显著影响的概率是 $1 - 0.9025 = 0.0975$ ，而不是 $0.05$ 。通过检验两个假设，研究人员可能会得出不同的结论。在进行推理时，需要考虑到这一事实。

多重假设可能会导致多重比较谬误 (family-wise error rate，简称 FWER) 或错误发现率 (false discovery rate，简称 FDR)。FWER 是指在一系列检验中至少有一次发现不显著的概率，FDR 是指在已经完成的检验中至少有一个发现不显著的概率。以下校正以及 Stata 命令能帮助人们更好地解决多重假设问题：

2.3.1 wyoung

Westfall 和 Young (1993) 提出了一种降重方法，它允许 $p$ 值之间的依赖关系，使用自举法来估计假设之间的依赖关系。然后，这种依赖结构被用来调整降重过程，假设子集枢轴性，以减少必要假设之间的比较次数，从而简化计算过程。Westfall-Young 过程是渐近最优的。

* 命令安装
ssc install wyoung, replace
* 命令语法
wyoung varlist, cmd(model) familyp(varlist) bootstraps(#) 
       [seed(#) strata(varlist) cluster(varlist) force
       subgroup(varname) singlestep detail noresampling familypexp replace]

这里我们以 auto.dta 数据集为例，具体命令如下：

. sysuse auto.dta, clear
. wyoung mpg headroom turn, cmd(regress OUTCOMEVAR displacement length,  ///
>     cluster(rep78)) cluster(rep78) familyp(displacement)               ///
>     bootstraps(100) seed(20)

familyp: displacement
    regress mpg displacement length, cluster(rep78)
    regress headroom displacement length, cluster(rep78)
    regress turn displacement length, cluster(rep78)

     +--------------------------------------------------------------------------+
  1. | k |           model   |  outcome |      familyp |       coef |    stderr |
     | 1 |    regress mpg~   |      mpg | displacement | -.01239323 | .00449325 |
     |--------------------------------------------------------------------------|
     |            p   |    pwyoung    |        pbonf    |           psidak      |
     |    .05094428   |        .45    |    .15283285    |        .14517911      |
     +--------------------------------------------------------------------------+
     +--------------------------------------------------------------------------+
  2. | k |             model |  outcome |      familyp |       coef |    stderr |
     | 2 | regress headroom~ | headroom | displacement |  .00087747 |  .0007628 |
     |--------------------------------------------------------------------------|
     |         p      |     pwyoun    |      pbonf      |        psidak         |
     | .31410222      |        .53    |  .31410222      |     .31410222         |
     +--------------------------------------------------------------------------+
     +--------------------------------------------------------------------------+
  3. | k |            model |  outcome |      familyp |       coef |     stderr |
     | 3 |    regress turn~ |     turn | displacement |  .01268401 |  .00635871 |
     |--------------------------------------------------------------------------|
     |         p      |    pwyoung     |     pbonf      |         psidak        |
     | .11681525      |        .53     |  .2336305      |       .2199847        |
     +--------------------------------------------------------------------------+

2.3.2 rwolf

Romano 和 Wolf (2005a，2005b) 提出了对 Westfall-Young 修正的推广。和 Westfall-Young (1993) 一样，Romano-Wolf 校正使用了一个自举重采样程序，利用了 $p$ 值之间的依赖性，能更有效地检测真正错误的零假设。该过程可以通过用户编写的命令 rwolf 在 Stata 中实现 (Clarke，2016)。

* 命令安装
ssc install rwolf, replace 
* 命令语法
rwolf depvars [if] [in] [weight], [options]

. sysuse auto.dta, clear
. rwolf headroom turn price rep78, indepvar(weight) controls(trunk mpg) reps(250)

Independent variable:  weight
Outcome variables:   headroom turn price rep78
Number of resamples: 250
------------------------------------------------------------------------------
   Outcome Variable | Model p-value    Resample p-value    Romano-Wolf p-value
--------------------+---------------------------------------------------------
           headroom |    0.6719             0.6494              0.6494
               turn |    0.0000             0.0040              0.0040
              price |    0.0075             0.0398              0.0598
              rep78 |    0.0998             0.1155              0.1673
------------------------------------------------------------------------------

2.3.3 mhtreg

List 等人 (2019) 开发了一个依赖于类似的自举重采样程序的 FWER 校正，提出了一个额外的假设，即检验过程是在个体水平上，通过简单的随机抽样完成的。Seidel 和 Xu (2016) 提供了 mhtexp 命令来实现 List-Shaikh-Xu 校正。

但是，该命令不允许在不同的回归中使用不同的控件。所以 Steinmayr (2020) 提供了一个更通用的 mhtexp 命令版本，称为 mhtreg，它允许在不同的回归中包含不同的控制，以及集群随机化。

* 命令安装
ssc install mhtreg, replace
* 命令语法
mhtreg (depvar indepvars [if]) (depvar indepvars [if]) [(depvar indepvars if) ] ... [, options]

. sysuse auto, clear
. mhtreg (price length trunk) (mpg length trunk)

             |       hyp       coef          p    pthm3_1      pbonf      pholm 
-------------+------------------------------------------------------------------
          r1 |         1   57.08706   .0036667   .0036667   .0073333   .0036667 
          r2 |         2   -.205419   .0003333   .0003333   .0006667   .0006667

3. 数值方法的渐进改进

随着计算成本的降低，基于数值方法的渐近改进变得越来越普遍 (Mackinnon，2006)。其中最主要的是自举法，这是由 Efron (1979，1982) 创造的一个术语。这个词来源于成语 “通过自己的引导自己”，意思是使用现有的资源来改善自己的处境。通过从现有数据中重新采样，以模拟底层数据，自举程序可以提供对人口统计数据的细化。

自举提供了一种替代传统推理的方法，它依赖于渐近公式，在许多情况下，自举改进了未知精确有限样本分布的渐近近似 (Horowitz，2019)。一个自举样本是从自己的数据中提取的，其中原始样本被视为它作为总体处理，再重复抽取一个新的随机样本，即自举样本。我们期望从数据中提取并构建的自举分布能提供一个很好的近似抽样分布。

3.1 经典自举

经典的自举法 (Classic Bootstrapping) 对函数形式的分布不做任何假设。从 $n$ 个观测值的原始样本开始，我们随机抽取一对 $(y_{i}, x_{i})$ ，即因变量及其相应的协变量。将这对值保存为新引导样本中的第一个观察，并替代原始样本，再随机抽取一个新的值并放置在自举样本中。这个过程重复了 $n$ 次。因为引导样本从原始样本中随机抽取并进行替换，所以一些原始观察结果会出现多次，而另一些则根本不会出现。

然后，计算我们感兴趣的系数的平均值和方差，以及其他的检验统计量，并保存它们。然后反复重复此过程，以构建复制统计信息的自举分布。

自举标准误被计算为：

\hat{σ} = {\frac{1}{r - 1} \sum ({\hat{β}}_{i} - \bar{β})^{2}}^{\frac{1}{2}}

$r$ 是重复次数， ${\hat{β}}_{i}$ 是第 $i$ 个自举样本的估计系数， $\bar{β}$ 是自举分布中该系数的平均值。要生成引导的 $t$ 统计量，可以计算：

t_{i} = \frac{{\hat{β}}_{i} - \bar{β}}{{\hat{σ}}_{i}}

其中 $t_{i}$ 是第 $i$ 个引导样本的 $t$ 统计量。按绝对值排序 $t_{i}$ ，使 | $t_{1}$ |≤···≤| $t_{r}$ |，确保对称临界值。然后从这些有序的元素中选择适当的临界值。所以，对于双面 95% 检验， $r$ =10000，有序元素是第 9500 个。还可以使用这个过程为假设检验创建自举置信区间，而不是标准表中的置信区间。

经典的自举法对于使用 Stata 进行估计相对简单，具体如下：

3.1.1 vce(bootstrap)

生成自举标准误的最简单方法是使用 vce(bootstrap)，它会自动处理集群和其他模型细节。我们使用 nlsw88.dta 数据集，用 vce(bootstrap) 自举，相关命令如下：

. sysuse "nlsw88.dta", clear
. global x "hours ttl_exp tenure south collgrad married" 
. reg wage $x, vce(bootstrap)

Linear regression                                       Number of obs =  2,227
                                                        Replications  =     50
                                                        Wald chi2(6)  = 506.92
                                                        Prob > chi2   = 0.0000
                                                        R-squared     = 0.1517
                                                        Adj R-squared = 0.1494
                                                        Root MSE      = 5.3183
------------------------------------------------------------------------------
             |   Observed   Bootstrap                         Normal-based
        wage | coefficient  std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
       hours |      0.053      0.011     4.73   0.000        0.031       0.075
     ttl_exp |      0.252      0.025    10.13   0.000        0.203       0.300
      tenure |      0.022      0.020     1.10   0.272       -0.017       0.061
       south |     -1.528      0.232    -6.60   0.000       -1.982      -1.075
    collgrad |      3.167      0.292    10.84   0.000        2.595       3.739
     married |     -0.200      0.282    -0.71   0.478       -0.752       0.352
       _cons |      2.546      0.459     5.55   0.000        1.646       3.446
------------------------------------------------------------------------------

3.1.2 bsample

bsample 从数据集替换和自举出新样本，其语法格式如下：

* 命令语法
bsample [exp] [if] [in] [, options]

其中，[exp] 指定了样本的大小，使用命令如下：

bsample 200

3.2 残差和原始自举

在残差自举 (Residual Bootstrapping) 中，协变量值是固定的，并从残差分布中抽取样本，根据每个观测值的预测值和残差提取生成因变量的新估计值。与经典自举一样，重复这个过程以建立一个残差的自举分布。该方法模拟了一个具有非随机误差的样本，并确保了协变量和回归残差的同等分布。

而原始自举 (Wild Bootstrapping) 是一种放宽独立性假设以允许异方差性的方法，适用于残差序列相关或存在簇内相关时使用。原始自举可以在有限的矩下有效地匹配未知的精确有限样本分布，使推理更加可靠。

boottest 首先去掉样本中每个观测值的平均值，从而使中心分布的平均值为零。接下来，从这个中心分布中选取一个替换的样本，对于每个随机样本，计算一个 $t$ 统计量，就像经典自举一样。由此，在假设总体均值为零的前提下，得到随机抽样期望的 $t$ 统计量的自举分布。

* 命令安装
ssc install boottest, replace

. sysuse nlsw88.dta, clear
. regress wage tenure ttl_exp collgrad, cluster(industry)

Linear regression                               Number of obs     =      2,217
                                                F(3, 11)          =      56.79
                                                Prob > F          =     0.0000
                                                R-squared         =     0.1255
                                                Root MSE          =      5.402
                              (Std. err. adjusted for 12 clusters in industry)
------------------------------------------------------------------------------
             |               Robust
        wage | Coefficient  std. err.      t    P>|t|     [95% conf. interval]
-------------+----------------------------------------------------------------
      tenure |      0.030      0.028     1.08   0.304       -0.032       0.093
     ttl_exp |      0.273      0.047     5.81   0.000        0.170       0.376
    collgrad |      3.255      0.397     8.21   0.000        2.382       4.129
       _cons |      3.420      0.285    11.98   0.000        2.792       4.048
------------------------------------------------------------------------------

. boottest tenure, svmat

Wild bootstrap-t, null imposed, 999 replications, Wald test, 
bootstrap clustering by industry, Rademacher weights:
  tenure
                           t(11) =     1.0772
                        Prob>|t| =     0.3053
95% confidence set for null hypothesis expression: [-.04135, .1302]

3.3 随机推断

随机推断 (Randomization Inference) 被认为是分析随机实验数据的一种较好的方法。尤其在样本量较少、具有聚类随机性、以及高杠杆点 (潜在离群值) 等情形下，其优势更加明显。

它的基本步骤是：(1) 保留样本的原始实验分配；(2) 根据原始实验组分配方式，重新随机分配实验组；(3) 用 “假的”，或者说 “安慰剂” 实验组作为一个附加项，重新估计原始回归方程；(4) 重复以上 (2)-(3) 过程；(5) 随机推断 $p$ 值是安慰剂处理效果大于估计实验效果的次数的比例。

* 命令安装
ssc install ritest, replace
* 命令语法
ritest resampvar exp_list [, options] : command

该语句随机排列 resampvar $n$ 次，并且每次都执行语句并收集 exp_list 中变量的实现值。options 中有分层 strata() 和聚类 cluster()，不设置分层即假设所有观测值在同一个单一分层中，不设置聚类即假设每个观测值为单一聚类。

. sysuse nlsw88.dta, clear
. reg wage $x

      Source |       SS           df       MS      Number of obs   =     2,227
-------------+----------------------------------   F(6, 2220)      =     66.17
       Model |  11229.3132         6   1871.5522   Prob > F        =    0.0000
    Residual |  62791.4974     2,220  28.2844583   R-squared       =    0.1517
-------------+----------------------------------   Adj R-squared   =    0.1494
       Total |  74020.8106     2,226  33.2528349   Root MSE        =    5.3183
------------------------------------------------------------------------------
        wage | Coefficient  Std. err.      t    P>|t|     [95% conf. interval]
-------------+----------------------------------------------------------------
       hours |      0.053      0.011     4.73   0.000        0.031       0.075
     ttl_exp |      0.252      0.030     8.26   0.000        0.192       0.311
      tenure |      0.022      0.025     0.88   0.381       -0.027       0.071
       south |     -1.528      0.229    -6.66   0.000       -1.978      -1.079
    collgrad |      3.167      0.267    11.86   0.000        2.643       3.691
     married |     -0.200      0.238    -0.84   0.401       -0.666       0.267
       _cons |      2.546      0.530     4.80   0.000        1.506       3.586
------------------------------------------------------------------------------

. ritest tenure _b[tenure]: reg wage $x

      Source |       SS           df       MS      Number of obs   =     2,227
-------------+----------------------------------   F(6, 2220)      =     66.17
       Model |  11229.3132         6   1871.5522   Prob > F        =    0.0000
    Residual |  62791.4974     2,220  28.2844583   R-squared       =    0.1517
-------------+----------------------------------   Adj R-squared   =    0.1494
       Total |  74020.8106     2,226  33.2528349   Root MSE        =    5.3183
------------------------------------------------------------------------------
        wage | Coefficient  Std. err.      t    P>|t|     [95% conf. interval]
-------------+----------------------------------------------------------------
       hours |      0.053      0.011     4.73   0.000        0.031       0.075
     ttl_exp |      0.252      0.030     8.26   0.000        0.192       0.311
      tenure |      0.022      0.025     0.88   0.381       -0.027       0.071
       south |     -1.528      0.229    -6.66   0.000       -1.978      -1.079
    collgrad |      3.167      0.267    11.86   0.000        2.643       3.691
     married |     -0.200      0.238    -0.84   0.401       -0.666       0.267
       _cons |      2.546      0.530     4.80   0.000        1.506       3.586
------------------------------------------------------------------------------

      Command: regress wage hours ttl_exp tenure south collgrad married
        _pm_1: _b[tenure]
  res. var(s):  tenure
   Resampling:  Permuting tenure
Clust. var(s):  __000001
     Clusters:  2246
Strata var(s):  none
       Strata:  1
------------------------------------------------------------------------------
T            |     T(obs)       c       n   p=c/n   SE(p) [95% Conf. Interval]
-------------+----------------------------------------------------------------
       _pm_1 |   .0219466      27     100  0.2700  0.0444  .1860664   .3680163
------------------------------------------------------------------------------
Note: Confidence interval is with respect to p=c/n.
Note: c = #{|T| >= |T(obs)|}

4. 相关推文

Note：产生如下推文列表的 Stata 命令为：
lianxh 标准误随机推断, m
安装最新版 lianxh 命令：
ssc install lianxh, replace

专题	嘉宾	直播/回看视频
⭐ 最新专题		文本分析、机器学习、效率专题、生存分析等
研究设计	连玉君	我的特斯拉-实证研究设计，-幻灯片-
面板模型	连玉君	动态面板模型，-幻灯片-
面板模型	连玉君	直击面板数据模型 [免费公开课，2小时]