分組回歸后的組間系數差異

問題：實證分析中，經常需要對比分析兩個子樣本組的系數是否存在差異。
例如，在公司金融領域，研究薪酬激勵是否有助於提升業績時，模型設定為：

關注的重點是系數  。
我們經常把樣本組分成“國有企業(SOE)”和“民營企業(PRI)”兩個樣本組，繼而比較  和  是否存在差異。通常認為，民營企業的薪酬激勵更有效果，即  。

如果兩個樣本組中的模型設定是相同的，則兩組之間的系數大小是可以比較的，而且這種比較在多數實證分析中都是非常必要的。

下面我們介紹三種檢驗組間系數差異的方法：

1. 引入交叉項（Chow 檢驗）

2. 基於似無相關模型的檢驗方法 (suest)

3. 費舍爾組合檢驗（Permutation test）

方法 1: 引入交叉項

這是文獻中最常用的方法，執行起來也最簡單。以檢驗 ttl_exp 在兩組之間的系數是否存在顯著差異為例。引入一個虛擬變量 ，若某個婦女是黑人，則  ，否則 。在如下命令中，black 變量即為這里的  。模型設定為：

這是最基本的包含虛擬變量，以及虛擬變量與一個連續變量交乘項的情形。

顯然，對於白人組而言，  ，則 (1) 式可以寫為：

 對於黑人組， (1) 式可以寫為：

 參數  和  分別反映了黑人組相對於白人組的截距和斜率差異。我們關注的是參數  ，它反映了 ttl_exp 這個變量在兩個樣本組中的系數差異。因此，檢驗 ttl_exp 在兩組之間的系數是否存在顯著差異就轉變為 。相應的估計命令如下：

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   //全樣本回歸+交乘

 

交乘項 [ttl_x_black] 的系數為 -.01818, 對應的 p-value 為 0.756，表明 [ttl_exp] 的系數在兩組之間並不存在顯著差異。

我們也可以不事先生成交乘項，而直接采用 stata 的因子變量表達式，得到完全相同的結果：

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

然而，需要特別強調的是，在上述檢驗過程中，我們無意識中施加了一個非常嚴格的假設條件：只允許變量 [ttl_exp] 的系數在兩組之間存在差異，而其他控制變量(如 married, south, hours 等) 的系數則不隨組別發生變化。

這顯然是一個非常嚴格的假設。因為，從 -Table 1- 的結果來看, married， south, hours 等變量在兩組之間的差異都比較明顯。

為此，我們放松上述假設，允許所有的變量在兩組之間都存在系數差異（注意：所有離散變量前都要加 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)    

方法 2: SUEST (基於似無相關模型SUR的檢驗)

顧名思義，所謂的似無相關模型（seemingly unrelated regression）其實就是表面上看起來沒有關系，但實質上有關系的兩個模型。這聽起來有點匪夷所思。這種“實質上”的關系其實是假設白人組和黑人組的干擾項彼此相關。為了表述方便，將白人和黑人組的模型簡寫如下：

若假設  ，則我們可以分別對白人組和黑人組進行 OLS 估計。

然而，雖然白人和黑人種族不同，但所處的社會和法律環境，面臨的勞動法規都有諸多相似之處，使得二者的干擾項可能相關，即 。此時，對兩個樣本組執行聯合估計（GLS）會更有效率（詳見 Greene (2012, Econometric analysis, 7th ed, 292–304)）。

執行完 SUR 估計后，我們就可以對兩組之間的系數差異進行檢驗了。

從上面的原理介紹，可以看出，基於 SUR 估計進行組間系數差異檢驗時，假設條件比第一種方法要寬松一些：

1. 其一，在估計過程中，並未預先限定白人組和黑人組各個變量的系數一定要相同，因此在 (2) 式中，我們分別用  和  表示白人組和黑人組各個變量的系數向量；

2. 其二，兩個組的干擾項可以有不同的分布，即 可以不同，即  ,  ，且允許二者的干擾項相關，  。

