二元选择模型：probit 还是 logit？

作者：黄彩虹 [编译] (知乎 | 简书 | 码云)  
     
本期责任编辑：王俊

---

Source:  

• Probit or Logit: ladies and gentlemen, pick your weapon
  

• regress, Probit, or Logit?

• Stata现场班 (2019.1.17-27) 课程详情

• Stata寒假精品现场班（2019.1月17日-27日）报名主页

开始接受报名

1 解读 Probit 模型和 Logit 模型的边际效应估计

对于二元被解释变量的分析常采用 Probit 或 Logit 模型。就模型设定而言，Logit 模型更简单。不过，借助 Stata 中的 margins  命令，上述两个模型都易于实现。两种模型估计将返回同样的结果，可以根据习惯或偏好来选用。

下面通过研究人员感兴趣的一系列效应来展示 Probit 和 Logit 模型的相同拟合结果。主要考虑连续协变量与离散协变量的变化对正向估计结果的影响，计算作用效应的平均值及协变量的均值。换句话说，我研究了平均边际效应（AME），平均处理效应（ATE），协变量均值处的边际效应（MEM）、协变量均值处的处理效应（TEM）。

1.1 估计结果

表1给出了真实数据生成过程符合 Probit 模型假定的估计结果，并给出 AME 和 ATE 估计的均值以及在原假设 5% 拒绝概率下的 Probit 和 Logit 估计。同时使用 2 千万观测样本，计算了 AME 和 ATE 的近似真实值，以及系数的真实值。

表1： AME和ATE ：真实 DGP Probit    

N=10000, k=4000 次重复的拟合结果

DGP: Data Generation Process (数据生成过程)

统计量	Approximate True Value	Probit	Logit
x1的AME	-.1536	-.1537	-.1537
5%拒绝概率		.050	.052
x2的ATE	.1418	.1417	.1417
5%拒绝概率		.050	.049

表2： MEM 和 TEM ：真实 DGP Probit 

N=10000, k=4000 次重复的拟合结果

统计量	Approximate True Value	Probit	Logit
x1的MEM	-.1672	-.1673	-.1537
5%拒绝概率		.056	.06
x2的TEM	.1499	.1498	.1471
5%拒绝概率		.053	.058

Logit 估计值接近真实值，拒绝率接近 5% 。当真实 DGP 满足 Probit 模型的假设时，用 Logit 对模型参数进行拟合不会带来偏误。如果真实 DGP 满足 Logit 模型的假设，结论是一致的，表 3 -表 4 给出了估计结果。

表3：平均边际和处理效应：真实 DGP Logit

N=10000, k=4000 次重复的拟合结果

统计量	Approximate True Value	Probit	Logit
x1的AME	-.1090	-.1088	-.1089
5%拒绝概率		.052	.052
x2的ATE	.1046	.1044	.1045
5%拒绝概率		.053	.051

表4： MEM 和 TEM ：真实 DGP Logit

N=10000, k=4000次重复的拟合结果

统计量	Approximate True Value	Probit	Logit
x1的MEM	-.1146	-.1138	-.1146
5%拒绝概率		.050	.051
x2的TEM	.1086	.1081	.1085
5%拒绝概率		.058	.058

1.2 原因是什么？

从上面的蒙特卡洛模拟分析结果可以看出，Logit 和 Probit 估计得到的结果非常接近。

极大似然估计为我们寻找使得数据符合分布假设的参数。所选择的似然是对真实似然的近似，如果二者接近的话，这就是一种有效的近似。将基于似然的模型视为有用的近似，而不是作为真实似然的模型，是拟似然理论的基础。详见 White(1996) 和 Wooldridge(2010) 。假设 Probit 模型和 Logit 模型中的不可观测随机变量分别来自标准正态分布和 Logistic 分布。这两种情况下的累积分布函数 (CDFs) 相互接近，特别是在平均值附近。因此，这两组假设下的估计量产生了相似的结果。为了说明这些论点，我们可以画出这两支 CDF 及其不同之处：

从左侧或右侧靠近均值时，两种分布函数的差异接近于 0 且始终小于 0.15 。

1.3 模特卡洛模拟分析过程解析

下面是进行模拟的代码。

• 第一部分，即 4-12 行，生成满足 Probit 模型假设的结果变量 y1 ,及满足 Logit 模型假设的结果变量 y2 。

