Stata：面板聚类标准误-自动确定最优聚类层级和数量-xtregcluster

发布时间：2021-11-06 阅读 7833

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

作者：吕卓阳 (厦门大学)
邮箱：lvzy20@163.com

1. 研究背景

在针对上市公司的研究中，大家经常在如下两种估计命令之间摇摆不定：第一种是我们控制了时间固定效应、以及公司或行业固定效应，进行面板回归：

. xtreg y x1 x2 i.year, fe 
. reg   y x1 x2 i.year i.industry

第二种是在第一种的基础上，考虑标准误的聚类：

. xtreg y x1 x2 i.year, fe vce(cluster firmID)
. reg   y x1 x2 i.year i.industry, vce(cluster industry)

回归的方程如下：

y_{i t} = α_{i} + λ_{t} + x_{i t} β + ε_{i t}

y_{i t} = θ_{j} + λ_{t} + x_{i t} β + ε_{i t}

两种回归的方式不同在于第二种回归考虑了标准误的聚类问题，即假设标准误在公司或行业的层面聚类。这种假设相比于假定标准误差满足 $i . i . d$ 的基础上更加合理，因为实际问题中标准误在某一层面上允许其相关是合理的，譬如行业内的标准误存在相关性，而行业间标准误不相关。

我们平时常常是基于经济学的意义预先设定聚类的数量和分块大小，鲜少从统计意义上去估计最优聚类方法，本文将为大家介绍一种新的「聚类分析」方法及其 Stata 实现命令 xtregcluster。其优点在于不需要先验地假设聚类的数量和分块的大小，而是根据数据建立惩罚函数来降低总残差平方和，本着「让数据说话」的理念寻找最优聚类数量和分块大小。

2. 方法介绍

对于标准的面板数据分析，一般施加所有个体具有相同斜率系数的限制，个体不可观测的异质性取决于个体不随时间改变的效应，但从理论与实证上都很难证明该限制 (Burnside，1996；Baltagi 和 Griffin，1997；Pesaran 等，1999)。

由此，Sarafidis 和 Weber (2015) 开发了一种回归聚类方法，将个体分为不同的聚类。这种聚类不依赖于参数分布的先验假设。譬如，我们想探究学生成绩的影响因素，若在计算标准误时，选择在班级层面上进行聚类调整，则意味着我们假定同一个班级的同学的成绩更具有相关性 (可能源于班风、班主任个人特征等因素)。但这类假设往往是基于研究者的理论分析给出的假设，而非通过数据分析得到的。

Sarafidis 开发的回归聚类方法是基于数据确定的，假设聚类数量与分块的大小都是未知，基于最小化总残差平方和的最优化方法并设置严格增加聚类数量的惩罚函数，从而将个体分为不同的聚类。在每个聚类中，斜率系数是同质的，而在组间斜率的系数是异质的。在没有参数结构的先验信息的情况下，xtregcluster 可以用于探索斜率参数的异质性，适用于分析大样本个体和短时间序列样本。此外，从统计意义上也可以检验单个实体先验分类是否最优。

3. 命令介绍

3.1 理论部分

考虑如下面板回归模型：

