Stata：分组回归系数比较的新思路

发布时间：2021-09-15 阅读 17134

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

温馨提示： 定期清理浏览器缓存，可以获得最佳浏览体验。

New！ lianxh 命令发布了：
随时搜索推文、Stata 资源。安装命令如下：
. ssc install lianxh
详情参见帮助文件 (有惊喜)：
. help lianxh

✌ 课程详情： https://gitee.com/lianxh/Course

⛳ 课程主页： https://gitee.com/lianxh/Course

⛳ Stata 系列推文：

全部 | Stata入门 | Stata教程 | Stata资源 | Stata命令
计量专题 | 论文写作 | 数据分享 | 专题课程
结果输出 | Stata绘图 | 数据处理 | Stata程序
回归分析 | 面板数据 | 交乘项-调节 | IV-GMM
内生性-因果推断 | 倍分法DID | 断点回归RDD | PSM-Matching | 合成控制法
Probit-Logit | 时间序列 | 空间计量 | 分位数回归 | 生存分析 | SFA-DEA
文本分析-爬虫 | Python-R-Matlab | 机器学习
Markdown | 工具软件 | 其它

作者：李适源 (北京大学)
邮箱：shiyuanli@pku.edu.cn

目录 [TOC]

1. 从 bdiff 命令的 “美中不足” 说起

在往期推文中，连老师为大家介绍了在 Stata 当中使用 bdiff 命令检验分组回归后的组间系数差异是否显著，详见推文「如何检验分组回归后的组间系数差异？」。但略有遗憾的是，在使用 bdiff 命令进行系数组间差异检验的操作过程中，大家可能会产生如下两点困惑。

1.1 两组变量对应相同的约束

困惑之一是两组变量对应相同的约束。我们需要保证参与系数比较的两组样本 (已婚组=1、单身组=0)，在分组回归中的自变量是对应相同的。在某些情形下，两组样本中自变量对应相同的要求其实很难达到。特别是，当我们在分组回归中引入了大量虚拟变量 (如地区固定效应、年份固定效应、乃至地区与年份的交叉固定效应) 的时候，“已婚组” 和 “单身组” 这两组样本在某些虚拟变量中可能并不存在观测值。

具体来说，设想我们设置了区县固定效应，在回归方程中引入了 29 个区县虚拟变量 (假设全样本中有 30 个区县)。“已婚组” 样本中可能没有受访者居住在 “区县 2 和区县 3” (这意味着，对于已婚组样本而言，区县 2 和区县 3 这两个虚拟变量对应系数无法估计)；而 “单身组” 样本中可能没有受访者居住在区县 4 和区县 5 (这意味着对于未婚组样本而言，区县 4 和区县 5 这两个虚拟变量对应系数无法估计)。

此时，我们就违反了 “两组变量对应相同” 的要求，无法直接运行 bdiff 命令来检验系数的组间差异。在上述的研究情形中，即使在每个区县中都有 “未婚组” 样本的受访者居住，也仍然无法直接运行 bdiff 命令。因为在这个时候，已婚组样本只有 27 个虚拟变量参与系数比较，而未婚组样本却有 29 个虚拟变量参与系数比较，我们仍会遭遇 “两组自变量个数不相同” 的问题。

为了更直观地了解这个问题，我们使用 Stata 的系统数据集 nlsw88 进行演示。这份数据是美国 1988 年收集的 2246 位女性的相关资料，核心变量包括：小时工资 wage，是否接受过大学教育 collgrad，每周工作时数 hours，工龄 ttl_exp，婚姻状态 married (已婚=1，单身=0)。

. sysuse nlsw88.dta,clear
. set linesize 80     //设置表格展示的格式
. set cformat  %5.4f  //回归结果中系数和标准误的显示格式
. set sformat  %4.2f  //回归结果中 t 值的显示格式
. set pformat  %4.3f  //回归结果中 p 值的显示格式
. sum wage collgrad hours ttl_exp  married //描述性统计

    Variable |        Obs        Mean    Std. dev.       Min        Max
-------------+---------------------------------------------------------
        wage |      2,246    7.766949    5.755523   1.004952   40.74659
    collgrad |      2,246    .2368655    .4252538          0          1
       hours |      2,242    37.21811    10.50914          1         80
     ttl_exp |      2,246    12.53498    4.610208   .1153846   28.88461
     married |      2,246    .6420303    .4795099          0          1

