作者：盧梅 | Stata連享會 (知乎 | 簡書 | 碼雲)

Stata連享會現場班課程

1、前言

倍差法（ DID ）是政策評價中被廣泛應用的方法，DID 的結果表明某項政策的實施將會對因變數平均產生多少影響，但是很多時候我們都想知道，在某一個自變數處於某一特定水平，或者是某幾個自變數處於特定水平時，政策效應到底有多大。

比如說，某項政策對於居民健康有正的影響，但是我們想知道這項政策對於 20 歲、 30 歲、 60 歲的居民影響有多大。將邊際分析（margin)應用到 DID 中，可以觀察到在任何一個想要的位置上政策衝擊的所帶來的邊際效應是多少。margins 命令可以實現這一功能。 接下來讓我們近距離觀察一下 margin。

2、margins 命令介紹

在 Stata 命令窗口中輸入 help margins 可以查看完整的幫助文件。事實上，它與 Stata 中的 predict，lincom 等命令有密切的關係。 - predict 的功能是用現有的估計結果（主要是估計所得的係數）來估計觀測對象的統計指標。 - lincom、nlcom的功能在於估計已有係數的線性或非線性表達式。 - margin 命令綜合了 predict 與 lincom、nlcom 的功能，來估計邊際效果。 margin 是指在基於數據集計算的擬合模型基礎上，看在不同的協變數值的情況下，邊際效應是多少。

簡單來說，margin 其實就是，先利用已有的數據估計出模型係數，然後再用這些係數算出線性或非線性表達式的值，接著預測協變數處於一定水平上時的邊際效果。

margin 本身就是一個可以實現眾多功能的命令，關於它此處不做過多介紹，大家有興趣的話可以 help 一下。

3、DID 的邊際分析

Andrew P Wheeler 教授生成一份數據來具體講述 margin 運用到 DID 中的示例，我們也使用 Andrew P Wheeler 生成的數據來給大家說明。

數據結構：Y 代表被解釋變數，Post 和 Exper 可以理解為政策衝擊的虛擬變數和是否為政策實施地區的虛擬變數，它們二者的交乘項係數表示政策效果，還有一個 Ord 序列變數就相當於我們經常用的 id。

很簡單的數據，但是要注意，雖然 DID 及後續分析中用到了 i.month，但是這份數據並不是多期 DID，不同的月份只是代表不同的分類， i.month 表示在全年 12 個月中每個月都有一個特有的 DID 回歸係數，可以將 month 類比為對於 20 多歲組， 30 多歲組， 40 多歲組這種，一月份組，二月份組，十二月份組。

DID 回歸代碼如下，之所以用 xtgee、poisson 廣義估計泊松模型是由這份數據的性質決定的，這個我們可以先不用管它，就把它當作一般的回歸就好。

tsset Exper Ord
asdoc xtgee Y Exper#Post i.Month, family(poisson) corr(ar1) replace

回歸結果為：

GEE population-averaged model Number of obs = 264
Group and time vars: Exper Ord Number of groups = 2
Link: log Obs per group:
Family: Poisson min = 132
Correlation: AR(1) avg = 132.0
max = 132
Wald chi2(14) = 334.52
Scale parameter: 1 Prob > chi2 = 0.0000

------------------------------------------------------------------------------
Y | Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
Exper#Post |
0 1 | .339348 .0898574 3.78 0.000 .1632307 .5154653
1 0 | 1.157615 .0803121 14.41 0.000 1.000206 1.315023
1 1 | .7602288 .0836317 9.09 0.000 .5963138 .9241439
|
Month |
2 | .214908 .1243037 1.73 0.084 -.0287228 .4585388
3 | .2149183 .1264076 1.70 0.089 -.0328361 .4626726
4 | .3119988 .1243131 2.51 0.012 .0683497 .5556479
5 | .4416236 .1210554 3.65 0.000 .2043594 .6788877
6 | .556197 .1184427 4.70 0.000 .3240536 .7883404
7 | .4978252 .1197435 4.16 0.000 .2631322 .7325182
8 | .2581524 .1257715 2.05 0.040 .0116447 .50466
9 | .222813 .1267606 1.76 0.079 -.0256333 .4712592
10 | .430002 .121332 3.54 0.000 .1921956 .6678084
11 | .0923395 .1305857 0.71 0.479 -.1636037 .3482828
12 | -.1479884 .1367331 -1.08 0.279 -.4159803 .1200035
|
_cons | .9699886 .1140566 8.50 0.000 .7464418 1.193536
------------------------------------------------------------------------------

