# Stata: 如何检验分组回归后的组间系数差异？

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

## 1. 问题背景

### 1.2 请问：2 > 1 吗？

``````. sysuse "nlsw88.dta", clear
. gen agesq = age*age
*-分组虚拟变量
. drop if race==3
. gen black = 2.race
. tab black
*-删除缺漏值
. global xx "ttl_exp married south hours tenure age* i.industry"
. reg wage \$xx i.race
. keep if e(sample)
*-分组回归
. global xx "ttl_exp married south hours tenure age* i.industry"
. reg wage \$xx if black==0
. est store White
. reg wage \$xx if black==1
. est store Black
*-结果对比
. local m "White Black"
. esttab `m', mtitle(`m') b(%6.3f) nogap drop(*.industry) ///
s(N r2_a) star(* 0.1 ** 0.05 *** 0.01)
``````

Table1: 白人组和黑人组工资影响因素差异对比

``````------------Table 1-------------------
(1)             (2)
White           Black
--------------------------------------
ttl_exp       0.251***     0.269***
(6.47)          (4.77)
married      -0.737**      0.091
(-2.31)          (0.23)
... (output ommited) ...
--------------------------------------
N          1615.000         572.000
r2_a          0.112           0.165
--------------------------------------
t statistics in parentheses
* p<0.1, ** p<0.05, *** p<0.01
``````

``````----------------------------------------
White         Black
----------------------------------------
ttl_exp
---------
beta     0.251***    0.269***
95% CI  [0.17, 0.33]   [0.16, 0.38]
----------------------------------------
married
---------
beta     -0.737**      0.091
95% CI  [-1.36, -0.11]  [-0.69, 0.87]
----------------------------------------
``````

``````. global xline "xline(0,lp(dash) lc(red*0.5))"
. coefplot White Black, keep(ttl_exp married) ///
nolabels \$xline ciopt(recast(rcap))
. graph export "Fig01.png", replace
``````

Stata连享会 专题系列 || 精彩推文

## 2. 组间系数差异的检验方法

• 方法 1： 引入交叉项（Chow 检验）
• 方法 2： 基于似无相关模型的检验方法 (suest)
• 方法 3： 费舍尔组合检验（Permutation test）

### 2.1 方法 1: 引入交叉项

#### B. Stata 实现

``````dropvars ttl_x_black marr_x_black
global xx "ttl_exp married south hours tenure age* i.industry" //Controls
gen ttl_x_black = ttl_exp*black  //交乘项
reg wage black ttl_x_black \$xx   //全样本回归+交乘项
``````

``````. reg wage black ttl_x_black \$xx   //全样本回归+交乘
Number of obs =  2,187
R-squared     = 0.1379
--------------------------------------------------
wage |      Coef.      SE     t    P>|t|
-------------+------------------------------------
black |     -0.819    0.802  -1.02  0.307
ttl_x_black |     -0.018    0.059  -0.31  0.756
ttl_exp |      0.254    0.035   7.23  0.000
married |     -0.465    0.253  -1.84  0.066
... (output omitted) ...
--------------------------------------------------
``````

``````reg wage i.black ttl_exp i.black#c.ttl_exp \$xx
``````

``````reg wage i.black##c.ttl_exp \$xx
``````

``````. reg wage i.black i.black#(c.ttl_exp i.married i.south c.hours) \$xx    //Model 2
``````

``````. global xx "c.ttl_exp married south c.hours c.tenure c.(age*) i.industry"
. reg wage i.black##(\$xx)                                               //Model 3
``````

#### C. 引入交叉项法的假设条件

• A1: 所有控制变量的系数在两组之间无差异，即 ${\gamma }^{White}={\gamma }^{Black}$;
• A2: 两组的干扰项具有相同的分布（因为估计时是将两组样本混合在一起进行估计的），即 ${ϵ}^{Wight}\sim N\left(0,{\sigma }^{2}\right)$ , 且 ${ϵ}^{Black}\sim N\left(0,{\sigma }^{2}\right)$，换言之，${\sigma }_{White}^{2}={\sigma }_{Black}^{2}$