我们的研究问题是，接受大学教育对于女性工资的影响 (高等教育回报率)，是否在 “已婚组” 和 “单身组” 之间存在显著差异。首先，我们对工资取对数 lnwage，并将 lnwage 作为回归方程的结果变量，将是否接受大学教育 collgrad 作为核心解释变量，并加入工时 hours、工龄 ttl_exp 等作为控制变量，进行 OLS 分组回归。此时，我们可直接使用 bdiff 命令，选择自助法 (Bootstrap) 重复抽样的方式，检验 collgrad 对应的回归系数在 “已婚组” 和 “单身组” 之间的差异。

. gen lnwage=ln(wage)       //工资对数
. global xx "ttl_exp hours" //控制变量
. bdiff, group(married) model(reg lnwage collgrad $xx) reps(500) bsample seed(123)

  Variables |      b0-b1    Freq     p-value
-------------+-------------------------------
    collgrad |      0.037     118       0.236
     ttl_exp |     -0.003     379       0.242
     married |      0.000     232       0.464
       south |     -0.043     401       0.198
       hours |      0.004      43       0.086
         age |     -0.009     449       0.102
       _cons |      0.221     109       0.218
---------------------------------------------

从上表的第一行系数可以看到，collgrad 的分组回归系数并不存在显著差异。现在，我们引入一组 “区县虚拟变量” countyid。由于原始数据中并没有记录受访者居住的区县信息，本文出于演示的目的，使用均匀分布的随机数，生成 30 个区县代码 (1、2……29、30)。并在上述回归当中引入 29 个区县虚拟变量 i.countyid 作为控制变量 (即设置 “区县固定效应”)。

为了说明 bdiff 对分组回归变量个数需要对应相同的约束，我们人为地设定，在 “已婚样本” 中，没有受访者居住在 “区县 10、区县 20、区县 30”；同时也人为设定，在 “单身样本” 中，没有受访者居住在 “区县 1、区县 11、区县 21”。

然后，我们可以再尝试使用 bdiff 来检验 collgrad 对应系数在已婚组和未婚组的差异。容易发现，由于两组回归的自变量对应 “不一致”，运行结果将会报错 conformability error。

. set seed 10101
. gen countyid=int(31*runiform()) //使用随机数生成区县代码 (1~30)
. recode countyid 10=. 20=. 30=. if married==1
. recode countyid 11=. 21=. 29=. if married==0
. bdiff, group(married) model(reg lnwage collgrad $xx i.countyid) reps(500) bsample

conformability error
r(503);

接下来，我们还可以考察，在两组回归的变量个数不同时的运行问题。首先，drop countyid，重新用随机数生成countyid。然后，我们人为地设定 “已婚样本” 中，没有受访者居住在 “区县 10、区县 20、区县 30”，而 “单身样本” 却在每一个区县中都存在受访者。然后我们再尝试使用 bdiff 来检验 collgrad 对应系数在已婚组和单身组的差异。容易发现，由于两组回归的自变量个数 “不同”，运行结果也会报错。

. drop countyid
. set seed 10101
. gen countyid=int(31*runiform()) //使用随机数生成区县代码 (1~30)
. recode countyid 10=. 20=. 30=. if married==1
. bdiff, group(married) model(reg lnwage collgrad $xx i.countyid) reps(500) bsample

Note: The numbers of effective independent variables in Group1 (married=0) and g
> roup2 (married=1) are not equal
This may occur because you use the factor variable to indicate the dummy varaibl
> es

为了化解上述问题，我们可能需要在分组回归前，检查两组样本的自变量是否对应相等。如果某组存在 “虚拟变量缺失”，则需要在回归方程中手动删去 “缺失的虚拟变量”。这种做法在思路上容易理解，但在实际操作中却比较麻烦。特别是当方程中引入了许多高维固定效应后，可能产生大量的虚拟变量，“事前检查和清理” 的工作很可能变得繁琐许多。

1.2 非线性模型的系数比较问题

根据标准的计量理论，非线性模型 (如 Logit、Probit) 的分组回归系数不能直接比较，因为存在 “误差项的方差设定不同、系数测度标准不同、遗漏无关变量对系数测度的干扰” 等问题 (Mize, 2019；Wooldridge, 2010；洪岩璧, 2015)，并且非线性模型的回归系数一般也不能直接度量某个解释变量对结果变量的边际影响 (Marginal Effect)。学者一般建议将非线性模型的回归系数转化为平均边际效应 (Average Marginal Effect, AME)。平均边际效应又称为平均偏效应 (Average Partial Effect, APE)。

