Stata：面板聚类标准误-自动确定最优聚类层级和数量-xtregcluster
👇 连享会 · 推文导航 | www.lianxh.cn
🍎 Stata：Stata基础 | Stata绘图 | Stata程序 | Stata新命令 📘 论文：数据处理 | 结果输出 | 论文写作 | 数据分享 💹 计量：回归分析 | 交乘项-调节 | IV-GMM | 时间序列 | 面板数据 | 空间计量 | Probit-Logit | 分位数回归 ⛳ 专题：SFA-DEA | 生存分析 | 爬虫 | 机器学习 | 文本分析 🔃 因果：DID | RDD | 因果推断 | 合成控制法 | PSM-Matching 🔨 工具：工具软件 | Markdown | Python-R-Stata 🎧 课程：公开课-直播 | 计量专题 | 关于连享会
🎦 效率分析前沿模型 最新模型和工具
📅 2021 年 11 月 13-14；20-21 日
🔑 张宁教授 (山东大学)；杜克锐副教授 (厦门大学)
🍓 课程主页：https://gitee.com/lianxh/TE
实现经济增长的途径主要有两个：一是增加投入，二是提高效率。后者不仅能增加经济总量，还能提升经济质量。因此，效率分析是经济学研究中的核心。
显然，若想发表高水平的论文，就必须紧跟前沿，引领趋势，这正是本期「效率分析前沿模型」课程的主旨。
张宁教授在环境效率和非期望产出领域学者中排名第二，他在 Science, Nature, Cell 和 Lancet 上均有发表。杜克锐副教授已在 SSC 发布了十余个 Stata 命令，有三篇论文发表于 Stata Journal，获得 2020 年 Stata 中国用户奖。
课程包括参数方法和非参数方法两个部分：
参数方法主要包括参数线性规划（Parametric LP）、计量经济方法和随机前沿方法（SFA）三大方法。 非参数方法，即 DEA 方法，主要包括 DDF, NDDF 和 SBM 等模型。
重要福利： 本课程所涉及的参数和非参数效率模型，老师们都将会提供封装后的 Stata 估计命令，只需几条简单的命令即可实现最新的模型。老师们也会分享编写过程心得，让大家有能力对这些封装程序进行修改和扩展。
扫码直达课程主页：
作者：吕卓阳 (厦门大学)
邮箱：lvzy20@163.com
目录
1. 研究背景
2. 方法介绍
3. 命令介绍
3.1 理论部分
3.2 语法结构
4. Stata 实例
5. 参考文献
6. 相关推文
温馨提示： 文中链接在微信中无法生效。请点击底部「阅读原文」。或直接长按/扫描如下二维码，直达原文：
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)
回归的方程如下：
两种回归的方式不同在于第二种回归考虑了标准误的聚类问题，即假设标准误在公司或行业的层面聚类。这种假设相比于假定标准误差满足 的基础上更加合理，因为实际问题中标准误在某一层面上允许其相关是合理的，譬如行业内的标准误存在相关性，而行业间标准误不相关。
我们平时常常是基于经济学的意义预先设定聚类的数量和分块大小，鲜少从统计意义上去估计最优聚类方法，本文将为大家介绍一种新的「聚类分析」方法及其 Stata 实现命令 xtregcluster。其优点在于不需要先验地假设聚类的数量和分块的大小，而是根据数据建立惩罚函数来降低总残差平方和，本着「让数据说话」的理念寻找最优聚类数量和分块大小。
2. 方法介绍
对于标准的面板数据分析，一般施加所有个体具有相同斜率系数的限制，个体不可观测的异质性取决于个体不随时间改变的效应，但从理论与实证上都很难证明该限制 (Burnside，1996；Baltagi 和 Griffin，1997；Pesaran 等，1999)。
由此，Sarafidis 和 Weber (2015) 开发了一种回归聚类方法，将个体分为不同的聚类。这种聚类不依赖于参数分布的先验假设。譬如，我们想探究学生成绩的影响因素，若在计算标准误时，选择在班级层面上进行聚类调整，则意味着我们假定同一个班级的同学的成绩更具有相关性 (可能源于班风、班主任个人特征等因素)。但这类假设往往是基于研究者的理论分析给出的假设，而非通过数据分析得到的。
Sarafidis 开发的回归聚类方法是基于数据确定的，假设聚类数量与分块的大小都是未知，基于最小化总残差平方和的最优化方法并设置严格增加聚类数量的惩罚函数，从而将个体分为不同的聚类。在每个聚类中，斜率系数是同质的，而在组间斜率的系数是异质的。在没有参数结构的先验信息的情况下，xtregcluster 可以用于探索斜率参数的异质性，适用于分析大样本个体和短时间序列样本。此外，从统计意义上也可以检验单个实体先验分类是否最优。
3. 命令介绍
3.1 理论部分
考虑如下面板回归模型：
其中， 表示个体 在 标准误聚类情形下 时刻被解释变量的观测值， 表示 阶在标准误聚类情形下的解释变量， 为 阶的估计系数，并且由于在 情形下聚类，系数在组内相同但在组间不同。残差项可以分解如下：
其中， 为个体不可观测的残差项， 表示时间固定效应， 可能与 相关。譬如背景介绍中的例子，我们想要在公司层面研究 对 的影响，但在残差项满足 的情形下我们忽视了残差项中存在相关性的问题，譬如同一行业的公司存在相关性，从而导致模型的偏误。
我们过往的聚类方法往往通过经济学角度进行预先设定，譬如我们假定公司内部或行业内部存在组内相关，从而我们在公司或行业的级别上聚类，但我们还可以通过统计的角度，通过最小化 RSS 来确定最优聚类分组。
由于在整个时间的某个个体只能属于一个集群，在 未知以及个体聚类关系未知的情况下，Sarafidis 和 Weber (2015) 开发了一种新的方法以估计模型系数的估计值，他们得到了 的强一致估计，即对于任意固定的 ，当 且 ，。Sarafidis 和 Weber 最优估计聚类方式是最小化残差平方和，即：
我们定义模型信息准则 (MIC)：
其中， 是关于 的严格增函数， 是取决于 的常数， 表示用于限制聚类数量的惩罚项。由此，我们可以通过最小化 MIC 来估计聚类数量和分区。
最优值具体的估计策略如下：
给定 初始值，其中 ，且 ，，我们常常假定 且 ； 随机生成初始识别值； 可以使用如下三种方式获得初始分区：i) 标准均匀分布随机生成，ii) 基于先验分布，iii) 基于给定观测使用 k-means 聚类； 将初始个体再分配给现有的聚类，每次保存 RSS 值，迭代得到使 最小的值； 每个个体重复第四步； 重复第三步与第四步使 RSS 无法再进一步减小； 给定不同的 重复第一步到第六步，最优的 是使模型信息准则 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)：定义 的数值范围；prevarsopt(kmeansopt)：与prevars(varlist)同时指定；theta(#)：惩罚函数的参数 ；seed：随机数种子，默认为seed(123)；name：存储识别最优聚类分区数量的变量名称，默认为name(omega#)；iterate：最小化RSS的最大迭代数量，默认值为 100；tolerance：收敛于 RSS 的容忍度，默认值为tolerance(1e-6;nolog：取消迭代日志；graph：不同聚类方法的散点图与线性拟合；tolerance：显示每个聚类的斜率系数，回报基于 使用聚类标准误的固定效应回归结果。
4. Stata 实例
我们使用西班牙奶牛生产产量估计作为示例，相关数据来源于网站「William Greene 的面板数据计量经济学」。面板数据是 1993-1998 年 247 个农场产出与若干投入的观测值，包括牛乳、奶牛数量、土地大小、劳动力和饲料。
我们首先需要使用 xtset 设定面板数据。其中， 表示的是产出 (牛乳) 的对数值， 是投入的 个变量的对数值。我们先基于 使用标准均匀分布的初始分布进行估计。
. 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，这是基于初始分布为标准均匀分布的情形下。当然，我们也可以使用自变量 来作为特定观测。
. 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
专题：Stata命令 Stata：聚类分析-cluster 专题：回归分析 acreg：允许干扰项随意相关的稳健性标准误 Stata：聚类调整标准误笔记 Stata：聚类调整后的标准误-Cluster-SE 小样本下OLS估计的纠偏聚类标准误 专题：面板数据 wcbregress：面板聚类标准误 专题：倍分法DID 倍分法(DID)的标准误：不能忽略空间相关性
New！ Stata 搜索神器：
lianxh和songblGIF 动图介绍
搜： 推文、数据分享、期刊论文、重现代码 ……
👉 安装：
. ssc install lianxh
. ssc install songbl
👉 使用：
. lianxh DID 倍分法
. songbl all
🍏 关于我们
连享会 ( www.lianxh.cn，推文列表) 由中山大学连玉君老师团队创办，定期分享实证分析经验。 直通车： 👉【百度一下：连享会】即可直达连享会主页。亦可进一步添加 「知乎」,「b 站」,「面板数据」,「公开课」 等关键词细化搜索。