#### D. 解决办法

• 对于 A1，实际操作过程中，可以通过引入更多的交乘项来放松 A1，如上文提到的 Model 2 或 Model 3。
• 对于 A2， 则可以在上述回归分析过程中加入 `vce(robust)` 选项，以便允许干扰项存在异方差；或加入 `vce(cluster varname)` 以便得到聚类调整后的稳健型标准误。当然，也可以在模型中允许二维聚类标准误，此时可以使用 `vce2way` 等命令。

### 2.2 方法 2: SUEST (基于似无相关模型 SUR 的检验)

#### A. 基本思想

• 其一，在估计过程中，并未预先限定白人组和黑人组各个变量的系数一定要相同，因此在 (2) 式中，我们分别用 ${\beta }_{1}$ 和 ${\beta }_{2}$ 表示白人组和黑人组各个变量的系数向量；
• 其二，两个组的干扰项可以有不同的分布，即 ${ϵ}_{1}\sim N\left(0,{\sigma }_{1}^{2}\right)$${ϵ}_{2}\sim N\left(0,{\sigma }_{2}^{2}\right)$ ，且允许二者的干扰项相关，$corr\left({ϵ}_{1},{ϵ}_{2}\right)\ne 0$

#### B. Stata 实现方法

##### B1. 基于 `suest` 命令的检验

• Step 1： 分别针对白人组和黑人组进行估计（不限于 OLS 估计，可以执行 Logit, Tobit 等估计），存储估计结果；
• Step 2： 使用 `suest` 命令执行 SUR 估计；
• Step 3： 使用 `test` 命令检验组间系数差异。

``````*-Step1: 分别针对两个样本组执行估计
reg wage \$xx if black==0
est store w  //white
reg wage \$xx if black==1
est store b  //black
*-Step 2: SUR
suest w b
*-Step 3: 检验系数差异
test [w_mean]ttl_exp = [b_mean]ttl_exp
test [w_mean]married = [b_mean]married
test [w_mean]south = [b_mean]south
``````

Step 2 的结果如下（为便于阅读，部分变量的系数未呈现）：

``````.  suest w b

Simultaneous results for w, b
Number of obs  = 2,187
---------------------------------------------------
|             Robus
|    Coef.   Std. Err.      z    P>|z|
------------+--------------------------------------
w_mean
ttl_exp |    0.251      0.036     6.92   0.000
married |   -0.737      0.349    -2.11   0.035
south |   -0.813      0.289    -2.81   0.005
... (output omitted) ..
------------+--------------------------------------
b_mean
ttl_exp |    0.269      0.051     5.25   0.000
married |    0.091      0.405     0.23   0.822
south |   -2.041      0.425    -4.80   0.000
... (output omitted) ..
------------------+--------------------------------
``````

• 白人组和黑人组的估计结果分别存储于 `w``b` 两个临时性文件中；
• 执行 `suest w b` 命令时，白人组和黑人组的被视为两个方程，即文的 (2a) 和 (2b) 式。Stata 会自动将两个方程对应的样本联合起来，采用 GLS 执行似无相关估计（SUR）；
• 由于 SUR 属于多方程模型，因此需要指定每个方程的名称，在下面呈现的回归结果中，`[w_mean]``[b_mean]` 分别是白人组和黑人组各自对应的方程名称。因此，`[w_mean]ttl_exp` 表示白人组方程中 `ttl_exp` 变量的系数，而 `[b_mean]ttl_exp` 则表示黑人组中 `ttl_exp` 变量的系数。

``````. *-Step 3: 检验系数差异

. test [w_mean]ttl_exp = [b_mean]ttl_exp
(1)  [w_mean]ttl_exp - [b_mean]ttl_exp = 0
chi2(  1) =    0.08
Prob > chi2 =    0.7735

. test [w_mean]married = [b_mean]married
(2)  [w_mean]married - [b_mean]married = 0
chi2(  1) =    2.40
Prob > chi2 =    0.1213