与非线性模型的回归系数不同，平均边际效应系数在不同组别间、或在不同模型设定下均具有较好的可比性。因此，如果想要比较分组回归系数的组间差异，学者往往建议比较不同组的平均边际效应 (而非回归系数) 是否存在组间差异。然而，从 bdiff 的说明性文件来看，针对非线性模型，bdiff 命令主要检验的是 probit 或 logit 回归系数的组间差异，不能直接用来检验平均边际效应系数的组间差异。

2. 新思路：从 “面面俱到” 转向 “一心一意”

应对 bdiff 的两处 “美中不足”，本文的解决思路是转换组间差异分析的重点，聚焦于某一个核心解释变量对应系数 (平均边际效应) 的组间差异。这样，在很大程度上可以克服 bdiff 对两组自变量 “对应相同” 的约束，同时也将关注重点从回归系数转向了更具可比性的平均边际效应。

诚然，上述方法并不是万全之策，我们无法像 bdiff 命令那样同时检验其他自变量对应系数的组间差异。但是，在社会科学实证研究中，我们往往只关心某一个核心解释变量 $D$ (也称原因变量或处理变量) 对于被解释变量 (结果变量) $Y$ 的边际影响，并不太关心控制变量们所对应的系数。

在组间系数比较的研究情境中，我们可能只需比较原因变量 $D$ 所对应的平均边际效应系数 $A M E (D)$ 是否在各个组别间存在显著差异即可。

2.1 边际效应 (ME) 与平均边际效应 (AME)

考虑下面这个总体回归方程，其中 $X$ 是一组控制变量，我们关心的是核心解释变量 (下面简称原因变量) $D$ 对结果变量 $Y$ 的边际效应 $M E (D)$ 。边际效应的计算方式有如下两种：如果 $D$ 是连续型变量，则用 $E (Y | D, X)$ 对 $D$ 求偏导数；如果 $D$ 是离散型变量，则对 $D$ 求差分。容易看出，在线性模型设定下 (且不存在高次项和交互项)，边际效应就是总体回归系数 $β_{1}$ 本身。

E (Y | D, X) = β_{0} + β_{1} D + X^{T} δ

M E (D) = \frac{\partial E (Y | D, X)}{\partial D}

M E (D) = E (Y | D = 1, X) - E (Y | D = 0, X)

当线性模型中引入了解释变量 $D$ 的高次项或与其他变量的交互项，或者是在非线性模型的设定下，边际效应往往 “因个体而异”，此时需要计算平均边际效应 (AME)，也即将样本中每位个体的边际效应进行加总平均。原因变量 $D$ 对于结果变量 $Y$ 的平均边际效应 (AME) 可解读为，在给定控制变量 $X$ 的情形下， $D$ 变动 (微小的) 一单位将会给条件期望 $E (Y | D, X)$ 带来的平均影响。

2.2 比较 AME 系数的组间差异

那么，应如何比较原因变量 $D$ 对应的平均边际效应系数是否在各个组别间存在显著差异呢？我们可以通过自助法 (Bootstrap) 获取平均边际效应系数之差 ${\hat{α}}_{d i f} = [A M E (G r o u p 1) - A M E (G r o u p 0)]$ 的标准误 $S E ({\hat{α}}_{d i f})$ 。然后, 构造统计量 $T = \frac{{\hat{α}}_{d i f}}{S E ({\hat{α}}_{d i f})}$ 检验原因变量 $D$ 的平均边际效应系数在两组间差异是否显著区别于 0。

传统文献中，多使用结构方程模型以及 Delta Method 来估计 ${\hat{α}}_{d i f}$ 的标准误。但相比之下，基于 Bootstrap 方法获取的标准误可能更接近经验分布。尤其是样本量偏小的时候，Bootstrap 得出的标准误往往具有更好的有限样本性质 (陈强, 2014)。

因此，本文下列内容主要介绍在常见的线性模型和非线性模型设定下，如何使用 Bootstrap 方法来检验原因变量 $D$ 对应的平均边际效应系数 (AME) 在各个组别间是否存在显著差异，主要内容如下：

截面 OLS 模型设定下，AME 系数的组间差异显著性检验；
截面 Probit 或 Logit 模型设定下，AME 系数的组间差异显著性检验；
工具变量的两阶段最小二乘 (IV-TSLS) 设定下，AME 系数的组间差异显著性检验；
面板数据的线性固定效应 (FE-OLS) 模型设定下，AME 系数的组间差异显著性检验。

3. AME 系数的组间差异检验：截面 OLS

本节选择了最为简单的模型设定，我们以工资对数 lnwage 作为结果变量，将是否接受大学教育 collgrad 作为原因变量 (核心解释变量)，以工时 hours，工龄 ttl_exp 作为一组控制变量，并设置区县固定效应 i.countyid。

