Stata:一文读懂计数模型和因变量受限模型

发布时间:2021-09-10 阅读 6482

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下载 - 推文合集

作者: 杨雨萱 (中山大学)
邮箱:yangyuxuanxuan1994@163.com

 ​

编者按:本文主要摘译自下文,特此致谢!
Source:Data Analysis Examples -Link-


目录


1. 计数模型

1.1 泊松回归

有些被解释变量只能取非负整数,例如专利个数、子女人数等。对于这一类计数数据的分析,经常会使用到泊松回归 (poisson regression)。接下来,我们以一个具体的例子来讲解。在该模型中,因变量 num_awards 为一所高中学生获得的奖励数量,自变量包括学生类型 prog (如职业、普通或学术) 以及他们的数学期末考试成绩 math

. cnssc install lxhuse, lianxh replace
. lxhuse poisson_sim.dta, clear
. sum num_awards math

    Variable |  Obs    Mean    Std. dev.   Min   Max
-------------+--------------------------------------
  num_awards |  200     .63    1.052921      0     6
        math |  200  52.645    9.368448     33    75

. tabstat num_awards, by(prog) stats(mean sd n)

Summary for variables: num_awards
Group variable: prog (type of program)

    prog |   Mean        SD         N
---------+---------------------------
 general |     .2  .4045199        45
academic |      1  1.278521       105
vocation |    .24  .5174506        50
---------+---------------------------
   Total |    .63  1.052921       200
-------------------------------------

. histogram num_awards, discrete freq scheme(s1mono)
. poisson num_awards i.prog math, vce(robust)

Iteration 0:   log pseudolikelihood = -182.75759  
Iteration 1:   log pseudolikelihood = -182.75225  
Iteration 2:   log pseudolikelihood = -182.75225  
Poisson regression                                Number of obs =    200
                                                  Wald chi2(3)  =  80.15
                                                  Prob > chi2   = 0.0000
Log pseudolikelihood = -182.75225                 Pseudo R2     = 0.2118
------------------------------------------------------------------------
           |               Robust                                    
num_awards | Coefficient  std. err.    z    P>|z|   [95% conf. interval]
-----------+------------------------------------------------------------
      prog |                                                   
 academic  |      1.084      0.322   3.37   0.001      0.453       1.715
 vocation  |      0.370      0.401   0.92   0.357     -0.417       1.157
           |                                                   
      math |      0.070      0.010   6.71   0.000      0.050       0.091
     _cons |     -5.247      0.648  -8.10   0.000     -6.516      -3.978
------------------------------------------------------------------------

math 的系数是 0.070,表明数学分数每增加一单位,得奖的平均次数是原来的 e0.07 倍。academic 的系数为 1.084,表明与 general 组相比,academic 组的 log(奖励次数) 平均增加了约 1.084。类似地,vocation 组比 general 组的 log(奖励次数) 平均增加了约 0.370。

为了帮助评估模型拟合的效果,可以使用 estat gof 命令来获得拟合优度卡方检验。这不是对模型系数的检验,而是对模型形式的测试,即泊松模型形式是否适合我们的数据?

. estat gof

         Deviance goodness-of-fit =  189.4496
         Prob > chi2(196)         =    0.6182

         Pearson goodness-of-fit  =  212.1437
         Prob > chi2(196)         =    0.2040

可以看出,模型拟合得相当好,因为拟合优度卡方检验在统计上不显著。如果该检验在统计上具有显著性,则表明该数据与模型不很吻合。在这种情况下,我们可以尝试确定是否有遗漏的自变量、线性假设是否成立、以及否存在过度分散的问题。

有时,我们可能希望将回归结果表示为事件发生比率,我们可以使用 irr 选项。这些 IRR 值与上面输出的系数之间存在如下换算关系:IRR=ecoef

. poisson, irr

Poisson regression                           Number of obs =    200
                                             Wald chi2(3)  =  80.15
                                             Prob > chi2   = 0.0000
Log pseudolikelihood = -182.75225            Pseudo R2     = 0.2118
-------------------------------------------------------------------
           |           Robust                                      
num_awards |    IRR   std. err.    z    P>|z|  [95% conf. interval]
-----------+-------------------------------------------------------
      prog |                                                   
 academic  |  2.956      0.951   3.37   0.001     1.573       5.555
 vocation  |  1.447      0.581   0.92   0.357     0.659       3.179
           |                                                  
      math |  1.073      0.011   6.71   0.000     1.051       1.095
     _cons |  0.005      0.003  -8.10   0.000     0.001       0.019
-------------------------------------------------------------------
Note: _cons estimates baseline incidence rate.