. test [w_mean]south = [b_mean]south
(3)  [w_mean]south - [b_mean]south = 0
chi2(  1) =    5.70
Prob > chi2 =    0.0169
``````

##### B2. 基于 `bdiff` 命令的检验

``````preserve
drop if industry==2  // 白人组中没有处于 Mining (industry=2) 的观察值
tab industry, gen(d)  //手动生成行业虚拟变量
local dumind "d2 d3 d4 d5 d6 d7 d8 d9 d10 d11" //行业虚拟变量
global xx "c.ttl_exp married south c.hours c.tenure c.age c.agesq `dumind'"
bdiff, group(black) model(reg wage \$xx) surtest
restore
``````

``````* -SUR- Test of Group (black 0 v.s 1) coeficients difference

Variables |      b0-b1    Chi2     p-value
-------------+-------------------------------
ttl_exp |     -0.017    0.07       0.788
married |     -0.814    2.32       0.128
south |      1.238    5.80       0.016
hours |      0.014    0.28       0.597
tenure |      0.030    0.32       0.571
age |     -1.027    0.23       0.629
agesq |      0.015    0.31       0.578
d2 |     -2.732    1.63       0.202
d3 |     -1.355    1.45       0.228
d4 |     -2.708    2.23       0.135
d5 |     -1.227    1.00       0.317
d6 |      0.087    0.00       0.950
d7 |     -0.534    0.07       0.785
d8 |     -1.316    1.26       0.261
d9 |      0.346    0.06       0.807
d10 |     -1.105    0.94       0.333
d11 |     -1.689    1.81       0.179
_cons |     18.770    0.20       0.652
---------------------------------------------
``````

• 使用 `suest` 时，允许两个样本组的解释变量个数不同。但由于一些技术上的问题尚未解决(很快可以解决掉)，`bdiff` 命令要求两个样本组中的解释变量个数相同。在上例中，白人组在 Mining 行业的观察值个数为零（输入`tab industry black` 可以查看），导致我们加入行业虚拟变量时，白人组只有 10 个行业虚拟变量，而黑人组则有 11 个行业虚拟变量。为此，在上述命令中，我使用 `drop if industry==2` 命令删除了 Mining 行业的观察值。
• 目前，`bdiff` 还不能很好地支持因子变量的写法 (`help fvvarlist`)，因此上例中的行业虚拟变量不能通配符方式写成 `d*`，而必须写成原始模样： `d2 d3 d4 d5 d6 d7 d8 d9 d10 d11`

#### C. 面板数据的处理方法

• `2018/5/3 8:04 更新` 目前，可以使用 `Federico Belotti.` 更新后的 `suest.ado` 文档替换 Stata 官方提供的 `suest.ado` 文档。前者支持 `xtreg` 命令。替换方法为：`net install suest, replace`。此时，`suest.ado` 文件被自动安装在 stata15\ado\plus\s 文件夹中。一劳永逸的做法是用该文件替换掉 stata15\ado\base\s 文件夹中的同名文件。

• (`2017/9/3 8:05`) `suest` 不支持 `xtreg` 命令，因此无法直接将该方法直接应用于面板数据模型，如 FE 或 RE。此时，可以预先手动去除个体效应，继而对变换后的数据执行 OLS 估计，步骤如下：

• step 1：对于固定效应模型而言，可以使用 `center``xtdata` 命令去除个体效应；对于随机效应模型而言，可以使用 `xtdata` 命令去除个体效应。
• step 2：按照截面数据的方法对处理后的数据进行分组估计，并执行 `suest` 估计和组间系数检验。
• 举个例子：