stata 中執行上述檢驗的步驟為：

Step 1: 分別針對白人組和黑人組進行估計（不限於OLS估計，可以執行 Logit, Tobit 等估計），存儲估計結果；

Step 2：使用 suest 命令執行 SUR 估計；

Step 3: 使用 test 命令檢驗組間系數差異。

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

對上述命令和結果的簡要解釋如下：

1. 白人組和黑人組的估計結果分別存儲於 w 和 b 兩個臨時性文件中；

2. 執行 - suest w b - 命令時，白人組和黑人組的被視為兩個方程，即文的 (2a) 和 (2b) 式。Stata 會自動將兩個方程對應的樣本聯合起來，采用 GLS 執行似無相關估計（SUR）；

3. 由於 SUR 屬於多方程模型，因此需要指定每個方程的名稱，在下面呈現的回歸結果中，[w_mean] 和 [b_mean] 分別是白人組和黑人組各自對應的方程名稱。因此，[w_mean]ttl_exp 表示白人組方程中 ttl_exp 變量的系數，而 [b_mean]ttl_exp 則表示黑人組中 ttl_exp 變量的系數。

使用biff命令

上述過程可以使用我編寫的 - bdiff - 命令非常快捷的加以實現，結果的輸出方式也更為清晰（在 stata 命令窗口中輸入 - ssc install bdiff, replace- 可以下載最新版命令包，進而輸入 - help bdiff - 查看幫助文件）：

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

幾點說明：

• 使用 -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。

面板數據的處理方法

- 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 - 中的類似范例。

新方法: 安裝新版 suest 命令

目前，可以使用 Federico Belotti. 更新后的 suest.ado 文檔替換 Stata 官方提供的 suest.ado 文檔。前者支持 xtreg 命令。

替換方法為：

• Step 1： 執行 net install suest, replace 命令，suest.ado 文件被自動安裝在 stata15\ado\plus\s 文件夾中。
• Step 2： 用 stata15\ado\plus\s 文件夾中的 suest.ado 文件替換掉 stata15\ado\base\s 文件夾中的同名文件。

然后就可以在完成 xtreg …… 估計后，使用 suest 命令進行組間系數差異檢驗了。

help suest xtreg 出現這個ado文件， 旁邊就有 click here to install

方法 3：費舍爾組合檢驗 （Fisher's Permutation test)

A. 基本思想

將二者的系數差異定義為 ，檢驗的原假設為：  。

我們仍然關注 ttl_exp 變量在兩組之間的系數差異。以 -Table 1- 中的結果為例，可以看到，ttl_exp 變量在 [white 組] 和 [black 組] 的系數估計值分別為 0.251 (  ) 和 0.269 (  )，因此，實際觀察到的系數差異為  。

這里，  是一個統計量，若能知道其分布特征，便可通過分析  在  的分布中的相對位置來判斷我們實際觀察到  的概率。若概率很小，則表明  是小概率事件，此時拒絕原假設，反之則無法拒絕原假設。

例如，若假設  服從標准正態分布，即  ，則基於實際觀察到的  ，我們很容易得出結論：無法拒絕原假設，即兩組之間的 ttl_exp 的系數不存在顯著差異。p-value 很容易計算 (當然，也可以查表得到)：

. dis normal(-0.018)    //單尾檢驗
.49281943

我們並不知道 d 的分布特征。此時，可以對現有樣本進行重新抽樣，以得到經驗樣本 (empirical sample)，進而利用經驗樣本構造出組間系數差異統計量 d 的經驗分布 (empirical distribution)，從而最終得到經驗 p 值 (empirical p-value)。

下面先通過一個小例子說明 “經驗 p 值” 和 “經驗分布” 的概念，進而介紹使用組合檢驗獲得 “經驗 p 值” 的流程。

B. 經驗P值