上述结果表明,保持其他变量不变,academic 的平均得奖率是参考组 general 的 2.956 倍。同样,保持其他变量不变, vocation 的平均得奖率是参考组 general 的 1.447 倍。数学成绩每增加一个单位,得奖事件发生率的百分比就增加 7.3%。

回顾我们的模型:

上述模型可以该写为:

可以看出,系数在 log(y) 尺度上具有 加法效应,IRR 在 y 尺度上具有乘法效应。为了更好地理解这个模型,下面使用 margin 命令来计算 prog 的每一层级的预测计数,在模型中保留所有其他变量的平均值。

. margins prog, atmeans

Adjusted predictions                              Number of obs = 200
Model VCE: Robust
Expression: Predicted number of events, predict()
At: 1.prog =   .225 (mean)
    2.prog =   .525 (mean)
    3.prog =    .25 (mean)
    math   = 52.645 (mean)
---------------------------------------------------------------------
           |        Delta-method                                     
           | Margin   std. err.      z    P>|z|  [95% conf. interval]
-----------+---------------------------------------------------------
      prog |                                                    
  general  |  0.211      0.063     3.37   0.001     0.088       0.334
 academic  |  0.625      0.089     7.05   0.000     0.451       0.799
 vocation  |  0.306      0.083     3.69   0.000     0.144       0.468
---------------------------------------------------------------------

可以看出,在保持数学成绩在均值水平和 prog 取值为 1 水平上时,预测得到的平均得奖数大约是 0.211。同理,在保持数学成绩在均值水平和 prog 取值为 2 水平上时,预测得到的平均得奖数大约是为 0.625。在保持数学成绩在均值水平时和 prog 取值为 3 水平上时,预测得到的平均得奖数大约是 0.306。并且,prog 取值为 2 水平上的预测计数是 prog 取值为 1 水平上的 (.6249446/.211411)= 2.96 倍,这与我们在 IRR 输出表中看到的结果相一致。

下面我们将获得数学成绩从 35 到 75 (以 10 为增量) 的预测平均得奖数。

. margins, at(math=(35(10)75)) vsquish

Predictive margins                         Number of obs = 200
Model VCE: Robust
Expression: Predicted number of events, predict()
1._at: math = 35
2._at: math = 45
3._at: math = 55
4._at: math = 65
5._at: math = 75
--------------------------------------------------------------
     |        Delta-method                                    
     | Margin   std. err.     z    P>|z|  [95% conf. interval]
-----+--------------------------------------------------------
 _at |                                                   
  1  |  0.131      0.036    3.66   0.000     0.061       0.201
  2  |  0.264      0.048    5.57   0.000     0.171       0.358
  3  |  0.533      0.058    9.27   0.000     0.421       0.646
  4  |  1.076      0.122    8.82   0.000     0.837       1.315
  5  |  2.170      0.412    5.27   0.000     1.363       2.976
--------------------------------------------------------------

上表结果显示,当 prog 变量观测值固定时,当数学成绩是 35 时,平均得奖数约为 0.131;当数学成绩是 75 分时,平均得奖数约为 2.170;其他同理。如果我们比较 math = 35 和 math = 45的平均得奖数时,计算得到的比率是 (0.2644714/0.1311326) = 2.017,即成绩每增加 10 单位,后者对应的得奖数是前者的 2.017 倍,与 IRR 值计算的结果 1.0727^10 = 2.017 一致。

如果想多获得一些其他指标,可以使用用户编写的 fitstat 命令,如果想比较模型,这些信息可能会有所帮助。

. fitstat

Measures of Fit for poisson of num_awards
Log-Lik Intercept Only:     -231.864     Log-Lik Full Model:  -182.752
D(195):                      365.505     LR(3):                 98.223
                                         Prob > LR:              0.000
McFadden's R2:                 0.212     McFadden's Adj R2:      0.190
Maximum Likelihood R2:         0.388     Cragg & Uhler's R2:     0.430
AIC:                           1.878     AIC*n:                375.505
BIC:                        -667.667     BIC':                 -82.328

可以使用下面的命令来绘制事件的预测数量。

. twoway scatter c1 c2 c3 math, connect(l l l) sort ///
>        ytitle("Predicted Count") ylabel( ,nogrid) legend(rows(3)) ///
>        legend(ring(0) position(10)) scheme(s1mono)
=

该图表明,大多数奖项预计被学术计划 (prog=2) 中的学生拿走,尤其是在学生有更高数学成绩时。预计获奖人数最少的是普通计划课程的学生 (prog=1)。

1.2 负二项回归

负二项回归用于计数变量建模,通常用于过分散的计数结果变量。泊松回归的局限是泊松分布的期望与方差相等,被称为 “均等分散”,但这个特征与实际数据不符。

