盈余管理、过度投资怎么算？分组回归获取残差
作者：吴世飞 | 知乎 | 简书 | 码云 | github
摘要：实证分析中，一般人很少关注残差，其实，残差是个神秘的黑箱，里面蕴含着丰富的信息。本文介绍如何通过分组回归获取的残差来构造衡量“盈余管理程度”、“过度投资”以及“过度消费”等指标。Stata 中的 asreg 命令可以快捷地完成这一任务！
Stata现场研讨班即将开班（2018年1月中旬，北京）
Stata现场班-初级｜Stata现场班-高级
1. 引言——福尔摩斯与怀表
相信很多人都读过英国侦探小说家阿瑟·柯南·道尔的中篇小说《四签名》。
故事中的主角大侦探福尔摩斯的名声可谓如雷贯耳，甚至在日本漫画家青山刚昌创作的侦探漫画《名侦探柯南》中，主角柯南也将福尔摩斯视作偶像。
在《四签名》中，福尔摩斯曾经准确地通过一块怀表的特征推断出了其旧主人的性格特点和生活习惯：“他是一个放荡不羁的人……最后因为好酒而死……”。
下面是相关的证据：
证据一: 这是一块价值五十多英镑的表，本应该很精心地保护和使用，但整个表的上面有无数的伤痕，这是因为惯于把表放在有钱币、钥匙一类硬东西的衣袋里的缘故。故知其主人“放荡不羁”。证据二： 表的里盖钥匙孔附近有上千的伤痕，这是因被钥匙摩擦而造成的。一般清醒的人插钥匙，一插就能够进去，只有醉汉才会因手腕颤抖而留下这些痕迹。故知其主人“好酒”。
上面的推理充分展现了福尔摩斯敏锐的观察力与严谨的分析能力。
然而，隐含在上述推理背后的逻辑是：
掌握一般情况下本应怎样，观察实际情况下却是怎样，将二者作一比对，便可以发现隐藏在后面的真相。
比如：一般情况下，对如此贵重的怀表本应精细呵护才是，
然而，实际上却在其表面发现许多伤痕，这些伤痕必然暗藏玄机……
智慧的真谛往往会得到跨领域的呼应。福尔摩斯这一推理方法也被广泛地运用于社会科学的学术研究中。
2. 实证中的黑匣子——残差
在基于大样本的社会科学实证研究中，我们往往会构建某种模型来描述事情在一般情况下本应怎样（回归模型的拟合值），然后再看一看事情的实际情况却是怎样（被解释变量的实际观测值），而上述二者之间的差异（残差）所隐藏的信息便成为了我们需要解释的重点。
我们可以看到，在很多领域，常常采用回归得到的残差作为某些重要变量的衡量指标，而这些指标继而会作为随后分析中的被解释变量。这便是上述的福尔摩斯推理思想的一种运用。例如：
在会计和公司金融领域，研究
盈余管理或盈余质量时，便是基于“分行业-分年度”回归得到的残差作为异常盈余的指标；在
事件研究法中，我们针对每家公司的日交易资料来估计市场模型(CAPM 的一个简化版本)，并取其残差用以衡量异常收益率 (Abnormal Return, AR)；在公司
投资行为的研究中，用实际投资支出与可能影响投资行为的变量进行回归，得到的残差往往被视为非预期投资，残差为负者视为投资不足，为正者则视为过度投资；在
个体消费行为的研究中，也会采取相同的思路，估计消费率方程后，用残差来衡量异常消费。
类似的例子和应用还有很多。这里我们以公司投资行为的研究为例，介绍此类模型的基本思想，并应用Stata范例来展现这一思想的实现过程。
3. 研究实例：过度投资与投资不足
Richardson（2006，“Over-Investment of Free Cash Flow”） 对公司的投资行为进行了研究，文中构建了如下形式的投资模型：
Inv[i,t] = α + βX[i,t-1] + u[i,t]模型中的被解释变量
Inv[i,t]是第i家公司在第t年的新增投资支出，解释变量 X 主要包括：公司的增长机会、杠杆率、公司规模、公司年龄、现金存量、股票回报率、新增投资支出，以及年度和行业虚拟变量等。干扰项 u 里则包含了各种无法观测，同时有可能影响公司投资支出的因素。比如，管理者的个人风格、公司文化、行业层面或宏观层面受到的各种冲击等等。
这些解释变量决定了公司新增投资支出的正常水平，因此上述回归模型的拟合值便是对公司“预期投资支出”的衡量。
若用 Inv_fit 表示拟合值，即 Inv_fit = α^ + β^*x[i,t] (这里，α^ 和 β^ 分别表示参数 α 和 β 的估计值)，则模型残差定义为： e = Inv - Inv_fit。它反映了公司的非预期投资支出。 如果 e 为正值，表示公司倾向于过度投资；反之，则意味着存在投资不足倾向 。
4. Stata 应用
4.1 残差和拟合值的估计
若假设参数 α 和 β 在全样本中为常数 （以 Leverage 变量为例，这意味着 A 公司的负债率增加一个单位对投资的边际影响与 B 公司完全相同；或者，2009 年 Leverage 增加一个单位对投资的边际影响与 2008 年和 2010 年也米有差别，这显然是一个很严格的假设条件），则估计残差和拟合值是非常简单的事情：只需在完成回归后执行 predict 命令即可。
为便于各位读者演练，这里使用 Stata 手册中的一份范例数据来说明。该数据源于 Grunfeld & Griliches (1960，“Is aggregation necessarily bad？”)。他们用公司前期市场价值和固定资产价值两个因素解释了公司的总投资支出。样本包括了 10 家公司 1935-1954 年 20 年间的数据。
虽然他们的模型与 Richardson (2006) 在设定上存在差异，但求取模型残差的思路是相通的。
webuse grunfeld, clear // 调入数据
reg invest mvalue kstock
predict inv_fit // invest 的拟合值
predict E0, res // 残差，需要附加 residual 选项，可以简写为 res
label var E0 "E0"
list comp year inv* E0 if mod(year,5)==0, sep(4)
*-Note: invest-总投资支出；
* mvalue-前期市场价值；
* kstock-前期固定资产价值部分结果呈现如下：
+-------------------------------------------------+
| company year invest inv_fit E0 |
|-------------------------------------------------|
1. | 1 1935 317.6 313.6896 3.910378 |
6. | 1 1940 461.2 541.7413 -80.54128 |
11. | 1 1945 561.2 577.8403 -16.64025 |
16. | 1 1950 642.9 644.8065 -1.906509 |
|-------------------------------------------------|
21. | 2 1935 209.9 127.138 82.76198 |
26. | 2 1940 361.6 270.496 91.10404 |
31. | 2 1945 258.7 220.4178 38.28222 |
36. | 2 1950 418.8 233.6664 185.1336 |
|-------------------------------------------------|
41. | 3 1935 33.1 115.123 -82.02305 |
46. | 3 1940 74.4 246.7319 -172.3319 |
51. | 3 1945 93.6 263.0246 -169.4246 |
56. | 3 1950 93.5 292.7397 -199.2397 |
|-------------------------------------------------|使用 xtline 命令可以很方便地实现上述结果的可视化：
xtline E0 if comp<=6, yline(0, lc(green) lp(dash))输出图形为：
4.2 分组回归的残差
然而，假设参数 α 和 β 为常数其实是一个非常严格的设定，也缺乏合理性。
我们以 Leverage 变量为例来说明。上述假设意味着 A 公司的负债率增加一个单位对投资的边际影响与 B 公司完全相同；或者，2009 年 Leverage 增加一个单位对投资的边际影响与 2008 年和 2010 年也没有差别；制造业公司的 Leverage 对投资的边际影响与零售业或金融业也完全相同。
因此，无论是在公司金融领域还是消费领域，学者们通常会放松上述假设，比如允许不同行业的 α 和 β 可以有所差异，甚至同一个行业不同年度上的 α 和 β 也可以变化。这就需要“分行业-分年度” 进行回归，并分别计算对应的残差。
以分年度计算残差为例，我们可以用循环语句来完成上述任务：
*-分年度变参数模型
*webuse grunfeld, clear
egen t = group(year) //生成 1,2, T 年份标示变量，防止原始年份数据不连续
sum t
local T = r(max) // 最后一年
gen Et = . // 用于记录残差的变量
forvalues i=1/`T'{
qui reg invest mvalue kstock if t==`i' // 分年度回归
qui predict e_i if e(sample), res // 第 t 年的残差
qui replace Et = e_i if e(sample) // 将第 t 年的残差计入变量 E
drop e_i
}
*-与参数不变模型的对比
xtline E0 Et if comp<=6, yline(0, lc(pink*0.6) lp(dash))最后一行中，仍然使用 xtline 命令同时绘制两个时序图，以作对比：
5. 便捷的 asreg 命令
上例中，我们只在年度层面上进行了分组回归，但实际应用过程中，可能还有更复杂的需求，例如：
有时需要在二维 (分年度-分行业)，甚至多维层面上进行分组回归；这就需要编写二维嵌套循环语句。
部分细分组中的观察值个数可能很少，以至于无法执行回归；此时，Stata 可能会报错，我们需要预先删除这些细分组。
有些研究中为了充分保证参数的时变性，还会进行滚动窗口回归；此时，程序会变得很复杂。
值得庆幸的是，借助 Stata 外部命令 asreg，我们可以很方便地实现上述需求。
asreg命令可以通过三种方式对样本进行分组并分别执行线性回归，最后以生成新变量的形式存储各组回归的对应统计量。
其中，三种分组方式分别为：
滚动窗口分组（rolling window）
递归窗口分组（recursive window）
一般分组
前两种分组方式的含义可参见 Stata 命令rolling的帮助文件（help rolling）,我们在此处主要运用的是第三种分组方式，即一般分组。
5.1 安装 asreg 命令
可以在 Stata 的命令窗口输入如下语句：
net install asreg, replace5.2 asreg 命令的基本语法
asreg depvar indepvars [if] [in] [,
window([rangevar] # )
recursive
minimum( # )
by(varlist)
statistics_options]window选项:设定滚动窗口或递归窗口分组中的窗口长度。recursive选项：指定采用递归窗口分组进行回归（否则，在设定了window选项的情况下，默认进行滚动窗口分组的回归）。by选项：指定进行一般分组的组别变量。minimum(#)选项：指定在分组回归中，执行每组回归所需的最小观察值数量。当设定了 min(#) 选项时，回归中使用的观察值个数将是模型中解释变量个数（包括常数项）与 # 中的较大者。如果某一组中包含的观察值数小于这里指定的数量，则与其对应的新生成变量数值为缺失值。
5.3 asreg 生成的新变量
asreg 命令会自动生成一系列 _ 开头的新变量，用于存储分组回归得到的统计量，相关统计量内容及其对应变量的命名规则如下：
| 统计量 | 命名规则 |
|---|---|
| 观察值数 | 存储各组回归观察值数量的变量为： |
| 回归系数 | 存储回归系数的变量名称将以原解释变量名称加前缀 “_b_” 组成。 |
| 常数项 | 存储常数项的变量名称为： _b_cons |
| R2 | R2和调整R2将分别被存储在名为 _R2 和 _adjR2 的变量中。 |
| 系数标准误 | 存储回归系数标准误的变量名称将以原解释变量名称加前缀 “_se_” 组成。 |
| 残差 | 存储残差的变量为： _residuals |
| 拟合值 | 存储拟合值的变量为： _fitted |
注：如果不进行额外设定，asreg命令默认存储的统计量包括：观察值数、回归系数、常数项、R2 和调整 R2 。因此，若要获得系数标准误、残差和拟合值，则需要分别添加se和fit选项(fit选项同时获取残差和拟合值)。
5.4 asreg 命令应用之 Stata 范例
掌握了这一“武器”的用法后，我们通过两个例子来实战操作一下 asreg 的用法。
范例 1：过度投资与投资不足
承接上例，我们先用 asreg 来实现分年度计算残差。上例中繁杂的循环语句此时简化为一条命令：
*- 验证 asreg 命令 -- 按年度分组回归取残差
asreg invest mvalue kstock, by(year) fit我们将命令自动生成的残差变量另存一份为 Et_as，并与此前手动编写代码得到的 Et 进行对比 —— 完全一样! 所以，日后可以放心地使用 asreg 了。
. gen Et_as = _residuals
*-对比
. list comp year Et* if mod(year,5)==0, sep(4)
+----------------------------------------+
| company year Et Et_as |
|----------------------------------------|
1. | 1 1935 1.709904 1.709904 |
6. | 1 1940 3.309084 3.309084 |
11. | 1 1945 33.3144 33.3144 |
16. | 1 1950 17.55818 17.55818 |
|----------------------------------------|
21. | 2 1935 70.00819 70.00819 |
26. | 2 1940 127.0704 127.0704 |
31. | 2 1945 57.78962 57.78962 |
36. | 2 1950 143.2945 143.2945 |
|----------------------------------------|
41. | 3 1935 -87.04494 -87.04494 |
46. | 3 1940 -139.7101 -139.7101 |
51. | 3 1945 -129.7079 -129.7079 |
56. | 3 1950 -163.8177 -163.8177 |
|----------------------------------------|前面已经提到，asreg 可以自动生成多个统计量，我们来查看一下：
format _* %4.3f
format _Nobs %2.0f
list comp year _Nobs _adjR2 /// // 样本数，R2-adj
_b_mvalue _b_kstock /// // 变量的系数估计值
_fitted _residuals /// //拟合值和残差
if mod(year,5)==0, sep(4) noobs结果如下：
+----------------------------------------------------------------------------+
| company year _Nobs _adjR2 _b_mva~e _b_kst~k _fitted _resid~s |
|----------------------------------------------------------------------------|
| 1 1935 10 0.827 0.102 -0.002 315.890 1.710 |
| 1 1940 10 0.793 0.095 0.202 457.891 3.309 |
| 1 1945 10 0.880 0.108 0.050 527.886 33.314 |
| 1 1950 10 0.817 0.176 -0.022 625.342 17.558 |
|----------------------------------------------------------------------------|
| 2 1935 10 0.827 0.102 -0.002 139.892 70.008 |
| 2 1940 10 0.793 0.095 0.202 234.530 127.070 |
| 2 1945 10 0.880 0.108 0.050 200.910 57.790 |
| 2 1950 10 0.817 0.176 -0.022 275.505 143.295 |
|----------------------------------------------------------------------------|
| 3 1935 10 0.827 0.102 -0.002 120.145 -87.045 |
| 3 1940 10 0.793 0.095 0.202 214.110 -139.710 |
| 3 1945 10 0.880 0.108 0.050 223.308 -129.708 |
| 3 1950 10 0.817 0.176 -0.022 257.318 -163.818 |
|----------------------------------------------------------------------------|minimum(#)选项的使用
上例中的数据是平行面板，每家公司均有 20 年的观察值，对于一个只有三个未知参数的回归模型而言，每各细分组都有足够的样本数。但有些情况下，个别细分组中的样本数很少，此时可以用 min() 选项预先删除这些细分组。
我们首先随机删除一些观察值，虚构一份非平行面板，进而以公司为单位进行分组回归，并要求每家公司至少要有 10 年的数据。
. webuse grunfeld, clear
. set seed 13579 //设定种子值，保证结果可重现
. sample 50 //随机抽取 50% 的观察值
. xtdes
. bysort comp: gen Ni = _N //每家公司的年数
. tab comp, sort
company | Freq. Percent Cum.
------------+-----------------------------------
4 | 14 14.00 14.00
9 | 12 12.00 26.00
10 | 12 12.00 38.00
5 | 11 11.00 49.00
6 | 10 10.00 59.00
7 | 9 9.00 68.00
8 | 9 9.00 77.00
1 | 8 8.00 85.00
3 | 8 8.00 93.00
2 | 7 7.00 100.00
------------+-----------------------------------
Total | 100 100.00可以看出，有 5 家公司的样本数都不足 10 年，它们在后续分组回归中将被忽略——通过设定 min(10) 选项来实现：
. asreg invest mvalue kstock, by(comp) min(10) fit
. format _res %4.2f
. list comp year Ni _res if mod(year,4)==0, sepby(comp)
+--------------------------------+
| company year Ni _resid~s |
|--------------------------------|
1. | 1 1944 8 . |
3. | 1 1936 8 . |
4. | 1 1940 8 . |
5. | 1 1948 8 . |
|--------------------------------|
11. | 2 1948 7 . |
15. | 2 1940 7 . |
|--------------------------------|
16. | 3 1944 8 . |
18. | 3 1936 8 . |
19. | 3 1952 8 . |
21. | 3 1948 8 . |
|--------------------------------|
26. | 4 1940 14 -4.08 |
28. | 4 1952 14 -0.36 |
30. | 4 1944 14 -9.04 |
33. | 4 1936 14 7.28 |
34. | 4 1948 14 5.11 |
|--------------------------------|
40. | 5 1944 11 -1.13 |
41. | 5 1936 11 5.37 |
|--------------------------------|
51. | 6 1940 10 -6.31 |
53. | 6 1948 10 7.21 |
57. | 6 1936 10 5.76 |
|--------------------------------|
66. | 7 1944 9 . |
67. | 7 1936 9 . |
|--------------------------------|
68. | 8 1940 9 . |
72. | 8 1944 9 . |
|--------------------------------|
77. | 9 1944 12 16.23 |
86. | 9 1948 12 -4.23 |
|--------------------------------|
90. | 10 1936 12 -0.23 |
97. | 10 1952 12 0.68 |
98. | 10 1944 12 -0.39 |
+--------------------------------+得到了各公司的回归残差之后，残差为正的公司可以视为存在过度投资，残差为负的公司可以视为存在投资不足。我们将各公司的残差状况分年度绘制出来。
webuse grunfeld, clear
asreg invest mvalue kstock, fit by(company)
xtline _residuals, yline(0,lpattern(dot))上图中，残差位于虚线上方的为过度投资部分，下方的为投资不足部分。
可以发现，大部分公司的回归残差是很接近 0 的，说明这些公司的投资水平是符合预期的。
只有前 3 家公司的回归残差出现了明显偏离 0 的现象，表明这些公司存在过度投资或投资不足。为了揭示哪些因素导致了这些公司的投资偏差，接下来要做的就是把这些回归残差作为被解释变量，去探寻其他因素对其产生的影响。
范例 2 ：超额工资
在上面的例子中，我们是在一个层面上（即，公司层面）进行分组并计算残差的。下面看一个在两个层面上分组计算残差的例子。
我们考察工资决定因素模型，并在种族-职业两个层面上分组计算超额工资（可正可负）。
sysuse "nlsw88.dta", clear // 调入数据
asreg wage age hours tenure collgrad married south, fit by(race occupation) // 在种族、职业两个维度上进行分组回归，并分别求取拟合值与残差得到超额工资的数据后，我们看一下超额工资在各分组层面上所表现出的特征。
graph hbox _residuals, over(occupation,) nooutsides //按职业绘制箱形图
从上图可见，在各种职业中，工资水平较为稳定的是交通运输（Transport）和服务业（Service）的从业人员，其回归残差均值很接近 0 ，且方差也很小。
工资水平内部差距较大的是专业技术人员（Professional/technical）、高管（ Managers/admin），以及办公人员或对专业技不作要求的职业人员（Clerical/unskilled）。在这些职业中，工资的方差很大，说明其中一部分人获得了明显的超额工资。
（注：在某些职业组别中没有数据，是由于在对应组中观察值数小于解释变量数，因而无法进行回归估计造成的。）
twoway (kdensity _residuals if race==1,color(red) legend(label(1 "white"))) /// //按种族绘制密度函数图
(kdensity _residuals if race==2,color(orange) legend(label(2 "black"))) ///
(kdensity _residuals if race==3,color(blue) legend(label(3 "other"))) ///
, legend(col(1))
从上图可见，在白种人和黑种人中，有部分人获得了明显的超额工资，而其他人种的工资水平整体比较正常，残差集中在 0 附近。
此外，在获得了超额工资的人中，获得超高工资（残差为正）的金额明显大于获得超低工资（残差为负）的金额，这可能是受到了“最低工资法”的影响，使得工资水平存在一个法定下限。
最后，再来看一下加入工会和未加入工会人员的工资差异。
twoway (histogram _residuals if union==0, legend(label(1 "nonunion"))) ///
(kdensity _residuals if union==1, legend(label(2 "union")))
从上图可见，加入工会人员（union）整体上获得了明显更多的超高工资和更少的超低工资，这是工会组织议价能力的一个体现。
6. 结语
正如本文引言部分所述，在大样本的实证研究中，对模型残差的考察已经被应用在越来越多的领域。
这其中蕴含的逻辑就是，通过对比正常情况与事实情况的差异而发现其中的异常情况。
这种方法无疑为判断某些因素对研究对象的影响提供了极大的帮助。但是，这种方法的运用也是需要一定前提的，那就是首先需要对正常情况具有充足的认识。
比如，在引言的例子中，福尔摩斯正确地认识到，一般清醒的人插钥匙，一插就能够进去这一事实，继而才能够准确判断怀表的主人是一个好酒之人。而如果福尔摩斯认为大多数人插钥匙均存在插不准的可能，则其不可能推断出怀表的主人与旁人有何不同。
甚至，如果福尔摩斯认为一般人插钥匙均不会用眼睛注视钥匙孔，因而钥匙孔附近应有更多划痕，那么他可能错误地推断该怀表的主人是一个习惯于注视钥匙孔的谨慎之人。
因此，在应用上述方法进行实证研究时，一定要对模型进行合理、完整的设定，使得模型能够充分地对研究对象作出解释，从而如实地反映正常情况。这需要我们对前期的研究成果具有充分的掌握，并在借鉴前期研究的基础上，加入自己的深思熟虑，这样才能使得模型残差确实反映了异常情况，进而保证我们从中得出准确的推断。
正如福尔摩斯所言：“在你得到所有证据之前就进行推理是个致命的错误，这会使结果带有偏见。”引述于此，与君共勉！
阅读资料
盈余管理
Jones, J. J. Earnings management during import relief investigations [J]. Journal of Accounting Research, 1991, 29 (2): 193-228.
Dechow, P. M., Sloan, R. G., Sweeney, A. P. Detecting earnings management [J]. The Accounting Review, 1995, 70 (2): 193-225.
Dechow, P. M., Hutton, A. P., Kim, J. H., Sloan, R. G. Detecting Earnings Management: A New Approach [J]. Journal of Accounting Research, 2012, 50 (2): 275-334.
Roychowdhury, S. Earnings management through real activities manipulation [J]. Journal of Accounting and Economics, 2006, 42 (3): 335-370.
事件研究
Binder, J. The event study methodology since 1969 [J]. Review of Quantitative Finance and Accounting, 1998, 11 (2): 111-137.
Slovin M B, Sushka M E. Ownership Concentration, Corporate Control Activity, and Firm Value: Evidence from the Death of Inside Blockholders[J]. Journal of Finance, 1993, 48(4):1293-1321.
张鸣, 郭思永. 大股东控制下的定向增发和财富转移——来自中国上市公司的经验证据[J]. 会计研究, 2009(5):80-88+99.
投资行为
Richardson, S. Over-investment of free cash flow [J]. Review of Accounting Studies, 2006, 11 (2): 159-189.
张功富, 宋献中. 我国上市公司投资:过度还是不足?——基于沪深工业类上市公司非效率投资的实证度量[J]. 会计研究, 2009(5):71-79+99.
杨兴全, 张照南, 吴昊旻. 治理环境、超额持有现金与过度投资——基于我国上市公司面板数据的分析[J]. 南开管理评论, 2010, 13(5):61-69.
消费行为
叶德珠, 连玉君, 黄有光, 李东辉. 消费文化、认知偏差与消费行为偏差 [J]. 经济研究, 2012, (2): 80-92.
Stata 寒假班 报名中……
连玉君主讲，2018年1月13日-21日，北京）
Stata初级班 ｜ Stata高级班 | Stata全程班
关于我们
【Stata 连享会(公众号：StataChina)】由中山大学连玉君老师团队创办，旨在定期与大家分享 Stata 应用的各种经验和技巧。
公众号推文同步发布于 【简书-Stata连享会】 和 【知乎-连玉君Stata专栏】。可以在简书和知乎中搜索关键词
Stata或Stata连享会后关注我们。推文中的相关数据和程序，以及 Markdown 格式原文 可以在 【Stata连享会-码云】 中获取。【Stata连享会-码云】 中还放置了诸多 Stata 资源和程序。如 Stata命令导航 || stata-fundamentals || Propensity-score-matching-in-stata || Stata-Training 等。
联系我们
欢迎赐稿： 欢迎将您的文章或笔记投稿至
Stata连享会(公众号: StataChina)，我们会保留您的署名；录用稿件达五篇以上，即可免费获得 Stata 现场培训 (初级或高级选其一) 资格。意见和资料： 欢迎您的宝贵意见，您也可以来信索取推文中提及的程序和数据。
招募英才： 欢迎加入我们的团队，一起学习 Stata。合作编辑或撰写稿件五篇以上，即可免费获得 Stata 现场培训 (初级或高级选其一) 资格。
联系邮件： StataChina@163.com