Stata因子变量：虚拟变量-交乘项批量处理

发布时间：2020-06-16 阅读 13904

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

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

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

✨ 课程主页： https://gitee.com/arlionn/Course

作者：连玉君 | 杨柳

Email： arlionn@163.com | philoyl@163.com

注：该文已发表: 连玉君, 杨柳.《郑州航空工业管理学院学报》, 2018, 36(2): 90-103. 【 -点击下载-】 || 本文PDF版本2

事倍功半 vs. 事半功倍

当需要控制公司个体效应、面对成千上万家公司的数据资料时，你会如何处理？或许你会首先想到对每家公司生成虚拟变量，但是这种做法的 工作量实在是太大了！ 那么，有没有 事半功倍 的方法呢？下面，小编就带你学习 Stata 软件中 因子变量 的使用方法。

1. 问题背景

实证分析中，我们经常需要在模型中加入反映类别的 虚拟变量，以便控制不可观测的 组间差异。而在另一些分析中，为了刻画 调节效应，尚需在模型中加入变量的 交乘项 或 平方项。传统的做法是，预先生成虚拟变量或交乘项，进而将它们加入模型。然而，当虚拟变量或交乘项的数目较多时，上述方法就显得尤为不便，不但浪费计算机内存，也会严重降低我们的工作效率。

在 Stata 中，我们可以使用 因子变量 (Factor Variable) 简化操作步骤、快捷地在回归模型中加入虚拟变量、交乘项、平方项或高次项。值得注意的是，由于引入交乘项或平方项后，解释变量对被解释变量的边际影响不再是常数，而是某个变量（调节变量）的函数，在有些模型设定下，这种关系可能是非线性的。此时，若使用 因子变量 ，并配合 Stata 中的 margins 和 marginsplot 命令，可以非常便捷、直观地分析关键变量的 边际效应 并进行图形化展示。

2. 什么是因子变量

因子变量 (Factor Variable) 是对现有变量的延伸，是从类别变量中生成虚拟变量、设定类别变量之间的交乘项、类别变量与连续型变量之间的交乘项或连续变量之间的交乘项（或多项式）。在 Stata 中的大多数回归命令和回归后的估计命令中都可以使用这些因子变量 (详情参阅 Stata 帮助文件 help fvvarlist)。

因子变量的五种运算符及其含义如下表 1 所示：

表 1. 因子变量的运算符及含义

运算符	含义
`i.`	标示为类别变量的一元运算符
`c.`	标示为连续变量的一元运算符
`o.`	省略一个变量或指示变量的一元运算符
`#`	交乘项的二元运算符（在模型中仅包含交乘项）
`##`	两个变量及其交乘项的二元运算符（在模型中既包含交乘项，也包含原始变量）

注：
(a). 由因子变量运算符生成的指示变量和交乘项是实际存在的变量，它们与数据表中的变量是一样的，但是在数据表中并不显示出来；
(b). 类别变量的值必须是非负整数，范围介于 [0, 32740] 之间；
(c). 因子变量的运算符有时可与时间序列的运算符 L. 和 F. 组合在一起使用。

接下来，我们以研究 妇女工资的决定因素 为例，使用 Stata 软件自带的数据文件 nlsw88.dta 介绍因子变量的使用。该数据包含了 1988 年采集的 2246 个美国妇女的资料，包括：小时工资 wage、每周工作时数 hours、种族 race、职业 occupation、年龄 age、是否大学毕业 collgrad、当前职业的工作年数 tenure、是否结婚 married、是否居住在南部地区 south、合计工作年数 ttl_exp 等变量。其中，小时工资 wage、每周工作时数 hours、年龄 age、当前职业的工作年数 tenure、合计工作年数 ttl_exp 为连续型变量；种族 race 为类别变量（ 1 代表白种人 white，2 代表黑种人black，3 代表其他人种 other）、职业 occupation 为类别变量（13 个职业类别）；是否大学毕业 collgrad、是否结婚 married、是否居住在南部地区 south 为虚拟变量。

在这份数据中有一个表示种族的类别变量 race，取值为 1、2、3，分别对应白人、黑人和 其他人种。假设我们想在模型中加入一个反映种族的虚拟变量 black，当某个妇女是黑人时，black 取值为 1，否则为 0。则传统的做法如下：
(为了便于说明，后续多数回归命令中都省略了控制变量)