我们关心的是原因变量 collgrad 对结果变量 lnwage 的平均边际效应 $A M E (D)$ ，是否在已婚组和单身组之间存在显著的组间差异。方便起见，假定在给定控制变量的情形下，原因变量 collgrad 是外生的 (总体模型的误差项 $u$ 均值独立于原因变量 collgrad)。因此，原因变量 collgrad 对应的平均边际效应 $A M E (D)$ ，可以解读为大学教育的收入回报率。

sysuse nlsw88.dta,clear          //导入示例数据
gen lnwage=ln(wage)              //生成工资对数
set seed 10101
gen countyid=int(31*runiform())  //使用随机数生成区县代码 (1~30)
recode countyid 10=. 20=. 30=. if married==1 //人为设置已婚组中缺失的虚拟变量
recode countyid 11=. 21=. 29=. if married==0 //人为设置未婚组中缺失的虚拟变量
global xx "ttl_exp hours i.countyid"        //控制变量的宏

*手动编写AME系数组间差异检验的命令
capture program drop AME_Linear      //清空历史命令
program define AME_Linear, eclass    //命令开始，将这段命令的名字定义为"AME_Linear"
  preserve
  tempname b b1 b0                   //定义三个用于临时存储的名字
	reg lnwage collgrad $xx i.countyid if married==1 //已婚组进行OLS回归
	margins,dydx(collgrad) post
	matrix `b1'=e(b)                //获取已婚组的平均边际效应b1(即已婚组的教育回报率)
	reg lnwage collgrad $xx i.countyid if married==0 //单身组进行OLS回归
	margins,dydx(collgrad) post
	matrix `b0'=e(b)                //获取单身组的平均边际效应b0 (即单身组的教育回报率)
	matrix `b'=`b1'-`b0'            //获取两组样本的平均边际效应之差 (b1-b0)
	ereturn post `b'
end

. *使用Bootstrap方法运行上述命令"AME_Linear"，重复进行500次
. bootstrap _b, nowarn reps(500) seed(10101): AME_Linear

Bootstrap results                                        Number of obs = 2,246
                                                         Replications  =   500

------------------------------------------------------------------------------
             |   Observed   Bootstrap                         Normal-based
             | coefficient  std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
    collgrad |    -0.0849     0.0598    -1.42   0.156      -0.2022      0.0323
------------------------------------------------------------------------------

. est store AME_diff_OLS //将组间差异检验的结果存储到"AME_diff_OLS"

在上面的操作步骤中，最关键的是要自己编写一个简短的命令 AME_Linear，这个命令主要做了以下三件工作：(1) 对已婚组进行 OLS 回归，并获取获取已婚组的平均边际效应 b1；(2) 对单身组进行 OLS 回归，并获取获取单身组的平均边际效应 b0；(3) 获取已婚组和单身组的平均边际效应之差 (b1-b0)。在编写完成 AME_Linear 命令之后，即可使用 bootstrap 命令，命令格式为：

bootstrap _b，reps(规定重复抽样的次数) seed(设置随机数种子): 在冒号后填写自行编写的命令-例如"AME_Linear"

bootstrap 命令将会对通过 AME_Linear 获得的平均边际效应之差 (b1-b0) 进行 500 次重复抽样，从而得到 (b1-b0) 的 Bootstrap 标准误、及其对应的 Z 值、P 值、95% 正态置信区间。我们据此可检验，collgrad 对应的 AME 系数是否在已婚组和未婚组之间存在显著差异。从 bootstrap 的运行结果看，检验统计量 Z 值为 -1.4，P 值大于 0.1，可认为 AME 系数在 10% 水平上并不存在显著差异。

为了更直观地展示 collgrad 在已婚样本、单身样本中的平均边际效应及其组间差异，我们可以使用 est store 和 esttab 命令。

. reg lnwage collgrad $xx i.countyid if married==1
. margins,dydx(collgrad) post
. est store AME_Group1_OLS     //保存已婚组的AME系数
. reg lnwage collgrad $xx i.countyid if married==0
. margins,dydx(collgrad) post
. est store AME_Group0_OLS     //保存单身组的AME系数
. esttab AME_Group1_OLS AME_Group0_OLS AME_diff_OLS,b(%9.4f) ///
>        se mtitle star(* 0.1 ** 0.05 *** 0.01) modelwidth(16)

------------------------------------------------------------------------
                          (1)                 (2)                 (3)   
               AME_Group1_OLS      AME_Group0_OLS        AME_diff_OLS   
------------------------------------------------------------------------
collgrad               0.3679***           0.4528***          -0.0849   
                     (0.0314)            (0.0479)            (0.0598)   
------------------------------------------------------------------------
N                        1325                 726                2246   
------------------------------------------------------------------------
Standard errors in parentheses
* p<0.1, ** p<0.05, *** p<0.01

4. AME 系数的组间差异检验：截面 Probit 与 Logit

根据上述思路，我们也可以在非线性模型设定下，检验 AME 系数的组间差异。出于演示的目的，我们将之前的结果变量lnwage (对数工资) 转化为二值变量 highwage (是否为高收入)。如果对数工资 lnwage 处于 75% 分位数以下，则 highwage=0；如果对数工资 lnwage 处于 75% 分位数及以上，则 highwage=1。

. sum lnwage, d   //75%分位数为2.261
. gen highwage=(lnwage>=2.261)
. tab highwage    //得到新的二值结果变量，高收入=1，不是高收入=0

   highwage |      Freq.     Percent        Cum.
------------+-----------------------------------
          0 |      1,684       74.98       74.98
          1 |        562       25.02      100.00
------------+-----------------------------------
      Total |      2,246      100.00

接下来，我们以受访者是否为高收入 highwage 作为结果变量，将是否接受大学教育 collgrad 作为原因变量 (核心解释变量)，以工时 hours，工龄 ttl_exp 作为一组控制变量，并设置区县固定效应 i.countyid。我们关心的是原因变量 collgrad 对结果变量 highwage 的平均边际效应 $A M E (D)$ ，是否在已婚组和单身组之间存在显著的组间差异。

在这里， $A M E (D)$ 可以解读为，”大学教育” (相比未接受过大学教育) 对于 “获得高收入的概率” 的平均影响，也可近似地理解为大学教育的收入回报率 (但此处的收入变量为二值变量)。由于结果变量是二值变量，我们可以选择 Probit 或 Logit 模型来估计平均边际效应 $A M E (D)$ (当然 OLS 在很多情况下也可以提供较好的线性近似)。

*手动编写AME系数组间差异检验的命令
capture program drop AME_Probit     //清空历史命令
program define AME_Probit, eclass   //命令开始，将这段命令的名字定义为"AME_Probit"
  preserve
  tempname b b1 b0                  //定义三个用于临时存储的名字
  probit highwage collgrad $xx i.countyid if married==1 //已婚组进行Probit回归
  margins,dydx(collgrad) post
  matrix `b1'=e(b)                 //获取已婚组的平均边际效应b1 (即已婚组的教育回报率)
  probit highwage collgrad $xx i.countyid if married==0 //单身组进行Probit回归
  margins,dydx(collgrad) post
  matrix `b0'=e(b)                 //获取单身组的平均边际效应b0 (即单身组的教育回报率)
  matrix `b'=`b1'-`b0'             //获取两组样本的平均边际效应之差 (b1-b0)
  ereturn post `b'
