Stata:聚类标准误的纠结

发布时间:2021-11-19 阅读 251

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 系列推文:

PDF下载 - 推文合集

作者: 白锐(厦门大学)
邮箱: ruibai0608@126.com

编者按: 标准误在统计推断中发挥着至关重要的作用,直接影响着系数的显著性和置信区间,并最终影响到假设检验的结论。简单的说,聚类标准差是当误差项之间存在相关性时调整标准误的一种有效的方法。


目录


本文主要编译自下文:
Source: To cluster or not to cluster

1.背景介绍

1.1 聚类还是不聚类,这是一个问题

在开始讨论前,有必要给出一些有趣的文献,它们对什么时候应该和不应该使用聚集标准误给出了明确的建议(参见Abadie等人2017)。事实上,就在最近(撰写本文时,已是2021年3月12日),Wooldridge教授就聚类的使用提供了一些强有力的建议,建议至少应该给出稳健标准误。这是大多数计量经济学中都可以使用的选项(包括我最喜欢的…Stata)。相比给出大量的数学证明,在这里我会针对具体问题给出通过模拟得到的证据。

1.2 如果本应该聚类,但没有聚类会怎样

此处给出的模拟数据是有限的,一些结论可能是有局限性的。 无论如何,我相信你可以从中学到以下三点。

a. 给定数据假设,如何通过模拟来分析。 b. 如何解释模拟结果。 c. 本应该聚类的却没有聚类会带来怎样的后果

2.实操准备

2.1 模拟设置

当考虑聚类标准错误时,脑海中最重要的例子是面板数据。 面板数据是微观经济学家和宏观经济学家能够接触到的一种特殊类型的数据。面板数据最重要的特点是既有横截面的维度,又有时间的维度,所以如果个体未观察到的特征信息在时间上对每个个体是固定的,就可以控制这一效应。 当然,这种结构并非面板数据所独有。如果你收集家庭成员和他们的家庭的数据,或者县和州的数据,也可以有类似的结论。 现在回到我们的问题,让我们考虑一个非常简单的面板结构,其中数据生成过程(Data Generating Process,简记DGP)定义为:

其中ei是个体无法观测但在时间上固定的因素,vit是随个体与时间而改变的扰动项 为了简单起见,假设这两个扰动项与解释变量xitzit无关,所以可以将它们视为外生的。我们可以假设xit随时间而变化,但zit是固定的。为了简单起见,我们构建一个平衡面板数据集,因此所有N个个体在T周期内都能被观察到。在这种结构下,我们可以进行如下模拟:

 . * Consider 100 individuals:
. set obs 100
number of observations (_N) was 0, now 100
. * For which I'll create an identifier.
. gen id=_n
. * And I will assume that each individual is observed for 10 periods
. expand 10
(900 observations created)
. * Lets assume that X and Z are independent from each other.
. * but that both are "normal"
. gen x = rnormal()
. gen z = rnormal()
. * And Z is fixed for each individual
. bysort id:replace z=z[1]
(900 real changes made)
. * I use the same code for the unobserved factors. 
. gen e = rnormal()
. gen v = rnormal()
. * And Z is fixed for each individual
. bysort id:replace e=e[1]
(900 real changes made)
. * Finally, we construct our dependent variable Y
. gen y=1+x+z+e+v

对于这个简单的DGP,我们将对比5种不同估算方法的结果:

(1)简单的OLS (2)基于稳健标准误差的OLS (3)基于聚类标准误的OLS (4)随机效应面板 (5)固定效应面板

你可能会注意到第1、2种和第3种并不明确地控制个体固定效应。而第5种,尝试更明确地控制这种效果。第4种还会构建标准误差来控制未观察到的影响。第3种并没有做更多的假设。

因为xitzit是独立于扰动项的外生变量,所以我们可以得到无偏估计。然而,并不是所有的估计都能做得到。那我们怎么考虑这个问题呢?我们怎么知道系数是无偏的,但标准误可能是有误的?基于以上问题,我们进行了多次模拟并给出了结果

2.2 主要程序

首先,创建一个程序来模拟数据(如前所述),并估计前面的列表中的模型。其次,需要大量重复这个过程(通常是数万次),并收集结果进行评估。

 . . * I'm writing this as an "eclass" estimator, so it is easy to collect the estimation results.