``````*-SUEST test for panel data
*-数据概况
webuse "nlswork", clear
xtset idcode year
xtdes
*-对核心变量执行组内去心：去除个体效应
help center   //外部命令, 下载命令为 ssc install center, replace
local y "ln_wage"
local x "hours tenure ttl_exp south"
bysort id: center `y', prefix(cy_)   //组内去心
bysort id: center `x', prefix(cx_)
*-分组回归分析
reg cy_* cx_* i.year if collgrad==0  // 非大学组
est store Yes
reg cy_* cx_* i.year if collgrad==1  //   大学组
est store No
*-列示分组估计结果
esttab Yes No, nogap mtitle(Yes_Coll No_Coll) ///
star(* 0.1 ** 0.05 *** 0.01) s(r2 N)
*-似无相关估计
suest Yes No
*-组间差异检验
test [Yes_mean]cx_ttl_exp = [No_mean]cx_ttl_exp
test [Yes_mean]cx_hours = [No_mean]cx_hours
``````

#### D. 小结

• 相对于方法1（引入交乘项），基于 SUR 的方法更为灵活一些。在上例中，白人组和黑人组的被解释变量相同 （均为 wage），此时方法 1 和方法 2 都能用。有些情况下，两个组中的被解释变量不同，此时方法 1 不再适用，而方法 2 则可以。
• 对于面板数据而言，可以预先使用 `center``xtdata` 命令去除个体效应，变换后的数据可以视为截面数据，使用 `regress` 命令进行估计即可。
• 为了便于呈现结果，可以使用`estadd` 命令将上述检验结果(chi2 值或 p值) 加入内存，进而使用 `esttab` 命令列示出来。可以参考 `help bdiff` 中的类似范例。

### 2.3 方法 3：费舍尔组合检验 （Fisher's Permutation test)

#### A. 基本思想

``````. dis normal(-0.018)    //单尾检验
.49281943
``````

#### B. 经验 p 值 (empirical p-value)

``````. preserve
. clear
. set obs 10000
. set seed 1357
. gen d = rnormal()  // d~N(0,1) 服从标准正态分布的随机数
. sum d, detail
. count if d<-0.018
. dis "Empirical p-value = "  4963/10000
. restore
``````

#### D. 费舍尔组合检验的步骤