y_{i_{ω} t} = β_{ω}^{'} x_{i_{ω} t} + ϵ_{i_{ω} t}

其中， $y_{i_{ω} t}$ 表示个体 $i$ 在 $w$ 标准误聚类情形下 $t$ 时刻被解释变量的观测值， $X_{i_{ω} t}$ 表示 $K \times 1$ 阶在标准误聚类情形下的解释变量， $β_{ω}$ 为 $K \times 1$ 阶的估计系数，并且由于在 $w$ 情形下聚类，系数在组内相同但在组间不同。残差项可以分解如下：

ϵ_{i_{ω} t} = η_{i_{ω}} + τ_{t} + e_{i_{ω} t}

其中， $e_{i_{ω} t}$ 为个体不可观测的残差项， $τ_{t}$ 表示时间固定效应， $e_{i_{ω} t}$ 可能与 $η_{i_{ω}}$ 相关。譬如背景介绍中的例子，我们想要在公司层面研究 $X$ 对 $Y$ 的影响，但在残差项满足 $i . i . d$ 的情形下我们忽视了残差项中存在相关性的问题，譬如同一行业的公司存在相关性，从而导致模型的偏误。

我们过往的聚类方法往往通过经济学角度进行预先设定，譬如我们假定公司内部或行业内部存在组内相关，从而我们在公司或行业的级别上聚类，但我们还可以通过统计的角度，通过最小化 RSS 来确定最优聚类分组。

由于在整个时间的某个个体只能属于一个集群，在 $ω_{0}$ 未知以及个体聚类关系未知的情况下，Sarafidis 和 Weber (2015) 开发了一种新的方法以估计模型系数的估计值，他们得到了 $ω_{0}$ 的强一致估计，即对于任意固定的 $T$ ，当 $prob \to 1$ 且 $N \to \infty$ ， $\hat{Ω} \to Ω_{0}$ 。Sarafidis 和 Weber 最优估计聚类方式是最小化残差平方和，即：

R S S_{ω} = {∥ {\tilde{Y}}_{C_{Ω ω}} - {\tilde{X}}_{C_{Ω ω}} {\hat{β}}_{Ω ω} ∥}^{2} R S S = RSS (Ω) = \sum_{ω = 1}^{Ω} {R S S}_{ω}

我们定义模型信息准则 (MIC)：

M I C = N \ln (\frac{R S S}{N \bar{T}}) + Ω θ_{N}

其中， $f (Ω)$ 是关于 $Ω$ 的严格增函数， $θ_{N}$ 是取决于 $N$ 的常数， $Ω θ_{N}$ 表示用于限制聚类数量的惩罚项。由此，我们可以通过最小化 MIC 来估计聚类数量和分区。

$Ω$ 最优值 $\hat{Ω}$ 具体的估计策略如下：

给定 $Ω$ 初始值，其中 $Ω \leq ξ \leq Ω_{0}$ ，且 $ξ > 1$ ， $ξ \in Z^{+}$ ，我们常常假定 $f (Ω) = Ω$ 且 $θ_{N} = \sqrt{N}$ ；
随机生成初始识别值；
可以使用如下三种方式获得初始分区：i) 标准均匀分布随机生成，ii) 基于先验分布，iii) 基于给定观测使用 k-means 聚类；
将初始个体再分配给现有的聚类，每次保存 RSS 值，迭代得到使 $R S S = \sum_{ω = 1}^{Ω} R S S_{ω}$ 最小的值；
每个个体重复第四步；
重复第三步与第四步使 RSS 无法再进一步减小；
给定不同的 $Ω$ 重复第一步到第六步，最优的 $\hat{Ω}$ 是使模型信息准则 MIC 最小的值。

3.2 语法结构

与所有的 xt 命令相同，xtregcluster 要求数据先加以设定。首先通过随机分类获得初始分区、预分类或根据特定观测分类，再将个体重新聚类，从而使 RSS 最小化。

*命令安装
cnssc install xtregcluster, replace

