聚类异质性：使用summclust进行统计推断

发布时间：2022-10-18 阅读 1833

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

作者：王本丞 (中国人民大学)
邮箱：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. 引言

在计量经济学的统计推断中，标准误扮演着重要角色。不过在当前的研究中，误差项满足独立同分布 (iid) 的假定并不能很好的反映真实情况，因此越来越多的研究者认为需要对标准误进行调整。其中聚类标准误 (cluster SE) 是一种最为常见的方法。

聚类稳健标准误放松了误差项满足独立同分布的假定，允许聚类内部个体间误差项存在相关性，但是聚类之间个体误差项不存在相关性。当上述条件得到满足，采用聚类调整的 OLS 估计值是无偏的。然而，聚类层级的选择会影响系数显著性和统计推断的结果，如何判断实证研究中聚类是否有效？

MacKinnon 等 (2022) 为聚类推断的有效性提供了检验方法。他们认为通过一系列特征统计量的分布特征可以分析聚类的有效性，并提供了 Stata 命令 summclust。接下来，本文将介绍 MacKinnon 等 (2022) 提出的聚类有效性分析思路以及 Stata 实操过程，以便加深读者对于相关内容的理解，提高实证研究的可靠性。

2. 聚类稳健标准误

聚类稳健标准误是当前使用最多的标准误调整方法。该方法在聚类内部允许个体之前存在任意形式的相关性和异质性，但是在聚类之间不存在任何相关性。给定采用聚类调整的回归方程：

y_{g} = X_{g} β + u_{g}, g = 1, \cdot \cdot \cdot \cdot, G

其中，全部样本划分为 $G$ 组，第 $g$ 组中包含 $N_{g}$ 个样本，总样本空间 $N = \sum_{g = 1}^{G} N_{g}$ 。采用 OLS 得到估计系数的 “三明治” 标准形式： $β = (X^{'} X)^{- 1} (\sum_{g = 1}^{G} \sum_{g}) (X^{'} X)^{- 1}$ 。在实际中，我们经常采用聚类稳健标准误 (CRVE) 如下所示：

