連玉君，楊 柳

(1. 中山大學 嶺南學院，廣東 廣州 510275； 2.西北大學 城市與環(huán)境學院，陜西 西安 710127)

一、問題背景

實證分析中，我們經(jīng)常需要在模型中加入反映類別的虛擬變量，以便控制不可觀測的組間差異。而在另一些分析中，為了刻畫調(diào)節(jié)效應，尚需在模型中加入變量的交乘項或平方項。傳統(tǒng)的做法是，預先生成虛擬變量或交乘項，進而將它們加入模型。然而，當虛擬變量或交乘項的數(shù)目較多時，上述方法就顯得尤為不便，不但浪費計算機內(nèi)存，也會嚴重降低我們的工作效率。

在Stata中，我們可以使用因子變量(Factor Variable)簡化操作步驟、快捷地在回歸模型中加入虛擬變量、交乘項、平方項或高次項。更為重要的是，由于引入交乘項或平方項后，解釋變量對被解釋變量的邊際影響不再是常數(shù)，而是某個變量(調(diào)節(jié)變量)的函數(shù)，在有些模型設(shè)定下，這種關(guān)系可能是非線性的。此時，若使用因子變量，并配合Stata中的margins和marginsplot命令，可以非常便捷、直觀地分析關(guān)鍵變量的邊際效應。

本文將以文獻中廣泛使用的一些模型設(shè)定形式為例，來說明Stata中因子變量的使用方法。后文結(jié)構(gòu)安排如下：第二部分介紹因子變量的基本用法，包括使用因子變量表示虛擬變量、交乘項和高次項；第三部分舉例說明因子變量在常用回歸模型中的應用；第四部分重點介紹如何使用因子變量進行邊際效應分析。

二、什么是因子變量

因子變量(Factor Variable)是對現(xiàn)有變量的延伸，是從類別變量中生成虛擬變量、設(shè)定類別變量之間的交乘項、類別變量與連續(xù)型變量之間的交乘項或連續(xù)變量之間的交乘項(或多項式)。在Stata中的大多數(shù)回歸命令和回歸后的估計命令中都可以使用這些因子變量。①

因子變量的五種運算符及其含義如表1所示：

表1 因子變量的運算符及含義

注：(a) 由因子變量運算符生成的指示變量和交乘項是實際存在的變量，它們與數(shù)據(jù)表中的變量是一樣的，但是在數(shù)據(jù)表中并不顯示出來；(b) 類別變量的值必須是非負整數(shù)，范圍介于0至32740之間；(c) 因子變量的運算符有時可與時間序列的運算符L.和F.組合在一起使用.②

以研究婦女工資的決定因素為例，使用Stata軟件自帶的數(shù)據(jù)文件(nlsw88.dta)。該數(shù)據(jù)包含了1988年采集的2246個美國婦女的資料，包括：小時工資wage、每周工作時數(shù)hours、種族race、職業(yè)occupation、年齡age、是否大學畢業(yè)collgrad、當前職業(yè)的工作年限tenure、是否結(jié)婚married、是否居住在南部地區(qū)south、合計工作年限ttl_exp等變量。其中，小時工資wage、每周工作時數(shù)hours、年齡age、當前職業(yè)的工作年限tenure、合計工作年限ttl_exp為連續(xù)型變量；種族race為類別變量(1代表白種人white，2代表黑種人black，3代表其他人種other)、職業(yè)occupation為類別變量(13個職業(yè)類別)；是否大學畢業(yè)collgrad、是否結(jié)婚married、是否居住在南部地區(qū)south為虛擬變量。

在這份數(shù)據(jù)中有一個表示種族的類別變量race，取值為1、2、3，分別對應“白人”、“黑人”和“其他人種”。假設(shè)我們想在模型中加入一個反映種族的虛擬變量black，當某個婦女是黑人時，black取值為1，否則為零。則傳統(tǒng)的做法如下：③

L1

. gen black=1

. replace black=0 if race!=2

. reg wage black

若延續(xù)這一思路，但使用因子變量來生成black變量，則命令為：

L2

. gen black = 2.race