• Step 0: 分别针对白人组和黑人组估计模型 (3a) 和 (3b)，得到系数估计值  ${\stackrel{^}{\beta }}_{1}$ 和 ${\stackrel{^}{\beta }}_{2}$，以及二者的系数差异 ${\stackrel{^}{d}}_{0}={\stackrel{^}{\beta }}_{1}-{\stackrel{^}{\beta }}_{2}$
• Step 1: 将白人组和黑人组的样本混合起来，得到 ${n}_{1}+{n}_{2}$ 个观察值构成的样本 $S$
• Step 2: 获得经验样本 —— 从 $S$ 中随机抽取 (无放回) ${n}_{1}$ 个观察值，将其视为“白人组”(记为 $Sw$)，剩下的 ${n}_{2}$ 观察值可以视为“黑人组” (记为 Sb)；
• Step 3: 分别针对经验样本 $Sw$ 和 $Sb$，估计模型 (3a) 和 (3b)，得到 ${\stackrel{^}{\beta }}_{1}^{{S}_{1}}$ 和 ${\stackrel{^}{\beta }}_{2}^{S1}$ （上标 ${}^{{S}_{1}}$ 表示利用第一笔经验样本得到的估计值），以及二者的差异 ${d}^{{S}_{1}}={\stackrel{^}{\beta }}_{1}^{{S}_{1}}-{\stackrel{^}{\beta }}_{2}^{{S}_{1}}$
• Step 4: 获得统计量 $d$ 的经验分布 —— 将 Step 2 和 Step 3 重复执行 $K$ 次 (如 $K=1000$)，则可以得到  ，亦可简记为 ${d}^{{S}_{j}}\phantom{\rule{thinmathspace}{0ex}}\left(j=1,2,\cdots ,K\right)$
• Step 5: 计算经验 p 值：$\stackrel{^}{p}=\frac{\mathrm{#}\left\{{d}^{{S}_{j}}>{\stackrel{^}{d}}_{0}\right\}}{K}$ 其中， $\mathrm{#}\left\{{d}^{{S}_{j}}>{\stackrel{^}{d}}_{0}\right\}$ 表示 Step 4 中得到的 $K$ 个 ${d}^{{S}_{j}}$ 中大于我们实际观测到的 ${\stackrel{^}{d}}_{0}$ 的个数。若 $\stackrel{^}{p}<0.05$ ，则可以在 5% 水平上拒绝原假设，表明两组的系数差异是显著的。

Stata连享会 专题系列 || 精彩推文

#### E. Stata 实现

``````* 数据处理
. sysuse "nlsw88.dta", clear
. gen agesq = age*age
. drop if race==3
. gen black = 2.race
. global xx "ttl_exp married south hours tenure age agesq"
*-检验
. bdiff, group(black) model(reg wage \$xx) reps(1000) detail
``````

• 选项 `group()` 中填写用于区分组别的类别变量（若有多个组，可以预先删除不参与比较的组，类似于上面的 `drop if race==3` 命令）；
• 选项 `model()` 用于设定回归模型，即上面提到的模型 (3a) 或 (3b)；
• 选项 `reps(#)` 用于设定抽样次数，即上文提到的 K ，通常设定为 1000-5000 次即可；
• 附加 `detail` 选项，可以进一步列表呈现两组的实际估计系数 ${\stackrel{^}{\beta }}_{1}$ 和 ${\stackrel{^}{\beta }}_{2}$

``````- Permutaion (1000 times)- Test of Group (black 0 v.s 1) coeficients difference

Variables |      b0-b1    Freq     p-value
-------------+-------------------------------
ttl_exp |      0.007     490       0.490
married |     -0.824     920       0.080
south |      1.411       5       0.005
hours |      0.010     344       0.344
tenure |     -0.006     512       0.488
age |     -1.579     751       0.249
agesq |      0.022     218       0.218
_cons |     28.051     267       0.267
---------------------------------------------
``````

``````*-数据预处理(可以忽略)
sysuse "nlsw88.dta", clear
gen agesq = age*age
*-分组虚拟变量
drop if race==3
gen black = 2.race
*-删除缺漏值
global xx "ttl_exp married south hours tenure age* i.industry"
qui reg wage \$xx i.race
keep if e(sample)
*-生成行业虚拟变量
drop if industry==2  // 白人组中没有处于 Mining (industry=2) 的观察值
tab industry, gen(d)  //手动生成行业虚拟变量
local dumind "d2 d3 d4 d5 d6 d7 d8 d9 d10 d11" //行业虚拟变量
global xx "c.ttl_exp married south c.hours c.tenure c.age c.agesq `dumind'"
*-permutation test
bdiff, group(black) model(reg wage \$xx) reps(1000) detail
``````

#### F. 面板数据

``````*-Panel Data (sample by cluster(id))
*-数据预处理
. webuse "nlswork.dta", clear
. xtset id year                  //声明为面板数据，否则视为截面数据
. gen agesq = age*age
. drop if race==3
. gen black = 2.race
*-检验
. global x "ttl_exp hours tenure south age agesq"
. local m "xtreg ln_wage \$x, fe" //模型设定
. bdiff, group(black) model(`m') reps(1000) bs first detail
``````

``````-Bootstrap (1000 times)- Test of Group (black 0 v.s 1) coeficients difference

Variables |      b0-b1    Freq     p-value
-------------+-------------------------------
ttl_exp |      0.011      47       0.047
hours |      0.003      58       0.058
tenure |      0.004     116       0.116
south |      0.093      14       0.014
age |     -0.022     984       0.016
agesq |      0.000      45       0.045
_cons |      0.302      22       0.022
---------------------------------------------
Ho: b0(ttl_exp) = b1(ttl_exp)
Observed difference =  0.011
Empirical p-value =  0.047
``````

• `bdiff` 命令中设定 `bs` 选项，则随机抽样为可重复抽样 (bootstrap)；
• `first` 选项便于将组间系数差异检验结果保存在内存中，方便后续使用 `esttab` 合并到回归结果表格中。具体使用方法参见 `help bdiff`
• 由于抽样过程具有随机性，因此每次检验的结果都有微小差异。在投稿之前，可以附加 `seed()` 选项，以保证检验结果的可复制性。
• 其他选项和使用方法参阅 `help bdiff` 的帮助文件。

#### G. 延伸阅读

• Cleary, S., 1999, The relationship between firm investment and financial status, Journal of Finance, 54 (2): 673-692. (pp.684-685) -PDF-

• 连玉君, 彭方平, 苏治, 2010, 融资约束与流动性管理行为, 金融研究, (10): 158-171. （pp.164）-PDF-

### H. 其它基于 Permutation test 的检验

• `mtest` 命令基于组合检验的思想来检验两个样本组是否具有相同的分布。
• Tebaldi P, Bonetti M, Pagano M. M statistic commands: Interpoint distance distribution analysis. Stata Journal, 2011, 11(2):271-289. - PDF -
• `tsrtest`, `mwtest` 命令用于检验两个样本组的均值，中位数，方差等多个维度的差异。
• Kaiser J, Lacy M G. A general-purpose method for two-group randomization tests[J]. Stata Journal, 2009, 9(1):70-85. - PDF -
• 上述命令都可以使用 `findit` 命令搜索到，或直接用 `ssc install` 命令下载安装。

## 3. 小结

• 方法1（加入交乘项） 在多数模型中都可以使用，但要注意其背后的假设条件是比较严格的。若在混合回归中，只引入你关心的那个变量（`ttl_exp`）与分组变量 (`black`) 的交乘项（`ttl_exp*black`），则相当于假设其他控制变量在两组之间的不存在系数差异。相对保守的处理方法是：在混合估计时，引入所有变量与分组变量的交乘项，同时附加 `vce(robust)` 选项，以克服异方差的影响。
• 方法2（基于 SUR 模型的检验） 执行起来比较方便，假设条件也比较宽松：允许两组中所有变量的系数都存在差异，也允许两组的干扰项具有不同的分布，且彼此相关。局限在于，有些命令无法使用 `suest` 执行联合估计，如几乎所有针对面板数据的命令都不支持（`xtreg``xtabond` 等）。
• 方法3（组合检验） 是三种方法中假设条件最为宽松的，只要求原始样本是从母体中随机抽取的（看似简单，但很难检验，只能靠嘴说了），而对于两个样本组中干扰项的分布，以及衡量组间系数差异的统计量 $d={\beta }_{1}-{\beta }_{2}$ 的分布也未做任何限制。事实上，在获取 经验 p 值 的过程中，我们采用的是“就地取材”、“管中窥豹”的思路，并未假设 $d$ 的分布函数；另一个好处在于可以适用于各种命令，如 `regress``xtreg``logit`, `ivregress' 等，而`suest` 则只能应用于部分命令。

## 相关课程

http://lianxh.duanshu.com

### 课程一览

Stata数据清洗 游万海 直播, 2 小时，已上线

Note: 部分课程的资料，PPT 等可以前往 连享会-直播课 主页查看，下载。

#### 关于我们

• Stata连享会 由中山大学连玉君老师团队创办，定期分享实证分析经验。直播间 有很多视频课程，可以随时观看。
• 连享会-主页知乎专栏，300+ 推文，实证分析不再抓狂。
• 公众号推文分类： 计量专题 | 分类推文 | 资源工具。推文分成 内生性 | 空间计量 | 时序面板 | 结果输出 | 交乘调节 五类，主流方法介绍一目了然：DID, RDD, IV, GMM, FE, Probit 等。
• 公众号关键词搜索/回复 功能已经上线。大家可以在公众号左下角点击键盘图标，输入简要关键词，以便快速呈现历史推文，获取工具软件和数据下载。常见关键词：`课程, 直播, 视频, 客服, 模型设定, 研究设计, stata, plus, 绘图, 编程, 面板, 论文重现, 可视化, RDD, DID, PSM, 合成控制法`

✏ 连享会学习群-常见问题解答汇总：
https://gitee.com/arlionn/WD