为此,Stata 命令提供了一个 LR 检验,其原假设为 “不存在过度分散,应使用泊松回归”。在完成负二项回归估计后,Stata 将自动输出检验结果 “LR test of alpha=0”,其中过度分散参数 “alpha=0” 对应于 “泊松回归”;即如果原假设 “alpha=0” 成立,则可使用泊松回归。如果使用选择项 “r” (稳健标准误),则 Stata 将只输出 alpha 的 95% 置信区间。根据这个置信区间,同样可以检验原假设 “alpha=0”。

下面我们以一个具体的例子,来看一下 Stata 具体操作流程,但由于和泊松回归过程类似,下面仅给出 Stata 命令。根据该数据构建的模型为,因变量为一两所高中学生的出勤天数 daysabs,自变量包括学生类型 prog (如职业、普通或学术) 以及他们的数学期末考试成绩 math

lxhuse nb_data.dta, clear
summarize daysabs math
histogram daysabs, discrete freq scheme(s1mono)
nbreg daysabs math i.prog
test 2.prog 3.prog
nbreg, irr
margins prog, atmeans
margins, at(math=(0(20)100)) vsquish
fitstat
predict c
sort math
twoway (line c math if prog==1) ///
       (line c math if prog==2) ///
       (line c math if prog==3), ///
       ytitle("Predicted Mean Value of Response") ylabel(0(3)15 ,nogrid) ///
       xtitle("Math Score") legend(order(1 "prog=1" 2 "prog=2" 3 "prog=3")) scheme(s1mono)

1.3 零膨胀泊松回归与零膨胀负二项回归

如果计数数据中含有大量的 0 值,则可考虑使用 “零膨胀泊松回归” 或 “零膨胀负二项回归”。此外,理论表明,多余的零可以看作是由一个单独的过程产生的计数值,多余的零可以独立建模。因此,零膨胀模型有两部分,泊松计数模型和用于预测多余零的 logit 模型。

究竟应该使用标准的泊松回归还是零膨胀泊松回归?Stata 提供了一个 “Vuong” 统计量,其渐近分布为标准正态。如果 “Vuong” 统计量很大 (为正数),则应该选择零膨胀泊松回归 (或零膨胀负二项回归);反之,如果 “Vuong” 统计量很小 (为负数),则应该选择标准的泊松回归 (或标准的负二项回归)。

需要注意的是,不建议将零膨胀模型应用于小样本。另外,零膨胀负二项回归的 Stata 操作流程与零膨胀泊松回归类似,读者可自行研究。

2.受限被解释变量

2.1 断尾回归

断尾回归主要用于,部分因变量的观测值由于某些原因未被纳入分析中。例如,yi 的总体为某地区所有企业的年销售收入,而统计局只收集有一定规模以上企业的数据,比如 yi100000。这样被解释变量在 100000 处就存在 “左边断尾”。接着,我们仍以一个简单例子来说明。

. lxhuse truncreg.dta, clear         
. summarize achiv langscore

  Variable |  Obs      Mean   Std. dev.  Min  Max
-----------+-------------------------------------
     achiv |  178  54.23596    8.96323    41   76
 langscore |  178  54.01124   8.944896    31   67

. tabstat achiv, by(prog) stats(n mean sd)

Summary for variables: achiv
Group variable: prog (type of program)
    prog |         N      Mean        SD
---------+------------------------------
 general |        40    51.575   7.97074
academic |       101  56.89109  9.018759
vocation |        37  49.86486  7.276912
---------+------------------------------
   Total |       178  54.23596   8.96323
----------------------------------------

. histogram achiv, bin(6) start(40) freq normal

下面我们使用 truncreg 命令来估计一个截断的回归模型。truncreg 命令中的 ll() 选项表示发生左截断的值,ul() 选项用于指示右截断值。

. truncreg achiv langscore i.prog, ll(40)
(0 obs truncated)

Fitting full model:
Iteration 0:   log likelihood = -598.11669  
Iteration 1:   log likelihood = -591.68358  
Iteration 2:   log likelihood = -591.31208  
Iteration 3:   log likelihood = -591.30981  
Iteration 4:   log likelihood = -591.30981  
Truncated regression
Limit: Lower =   40                        Number of obs     =        178
       Upper = +inf                        Wald chi2(3)      =      54.76
Log likelihood = -591.30981                Prob > chi2       =     0.0000
-------------------------------------------------------------------------
     achiv | Coefficient  Std. err.      z    P>|z|  [95% conf. interval]
-----------+-------------------------------------------------------------
 langscore |      0.713      0.114     6.22   0.000     0.488       0.937
           |                                                             
      prog |                                                             
 academic  |      4.065      2.055     1.98   0.048     0.038       8.093
 vocation  |     -1.136      2.670    -0.43   0.671    -6.369       4.097
           |                                                             
     _cons |     11.302      6.773     1.67   0.095    -1.973      24.576
