全面解读Logit模型

发布时间:2022-05-28 阅读 640

Stata连享会   主页 || 视频 || 推文 || 知乎 || Bilibili 站

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

New! lianxh 命令发布了:
随时搜索推文、Stata 资源。安装:
. ssc install lianxh
详情参见帮助文件 (有惊喜):
. help lianxh
连享会新命令:cnssc, ihelp, rdbalance, gitee, installpkg

课程详情 https://gitee.com/lianxh/Course

课程主页 https://gitee.com/lianxh/Course

⛳ Stata 系列推文:

PDF下载 - 推文合集

作者:陈卓然(中山大学)
邮箱chenzhr25@mail2.sysu.edu.cn

编者按:本文主要摘译自下文,特此致谢!
Source:Uberti L J . Interpreting logit models:[J]. The Stata Journal, 2022, 22(1):60-76. -PDF-


目录


1. 简要回顾

给定一组协变量 X,一个二元被解释变量 Y 的期望值等于以下概率:

在线性概率模型 (LPM) 中,这一条件概率的表达式为 P(Yi=1Xi)=Xiβ,其中 β 可以通过 OLS 来估计。但是线性概率模型的最大问题是估计出来的概率可能会超过 [0,1] 区间。因此对于二元被解释变量,我们往往会采用 Logistic 方程来建模,具体而言:

在 Stata 中,β 的极大似然估计值可以通过 logit 命令估计。为了更好理解,我们使用一份包含 753 个已婚女性的调查数据 (Mroz 1987) 来估计劳动方程。其中,结果变量 inlf 是一个二元变量,取值为 1 意味着参加了工作,取值为 0 意味着未参加工作。协变量 nokid 是一个二元变量,取值为 1 意味着这一妇女在调查时没有小于 6 岁的孩子;nwifeinc 是一个连续变量,衡量了这一妇女的丈夫的年薪 (单位为千欧元)。

. webuse mroz, clear
. gen nokid = (kidslt6 == 0)
. logit inlf nwifeinc i.nokid, robust
. est sto logit 
. reg inlf nwifeinc i.nokid, robust
. est sto OLS
. esttab logit OLS 

--------------------------------------------
                      (1)             (2)   
                     inlf            inlf   
--------------------------------------------
main                                        
nwifeinc          -0.0211**      -0.00479** 
                  (-2.96)         (-3.12)   
0.nokid                 0               0   
                      (.)             (.)   
1.nokid             1.060***        0.255***
                   (5.51)          (5.78)   
_cons              -0.150           0.460***
                  (-0.68)          (8.97)   
--------------------------------------------
N                     753             753   
--------------------------------------------
t statistics in parentheses
* p<0.05, ** p<0.01, *** p<0.001

2. 连续变量的边际影响

连续变量的边际影响可以表示为如下形式:

其中,β 是 xi 的系数。需要注意的是,首先在 LPM 模型中,xi 的边际影响是 Pxi=β,但这 Logit 模型中,xi 的边际影响并不等于 β。其次,Pxi 并不是一个常数,它随着 xi 和 Xi 的变化而变化。因此在报告结果时,我们往往采用如下两种方式:一是取在 xi 和 Xi 均值处的边际影响 (marginal effects at the means),主要通过以下命令实现:

. logit inlf nwifeinc i.nokid, robust
. margins, dydx(nwifeinc) atmeans

Conditional marginal effects                               Number of obs = 753
Model VCE: Robust
Expression: Pr(inlf), predict()
dy/dx wrt:  nwifeinc
At: nwifeinc = 20.12896 (mean)
    0.nokid  = .1952191 (mean)
    1.nokid  = .8047809 (mean)
------------------------------------------------------------------------------
             |            Delta-method
             |      dy/dx   std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
    nwifeinc |     -0.005      0.002    -2.96   0.003       -0.009      -0.002
------------------------------------------------------------------------------

