Stata: gen 命令中的 group() 函数的潜在风险

查看原文

其他

Stata: gen 命令中的 group() 函数的潜在风险

Stata连享会 Stata连享会 2020-02-10

作者：连玉君 (知乎 | 简书 | 码云)

python 爬虫与文本分析专题-现场班

连享会-Python爬虫与文本分析现场班-山西大学 2019.5.17-19

一起学空间计量……

空间计量专题-西安 2019.6.27-30

1. 问题背景

我们经常使用 generate (后文简称 gen) 命令提供的 group() 函数对某个变量进行分组，产生分组变量 gg，继而基于 gg 变量进行后续的分组回归分析。

例如，在公司金融中，常用如下代码产生融资约束的分组指标：

*-用公司规模衡量融资约束，分成三组
. bysort code: egen av_size = mean(size)
. sort av_size
. gen gg = group(3)
. gen FC = (gg==1) //小规模公司定义为 FC 组
. replace FC=. if gg==2 //丢弃中间组
. reg y x if FC==0
. est store m_FC
. reg y x if FC==1
. est store m_NFC
. esttab m_FC m_NFC

恐怖的事情就要发生了！

后续做分组回归时，你会发现：执行相同的代码，但两次得到的结果居然不同！

原因何在？？

一个小例子

我先虚构一份数据，让各位了解 group() 函数的工作原理，搞明白这件事情后，上面的问题就迎刃而解了。

这份数据很简单，只有 4 行观察值。我们对变量 x 排序后再执行 gen g = group(2) 命令，以便将样本分成两组。为了测试分组结果是否唯一\稳定，我进一步使用 tatstat 命令计算了两组的均值。

clear
input x y
3 13
2 10
1 1
2 8
end
sort x \\ 由小到大排序
gen g = group(2) \\ 等分两组
tabstat x y, by(g) f(%3.1f)
sort g x y
list, sepby(g) noobs

第一轮执行的结果如下：

. tabstat x y, by(g) f(%3.1f)
Summary statistics: mean
by categories of: g
g | x y
---------+--------------------
1 | 1.5 5.5
2 | 2.5 10.5
---------+--------------------
Total | 2.0 8.0
------------------------------
. sort g x y
. list, sepby(g) noobs
+------------------+
| id x y g |
|------------------|
| 301 1 1 1 |
| 201 2 10 1 |
|------------------|
| 401 2 8 2 |
| 101 3 13 2 |
+------------------+

可以看到，g 取值为 1 和 2 时两组的均值分别为 5.5 和 10.5，差为 5。

我又连续执行了两遍上述代码，结果都没有任何变化，似乎表明上述分组结果是稳定的。

然而，当我执行第四次时，得到了如下结果：

. tabstat x y, by(g) f(%3.1f)
Summary statistics: mean
by categories of: g
g | x y
---------+--------------------
1 | 1.5 4.5
2 | 2.5 11.5
---------+--------------------
Total | 2.0 8.0
------------------------------
. sort g x y
. list, sepby(g) noobs
+------------------+
| id x y g |
|------------------|
| 301 1 1 1 |
| 401 2 8 1 |
|------------------|
| 201 2 10 2 |
| 101 3 13 2 |
+------------------+

此时，两组的均值分比为 4.5 和 11.5，差为 7。

留给诸位 2 分钟，对比一下两组结果，然后想想为什么会出现结果不一致的现象？

2. 揭秘：group() 函数的工作原理

细心的读者已经发现了结果发生漂移的原因：id=201 和 id=401 的两个观察值的 x 变量具有相同的取值，x=2。然而，他们的 y 变量取值却不同。当我们执行 sort x 命令对 x 变量进行排序时，id=201 会随机地被排在第二位或第三位。虽然这对 x 的排序结果没有任何影响，但 y 变量中各个观察值出现的顺序却存在差异。

让我们更细致地解读一下 group() 函数的工作原理。

上面的分析中我们用到了如下两条核心命令：

sort x \\ 由小到大排序
gen g = group(2) \\ 等分两组

也就是说，gen 命令下的 group() 函数按如下两个步骤工作：

Step 1： 使用 sort 命令对样本进行排序；
Step 2： 将观察值平均分配为两组 (group(2))。比如，上例中样本数 N=4，那么排序后的第 1 个和第 2 个观察值就被分到第一组 (gg=1)，第 3 和 4 个观察值被分到第二组 (gg=2)。

显然，某个观察值被分到第一组还是第二组，并不是由 group() 函数决定的，而是它之前的 sort 命令决定的。

如果诸位此时输入 help sort 命令查看其帮助文件，就会发现他有一个唯一的选项——stable（我用了十几年的 Stata，今天才发现原来 sort 命令还有选项！）

执行 sort x 时，如果不加入 stable 选项，那么就会随机出现如下两种情形中的一种：

CaseA Case_B
----------- -----------
x y | x y
1 1 | 1 1
2 10 | 2 8
2 8 | 2 10
3 13 | 3 13
----------- -----------

这就是我们此前看到的两种情形。

加入 stable 选项后，每次执行 sort x 后得到的结果都是相同的，即 Case A (它维持了原始数据中 y 变量各个观察值的相对顺序)。

其实，上述表述并不严谨，sort x, stable 其实上只是帮助我们维持在执行 sort 命令前一刻内存中数据的相对状态。然而，如果在 sort 命令之前执行过会打乱样本中各个观察值相对位置的命令，即使在 sort 命令中附加 stable 选项，仍然无法保证排序结果的稳定性，也就无法保证 group() 函数分组的唯一性。

有兴趣的读者可以执行如下命令测试一下。