end

. *使用Bootstrap方法运行上述命令"AME_Probit"，重复进行500次
. bootstrap _b, nowarn reps(500) seed(10101): AME_Probit

Bootstrap results                                        Number of obs = 2,246
                                                         Replications  =   500
------------------------------------------------------------------------------
             |   Observed   Bootstrap                         Normal-based
             | coefficient  std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
    collgrad |    -0.0338     0.0343    -0.98   0.325      -0.1010      0.0334
------------------------------------------------------------------------------

. est store AME_diff_Probit //将组间差异检验的结果存储到"AME_diff_Probit"

从 bootstrap 的运行结果看，检验统计量 Z 值为-0.985，P 值大于 0.1，可认为 AME 系数在 10% 水平上并不存在显著差异。类似地，为了更直观地展示 collgrad 在已婚样本、单身样本中的平均边际效应及其组间差异，我们可以使用 est store 和 esttab 命令。

. probit highwage collgrad $xx i.countyid if married==1
. margins,dydx(collgrad) post
. est store AME_Group1_Probit     //保存已婚组的AME系数
. probit highwage collgrad $xx i.countyid if married==0
. margins,dydx(collgrad) post
. est store AME_Group0_Probit     //保存单身组的AME系数
. esttab AME_Group1_Probit AME_Group0_Probit AME_diff_Probit, ///
>        b(%9.4f) se mtitle star(* 0.1 ** 0.05 *** 0.01) modelwidth(16)

------------------------------------------------------------------------
                          (1)                 (2)                 (3)   
             AME_Group1_Pro~t    AME_Group0_Pro~t     AME_diff_Probit   