-----------+-------------------------------------------------------------
    /sigma |      8.755      0.667    13.13   0.000     7.448      10.062
-------------------------------------------------------------------------

可以看出,langscore 系数为 0.713,表明语言成绩每增加一个单位,预期成绩就会增加 0.713 分。academic 系数为 4.065,表明在保持其他条件一定时,与 general 项目的同学相比,academic 项目的同学的预测成绩高了约 4.07 分。为了确定 prog 变量本身是否具有统计意义,我们可以使用 test 命令来获得该变量的二自由度检验。

. test 2.prog 3.prog

 (1)  [eq1]2.prog = 0
 (2)  [eq1]3.prog = 0
           chi2(  2) =    7.19
         Prob > chi2 =    0.0274

卡方检验表明,拒绝原假设,表明 prog 是具有统计学意义的预测因子。我们还可以使用 margins 命令来获得期望的边际效应。

. margins prog

Predictive margins                              Number of obs = 178
Model VCE: OIM
Expression: Linear prediction, predict()
-------------------------------------------------------------------
          |        Delta-method                                   
          | Margin   std. err.      z    P>|z| [95% conf. interval]
----------+--------------------------------------------------------
     prog |                                                        
 general  | 49.789      1.897    26.24   0.000   46.070      53.507
academic  | 53.854      1.150    46.83   0.000   51.600      56.108
vocation  | 48.653      2.140    22.73   0.000   44.458      52.848
--------------------------------------------------------------------

. marginsplot

在上表中的结果,我们可以看到,在 prog 等于 general 时,avchiv 的期望均值为 49.789。同理,在 prog 等于 academic 时,avchiv 的期望均值为 53.854。在 prog 等于 vocation 时,avchiv 的期望均值为 48.65。

truncreg 命令输出既不包含 R2 也不包含伪 R2。通过将 achiv 与预测值关联起来,然后将结果平方,从而计算出关联程度的粗略估计值。

. predict p
. correlate p achiv
. display r(rho)^2

.30519203

计算结果为 0.305,这个值是对回归中 R2 的粗略估计。观察到的 achiv 和预测的 achiv 之间的相关系数的平方约为 0.305,表明这些预测因素占结果变量的可变性的 30.5% 左右。

2.2 Tobit模型

tobit 模型又称归并回归 (censored regression) 模型,用于估计因变量中存在左归并或右归并时变量之间的线性关系(也称为上下归并)。具体来说,当某个值大于或等于某一阈值时,就会出现上述归并,因此真实值可能等于某一阈值,但也可能更高。那些高于或低于某个阈值的值将被归并。例如,在 20 世纪 80 年代,联邦法律限制速度计读数不超过每小时 85 英里。所以,如果你想通过马力和引擎大小来预测一辆车的最高时速,你得到的读数不会高于 85,即使这辆车的实际行驶速度更快。这是一个典型的数据右截 (右截) 的例子。我们唯一能确定的是这些车辆的时速至少为 85 英里。

我们以某个具体例子来说明,其中因变量为学业能力分数 apt,自变量为阅读成绩 read、数学成绩 math、以及学生类型 prog

. lxhuse tobitreg.dta, clear
. summarize apt read math

 Variable |   Obs      Mean    Std. dev.    Min    Max
----------+-------------------------------------------
      apt |   200   640.035    99.21903     352    800
     read |   200     52.23    10.25294      28     76
     math |   200    52.645    9.368448      33     75

. tabulate prog

    type of |
    program |      Freq.     Percent        Cum.
------------+-----------------------------------
   academic |         45       22.50       22.50
    general |        105       52.50       75.00
 vocational |         50       25.00      100.00
------------+-----------------------------------
      Total |        200      100.00

. histogram apt, normal bin(10) xline(800)

从上面的直方图中,我们可以看到数据中的归并 (censored)。需要注意的是,重点考察 y 在哪里被截断的,要有明确的来源和根据。接下来我们使用 tobit 命令进行回归分析,在 tobit 命令中的 ul() 选项表示右归并开始的值 (即上限),ll() 选项用于指示左归并值 (下限),在本例中不需要这个值。具体模型解释和 Truncated Regression 类似,下面仅提供代码,不做详细解释。

lxhuse tobitreg.dta, clear
tobit apt read math i.prog, ul(800)
test 2.prog 3.prog
test 2.prog = 3.prog
predict yhat
correlate apt yhat
fitstat

3. 参考资料

  • Data Analysis Examples -Link-
  • 陈强. 高级计量经济学及 Stata 应用[M]. 高等教育出版社, 2014.

4. 相关推文

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

相关课程

免费公开课

最新课程-直播课

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

课程主页

课程主页

关于我们

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

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

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

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

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