常用命令

数据相关

生成新数据

利用旧数据定义新数据，注意，如果使用的旧数据为两个及以上，只要其中一个为空值，则新数据也为空值（即空值不会自动转化为0）：
gen new_valuename = f(value1，value2,……)
如果要生成时间序列，可以定义原有的时间序列数据：
tsset valueame
如果原本没有，可以定义一列新的时间序列：
gen t=_n
tsset t

删除和修改数据

删除数据：
drop value1 value2 value3 ……
保留数据（除保留的以外全部删除）：
keep value1 value2 value3 ……
修改数据：
replace valuename = f（value1 , value2 ,……）
修改数据时也可以将自己作为自变量，比如：
replace age = age -1
则age数据全部减少1
注意，请谨慎使用 drop 命令，该命令不可撤销，数据一旦丢失将无法追回，建议使用do文件减少风险（找不回，只是可以选择不保存，然后不用重输命令）

一键删除全部含有缺失值的对象：
egen mis = rowmiss(_all)
drop if mis
删除重复值（删除重复值后会保留第一个数据）：
duplicates drop valuename , force

改变数据类型

excel导入的时候，数值型的数据如果有缺失，可能会被导成数值型，可以利用以下操作强制改为数值型
destring value ,force

条件命令

条件命令为 if ，并且的符号为&，或者的符号为|，否认的符号为！，以下为示例：
drop value1 if value 1 >=1 & value 2>=3
drop value1 if value 1 >=1 | value 2>=3
drop value1 if value 4 ！= “error”

取对数命令

底数为e取对数：
log(value)
底数为10取对数：
log10(value)
取对数后新数值的系数表示百分比变化，如果回归结果不好，可以尝试对其中相对大的数据列取对数。

输出相关

将不同回归结果输出到同一张表中，其中的省略号为所作的操作，注意nest与逗号之间要有空格，不然会失败（第一条命令为安装命令，已经安装的可以跳过）：
ssc install asdoc, replace
asdoc …… , nest replace
asdoc …… , nest append
如果要输出描述统计或者其他结果，则可以用：
asdoc …… , save(docname.doc) title(titlename)

常用回归

非线性选择回归

logit回归

二元选择回归，因变量y为虚拟变量：
logistic y x1 x2 x3……
此外逻辑回归也可以做到分组回归，多元选择回归，在学习完成后在这个博客上将继续补充

probit回归

probit y x1 x2 x3……
与逻辑回归的用法相同，虽然思想不一样，但通常情况下，probit回归估计出的参数值乘以1.814，大致会等于logistic回归中的参数值。

线性回归

OLS

reg y x1 x2 x3……
如果有虚拟变量，可以使用xi命令快捷生成：
xi:reg y x1 x2 x3 i.x4……

Heckman 回归

选择样本模型的一种，这里给出heckman二步法的命令：
heckman y x1 x2 x3 …… , select(xi = z1 z2) twostep
后面的select变量中必须含有至少一种影响xi但对y没有影响的变量，xi为选择变量，这种方法会直接输出第二步回归时的R^2，但无第一步的信息，且默认第二步为OLS回归。可以使用下列方法分布计算：
dprobit xi x1 x2 z1 z2
predict zg if e(sample),xb
gen lambda = normalden(zg) / normal(zg)
reg y x1 x2 x3 …… lambda
符号含义同上，新增lambda变量意为“反Mills比例”，若其为负数则说明确实存在选择问题，使用该方法是合理的。另外，对最后一步可以替换为其他回归，例如2SLS等。

Tobit回归

选择样本模型的一种，只给出一般的tobit命令：
tobit y x1 x2 x3……

2SLS回归（工具变量法）

被解释变量：y

解释变量：x1

控制变量：x2 x3 x4

工具变量：z1 z2

reg x1 z1 x2 x3 x4, robust
ivregress 2sls y x1 x2 x3 x4 (x1= z1 z2), robust

常用检验

在做完回归后进行输入命令：

异方差检验

estat hettest
假如prob > chi2 的值小于0.05 则检验通过

多重共线性检验

estat vif
假如整体的vif小于10则通过

自相关检验