------------------------------------------------------------------------
collgrad               0.2568***           0.2906***          -0.0338   
                     (0.0194)            (0.0268)            (0.0343)   
------------------------------------------------------------------------
N                        1325                 726                2246   
------------------------------------------------------------------------
Standard errors in parentheses
* p<0.1, ** p<0.05, *** p<0.01

基于 Logit 模型的操作步骤是高度相似的，我们只需要将自己编写的命令名字从 AME_Probit 改为 AME_Logit，同时将所有涉及回归的命令从 probit 改成 logit 即可。在此不再赘述。

5. AME 系数的组间差异检验：引入工具变量

在上面的系数比较中，我们都假定了原因变量 collgrad 是外生的 (总体模型的误差项 $u$ 均值独立与原因变量collgrad)。如果原因变量 collgrad 具有内生性 (例如存在未观测到的遗漏变量，导致 collgrad 与误差项 $u$ 具有相关性) ，我们可以使用工具变量法，获得对 collgrad 的平均边际效应的一致估计。简便起见，我们暂不讨论因果效应具有个体异质性时的局部平均因果效应 (LATE)。

本节的模型设定是，我们以工资对数 lnwage 作为结果变量 (连续型变量)，是否接受大学教育 collgrad 作为原因变量 (核心解释变量)，以工时 hours，工龄 ttl_exp 作为一组控制变量，并设置区县固定效应 i.countyid。

与上文不同的是，我们首先为 collgrad 生成一个有效的工具变量 $Z$ ,使得工具变量 $Z$ 满足外生性假定和相关性假定，也即在给定控制变量的情形下，工具变量 $Z$ 与误差项 $u$ 不相关，而且工具变量 $Z$ 与原因变量 $u$ 具有偏相关关系。工具变量 $Z$ 的生成步骤如下：

. gen e=rchi2(3)-3
. gen Z_instrument=0.3+3.5*collg+e  //生成有效的工具变量，与原因变量相关，与误差项无关

接下来，通过分组回归的方式，分别在已婚组和单身组当中使用两阶段最小二乘法 (Two Stage Least Squares, TSLS) 估计 collgrad 的平均边际效应 AME 系数，进而检验 AME 系数是否存在显著的组间差异。

