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 系列推文：

全部 | Stata入门 | Stata教程 | Stata资源 | Stata命令
计量专题 | 论文写作 | 数据分享 | 专题课程
结果输出 | Stata绘图 | 数据处理 | Stata程序
回归分析 | 面板数据 | 交乘项-调节 | IV-GMM
内生性-因果推断 | 倍分法DID | 断点回归RDD | PSM-Matching | 合成控制法
Probit-Logit | 时间序列 | 空间计量 | 分位数回归 | 生存分析 | SFA-DEA
文本分析-爬虫 | Python-R-Matlab | 机器学习
Markdown | 工具软件 | 其它

☝ 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，表明数学分数每增加一单位，得奖的平均次数是原来的 $e^{0.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 值与上面输出的系数之间存在如下换算关系： $I R R = e^{c o e f}$ 。

. 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%。

回顾我们的模型：

l o g (n u m_a w a r d s) = I n t e r c e p t + b 1 (p r o g = 2) + b 2 (p r o g = 3) + b 3 m a t h

上述模型可以该写为：

\begin{aligned} n u m_a w a r d s & = e x p (I n t e r c e p t + b 1 (p r o g = 2) + b 2 (p r o g = 3) + b 3 m a t h) \\ = e x p (I n t e r c e p t) * e x p (b 1 (p r o g = 2)) * e x p (b 2 (p r o g = 3)) * e x p (b 3 m a t h) \end{aligned}

可以看出，系数在 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 断尾回归

断尾回归主要用于，部分因变量的观测值由于某些原因未被纳入分析中。例如， $y_{i}$ 的总体为某地区所有企业的年销售收入，而统计局只收集有一定规模以上企业的数据，比如 $y_{i} \geq 100000$ 。这样被解释变量在 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 命令输出既不包含 $R^{2}$ 也不包含伪 $R^{2}$ 。通过将 achiv 与预测值关联起来，然后将结果平方，从而计算出关联程度的粗略估计值。

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

.30519203

计算结果为 0.305，这个值是对回归中 $R^{2}$ 的粗略估计。观察到的 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小时]

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

1. 计数模型

1.1 泊松回归

1.2 负二项回归

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

2.受限被解释变量

2.1 断尾回归

2.2 Tobit模型

3. 参考资料

4. 相关推文

相关课程

免费公开课

最新课程-直播课

⛳ 课程主页

⛳ 课程主页

关于我们

热门资讯