Stata：聚类标准误的纠结

发布时间：2021-11-19 阅读 7315

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下载 - 推文合集

作者： 白锐（厦门大学）
邮箱： ruibai0608@126.com

编者按： 标准误在统计推断中发挥着至关重要的作用，直接影响着系数的显著性和置信区间，并最终影响到假设检验的结论。本文通过数据模拟分析，指出了使用聚类标准误的必要性。

本文主要编译自下文：

Source： To cluster or not to cluster

1. 背景介绍

1. 1 背景参考

关于是否需要聚类这一问题，已经有一些有趣的文献进行了讨论，这些文献对何时应该使用聚集标准误给出了明确的建议(参见 Abadie 等人 2017)。特别地，如果存在序列自相关、截面自相关以及异方差等问题时，不使用聚类标准误就会导致得到错误的统计推断结果。不同于给出大量的数学证明，本文作者试图通过数据模拟分析来验证关于是否需要使用聚类标准误的问题。

1. 2 研究问题

本文的研究问题聚焦在：如果应该使用聚类标准误，却没有使用会产生怎样的后果？当然本文给出的模拟数据是有限的，所以结论可能有一定局限性。但本文仍然有以下三点有益的启示：

给定数据假设，如何进行模拟分析
如何解释模拟结果
本应该使用聚类标准误却没有使用会带来怎样的后果

2. 实操准备

2. 1 模拟设置

讨论聚类标准误时，面板数据是很好的例子。面板数据包含两个维度的数据，截面维度和时间维度，所以可以控制不随时间变化的不可观测的个体异质性。本文考虑了一个非常简单的面板结构，其中数据生成过程（Data Generating Process,简记 DGP）如下:

y_{i t} = a_{0} + a_{1} x_{i t} + a_{2} z_{i t} + e_{i} + v_{i t} (1)

其中， $e_{i}$ 是个体不随时间改变的无法观测的影响因素， $v_{i t}$ 是随个体与时间而改变的不可观测的扰动项。

为了简单起见，假设二者与解释变量 $x_{i t}$ 和 $z_{i t}$ 均不相关，所以可以将它们视为外生的。进一步，本文假设 $x_{i t}$ 随时间而变化，而 $z_{i t}$ 是不随时间变化的。至此，本文构建了一个平衡面板数据集，它既有横截面的维度(N 个个体)，又有时间维度(T 个时期)。基于上述假定，本文的模拟设置如下:

. * 考虑有100个个体:
. set obs 100
number of observations (_N) was 0, now 100
. * 生成id变量
. gen id=_n
. * 假设观测10个时期
. expand 10
(900 observations created)
. * 假设解释变量X和Z相互独立
. * 假设解释变量为正态分布
. gen x = rnormal()
. gen z = rnormal()
. * 假设Z不随时间变化
. bysort id:replace z=z[1]
(900 real changes made)
. * 同理构造不可观测的影响因素
. gen e = rnormal()
. gen v = rnormal()
. * 假设Z不随时间变化
. bysort id:replace e=e[1]
(900 real changes made)
. * 构造因变量Y
. gen y=1+x+z+e+v

基于这一简单的 DGP，本文接下来将对比 5 种不同估计方法的结果，包括:

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

因为前文假设 $x_{i t}$ 和 $z_{i t}$ 是独立于扰动项的外生变量，所以模型可以得到无偏估计，然而，标准误不一定是正确的。基于以上问题，本文展开了多次模拟分析并给出了结果。

2. 2 主要程序

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

. program panel_re, 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
* 对模型进行估计
* 设置2个局部宏: 1 and 2.
* 在特定情况下，需将数据设置为面板数据
 13.         xtset id
 14.         `1' y x z, `2'
 15.         ** 设置完成
. end

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

*panel_re reg  cluster(id)
*          ^       ^
*          |       |
*         '1'     '2'
*其中reg对应'1'，cluster(id)对应'2'

. panel_re reg  cluster(id)

       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 给出了无偏估计量。 $x$ 和 $z$ 的系数几乎与总体系数 1 相同。常数项也是如此。

本文给出 2 个标准误的估计量：一个是基于模拟结果(_b),另一个是基于模拟的平均值(_se)。

结果显示，对于时变变量 $x$ 的标准误的估计是正确的，因为 x (_se_x)的“平均”标准误接近_b_x 的标准误。然而， $z$ 和常数项的估计结果就没有那么好了。在这种情况下，我们无法从这个模型中做出正确的统计推断。以零假设 $H{0}:\beta_{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%，远高于 5%的理论水平。

进一步，在其他条件不变的情况下，本文给出将标准误替换为稳健标准误差后的结果。

. 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%左右。

接着，本文明确地控制个体效应，考虑面板固定效应。

. 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$ 的系数。

上述模拟分析结果显示，面板固定效应下的估计量的标准误比 OLS 聚类标准误要小一些。

随后，本文考虑使用面板随机效应模型：

. 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）其标准误与基于固定效应模型估计的一样小。

（2）能够得到 $z$ 的估计，其标准误的大小与考虑聚类标准误差的 OLS 情形结果相似。

（3）能够得到常数项的假设检验结果。

接下来还有最后一种方法，被称为相关随机效应模型。前文假设不随时间变化的误差项与解释变量 $x_{i t}$ 和 $z_{i t}$ 均不相关，而在分析真实数据时，这是未知的，所以要应用相关随机效应模型，需要放宽这一假设并修改前文模拟程序：

program define panel_cre, eclass
        clear
        set obs 100
        gen id=_n
        expand 10
        gen x = rnormal()
        gen z = rnormal()
        sort id, stable
        by id:replace z=z[1]
        gen e = rnormal()
        gen v = rnormal()
        by id:replace e=e[1]
        gen y=1+x+z+e+v
        xtset id
        by id:egen m_x=mean(x) // 此处估计了组内均值
        xtreg y x z m_x, re    // 考虑随机效应
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

最后一种数据模拟分析方法也表现良好。效果可与面板固定效应和随机效应情形下的估计量相媲美。当然，如果 $e$ 与模型中的任何解释变量相关，它会比考虑随机效应情形的结果更好。

3. 2 基本结论

是否需要使用聚类标准误?回到本文最初提出的问题，我们得到的答案是最好使用聚集标准误来估计模型。

从实证的方法，本文已经表明，如果使用面板数据，考虑存在个体不随时间改变的无法观测的影响因素，结合数据生成过程来看，使用普通标准误将得到错误结果，甚至使用稳健标准误也只能得到错误结果。如果不使用聚类标准误，就会得到上述错误的标准误结果，从而得到错误的统计推断结果，而通过使用聚类标准误，可以修正得到的标准误结果。

上述讨论仍有未尽之处，并不能概括所有情况。不过本文讨论的范式值得借鉴，可参照本文构建不同的模拟设置，以测试特定假设及特定估计方法会如何影响结果。进一步的研究可以参考 Abadie 等人(2017)的工作。

4. 参考资料

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-
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 标准误稳健性, m
安装最新版 lianxh 命令：
ssc install lianxh, replace

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

Stata：聚类标准误的纠结

1. 背景介绍

1. 1 背景参考

1. 2 研究问题

2. 实操准备

2. 1 模拟设置

2. 2 主要程序

3. 主要结果

3. 1 不同情况下的结果对比

3. 2 基本结论

4. 参考资料

5. 相关推文

相关课程

免费公开课

最新课程-直播课

⛳ 课程主页

⛳ 课程主页

关于我们

热门资讯