Stata：聚类标准误的纠结
👇 连享会 · 推文导航 | www.lianxh.cn
🍎 Stata：Stata基础 | Stata绘图 | Stata程序 | Stata新命令 📘 论文：数据处理 | 结果输出 | 论文写作 | 数据分享 💹 计量：回归分析 | 交乘项-调节 | IV-GMM | 时间序列 | 面板数据 | 空间计量 | Probit-Logit | 分位数回归 ⛳ 专题：SFA-DEA | 生存分析 | 爬虫 | 机器学习 | 文本分析 🔃 因果：DID | RDD | 因果推断 | 合成控制法 | PSM-Matching 🔨 工具：工具软件 | Markdown | Python-R-Stata 🎧 课程：公开课-直播 | 计量专题 | 关于连享会
连享会课程 · 基于机器学习的因果推断方法
作者： 白锐（厦门大学）
邮箱： ruibai0608@126.com
编者按： 标准误在统计推断中发挥着至关重要的作用，直接影响着系数的显著性和置信区间，并最终影响到假设检验的结论。本文通过数据模拟分析，指出了使用聚类标准误的必要性。
目录
1. 背景介绍
1. 1 背景参考
1. 2 研究问题
2. 实操准备
2. 1 模拟设置
2. 2 主要程序
3. 主要结果
3. 1 不同情况下的结果对比
3. 2 基本结论
4. 参考资料
5. 相关推文
温馨提示： 文中链接在微信中无法生效。请点击底部「阅读原文」。或直接长按/扫描如下二维码，直达原文：
本文主要编译自下文：
Source： To cluster or not to cluster
1. 背景介绍
1. 1 背景参考
关于是否需要聚类这一问题，已经有一些有趣的文献进行了讨论，这些文献对何时应该使用聚集标准误给出了明确的建议(参见 Abadie 等人 2017)。特别地，如果存在序列自相关、截面自相关以及异方差等问题时，不使用聚类标准误就会导致得到错误的统计推断结果。不同于给出大量的数学证明，本文作者试图通过数据模拟分析来验证关于是否需要使用聚类标准误的问题。
1. 2 研究问题
本文的研究问题聚焦在：如果应该使用聚类标准误，却没有使用会产生怎样的后果？当然本文给出的模拟数据是有限的，所以结论可能有一定局限性。 但本文仍然有以下三点有益的启示：
给定数据假设，如何进行模拟分析 如何解释模拟结果 本应该使用聚类标准误却没有使用会带来怎样的后果
2. 实操准备
2. 1 模拟设置
讨论聚类标准误时，面板数据是很好的例子。面板数据包含两个维度的数据，截面维度和时间维度，所以可以控制不随时间变化的不可观测的个体异质性。 本文考虑了一个非常简单的面板结构，其中数据生成过程（Data Generating Process,简记 DGP）如下:
其中， 是个体不随时间改变的无法观测的影响因素， 是随个体与时间而改变的不可观测的扰动项。
为了简单起见，假设二者与解释变量 和 均不相关，所以可以将它们视为外生的。进一步，本文假设 随时间而变化，而 是不随时间变化的。至此，本文构建了一个平衡面板数据集，它既有横截面的维度(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）固定效应面板
因为前文假设 和 是独立于扰动项的外生变量，所以模型可以得到无偏估计，然而，标准误不一定是正确的。基于以上问题，本文展开了多次模拟分析并给出了结果。
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 给出了无偏估计量。 和 的系数几乎与总体系数 1 相同。常数项也是如此。
本文给出 2 个标准误的估计量：一个是基于模拟结果(_b),另一个是基于模拟的平均值(_se)。
结果显示，对于时变变量 的标准误的估计是正确的，因为 x (_se_x)的“平均”标准误接近_b_x 的标准误。然而， 和常数项的估计结果就没有那么好了。在这种情况下，我们无法从这个模型中做出正确的统计推断。以零假设 $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 回归，对 而言，拒绝原假设的平均值为 2%，但对 和常数项而言，拒绝率分别为是 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
基于上述结果， 和常数的平均标准误与基准模型模拟结果相似。此外，现在所有的零假设的拒绝率都和预期的一样，在 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
考虑了固定效应的结果与前三个稍有不同。因为模型中明确地包含个体固定效应，其不能估计不随时间而变的变量，所以无法再估计 的系数。
上述模拟分析结果显示，面板固定效应下的估计量的标准误比 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）能够得到 的估计，其标准误的大小与考虑聚类标准误差的 OLS 情形结果相似。
（3）能够得到常数项的假设检验结果。
接下来还有最后一种方法，被称为相关随机效应模型。前文假设不随时间变化的误差项与解释变量 和 均不相关，而在分析真实数据时，这是未知的，所以要应用相关随机效应模型，需要放宽这一假设并修改前文模拟程序：
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
最后一种数据模拟分析方法也表现良好。效果可与面板固定效应和随机效应情形下的估计量相媲美。当然，如果 与模型中的任何解释变量相关，它会比考虑随机效应情形的结果更好。
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
专题：回归分析 稳健性检验！稳健性检验！ Stata：聚类标准误的纠结 acreg：允许干扰项随意相关的稳健性标准误 Stata：聚类调整标准误笔记 Stata：聚类调整后的标准误-Cluster-SE 小样本下OLS估计的纠偏聚类标准误 专题：面板数据 Stata：面板聚类标准误-自动确定最优聚类层级和数量-xtregcluster wcbregress：面板聚类标准误 专题：倍分法DID 倍分法(DID)的标准误：不能忽略空间相关性 专题：内生性-因果推断 Stata新命令：konfound - 因果推断的稳健性检验 专题：时间序列 gcrobustvar：基于VAR的稳健性Granger因果检验
课程推荐：因果推断实用计量方法
主讲老师：邱嘉平教授
🍓 课程主页：https://gitee.com/lianxh/YGqjp
New！ Stata 搜索神器：
lianxh和songblGIF 动图介绍
搜： 推文、数据分享、期刊论文、重现代码 ……
👉 安装：
. ssc install lianxh
. ssc install songbl
👉 使用：
. lianxh DID 倍分法
. songbl all
🍏 关于我们
连享会 ( www.lianxh.cn，推文列表) 由中山大学连玉君老师团队创办，定期分享实证分析经验。 直通车： 👉【百度一下： 连享会】即可直达连享会主页。亦可进一步添加 「知乎」,「b 站」,「面板数据」,「公开课」 等关键词细化搜索。