. program panel_re, eclass
  1. * a new line. "clear" to start from a clean dataset each time the program is called:
.         clear
  2.         set obs 100
  3.         gen id=_n
  4.         expand 10
  5.         gen x = rnormal()
  6.         gen z = rnormal()
  7.         sort id, stable
  8.         by id:replace z=z[1]
  9.         gen e = rnormal()
 10.         gen v = rnormal()
 11.         by id:replace e=e[1]
 12.         gen y=1+x+z+e+v
* here I estimate the model
* I will use two "locals": 1 and 2. 
* This will be used as place holders for the model
* and options for the estimation.
* I'll also set this data as panel data (it will be needed for some cases)
 13.         xtset id
 14.         `1' y x z, `2'
 15.         ** and that is it! 
. end

基本程序如上。例如,如果想运行简单的OLS模型并考虑聚类标准误的模拟,可以按照以下操作。

. *panel_re reg  cluster(id)
. *          ^       ^
. *          |       |
. *         '1'     '2'
. *Here reg will be argument 1
. *and cluster(id) argument 2
. panel_re reg  cluster(id)
number of observations (_N) was 0, now 100
(900 observations created)
(900 real changes made)
(900 real changes made)
       panel variable:  id (balanced)

Linear regression                               Number of obs     =      1,000
                                                F(2, 99)          =     355.71
                                                Prob > F          =     0.0000
                                                R-squared         =     0.4973
                                                Root MSE          =     1.4754

                                   (Std. Err. adjusted for 100 clusters in id)
------------------------------------------------------------------------------
             |               Robust
           y |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
           x |   1.033408   .0409937    25.21   0.000     .9520676    1.114748
           z |   .9983844   .0895856    11.14   0.000     .8206271    1.176142
       _cons |   .9965414   .1162236     8.57   0.000     .7659286    1.227154
------------------------------------------------------------------------------

3. 主要结果

3.1 不同情况下的结果对比

根据设置好的程序,我们可以执行几次模拟,看看如果在估计这个模型时使用或忽略聚类会发生什么。虽然通常建议运行10000次迭代以正确地评估结果,但我将只使用100次迭代,因为这足以得到结果。当然,欢迎你使用更多的重复,看看会发生什么。 对于模拟,我们将使用Stata命令simulate,这使分析变得简单。我将使用seed(101)使结果可重复。 让我们从只考虑简单标准误的OLS回归开始。每次我们会收集点估计和标准误结果:

. simulate _b _se, reps(100) seed(101):panel_re reg  

      command:  panel_re reg

Simulations (100)
----+--- 1 ---+--- 2 ---+--- 3 ---+--- 4 ---+--- 5 
..................................................    50
..................................................   100

. sum, sep(3)

    Variable |        Obs        Mean    Std. Dev.       Min        Max
-------------+---------------------------------------------------------
        _b_x |        100    1.000355    .0381124   .8919857   1.085083
        _b_z |        100    .9876579    .0994754   .7406453     1.2002
     _b_cons |        100    .9984308    .1107775   .7436068   1.238203
-------------+---------------------------------------------------------
       _se_x |        100    .0444944    .0023557   .0397996   .0505446
       _se_z |        100    .0454697    .0039679   .0382312   .0611115
    _se_cons |        100    .0446636    .0022492   .0396397   .0515558

首先,正如预期,OLS提供了无偏估计量。xz的系数几乎与1的总体系数相同。常数项也是一样。

我们现在有两个对于标准误的估计量。一个是基于模拟结果(__b)_,另一个是基于模拟的平均值(_se)。

结果显示,对于时变变量x的标准误的估计是正确的,因为x (__se_x)的“平均”标准误接近_b_x的标准偏差。然而,z和常数项的对比结果就没有那么好了。在这种情况下,我们无法从这个模型中做出正确的统计推断。例如,拒绝零假设 H0:βz=1 的频率“太频繁”,如下所示。

. gen is_bx_1 = !inrange(1,_b_x-_se_x*1.96,_b_x+_se_x*1.96)

. gen is_bz_1 = !inrange(1,_b_z-_se_z*1.96,_b_z+_se_z*1.96)

