[旧版]Stata: 如何检验分组回归后的组间系数差异？

2018-09-29 本文已影响543人 stata连享会

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

Stata连享会精彩推文1 || 精彩推文2

Note：该文已发表。连玉君, 2017, 如何检验分组回归后的组间系数差异?, 郑州航空工业管理学院学报 35, 97-109. PDF 原文下载

问题背景

如果两个样本组中的模型设定是相同的，则两组之间的系数大小是可以比较的，而且这种比较在多数实证分析中都是非常必要的。

举几个例子，让诸位对这类问题有点感觉：

Cleary, S., 1999, The relationship between firm investment and financial status, Journal of Finance, 54 (2): 673-692. Tabel IV

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

下面使用我在stata初级班讲座 (连玉君课程）中的例子，列举几种方法。

调入 stata 自带的数据集 nlsw88.dta。

这份数据包含了 1988 年采集的 2246 个妇女的资料，包括：小时工资 wage，每周工作时数 hours，种族 race 等变量。
我们想研究的是妇女的工资决定因素是否存种族差异。
最为关注的是白人和黑人（相当于把原始数据分成了两个样本组：白人组和黑人组）的工资决定因素是否存在差异。
分析的重点集中于工龄（ttl_exp）和婚姻状况(married) 这两个变量的系数在两组之间是否存在显著差异。

下面是分组执行 OLS 回归的命令和结果：

  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)

结果：

------------------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)   
south              -0.813***       -2.041***
                  (-2.71)         (-4.92)   
hours               0.051***        0.037   
                   (3.81)          (1.39)   
tenure              0.025          -0.004   
                   (0.77)         (-0.09)   
age                 0.042           0.995   
                   (0.03)          (0.54)   
agesq              -0.001          -0.015   
                  (-0.09)         (-0.66)   
_cons               3.333         -14.098   
                   (0.14)         (-0.39)   
--------------------------------------------
N                1615.000         572.000   
r2_a                0.112           0.165   
--------------------------------------------
t statistics in parentheses
* p<0.1, ** p<0.05, *** p<0.01

可以看到，ttl_exp 变量在 [white 组] 和 [black 组] 的系数分别为 0.251 和 0.269，二者都在 1% 水平上显著异于零。

问题在于：我们能说 0.269 比 0.251 大吗？

从统计意义上来看，答案显然没有那么明确（小学五年级的小朋友会觉得这根本不是个问题！）。

相对而言，若把注意力放在 married 这个变量上，或许更容易判断二者的差异是否显著。因为，_b[married]_white (白人组的 married 估计系数) 为 -0.737**，而 _b[married]_black 为 0.091 —— 前者在 5% 水平上显著为负，而后者不显著。

即便如此，我们仍然无法直接作出结论：_b[married]_white < _b[married]_black，因为二者的置信区间尚有重叠：

      *----------------------------------------
      *             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] 
      *----------------------------------------

下面我们介绍三种检验组间系数差异的方法：

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

方法 1: 引入交叉项

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   //全样本回归+交乘

     Source |     SS       df       MS        Number of obs =    2187
------------+----------------------------     F( 20,  2166) =   17.32
      Model | 10074.761    20  503.738052     Prob > F      =  0.0000
   Residual |63000.2591  2166  29.0859922     R-squared     =  0.1379
------------+----------------------------     Adj R-squared =  0.1299
      Total |73075.0201  2186   33.428646     Root MSE      =  5.3931
----------------------------------------------------------------------
       wage |    Coef.   Std. Err.    t   P>|t|   [95% Conf. Interval]
------------+---------------------------------------------------------
      black | -.818647   .8015272  -1.02  0.307   -2.39049    .7531957
ttl_x_black |-.0181844   .0585517  -0.31  0.756  -.1330077    .0966389
    ttl_exp | .2537358   .0351178   7.23  0.000   .1848676     .322604
    married |-.4646649   .2530829  -1.84  0.066  -.9609756    .0316458
      south |-1.127138   .2453123  -4.59  0.000   -1.60821   -.6460662
      hours | .0516672   .0116886   4.42  0.000   .0287451    .0745894
     tenure | .0198005   .0260971   0.76  0.448  -.0313775    .0709786
        age | .1685498   1.035721   0.16  0.871  -1.862562    2.199661
      agesq |-.0034668   .0131097  -0.26  0.791  -.0291756     .022242
                           ...
      _cons | .8448605   20.39973   0.04  0.967  -39.16022    40.84994
----------------------------------------------------------------------

可以看出，交乘项 [ttl_x_black] 的系数为 −0.01818, 对应的 p-value 为 0.756，表明 [ttl_exp] 的系数在两组之间并不存在显著差异。

我们也可以不事先生成交乘项，而直接采用 stata 的因子变量表达式，得到完全相同的结果：

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

或如下更为简洁的方式 (详情参见 help fvvarlist)：

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

然而，需要特别强调的是，在上述检验过程中，我们无意识中施加了一个非常严格的假设条件：只允许变量 [ttl_exp] 的系数在两组之间存在差异，而其他控制变量(如 married, south, hours 等) 的系数则不随组别发生变化。

这显然是一个非常严格的假设。因为，从 -Table 1- 的结果来看, married， south, hours 等变量在两组之间的差异都比较明显。

为此，我们放松上述假设，允许 married， south, hours 等变量在两组之间的系数存在差异：

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

在这种相对灵活的设定下，[ttl_exp] 的系数为 −0.016147，相应的 p-value = 0.787，依然不显著。