接下來對一月份組進行的邊際分析，就是說此項政策實施對於一月份組來說效果有多大（base 默認為一月），代碼如下：

margins Post#Exper, at( (base) Month )

結果如下：

Adjusted predictions Number of obs = 264
Model VCE : Conventional

Expression : Exponentiated linear prediction, predict()
at : Month = 1

------------------------------------------------------------------------------
| Delta-method
| Margin Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
Post#Exper |
0 0 | 2.637914 .3008716 8.77 0.000 2.048217 3.227612
0 1 | 8.394722 .8243735 10.18 0.000 6.77898 10.01046
1 0 | 3.703716 .3988067 9.29 0.000 2.922069 4.485363
1 1 | 5.641881 .5783881 9.75 0.000 4.508261 6.775501
------------------------------------------------------------------------------

以圖形展示 margin 結果，代碼如下：

marginsplot, recast(scatter)

圖形如下：

如果不想要兩邊邊的 end cross hairs（還是不要讓我翻譯了，哈哈哈，我的第一反應是邊邊起橫長的毛髮）可以輸入如下代碼，然後看一下效果，就知道差別在哪裡了。

marginsplot, recast(scatter) ciopts(recast(rspike))

去掉 end cross hairs 的圖形：

直觀來看，對照組因變數上升了 1 個單位，實驗組因變數下降了 2 個單位。

但是Andrew教授的同事提出，設置對照組的意義在於看到如果沒有實行政策時會是怎樣的情形，對照組因變數是下降的，那麼實驗組因變數也應該 是下降的，我們更應該看的是實驗組因變數的相對變動。

我們需要藉助lincom、nlcom等指令計算已估係數的表達式，然後還原出真正實驗組對於對照組的相對變化。

4、對margin結果的修正

lincom、nlcom之前，我們可以先用test檢驗一下相對下降是否存在，代碼如下：

test 1.Exper#1.Post = (1.Exper#0.Post + 0.Exper#1.Post)

結果為：

chi2( 1) = 46.50
Prob > chi2 = 0.0000

結果顯示，相對下降顯著。 當然也可以對各分組係數做聯合檢驗，代碼如下：

testparm i.Month

結果如下：

(1) 2.Month = 0
(2) 3.Month = 0
(3) 4.Month = 0
(4) 5.Month = 0
(5) 6.Month = 0
(6) 7.Month = 0
(7) 8.Month = 0
(8) 9.Month = 0
(9) 10.Month = 0
(10) 11.Month = 0
(11) 12.Month = 0

chi2( 11) = 61.58
Prob > chi2 = 0.0000

結果依然顯著。接下來我們使用lincom、nlcom來看真實的相對下降是多少。

lincom 1.Exper#1.Post - 0.Exper#1.Post - 1.Exper#0.Post

結果如下：

------------------------------------------------------------------------------
Y | Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
(1) | -.7367337 .1080418 -6.82 0.000 -.9484917 -.5249757
------------------------------------------------------------------------------

這一結果與在回歸模型中，直接加入i.Exper i.Post的效果是一樣的。輸入如下代碼：

xtgee Y i.Exper i.Post Exper#Post i.Month, family(poisson) corr(ar1)

結果如下：

GEE population-averaged model Number of obs = 264
Group and time vars: Exper Ord Number of groups = 2
Link: log Obs per group:
Family: Poisson min = 132
Correlation: AR(1) avg = 132.0
max = 132
Wald chi2(14) = 334.52
Scale parameter: 1 Prob > chi2 = 0.0000