. gen is_bc_1 = !inrange(1,_b_c-_se_c*1.96,_b_c+_se_c*1.96)

. sum is*

    Variable |        Obs        Mean    Std. Dev.       Min        Max
-------------+---------------------------------------------------------
     is_bx_1 |        100         .02    .1407053          0          1
     is_bz_1 |        100         .43    .4975699          0          1
     is_bc_1 |        100         .42     .496045          0          1

因此,使用OLS回归,对x而言,在2%的显著性水平上拒绝原假设(正如你所期望的),但对z和常数项而言,这个结果是43%和42%。

也许如果我们使用稳健标准误差呢?这可能不会有任何影响。虽然这个模型已经假设误差项是同方差的,但尝试一下也无妨:

. simulate _b _se, reps(100) seed(101):panel_re reg robust 

      command:  panel_re reg robust

Simulations (100)
----+--- 1 ---+--- 2 ---+--- 3 ---+--- 4 ---+--- 5 
..................................................    50
..................................................   100

. gen is_bx_1 = !inrange(1,_b_x-_se_x*1.96,_b_x+_se_x*1.96)
. gen is_bz_1 = !inrange(1,_b_z-_se_z*1.96,_b_z+_se_z*1.96)
. gen is_bc_1 = !inrange(1,_b_c-_se_c*1.96,_b_c+_se_c*1.96)
. sum, sep(3)

    Variable |        Obs        Mean    Std. Dev.       Min        Max
-------------+---------------------------------------------------------
        _b_x |        100    1.000355    .0381124   .8919857   1.085083
        _b_z |        100    .9876579    .0994754   .7406453     1.2002
     _b_cons |        100    .9984308    .1107775   .7436068   1.238203
-------------+---------------------------------------------------------
       _se_x |        100    .0444072    .0027773   .0391417    .052421
       _se_z |        100    .0453708    .0049049   .0347456    .064886
    _se_cons |        100     .044696     .002205   .0396771   .0515227
-------------+---------------------------------------------------------
     is_bx_1 |        100         .02    .1407053          0          1
     is_bz_1 |        100          .4     .492366          0          1
     is_bc_1 |        100         .42     .496045          0          1

情况和之前一样。

下一个是考虑聚类标准误差的回归。

. . simulate _b _se, reps(100) seed(101):panel_re reg cluster(id) 

      command:  panel_re reg cluster(id)

Simulations (100)
----+--- 1 ---+--- 2 ---+--- 3 ---+--- 4 ---+--- 5 
..................................................    50
..................................................   100

. gen is_bx_1 = !inrange(1,_b_x-_se_x*1.96,_b_x+_se_x*1.96)
. gen is_bz_1 = !inrange(1,_b_z-_se_z*1.96,_b_z+_se_z*1.96)
. gen is_bc_1 = !inrange(1,_b_c-_se_c*1.96,_b_c+_se_c*1.96)
. sum, sep(3)

    Variable |        Obs        Mean    Std. Dev.       Min        Max
-------------+---------------------------------------------------------
        _b_x |        100    1.000355    .0381124   .8919857   1.085083
        _b_z |        100    .9876579    .0994754   .7406453     1.2002
     _b_cons |        100    .9984308    .1107775   .7436068   1.238203
-------------+---------------------------------------------------------
       _se_x |        100    .0441533    .0051071   .0332627   .0670466
       _se_z |        100    .1049323    .0155818   .0745955   .1655501
    _se_cons |        100    .1045397    .0085235   .0841128   .1316457
-------------+---------------------------------------------------------
     is_bx_1 |        100         .04    .1969464          0          1
     is_bz_1 |        100         .05    .2190429          0          1
     is_bc_1 |        100         .06    .2386833          0          1

基于上述结果,z和常数的平均标准误与基准模型模拟结果相似。此外,现在所有的零假设的拒绝率都和预期的一样,在5%左右。

如果我们明确地控制这些个体效应会怎样呢?我们可以用FE面板来做

. simulate _b _se, reps(100) seed(101):panel_re xtreg fe 

      command:  panel_re xtreg fe

Simulations (100)
----+--- 1 ---+--- 2 ---+--- 3 ---+--- 4 ---+--- 5 
..................................................    50
..................................................   100