• 第二部分，即 13-16 行，计算 Logit 模型和 Probit 模型的边际效应。对于离散协变量，边际效应是一种处理效应。

• 第三部分，即 17-25 行，计算均值处估计的边际效应。后面将用这些估计值来对这些影响的真实值进行近似估计。

1. program define mkdata

2.        syntax, [n(integer 1000)]

3.        clear

4.        quietly set obs `n'

5.        // 1. Generating data from Probit, Logit, and misspecified

6.        generate x1    = rnormal()

7.        generate x2    = rbeta(2,4)>.5

8.        generate e1    = rnormal()

9.        generate u     = runiform()

10.        generate e2    = ln(u) -ln(1-u)

11.        generate xb    = .5*(1 -x1 + x2)

12.        generate y1    =  xb + e1 > 0

13.        generate y2    =  xb + e2 > 0

14.        // 2. Computing Probit & Logit marginal and treatment effects

15.        generate m1    = normalden(xb)*(-.5)

16.        generate m2    = normal(1 -.5*x1 ) - normal(.5 -.5*x1)

17.        generate m1l   = exp(xb)*(-.5)/(1+exp(xb))^2

18.        generate m2l   = exp(1 -.5*x1)/(1+ exp(1 -.5*x1 )) - ///

19.                         exp(.5 -.5*x1)/(1+ exp(.5 -.5*x1 ))

20.        // 3. Computing Probit & Logit marginal and treatment effects at means

21.        quietly mean x1 x2

22.        matrix A         = r(table)

23.        scalar a         = .5 -.5*A[1,1] + .5*A[1,2]

24.        scalar b1        =  1 -.5*A[1,1]

25.        scalar b0        = .5 -.5*A[1,1]

26.        generate mean1   = normalden(a)*(-.5)

27.        generate mean2   = normal(b1) - normal(b0)

28.        generate mean1l  = exp(a)*(-.5)/(1+exp(a))^2

29.        generate mean2l  = exp(b1)/(1+ exp(b1)) - exp(b0)/(1+ exp(b0))

30. end

我用 2000 万样本观测数据近似估计了真实边际效应，在这个例子中是一种合理策略。例如，对 Probit 模型中连续的协变量 xk 取平均边际效应：

上式是对的近似估计，为了获得这个期望值，需要对所有协变量的分布进行整合。这并不实际，且会限制协变量的选择。因此，采用2000万样本观测数据的来代替，把它视作真实值。对于其他边际效应，遵循相同的逻辑。下面是计算近似真实边际效应的代码。调用2000万的观测结果，计算了在模拟中将使用的平均值，并为每一个近似的真实值创建了暂元。

1. . mkdata, n(20000000)

2. . local values "m1 m2 m1l m2l mean1 mean2 mean1l mean2l"

3. . local means  "mx1 mx2 mx1l mx2l meanx1 meanx2 meanx1l meanx2l"

4. . local n : word count `values'

5. . forvalues i= 1/`n' {

6.  2.         local a: word `i' of `values'

7.  3.         local b: word `i' of `means'

8.  4.         sum `a', meanonly

9.  5.         local `b' = r(mean)

10.  6. }

现在可以进行以上部分用于生成结果的所有估计了。当真实 DGP 符合 Probit 模型时估计 ATE 和 AM 的命令是：

1. . postfile mProbit y1p y1p_r y1l y1l_r y2p y2p_r y2l y2l_r ///

2.                   using simsmProbit, replace