只需要一條命令，而且含義非常明確。這里“2.race”本質(zhì)上是一個條件判斷語句：判斷某一行觀察值中的race變量取值是否為2，若是，則返回1到變量black中，否則返回0。

然而，在多數(shù)情況下，我們的目的只是希望得到虛擬變量black的估計系數(shù)，而不希望生成或存儲這個變量。④Stata中的因子變量語法完全注意到了這個問題。使用因子變量的標準做法如下：

L3

. sysuse "nlsw88.dta", clear

. reg wage 2.race

注意，我們無需預先生成black變量，而是直接在回歸模型中加入了2.race因子變量。

有些讀者注意到，race變量有三個取值，因此，我們可以在模型中放入兩個虛擬變量，此時可以書寫如下命令：

L4

. reg wage i.race

回歸結(jié)果如圖1所示：

圖1 模型回歸結(jié)果

從圖1的結(jié)果中可以看到，回歸模型中加入了black虛擬變量和other虛擬變量，分別對應race變量的第二個和第三個類別，而第一個類別white被Stata默認作為基準組，目的在于防止完全共線性。black變量的系數(shù)值為-1.238，表示黑人的平均工資比白人低1.238個單位，并在統(tǒng)計上顯著；other變量的系數(shù)值為0.468，表示其他人種的平均工資比白人高0.468個單位，但是在統(tǒng)計上不顯著。

在實證分析中，有時會根據(jù)研究內(nèi)容的需要改變基準組的設(shè)定，此時可使用ib.或b.的前綴，具體寫法見表2中描述。

在實證分析中，變量的交乘項或高次項往往是重要的解釋變量。以研究婦女工資的決定因素為例，若我們想在模型中加入黑人每周工作時數(shù)的變量，則傳統(tǒng)做法是預先生成一個虛擬變量(black)代表是否黑人，再生成一個新變量(black_x_hours)表示黑人與每周工作時數(shù)的交乘項，然后再將這個新變量放入回歸模型中。Stata命令為：

L5

. gen black=1

. replace black=0 if race!=2

. gen black_x_hours =black * hours

. reg wage black_x_hours

若在Stata中使用因子變量實現(xiàn)上述過程，則命令十分簡潔：

L6

. reg wage 2.race#c.hours

表2 基準組的運算符及含義