在這個小例子中，我們先隨機生成一個服從標准正態分布的隨機數 d，共有 10000 個觀察值。這些觀察值是通過模擬產生的。如果這些觀察值構成的樣本是通過從原始樣本（原始樣本是從母體中一次隨機抽樣，稱為 “抽樣樣本，sample”）中二次抽樣得到的，則稱為 “經驗樣本 (empirical sample)”。

然后，我們數一下在這 10000 個隨機數中，有多個是大於  (我們實際觀察到的數值)，命令為 count if d<-0.018。一共有 4963 個觀察值大於  ，由此可得，經驗 p 值 = 4963/10000 = 0.4963。這與采用 normal() 函數得到的結果 (0.4928) 非常接近。

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 

C. 經驗樣本

上例中，我們假設 d 服從標准正態分布，從而可以通過 monte carlo 模擬的方式產生 10000 個觀察值，這事實上是構造了一個經驗樣本。但多數情況下，我們並不知道 d 的分布特征，此時無法使用 monte carlo 模擬。然而，若假設抽樣樣本 (sample) 是從母體 (population) 中隨機抽取的，則可以通過抽樣樣本中二次抽樣得到經驗樣本 (empirical sample)，這些經驗樣本也可以視為對母體的隨機抽樣。

若抽樣過程中為無放回抽樣 (sampling with no replacement)，則相應的檢驗方法稱為 “組合檢驗（permutation test）”；若抽樣過程中為有放回抽樣 (也稱為可重復抽樣，sampling with replacement)，則對應的檢驗方法稱為 “基於 Bootstrap 的檢驗”。

D. 費舍爾組合檢驗的步驟

若  是正確的，則對於任何一個婦女而言（不論她是白人還是黑人），其 x 對 y 的邊際影響都是相同的。因此，我們可以將白人組和黑人組的觀察值混合起來，從中隨機抽取 n1 個觀察值，並將其視為"白人組"，剩下的 n2 個觀察值可以視為“黑人組”。

• Step 0: 分別針對白人組和黑人組估計模型 (3a) 和 (3b)，得到系數估計值  和 ，以及二者的系數差異  ；
• Step 1: 將白人組和黑人組的樣本混合起來，得到 n1+n2 個觀察值構成的樣本 S；
• Step 2: 獲得經驗樣本 —— 從 S 中隨機抽取 (無放回) n1 個觀察值，將其視為“白人組”(記為 Sw)，剩下的 n2 個觀察值可以視為“黑人組” (記為 Sb)；
• Step 3: 分別針對經驗樣本 Sw 和 Sb，估計模型 (3a) 和 (3b)，得到  和  （上標  表示利用第一筆經驗樣本得到的估計值），以及二者的差異  ；
• Step 4: 獲得統計量  的經驗分布 —— 將 Step 2 和 Step 3 重復執行  次 (如  )，則可以得到  ，亦可簡記為  ；
• Step 5: 計算 經驗 p 值，  ，其中  表示 Step 4 中得到的  個  中大於我們實際觀測到的  的個數。若  ，則可以在5%水平上拒絕原假設，表明兩組的系數差異是顯著的。
• 需要說明的是，由於  的分布未必是對稱的，因此， 與  都可以視為在 5% 水平上拒絕原假設的證據。因為，前者意味着  在 1000 個  中屬於非常大的數值，而后者意味着它是非常小的數值。無論如何，在原假設  下觀察到  都是小概率事件，也就意味着原假設是不合理的。
• 此外，該方法的並不局限於普通的線性回歸模型(-regress-命令)，可以應用於各種模型的估計命令，如 -xtreg-, -xtabond-, -logit-, -ivregress- 等。

E. stata實現

上述過程可以使用連玉君編寫的 -bdiff- 命令來實現。在命令窗口中輸入 -ssc install bdiff, replace- 可以自動安裝該命令。幫助文件中提供了多個范例。

先使用一個簡單的例子，不考慮行業虛擬變量：

 *-數據處理
        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