3. . forvalues i=1/4000 {

4.  2.  quietly {

5.  3.     mkdata, n(10000)

6.  4.     Probit y1 x1 i.x2, vce(robust)

7.  5.     margins, dydx(*) atmeans post

8.  6.     local y1p = _b[x1]

9.  7.     test _b[x1] = `meanx1'

10.  8.     local y1p_r   = (r(p)<.05)

11.  9.     local y2p = _b[1.x2]

12. 10.     test _b[1.x2] = `meanx2'

13. 11.     local y2p_r   = (r(p)<.05)

14. 12.     Logit  y1 x1 i.x2, vce(robust)

15. 13.     margins, dydx(*) atmeans post

16. 14.     local y1l = _b[x1]

17. 15.     test _b[x1] = `meanx1'

18. 16.     local y1l_r   = (r(p)<.05)

19. 17.     local y2l = _b[1.x2]

20. 18.     test _b[1.x2] = `meanx2'

21. 19.     local y2l_r   = (r(p)<.05)

22. 20.     post mProbit (`y1p') (`y1p_r') (`y1l') (`y1l_r') ///

23.                 (`y2p') (`y2p_r') (`y2l') (`y2l_r')

24. 21.  }

25. 22. }

1. . use simsProbit

2. . summarize

3.    Variable |        Obs        Mean    Std. Dev.       Min        Max

4. -------------+---------------------------------------------------------

5.         y1p |      4,000   -.1536812    .0038952  -.1697037  -.1396532

6.       y1p_r |      4,000         .05    .2179722          0          1

7.         y1l |      4,000   -.1536778    .0039179  -.1692524  -.1396366

8.       y1l_r |      4,000      .05175    .2215496          0          1

9.         y2p |      4,000     .141708    .0097155   .1111133   .1800973

10. -------------+---------------------------------------------------------

11.       y2p_r |      4,000       .0495    .2169367          0          1

12.         y2l |      4,000    .1416983    .0097459   .1102069   .1789895

13.       y2l_r |      4,000        .049     .215895          0          1

对于真实 DGP 符合 Probit 模型的 MEM 和 TEM ，使用包含 atmeans 选项的 margins 命令。使用稳健标准误以确保最逼近真实似然，并且使用选项vce以说明使用了两步 M 估计。两步 M 估计详见 Wooldridge (2010) 。

1.4 总结

给出了 Probit 与 Logit 的估计结果差异可以忽略不计的模拟运算证据，原因在于拟似然理论，特别是， Probit 与 Logit 模型的累积分布函数是相似的，尤其是在均值附近。
相关code

2. Probit/Logit 与线性概率模型对比

也有很多人直接使用 regress 估计二值选择模型，这其实是在估计 线性概率模型，其实就是一个简单的 OLS 回归。由于模型的被解释变量是概率值 (取值介于 0 和 1 之间)，而右侧的线性拟合值的取值范围未必在 0 和 1 之间，这个模型的局限也就显而易见了。

当真实模型是 Probit 或 Logit 时，运算结果显示，对于研究者所关注的边际效应，使用线性概率模型将生成不一致的估计。结果取决于真实模型究竟是 Probit 模型还是 Logit 模型。

2.1 估计结果

以下模拟使用 1000 个观察值，执行 5000 次重复（replications）。真实的数据生成过程 （DGPs） 是用一个离散的协变量和一个连续的协变量来构造的。

我研究了连续变量变化对条件概率 (AME) 的平均影响和离散协变量变化对条件概率 (ATE) 的平均影响。根据协变量的平均值计算，我还研究了连续变量的变化对条件概率的影响 (MEM) ，以及离散协变量的变化对条件概率的影响 (TEM) 。

表 1 给出了真实 DGP 满足 Logit 模型假设时的模拟结果。给出了 AME 和 ATE 估计的平均值，以及真实的零假设的 5% 的拒绝率，还提供了 AME 和 ATE 的近似真值。在系数的真实值处利用 2000 万个观测数据，通过计算 ATE 和 AME ，得到了近似的真实值。

表1： AME和ATE ：真实 DGP Probit

N=10000, k=5000 次重复的拟合结果

统计量	Approximate True Value	Logit	Regress (LPM)
x1的AME	-.084	-.084	-0.1537
5%拒绝概率		.050	.99
x2的ATE	.092	.091	.091
5%拒绝概率		.058	.058

从表1中可以看出， Logit 模型估计接近真实值，真零假设的拒绝率接近 5% 。对于线性概率模型， AME 的拒绝率为 99% 。对于 ATE ，拒绝率和点估计接近使用 Logit 估计的值。对于 MEM 和 TEM ，有:

表2： MEM 和 TEM ：真实 DGP Probit

N=10000, k=5000次重复的拟合结果

统计量	Approximate True Value	Logit	Regress (LPM)
x1的MEM	-.099	-.099	-.094
5%拒绝概率		.054	.618
x2的TEM	.109	.109	.092
5%拒绝概率		.062	.073

同样， Logit 估计与预期一致。对于线性概率模型，真实零假设的拒绝率为 62% 。对 TEM 的拒绝率为 7.3% ，估计效应小于真实效应。对于 AME 和 ATE ，当真实 DGP 为 Probit 模型时，结果如下：

表3：平均边际和处理效应：真实 DGP Probit

N=10000, 5000次重复的拟合结果

统计量	Approximate True Value	Probit	Regress (LPM)
x1的AME	-.094	-.094	-.121
5%拒绝概率		.047	1
x2的ATE	.111	.111	.111
5%拒绝概率		.065	.061

Probit 模型估计接近真实值，真实原假设的拒绝率接近 5% 。对于线性概率模型， AME 的拒绝率为 100% 。对于 ATE ，拒绝率和点估计与 Probit 估计的值近似。对于 MEM 和 TEM，结果如下：

表4： MEM 和 TEM ：真实 DGP Probit

N=10000, 5000次重复的拟合结果

统计量	Approximate True Value	Probit	Regress (LPM)
x1的MEM	-.121	-.122	-.121
5%拒绝概率		.063	.054
x2的TEM	.150	.150	.110
5%拒绝概率		.059	.158

对于 MEM ， Probit 模型和线性概率模型给出了可靠的推论。对于 TEM ， Probit 边际效应表现与预期一致，但线性概率模型的拒绝率为 16% ，且点估计与真实值不接近。

2.2 模特卡洛模拟分析过程解析

• 在第一部分，第 6 至 13 行，生成满足 Logit 模型和 Probit 模型假设的结果变量， y 和 yp 。

• 在第二部分，第 15 至 19 行，计算了 Logit 和 proit 模型的边际效应。对于离散协变量，边际效应是一种处理效应。

• 在第三部分，第 21 至 29 行，计算了均值处的边际效应。稍后将使用这些估计来计算真实效应的近似值。

1. program define mkdata

2.    syntax, [n(integer 1000)]

3.    clear

4.    quietly set obs `n'

5.    // 1. Generating data from Probit, Logit, and misspecified

6.    generate x1    = rchi2(2)-2

7.    generate x2    = rbeta(4,2)>.2

8.    generate u     = runiform()

9.    generate e     = ln(u) -ln(1-u)

10.    generate ep    = rnormal()

11.    generate xb    = .5*(1 - x1 + x2)

12.    generate y     =  xb + e > 0

13.    generate yp    = xb + ep > 0

14.    // 2. Computing  Probit & Logit marginal and treatment effects

15.    generate m1   = exp(xb)*(-.5)/(1+exp(xb))^2

16.    generate m2   = exp(1 -.5*x1)/(1+ exp(1 -.5*x1 )) - ///

17.                  exp(.5 -.5*x1)/(1+ exp(.5 -.5*x1 ))

18.    generate m1p  = normalden(xb)*(-.5)

19.    generate m2p  = normal(1 -.5*x1 ) - normal(.5 -.5*x1)

20.    // 3. Computing marginal and treatment effects at means

21.    quietly mean x1 x2

22.    matrix A        = r(table)

23.    scalar a        = .5 -.5*A[1,1] + .5*A[1,2]

24.    scalar b1       =  1 -.5*A[1,1]

25.    scalar b0       = .5 -.5*A[1,1]

26.    generate mean1  = exp(a)*(-.5)/(1+exp(a))^2

27.    generate mean2  = exp(b1)/(1+ exp(b1)) - exp(b0)/(1+ exp(b0))

28.    generate mean1p = normalden(a)*(-.5)

29.    generate mean2p = normal(b1) - normal(b0)

30. end

用 2000 万个观测数据来近似真实的边际效应。在这种情况下，这是一个合理的策略。例如，以连续协变量 xk 的平均边际效应为例，在 Probit 模型中：

上式是对的近似估计，为了获得这个期望值，需要对所有协变量的分布进行整合。这是不实际的，将限制协变量的选择。相反，我抽取了 2000 万个观测数据，计算，将其作为真实值。其他边际效应的估计遵循同样的逻辑。下面是计算近似真实边际效应的代码。采用 2000 万个观测数据，计算了模拟中使用的平均值，并为每一个近似的真实值创建了暂元。

1. . mkdata, n(`L')

2. (2 missing values generated)

3. . local values "m1 m2 mean1 mean2 m1p m2p mean1p mean2p"

4. . local means  "mx1 mx2 meanx1 meanx2 mx1p mx2p meanx1p meanx2p"

5. . local n : word count `values'

6. .

7. . forvalues i= 1/`n' {

8.  2.         local a: word `i' of `values'

9.  3.         local b: word `i' of `means'

10.  4.         sum `a', meanonly

11.  5.         local `b' = r(mean)

12.  6. }  

