Stata：聚类调整后的标准误-Cluster-SE

发布时间：2020-06-11 阅读 30717

Stata 连享会主页 || 视频 || 推文

温馨提示： 定期清理浏览器缓存，可以获得最佳浏览体验。

✌ 课程详情： https://gitee.com/arlionn/Course | lianxh.cn

✨ 课程主页： https://gitee.com/arlionn/Course

杨鑫 (南京大学)，njuyangxin@smail.nju.edu.cn
秦利宾 (厦门大学)，qlb150@163.com
连玉君 (中山大学)，arlionn@163.com

1. 引言

标准误在统计推断中发挥着至关重要的作用，直接影响着系数的显著性和置信区间，并最终影响到假设检验的结论。因此，正确地估计标准误在实证分析的过程中显得尤为重要。当干扰项满足「独立同分布 (iid)」条件时， OLS 所估计的标准误是无偏的。但是当误差项之间存在相关性时，OLS 所估计的标准误是有偏的，不能很好地反映估计系数的真实变异性 (Petersen, 2009)，故需要对标准误进行调整。在多种调整标准误的方式中，「聚类调整标准误 (cluster)」是一种有效的方法 (Petersen, 2009)。

本文主要对聚类调整标准误的原理及其在 Stata 中的具体应用进行简要介绍，包括不同类型的模型中进行「一维聚类调整标准误」和「二维聚类调整标准误」的操作方法。对于该方法更深入的了解，可参考 Petersen (2009)、Thompson (2011)、 Cameron and Miller (2015)、 Abadie et al. (2017) 、Gu and Yoo (2019)等文献。在文章末尾，还对常见的与标准误相关的问题进行了探讨，以便加深对相关内容的理解。

2. 认识标准误

2.1 什么是标准误

为了简便，以仅含有一个非随机解释变量，且不含有截距项回归模型为例予以说明，具体如下：

y_{i} = β x_{i} + u_{i} (1)

其中， $i = 1, \dots, N$ ， $E [u_{i}] = 0$ 。

采用 OLS 方法进行估计，系数的估计量可表示为：

\hat{β} = \sum_{i} x_{i} y_{i} / \sum_{i} x_{i}^{2} (2)

将式 (2) 中的 $y_{i}$ 用式 (1) 替换，整理得：

\hat{β} - β = \sum_{i} x_{i} u_{i} / \sum_{i} x_{i}^{2} (3)

系数方差的一般形式可以表示为：

V [\hat{β}] = E [(\hat{β} - β)^{2}] = V [\sum_{i} x_{i} u_{i}] / {(\sum_{i} x_{i}^{2})}^{2} (4)

若误差项间不相关，则 $V [Σ_{i} x_{i} u_{i}]$ 可以表示为：

V [\sum_{i} x_{i} u_{i}] = \sum_{i} V [x_{i} u_{i}] = \sum_{i} x_{i}^{2} V [u_{i}] (5)

进一步，若「同方差」，则 $V [u_{i}] = σ^{2}$ ，式 (4) 可以表示为：

V [\hat{β}] = σ^{2} / \sum_{i} x_{i}^{2} (6)

若「异方差」，由于 $E [u_{i}] = 0$ ，则 $V [u_{i}] = E [u_{i}^{2}]$ ，式 (4) 可以表示为：

V [\hat{β}] = (\sum_{i} x_{i}^{2} E [u_{i}^{2}]) / {(\sum_{i} x_{i}^{2})}^{2} (7)

White (1980) 认为当 $N \to \infty$ 时， $Σ_{i} x_{i}^{2} E [u_{i}^{2}]$ 可以由 $Σ_{i} x_{i}^{2} {\hat{u}}_{i}^{2}$ 表示，其中， ${\hat{u}}_{i} = y_{i} - \hat{β} x_{i}$ 。

\hat{V} [\hat{β}] = (\sum_{i} x_{i}^{2} {\hat{u}}_{i}^{2}) / {(\sum_{i} x_{i}^{2})}^{2} (8)

这里 $\hat{β}$ 的标准误就是稳健标准误 (robust standard error)，更为准确的表述为异方差稳健标准误 (heteroskedastic-robust standard error)。

若误差项间存在自相关，则 $V [Σ_{i} x_{i} u_{i}]$ 可以表示为：

V [\sum_{i} x_{i} u_{i}] = \sum_{i} \sum_{j} Cov [x_{i} u_{i}, x_{j} u_{j}] = \sum_{i} \sum_{j} x_{i} x_{j} E [u_{i} u_{j}] (9)