. sysuse "nlsw88", clear
. gen black=1
. replace black=0 if race!=2
. reg wage black

若延续这一思路，但使用因子变量来生成 black 变量，则命令为：

. gen black = 2.race

只需要一条命令，而且含义非常明确。这里 2.race 本质上是一个条件判断语句：判断某一行观察值中的 race变量取值是否为 2，若是，则返回 1 到变量 black 中，否则返回 0。

然而，在多数情况下，我们的目的只是希望得到虚拟变量 black 的估计系数，不希望生成或存储这个变量 (事实上，只要你的数据中存储了 race 变量，我们只需要保存好 dofile 文件，就无需生成 black 这个中间变量)。Stata 中的因子变量语法完全注意到了这个问题，使用因子变量的标准做法如下：

. sysuse "nlsw88.dta", clear
. reg wage 2.race

注意，我们无需预先生成 black 变量，而是直接在回归模型中加入了 2.race 因子变量。有些读者注意到，race 变量有三个取值，因此，我们可以在模型中放入两个虚拟变量，此时可以书写如下命令：

. reg wage i.race

回归结果如下所示：

---------------------------------------------------
   wage |      Coef.   Std. Err.      t    P>|t|   
--------+------------------------------------------
   race |
 black  |  -1.238442   .2764488    -4.48   0.000   
 other  |   .4677818   1.133005     0.41   0.680   
        |
  _cons |   8.082999   .1416683    57.06   0.000   
---------------------------------------------------

从上述结果中可以看到，回归模型中加入了 black 虚拟变量和 other 虚拟变量，分别对应 race 变量的第二个和第三个类别，而第一个类别 white 被 Stata 默认作为基准组，目的在于防止完全共线性。black 变量的系数值为 -1.238，表示黑人的平均工资比白人低 1.238 个单位，并在统计上显著；other 变量的系数值为 0.468，表示其他人种的平均工资比白人高 0.468 个单位，但是在统计上不显著。

在实证分析中，有时会根据研究内容的需要改变基准组的设定，此时可使用 ib. 或 b. 的前缀，具体写法如下表 2 中所描述：

表 2. 设定基准组的运算符及含义

基准组运算符	含义
`ib*.`	使用作为基准组，为变量中其中一类的值
`ib(#*).`	使用变量值中的第*位排序的值所对应的类别作为基准组
`ib(first).`	使用变量的最小值所对应的类别作为基准组（该项为 Stata 默认选项）
`ib(last).`	使用变量的最大值所对应的类别作为基准组
`ib(freq).`	使用变量值的频数最大的类别作为基准组
`ibn.`	不设基准组