注：(a) i可以省略不寫。例如：ib2.group與b2.group的寫法等價；(b) ib(#2)指使用變量值中的第二位排序的值所對應的類別作為基準組；(c) 如果想在線性回歸模型中將group變量中值為3的類別設(shè)置為基準組，則命令可寫為：reg y i.sex ib3.group.

下面，我們介紹如何在Stata中使用因子變量表示變量的交乘項或高次項。我們以研究婦女工資的決定因素為例進行說明。

1.兩個類別變量的交乘項

在回歸模型中加入種族和職業(yè)類別的交乘項，Stata命令為：

L7

. reg wage i.race#i.occupation

若在回歸模型中既要放入種族與職業(yè)的虛擬變量，又需要同時放入這兩個變量的交乘項，則在回歸命令中使用“i.race##i.occupation”，相應的命令為：

L8

. reg wage i.race##i.occupation

2.類別變量與連續(xù)變量的交乘項

在回歸模型中加入種族和每周工作時數(shù)的交乘項，Stata命令為：

L9

. reg wage i.race#c.hours

需要注意的是，在上例中，由于我們把hours變量視為連續(xù)變量，因此，需要在其前面加上c.符號以便告知Stata該變量是連續(xù)變量。

3.連續(xù)變量與連續(xù)變量的交乘項(高次項)

在回歸模型中加入年齡age變量，以及其平方項，Stata命令如下：

L10

. reg wage c.age##c.age

上述命令中c.age表示年齡age變量被當成連續(xù)型變量。如果我們在Stata命令中使用i.age，則年齡age變量被當成類別變量處理，此時，類別的個數(shù)為年齡age變量中不同取值的個數(shù)。

三、常用回歸模型的因子變量表述

范例1：鄒氏檢驗

由于不同組別之間可能會存在差異(截距項或斜率項存在差異)，因此，我們需要檢驗這些差異在統(tǒng)計上是否顯著。這時，我們可以使用鄒氏檢驗。以研究婦女工資的影響因素為例，我們可以使用chowtest命令檢驗工會成員與非工會成員兩個樣本組中工資影響因素是否存在差異(或稱之為存在結(jié)構(gòu)變化)：

L11

. global xx "hours age tenure ttl_exp married"

. chowtest wage $xx, group(union) detail

事實上，我們也可以使用因子變量的語法，在回歸模型中加入分組變量與其他控制變量的交乘項，然后再聯(lián)合檢驗分組變量的系數(shù)以及所有交乘項的系數(shù)是否都等于0。此時，即使我們不使用chowtest命令，也可以輕松實現(xiàn)鄒氏檢驗：

L12

. global xx "hours age tenure ttl_exp married"

. reg wage $xx i.unioni.union#c.($xx)

. testparmi.unioni.union#c.($xx)

范例2：雙向固定效應模型

在實證分析過程中，經(jīng)常需要在模型中加入反映年度、公司或行業(yè)特征的虛擬變量。當虛擬變量的數(shù)目眾多時，采用手動輸入變量的方式會非常耗時。例如，下述模型(1)是文獻中廣泛應用的“雙向固定效應模型”：

(1)

上式等價于

(2)

其中，Di是對應于每家公司的虛擬變量，對于公司i，Di=1，否則Di=0。對于一個N=1000家公司的面板數(shù)據(jù)而言，共有1000個虛擬變量。Wt是T-1個年度虛擬變量，⑤定義方式與Di相似。

若使用因子變量，則語法很簡單：⑥

L13

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

其中，i.id表示N個公司虛擬變量，即(2)式中的αi。

當然，我們也可以用xtreg命令自帶的fe選項來控制個體效應αi，同時使用因子變量來加入年度虛擬變量λt，命令如下：

L14

. xtreg y x1 x2 x3 i.year, fe

范例3：DID模型

在實證分析中，常常需要分析政策實施后帶來的效果。這時，我們就需要采集實驗組與控制組兩組樣本(實驗組的樣本代表政策實施前后的情況，控制組的樣本代表不實施政策的情況)，再把這兩組樣本合并為一份數(shù)據(jù)后進行回歸分析。在建模時，我們需要設(shè)定一個是否為實驗組或控制組的處理虛擬變量放入模型，還需要設(shè)定一個反映政策實施時間前與實施后的時間虛擬變量放入模型，此外，還需要在模型中加入變量與變量的交乘項，該交乘項的系數(shù)估計值就是政策實施后帶來的效果，即實驗組在政策實施后相對于政策實施前的差異。回歸模型如式(3)所示：

y=α0+α1Treat+α2Time+α3Treat×Time+xβ+ε

(3)

在Stata中的命令寫法如下：

L15

. reg Treat Time Treat#Time x1 x2 x3

*-或者寫為

. reg Treat##Time x1 x2 x3

在多期DID分析中，我們常常需要加入年度虛擬變量與處理變量的交乘項來檢驗“共同趨勢假設(shè)(common trend)”以及政策效果。例如，在Acemoglu and Angrist(2001)文中，作者采集了1988年至1997年的人口調(diào)查數(shù)據(jù)，將樣本分為殘疾人與非殘疾人，使用多期DID模型研究了1992年美國對殘疾人工作保護法案(ADA)的實施效果。文中模型設(shè)定如下：

其中，i為殘疾人個體；t為年份；yit為工作周數(shù)或周平均工資；xi為一系列控制變量，πt為時間趨勢項；di為是否殘疾，其系數(shù)為δ；αt為隨年份變化的實施ADA與未實施ADA對殘疾人工作保護效應的系數(shù)。當t≥1992年時，ai量化了實施ADA之后對殘疾人工作保護的效應(使用非殘疾人作為控制組)；當t<1992年時，ai檢驗了實施ADA之前的各年份中殘疾人和非殘疾人的工作周數(shù)或周平均工資是否具有統(tǒng)計上顯著的差異，相當于同時進行了DID模型的“共同趨勢假設(shè)”的檢驗。使用如下命令即可得到Acemoglu、Angrist(2001)文中表2第(1)列中的結(jié)果。⑦

L16

. use "ABA_JPE2001.dta", clear

. global controls "i.age_Gi.edu_Gi.race_Gi.region"

. reg y i.year##($controls) Treat i.yeari.Treat#i.year

范例4：超越對數(shù)生產(chǎn)函數(shù)

很多文獻使用超越對數(shù)生產(chǎn)函數(shù)估計生產(chǎn)效率，模型中會包含投入要素K和L的高階項，包括平方項和二者的交乘項⑧。例如在Kumbhakar(1989)中，作者通過在模型中加入二次項來捕捉要素的交互影響和潛在的非線性關(guān)系，模型如式(5)所示。

此時，使用因子變量會讓Stata中的命令變得異常簡潔，如下所示：

L17

. webuse frontier1.dta, clear

L18

. global y "lnoutput" // y

. global x "lnlaborlncapital" // x1 x2

*-超越對數(shù)生產(chǎn)函數(shù)

. sfcross $y c.($x)##c.($x) // Eq.

(5)

四、邊際效應分析和圖形化呈現(xiàn)

一些實證研究的文獻中常常加入變量的交乘項以反映可能存在的調(diào)節(jié)效應，例如Faulkender and Wang(2006)，戴魁早與劉友金(2016)，張?zhí)K與高揚(2012)⑨，此時最重要的是分析x對y的邊際效應。