V_{c o r} [\hat{β}] = (\sum_{i} \sum_{j} x_{i} x_{j} E [u_{i} u_{j}]) / {(\sum_{i} x_{i}^{2})}^{2} (10)

一个直接的想法是对 White (1980) 扩展，采用 ${\hat{u}}_{i} {\hat{u}}_{j}$ 替代 $E [u_{i} u_{j}]$ ，但是由于 $Σ_{i} x_{i} {\hat{u}}_{i} = 0$ ，使得 $\hat{V} [\hat{β}] = (Σ_{i} Σ_{j} x_{i} x_{j} {\hat{u}}_{i} {\hat{u}}_{j}]) / {(Σ_{i} x_{i}^{2})}^{2}$ 也为 0。

对于时间序列数据，假设误差项在间隔 m 期存在自相关和异方差问题，那么 White (1980) 可以扩展产生异方差自相关一致性估计 (heteroskedastic - and autocorrelation-consistent, HAC)，详见 Newey and West (1987)。

与上述解决同时存在自相关和异方差问题思路类似，聚类标准误 (cluster errors) 假设样本 i 和 j 不在同一组时， $E [u_{i} u_{j}] = 0$ ，可得：

V_{c l u} [\hat{β}] = (\sum_{i} \sum_{j} x_{i} x_{j} E [u_{i} u_{j}] 1 [i, j in same cluster]) / {(\sum_{i} x_{i}^{2})}^{2} (11)

进一步，用 ${\hat{u}}_{i} {\hat{u}}_{j}$ 替代 $E [u_{i} u_{j}]$ ，可得：

{\hat{V}}_{c l u} [\hat{β}] = (\sum_{i} \sum_{j} x_{i} x_{j} {\hat{u}}_{i} {\hat{u}}_{j} 1 [i, j in same cluster]) / {(\sum_{i} x_{i}^{2})}^{2} (12)

其中， $1 [A]$ 为指示函数。在事件 $A$ 发生时，等于 $1$ ，反之为 $0$ 。这里 $\hat{β}$ 的标准误就是聚类稳健标准误 (cluster-robust standard error)。

系数估计量的标准差和标准误是既有联系又有区别的两个统计量：

系数估计量 $\hat{β}$ 的标准差 (standard deviation) 为其方差的平方根：

s d (\hat{β}) = {[V (\hat{β})]}^{1 / 2} (13)

系数估计量 $\hat{β}$ 的标准误 (standard error) 为其方差估计量的平方根：

s e (\hat{β}) = {[\hat{V} (\hat{β})]}^{1 / 2} (14)

2.2 标准误的作用

标准误在统计推断中的作用主要有以下两个方面：

构建 t 统计量。在进行统计推断时，需要构建 t 统计量来对单个参数进行假设检验， ${\hat{β}}_{i}$ 所对应的 t 统计量为：

t_{{\hat{β}}_{i}} = {\hat{β}}_{i} / s e ({\hat{β}}_{i}) (15)

构建置信区间。利用 ${\hat{β}}_{i}$ 的标准误还可以构建总体参数 $β_{i}$ 的置信区间，如 $β_{i}$ 的 95% 置信区间为：

{\hat{β}}_{i} \pm 1.96 \cdot s e ({\hat{β}}_{j}) (16)

3. 聚类调整标准误的基本思想

使用聚类方法调整标准误时，放宽了随机误差项「独立同分布」的假定，要点如下：

允许组内个体的干扰项之间存在相关性；
不同组之间个体的干扰项之间彼此不相关；
系数估计值仍然采用 OLS 估计值，因为它是无偏的。

一维聚类调整标准误

对式 (12) 做进一步处理，可表示为 (Cameron and Miller, 2015)：

{\hat{V}}_{c l u} [\hat{β}] = (\sum_{g = 1}^{G} \sum_{i = 1}^{N_{g}} \sum_{j = 1}^{N_{g}} x_{i g} x_{j g} ω_{i g, j g}) / {(\sum_{i} x_{i}^{2})}^{2} (17)

其中， $G$ 为聚类分组的数量「如行业数量」， $N_{g}$ 为第 $g$ 组样本数量「如某个行业样本量」， $ω_{i g, j g}$ 为 $i g^{t h}$ 和 $j g^{t h}$ 样本的协方差。

二维聚类调整标准误