• 選項 group() 中填寫用於區分組別的類別變量（若有多個組，可以預先刪除不參與比較的組，類似於上面的 drop if race==3 命令）；
• 選項 model() 用於設定回歸模型，即上面提到的模型 (3a) 或 (3b)；
• 選項 reps(#) 用於設定抽樣次數，即上文提到的  ，通常設定 1000-5000 次即可；
• 附加 detail 選項，可以進一步列表呈現兩組的實際估計系數  和  ；

上述過程大約用時 13 秒，結果如下：

 

 可以看到，ttl_exp 的經驗 p 值為 0.49，表明白人和黑人組的 ttl_exp 系數不存在顯著差異；married 變量的 p 值為 0.08，我們可以在 10% 水平上拒絕原假設。細心的讀者會發現，該變量對應的 Freq = 920，為什么？（答案在上面 Step 5 處）。

添加行業虛擬變量：若需在模型中加入虛擬變量，處理過程會稍微復雜一些。需要手動生成行業虛擬變量，並保證兩個樣本組中參與回歸的行業虛擬變量個數相同。此外，書寫命令時，不能使用通配符。

*-數據預處理(可以忽略)
  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

面板數據

若原始數據為面板數據，通常會采用 -xtreg-, -xtabond- 等考慮個體效應的方法進行估計。抽樣過程必須考慮面板數據的特征。在執行 -bdiff- 命令之前，只需設定 -xtset id year-，聲明數據為面板數據格式，則抽樣時便會以 id (公司或省份代碼) 為單位，以保持 id 內部的時序特征。

*-Panel Data (sample by cluster(id))
    *-數據預處理
        webuse "nlswork.dta", clear
    xtset id year                  //聲明為面板數據，否則視為截面數據
    gen agesq = age*age
    drop if race==3
    gen black = 2.race
    *-檢驗
    global x "ttl_exp hours tenure south age agesq"
    local m "xtreg ln_wage $x, fe" //模型設定
    bdiff, group(black) model(`m') reps(1000) bs first detail    

• -bdiff- 命令中設定 -bs- 選項，則隨機抽樣為可重復抽樣 (bootstrap)；
• -first- 選項便於將組間系數差異檢驗結果保存在內存中，方便后續使用 esttab 合並到回歸結果表格中。具體使用方法參見 -help bdiff-。
• 由於抽樣過程具有隨機性，因此每次檢驗的結果都有微小差異。在投稿之前，可以附加 -seed()- 選項，以保證檢驗結果的可復制性。
• 其他選項和使用方法參閱 -help bdiff- 的幫助文件。

小結

• 方法1（加入交乘項）在多數模型中都可以使用，但要注意其背后的假設條件是比較嚴格的。若在混合回歸中，只引入你關心的那個變量（ttl_exp）與分組變量 (black) 的交乘項（ttl_exp*black），則相當於假設其他控制變量在兩組之間的不存在系數差異。相對保守的處理方法是：在混合估計時，引入所有變量與分組變量的交乘項，同時附加 vce(robust) 選項，以克服異方差的影響。
• 方法2（基於 SUR 模型的檢驗）執行起來也比較方便，假設條件也比較寬松：允許兩組中所有變量的系數都存在差異，也允許兩組的干擾項具有不同的分布，且彼此相關。局限在於，有些命令無法使用 -suest- 執行聯合估計，如幾乎所有針對面板數據的命令都不支持（-xtreg-， -xtabond- 等）。
• 方法3（組合檢驗）是三種方法中假設條件最為寬松的，只要求原始樣本是從母體中隨機抽取的（看似簡單，但很難檢驗，只能靠嘴說了），而對於兩個樣本組中干擾項的分布，以及衡量組間系數差異的統計量  的分布也未做任何限制。事實上，在獲取經驗 p 值的過程中，我們采用的是“就地取材”、“管中窺豹”的思路，並未假設  的分布函數；另一個好處在於可以適用於各種命令，如 regress，xtreg，logit, ivregress 等，而 -suest- 則只能應用於部分命令。
• 方法無優劣。無論選擇哪種方法，都要預先審視一下是否符合這些檢驗方法的假設條件。