estat dwatson
一般DW值在2附近则通过（也可以用DW表对照自己判断）

格兰杰因果检验

存在滞后项时使用
先下载格兰杰因果检验程序gcause
ssc install gcause
检验x对y的因果关系，以下命令为滞后1、2期，显示其AIC、BIC的取值，一般的检验到滞后三期为止，之后进行y对x的检验，命令类似，不予贴出：
gcause y x,lags(1)
estat ic
gcause y x,lags(2)
estat ic

稳健性检验

通常是通过改变工具变量或是使用不同的计量方法，如果回归结果类似，则说明结果具有稳健性

内生性检验

通常意味着变量之间相互影响，可以考虑改变工具变量或是使用Durbin－Wu－Hausman 方法进行检验
代码为：
estat endogenous
输出结果示例：

原假设为“变量是外生的”，若p值较大则接受原假设。这里可采用较大的p值以求结果的可靠性。

建立新的变量——generate

generate [type] newva = exp [if] [in]

type是可选项 用于指定创建的变量的类型 类型包括byte、int、long、float、double，

Newva表明新创建的变量名称，exp是关于变量赋值的表达式。

例：

gen educ2 = educ^2

gen educ_exper = educ*exper if educ >= 9  //生成交叉项

gen logexper = log(exper)  //变量exper取对数

gen educ3 = round (sqrt(educ))  //round的含义是取整，sqrt()表示开根号，若（）内为负数该函数将返回缺失值

更改已有的变量——replace

replace oldvar = exp [if] [in] [,nopromate] 改变已经存在的变量的赋值

[,nopromate]：阻止Stata 改变变量的类型来适应变量新的赋值

gen eudcat = 0

replace educat = 1 if educ==4 | educ==5 | educ==6

replace educat = 2 if educ >= 7 & educat <= 9

创建变量——egen

egen [type] newvar = fcn(arguments) [if] [in] [,options]

其中，fcn(arguments)指特有函数

egen命令函数清单：

sd(exp) 创建标准差

skew(varname) 创建偏度

std(exp) 创建标准化的赋值，默认均值为0，标准差为1

total(exp) 数值为表达式exper的加总

iqr(exp) 创建含有四分位距的变量

kurt(varname) 创建峰度新变量

mean(exp) 创建平均数变量

median(exp) 创建中位数变量

egen educavg = mean(educ)
egen educed = educ - educavg
sum educ educavg educed  //总结这三个变量的关键统计指标，呈现基本情况
bysort female: egen educmed = median(educ)  //分性别计算变量，产生educ中位数
egen stdeduc = std(educ)
egen higheduc = anyvalue(educ), v(13/18)  //若educ的取值为13-18，则生成新变量higheduc赋值为educ；否则higheduc取值为缺失值“.”
list higheduc educ in 1/20  //列举前20个观测值，对比这两个变量。
egen sexmar = group(female married) //使用了group函数，说明用于分组的标准为female 和married

数值和字符串的转换——encode、decode

字符型到数值型:为已经存在的字符串变量添加一个去了标签的数值型变量

encode varname [if] [in], generate (newvar) [label(name)]

数值型到字符型：为已经存在的数值型变量和它的标签生成一个字符串变量