现在可以进行以上部分用于生成结果的所有估计了。当真实 DGP 符合 Logit 模型时估计 TEM 和 MEM 的命令是：

1. . postfile lpm y1l y1l_r y1lp y1lp_r y2l y2l_r y2lp y2lp_r ///

2. >                 using simslpm, replace

3. . forvalues i=1/`R' {

4.  2.         quietly {

5.  3.                 mkdata, n(`N')

6.  4.                 Logit  y x1 i.x2, vce(robust)

7.  5.                 margins, dydx(*) atmeans post  vce(unconditional)

8.  6.                 local y1l = _b[x1]

9.  7.                 test _b[x1] = `meanx1'

10.  8.                 local y1l_r   = (r(p)<.05)

11.  9.                 local y2l = _b[1.x2]

12. 10.                 test _b[1.x2] = `meanx2'

13. 11.                 local y2l_r   = (r(p)<.05)

14. 12.                 regress  y x1 i.x2, vce(robust)

15. 13.                 margins, dydx(*) atmeans post  vce(unconditional)

16. 14.                 local y1lp = _b[x1]

17. 15.                 test _b[x1] = `meanx1'

18. 16.                 local y1lp_r   = (r(p)<.05)

19. 17.                 local y2lp = _b[1.x2]

20. 18.                 test _b[1.x2] = `meanx2'