y=α+xβ1+(x·z)β2+ε

(6)

x對y的邊際影響為：

顯然，?y/?x取決于z的取值。這種邊際效應分析若是采用手動計算會非常煩瑣，而后續(xù)的圖形呈現(xiàn)則更為復雜。然而，若是借助因子變量，并配合Stata中的margins和marginsplot命令，對于包含交乘項的模型中的邊際效應的分析和圖形化展示都變得異常輕松。

下面，我們以研究婦女工資的決定因素為例進行說明。使用Stata軟件的自帶數(shù)據(jù)nlsw88.dta(1988年美國婦女小時工資)，以wage(婦女的小時工資)作為被解釋變量，以race(種族類別)、collgrad(是否大學畢業(yè))、race與collgrad的交乘項作為解釋變量建立線性回歸模型，Stata中的命令如下：

L19

. sysuse"nlsw88.dta", clear

. reg wage i.race##collgrad

回歸結(jié)果如圖2所示：

圖2 模型回歸結(jié)果

我們使用margins命令來計算種族(race)與是否大學畢業(yè)(collgrad)的交乘項的各個類別對婦女小時工資(wage)的平均邊際效應，Stata命令和結(jié)果如下所示：

L20

. margins i.race#collgrad

計算結(jié)果如圖3所示：

圖3 邊際效應計算結(jié)果

使用marginsplot命令將計算結(jié)果用圖4的形式表示，Stata命令和結(jié)果如下所示：

L21

. marginsplot

圖4 種族與是否大學畢業(yè)的交乘項對婦女工資的平均邊際效應

從圖4中可以直觀地看到，不論種族類別，大學畢業(yè)的婦女的平均工資高于非大學畢業(yè)的婦女，該結(jié)果符合我們的一般認知；另外，我們從圖中還發(fā)現(xiàn)一個有趣的結(jié)果，當婦女為非大學畢業(yè)時，白人的平均工資高于黑人，而當婦女為大學畢業(yè)時，白人的平均工資低于黑人。

下面，我們使用margins命令附加atmeans選項來計算當其他變量取均值時不同種族類別對婦女工資的平均邊際效應，Stata命令和結(jié)果如下所示：

L22

. margins i.race, atmeans

*-或者寫為：

. margins race, atmeans

計算結(jié)果如圖5所示。

圖5 邊際效應計算結(jié)果

使用marginsplot命令將計算結(jié)果用圖6的形式表示，Stata命令和結(jié)果如下所示：

L23

. marginsplot

圖6 種族對婦女工資的平均邊際效應

從圖6中可以直觀地看到，白人與其他種族婦女的平均工資高于黑人婦女。

