Stata：一文读懂计数模型和因变量受限模型
👇 连享会 · 推文导航 | www.lianxh.cn
🍎 Stata：Stata基础 | Stata绘图 | Stata程序 | Stata新命令 📘 论文：数据处理 | 结果输出 | 论文写作 | 数据分享 💹 计量：回归分析 | 交乘项-调节 | IV-GMM | 时间序列 | 面板数据 | 空间计量 | Probit-Logit | 分位数回归 ⛳ 专题：SFA-DEA | 生存分析 | 爬虫 | 机器学习 | 文本分析 🔃 因果：DID | RDD | 因果推断 | 合成控制法 | PSM-Matching 🔨 工具：工具软件 | Markdown | Python-R-Stata 🎧 课程：公开课-直播 | 计量专题 | 关于连享会
连享会视频课 · 因果推断实用计量方法
作者： 杨雨萱 (中山大学)
邮箱：yangyuxuanxuan1994@163.com
编者按：本文主要摘译自下文，特此致谢！
Source：Data Analysis Examples -Link-
温馨提示： 文中链接在微信中无法生效。请点击底部「阅读原文」。或直接长按/扫描如下二维码，直达原文：
目录
1. 计数模型
1.1 泊松回归
1.2 负二项回归
1.3 零膨胀泊松回归与零膨胀负二项回归
2.受限被解释变量
2.1 断尾回归
2.2 Tobit模型
3. 参考资料
4. 相关推文
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，表明数学分数每增加一单位，得奖的平均次数是原来的 倍。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 值与上面输出的系数之间存在如下换算关系：。
. 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 断尾回归
断尾回归主要用于，部分因变量的观测值由于某些原因未被纳入分析中。例如， 的总体为某地区所有企业的年销售收入，而统计局只收集有一定规模以上企业的数据，比如 。这样被解释变量在 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 命令输出既不包含 也不包含伪 。通过将 achiv 与预测值关联起来，然后将结果平方，从而计算出关联程度的粗略估计值。
. predict p
. correlate p achiv
. display r(rho)^2
.30519203
计算结果为 0.305，这个值是对回归中 的粗略估计。观察到的 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
专题：Stata命令 Stata新命令：面板-LogitFE-ProbitFE 专题：交乘项-调节 Logit-Probit中的交乘项及边际效应图示 专题：Probit-Logit Logit-Probit：非线性模型中交互项的边际效应解读 秒懂小罗肥归：logit与mlogit详解 reg2logit：用OLS估计Logit模型参数 feologit：固定效应有序Logit模型 Stata：多元 Logit 模型详解 (mlogit) Stata：Logit模型一文读懂 详解 Logit/Probit 模型中的 completely determined 问题 Stata：Logit 模型评介 二元选择模型：Probit 还是 Logit？ Stata：何时使用线性概率模型而非Logit？ Stata：嵌套 Logit 模型 (Nested Logit) Stata：二元Probit模型 动态 Probit 模型及 Stata 实现
课程推荐：因果推断实用计量方法
主讲老师：邱嘉平教授
🍓 课程主页：https://gitee.com/lianxh/YGqjp
New！ Stata 搜索神器：
lianxh和songblGIF 动图介绍
搜： 推文、数据分享、期刊论文、重现代码 ……
👉 安装：
. ssc install lianxh
. ssc install songbl
👉 使用：
. lianxh DID 倍分法
. songbl all
🍏 关于我们
连享会 ( www.lianxh.cn，推文列表) 由中山大学连玉君老师团队创办，定期分享实证分析经验。 直通车： 👉【百度一下： 连享会】即可直达连享会主页。亦可进一步添加 「知乎」,「b 站」,「面板数据」,「公开课」 等关键词细化搜索。