{CV}_{1} : \frac{G (N - 1)}{(G - 1) (N - k)} (X^{'} X)^{- 1} (\sum_{g = 1}^{G} {\hat{s}}_{g} {\hat{s}}_{g}^{'}) (X^{'} X)^{- 1}

其中， ${\hat{s}}_{g} = X_{g}^{'} {\hat{u}}_{g}$ 。Stata 中汇报的聚类稳健标准误如上式所示，当 $G = N$ 时，聚类稳健标准误退化到常规的 $H C_{1}$ 标准误 (针对未知形式的异方差稳健标准误)。当 $G \to \infty$ 时，聚类稳健标准误具有良好的渐进性质 (Hansen 和 Lee，2019)。

一般意义上讲，估计量的渐进性质取决于聚类数量 $G$ 和聚类的异质性程度。当不同聚类之间异质性程度越高，聚类稳健标准误的渐进性质越难得到满足。

3. 聚类标准误异质性识别

在样本数据的分析中，我们常用离群值 (outlier)、杠杆点 (leverage)、偏杠杆点 (partial leverage) 和强重要点 (influence) 进行异质性分析。为了更好的理解 MacKinnon 等 (2022) 的处理，将上述概念展开如下：

outlier：因变量分布中远离样本空间中心的点。
high leverage：自变量分布中远离样本空间中心的点。常规测度方法是 “帽子” 统计量大于平均 “帽子” 统计量的 3 倍。其中 “帽子” 矩阵表示为 $H = X (X^{'} X)^{- 1} X^{'}$ 。
partial leverage：测度单个自变量对于每个观测值杠杆的影响，计算方式为 $(P L_{j})_{i} = \frac{(X_{j \cdot [j]})_{i}^{2}}{\sum_{k}^{n} (X_{j \cdot [j]})_{k}^{2})}$ 。其中 $j$ 表示自变量， $i$ 表示观测样本点， $X_{j \cdot [j]}$ 表示 $X_{j}$ 回归到其余自变量上的残差。
influence：对模型有较大的影响的点，删除该点后显著的改变拟合回归方程。

上述特殊点的存在，表征了样本分布存在较强的异质性。在聚类调整估计中，我们的估计单位不再是样本个体，而是聚类使用的分组 $g$ 。因此可以针对聚类的异质性进行考察，从而明确聚类估计的可信度。如果聚类空间 $G$ 的分布具有较强的异质性，则当前使用的聚类估计可信度较差，应该采用更稳健的估计量进行调整。同时在小样本情况下应该采取 Jackknife 方法或 Bootstrap 方法。

3.1 高杠杆类 (high leverage)

从聚类的角度考察单个聚类与聚类空间中心的偏离程度。在估计 $β$ 的时候，扔掉聚类 $g^{t h}$ 会导致第 $g^{t h}$ 的回归残差从 ${\hat{u}}_{g}$ 变为 $(I - H_{g})^{- 1} {\hat{u}}_{g}$ ，其中 $H_{g} = X_{g} (X^{'} X)^{- 1} X_{g}^{'}$ 。直接汇报矩阵 $H_{g}$ 并不方便，我们通过汇报矩阵 $H_{g}$ 来测度杠杆，即：

L_{g} = T r (H_{g}) = T r (X_{g} (X^{'} X)^{- 1} X_{g}^{'})

高杠杆的聚类可以通过比较 $L_{g}$ 和 $k / G$ 来识别。如果对于某些聚类 $h$ 有 $L_{h} > k / G$ ，则聚类 $h$ 有可能是高杠杆聚类，在聚类标准误的分析中要谨慎。

3.2 偏杠杆类 (partial leverage)

聚类层面的偏杠杆计算如下：

L_{g j} = \frac{{\dot{x}}_{g j}^{'} {\dot{x}}_{g j}}{{\dot{x}}_{j}^{'} {\dot{x}}_{j}}

其中， ${\dot{x}}_{g j}$ 是 $g^{t h}$ 聚类对应的 $x_{j}$ 中的子向量，通过上式可以简便的计算任何感兴趣的估计系数。如果聚类 $h$ 存在 $L_{h j} >> 1 / G$ ，则表明该聚类对于系数估计具有较高的偏杠杆影响。

3.3 强影响类 (influence)

强影响类聚类是指在删除该聚类后，标准误估计出现较大变化的类，删除聚类 $g$ 之后的 OLS 估计为：

{\hat{β}}^{(g)} = (X^{'} X - X_{g}^{'} X_{g})^{- 1} (X^{'} y - X_{g}^{'} y_{g})

汇报特定参数估计的聚类 $g$ 对应的 ${\hat{β}}^{(g)}$ ，可以明确单个聚类对于参数估计的影响。如果某个聚类对应的 ${\hat{β}}^{(g)}$ 与其他聚类对应的参数 ${\hat{β}}^{(g)}$ 差异较大，则表明聚类 $g$ 可能是强影响类。

通过上述三种表征聚类异质性的方法，MacKinnon 等 (2022) 提供了分析实证分析中聚类可靠性的思路。当 $G$ 较小的时候，应当计算 $L_{g}$ 以及少数参数对应的 $L_{g j}, {\hat{β}}^{(g)}$ 。当 $G$ 较大的时候应当汇报上述三类参数的分布情况 (描述性统计)，从而对于当前的聚类方法进行可靠性评估。

3.4 小样本聚类标准误

聚类稳健标准误 (CV_1) 在样本较小的情况下不具备良好的性质，因此可以采用 Jacknife 方法进行小样本聚类推断。实践中，采用 Jacknife 方法估计 $var (\hat{β})$ 的计算如下：

{CV}_{3 J} : \frac{G - 1}{G} \sum_{g = 1}^{G} ({\hat{β}}^{(g)} - \bar{β}) ({\hat{β}}^{(g)} - \bar{β})^{'}

同时，聚类情况下的 CV_3 标准误计算如下：

{CV}_{3} : \frac{G - 1}{G} \sum_{g = 1}^{G} ({\hat{β}}^{(g)} - \hat{β}) ({\hat{β}}^{(g)} - \hat{β})^{'}

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 计算有效聚类的 $G^{*} ()$ 和 $G^{*} (1)$ 数量；
table 汇报每个聚类的描述性统计结果；
jackknife 汇报 ${CV}_{3 J}$ 标准误等。

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， $p$ 值为 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， $p$ 值为 0.0763，在 95% 水平上不显著，这表明使用聚类稳健标准误 (CV_1) 进行统计推断可能存在可信度的问题。下表汇报了聚类统计量的描述性统计结果：

Ng 表示聚类内部包含样本个数的分布，从聚类样本分布上来看，不同类内部存在较大的异质性。
Leverage 表示杠杆 $L_{g}$ ，从杠杆分布来看最大值和最小值之间存在非常大的差异 (217 倍)，表明不同聚类间处在较强的异质性，可能存在高杠杆类对于统计推断产生影响。
Partial L. 表示偏杠杆 $L_{g j}$ ，从偏杠杆的分布来看极值之间的差异较大 (193 倍)。
beta no g 表示剔除聚类 $g$ 之后的系数估计 ${\hat{β}}^{(g)}$ ，从分布来看相对均衡，剔除某一聚类后并未改变估计系数的符号。

上述分析表明行业层面聚类存在类之间较强的异质性，对于统计推断产生了影响，为了增加实证研究的可信度，应当汇报 ${CV}_{3}$ 或 ${CV}_{3 J}$ 等标准误。


. 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 命令汇报 $L_{g}$ 、 $L_{g j}$ 、 ${\hat{β}}^{(h)}$ 的分布情况，分析在不同类之间是否存在显著的异质性。如果存在的话当前的聚类标准可能存在统计推断的问题，为了提高研究的可信度，应当汇报 ${CV}_{3}$ 或 ${CV}_{3 J}$ 的标准误。

5. 相关推文

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

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

聚类异质性：使用summclust进行统计推断

1. 引言

2. 聚类稳健标准误

3. 聚类标准误异质性识别

3.1 高杠杆类 (high leverage)

3.2 偏杠杆类 (partial leverage)

3.3 强影响类 (influence)

3.4 小样本聚类标准误

4. Stata 实操

4.1 命令介绍

4.2 案例演示

5. 相关推文

相关课程

免费公开课

最新课程-直播课

⛳ 课程主页

⛳ 课程主页

关于我们

热门资讯