使用二维聚类调整时，方差估计量由三个方差矩阵计算得来 (Cameron and Miller, 2015 p337) ，其公式的一般表达式为:

{\hat{V}}_{2 w a y} [\hat{β}] = {\hat{V}}_{1} [\hat{β}] + {\hat{V}}_{2} [\hat{β}] - {\hat{V}}_{1 \cap 2} [\hat{β}] (18)

上述公式表明，二维聚类调整的本质即是在两个维度上分别进行一维的聚类调整，再将有交叉的部分去掉。

以行业、年度二维聚类调整为例，标准误调整的步骤可分解为：

估计模型，在行业层面进行聚类调整，计算得到方差矩阵 ${\hat{V}}_{1} [\hat{β}]$ ；
估计模型，在年度层面进行聚类调整，计算得到方差矩阵 ${\hat{V}}_{2} [\hat{β}]$ ；
估计模型，在行业和年度的交互层面（使用行业和年度虚拟变量进行交乘，生成一个新的分组变量）进行聚类调整，计算得到方差矩阵 ${\hat{V}}_{1 \cap 2} [\hat{β}]$ ；
用式 (18) 计算得到二维聚类调整下的方差，并进一步求得标准误。

对于面板数据，在公司、年度层面进行二维聚类调整，方差的估计量还可表示为以下形式 (Thompson, 2011) ：

\hat{V} (\hat{β}) = {\hat{V}}_{f i r m} + {\hat{V}}_{t i m e} - {\hat{V}}_{w h i t e} (19)

其中， ${\hat{V}}_{f i r m}$ 和 ${\hat{V}}_{t i m e}$ 分别表示在公司层面和年度层面进行一维聚类调整的方差， ${\hat{V}}_{w h i t e}$ 表示进行 White (1980) 异方差调整的方差。公式 (19) 可以看作是公式 (18) 的一种特殊形式。

多维聚类调整标准误

多维聚类与式 (18) 类似，以三维聚类为例 (Gu and Yoo, 2019)：

{\hat{V}}_{3 w a y} [\hat{β}] = {\hat{V}}_{1} [\hat{β}] + {\hat{V}}_{2} [\hat{β}] + {\hat{V}}_{3} [\hat{β}] - {\hat{V}}_{1 \cap 2} [\hat{β}] - {\hat{V}}_{1 \cap 3} [\hat{β}] - {\hat{V}}_{2 \cap 3} [\hat{β}] + {\hat{V}}_{1 \cap 2 \cap 3} [\hat{β}] (20)

4. Stata 实操

4.1 一维聚类调整

对标准误进行一维聚类调整时， Stata 命令有如下几种表现形式：

/*
*-截面数据，在公司层面进行聚类，以下两种写法等价
  reg y x, cluster(id) 
  reg y x, vce(cluster id)  
  
*-面板数据，在公司层面进行聚类，以下三种写法等价
  xtset id year
  xtreg y x, fe cluster(id)  
  xtreg y x, fe vce(cluster id)
  xtreg y x, fe robust  // If you specify -xtreg, fe robust-, Stata will automatically, and without even telling you, use vce(cluster panel_variable) instead. (This is true since version 13.)
*/

以 nlswork.dta 为例，对 Stata 相关命令和结果予以说明。

*-调入数据
  *copy http://www.stata-press.com/data/r9/nlswork.dta nlswork.dta, replace
  use nlswork.dta, clear
  
*-定义全局暂元
  global x "age grade"

*-回归结果
  reg ln_wage $x          //干扰项同方差
  est store m1
  reg ln_wage $x, robust  //干扰项异方差
  est store m2
  reg ln_wage $x, vce(cluster idcode)
  est store m3

4.2 二维聚类调整

在对标准误进行二维聚类调整时， Stata 命令有以下几种不同形式：

*-cluster2 (Petersen-2009, RFS) 
  cluster2 ln_wage $x, fcluster(idcode) tcluster(year) 
  *该命令没有帮助文件，所有功能都可以用 cgmreg 和 vce2way 代替
  *因此，建议日后不必使用该命令