*命令语法
xtregcluster depvar indepvars if in weight , 
    {random|preclass(varname)|prevars(varlist|X|b)} omega(numlist)
    [prevarsopt(kmeansopt) theta(#) seed(#) name(varname) iterate(#)
     tolerance(#) nolog graph table]

random：从标准均匀分布中随机选取初始分区；
preclass(varname)：根据分类变量获得初始分区；
prevars(varlist)：使用 Stata 官方命令 kmeans 聚类方法获得初始分区，prevars(X) 包含所有的回归自变量，prevars(b) 使用所有个体的估计斜率；
omega(numlist)：定义 $Ω > 1$ 的数值范围；
prevarsopt(kmeansopt)：与 prevars(varlist) 同时指定；
theta(#)：惩罚函数的参数 $θ_{N}$ ；
seed：随机数种子，默认为 seed(123)；
name：存储识别最优聚类分区数量的变量名称，默认为 name(omega#)；
iterate：最小化RSS的最大迭代数量，默认值为 100；
tolerance：收敛于 RSS 的容忍度，默认值为 tolerance(1e-6;
nolog：取消迭代日志；
graph：不同聚类方法的散点图与线性拟合；
tolerance：显示每个聚类的斜率系数，回报基于 $\hat{Ω}$ 使用聚类标准误的固定效应回归结果。

4. Stata 实例

我们使用西班牙奶牛生产产量估计作为示例，相关数据来源于网站「William Greene 的面板数据计量经济学」。面板数据是 1993-1998 年 247 个农场产出与若干投入的观测值，包括牛乳、奶牛数量、土地大小、劳动力和饲料。

我们首先需要使用 xtset 设定面板数据。其中， $y_{i t}$ 表示的是产出 (牛乳) 的对数值， $x_{k l}$ 是投入的 $K$ 个变量的对数值。我们先基于 $Ω = 2, 3$ 使用标准均匀分布的初始分布进行估计。

. lxhuse dairy_wg, clear
. xtset farm year
. xtregcluster yit x1-x34, random omega(2/3)

Initial partition via randomized classification and seed 123
Omega = 2
Iteration 0:   Total RSS =           7.654266
...省略
Iteration 5:   Total RSS =           6.316776

Omega = 3
Iteration 0:   Total RSS =           7.579113
...省略
Iteration 9:   Total RSS =           5.461213
+-----------------------------------------------+
|     Omega     |    Total RSS  |      MIC      |
|---------------+---------------+---------------|
|       1       |         7.887 |    -1280.962  |
|       2       |         6.317 |    -1323.483  |
|       3       |         5.461 |    -1347.117  |
+-----------------------------------------------+
Proceed with xtreg if omega3==`i',fe  where `i'=1,2,3

由估计结果可知，在此模型设定下本数据存在至少三个聚类分区，此外，omega2 和 omega3 指示了所有个体的聚类指示。而 MIC 最小值的 $Ω$ 等于 omega() 的最大值。由此，我们需要进一步探索最优解，可以通过nolog 来节约空间。我们同样可以指定 om 名称来存储最优的聚类分区数量。

. xtregcluster yit x1-x34,random omega(2/10) name(om) nolog

Initial partition via randomized classification and seed 123
+-----------------------------------------------+
|     Omega     |    Total RSS  |      MIC      |
|---------------+---------------+---------------|
|       1       |         7.887 |    -1280.962  |
|       2       |         6.317 |    -1323.483  |
|       3       |         5.461 |    -1347.117  |
|       4       |         4.999 |    -1356.638  |
|       5       |         4.626 |    -1363.494  |
|       6       |         4.483 |    -1358.955  |
|       7       |         4.001 |    -1374.719  |
|       8       |         3.820 |    -1373.825  |
|       9       |         3.781 |    -1364.031  |
|       10      |         3.580 |    -1365.253  |
+-----------------------------------------------+
Proceed with xtreg if om7==`i',fe  where `i'=1,2,3,4,5,6,7

由结果可知，最优的 $Ω$ 为7，这是基于初始分布为标准均匀分布的情形下。当然，我们也可以使用自变量 $X$ 来作为特定观测。

. xtregcluster yit x1-x34,prevars(X) omega(2/10) name(omX) nolog

Initial partition via the variation in x1 x2 x3 x4 x11 x22 x33 x44 x12 x13 x14 x23 x24 
> x34 and seed 123
+-----------------------------------------------+
|     Omega     |    Total RSS  |      MIC      |
|---------------+---------------+---------------|
|       1       |         7.887 |    -1280.962  |
|       2       |         6.155 |    -1329.891  |
|       3       |         5.658 |    -1338.374  |
|       4       |         4.976 |    -1357.779  |
|       5       |         4.684 |    -1360.402  |
|       6       |         4.374 |    -1364.992  |
|       7       |         4.132 |    -1366.781  |
|       8       |         3.863 |    -1371.059  |
|       9       |         3.805 |    -1362.471  |
|       10      |         3.574 |    -1365.671  |
+-----------------------------------------------+
Proceed with xtreg if omX8==`i',fe  where `i'=1,2,3,4,5,6,7,8

由结果可知，最优的 $Ω$ 为 8，且 MIC 要大于标准均匀分布的 MIC。由此，使用不同初始分布设定的结果是不同的。此外，我们可以使用 tabulate 查看每个聚类中分配了多少人。我们也可以使用 seed() 设定不同的随机数种子。table 和 graph 参数可以提供更多的模型信息，table 提供模型的估计，而 graph 包含一系列的散点拟合图。

* Pooled fitted values
quietly xtreg yit x1-x34,fe
predict xb,xb
* Fitted values by cluster
forvalues i = 1/7{
  quietly xtreg yit x1-x34 if om7==`i', fe
  estimates store om`i'
  quietly predict xb`i' if om7==`i ', xb
  twoway (lfit yit xb, lwidth(*3) lpattern(solid) lcolor(gs8)) ///
    (scatter yit xb`i' if om7==`i', msymbol(oh) mlwidth(*.3) mcolor(gs0))   ///
    (lfit yit xb`i' if om7==`i', lwidth(*2.25) lpattern(dash) lcolor(gs0)), ///
    aspect(1) ysize(1) xsize(1) scheme(sj) legend(off) ///
    ytitle("Log of milk production (output)")          ///
    title("{&omega} = `i'", ring(0) pos(11) margin(medium)) name(g`i', replace)
}
graph combine g1 g2 g3 g4 g5 g6 g7, row(2) ysize(2) xsize(4) imargin(small) scheme(sj)

实线拟合线表示混合面板回归的拟合线，而虚线表示每个 omega 值的线性拟合，除了 omega=2 的情形，其余的线性拟合均相较于混合回归斜率更小。

5. 参考文献

Christodoulou D, Sarafidis V. Regression clustering for panel-data models with fixed effects[J]. The Stata Journal, 2017, 17(2): 314-329. -PDF-
Sarafidis V, Weber N. A partially heterogeneous framework for analyzing panel data[J]. Oxford Bulletin of Economics and Statistics, 2015, 77(2): 274-296. -PDF-

6. 相关推文

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

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

Stata：面板聚类标准误-自动确定最优聚类层级和数量-xtregcluster

1. 研究背景

2. 方法介绍

3. 命令介绍

3.1 理论部分

3.2 语法结构

4. Stata 实例

5. 参考文献

6. 相关推文

相关课程

免费公开课

最新课程-直播课

⛳ 课程主页

⛳ 课程主页

关于我们

热门资讯