当然，我们也可以采用更为灵活的方式：允许所有的变量在两组之间都存在系数差异（注意：所有离散变量前都要加 i. 前缀，否则将被视为连续变量进行处理（对于取值为0/1的虚拟变量，可以省略前缀 i.）；连续变量则需加 c. 前缀）：

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

这其实就是大名鼎鼎的 Chow test（邹检验），可以用 howtest命令快捷地完成。

解决办法：

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

上述范例中，是以基于截面数据的 OLS 回归为例的，但这一方法也适用于其他命令，如针对面板数据的 xtreg, 针对离散数据的 logit, probit 等。

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

*-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   =      2187
------------------------------------------------------------------------
          |   Robu
          |     Coef.   Std. Err.      z    P>|z|   [95% Conf. Interval]
----------+-------------------------------------------------------------
w_mean   
  ttl_exp |  .2505707   .0362302     6.92   0.000   .1795609    .3215805
  married | -.7367238   .3486479    -2.11   0.035  -1.420061   -.0533865
    south | -.8125914   .2892359    -2.81   0.005  -1.379483   -.2456994
    hours |  .0507118   .0126268     4.02   0.000   .0259637    .0754599
   tenure |  .0246063   .0289939     0.85   0.396  -.0322207    .0814333
      age |  .0415616   1.107902     0.04   0.970  -2.129887    2.213011
    agesq | -.0014454   .0138965    -0.10   0.917   -.028682    .0257911
                       ... .
    _cons |  3.333308   22.03159     0.15   0.880  -39.84781    46.51442
----------+-------------------------------------------------------------
w_lnvar  
    _cons |    3.4561   .0971137    35.59   0.000   3.265761     3.64644
----------+-------------------------------------------------------------
b_mean   
  ttl_exp |  .2686185   .0511831     5.25   0.000   .1683015    .3689355
  married |  .0913607   .4050685     0.23   0.822  -.7025589    .8852804
    south | -2.040701   .4252889    -4.80   0.000  -2.874252    -1.20715
    hours |  .0367536   .0229076     1.60   0.109  -.0081444    .0816516
   tenure |  -.003914   .0450917    -0.09   0.931  -.0922921    .0844641
      age |  .9952064   1.814351     0.55   0.583  -2.560856    4.551269
    agesq | -.0153824   .0229045    -0.67   0.502  -.0602744    .0295096
    _cons | -14.09831   35.40206    -0.40   0.690  -83.48508    55.28846
                             ... .
----------+-------------------------------------------------------------
b_lnvar  
    _cons |  3.077588   .2090237    14.72   0.000   2.667909    3.487267
------------------------------------------------------------------------

. *-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

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。

面板数据的处理方法

2018/5/3 8:04 更新 目前，可以使用 Federico Belotti. 更新后的 suest.ado 文档替换 Stata 官方提供的 suest.ado 文档。前者支持 xtreg 命令。替换方法为：net install suest, replace。
(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

小结

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

Stata连享会精彩推文1 || 精彩推文2

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

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

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

*-数据处理
  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

-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

耗时 608 秒才完成，结果如下：

-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- 的帮助文件。

F、延伸阅读

关于这一方法的更为一般化的介绍参见 Efron, B., R. Tibshirani. An introduction to the bootstrap[M]. New York: Chapmann & Hall, 1993 （Section 15.2, pp.202）.

如下论文使用了这一方法检验了 “投资-现金流敏感性” 分析中的组间系数差异：

Cleary, S., 1999, The relationship between firm investment and financial status, Journal of Finance, 54 (2): 673-692. （pp.684-685）
连玉君, 彭方平, 苏治, 2010, 融资约束与流动性管理行为, 金融研究, (10): 158-171. （pp.164）

其它基于 Permutation test 的检验

mtest 命令基于组合检验的思想来检验两个样本组是否具有相同的分布（Stata Journal, 11-2, http://www.stata-journal.com/sjpdf.html?articlenum=st0228）；
tsrtest， mwtest 命令用于检验两个样本组的均值，中位数，方差等多个维度的差异（Stata Journal, 9-1, http://www.stata-journal.com/sjpdf.html?articlenum=st0158）
上述命令都可以使用 - findit - 命令搜索到，或直接用 - ssc install 命令名称 - 下载安装。

小结

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

方法无优劣。 无论选择哪种方法，都要预先审视一下是否符合这些检验方法的假设条件。

image

往期回顾

¤ Stata兄弟：提问时能否有点诚意？用 dataex 吧
¤ 一个博士生该掌握哪些基本工具(武器)？
¤ 第一届stata中国用户大会嘉宾PDF讲稿下载
¤ 彼此不再煎熬-如何做好毕业答辩陈述？
¤ 连玉君 Markdown 笔记
¤ Stata快捷键GIF：键盘就是你的武器
¤ 如何处理时间序列中的日期间隔 (with gaps) 问题？
¤ 使用 -import fred 命令导入联邦储备经济数据库 (FRED)
¤ Fuzzy Differences-in-Differences （模糊倍分法）

关于我们

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

Stata连享会精彩推文1 || 精彩推文2

联系我们

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

往期精彩推文

Stata连享会推文列表

Stata连享会精彩推文1 || 精彩推文2

[旧版]Stata: 如何检验分组回归后的组间系数差异？

问题背景

方法 1: 引入交叉项

解决办法：

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

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

小结

往期回顾

关于我们

联系我们

往期精彩推文

猜你喜欢

热点阅读

[旧版]Stata: 如何检验分组回归后的组间系数差异？

问题背景

方法 1: 引入交叉项

解决办法：

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

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

小结

往期回顾

关于我们

联系我们

往期精彩推文

猜你喜欢

热点阅读

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