兩個連續(xù)變量的交乘項對被解釋變量的邊際效應也可以使用margins命令來計算。我們?nèi)砸匝芯繈D女工資的決定因素為例進行說明。在回歸模型中加入ttl_exp(合計工作年限)及其平方項，并將hours(每周工作時數(shù))、age(婦女年齡)、tenure(當前職業(yè)的工作年限)、married(是否結(jié)婚)、south(是否居住在南部地區(qū))、race(種族類別)作為控制變量，Stata命令如下所示：

L24

. global xx "hours age tenure married south i.race"

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

由于在模型中加入了ttl_exp的平方項，因此，ttl_exp對wage的邊際效應將會受到ttl_exp取值的影響。上述Stata命令對應的模型設(shè)定如下：

wage=α+ttl_expβ1+(ttl_exp·ttl_exp)β2+β3x3+…+β8x8+ε

(7)

ttl_exp對wage的邊際影響為：

顯然，當ttl_exp取值不同時，ttl_exp對wage的邊際效應是不相同的。因此，需要先使用描述性統(tǒng)計分析的命令(如sum ttl_exp)查看ttl_exp取值的范圍，然后再計算當ttl_exp取不同的值所對應的邊際效應，Stata命令如下：

L25

. sum ttl_exp

. margins, dydx(ttl_exp) at(ttl_exp=(0.1 1(3)28 28.9))

計算結(jié)果見附錄1。

使用marginsplot命令將計算結(jié)果用圖7的形式表示，Stata命令和結(jié)果如下所示：

L26

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

圖7 合計工作年限對婦女工資的平均邊際效應

Stata中的多數(shù)命令都支持margins和marginsplot命令。因此，即使對于非線性模型，如logit, tobit等，我們?nèi)匀豢梢越柚@兩個命令很方便地分析邊際效應。

五、輸出回歸結(jié)果時的問題及解決辦法

以研究婦女工資的決定因素為例。wage(婦女的小時工資)作為回歸模型的被解釋變量，race(種族類別)、collgrad(是否大學畢業(yè))、race與collgrad的交乘項作為解釋變量，并將回歸結(jié)果輸出到Excel中，Stata命令和結(jié)果如下：

L27

. reg wage i.race##collgrad

. est store R1

. esttab R1 using D:/Table_factor_1.csv, nogap replace

表3 1988年美國婦女工資模型結(jié)果

注：括號內(nèi)為t統(tǒng)計量；* p<0.05, ** p<0.01, *** p<0.001.

我們發(fā)現(xiàn)在表3中有很多變量的系數(shù)值為0，并缺失t統(tǒng)計量。造成這種情況的原因有二：一方面是由于有些虛擬變量作為基準組，例如：1.race，0.collgrad與1.race#0.collgrad，Stata默認將它們作為基準組，所以就缺失這些基準組的估計系數(shù)值和標準誤；另一方面是由于有些交乘項中的其中一個因子變量是基準組，其他變量與這個作為基準組的因子變量交乘后的交乘項就被忽略了，所以其估計系數(shù)值和標準誤就會缺失，例如：1.race#1.collgrad，2.race#0.collgrad與3.race#0.collgrad。此時，可以使用esttab命令的drop()選項來屏蔽這些系數(shù)的顯示，還可以使用nobase和noomit的選項，Stata命令和結(jié)果如下所示：

L28

*-輸出結(jié)果(不顯示基準組和忽略組的系數(shù))

. esttabR1 using D:/Table_factor_2.csv, nogap///