. gen is_bx_1 = !inrange(1,_b_x-_se_x*1.96,_b_x+_se_x*1.96)
. gen is_bc_1 = !inrange(1,_b_c-_se_c*1.96,_b_c+_se_c*1.96)
. sum, sep(3)

    Variable |        Obs        Mean    Std. Dev.       Min        Max
-------------+---------------------------------------------------------
        _b_x |        100    .9956048     .029711   .9376439   1.062583
      _sim_2 |        100           0           0          0          0
     _b_cons |        100    .9955147    .1434849   .6446086   1.326498
-------------+---------------------------------------------------------
       _se_x |        100    .0333269    .0011717   .0301809   .0364448
      _sim_5 |        100           0           0          0          0
    _se_cons |        100    .0316341     .000825   .0287964   .0339852
-------------+---------------------------------------------------------
     is_bx_1 |        100           0           0          0          0
     is_bc_1 |        100         .68    .4688262          0          1

所以事情现在看起来不一样了。因为模型中明确地包含个体固定效应,所以我们不再能z的系数。

有点奇怪的是,我们从不拒绝 H0:βx=1,虽然这可能是由于有限的重复次数,但如果我们重复更多次(比如1000次),也观察到这个值小于5%。同样重要的是,要注意到面板FE估计量的标准误比OLS聚集标准误差要小一些。

这可以看作是一种估计βx效率上的增益。

不幸的是,天下没有免费的午餐!虽然我们似乎可以更有效地估计βx,但我们不能再估计常数项了。事实上,基于模拟的标准误增加,零假设的拒绝率也会增加。

最后,我们可以使用面板随机效应模型:

. simulate _b _se, reps(100) seed(101):panel_re xtreg re 

      command:  panel_re xtreg re

Simulations (100)
----+--- 1 ---+--- 2 ---+--- 3 ---+--- 4 ---+--- 5 
..................................................    50
..................................................   100

. gen is_bx_1 = !inrange(1,_b_x-_se_x*1.96,_b_x+_se_x*1.96)
. gen is_bz_1 = !inrange(1,_b_z-_se_z*1.96,_b_z+_se_z*1.96)
. gen is_bc_1 = !inrange(1,_b_c-_se_c*1.96,_b_c+_se_c*1.96)
. sum, sep(3)

    Variable |        Obs        Mean    Std. Dev.       Min        Max
-------------+---------------------------------------------------------
        _b_x |        100    .9961247    .0289081   .9367011   1.064042
        _b_z |        100    .9877817    .0993992   .7413588   1.199385
     _b_cons |        100    .9984448    .1109847   .7435786   1.240776
-------------+---------------------------------------------------------
       _se_x |        100    .0331599     .001161    .030052   .0362103
       _se_z |        100    .1067879    .0120344   .0871601   .1579433
    _se_cons |        100    .1048667    .0087699   .0846491   .1327749
-------------+---------------------------------------------------------
     is_bx_1 |        100           0           0          0          0
     is_bz_1 |        100         .03    .1714466          0          1
     is_bc_1 |        100         .06    .2386833          0          1

在DGP的严格假设下,这个模型似乎是所有策略中最好的。

(1)其标准误与基于FE估计的一样小。

(2)我们仍然能够得到z的估计,标准误的大小与考虑聚类标准误差的OLS相似。

(3)我们还可以得到常数项的显著性水平

接下来还有第六种方法,被称为相关随机效应模型。基于前文假设,其实这里不需要这个模型,因为我们知道已经设定时间固定误差与其他特征不相关。但在分析真实数据时,这是未知的。

要应用这种方法,需要稍微修改模拟程序:

. program panel_cre, eclass
  1.         clear
  2.         set obs 100
  3.         gen id=_n
  4.         expand 10
  5.         gen x = rnormal()
  6.         gen z = rnormal()
  7.         sort id, stable
  8.         by id:replace z=z[1]
  9.         gen e = rnormal()
 10.         gen v = rnormal()
 11.         by id:replace e=e[1]
 12.         gen y=1+x+z+e+v
 13.         xtset id
 14.         by id:egen m_x=mean(x)  // <-- Here we estimate the within person mean characteristic
 15.         xtreg y x z m_x, re   // and control for it in the RE effects model.
 16. end