*-cgmreg (CGM2011, Mitchell Petersen's -cluster2.ado- 的升级版)
  *需手动下载：https://sites.google.com/site/judsoncaskey/data
  *help cgmreg  
  cgmreg ln_wage $x, cluster(idcode year)	
  est store m4

*-vce2way (CGM2011, 支持 Panel data, xtreg 等命令)
  *ssc install vce2way
  *help vce2way
  vce2way reg ln_wage $x, cluster(idcode year)	
  est store m5
  
*-vcemway (Gu and Yoo-2019, 该命令在 vce2way 的基础上扩展到多维)
  *ssc install vcemway
  *help vcemway
   vcemway reg ln_wage $x, cluster(idcode year)	
   est store m6

*-结果对比
  local m  "m1 m2 m3 m4 m5 m6"
  local mt "OLS Robust 1Clus 2_cgmreg 2_vce2way 2_vcemway"
  esttab `m', mtitle(`mt') nogap b(%4.3f) se(%6.4f) brackets ///
	      star(* 0.1 ** 0.05 *** 0.01) s(N r2) compress

Note: 对于 Logit、Probit 和 Tobit 模型的二维聚类，可以使用 logit2、probit2 和 logit2 实现，其用法与 cluster2 类似。上述命令的 ado 文件，均可从 Mitchell A. Petersen 的个人主页进行下载。

----------------------------------------------------------------------------------------
                 (1)          (2)          (3)          (4)          (5)          (6)   
                 OLS       Robust        1Clus     2_cgmreg    2_vce2way    2_vcemway   
----------------------------------------------------------------------------------------
age            0.015***     0.015***     0.015***     0.015***     0.015***     0.015***
            [0.0006]     [0.0010]     [0.0004]     [0.0004]     [0.0011]     [0.0011]   
grade          0.083***     0.083***     0.083***     0.083***     0.083***     0.083***
            [0.0022]     [0.0030]     [0.0011]     [0.0011]     [0.0035]     [0.0035]   
_cons          0.218***     0.218***     0.218***     0.218***     0.218***     0.218***
            [0.0289]     [0.0345]     [0.0165]     [0.0165]     [0.0419]     [0.0419]   
----------------------------------------------------------------------------------------
N            2.9e+04      2.9e+04      2.9e+04      2.9e+04      2.9e+04      2.9e+04   
r2             0.233        0.233        0.233        0.233        0.233        0.233   
----------------------------------------------------------------------------------------
Standard errors in brackets
* p<0.1, ** p<0.05, *** p<0.01

4.3 二维聚类调整逻辑分析

根据式 (18) 和 (19) 手工分步计算二维聚类标准误：

*-分别进行 idcode、year、id_year 一维聚类
  *在 idcode 维度进行一维聚类
  reg ln_wage $x, cluster(idcode)
  est store m1
  
  *在 year 维度进行一维聚类
  reg ln_wage $x, cluster(year)
  est store m2
  
  *在 idcode 和 year 交互维度进行一维聚类
  egen id_year = group(idcode year)
  reg ln_wage $x, cluster(id_year)
  est store m3

*-异方差调整
  reg ln_wage $x, robust
  est store m4

*-vcemway 计算
  vcemway reg ln_wage $x, cluster(idcode year)
  est store m5

*-结果对比
  local m  "m1 m2 m3 m4 m5"
  local mt "Clu_id Clu_year Clu_id_year Robust 2_vcemway"
  esttab `m', mtitle(`mt') nogap b(%4.3f) se(%6.4f) brackets ///
	      star(* 0.1 ** 0.05 *** 0.01) s(N r2) compress
			  
*-手工计算二维聚类标准误
  est restore m1
  scalar se_idcode  = _se[age]
  est restore m2
  scalar se_year    = _se[age]
  est restore m3
  scalar se_id_year = _se[age]
  est restore m4
  scalar se_robust = _se[age]
  *式 (18)
  scalar se_2way1 = sqrt(se_idcode^2+se_year^2-se_id_year^2)
  *式 (19)
  scalar se_2way2 = sqrt(se_idcode^2+se_year^2-se_robust^2)
  scalar list se_2way1 se_2way1


---------------------------------------------------------------------------
                 (1)          (2)          (3)          (4)          (5)   
              Clu_id     Clu_year    Clu_id_~r       Robust    2_vcemway   
---------------------------------------------------------------------------
age            0.015***     0.015***     0.015***     0.015***     0.015***
            [0.0006]     [0.0010]     [0.0004]     [0.0004]     [0.0011]   
grade          0.083***     0.083***     0.083***     0.083***     0.083***
            [0.0022]     [0.0030]     [0.0011]     [0.0011]     [0.0035]   
_cons          0.218***     0.218***     0.218***     0.218***     0.218***
            [0.0289]     [0.0345]     [0.0165]     [0.0165]     [0.0419]   
---------------------------------------------------------------------------
N            2.9e+04      2.9e+04      2.9e+04      2.9e+04      2.9e+04   
r2             0.233        0.233        0.233        0.233        0.233   
---------------------------------------------------------------------------
Standard errors in brackets
* p<0.1, ** p<0.05, *** p<0.01

. scalar list se_2way1 se_2way1
  se_2way1 =  .00108169
  se_2way1 =  .00108169

我们在这仅列示了「age」的标准误，可以看出手工计算均为 .00108169，与我们采用命令 vcemway 计算的结果保持一致。

5. FAQs

如下问题基于连享会课程群学员提问和助教解答整理。
往期学员问答 (WD) 可前往 https://gitee.com/arlionn/WD 网站查看。
若需加入课程群，可以扫描如下二维码：

Q1： 我们什么时候需要使用聚类调整？

研究中应该评估采样过程和分配机制是否聚类的，如果两者回答都是“否”，则无论该调整是否会改变标准误差，都不应该进行聚类调整 (Abadie et al., 2017)。

Q2: 一维聚类和二维聚类如何选择？二维聚类一定比一维聚类更优吗？

更稳健的标准误会降低统计推断的偏差，但同时也会使方差增大，使结果更加不显著，增大犯第二类错误的概率 (Thompson, 2011)。因此，如何在二者之间进行权衡、选择何种聚类方式则应当根据具体的数据结构和逻辑来进行判断。

Q3: 进行聚类调整后， t 值达到了大样本条件下的临界值，但 p 值却未达到相应的显著性水平，这是何种原因？如 t 值为 1.67，但 p 值却大于 0.1 。

由于计算 p 值的公式中其中一个参数为自由度，而聚类调整会影响到模型的自由度，故而会影响到最终计算得到的 p 值。

Q4： Fixed Effect 和 cluster 的区别，如控制了企业固定效应，同时也在企业层面进行了聚类调整。

企业固定效应是控制了企业不随时间变化的特征，而企业层面的 cluster 调整则是认为误差项在企业层面存在相关性。

Q5： 聚类回归结果无 F 和 P值。

自变量的个数必须小于聚类的个数，否则没有 F 值及 P 值，需要重新完善模型。详见 Missing F-statistic when using xtreg with fe, vce(cluster) and adding time-fixed effects。

Q6： 聚类稳健标准误回归中，聚类只有20个，对结果是否有影响？聚类是否要不少于 50 时，使用聚类稳健标准误才有效？

是这样的，否则聚类标准误无效。详见 Problem with small number of clusters using reghdfe and vce suboptions、How misleading are clustered SEs in designs with few clusters?、Beware of studies with a small number of clusters。

Q7： 双维 cluster 修正标准误，是不是只在固定效应中使用？

随机效应也可以使用，详见 vcemway 命令。

Q8： 为什么换了 cluster 对象，系数也变了？

无论对标准误作何处理，该变的只有标准误，系数是不该变。如果发现调整 cluster 对象系数改变，很可能是样本发生改变。如 cluster(id) 和 cluster(industry) 不同的话，和可能是 id 或 industry 存在缺失值。

6. 参考文献

[1] Petersen, M. A. 2009. Estimating Standard Errors in Finance Panel Data Sets: Comparing Approaches. Review of Financial Studies, 22(1): 435-480. [PDF]

[2] Thompson, S. B., 2011, Simple formulas for standard errors that cluster by both firm and time, Journal of Financial Economics, 99 (1): 1-10. [PDF]

[3] Cameron, C. A., D. L. Miller, 2015, A practitioner’s guide to cluster-robust inference, Journal of Human Resources, 50 (2): 317-372. [PDF]

[4] Abadie, A., S. Athey, G. W. Imbens, J. Wooldridge, 2017, When should you adjust standard errors for clustering?, Working Paper. [PDF]

[5] Gu, A. and Yoo, H. I., 2019, Vcemway: A One-Stop Solution for Robust Inference with Multiway Clustering, The Stata Journal, 19(4): pp.900-912. [PDF]

专题	嘉宾	直播/回看视频
⭐ 最新专题 ⭐		DSGE, 因果推断, 空间计量等
⭕ Stata数据清洗	游万海	直播, 2 小时，已上线
研究设计	连玉君	我的特斯拉-实证研究设计，-幻灯片-
面板模型	连玉君	动态面板模型，-幻灯片-
面板模型	连玉君	直击面板数据模型 [免费公开课，2小时]