温馨提示: 定期 清理浏览器缓存,可以获得最佳浏览体验。
New!
lianxh
命令发布了:
随时搜索连享会推文、Stata 资源,安装命令如下:
. ssc install lianxh
使用详情参见帮助文件 (有惊喜):
. help lianxh
作者:郑浩文(中山大学)
邮箱: zhenghw25@mail2.sysu.edu.cn
目录
现实中,存在一类模型,其解释变量是离散的,如事件发生次数、物品个数等,这类模型被称为计数模型。
计数模型一般使用现有分布来拟合解释变量的频率,根据数据特征可使用正态分布、泊松分布、负二项分布等。而因为解释变量非连续,且当解释变量取值有明确上下界(如非负约束)时,正态分布通常不适用,而泊松分布等将有更好的拟合效果。计数模型的估计往往使用极大似然估计。
考虑以泊松分布来建模。根据泊松分布,随机变量取值为
为保证
以下 Stata 代码可以生成一个满足泊松分布的变量 poi。
. set obs 100000 // 设置 100000 的样本量
. set seed 1234 // 设置种子值
. gen poi = rpoisson(4) // λ 为 4 的泊松分布随机
. hist poi ,fraction xtitle("yi") ytitle("Pr") // 绘制直方图
分布如图:
观察上图可以发现,取值的概率随着取值的增大而先增大后减小。然而,在一些场景下,0 取值的频率很高,无法使用普通泊松分布解释,于是出现了所谓零膨胀问题。如在车险理赔次数的频率拟合中,存在许多的零赔付的保单。这可能由两部分原因组成:一是车祸发生的概率较低,二是车主对一些影响较小的擦碰选择不索赔,以免续保时保费增加。
为此,需要对泊松回归模型进行一定的修正。以车险为例,假设实际发生事故并可理赔的次数满足泊松分布,而车主在发生事故后有
其中,
同时,方程可以等价地表示为两个相互独立随机变量乘积的分布,一个变量满足 0-1 两点分布,另一个满足泊松分布:
因此可以将 ZIP 分为两个部分:
以 ZIP 的基本思想为基础,可以拓展出
通过以下 Stata 代码可以构造零膨胀的数据。方法是先构造基分布(两点分布和泊松分布),再对泊松分布进行零膨胀处理。
. set obs 100000 // 设置样本量
. set seed 1234 // 设置随机数种子
. gen poi = rpoisson(4) // 参数为 4 的泊松分布
. gen assign = runiform(0,100) // 分配变量 assign 满足 0-100 的均匀分布
. scalar separate = 12 // 设置分配变量阈值
. gen num = poi // 变量 num 满足零膨胀后的分布
. replace num = 0 if assign<=separate // 分配变量<=12的 变量设置为 0
. * 相当于两点分布随机变量取 0 值的概率为 12%
num 的直方图如下,曲线为泊松分布:
翻译并修改自 ZERO-INFLATED POISSON REGRESSION | STATA DATA ANALYSIS EXAMPLES。原文使用 Stata 12,本文使用 Stata 15 重现结果,差异之处通过译注形式说明。
离散的计数数据常用泊松回归模型来拟合。若数据集中包含过多的零时,数据的分布就与标准泊松分布有较大偏差,此时应该进行修正。理论表明,数据中零的产生是一个相对独立的过程,可以与其他数据分开建模。由此,零膨胀泊松回归模型( ZIP )应运而生。模型分为两个子模型:Logit模型和泊松计数模型。前者用于拟合零值的频率,后者用于拟合其他值的频率。关于这两个模型的使用,可以参考以下两篇数据分析样例介绍:泊松回归 、Logit回归 。
注:本文意在展示如何使用各种数据分析命令,并不完全覆盖研究者应有研究过程的各个方面。例如,文章不包含数据清洗和检查、假设检验、模型诊断和其他后续可能进行的分析过程。
适用的样例
例 1. 学校管理者研究了两所高中的高二学生在某一个学期中的出勤表现。出勤指标的定义是缺勤天数,并使用性别、三科(数学、语文和艺术)标准化后的成绩共四个变量进行预测。大多数学生整个学期都没有缺勤记录,意味着零值的频率很高。
例 2. 州野生生物学家想对一个州立公园内渔夫的捕鱼量进行建模。他们向游客收集以下信息:是否有野营车、团队中的人数、团队内是否有小孩和捕鱼数量。数据中样本捕鱼量为零原因有两种:一是游客没有去抓鱼,二是没有抓到鱼。前者是导致数据零值频率偏高的原因,而数据中并没有包含是否有去抓鱼的信息,因此无法将两者分开。
数据介绍
让我们用 ZIP 模型来探究上述的样例 2,使用的数据集为 Stata 手册附带数据 fish.dta。该数据集中一共有 250 群游客,每一群游客都被问及:捕鱼量(count)、有几个孩子(child)、总人数(persons)和是否有野营车(camper)。
在预测捕鱼量的同时,我们还可以预测超额零值的频率,即有多少零值不是因为抓不到鱼而产生的。
*. use http://www.stata-press.com/data/r10/fish, clear
. webuse "fish.dta", clear
. summarize count child persons camper
Variable | Obs Mean Std. Dev. Min Max
----------+---------------------------------------
count | 250 3.296 11.63503 0 149
child | 250 .684 .8503153 0 3
persons | 250 2.528 1.11273 1 4
camper | 250 .588 .4931824 0 1
. histogram count, discrete freq
从上面的直方图可以看出,count=0 频数接近 150,约占样本数的 60% (150/250),意味着 零膨胀 问题很明显。我们也可以使用 tabulate
命令更为细致地查看 count 变量的频数分布:
. tab count
number of |
fish caught | Freq. Percent Cum.
------------+-----------------------------------
0 | 142 56.80 56.80
1 | 31 12.40 69.20
2 | 20 8.00 77.20
3 | 12 4.80 82.00
4 | 6 2.40 84.40
5 | 10 4.00 88.40
6 | 4 1.60 90.00
7 | 3 1.20 91.20
8 | 2 0.80 92.00
9 | 2 0.80 92.80
10 | 1 0.40 93.20
11 | 1 0.40 93.60
13 | 1 0.40 94.00
14 | 1 0.40 94.40
15 | 2 0.80 95.20
16 | 1 0.40 95.60
21 | 2 0.80 96.40
22 | 1 0.40 96.80
29 | 1 0.40 97.20
30 | 1 0.40 97.60
31 | 1 0.40 98.00
32 | 2 0.80 98.80
38 | 1 0.40 99.20
65 | 1 0.40 99.60
149 | 1 0.40 100.00
------------+-----------------------------------
Total | 250 100.00
可见,有 142 位游客的捕鱼量 (count) 为 0,而在 捕到鱼 的人群中,count 的取值集中在 1-5 之间,占比约为 17.6%。count 的最大值为 149,导致 count 变量的分布极为分散 (从上面的直方图可以非常直观地看到这一特征)。简言之,捕鱼量 (count) 的分布兼具了 零膨胀 和 过度分散 的特征。
我们也可以看看其他特征变量的统计特征:
. tab1 child persons camper
-> tabulation of child
child | Freq. Percent Cum.
------------+-----------------------------------
0 | 132 52.80 52.80
1 | 75 30.00 82.80
2 | 33 13.20 96.00
3 | 10 4.00 100.00
------------+-----------------------------------
Total | 250 100.00
-> tabulation of persons
persons | Freq. Percent Cum.
------------+-----------------------------------
1 | 57 22.80 22.80
2 | 70 28.00 50.80
3 | 57 22.80 73.60
4 | 66 26.40 100.00
------------+-----------------------------------
Total | 250 100.00
-> tabulation of camper
camper | Freq. Percent Cum.
------------+-----------------------------------
0 | 103 41.20 41.20
1 | 147 58.80 100.00
------------+-----------------------------------
Total | 250 100.00
可考虑的分析方法
下列方法都可用于数据的分析。有的方法比较合理,而其余的不是不被研究者偏好,就是应用条件有限制。
零膨胀泊松回归
我们运行 zip 命令,以变量 child 和 camper 作 count 的解释变量,persons 作为过多零值的解释变量。我们使用了 vuong 选项,可以检验需要使用零膨胀模型还是普通计数模型。
. zip count child camper, inflate(persons) forcevuong
Fitting constant-only model:
Iteration 0: log likelihood = -1347.807
Iteration 1: log likelihood = -1315.5343
Iteration 2: log likelihood = -1126.3689
Iteration 3: log likelihood = -1125.5358
Iteration 4: log likelihood = -1125.5357
Iteration 5: log likelihood = -1125.5357
Fitting full model:
Iteration 0: log likelihood = -1125.5357
Iteration 1: log likelihood = -1044.8553
Iteration 2: log likelihood = -1031.8733
Iteration 3: log likelihood = -1031.6089
Iteration 4: log likelihood = -1031.6084
Iteration 5: log likelihood = -1031.6084
Zero-inflated Poisson regression Number of obs = 250
Nonzero obs = 108
Zero obs = 142
Inflation model = logit LR chi2(2) = 187.85
Log likelihood = -1031.608 Prob > chi2 = 0.0000
------------------------------------------------------------------------------
count | Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
count |
child | -1.042838 .0999883 -10.43 0.000 -1.238812 -.846865
camper | .8340222 .0936268 8.91 0.000 .650517 1.017527
_cons | 1.597889 .0855382 18.68 0.000 1.430237 1.76554
-------------+----------------------------------------------------------------
inflate |
persons | -.5643472 .1629638 -3.46 0.001 -.8837503 -.244944
_cons | 1.297439 .3738522 3.47 0.001 .5647022 2.030176
------------------------------------------------------------------------------
Vuong test of zip vs. standard Poisson: z = 3.57 Pr>z = 0.0002
Warning: The Vuong test is not appropriate for testing zero-inflation.
以上输出看起来很像 OLS 回归的输出:
Cameron 和 Trivedi( 2009 ) 建议使用泊松模型的稳健标准误。我们加入 vce(robust) 选项后重新运行模型。我们不能在上述命令后加入该选项,因为 robust 和 vuong 选项不能同时出现在一个模型中(见译注 2 )。
. zip count child i.camper, inflate(persons) vce(robust)
Fitting constant-only model:
Iteration 0: log pseudolikelihood = -1347.807
Iteration 1: log pseudolikelihood = -1315.5343
Iteration 2: log pseudolikelihood = -1126.3689
Iteration 3: log pseudolikelihood = -1125.5358
Iteration 4: log pseudolikelihood = -1125.5357
Iteration 5: log pseudolikelihood = -1125.5357
Fitting full model:
Iteration 0: log pseudolikelihood = -1125.5357
Iteration 1: log pseudolikelihood = -1044.8553
Iteration 2: log pseudolikelihood = -1031.8733
Iteration 3: log pseudolikelihood = -1031.6089
Iteration 4: log pseudolikelihood = -1031.6084
Iteration 5: log pseudolikelihood = -1031.6084
Zero-inflated Poisson regression Number of obs = 250
Nonzero obs = 108
Zero obs = 142
Inflation model = logit Wald chi2(2) = 7.25
Log pseudolikelihood = -1031.608 Prob > chi2 = 0.0266
------------------------------------------------------------------------------
| Robust
count | Coef. Std. Err. z P>|z| [95% Conf. Interval]
-------------+----------------------------------------------------------------
count |
child | -1.042838 .3893772 -2.68 0.007 -1.806004 -.2796731
1.camper | .8340222 .4076029 2.05 0.041 .0351352 1.632909
_cons | 1.597889 .2934631 5.44 0.000 1.022711 2.173066
-------------+----------------------------------------------------------------
inflate |
persons | -.5643472 .2888849 -1.95 0.051 -1.130551 .0018567
_cons | 1.297439 .493986 2.63 0.009 .3292445 2.265634
------------------------------------------------------------------------------
现在我们探讨结果的细节:
我们用 margins 命令( Stata 11 引入)来帮助理解我们的模型。首先我们计算类别变量 camper 的期望捕鱼量,计算时使用 atmeans 选项将 child 保持在其均值处。
. margins camper, atmeans
Adjusted predictions Number of obs = 250
Model VCE : Robust
Expression : Predicted number of events, predict()
at : child = .684 (mean)
0.camper = .412 (mean)
1.camper = .588 (mean)
persons = 2.528 (mean)
-------------------------------------------------------------------------
| Delta-method
| Margin Std. Err. z P>|z| [95% Conf. Interval]
--------+----------------------------------------------------------------
camper |
0 | 1.289132 .4393168 2.93 0.003 .4280866 2.150177
1 | 2.968305 .619339 4.79 0.000 1.754423 4.182187
-------------------------------------------------------------------------
结果显示,在 child persons 均值处,无野营车的游客的期望捕鱼量为 1.289 ,有野营车的为 2.968。
使用 dydx 选项可以计算 camper=0 和 camper=1 之间的期望捕鱼量差异。child 和 persons 仍保持在均值处,分别为 .684 和 2.528 。
. margins, dydx(camper) atmeans
Conditional marginal effects Number of obs = 250
Model VCE : Robust
Expression : Predicted number of events, predict()
dy/dx w.r.t. : 1.camper
at : child = .684 (mean)
0.camper = .412 (mean)
1.camper = .588 (mean)
persons = 2.528 (mean)
-----------------------------------------------------------------------
| Delta-method
| dy/dx Std. Err. z P>|z| [95% Conf. Interval]
----------+------------------------------------------------------------
1.camper | 1.679173 .7754611 2.17 0.030 .1592975 3.199049
-----------------------------------------------------------------------
Note: dy/dx for factor levels is the discrete change from the base level.
结果显示,有无野营车的差距为 1.679,是统计上显著的(译注:上一条命令结果中两类别差也是 1.679,使用 dydx 可以检验差异的统计显著性 )。
我们最后一条 margins 命令将给出 camper 分别取 0/1 、 child 变量分别取 0 到 3 时的期望捕鱼量。
. margins, at(child=(0(1)3) camper=(0/1)) vsquish
Predictive margins Number of obs = 250
Model VCE : Robust
Expression : Predicted number of events, predict()
1._at : child = 0
camper = 0
2._at : child = 0
camper = 1
3._at : child = 1
camper = 0
4._at : child = 1
camper = 1
5._at : child = 2
camper = 0
6._at : child = 2
camper = 1
7._at : child = 3
camper = 0
8._at : child = 3
camper = 1
-----------------------------------------------------------------------
| Delta-method
| Margin Std. Err. z P>|z| [95% Conf. Interval]
------+----------------------------------------------------------------
_at |
1 | 2.616441 .6470522 4.04 0.000 1.348242 3.88464
2 | 6.024516 2.159288 2.79 0.005 1.79239 10.25664
3 | .922172 .4142303 2.23 0.026 .1102954 1.734048
4 | 2.123358 .4771534 4.45 0.000 1.188154 3.058561
5 | .3250221 .2611556 1.24 0.213 -.1868335 .8368777
6 | .7483834 .3929987 1.90 0.057 -.0218798 1.518647
7 | .114555 .1351887 0.85 0.397 -.1504101 .37952
8 | .2637699 .2365495 1.12 0.265 -.1998587 .7273984
-----------------------------------------------------------------------
. marginsplot
Variables that uniquely identify margins: child camper
结果显示,期望捕鱼量对是否有野营车的游客来说,都会随着孩子数量的增加而减少。
使用 fitstat 命令可以展示许多的拟合统计量。此命令是 J. Scott Long 和 Jeremy Freese 编写的 spostado 包的一部分( 通过 search spostado 查找并获取 )。
. fitstat
Measures of Fit for zip of count
Log-Lik Intercept Only: -1127.023 Log-Lik Full Model: -1031.608
D(244): 2063.217 LR(4): 190.829
Prob > LR: 0.000
McFadden's R2: 0.085 McFadden's Adj R2: 0.079
Maximum Likelihood R2: 0.534 Cragg & Uhler's R2: 0.534
AIC: 8.301 AIC*n: 2075.217
BIC: 715.980 BIC': -168.743
注意事项:
更多信息
参考文献
译注 1 : 原文提供的是 vuong 选项。而在 Stata 15 中,命令报错,提示已经不建议使用 vuong 检验来选择模型。本文使用 forcevuong 强制进行 vuong 检验。同样的, forcevuong 和 vce(robust) 不能共存。
. zip count child camper, inflate(persons) vuong
Vuong test is not appropriate for testing zero-inflation. Specify option forcevuong to perform the test
anyway.
r(498);
Stata 提供的依据是 2015 年 Wilson 的论文 The misuse of the Vuong test for non-nested models to test for zero-inflation。Vuong(1989)提出的 Vuong 非嵌套检验适用于非嵌套模型。在此条件下,两个严格非嵌套模型的对数似然比分布满足正态性。而 Wilson 通过模拟发现, zip 模型与标准泊松回归的对数似然比分布不满足正态性,说明 Vuong 检验不适用于零膨胀的检验。其原因在于 ZIP 和 标准泊松回归是嵌套模型,即当 (1) 式中的
此外,Wilson 还关注了一个细节。根据 Logit 分布,
令
. help j_vuong##_new
Vuong test is not appropriate for testing zero inflation
Greene (1994) proposed using the Vuong test for nonnested models (Vuong 1989) to test for zero inflation. Despite many earlier citations, recent work by Wilson (2015) has shown that the Vuong test is inappropriate for testing zero inflation. Nesting occurs when the probability of zero inflation is 0, which is on the boundary, and this violates the regularity conditions of the Vuong test for nonnested models. So the distribution of the test statistic is not standard normal. The actual distribution is unknown and thus cannot be used for inference.
You may consider using information criteria to choose between the standard and the zero-inflated models; see example 2 in [R] zip.
...
Van den Broek(1995)提出用 Score 检验来判别数据是否需要使用零膨胀模型。在 Stata 官方帮助文档展示的 zip example 2 中,Stata 建议使用 AIC、BIC 准则进行模型选择,若 ZIP 模型的 AIC、BIC 比标准泊松回归小,则零膨胀模型更优。
以下代码将使用 AIC、BIC 准则进行模型选择。可以看出,零膨胀模型的指标值更小,因此零膨胀模型更优。
. * zip 回归 并储存结果
. cap zip count child camper, inflate(persons) vce(robust)
. est store zpvce
. * 标准泊松回归 并储存结果
. cap poisson count child camper,vce(robust)
. est store poi
. * 比较统计量
. est stats poi zpvce
Akaike's information criterion and Bayesian information criterion
--------------------------------------------------------------
Model | Obs ll(null) ll(model) df AIC BIC
--------+-----------------------------------------------------
poi | 250 -1647.716 -1358.593 3 2723.186 2733.75
zpvce | 250 -1125.536 -1031.608 5 2073.217 2090.824
--------------------------------------------------------------
Note: N=Obs used in calculating BIC; see [R] BIC note.
译注 2 : 也可以使用第一次 zip
命令的变量设定:
. zip count child camper, inflate(persons) vce(robust)
差异在于 camper 在回归时不作为虚拟变量。其运行结果完全一致,这是因为数据中 camper 变量本为 0/1 变量。然而从变量含义上看,使用 i.camper 更为规范。
陈强. 高级计量经济学及 Stata 应用[M]. 高等教育出版社, 2014.
徐昕.(2020). 零膨胀广义泊松模型的推广及其在精算中应用. 数理统计与管理(06),1010-1021. doi:10.13860/j.cnki.sltj.20200930-007. -Link-
Wilson, P. 2015. The misuse of the Vuong test for non-nested models to test for zero-inflation. Economics Letters 127:51-53. -Link-, -PDF-
Van den Broek, J. (1995). A Score Test for Zero Inflation in a Poisson Distribution. Biometrics, 51(2), 738-743. -Link-, -PDF-
Vuong, Q. H. (1989). Likelihood Ratio Tests for Model Selection and Non-Nested Hypotheses. Econometrica, 57(2), 307–333. -Link-, -PDF-
Note:产生如下推文列表的命令为:
lianxh logit probit, m
安装最新版lianxh
命令:
ssc install lianxh, replace
连享会-直播课 上线了!
http://lianxh.duanshu.com
免费公开课:
直击面板数据模型 - 连玉君,时长:1小时40分钟,课程主页 Stata 33 讲 - 连玉君, 每讲 15 分钟. Stata 小白的取经之路 - 龙志能,时长:2 小时,课程主页 部分直播课 课程资料下载 (PPT,dofiles等)
支持回看
专题 | 嘉宾 | 直播/回看视频 |
---|---|---|
⭐ 最新专题 | 因果推断, 空间计量,寒暑假班等 | |
⭕ 数据清洗系列 | 游万海 | 直播, 88 元,已上线 |
研究设计 | 连玉君 | 我的特斯拉-实证研究设计,-幻灯片- |
面板模型 | 连玉君 | 动态面板模型,-幻灯片- |
面板模型 | 连玉君 | 直击面板数据模型 [免费公开课,2小时] |
Note: 部分课程的资料,PPT 等可以前往 连享会-直播课 主页查看,下载。
关于我们
课程, 直播, 视频, 客服, 模型设定, 研究设计, stata, plus, 绘图, 编程, 面板, 论文重现, 可视化, RDD, DID, PSM, 合成控制法
等
连享会小程序:扫一扫,看推文,看视频……
扫码加入连享会微信群,提问交流更方便
✏ 连享会学习群-常见问题解答汇总:
✨ https://gitee.com/arlionn/WD
New!
lianxh
命令发布了:
随时搜索连享会推文、Stata 资源,安装命令如下:
. ssc install lianxh
使用详情参见帮助文件 (有惊喜):
. help lianxh