*手动编写AME系数组间差异检验的命令
capture program drop AME_TSLS     //清空历史命令
program define AME_TSLS, eclass   //命令开始，将这段命令的名字定义为"AME_TSLS"
  preserve
  tempname b b1 b0                 //定义三个用于临时存储的名字
  ivreg2 lnwage (collgrad=Z) $xx i.countyid if married==1 //已婚组进行TSLS回归
  margins,dydx(collgrad) post
  matrix `b1'=e(b)                 //获取已婚组的平均边际效应b1 (即已婚组的教育回报率)
  ivreg2 lnwage (collgrad=Z) $xx i.countyid if married==0 //单身组进行TSLS回归
  margins,dydx(collgrad) post
  matrix `b0'=e(b)                 //获取单身组的平均边际效应b0 (即单身组的教育回报率)
  matrix `b'=`b1'-`b0'             //获取两组样本的平均边际效应之差 (b1-b0)
  ereturn post `b'
end

. *使用Bootstrap方法运行上述命令"AME_TSLS"，重复进行500次
. bootstrap _b, nowarn reps(500) seed(10101): AME_TSLS

Bootstrap results                                        Number of obs = 2,246
                                                         Replications  =   500

------------------------------------------------------------------------------
             |   Observed   Bootstrap                         Normal-based
             | coefficient  std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
    collgrad |    -0.1623     0.1082    -1.50   0.134      -0.3744      0.0498
------------------------------------------------------------------------------

. est store AME_diff_TSLS //将组间差异检验的结果存储到"AME_diff_TSLS"

类似地，为直观地展示 collgrad 在已婚样本、单身样本中的平均边际效应及其组间差异，我们可以使用 est store 和 esttab 命令。

. ivreg2 lnwage (collgrad=Z) $xx i.countyid if married==1,r
. margins,dydx(collgrad) post
. est store AME_Group1_TSLS     //保存已婚组的AME系数
. ivreg2 lnwage (collgrad=Z) $xx i.countyid if married==0,r
. margins,dydx(collgrad) post
. est store AME_Group0_TSLS     //保存单身组的AME系数
. esttab AME_Group1_TSLS AME_Group0_TSLS AME_diff_TSLS,b(%9.4f) ///
>        se mtitle star(* 0.1 ** 0.05 *** 0.01) modelwidth(16)

------------------------------------------------------------------------
                          (1)                 (2)                 (3)   
              AME_Group1_TSLS     AME_Group0_TSLS       AME_diff_TSLS   
------------------------------------------------------------------------
collgrad               0.3708***           0.5331***          -0.1623   
                     (0.0615)            (0.0884)            (0.1082)   
------------------------------------------------------------------------
N                        1325                 726                2246   
------------------------------------------------------------------------
Standard errors in parentheses
* p<0.1, ** p<0.05, *** p<0.01

6. AME 系数的组间差异检验：面板数据的线性固定效应情形

上文介绍的都是截面数据中的系数比较问题，本节扩展至面板数据的固定效应模型 (Fixed Effect Model)。需要注意的是，面板数据中存在误差项自相关的问题。换言之，同一个受访者在不同时期的误差项是相关的，例如：

C o v (u_{i t,} u_{i t - 1}) \neq 0

因此，我们往往需要使用 “聚类稳健标准误” (以个体或者更大范围的单位为聚类)。如果要使用 Bootstrap 方法获取标准误，应该采取分块自助法 (Block Bootstrap)，而不是简单的自助法，否则将可能低估标准误。下面我们演示，如何在面板数据的线性固定效应模型设定下，检验某个原因变量的平均边际效应是否存在显著的组间差异。

我们使用 Wooldridge (2010) 提供的示例数据集 nls81_87，这份数据是 1981-1987 年美国女性劳动状况的面板数据，包括了工资 lnwage、工龄 exper、受教育程度 educ、种族 black、是否为制造业从业者 manuf、年份 year 以及个体 id 等信息。我们以工资对数 lnwage 作为结果变量，以工龄 exper 作为核心解释变量 (原因变量)，并设置个体固定效应、时间固定效应 i.year。

我们想要考察工龄 exper 对其工资对数 lnwage 的平均边际效应，是否在 “制造业劳动者” 和 “非制造业劳动者” 这两组群体间存在着显著差异。需要特别注意的是，我们首先要生成一个新的 id 变量 gen newid=id，并且在设置面板数据格式时，要以 newid (而不是原有的 id) 作为 panel name。

这是因为，之后我们将进行分块自助抽样 (Block Bootstrap)，每一次抽样都是以个体 id 为单位 (而不是以人-年为单位)，一旦抽中个体 “张三”，那么意味着 “张三” 在每一年的观测值将全部进入到样本中。如果有些个体在同一次自助抽样中被抽中了两次或多次，Stata 将会给这些幸运的个体指派新的标识码，这些标识码就存放在 newid 当中。简单地说，新的标识变量 newid 让分块自助抽样 (Block Bootstrap) 能够顺利地重复运行下去。

. cnssc install bcuse, replace
. bcuse nls81_87.dta, clear
. gen newid=id //生成新的个体标识变量
. xtset newid year // 设置面板数据格式，注意panel name应该用newid而不是原有的id变量
. gen lnwage=ln(wage)  //生成工资对数

*手动编写AME系数组间差异检验的命令
capture program drop AME_Panel     //清空历史命令
program define AME_Panel, eclass   //命令开始，将这段命令的名字定义为"AME_Probit"
  preserve
  tempname b b1 b0                  //定义三个用于临时存储的名字
  xtreg lnwage exper c.exper#c.exper i.year if manuf==1,fe  //制造业组进行FE回归
  margins,dydx(exper) post
  matrix `b1'=e(b)                 //获取制造业组的平均边际效应b1
  xtreg lnwage exper  c.exper#c.exper i.year if manuf==0,fe  //非制造业组进行FE回归
  margins,dydx(exper) post
  matrix `b0'=e(b)                 //获取非制造业组的平均边际效应b0
  matrix `b'=`b1'-`b0'             //获取两组样本的平均边际效应之差 (b1-b0)
  ereturn post `b'
end                                //自己编写的命令到此结束

. *使用Bootstrap方法运行上述命令"AME_Panel"，重复进行500次。
. bootstrap _b, cluster(id) idcluster(newid) nowarn reps(500) seed(10101): AME_Panel

Bootstrap results                                        Number of obs = 3,710
                                                         Replications  =   500

                                    (Replications based on 530 clusters in id)
------------------------------------------------------------------------------
             |   Observed   Bootstrap                         Normal-based
             | coefficient  std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
       exper |     0.0805     0.0716     1.12   0.261      -0.0599      0.2208
------------------------------------------------------------------------------

. est store AME_diff_FE //将组间差异检验的结果存储到"AME_diff_FE"