但是上述这种方式,在协变量中存在二元虚拟变量时并不是很符合直觉。因此我们更喜欢第二种方式,即对样本中所有 xi 和 Xi 计算 P/xi,然后对所有的边际影响取平均值。这一方式也是 Stata 中默认计算边际影响的方式。

. margins, dydx(nwifeinc)

Average marginal effects                                   Number of obs = 753
Model VCE: Robust
Expression: Pr(inlf), predict()
dy/dx wrt:  nwifeinc
------------------------------------------------------------------------------
             |            Delta-method
             |      dy/dx   std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
    nwifeinc |     -0.005      0.002    -3.03   0.002       -0.008      -0.002
------------------------------------------------------------------------------

从上述结果可以看出,平均而言,一位丈夫的年薪增加每 1000 欧元,会导致其妻子参加工作的概率降低 0.005。

3. 连续变量的预测概率

在估计完 Logit 模型之后,我们通常会想绘制出预测概率 P^(Yi=1|Xi)。例如,我们让 nwifeinc 从 0 增加到 100000 欧元,同时将妇女区分为有小孩和没有小孩。为了绘制 P^,我们需要使用 marginsplot 等官方命令。但是一个更为推荐的外部命令是 mcp,这一命令可以允许我们在给定协变量的取值范围下绘制 P^(Yi=1|Xi) (在命令选项中指定 at1(numlist)at2(numlist))。

. net install gr0056.pkg, replace 
. mcp nwifeinc nokid, at1(0(10)100) at2(0 1) ci plotopts(ycommon)

此外,我们也可以使用外部命令 mtable 将上图以表格的形式列示出来,具体如下:

. net install spost13_ado.pkg, replace 
. mtable, at(nokid = (0) nwifeinc = (0(20)100)) at(nokid = (1) nwifeinc = (0(20)100)) brief

Expression: Pr(inlf), predict()
           | nwifeinc     nokid     Pr(y)
 ----------+-----------------------------
         1 |        0         0     0.463
         2 |       20         0     0.361
         3 |       40         0     0.270
         4 |       60         0     0.195
         5 |       80         0     0.137
         6 |      100         0     0.095
         7 |        0         1     0.713
         8 |       20         1     0.620
         9 |       40         1     0.517
        10 |       60         1     0.412
        11 |       80         1     0.315
        12 |      100         1     0.232

4. Logit 与 LPM 的对比

. logit inlf nwifeinc i.nokid, robust
. eststo margin: margins, dydx(nwifeinc) post
. est sto logitmargins
. reg inlf nwifeinc i.nokid, r
. eststo margin: margins, dydx(nwifeinc) post
. est sto OLSmargins
. esttab logitmargins OLSmargins, mtitles("Logit" "OLS")

--------------------------------------------
                      (1)             (2)   
                    Logit             OLS   
--------------------------------------------
nwifeinc         -0.00489**      -0.00479** 
                  (-3.03)         (-3.12)   
--------------------------------------------
N                     753             753   
--------------------------------------------
t statistics in parentheses
* p<0.05, ** p<0.01, *** p<0.001

从图中不难看出,在 LPM 模型下丈夫工资的边际影响是 -0.00479,这与 Logit 模型的边际影响非常相似。这一点也可以从上图中看出,即丈夫工资的边际影响基本是线性的。此外如果我们对比二者的预测概率,不难发现结果也是极为相似的。

. logit inlf nwifeinc i.nokid,robust
. predict P_logit, pr
. regress inlf nwifeinc i.nokid, robust
. predict P_LPM, xb
. summarize P_LPM P_logit

    Variable |        Obs        Mean    Std. dev.       Min        Max
-------------+---------------------------------------------------------
       P_LPM |        753    .5683931    .1166285  -.0003462   .7147611
     P_logit |        753    .5683931     .116784   .1020306   .7132691

. correlate P_LPM P_logit

             |    P_LPM  P_logit
-------------+------------------
       P_LPM |   1.0000
     P_logit |   0.9987   1.0000

但是如果数据生成过程是高度非线性的,那么 Logit 和 LPM 模型的差异可能就会变得非常大。

