聚类异质性：使用summclust进行统计推断
👇 连享会 · 推文导航 | www.lianxh.cn
🍎 Stata：Stata基础 | Stata绘图 | Stata程序 | Stata新命令 📘 论文：数据处理 | 结果输出 | 论文写作 | 数据分享 💹 计量：回归分析 | 交乘项-调节 | IV-GMM | 时间序列 | 面板数据 | 空间计量 | Probit-Logit | 分位数回归 ⛳ 专题：SFA-DEA | 生存分析 | 爬虫 | 机器学习 | 文本分析 🔃 因果：DID | RDD | 因果推断 | 合成控制法 | PSM-Matching 🔨 工具：工具软件 | Markdown | Python-R-Stata 🎧 课程：公开课-直播 | 计量专题 | 关于连享会
连享会 · 2022 面板数据因果推断专题
作者：王本丞 (中国人民大学)
邮箱：wangbencheng@ruc.edu.cn
编者按：本文主要参考自下文，特此致谢！
Source：MacKinnon J G, Nielsen M Ø, Webb M D. Leverage, influence, and the jackknife in clustered regression models: Reliable inference using summclust[J]. arXiv preprint arXiv:2205.03288, 2022. -PDF- -Link-
温馨提示： 文中链接在微信中无法生效。请点击底部「阅读原文」。或直接长按/扫描如下二维码，直达原文：
目录
1. 引言
2. 聚类稳健标准误
3. 聚类标准误异质性识别
3.1 高杠杆类 (high leverage)
3.2 偏杠杆类 (partial leverage)
3.3 强影响类 (influence)
3.4 小样本聚类标准误
4. Stata 实操
4.1 命令介绍
4.2 案例演示
5. 相关推文
1. 引言
在计量经济学的统计推断中，标准误扮演着重要角色。不过在当前的研究中，误差项满足独立同分布 (iid) 的假定并不能很好的反映真实情况，因此越来越多的研究者认为需要对标准误进行调整。其中聚类标准误 (cluster SE) 是一种最为常见的方法。
聚类稳健标准误放松了误差项满足独立同分布的假定，允许聚类内部个体间误差项存在相关性，但是聚类之间个体误差项不存在相关性。当上述条件得到满足，采用聚类调整的 OLS 估计值是无偏的。然而，聚类层级的选择会影响系数显著性和统计推断的结果，如何判断实证研究中聚类是否有效？
MacKinnon 等 (2022) 为聚类推断的有效性提供了检验方法。他们认为通过一系列特征统计量的分布特征可以分析聚类的有效性，并提供了 Stata 命令 summclust。接下来，本文将介绍 MacKinnon 等 (2022) 提出的聚类有效性分析思路以及 Stata 实操过程，以便加深读者对于相关内容的理解，提高实证研究的可靠性。
2. 聚类稳健标准误
聚类稳健标准误是当前使用最多的标准误调整方法。该方法在聚类内部允许个体之前存在任意形式的相关性和异质性，但是在聚类之间不存在任何相关性。给定采用聚类调整的回归方程：
其中，全部样本划分为 组，第 组中包含 个样本，总样本空间 。采用 OLS 得到估计系数的 “三明治” 标准形式：。在实际中，我们经常采用聚类稳健标准误 (CRVE) 如下所示：
其中，。Stata 中汇报的聚类稳健标准误如上式所示，当 时，聚类稳健标准误退化到常规的 标准误 (针对未知形式的异方差稳健标准误)。当 时，聚类稳健标准误具有良好的渐进性质 (Hansen 和 Lee，2019)。
一般意义上讲，估计量的渐进性质取决于聚类数量 和聚类的异质性程度。当不同聚类之间异质性程度越高，聚类稳健标准误的渐进性质越难得到满足。
3. 聚类标准误异质性识别
在样本数据的分析中，我们常用离群值 (outlier)、杠杆点 (leverage)、偏杠杆点 (partial leverage) 和强重要点 (influence) 进行异质性分析。为了更好的理解 MacKinnon 等 (2022) 的处理， 将上述概念展开如下：
outlier：因变量分布中远离样本空间中心的点。 high leverage：自变量分布中远离样本空间中心的点。常规测度方法是 “帽子” 统计量大于平均 “帽子” 统计量的 3 倍。其中 “帽子” 矩阵表示为 。 partial leverage：测度单个自变量对于每个观测值杠杆的影响，计算方式为 。其中 表示自变量， 表示观测样本点， 表示 回归到其余自变量上的残差。 influence：对模型有较大的影响的点，删除该点后显著的改变拟合回归方程。
上述特殊点的存在，表征了样本分布存在较强的异质性。在聚类调整估计中，我们的估计单位不再是样本个体，而是聚类使用的分组 。因此可以针对聚类的异质性进行考察，从而明确聚类估计的可信度。如果聚类空间 的分布具有较强的异质性，则当前使用的聚类估计可信度较差，应该采用更稳健的估计量进行调整。同时在小样本情况下应该采取 Jackknife 方法或 Bootstrap 方法。
3.1 高杠杆类 (high leverage)
从聚类的角度考察单个聚类与聚类空间中心的偏离程度。在估计 的时候，扔掉聚类 会导致第 的回归残差从 变为 ，其中 。直接汇报矩阵 并不方便，我们通过汇报矩阵 来测度杠杆，即：
高杠杆的聚类可以通过比较 和 来识别。如果对于某些聚类 有 ，则聚类 有可能是高杠杆聚类，在聚类标准误的分析中要谨慎。
3.2 偏杠杆类 (partial leverage)
聚类层面的偏杠杆计算如下：
其中， 是 聚类对应的 中的子向量，通过上式可以简便的计算任何感兴趣的估计系数。如果聚类 存在 ，则表明该聚类对于系数估计具有较高的偏杠杆影响。
3.3 强影响类 (influence)
强影响类聚类是指在删除该聚类后，标准误估计出现较大变化的类，删除聚类 之后的 OLS 估计为：
汇报特定参数估计的聚类 对应的 ，可以明确单个聚类对于参数估计的影响。如果某个聚类对应的 与其他聚类对应的参数 差异较大，则表明聚类 可能是强影响类。
通过上述三种表征聚类异质性的方法，MacKinnon 等 (2022) 提供了分析实证分析中聚类可靠性的思路。当 较小的时候，应当计算 以及少数参数对应的 。当 较大的时候应当汇报上述三类参数的分布情况 (描述性统计)，从而对于当前的聚类方法进行可靠性评估。
3.4 小样本聚类标准误
聚类稳健标准误 (CV_1) 在样本较小的情况下不具备良好的性质，因此可以采用 Jacknife 方法进行小样本聚类推断。实践中，采用 Jacknife 方法估计 的计算如下：
同时，聚类情况下的 CV_3 标准误计算如下：
4. Stata 实操
MacKinnon 等 (2022) 为上述聚类异质性分析提供了 Stata 命令 summclust，接下来，我们将对该命令进行介绍。
4.1 命令介绍
summclust 命令安装：
ssc install summclust, replace
summclust 命令语法：
summclust depvar, yvar(varname) xvar(varlist) cluster(varname) [ options]
其中，
depavr：回归中主要关注的解释变量，也是汇报不同聚类标准误的变量；yvar(varname)：被解释变量；xvar(varlist)：其他解释变量；cluster(varname)：聚类层级，例如行业层面、城市层面等；fevar(varname)：控制固定效应，例如城市固定效应、时间固定效应等。
options 包括：
svars汇报其他的描述性统计结果；gstar计算有效聚类的 和 数量；table汇报每个聚类的描述性统计结果；jackknife汇报 标准误等。
4.2 案例演示
我们研究的问题是女性婚姻的工资回报率，使用的数据是 nlswork.dta。其中，被解释变量是工资水平 lnwage，核心解释变量为是否结婚 msp，控制变量包括 union 和 race。此外还控制了 grade、age、birth_yr 等固定效应，并聚类到行业 (ind) 层面。
. webuse nlswork, clear
. reg ln_wage msp union race i.grade i.age i.birth_yr, cluster(ind)
Linear regression Number of obs = 19,130
F(11, 11) = .
Prob > F = .
R-squared = 0.2586
Root MSE = .40341
(Std. err. adjusted for 12 clusters in ind_code)
------------------------------------------------------------------------------
| Robust
ln_wage | Coefficient std. err. t P>|t| [95% conf. interval]
-------------+----------------------------------------------------------------
msp | -0.028 0.009 -2.96 0.013 -0.048 -0.007
union | 0.204 0.061 3.33 0.007 0.069 0.339
race | -0.086 0.016 -5.34 0.000 -0.122 -0.051
------------------------------------------------------------------------------
聚类到行业层面的回归结果如上所示，msp 的估计系数为 -0.028，聚类稳健标准误为 0.009， 值为 0.013，在 95% 水平上显著，即结婚会显著的降低个体的工资水平。下面我们使用 summclust 命令进行估计，分析聚类到行业层面是否存在异质性。
. summclust msp, yvar(ln_wage) xvar(union race) fevar(grade age birth_yr) cluster(ind)
SUMMCLUST - MacKinnon, Nielsen, and Webb
Cluster summary statistics for msp when clustered by ind_code.
There are 19130 observations within 12 ind_code clusters.
Regression Output
s.e. | Coeff Sd. Err. t-stat P value CI-lower CI-upper
-------+----------------------------------------------------------------
CV1 | -0.027515 0.009293 -2.9608 0.0130 -0.047969 -0.007061
CV3 | -0.027515 0.014064 -1.9564 0.0763 -0.058470 0.003440
------------------------------------------------------------------------
Cluster Variability
Statistic | Ng Leverage Partial L. beta no g
-----------+-----------------------------------------------------
min | 38.00 0.093321 0.001622 -0.033200
q1 | 159.00 0.672924 0.008649 -0.029275
median | 995.50 3.515491 0.056682 -0.027765
mean | 1594.17 5.416667 0.083333 -0.026920
q3 | 2335.50 7.731883 0.120933 -0.025975
max | 6335.00 20.289183 0.312995 -0.015835
-----------+-----------------------------------------------------
coefvar | 1.19 1.152965 1.141326 0.162898
回归中汇报了 msp 的不同聚类标准误，CV_1 表示常规的聚类稳健标准误，结果和上表中 reg 回归结果一致。在这种情况下，结婚对于工资收入具有负显著 (95% 水平) 的影响。但是，当我们使用 CV_3 标准误时，标准误估计为 0.014， 值为 0.0763，在 95% 水平上不显著，这表明使用聚类稳健标准误 (CV_1) 进行统计推断可能存在可信度的问题。下表汇报了聚类统计量的描述性统计结果：
Ng 表示聚类内部包含样本个数的分布，从聚类样本分布上来看，不同类内部存在较大的异质性。 Leverage 表示杠杆 ，从杠杆分布来看最大值和最小值之间存在非常大的差异 (217 倍)，表明不同聚类间处在较强的异质性，可能存在高杠杆类对于统计推断产生影响。 Partial L. 表示偏杠杆 ，从偏杠杆的分布来看极值之间的差异较大 (193 倍)。 beta no g 表示剔除聚类 之后的系数估计 ，从分布来看相对均衡，剔除某一聚类后并未改变估计系数的符号。
上述分析表明行业层面聚类存在类之间较强的异质性，对于统计推断产生了影响，为了增加实证研究的可信度，应当汇报 或 等标准误。
. summclust msp, yvar(ln_wage) xvar(union race) fevar(grade age birth_yr) ///
> absorb(ind) cluster(ind) table svars jack rho(0.5)
SUMMCLUST - MacKinnon, Nielsen, and Webb
Cluster summary statistics for msp when clustered by ind_code.
There are 19130 observations within 12 ind_code clusters.
Regression Output
s.e. | Coeff Sd. Err. t-stat P value CI-lower CI-upper
-------+----------------------------------------------------------------
CV1 | -0.020895 0.007084 -2.9494 0.0132 -0.036488 -0.005302
CV3 | -0.020895 0.007931 -2.6345 0.0232 -0.038352 -0.003438
CV3J | -0.020895 0.007921 -2.6381 0.0231 -0.038328 -0.003462
------------------------------------------------------------------------
Cluster Variability
Statistic | Ng Leverage Partial L. beta no g
-----------+-----------------------------------------------------
min | 38.00 0.087112 0.001561 -0.023382
q1 | 159.00 0.656606 0.008621 -0.022428
median | 995.50 3.442673 0.056073 -0.021258
mean | 1594.17 5.333333 0.083333 -0.020770
q3 | 2335.50 7.605927 0.121546 -0.020189
max | 6335.00 20.011074 0.312377 -0.015001
-----------+-----------------------------------------------------
coefvar | 1.19 1.155829 1.141658 0.120094
Effective Number of Clusters
-----------------------------
G*(0) = 5.468
-----------------------------
G*(rho) and G*(1) are not available.
There are fixed effects at the cluster or subcluster level.
Alternative Sample Means and Ratios to Arithmetic Mean
| Ng Leverage Partial L. beta no g
----------------+------------------------------------------------------
Harmonic Mean | 227.315 0.644626 0.010404 .
Harmonic Ratio | 0.143 0.120867 0.124846 .
Geometric Mean | 687.061 2.334997 0.035552 .
Geometric Ratio | 0.431 0.437812 0.426622 .
Quadratic Mean | 2413.502 7.954736 0.123456 0.020907
Quadratic Ratio | 1.514 1.491513 1.481475 -1.006589
-----------------------------------------------------------------------
Cluster by Cluster Statistics
ind_code | Ng Leverage Partial L. beta no g
-----------+-----------------------------------------------------
1 | 130 0.592119 0.005974 -0.021222
2 | 38 0.087112 0.001561 -0.021028
3 | 185 0.721093 0.009621 -0.021496
4 | 3747 13.514215 0.201461 -0.015001
5 | 1069 3.346972 0.060079 -0.023382
6 | 2912 10.224230 0.151294 -0.021295
7 | 1759 4.987623 0.091798 -0.019349
8 | 572 2.773414 0.027859 -0.023240
9 | 922 3.792664 0.052068 -0.021211
10 | 133 0.411110 0.007621 -0.021984
11 | 6335 20.011074 0.312377 -0.017157
12 | 1328 3.538374 0.078286 -0.022872
-----------------------------------------------------------------
总之，MacKinnon 等 (2022) 的建议在于应当对于聚类标准误的使用更加谨慎。在使用聚类标准误的时候，可以通过 summclust 命令汇报 、、 的分布情况，分析在不同类之间是否存在显著的异质性。如果存在的话当前的聚类标准可能存在统计推断的问题，为了提高研究的可信度，应当汇报 或 的标准误。
5. 相关推文
Note：产生如下推文列表的 Stata 命令为：
lianxh 标准误, m
安装最新版lianxh命令：
ssc install lianxh, replace
专题：Stata命令 Stata：双重机器学习-多维聚类标准误的估计方法-crhdreg 专题：回归分析 Stata：标准误！标准误！ Stata：聚类标准误的纠结 acreg：允许干扰项随意相关的稳健性标准误 Stata：聚类调整标准误笔记 Stata：聚类调整后的标准误-Cluster-SE 小样本下OLS估计的纠偏聚类标准误 专题：面板数据 Stata：面板聚类标准误-自动确定最优聚类层级和数量-xtregcluster wcbregress：面板聚类标准误 专题：倍分法DID 倍分法(DID)的标准误：不能忽略空间相关性
课程推荐：因果推断实用计量方法
主讲老师：邱嘉平教授
🍓 课程主页：https://gitee.com/lianxh/YGqjp
New！ Stata 搜索神器：
lianxh和songblGIF 动图介绍
搜： 推文、数据分享、期刊论文、重现代码 ……
👉 安装：
. ssc install lianxh
. ssc install songbl
👉 使用：
. lianxh DID 倍分法
. songbl all
🍏 关于我们
连享会 ( www.lianxh.cn，推文列表) 由中山大学连玉君老师团队创办，定期分享实证分析经验。 直通车： 👉【百度一下： 连享会】即可直达连享会主页。亦可进一步添加 「知乎」,「b 站」,「面板数据」,「公开课」 等关键词细化搜索。