decode varname [if] [in], generate (newvar) [maxlength(#)]

real() 用于从合适的字符串表达式中得到相应的数值

encode sex, gen(gender)   //生成数值型变量gender

list sex gender in 1/4, nolabel   //不加标签的情况下罗列前四个观测值

label list gender  //查看gender标签的具体情况

decode gender, gen(gender2)  //根据gender变量重新生成一个字符型变量gender2

real("5.2")+1=6.2 //real()得到字符串的真实数值，如果不包含数字则返回缺失值

生成分类变量和虚拟变量

1.虚拟变量

gen college = 0

replace college=1 if educ>=12

list educ college in 1/10 //列表呈现前十个观测值


gen college = (educ>=12)

gen collede2 = (educ>12 & educ<17)

gen master = (educ >16 & educ<.)

tab1 edu0-master //查看新生成变量的频数表

2.分类变量

recode命令

recode oldvar (rule) [,gen(newvar)]

rule: #=#  3=1  3转换为1；##=#  2.=9  2和缺失值转化为9；#/#=#  1/5=9  1到5转化为9

nonmissing=#  Nonmissing=8  所有非缺失值转化为8

missing = #    miss = 9      所有缺失值转化为9

例：recode x (1 = 2) ,gen(nx)

recode x (1 = 2) (2 = 1) ,gen(nx1)

recode x (1 2 = 3) (4/7 = 3),gen(nx1)

(1)gen edu6 = 0

replace edu6=1 if educ>0 & educ<7

replace edu6=2 if educ>6 & educ<10

replace edu6=3 if educ>9 & educ<13

replace edu6=4 if educ>12 & educ<17

replace edu6=5 if educ>16 & educ<.

Tabulate edu6



(2)autocode(x,n,xmin,xmax) //把x自动分成间隔相等的n组

group(x) //等距分成x组

recode(x,x1,x2,...,xn) //当x缺失时求得缺失值


*gen exper1 = autocode(exper,5,1,51)

tabulate exper1

*sort exper

gen exper2 = group(5)

tabulate exper2

*gen exper3 = recode (exper,5,15,25,40,51)

tabulate exper3

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

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

Stata连享会   计量专题 || 精品课程 || 简书推文 || 公众号合集

 

连享会内生性专题……

 

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

1. 问题背景

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

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

2. 什么是因子变量

因子变量（Factor Variable）是对现有变量的延伸，是从类别变量中生成虚拟变量、设定类别变量之间的交乘项、类别变量与连续型变量之间的交乘项或连续变量之间的交乘项（或多项式）。在Stata中的大多数回归命令和回归后的估计命令中都可以使用这些因子变量^[1]。

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

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

以研究妇女工资的决定因素为例，使用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。则传统的做法如下^[3]：

. 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的估计系数，而不希望生成或存储这个变量^[4]。Stata中的因子变量语法完全注意到了这个问题。使用因子变量的标准做法如下：

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

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

. reg wage i.race

回归结果如下所示：

(部分回归结果省略)
------------------------------------------------------------------------------
        wage |      Coef.   Std. Err.      t    P>|t|     [95% Conf. Interval]
-------------+----------------------------------------------------------------
        race |
      black  |     -1.238      0.276    -4.48   0.000       -1.781      -0.696
      other  |      0.468      1.133     0.41   0.680       -1.754       2.690
             |
       _cons |      8.083      0.142    57.06   0.000        7.805       8.361
------------------------------------------------------------------------------

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

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

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

在实证分析中，变量的交乘项或高次项往往是重要的解释变量。以研究妇女工资的决定因素为例，若我们想在模型中加入黑人每周工作时数的变量，则传统做法是预先生成一个虚拟变量 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，相应的命令为：

. 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：邹氏检验 (Chow test)

由于不同组别之间可能会存在差异（截距项或斜率项存在差异），因此，我们需要检验这些差异在统计上是否显著。这时，我们可以使用邹氏检验。以研究妇女工资的影响因素为例，我们可以使用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_{it}={\alpha }_i+{\lambda }_t+{x_{it}}^{{}^{\prime }}\beta +{\epsilon }_{it}\text{\hspace{1em}(1)}$$
上式等价于
$$y_{it}=\sum _{i=1}^N{D}_i\cdot {\alpha }_i+\sum _{t=2}^T{W}_t\cdot {\lambda }_t+{x_{it}}^{{}^{\prime }}\beta +{\epsilon }_{it}\text{\hspace{1em}(2)}$$
其中，$D_i$ 是对应于每家公司的虚拟变量，对于公司 $i$，$D_i=1$，否则 $D_i=0$。对于一个有 $N=1000$ 家公司的面板数据而言，共有 1000 个虚拟变量。$W_t$ 是 $T-1$ 个年度虚拟变量^[5]， 定义方式与 $D_i$ 相似。

若使用因子变量，则语法很简单^[6]：

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

其中，i.id 表示 N-1 个公司虚拟变量（Stata 会默认将一家公司作为基准组），即(1)式中的 ${\alpha }_i$。

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

. xtreg y x1 x2 x3 i.year, fe

3.3 范例3：倍分法 DID 模型

在实证分析中，常常需要分析政策实施后带来的效果。这时，我们就需要采集实验组与控制组两组样本（实验组的样本代表政策实施前后的情况，控制组的样本代表不实施政策的情况），再把这两组样本合并为一份数据后进行回归分析。在建模时，我们需要设定一个是否为实验组或控制组的处理虚拟变量 Treat 放入模型，还需要设定一个反映政策实施时间前与实施后的时间虚拟变量 Time 放入模型，此外，还需要在模型中加入 Treat 变量与 Time 变量的交乘项，该交乘项的系数估计值 $a_3$ 就是政策实施后带来的效果，即实验组在政策实施后与假设未实施政策的差异。回归模型如下式(3)所示：
$$y=a_0+a_{1}Treat+a_{2}Time+a_{3}Treat\times Time+x^{{}^{\prime }}\beta +\epsilon \text{\hspace{1em}(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_{it}=a_0+{x_{it}}^{{}^{\prime }}\cdot \gamma +{p_{it}}^{{}^{\prime }}\cdot \eta +x_{it}^{{}^{\prime }}\cdot p_{it}\cdot \beta +d_{it}^{{}^{\prime }}\cdot \delta +d_{it}^{{}^{\prime }}\cdot p_{it}\cdot {\alpha }_t+{\epsilon }_{it})\text{\hspace{1em}(4)}$$

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

. 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.   Std. Err.      t    P>|t|     [95% Conf. Interval]
------------------+----------------------------------------------------------------
------------------|----------------------部分回归结果省略----------------------------
                  |
    disabled#year |
     Disabl#1988  |  -.6873656   .7197987    -0.95   0.340    -2.098154    .7234227
     Disabl#1989  |  -.5385104   .7056609    -0.76   0.445    -1.921589    .8445682
     Disabl#1990  |  -2.289956   .6945586    -3.30   0.001    -3.651274   -.9286372
     Disabl#1991  |  -2.158817   .6916205    -3.12   0.002    -3.514376    -.803257
     Disabl#1992  |  -1.386539   .6873483    -2.02   0.044    -2.733725   -.0393526
     Disabl#1993  |  -2.743325   .6953003    -3.95   0.000    -4.106097   -1.380553
     Disabl#1994  |  -3.917873   .7198205    -5.44   0.000    -5.328704   -2.507042
     Disabl#1995  |  -3.702671    .748066    -4.95   0.000    -5.168863    -2.23648
     Disabl#1996  |  -4.472484   .7373381    -6.07   0.000    -5.917649   -3.027319
                  |
            _cons |   41.05231   .4546278    90.30   0.000     40.16125    41.94337
-----------------------------------------------------------------------------------

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

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

$$lny=\alpha +\sum _{k=1}^k{\beta }_k\cdot ln{x}_k+\frac{1}{2}\sum _{k=1}^k\sum _{m=1}^m{\gamma }_{km}\cdot ln{x}_k\cdot ln{x}_m\text{\hspace{1em}(5)}$$

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

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

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

一些实证研究的文献中常常加入变量的交乘项以反映可能存在的调节效应，例如 Faulkender and Wang(2006) 、戴魁早和刘友金(2016) 与 张苏和高扬(2012) ^[9]，此时最重要的是分析 x 对 y 的边际效应。

$$y={\alpha }_1+x{\beta }_1+z{\beta }_2+(x\cdot z){\beta }_3+\epsilon \text{\hspace{1em}(6)}$$

x 对 y 的边际影响为：

$$\frac{\partial y}{\partial x}={\beta }_1+z{\beta }_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|     [95% Conf. Interval]
--------------------+----------------------------------------------------------------
               race |
             black  |     -1.442      0.297    -4.85   0.000       -2.026      -0.859
             other  |     -0.380      1.348    -0.28   0.778       -3.024       2.263
                    |
           collgrad |
      college grad  |      2.981      0.312     9.54   0.000        2.368       3.593
                    |
      race#collgrad |
black#college grad  |      2.502      0.676     3.70   0.000        1.177       3.827
other#college grad  |      1.678      2.297     0.73   0.465       -2.826       6.182
                    |
              _cons |      7.318      0.158    46.25   0.000        7.008       7.629
-------------------------------------------------------------------------------------

我们使用margins命令来计算种族(race)与是否大学毕业(collgrad)的交乘项的各个类别的妇女小时工资(wage)的预测边际值，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% Conf. Interval]
------------------------+----------------------------------------------------------------
          race#collgrad |
white#not college grad  |   7.318251   .1582176    46.25   0.000     7.007983     7.62852
    white#college grad  |   10.29895   .2693243    38.24   0.000     9.770797     10.8271
black#not college grad  |   5.875918   .2519298    23.32   0.000     5.381878    6.369958
    black#college grad  |   11.35861    .543853    20.89   0.000      10.2921    12.42512
other#not college grad  |   6.938195   1.338677     5.18   0.000     4.313019    9.563371
    other#college grad  |   11.59678   1.839835     6.30   0.000     7.988818    15.20474
-----------------------------------------------------------------------------------------

使用marginsplot命令将计算结果用图的形式表示，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% Conf. Interval]
-------------+----------------------------------------------------------------
        race |
      white  |   8.024276    .136558    58.76   0.000     7.756482    8.292069
      black  |   7.174578    .231424    31.00   0.000      6.72075    7.628406
      other  |   8.041653   1.110659     7.24   0.000     5.863625    10.21968
------------------------------------------------------------------------------

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

. marginsplot

图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|     [95% Conf. Interval]
------------------+----------------------------------------------------------------
           tenure |   .3015859   .0704148     4.28   0.000        .1635    .4396718
                  |
c.tenure#c.tenure |  -.0076366   .0037758    -2.02   0.043    -.0150411   -.0002321
                  |
            hours |   .0749604   .0115722     6.48   0.000     .0522669    .0976539
              age |  -.0850009   .0388991    -2.19   0.029    -.1612832   -.0087185
          married |  -.4619864   .2544385    -1.82   0.070     -.960949    .0369761
            south |  -1.283751   .2493374    -5.15   0.000     -1.77271   -.7947914
                  |
             race |
           black  |  -1.221256   .2870854    -4.25   0.000    -1.784241   -.6582719
           other  |   .3129332   1.098312     0.28   0.776    -1.840893     2.46676
                  |
            _cons |   8.183759   1.632092     5.01   0.000     4.983171    11.38435
-----------------------------------------------------------------------------------

由于在模型中加入了 tenure 的平方项，因此，tenure 对 wage 的边际效应将会受到 tenure 取值的影响。上述Stata命令对应的模型设定如下：

$$wage=\alpha +tenure{\beta }_1+(tenure\cdot tenure){\beta }_2+x_3{\beta }_3+\cdots +x_8{\beta }_8+\epsilon \text{\hspace{1em}(7)}$$

tenure 对 wage 的边际影响为：

$$\frac{\partial wage}{\partial tenure}={\beta }_1+2\cdot tenure\cdot {\beta }_2$$

显然，当 tenure 取值不同时，tenure 对 wage 的边际效应是不相同的。因此，需要先使用描述性统计分析的命令（如sum tenure）查看 tenure 取值的范围，然后再计算当 tenure 取不同的值所对应的边际效应，Stata命令如下：

. preserve
. keep if e(sample)
. sum tenure
. restore
. 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  |   .3015859   .0704148     4.28   0.000        .1635    .4396718
          2  |   .2863127   .0632764     4.52   0.000     .1622256    .4103998
          3  |   .2404932   .0427222     5.63   0.000     .1567136    .3242728
          4  |   .1946737   .0259414     7.50   0.000     .1438016    .2455457
          5  |   .1488542   .0233928     6.36   0.000     .1029801    .1947282
          6  |   .1030346   .0380523     2.71   0.007     .0284127    .1776566
          7  |   .0572151   .0580967     0.98   0.325    -.0567145    .1711447
          8  |   .0113956   .0795548     0.14   0.886    -.1446141    .1674053
          9  |  -.0344239   .1015342    -0.34   0.735    -.2335359    .1646881
         10  |  -.0802434   .1237575    -0.65   0.517    -.3229361    .1624492
         11  |  -.0939893   .1304534    -0.72   0.471    -.3498129    .1618343
------------------------------------------------------------------------------

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

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

图3 工作年限对妇女工资的平均边际效应

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

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

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

. sysuse "nlsw88.dta", clear
. reg wage i.race##collgrad
. est store R1
. esttab R1 using D:/Table_factor_1.csv, nogap replace

表3 1988年美国妇女工资模型结果

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

*-输出结果（不显示基准组和忽略组的系数，使用drop选项）
. esttab R1 using D:/Table_factor_2.csv, nogap      ///
	     drop(1.race 0.collgrad 1.race#0.collgrad   ///
         1.race#1.collgrad 2.race#0.collgrad        ///
         3.race#0.collgrad) replace
*-输出结果（不显示基准组和忽略组的系数，使用nobase与noomit选项）
. esttab R1 using D:/Table_factor_3.csv, nogap nobase noomit replace

表4 1988年美国妇女工资模型结果（删除缺失估计系数值的变量）

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. [PDF]
• Altunbas, Y., M.-H. Liu, P. Molyneux, R. Seth, 2000, Efficiency and Risk in Japanese Banking, Journal of Banking & Finance, 24 (10): 1605-1628. [PDF]
• Faulkender, M., R. Wang, 2006, Corporate Financial Policy and the Value of Cash, Journal of Finance, 61 (4): 1957-1990. [PDF]
• 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. [PDF]
• 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. [PDF]
• 戴魁早, 刘友金, 2016, 要素市场扭曲与创新效率——对中国高技术产业发展的经验分析, 经济研究, (7): 72-86.
• 王德祥, 李建军, 2009, 我国税收征管效率及其影响因素——基于随机前沿分析(SFA)技术的实证研究, 数量经济技术经济研究, (4): 152-160.
• 张苏, 高扬, 2012, 大学生学习行为与国家竞争力关联关系的实证研究, 管理世界, (4): 175-176.

---

注释：
[1] 详情参阅Stata帮助文件help fvvarlist。
[2] 详情参见Stata帮助文件help varlist。
[3] 为了便于说明，后续多数回归命令中都省略了控制变量。
[4] 事实上，只要你的数据中存储了race变量，我们只需要保存好dofile文件，就无需生成black这个中间变量。
[5] 之所以加入T - 1 个年度虚拟变量，是为了防止完全共线性。
[6] 此处设定i.year会自动加入 T 个年度虚拟变量，但Stata会自动删除一个，以防完全共线性。
[7] 文中实证分析所用的原始数据和相关程序可以从作者主页上下载：http://economics.mit.edu/faculty/acemoglu/data/aa2001。
[8] 例如，Altunbas, Liu, Molyneux and Seth(2000) 使用超越对数成本函数估算了日本银行的效率和风险；Wang(2007) 则使用超越对数生产函数研究30个国家R&D效率；王德祥和李建军(2009) 基于超越对数生产函数估算了我国的税收流失率。
[9] Faulkender and Wang(2006) 检验了由公司融资约束对现金持有边际市场价值的影响。戴魁早和刘友金(2016) 研究发现了要素市场扭曲对创新效率的影响存在着企业差异，企业规模在规避要素市场扭曲对创新效率的抑制效应中具有积极作用。张苏和高扬(2012) 的实证研究发现学生来源于城市或农村地区对国家竞争力的影响作用受到每周上网时间的影响，如果上网时间在每周8小时以下，城市大学生的学习行为落入“增进国家竞争力导向”上高效率区域而不是低效率区域的概率比农村大学生要高，否则要低。

 

关于我们

• 「Stata 连享会」 由中山大学连玉君老师团队创办，定期分享实证分析经验， 公众号：StataChina。
• 公众号推文同步发布于 CSDN 、简书 和 知乎Stata专栏。可在百度中搜索关键词 「Stata连享会」查看往期推文。
• 点击推文底部【阅读原文】可以查看推文中的链接并下载相关资料。
• 欢迎赐稿： 欢迎赐稿。录用稿件达 三篇 以上，即可 免费 获得一期 Stata 现场培训资格。
• E-mail： StataChina@163.com
• 往期推文：计量专题 || 精品课程 || 简书推文 || 公众号合集

---