5. 连续变量的偏效应

假设我们想要知道丈夫工资增加一个标准差 (δ=11635 欧元),那么妇女参加工作的概率会降低多少。一个很常见的做法是将边际效应 (0.005) 乘以 11.63 (千欧元),从而得到 P 将下降 5.7 个百分点。然而上述计算方法的前提是边际影响是线性的。实际上,正确的计算方法如下:

在 Stata 中,我们使用 mchange 命令实现上述计算,具体代码如下:

. logit inlf nwifeinc i.nokid, robust
. mchange nwifeinc, amount(sd) stats(change se zvalue pvalue ll ul) decimals(4) brief  
. // ll 表示置信区间的下界,ul 表示置信区间的上界

logit: Changes in Pr(y) | Number of obs = 753
Expression: Pr(inlf), predict(pr)
             |    Change    Std Err    z-value    p-value         LL         UL 
-------------+------------------------------------------------------------------
nwifeinc     |                                                                  
         +SD |   -0.0576     0.0192    -3.0074     0.0026    -0.0951    -0.0201 

如果不想改变一个标准差,则可以将上述命令中的 amount(sd) 换成 delta(#)。比如说 delta(4) 表示增加4。

6. 二元变量的边际影响

当协变量 xi 是一个二元变量时,我们可以通过两种方式去评估其边际影响的大小。首先,我们可以计算结果变量为 1 的预测概率。例如,有小孩妇女工作的概率 P^(Yi=1xi=1,Xi) 和没有小孩妇女工作的概率 P^(Yi=1xi=0,Xi);其次,我们可以计算这两种预测概率的差值。在 Stata 中,我们可以通过如下两种方式实现:

. * 方法一
. logit inlf nwifeinc i.nokid, robust
. margins i.nokid

Predictive margins                                         Number of obs = 753
Model VCE: Robust
Expression: Pr(inlf), predict()
------------------------------------------------------------------------------
             |            Delta-method
             |     Margin   std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
       nokid |
          0  |      0.362      0.039     9.18   0.000        0.285       0.440
          1  |      0.618      0.020    31.52   0.000        0.580       0.657
------------------------------------------------------------------------------
 
. * 方法二
. logit inlf nwifeinc i.nokid, robust
. margins r.nokid

Contrasts of predictive margins                            Number of obs = 753
Model VCE: Robust
Expression: Pr(inlf), predict()
------------------------------------------------
             |         df        chi2     P>chi2
-------------+----------------------------------
       nokid |          1       33.64     0.0000
------------------------------------------------
--------------------------------------------------------------
             |            Delta-method
             |   Contrast   std. err.     [95% conf. interval]
-------------+------------------------------------------------
       nokid |
   (1 vs 0)  |      0.256      0.044         0.169       0.342
--------------------------------------------------------------

可以看出,平均而言,在没有小孩的妇女中有 62% 的妇女参加了工作,而在拥有小孩的妇女中有 36% 的妇女参加了工作。相应地,没有小孩可以使一个妇女参加工作的概率提升大约 26 个百分点,这种影响是显著异于 0。当然你也可以使用和连续变量相同的命令来得到同样的结果,此时 Stata 会非常贴心地提醒你这一个离散变量,因此报告的是两种预测概率的差值。

. margins, dydx(nokid)

Average marginal effects                                   Number of obs = 753
Model VCE: Robust
Expression: Pr(inlf), predict()
dy/dx wrt:  1.nokid
------------------------------------------------------------------------------
             |            Delta-method
             |      dy/dx   std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
     1.nokid |      0.256      0.044     5.80   0.000        0.169       0.342
------------------------------------------------------------------------------
Note: dy/dx for factor levels is the discrete change from the base level.

7. 弹性

有时我们需要报告预测概率的相对变化而非绝对变化,此时就需要计算半弹性。如果 xi 是连续变量,那么半弹性的计算公式为:

在 Stata 中可以通过如下命令实现:

. logit inlf nwifeinc i.nokid,robust
. margins, eydx(nwifeinc)

Average marginal effects                                   Number of obs = 753
Model VCE: Robust
Expression: Pr(inlf), predict()
ey/dx wrt:  nwifeinc
------------------------------------------------------------------------------
             |            Delta-method
             |      ey/dx   std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
    nwifeinc |     -0.009      0.003    -2.93   0.003       -0.015      -0.003
------------------------------------------------------------------------------

可以看出,平均而言,丈夫的工资增加 1000 欧元可以导致妇女参加工作的概率降低 0.9%。

当 xi 为二元变量时,半弹性的计算公式为:

当 ΔP 很小时,采用对数近似的方法,我们可以将上式改写为:

不幸的是在 Stata 中,我们使用 logit 命令估计完模型之后,再使用 margins, eydx(i.nokid) 命令,计算的是改写式。

8. odds ratio 和 risk ratio

odds 定义为一件事情发生概率与这件事情没有发生概率的比值,即 P(Yi=1xi,Xi)/P(Yi=0xi,Xi),其中 xi 是一个二元变量。具体地,当 xi 从 0 变动到 1 时,

其中 β 是 xi 的系数。在 Stata 中,对于 odds ratio 的计算可以采用如下方式:

. logit inlf nwifeinc i.nokid,robust
. nlcom exp(_b[1.nokid])

       _nl_1: exp(_b[1.nokid])
------------------------------------------------------------------------------
        inlf | Coefficient  Std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
       _nl_1 |      2.887      0.556     5.20   0.000        1.798       3.977
------------------------------------------------------------------------------

可以看出,有小孩的妇女参加工作的概率是没有小孩的妇女参加工作概率的 3 倍左右。当然你也可以选择在 logit 命令中加入 or 的选项,Stata 会自动报告全部协变量的 odds ratio,不过当协变量连续时,odds ratio 不太好解释。我们通常也会报告风险比 risk ratio,具体计算如下:

不幸的是在 Stata 中并没有直接计算 risk ratio 的命令,我们需要借助 margins 命令的 expression 选项来计算。

. logit inlf nwifeinc i.nokid, robust or

Logistic regression                                     Number of obs =    753
                                                        Wald chi2(2)  =  37.78
                                                        Prob > chi2   = 0.0000
Log pseudolikelihood = -493.81541                       Pseudo R2     = 0.0409
------------------------------------------------------------------------------
             |               Robust
        inlf | Odds ratio   std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
    nwifeinc |      0.979      0.007    -2.96   0.003        0.966       0.993
     1.nokid |      2.887      0.556     5.51   0.000        1.980       4.210
       _cons |      0.861      0.189    -0.68   0.495        0.560       1.324
------------------------------------------------------------------------------
Note: _cons estimates baseline odds.

. margins, expression((1+exp(-nwifeinc*_b[nwifeinc]-0*_b[1.nokid]-_b[_cons])) ///
>     / (1+exp(-nwifeinc*_b[nwifeinc]-1*_b[1.nokid]-_b[_cons])))


Predictive margins                                         Number of obs = 753
Model VCE: Robust
Expression: (1+exp(-nwifeinc*_b[nwifeinc]-0*_b[1.nokid]-_b[_cons])) 
			/ (1+exp(-nwifeinc*_b[nwifeinc]-1*_b[1.nokid]-_b[_cons]))
------------------------------------------------------------------------------
             |            Delta-method
             |     Margin   std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
       _cons |      1.721      0.199     8.65   0.000        1.331       2.110
------------------------------------------------------------------------------

可以看出,没有孩子的妇女参加工作的概率是有孩子的妇女参加工作概率的 1.72 倍。

9. 交互项

通常来讲,一个接受过良好高等教育的女性往往更为独立,丈夫的工资对其影响也相对较小。因此我们可以认为,丈夫的工资对妇女参加工作概率的影响,受到妇女教育水平的调节。为此,我们在模型中加入 nwifeinceduc 的交乘项。

在 Stata 中,我们可以通过如下命令语句实现上述计算:

. logit inlf c.nwifeinc i.nokid c.nwifeinc#c.educ c.educ, robust nolog

Logistic regression                                     Number of obs =    753
                                                        Wald chi2(4)  =  75.22
                                                        Prob > chi2   = 0.0000
Log pseudolikelihood = -468.34427                       Pseudo R2     = 0.0904
-----------------------------------------------------------------------------------
                  |               Robust
             inlf | Coefficient  std. err.      z    P>|z|     [95% conf. interval]
------------------+----------------------------------------------------------------
         nwifeinc |     -0.048      0.050    -0.95   0.343       -0.146       0.051
          1.nokid |      1.259      0.202     6.24   0.000        0.864       1.655
c.nwifeinc#c.educ |      0.001      0.004     0.20   0.841       -0.007       0.008
             educ |      0.249      0.080     3.12   0.002        0.093       0.406
            _cons |     -3.015      1.021    -2.95   0.003       -5.017      -1.014
-----------------------------------------------------------------------------------

. testparm nwifeinc c.nwifeinc#c.educ

 ( 1)  [inlf]nwifeinc = 0
 ( 2)  [inlf]c.nwifeinc#c.educ = 0
           chi2(  2) =   20.11
         Prob > chi2 =    0.0000

不难看出,nwifeincc.nwifeinc#c.educ 是联合显著的,并且 c.nwifeinc#c.educ 的系数似乎告诉我们,多受一年教育可以削弱丈夫的工资对妇女参加工作积极性的负面影响。然而这样的想法在非线性模型中是错误的。

原因很简单,根据定义交互效应衡量了教育的边际改变是如何修正处理效应的,即

在 LPM 模型当中,这个二阶偏导等价于 P/(nwifeinc×educ),也就是这个交乘项的系数,但是在非线性模型中这个二阶偏导并非交互项的系数。具体在 Stata 当中,交互效应可以通过在使用 logit 后,采用 inteff 命令进行正确的计算。

. generate nwifeinc_educ = nwifeinc*educ
. quietly logit inlf nwifeinc educ nwifeinc_educ nokid, robust nolog
. net install st0063.pkg, replace 
. inteff inlf nwifeinc educ nwifeinc_educ nokid

    Variable |        Obs        Mean    Std. dev.       Min        Max
-------------+---------------------------------------------------------
   _logit_ie |        753    .0004244      .00062  -.0009878   .0010581
   _logit_se |        753    .0008147    .0002286   .0000674   .0016656
    _logit_z |        753     .547078    1.268062  -7.134073   5.891206

需要注意的是,交互效应 (educ)=2P/(nwifeinc)(educ) 的估计是随着解释变量取值的不同而不同的。从上表也可以看出,交互项的置信区间是非常大的 (从 -0.00099 到 0.00106)。另外一点值得注意的是 _logit_se_logit_z 分别报告了交互效应的标准误和 z 统计量,他们的 mean 并不是上表第二行中交互效应平均值的标准误和 z 统计量,而是平均的交互效应标准误和 z 统计量。

不过我们并不关心交互效应,而是更关心下式:

. quietly logit inlf c.nwifeinc c.educ c.nwifeinc#c.educ i.nokid, robust nolog
. quietly margins, dydx(nwifeinc) at(educ=(5(1)17) nokid=(0) nwifeinc=(60))
. quietly marginsplot, recast(line) recastci(rarea) name(kidhigh)  ///
>     scheme(white_tableau) title("nokid=0, nwifeinc=60")
. quietly margins, dydx(nwifeinc) at(educ=(5(1)17) nokid=(1) nwifeinc=(10))
. quietly marginsplot, recast(line) recastci(rarea) name(nokidlow) ///
>     scheme(white_tableau) title("nokid=1, nwifeinc=10")
. graph combine kidhigh nokidlow, ycommon xcommon

从上图可以看出,丈夫的工资对妇女参加工作概率的影响取决于妇女的教育水平的。对于没有孩子且丈夫的工资又很低的妇女,其受教育年限的提高可以增加参加工作的概率,这一点符合我们的常识。然而对于那些有孩子且丈夫的工资又很高的妇女,我们不禁吃惊地发现随着妇女教育水平的提高,其参加工作的积极性反而会下降。

如果想要计算某一特定点上的交互效应,我们需要采用如下的方式进行计算:

. quietly logit inlf c.nwifeinc c.educ c.nwifeinc#c.educ i.nokid, robust nolog
. global L "(1 / (1+exp(-_b[c.nwifeinc]*c.nwifeinc - _b[c.educ]*c.educ -_b[c.nwifeinc#c.educ]*
>	  c.nwifeinc#c.educ - _b[1.nokid]*nokid - _b[_cons])))"
. margins, expression( _b[c.nwifeinc#c.educ] * $L * (1 - $L ) + (_b[c.nwifeinc]         ///
>     + _b[c.nwifeinc#c.educ]*c.educ)*(_b[c.educ] + _b[c.nwifeinc#c.educ]*c.nwifeinc) * ///
>     $L * (1 - $L ) * (1 - 2*$L ) ) at(nokid=0 nwifeinc=60 educ=(5(3)17))


1._at: nwifeinc = 60
       educ     =  5
       nokid    =  0
2._at: nwifeinc = 60
       educ     =  8
       nokid    =  0
3._at: nwifeinc = 60
       educ     = 11
       nokid    =  0
4._at: nwifeinc = 60
       educ     = 14
       nokid    =  0
5._at: nwifeinc = 60
       educ     = 17
       nokid    =  0
------------------------------------------------------------------------------
             |            Delta-method
             |     Margin   std. err.      z    P>|z|     [95% conf. interval]
-------------+----------------------------------------------------------------
         _at |
          1  |     -0.000      0.000    -2.70   0.007       -0.000      -0.000
          2  |     -0.000      0.000    -5.34   0.000       -0.000      -0.000
          3  |     -0.001      0.000    -4.73   0.000       -0.001      -0.000
          4  |     -0.001      0.000    -3.84   0.000       -0.001      -0.000
          5  |     -0.001      0.001    -0.61   0.545       -0.003       0.002
------------------------------------------------------------------------------

10. 相关推文

Note:产生如下推文列表的 Stata 命令为:
lianxh probit logit, m
安装最新版 lianxh 命令:
ssc install lianxh, replace

相关课程

免费公开课

最新课程-直播课

专题 嘉宾 直播/回看视频
最新专题 文本分析、机器学习、效率专题、生存分析等
研究设计 连玉君 我的特斯拉-实证研究设计-幻灯片-
面板模型 连玉君 动态面板模型-幻灯片-
面板模型 连玉君 直击面板数据模型 [免费公开课,2小时]
  • Note: 部分课程的资料,PPT 等可以前往 连享会-直播课 主页查看,下载。

课程主页

课程主页

关于我们

  • Stata连享会 由中山大学连玉君老师团队创办,定期分享实证分析经验。
  • 连享会-主页知乎专栏,700+ 推文,实证分析不再抓狂。直播间 有很多视频课程,可以随时观看。
  • 公众号关键词搜索/回复 功能已经上线。大家可以在公众号左下角点击键盘图标,输入简要关键词,以便快速呈现历史推文,获取工具软件和数据下载。常见关键词:课程, 直播, 视频, 客服, 模型设定, 研究设计, stata, plus, 绘图, 编程, 面板, 论文重现, 可视化, RDD, DID, PSM, 合成控制法

连享会小程序:扫一扫,看推文,看视频……

扫码加入连享会微信群,提问交流更方便

✏ 连享会-常见问题解答:
https://gitee.com/lianxh/Course/wikis

New! lianxhsongbl 命令发布了:
随时搜索连享会推文、Stata 资源,安装命令如下:
. ssc install lianxh
使用详情参见帮助文件 (有惊喜):
. help lianxh