21. 19.                 local y2lp_r   = (r(p)<.05)

22. 20.                 post lpm (`y1l') (`y1l_r') (`y1lp') (`y1lp_r') ///

23. >                          (`y2l') (`y2l_r') (`y2lp') (`y2lp_r')

24. 21.         }

25. 22. }

26. . postclose lpm

27. . use simslpm, clear

28. . sum

29.    Variable |        Obs        Mean    Std. Dev.       Min        Max

30. -------------+---------------------------------------------------------

31.         y1l |      5,000   -.0985646      .00288  -.1083639  -.0889075

32.       y1l_r |      5,000       .0544     .226828          0          1

33.        y1lp |      5,000   -.0939211    .0020038  -.1008612  -.0868043

34.      y1lp_r |      5,000       .6182    .4858765          0          1

35.         y2l |      5,000    .1084959     .065586  -.1065291   .3743112

36. -------------+---------------------------------------------------------

37.       y2l_r |      5,000       .0618     .240816          0          1

38.        y2lp |      5,000    .0915894     .055462  -.0975456   .3184061

39.      y2lp_r |      5,000       .0732    .2604906          0          1

同样地，对于真实 DGP 符合 Logit 模型的 MEM 和 TEM ，使用包含atmeans选项的margins命令。似然模型逼近真实似然，因此使用稳健标准误，并且使用选项vce以说明使用了两步 M 估计。两步 M 估计详见 Wooldridge (2010)  。

2.3 总结

使用 Probit 或 Logit 模型可以产生等价的边际效应，但线性概率模型的边际效应与 Probit 及 Logit 模型不同。
相关code

White, H. 1996. Estimation, Inference, and Specification Analysis>. Cambridge: Cambridge University Press.

Wooldridge, J. M. 2010. Econometric Analysis of Cross Section and Panel Data. 2nd ed. Cambridge, Massachusetts: MIT Press.

---

关于我们

• 【Stata 连享会(公众号：StataChina)】由中山大学连玉君老师团队创办，旨在定期与大家分享 Stata 应用的各种经验和技巧。

• 公众号推文同步发布于 CSDN-Stata连享会 、简书-Stata连享会 和 知乎-连玉君Stata专栏。可以在上述网站中搜索关键词Stata或Stata连享会后关注我们。

• 点击推文底部【阅读原文】可以查看推文中的链接并下载相关资料。

• Stata连享会 精彩推文1  || 精彩推文2

联系我们

• 欢迎赐稿： 欢迎将您的文章或笔记投稿至Stata连享会(公众号: StataChina)，我们会保留您的署名；录用稿件达五篇以上，即可免费获得 Stata 现场培训 (初级或高级选其一) 资格。

• 意见和资料： 欢迎您的宝贵意见，您也可以来信索取推文中提及的程序和数据。

• 招募英才： 欢迎加入我们的团队，一起学习 Stata。合作编辑或撰写稿件五篇以上，即可免费获得 Stata 现场培训 (初级或高级选其一) 资格。

• 联系邮件： [e-mail]

往期精彩推文

• Stata现场班 (2019.1.17-27) 课程详情

• Stata寒假精品现场班（2019.1月17日-27日）报名主页

---