clear
input x y
3 13
2 10
1 1
2 8
end
rsort // 外部命令，随机排序，使用 ssc install rsort 下载
sort x, stable
gen g = group(2)
tabstat x y, by(g) f(%3.1f)
sort g x y
list, sepby(g) noobs

当然，实证分析中很少有人会使用上例中的 rsort 命令主动打乱观察值顺序 (一个典型的特例是在 PSM 分析中使用最近邻匹配时，需要预先打乱观察值顺序)。更一般的情形是如下命令 (你在毫不自知的情况下打乱了样本顺序)：

bysort industry year: egen av_x = mean(x)
sort av_x, stable
gen g_x = group(10)
……

2.1 同类问题

其实，当分组变量本身存在较多重复值时，egen 命令提供的 pctile() 函数，以及 quantiles 命令都存在上述问题，因为背后的道理都是相同的。

下面是针对 quantiles 命令的测试代码，也会存在结果不稳定的问题。

clear
input x y
3 13
2 10
1 1
2 8
end
quantiles x, gen(gg) n(2)
tabstat x y, by(g) f(%3.1f)
sort g x y
list, sepby(g) noobs

2.2 小结

如果用于分组变量存在诸多重复值，就非常容易导致上述问题；
上例中，x 变量只有一个重复值 (数值 2 重复出现了两次)，极端状况是 x 的所有观察值都集中于某一个取值，那么此时使用 group() 函数进行分组，就相当于随机分组，可能每次的结果都会不同。有兴趣的读者可以反复运行一下如下代码：

clear
input x y
2 13
2 10
2 1
1 8
2 12
end
sort x
gen g = group(2)
tabstat x y, by(g) f(%3.1f)
sort g x y
list, sepby(g) noobs

3. 解决方法

狭义而言，这个问题无解。

实证分析中的建议就是：

如果分组变量本身有很多重复值，那么使用 group() 函数进行分组是非常糟糕的选择。Note：类别变量 (如教育水平)、计数变量 (如教育年限, 专利个数等) 都会存在很多重复值。此时最好是人为地事先设定分组界点，比如教育年限低于 12 年定义为「Low」，高于 12 年定义为「High」。
如果分组变量是连续变量，如公司规模 (ln(总资产))、负债率等，出现重复值的概率很小，此时可以使用 generate 命令下的 group() 函数，但更好的办法是使用 quantiles 命令，基于分位数进行分组。但需要注意的是，quantiles 命令在样本数太小或分组个数太小时都不适用。
更好的办法是基于分位数分组 (但分组结果不再保证每组的样本数相当近似相等)，此时可以开始用 xtile，或 egen 命令下的 pctile() 函数。

使用 xtile 命令

由于 xtile 基于分位数分组，使用 sum x, detail 命令可知，x 变量的中位数是 2，因此，xitle gg=x, n(2) 会以 x=2 为分界点 (cutpoint) 将样本分成两组，用新生成的变量 gg 加以标记。此时，分组结果并不均等，这在小样本中是非常普遍的事情，对于大样本而言，组间样本数的差异基本上可以忽略。

对于前文提到的类别变量，或续别变量，我们可以使用 xtile 命令的 cutpoints() 选项人为设定分界点，此处不再赘述。也可以使用 egen 命令中的 cut() 函数实现相似的目的，缺陷是此时无法配合使用 bysort 前缀。

clear
input id x y
101 3 13
201 2 10
301 1 1
401 2 8
end
xtile gg = x, n(2) // New
tabstat x y, by(g) f(%3.1f)
sort g x y
list, sepby(g) noobs

结果如下：

. tabstat x y, by(g) f(%3.1f)
Summary statistics: mean
by categories of: gg (2 quantiles of x)
gg | x y
---------+--------------------
1 | 1.7 6.3
2 | 3.0 13.0
---------+--------------------
Total | 2.0 8.0
------------------------------
. sort g x y
. list, sepby(g) noobs
+-------------------+
| id x y gg |
|-------------------|
| 301 1 1 1 |
| 401 2 8 1 |
| 201 2 10 1 |
|-------------------|
| 101 3 13 2 |
+-------------------+

关于我们

【Stata 连享会(公众号：StataChina)】由中山大学连玉君老师团队创办，旨在定期与大家分享 Stata 应用的各种经验和技巧。
公众号推文同步发布于 CSDN-Stata连享会、简书-Stata连享会和知乎-连玉君Stata专栏。可以在上述网站中搜索关键词Stata或Stata连享会后关注我们。
点击推文底部【阅读原文】可以查看推文中的链接并下载相关资料。
Stata连享会精彩推文1 || 精彩推文2

联系我们

欢迎赐稿： 欢迎将您的文章或笔记投稿至Stata连享会(公众号: StataChina)，我们会保留您的署名；录用稿件达五篇以上，即可免费获得 Stata 现场培训 (初级或高级选其一) 资格。
意见和资料： 欢迎您的宝贵意见，您也可以来信索取推文中提及的程序和数据。
招募英才： 欢迎加入我们的团队，一起学习 Stata。合作编辑或撰写稿件五篇以上，即可免费获得 Stata 现场培训 (初级或高级选其一) 资格。
联系邮件： StataChina@163.com

往期精彩推文

欢迎加入Stata连享会(公众号: StataChina)

宾曰语云被法学教授投诉：严重侵权，“违法犯罪”！

京东Plus的隐藏特权，很多会员都没领取，白交了会员费...

呼吁四川大学澄清：1998年1月，川大有多少个“姜涛与爱人程月玲”？

二湘：朱令去世一周年，清华学子控诉清华在朱令案中的冷血和无耻

多长高8厘米！国内知名专家首次公开“追高秘笈”！担心孩子长不高的家长速来！