. simulate _b _se, reps(100) seed(101):panel_cre 

      command:  panel_cre

Simulations (100)
----+--- 1 ---+--- 2 ---+--- 3 ---+--- 4 ---+--- 5 
..................................................    50
..................................................   100

. gen is_bx_1 = !inrange(1,_b_x-_se_x*1.96,_b_x+_se_x*1.96)

. gen is_bz_1 = !inrange(1,_b_z-_se_z*1.96,_b_z+_se_z*1.96)

. gen is_bc_1 = !inrange(1,_b_c-_se_c*1.96,_b_c+_se_c*1.96)

. sum, sep(4)

    Variable |        Obs        Mean    Std. Dev.       Min        Max
-------------+---------------------------------------------------------
        _b_x |        100    .9956048     .029711   .9376439   1.062583
        _b_z |        100    .9863926    .1006587   .7307997   1.207344
      _b_m_x |        100     .051984    .3436732  -.6213726   1.080986
     _b_cons |        100    .9982991    .1093144   .7438877   1.214026
-------------+---------------------------------------------------------
       _se_x |        100    .0333269    .0011717   .0301809   .0364448
       _se_z |        100    .1072043    .0120987    .087363   .1584293
     _se_m_x |        100    .3380803    .0343612   .2572779   .4165666
    _se_cons |        100    .1055241    .0090165   .0849954   .1330222
-------------+---------------------------------------------------------
     is_bx_1 |        100           0           0          0          0
     is_bz_1 |        100         .04    .1969464          0          1
     is_bc_1 |        100         .06    .2386833          0          1

最后一种评估策略也表现良好。效果可与面板FE和随机FE估计量相媲美。然而,需要注意的是,只有在使用随机效应模型或使用聚集标准误差(当使用OLS时)时,它才会表现良好。否则,结果将和简单的OLS方法一样糟糕。当然,如果e与模型中的任何解释变量相关,它将比RE更好。

3.2 基本结论

聚类还是不聚类?

我从这个问题开始,是为了在你怀疑某个特定维度中存在无法观测到的因素时,促使你决定使用聚集标准误差来估计模型。

从实证的方法,我已经表明,如果使用面板数据,有个别未观察到的影响是随时间固定的,不使用聚类标准误将是一个错误。即使是稳健标准误的,也不会起多大作用。

在这种情况下,通过控制固定效应(例如使用面板FE)或随机效应或至少纠正存在的标准误(聚集标准误差),会更好。

这些讨论并不能概括所有情况。不过总的来说,它们是一个很好的起点。人们可以构建不同的模拟设置,以测试特定假设的变化会在多大程度上影响结果,以及某一估计方法会如何处理这个问题。

如果你对何时使用(或不使用)聚集标准误的的讨论有更进一步的兴趣,可以参考Abadie等人(2017)的工作。

4.参考资料

1、Alberto Abadie & Susan Athey & Guido W. Imbens & Jeffrey Wooldridge, 2017. "When Should You Adjust Standard Errors for Clustering?," NBER Working Papers 24003, National Bureau of Economic Research, Inc.-PDF-

2、Abowd, J. M., F. Kramarz, and S. Woodcock. 2008. Econometric analyses of linked employer-employee data. In The Econometrics of Panel Data: Fundamentals and Recent Developments in Theory and Practice, ed. L. M´aty´as and P. Sevestre, 3rd ed., 727–760. Berlin: Springer.-PDF-

5. 相关推文

Note:产生如下推文列表的 Stata 命令为:
lianxh 聚类
安装最新版 lianxh 命令:
ssc install lianxh, replace

相关课程

免费公开课

最新课程-直播课

专题 嘉宾 直播/回看视频
最新专题 文本分析、机器学习、效率专题、生存分析等
研究设计 连玉君 我的特斯拉-实证研究设计-幻灯片-
面板模型 连玉君 动态面板模型-幻灯片-
面板模型 连玉君 直击面板数据模型 [免费公开课,2小时]
  • Note: 部分课程的资料,PPT 等可以前往 连享会-直播课 主页查看,下载。

课程主页

课程主页

关于我们

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

连享会小程序:扫一扫,看推文,看视频……

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

✏ 连享会-常见问题解答:
https://gitee.com/lianxh/Course/wikis

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