drop(1.race 0.married 1.race#0.married ///

1.race#1.married 2.race#0.married ///

3.race#0.married) replace

*-輸出結(jié)果(不顯示基準組和忽略組的系數(shù)，使用nobase與noomit的選項)

. esttabR1 using D:/Table_factor_3.csv, nogapnobasenoomit replace

刪除缺失估計系數(shù)值后，模型結(jié)果呈表4所示。

表4 1988年美國婦女工資模型結(jié)果(刪除缺失估計系數(shù)值的變量)

注：括號內(nèi)為t統(tǒng)計量；* p<0.05, ** p<0.01, *** p<0.001.

六、小 結(jié)

本文介紹了Stata中因子變量產(chǎn)生虛擬變量與交乘項的使用方法，以常用經(jīng)典回歸模型為例，提供了它們的Stata命令，并進一步提供了因子變量與margins及marginsplot命令相配合分析邊際效應的示例。在Stata中的大多數(shù)命令中都可以使用因子變量的表述方法。該方法可以使Stata命令更加簡潔，并能夠大幅度提高實證分析的效率，但需要注意分析使用因子變量表述方法后得到的模型設(shè)定結(jié)構(gòu)。

注釋：

①詳情參閱Stata幫助文件(help fvvarlist).

②詳情參見Stata幫助文件help varlist.

③為了便于說明，后續(xù)多數(shù)回歸命令中都省略了控制變量.

④事實上，只要你的數(shù)據(jù)中存儲了race變量，我們只需要保存好dofile，就無需生成black這個中間比變量.

⑤之所以加入T-1個年度虛擬變量，是為了防止完全共線性.

⑥此處設(shè)定i.year會自動加入T個年度虛擬變量，但Stata會自動刪除一個，以防完全共線性.

⑦文中實證分析所用的原始數(shù)據(jù)和相關(guān)程序可以從作者主頁上下載：http://economics.mit.edu/faculty/acemoglu/data/aa2001.

⑧例如，Altunbas, Liu, Molyneux and Seth(2000)使用超越對數(shù)成本函數(shù)估算了日本銀行的效率和風險；Wang(2007)則使用超越對數(shù)生產(chǎn)函數(shù)研究30個國家R&D效率；王德祥、李建軍(2009)基于超越對數(shù)生產(chǎn)函數(shù)估算了我國的稅收流失率.

⑨Faulkender and Wang(2006)檢驗了由公司融資約束對現(xiàn)金持有邊際市場價值的影響。戴魁早和劉友金(2016)研究發(fā)現(xiàn)了要素市場扭曲對創(chuàng)新效率的影響存在著企業(yè)差異，企業(yè)規(guī)模在規(guī)避要素市場扭曲對創(chuàng)新效率的抑制效應中具有積極作用。張?zhí)K和高揚(2012)的實證研究發(fā)現(xiàn)了學生來源于城市或農(nóng)村地區(qū)對國家競爭力的影響作用受到每周上網(wǎng)時間的影響，如果上網(wǎng)時間在每周8小時以下時，城市大學生的學習行為落入“增進國家競爭力導向”上高效率區(qū)域而不是低效率區(qū)域的概率比農(nóng)村大學生要高，否則要低.

參考文獻：

[1]Acemoglu D, Joshua D Angrist.Consequences of employment protection? The case of the americans with disabilities act[J].Journal of Political Economy,2001,109 (5): 915-957.

[2]Altunbas Y, M-H Liu, P Molyneux, R Seth.Efficiency and risk in japanese banking[J].Journal of Banking & Finance,2000,24 (10): 1605-1628.

[3]Faulkender M,R Wang.Corporate financial policy and the value of cash[J].Journal of Finance,2006，61 (4): 1957-1990.

[4]Kumbhakar S C.Estimation of technical efficiency using flexible functional forms and panel data[J].Journal of Business and Economic Statistics,1989，(7):253-258.

[5]Wang E C.R&d efficiency and economic performance: A cross-country analysis using the stochastic[J].frontier approach,Journal of Policy Modeling,2007，29 (2): 345-360.

[6]戴魁早, 劉友金.要素市場扭曲與創(chuàng)新效率——對中國高技術(shù)產(chǎn)業(yè)發(fā)展的經(jīng)驗分析[J].經(jīng)濟研究,2016，(7): 72-86.

[7]連玉君，廖俊平.如何檢驗分組回歸后的組間系數(shù)差異[J].鄭州航空工業(yè)管理學院學報，2017，(6)：97-109.

[8]王德祥, 李建軍.我國稅收征管效率及其影響因素——基于隨機前沿分析(sfa)技術(shù)的實證研究[J].數(shù)量經(jīng)濟技術(shù)經(jīng)濟研究,2009，(4): 152-160.

[9]張 蘇,高 揚.大學生學習行為與國家競爭力關(guān)聯(lián)關(guān)系的實證研究[J].管理世界，2012，(4): 175-176.

附錄1邊際效應計算結(jié)果