注：(a). i 可以省略不写。例如：ib2.group 与 b2.group 的写法等价；(b). ib(#2). 指使用变量值中的第二位排序的值所对应的类别作为基准组；(c). 如果想在线性回归模型中将 group 变量中值为 3 的类别设置为基准组，则命令可写为：reg y i.sex ib3.group 。

在实证分析中，变量的 交乘项 或 高次项 往往是重要的解释变量。以研究妇女工资的决定因素为例，若我们想在模型中加入 黑人每周工作时数 的变量，则传统做法是预先生成一个虚拟变量 black 代表是否黑人，再生成一个新变量 black_x_hours 表示黑人与每周工作时数的交乘项，然后再将这个新变量放入回归模型中。Stata 命令为：

. gen black=1
. replace black=0 if race!=2
. gen black_x_hours = black * hours
. reg wage black_x_hours

若在 Stata 中使用因子变量实现上述过程，则命令十分简洁：

. reg wage 2.race#c.hours

下面，我们介绍如何在 Stata 中使用 因子变量 表示变量的 交乘项 或 高次项。我们以研究妇女工资的决定因素为例进行说明。

（1）两个类别变量的交乘项： 在回归模型中加入种族 race 和职业类别 occupation 的 交乘项，Stata 命令为：

. reg wage i.race#i.occupation

若在回归模型中既要放入种族 race 和职业类别 occupation 的虚拟变量，又需要同时放入这两个变量的 交乘项，则在回归命令中使用 i.race##i.occupation，相应的 Stata 命令为：

. reg wage i.race##i.occupation

（2）类别变量与连续变量的交乘项： 在回归模型中加入种族 race 和每周工作时数 hours 的交乘项，Stata 命令为：

. reg wage i.race#c.hours

需要注意的是，在上例中，由于我们把 hours 变量视为 连续变量，因此，需要在其前面加上 c. 符号以便告知 Stata 该变量是 连续变量。

（3）连续变量与连续变量的交乘项（高次项）： 在回归模型中加入年龄 age 变量，及其 平方项，Stata 命令如下：

. reg wage c.age##c.age

上述命令中 c.age 表示年龄 age 变量被当成 连续型变量。如果我们在 Stata 命令中使用 i.age，则年龄 age 变量被当成 类别变量 处理，此时，类别的个数为年龄 age 变量中不同取值的个数。

3. 常用回归模型的因子变量表述

3.1 范例 1：邹氏检验

由于不同组别之间可能会存在差异 (截距项或斜率项存在差异)，因此，我们需要检验这些差异在统计上是否显著。这时，我们可以使用 邹氏检验。以研究妇女工资的影响因素为例，我们可以使用 chowtest 命令检验工会成员与非工会成员两个样本组中工资影响因素是否存在差异（或称之为存在结构变化），Stata 命令如下：

. global xx "hours age tenure ttl_exp married"
. chowtest wage $xx, group(union) detail

事实上，我们也可以使用 因子变量 的语法，在回归模型中加入 分组变量 与 其他控制变量 的 交乘项，然后再 联合检验 分组变量的系数以及所有交乘项的系数是否都等于 0。此时，即使我们不使用 chowtest 命令，也可以轻松实现 邹氏检验，Stata命令如下：

. global xx "hours age tenure ttl_exp married"
. reg wage $xx i.union i.union#c.($xx)
. testparm i.union i.union#c.($xx)

3.2 范例 2：双向固定效应模型

在实证分析过程中，经常需要在模型中加入反映 年度、公司或行业特征 的 虚拟变量。当虚拟变量的数目众多时，采用手动输入变量的方式会非常耗时。例如，下述模型 (1) 是文献中广泛应用的 双向固定效应模型：

y_{i t} = α_{i} + λ_{t} + {x_{i t}}^{^{'}} β + ε_{i t} (1)

上式等价于

y_{i t} = \sum_{i = 1}^{N} D_{i} \cdot α_{i} + \sum_{t = 2}^{T} W_{t} \cdot λ_{t} + {x_{i t}}^{^{'}} β + ε_{i t} (2)

其中， $D_{i}$ 是对应于每家公司的虚拟变量，对于公司 $i$ ， $D_{i} = 1$ ，否则 $D_{i} = 0$ 。对于一个有 $N = 1000$ 家公司的面板数据而言，共有 1000 个虚拟变量。 $W_{t}$ 是 $T - 1$ 个年度虚拟变量 (之所以加入 $T - 1$ 个年度虚拟变量，是为了防止完全共线性)，定义方式与 $D_{i}$ 相似。

若使用因子变量，则语法很简单：

. reg y x1 x2 x3 i.id i.year

其中，i.id 表示 $N - 1$ 个公司的虚拟变量（Stata 会默认将一家公司作为基准组），即 (1) 式中的 $α_{i}$ 。此处设定 i.year 会自动加入 $T$ 个年度虚拟变量，但 Stata 会自动删除一个，以防完全共线性。

当然，我们也可以用 xtreg 命令自带的 fe 选项来控制个体效应，同时使用因子变量来加入年度虚拟变量，命令如下：

. xtreg y x1 x2 x3 i.year, fe

3.3 范例 3：DID 模型

在实证分析中，常常需要分析 政策实施后带来的效果。这时，我们就需要采集 实验组 与 控制组 两组样本（ 实验组 的样本代表政策实施前后的情况，控制组 的样本代表不实施政策的情况），再把这两组样本合并为一份数据后进行回归分析。在建模时，我们需要设定一个是否为实验组或控制组的处理虚拟变量 $T r e a t$ 放入模型，还需要设定一个反映政策实施时间前与实施后的时间虚拟变量 $T i m e$ 放入模型，此外，还需要在模型中加入 $T r e a t$ 变量与 $T i m e$ 变量的 交乘项，该交乘项的系数估计值 $a_{3}$ 就是政策实施后带来的效果，即实验组在政策实施后与假设未实施政策的差异。回归模型如下式 (3) 所示：

y = a_{0} + a_{1} T r e a t + a_{2} T i m e + a_{3} T r e a t \times T i m e + x^{^{'}} β + ε (3)

在 Stata 中的命令写法如下：

. reg y Treat Time Treat#Time x1 x2 x3
*-或者写为
. reg y Treat##Time x1 x2 x3

在多期 DID 分析中，我们常常需要加入 年度虚拟变量 与 处理虚拟变量 的 交乘项 来检验 共同趋势假设 (common trend) 以及 政策效果。例如，在 Acemoglu and Angrist(2001) 文中，作者采集了 1988 年至 1997 年的人口调查数据，将样本分为残疾人与非残疾人，使用多期 DID 模型研究了 1992 年美国对残疾人工作保护法案（ADA）的实施效果。文中模型设定如下 (pp. 925)：

y_{i t} = a_{0} + {x_{i t}}^{^{'}} \cdot γ + {p_{i t}}^{^{'}} \cdot η + x_{i t}^{^{'}} \cdot p_{i t} \cdot β + d_{i t}^{^{'}} \cdot δ + d_{i t}^{^{'}} \cdot p_{i t} \cdot α_{t} + ε_{i t} (4)

其中， $i$ 为残疾人个体； $t$ 为年份； $y_{i t}$ 为工作周数或周平均工资； $x_{i t}^{^{'}}$ 为一系列控制变量， $p_{i t}$ 为时间趋势项； $d_{i t}$ 为是否残疾，其系数为 $δ$ ； $α_{t}$ 为随年份变化的实施 ADA 与未实施 ADA 对残疾人工作保护效应的系数。当 $t \geq 1992$ 年时， $α_{t}$ 量化了实施 ADA 之后对残疾人工作保护的效果（使用非残疾人作为控制组）；当 $t < 1992$ 年时， $α_{t}$ 检验了实施 ADA 之前的各年份中残疾人和非残疾人的工作周数或周平均工资是否在统计上具有显著差异，相当于同时进行了DID模型的 共同趋势假设 的检验，若 $α_{t}$ 系数在 1992 年之前不显著，表明存在 共同趋势。Stata 命令与结果如下：

(文中实证分析所用的原始数据和相关程序可以从作者主页上下载：http://economics.mit.edu/faculty/acemoglu/data/aa2001)

. use "ABA_JPE2001.dta", clear
. global controls "i.age_G i.edu_G i.race_G i.region"
. reg wkswork1 i.year##($controls) i.disabled##i.year

Stata 回归结果如下：

----------------------------------------------------------------
      wkswork1 |  Coef.   S.E.     t     P>|t|      [95% CI]
---------------+------------------------------------------------
(部分回归结果省略
               
 disabled#year 
  Disabl#1988  | -0.687   0.720  -0.95   0.340    -2.098   0.723
  Disabl#1989  | -0.539   0.706  -0.76   0.445    -1.922   0.845
  Disabl#1990  | -2.290   0.695  -3.30   0.001    -3.651  -0.929
  Disabl#1991  | -2.159   0.692  -3.12   0.002    -3.514  -0.803
  Disabl#1992  | -1.387   0.687  -2.02   0.044    -2.734  -0.039
  Disabl#1993  | -2.743   0.695  -3.95   0.000    -4.106  -1.381
  Disabl#1994  | -3.918   0.720  -5.44   0.000    -5.329  -2.507
  Disabl#1995  | -3.703   0.748  -4.95   0.000    -5.169  -2.236
  Disabl#1996  | -4.472   0.737  -6.07   0.000    -5.918  -3.027
               
         _cons | 41.052   0.455  90.30   0.000    40.161  41.943
----------------------------------------------------------------

3.4 范例 4：超越对数生产函数

很多文献使用超越对数生产函数估计生产效率。例如，Altunbas, Liu, Molyneux and Seth (2000) 使用超越对数成本函数估算了日本银行的效率和风险；Wang (2007) 则使用超越对数生产函数研究30个国家R&D效率；王德祥和李建军(2009) 基于超越对数生产函数估算了我国的税收流失率。

超越对数生产函数模型中会包含投入要素 $K$ 和 $L$ 的高阶项，包括平方项和二者的交乘项。例如，在 Kumbhakar (1989) 中，作者通过在模型中加入二次项来捕捉要素的交互影响和潜在的非线性关系，模型如下式 (5) 所示。

l n y = α + \sum_{k = 1}^{k} β_{k} \cdot l n x_{k} + \frac{1}{2} \sum_{k = 1}^{k} \sum_{m = 1}^{m} γ_{k m} \cdot l n x_{k} \cdot l n x_{m} (5)

此时，使用因子变量会让 Stata 中的命令变得异常简洁，如下所示：

. webuse frontier1.dta, clear  
. global y "lnoutput"            // y
. global x "lnlabor lncapital"   // x1 x2
*-超越对数生产函数
. sfcross $y c.($x)##c.($x)      // Eq.(5)

4. 边际效应分析和图形化呈现

一些实证研究的文献中常常加入 变量的交乘项 以反映可能存在的 调节效应。例如， Faulkender and Wang (2006) 检验了由公司融资约束对现金持有边际市场价值的影响。戴魁早和刘友金 (2016) 研究发现了要素市场扭曲对创新效率的影响存在着企业差异，企业规模在规避要素市场扭曲对创新效率的抑制效应中具有积极作用。张苏和高扬 (2012) 的实证研究发现学生来源于城市或农村地区对国家竞争力的影响作用受到每周上网时间的影响，如果上网时间在每周 8 小时以下，城市大学生的学习行为落入 “增进国家竞争力导向” 上高效率区域而不是低效率区域的概率比农村大学生要高，否则要低。

在类似上述的研究中，最重要的是分析 $x$ 对 $y$ 的 边际效应。

y = α_{1} + x β_{1} + z β_{2} + (x \cdot z) β_{3} + ε (6)

$x$ 对 $y$ 的边际影响为：

\frac{\partial y}{\partial x} = β_{1} + z β_{3}

显然， $x$ 对 $y$ 的边际影响取决于 $z$ 的取值。这种边际效应分析若是采用手动计算会非常繁琐，而后续的图形呈现则更为复杂。然而，若是借助 因子变量，并配合 Stata 中的 margins 和 marginsplot 命令，对于包含 交乘项 的模型中的 边际效应 的分析和 图形化 展示都变得异常轻松。

下面，我们以研究妇女工资的决定因素为例进行说明。使用 Stata 软件的自带数据 nlsw88.dta(1988年美国妇女小时工资)，以 wage (妇女的小时工资) 作为被解释变量，以 race (种族类别)，collgrad (是否大学毕业)，race 与 collgrad 的交乘项 作为解释变量建立线性回归模型，Stata 中的命令如下：

. sysuse "nlsw88.dta", clear
. reg wage i.race##collgrad

回归结果如下所示：

(部分回归结果省略)
----------------------------------------------------------
               wage |  Coef.   Std. Err.      t    P>|t|  
--------------------+-------------------------------------
               race |
             black  | -1.442      0.297    -4.85   0.000  
             other  | -0.380      1.348    -0.28   0.778  
                    |
           collgrad |
      college grad  |  2.981      0.312     9.54   0.000  
                    |
      race#collgrad |
black#college grad  |  2.502      0.676     3.70   0.000  
other#college grad  |  1.678      2.297     0.73   0.465  
                    |
              _cons |  7.318      0.158    46.25   0.000  
----------------------------------------------------------

我们使用 margins 命令来计算 race 与 collgrad 的交乘项各类别的妇女工资的 预测边际值，Stata 命令和结果如下所示：

. margins i.race#collgrad

计算结果如下所示：

Adjusted predictions                            Number of obs     =      2,246
Model VCE    : OLS

Expression   : Linear prediction, predict()

--------------------------------------------------------------------------
                    |         Delta-method
                    |  Margin   Std. Err.      t    P>|t|      [95% CI]
--------------------+-----------------------------------------------------
      race#collgrad |
             white #|
  not college grad  |   7.318      0.158    46.25   0.000    7.008   7.629
white#college grad  |  10.299      0.269    38.24   0.000    9.771  10.827
             black #|
  not college grad  |   5.876      0.252    23.32   0.000    5.382   6.370
black#college grad  |  11.359      0.544    20.89   0.000   10.292  12.425
             other #|
  not college grad  |   6.938      1.339     5.18   0.000    4.313   9.563
other#college grad  |  11.597      1.840     6.30   0.000    7.989  15.205
--------------------------------------------------------------------------

使用 `marginsplot` 命令将计算结果用图的形式表示，Stata 命令和结果如下所示：

```stata
. marginsplot

图 1. 种族与是否大学毕业交乘项各类别的妇女小时工资的预测边际值

从图 1 中可以直观的看到，不论种族类别，大学毕业的妇女的平均工资高于非大学毕业的妇女，该结果符合我们的一般认知。我们从图中还发现一个有趣的结果，非大学毕业的白人妇女的平均工资高于非大学毕业的黑人妇女，而大学毕业的白人妇女的平均工资低于大学毕业的黑人妇女。

下面，我们使用 margins 命令附加 atmeans 选项来计算 当其他变量取均值 时 不同种族类别 的妇女小时工资的 预测边际值，Stata 命令和结果如下所示：

. margins i.race, atmeans
*-或者写为
. margins race, atmeans

计算结果如下所示：

Adjusted predictions                Number of obs     =      2,246
Model VCE    : OLS

Expression   : Linear prediction, predict()
at           : 1.race          =    .7288513 (mean)
               2.race          =    .2595726 (mean)
               3.race          =    .0115761 (mean)
               0.collgrad      =    .7631345 (mean)
               1.collgrad      =    .2368655 (mean)

-----------------------------------------------------------
        |        Delta-method
        | Margin   Std. Err.      t    P>|t|     [95% CI]
--------+--------------------------------------------------
   race |
 white  |  8.024      0.137    58.76   0.000   7.756  8.292
 black  |  7.175      0.231    31.00   0.000   6.721  7.628
 other  |  8.042      1.111     7.24   0.000   5.864 10.220
-----------------------------------------------------------

使用 marginsplot 命令将计算结果用图的形式表示，Stata 命令和结果如下所示：

. marginsplot

杨柳因子变量Fig02边际效应示例2 图 2. 不同种族类别的妇女小时工资的预测边际值

从图 2 中可以直观的看到，白人与其他种族的妇女的平均工资高于黑人妇女。

两个 连续变量 的 交乘项 对被解释变量的边际效应也可以使用 margins 命令来计算。我们仍以研究妇女工资的决定因素为例进行说明。在回归模型中加入 tenure (当前职业的工作年数) 及其平方项，并将 hours (每周工作时数)，age (妇女年龄)，married (是否结婚)，south (是否居住在南部地区)，race (种族类别) 作为控制变量，Stata 命令如下所示：

. global xx "hours age married south i.race"
. reg wage c.tenure##c.tenure $xx

回归结果如下：

(部分回归结果省略)
--------------------------------------------------------
             wage |  Coef.   Std. Err.      t    P>|t|  
------------------+-------------------------------------
           tenure |  0.302      0.070     4.28   0.000  
                  |
c.tenure#c.tenure | -0.008      0.004    -2.02   0.043  
                  |
            hours |  0.075      0.012     6.48   0.000  
              age | -0.085      0.039    -2.19   0.029  
          married | -0.462      0.254    -1.82   0.070  
            south | -1.284      0.249    -5.15   0.000  
                  |
             race |
           black  | -1.221      0.287    -4.25   0.000  
           other  |  0.313      1.098     0.28   0.776  
                  |
            _cons |  8.184      1.632     5.01   0.000  
--------------------------------------------------------

由于在模型中加入了 $t e n u r e$ 的平方项，因此， $t e n u r e$ 对 $w a g e$ 的边际效应将会受到 $t e n u r e$ 取值的影响。上述 Stata 命令对应的模型设定如下：

w a g e = α + t e n u r e β_{1} + (t e n u r e \cdot t e n u r e) β_{2} + x_{3} β_{3} + \dots + x_{8} β_{8} + ε (7)

$t e n u r e$ 对 $w a g e$ 的边际影响为：

\frac{\partial w a g e}{\partial t e n u r e} = β_{1} + 2 \cdot t e n u r e \cdot β_{2}

显然，当 $t e n u r e$ 取值不同时， $t e n u r e$ 对 $w a g e$ 的边际效应是不相同的。因此，需要先使用描述性统计分析的命令 (如 sum tenure) 查看 $t e n u r e$ 取值的范围，然后再计算当 $t e n u r e$ 取不同的值时

对应的边际效应，Stata 命令如下：

. preserve
. keep if e(sample)
. sum tenure
. restore
. margins, dydx(tenure) at(tenure=(0 1(3)25 25.9))

计算结果如下所示：

. sum tenure

    Variable |   Obs        Mean    Std. Dev.  Min        Max
-------------+-----------------------------------------------
      tenure | 2,227    5.981702    5.509708     0   25.91667


. margins, dydx(tenure) at(tenure=(0 1(3)25 25.9))

Average marginal effects                        Number of obs     =      2,227
Model VCE    : OLS

Expression   : Linear prediction, predict()
dy/dx w.r.t. : tenure

1._at        : tenure   =       0

2._at        : tenure   =       1

3._at        : tenure   =       4

4._at        : tenure   =       7

5._at        : tenure   =      10

6._at        : tenure   =      13

7._at        : tenure   =      16

8._at        : tenure   =      19

9._at        : tenure   =      22

10._at       : tenure   =      25

11._at       : tenure   =    25.9

-----------------------------------------------------------------
       |        Delta-method
       |  dy/dx   Std. Err.      t    P>|t|  [95% Conf. Interval]
-------+---------------------------------------------------------
tenure |
   _at |
    1  |  0.302      0.070     4.28   0.000     0.163       0.440
    2  |  0.286      0.063     4.52   0.000     0.162       0.410
    3  |  0.240      0.043     5.63   0.000     0.157       0.324
    4  |  0.195      0.026     7.50   0.000     0.144       0.246
    5  |  0.149      0.023     6.36   0.000     0.103       0.195
    6  |  0.103      0.038     2.71   0.007     0.028       0.178
    7  |  0.057      0.058     0.98   0.325    -0.057       0.171
    8  |  0.011      0.080     0.14   0.886    -0.145       0.167
    9  | -0.034      0.102    -0.34   0.735    -0.234       0.165
   10  | -0.080      0.124    -0.65   0.517    -0.323       0.162
   11  | -0.094      0.130    -0.72   0.471    -0.350       0.162
-----------------------------------------------------------------

使用 marginsplot 命令将计算结果用图的形式表示，Stata 命令和结果如下所示：

. marginsplot, xlabel(,format(%3.1f) angle(60))

杨柳因子变量Fig03边际效应示例3 图 3. 当前职业的工作年数对妇女工资的平均边际效应

Stata 中的多数命令都支持 margins 和 marginsplot 命令。因此，即使对于非线性模型，如 Logit, Tobit 等，我们仍然可以借助这两个命令很方便地分析边际效应。

5. 输出回归结果时的问题与解决办法

以研究妇女工资的决定因素为例。wage (妇女的小时工资) 作为回归模型的被解释变量，race (种族类别)、collgrad (是否大学毕业)、race 与 collgrad 的交乘项 作为解释变量，并将回归结果输出，Stata 命令和结果如下：

. sysuse "nlsw88.dta", clear
. reg wage i.race##collgrad
. est store R1
. esttab R1, nogap

----------------------------
                      (1)   
                     wage   
----------------------------
1.race                  0   
                      (.)   
2.race             -1.442***
                  (-4.85)   
3.race             -0.380   
                  (-0.28)   
0.collgrad              0   
                      (.)   
1.collgrad          2.981***
                   (9.54)   
1.race#0.c~d            0   
                      (.)   
1.race#1.c~d            0   
                      (.)   
2.race#0.c~d            0   
                      (.)   
2.race#1.c~d        2.502***
                   (3.70)   
3.race#0.c~d            0   
                      (.)   
3.race#1.c~d        1.678   
                   (0.73)   
_cons               7.318***
                  (46.25)   
----------------------------
N                    2246   
----------------------------
t statistics in parentheses
* p<0.05, ** p<0.01, *** p<0.001

我们发现在上述结果中有很多变量的系数值为 0，并缺失 $t$ 统计量。造成这种情况的原因有两个方面：一方面是由于 有些虚拟变量作为了基准组，例如： 1.race ， 0.collgrad 与 1.race#0.collgrad ，Stata 默认将它们作为基准组，所以就缺失这些基准组的估计系数值和 $t$ 统计量；另一方面是由于有些交乘项中的其中一个因子变量是基准组，其它变量与这个作为基准组的因子变量交乘后的交乘项就被忽略了，所以其估计系数值和 $t$ 统计量就会缺失，例如： 1.race#1.collgrad ， 2.race#0.collgrad 与 3.race#0.collgrad 。在这些情况下，可以使用 esttab 命令的 drop() 选项来屏蔽掉这些缺失系数的显示，还可以使用 nobase 和 noomit 的选项，Stata 命令和结果如下所示：

*-输出结果（不显示基准组和忽略组的系数，使用 drop 选项）

. esttab R1, nogap ///
   drop(1.race 0.collgrad 1.race#0.collgrad   ///
        1.race#1.collgrad 2.race#0.collgrad   ///
        3.race#0.collgrad)

*-输出结果（不显示基准组和忽略组的系数，使用nobase与noomit选项）
.    esttab R1, nogap nobase noomit

----------------------------
                      (1)   
                     wage   
----------------------------
2.race             -1.442***
                  (-4.85)   
3.race             -0.380   
                  (-0.28)   
1.collgrad          2.981***
                   (9.54)   
2.race#1.c~d        2.502***
                   (3.70)   
3.race#1.c~d        1.678   
                   (0.73)   
_cons               7.318***
                  (46.25)   
----------------------------
N                    2246   
----------------------------
t statistics in parentheses
* p<0.05, ** p<0.01, *** p<0.001

6. 小结

本文介绍了 Stata 中 因子变量产生虚拟变量与交乘项 的使用方法，以常用经典回归模型为例，提供了它们的 Stata 命令，并进一步提供了 因子变量 与 margins 与 marginsplot 命令相配合分析 边际效应 的示例。在 Stata 中的大多数命令中都可以使用 因子变量 的表述方法。该方法可以使 Stata 命令更加简洁，并能够大幅度提高实证分析的效率，但需要注意的是分析使用因子变量表述方法后得到的 模型设定。

参考文献：

Acemoglu, D., J. D. Angrist, 2001, Consequences of Employment Protection? The Case of the Americans with Disabilities Act, Journal of Political Economy, 109 (5): 915-957.
Altunbas, Y., M.-H. Liu, P. Molyneux, R. Seth, 2000, Efficiency and Risk in Japanese Banking, Journal of Banking & Finance, 24 (10): 1605-1628.
Faulkender, M., R. Wang, 2006, Corporate Financial Policy and the Value of Cash, Journal of Finance, 61 (4): 1957-1990.
Kumbhakar, S. C., 1989, Estimation of Technical Efficiency Using Flexible Functional Form and Panel Data, Journal of Business and Economic Statistics, 7 (2): 253-258.
Wang, E. C., 2007, R&D Efficiency and Economic Performance: A Cross-Country Analysis Using the Stochastic Frontier Approach, Journal of Policy Modeling, 29 (2): 345-360.
戴魁早, 刘友金, 2016, 要素市场扭曲与创新效率——对中国高技术产业发展的经验分析, 经济研究, (7): 72-86.
王德祥, 李建军, 2009, 我国税收征管效率及其影响因素——基于随机前沿分析(SFA)技术的实证研究, 数量经济技术经济研究, (4): 152-160.
张苏, 高扬, 2012, 大学生学习行为与国家竞争力关联关系的实证研究, 管理世界, (4): 175-176.

专题	嘉宾	直播/回看视频
⭐ 最新专题 ⭐		DSGE, 因果推断, 空间计量等
⭕ Stata数据清洗	游万海	直播, 2 小时，已上线
研究设计	连玉君	我的特斯拉-实证研究设计，-幻灯片-
面板模型	连玉君	动态面板模型，-幻灯片-
面板模型	连玉君	直击面板数据模型 [免费公开课，2小时]

Stata因子变量：虚拟变量-交乘项批量处理

1. 问题背景

2. 什么是因子变量

3. 常用回归模型的因子变量表述

3.1 范例 1：邹氏检验

3.2 范例 2：双向固定效应模型

3.3 范例 3：DID 模型

3.4 范例 4：超越对数生产函数

4. 边际效应分析和图形化呈现

5. 输出回归结果时的问题与解决办法

6. 小结

参考文献：

相关课程

课程一览

关于我们

热门资讯