------------------------------------------------------------------------------
Y | Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
1.Exper | 1.157615 .0803121 14.41 0.000 1.000206 1.315023
1.Post | .339348 .0898574 3.78 0.000 .1632307 .5154653
|
Exper#Post |
1 1 | -.7367337 .1080418 -6.82 0.000 -.9484917 -.5249757
|
Month |
2 | .214908 .1243037 1.73 0.084 -.0287228 .4585388
3 | .2149183 .1264076 1.70 0.089 -.0328361 .4626726
4 | .3119988 .1243131 2.51 0.012 .0683497 .5556479
5 | .4416236 .1210554 3.65 0.000 .2043594 .6788877
6 | .556197 .1184427 4.70 0.000 .3240536 .7883404
7 | .4978252 .1197435 4.16 0.000 .2631322 .7325182
8 | .2581524 .1257715 2.05 0.040 .0116447 .50466
9 | .222813 .1267606 1.76 0.079 -.0256333 .4712592
10 | .430002 .121332 3.54 0.000 .1921956 .6678084
11 | .0923395 .1305857 0.71 0.479 -.1636037 .3482828
12 | -.1479884 .1367331 -1.08 0.279 -.4159803 .1200035
|
_cons | .9699886 .1140566 8.50 0.000 .7464418 1.193536
--------------------------------

結果與通過lincom命令得到的交乘項係數是一樣的，因為模型是泊松模型，所以我們再使用nlcom還原到更直觀的數值變動。代碼如下：

nlcom exp(_b[1.Exper] + _b[1.Post] + _b[_cons])

結果如下：

------------------------------------------------------------------------------
Y | Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
_nl_1 | 11.78646 1.583091 7.45 0.000 8.683656 14.88926
------------------------------------------------------------------------------

因為看的是第一個月組的情況，month 前的係數為0，所以nlcom命令常數項後面不用加 month 的係數，如果想看2月份組、3月份組的政策效應的話，就應該把 month 的係數加上。

政策實施後因變數是 11.79 接近 12 的值，再想想之前那個長著 end cross hairs 的圖，實驗組從 8 變到 6 ，其實是不準確的，應該是從 12 變到 6，政策效應是使因變數減少了 8 個單位每月！

所以DID時，其實應該使用xtgee Y i.Exper i.Post Exper#Post i.Month, family(poisson) corr(ar1)這種形式的回歸，得到回歸結果後，再用nlcom計算出實際變動值。

相關鏈接

andrewpwheeler.wordpress.com

主要代碼

tsset Exper Ord
xtgee Y Exper#Post i.Month, family(poisson) corr(ar1)
margins Post#Exper, at( (base) Month )
marginsplot, recast(scatter)
test 1.Exper#1.Post = (1.Exper#0.Post + 0.Exper#1.Post)
testparm i.Month
lincom 1.Exper#1.Post - 0.Exper#1.Post - 1.Exper#0.Post
xtgee Y i.Exper i.Post Exper#Post i.Month, family(poisson) corr(ar1)
nlcom exp(_b[1.Exper] + _b[1.Post] + _b[_cons])

關於我們

  • 【Stata 連享會(公眾號：StataChina)】由中山大學連玉君老師團隊創辦，旨在定期與大家分享 Stata 應用的各種經驗和技巧。
  • 公眾號推文同步發佈於 CSDN-Stata連享會 、簡書-Stata連享會 和 知乎-連玉君Stata專欄。可以在上述網站中搜索關鍵詞Stata或Stata連享會後關注我們。
  • 點擊推文底部【閱讀原文】可以查看推文中的鏈接並下載相關資料。
  • Stata連享會 精彩推文1 || 精彩推文2

聯繫我們

  • 歡迎賜稿： 歡迎將您的文章或筆記投稿至Stata連享會(公眾號: StataChina)，我們會保留您的署名；錄用稿件達五篇以上，即可免費獲得 Stata 現場培訓 (初級或高級選其一) 資格。
  • 意見和資料： 歡迎您的寶貴意見，您也可以來信索取推文中提及的程序和數據。
  • 招募英才： 歡迎加入我們的團隊，一起學習 Stata。合作編輯或撰寫稿件五篇以上，即可免費獲得 Stata 現場培訓 (初級或高級選其一) 資格。
  • 聯繫郵件： [email protected]

往期精彩推文

  • Stata連享會推文列表1
  • Stata連享會推文列表2
  • Stata連享會 精彩推文1 || 精彩推文2


推薦閱讀：
相关文章