上述操作与截面数据的情形大同小异，需要注意最后一步，使用 bootstrap 命令时，需要在选项中写明 cluster (原有的个体标识变量) 以及 idcluster (新生成的个体标识变量)，这样才能让 “以受访者个体为单位的分块自助抽样” 顺利地进行下去。另外一个小细节是，我们在回归中通过因子变量的形式加入了工龄的二次项 c.exper#c.exper，以捕捉工龄对工资的非线性影响，高次项的设定并不会对平均边际效应的计算带来困扰，反而可能改善 AME 系数的估计。

最后，为了直观展示工龄 exper 在制造业组、非制造业组样本中的平均边际效应及其组间差异，我们可以用 est store 和 esttab 命令。从下表的三列系数可以看出，虽然 AME 系数在两组间存在较大的差异 (非制造业组的 AME 系数值仅为制造业组 AME 系数值的 50%) ，但是其组间差异并不统计显著。

. xtreg lnwage exper c.exper#c.exper i.year if manuf==1,fe vce(cl id)  //制造业组进行FE回归
. margins,dydx(exper) post
. est store AME_Group1_FE     //保存制造业组的AME系数
. xtreg lnwage exper c.exper#c.exper i.year if manuf==0,fe vce(cl id)  //非制造业组进行FE回归
. margins,dydx(exper) post
. est store AME_Group0_FE     //保存非制造业组的AME系数
. esttab AME_Group1_FE AME_Group0_FE AME_diff_FE,b(%9.4f) ///
>        se mtitle star(* 0.1 ** 0.05 *** 0.01) modelwidth(16)

------------------------------------------------------------------------
                          (1)                 (2)                 (3)   
                AME_Group1_FE       AME_Group0_FE         AME_diff_FE   
------------------------------------------------------------------------
exper                  0.1624***           0.0820*             0.0805   
                     (0.0493)            (0.0440)            (0.0716)   
------------------------------------------------------------------------
N                         650                1606                3710   
------------------------------------------------------------------------
Standard errors in parentheses
* p<0.1, ** p<0.05, *** p<0.01

7. 总结

本文聚焦于比较分组回归系数的组间差异，针对常用的 bdiff 命令可能遭遇的两处困惑，提出了可行的化解思路。本文认为，分析者可以聚焦于某一个 “核心解释变量” 对应系数 (平均边际效应) 的组间差异。这样在很大程度上可以克服 bdiff 要求两组变量 “对应相同” 的约束，同时也将关注重点从回归系数转向了更具可比性的平均边际效应 (AME)。

接下来，本文介绍了多种模型设定情境下的 AME 系数比较问题，这些情境包括截面 OLS，截面 Probit 或 Logit，截面工具变量回归 (TSLS) 以及面板固定效应 (FE)。

需要指出的是，本文的化解思路也并非万金油，它需要研究者自行编写一小段命令，并不如常用的 bdiff 命令那样简洁。另外，本文介绍的方法主要适用于检验某一个核心解释变量对应系数的组间差异，如果研究者关心更多变量对应系数的组间差异，那么 bdiff 命令可能会更加适用。

8. 参考资料

Cameron, A. Colin, and Douglas L. Miller. "A practitioner’s guide to cluster-robust inference." Journal of Human Resources 50.2 (2015): 317-372. -PDF
Mize, Trenton D., Long Doan, and J. Scott Long. "A general framework for comparing predictions and marginal effects across models." Sociological Methodology 49.1 (2019): 152-189. -PDF
Wooldridge, Jeffrey M. Econometric analysis of cross section and panel data. MIT press, 2010.
陈强, 2014. 高级计量经济学及 Stata 应用(第二版)[M]. 高等教育出版社.
洪岩璧. 2015. Logistic 模型的系数比较问题及解决策略：一个综述[J]. 社会, 35(4): 220-241.
连玉君, 廖俊平. 如何检验分组回归后的组间系数差异？[J]. 郑州航空工业管理学院学报, 2017, 35(6): 97-109.

9. 相关推文

Note：产生如下推文列表的 Stata 命令为：
lianxh 分组回归边际效应 cluster, m
安装最新版 lianxh 命令：
ssc install lianxh, replace

专题	嘉宾	直播/回看视频
⭐ 最新专题		文本分析、机器学习、效率专题、生存分析等
研究设计	连玉君	我的特斯拉-实证研究设计，-幻灯片-
面板模型	连玉君	动态面板模型，-幻灯片-
面板模型	连玉君	直击面板数据模型 [免费公开课，2小时]