Stata: 用负二项分布预测蚊子存活率

作者：崔颖（中央财经大学）

Stata 连享会：知乎 | 简书 | 码云 | CSDN

特别说明

文中包含的链接在微信中无法生效。请点击本文底部左下角的【阅读原文】，转入本文【简书版】。

Source: Survival Rates and Negative Binomial (Francis, 2012)

炎炎夏日，你是否也深受蚊虫叮咬困扰？

本篇推文将帮助你用 Stata 数值模拟预测在给定不同的死亡率(可能被捕食或遭遇意外等)情况下，蚊子的平均存活期限及存活概率。

若给定失败概率和淘汰所需失败次数，则预期存活天数服从负二项分布，这里就是第1次“死亡事件”出现时所需的天数服从负二项分布。根据负二项分布定义及概率密度函数可知，预期存活天数的期望为

其中，r=1.

给定若干死亡率：

死亡率=0.01时，p=0.99, ；

死亡率=0.015时，p=0.985,；

死亡率=0.02时，p=0.98, ；

死亡率=0.025时，p=0.975, ；

死亡率=0.03时，p=0.97, ；

死亡率=0.035时，, .

1. 数据生成

首先设置 12,000 个观测值。

1. clear

2. set obs 12000

将其分为6组，每组200个观测值。6组观测值的死亡率取值分别为0.01，0.015，0.02，0.025，0.03，0.035。

1. gen mort_rate=mod(_n,6)/200+0.01

2. tab mort_rate

生成一个哑变量表示蚊子尚且存活：

1. gen living=1

生成一个连续变量表示蚊子存活的天数：

1. gen survived=0

2. 数值模拟

假设观测了 500 天，看每天还剩多少蚊子存活。

1. forvalue i=1/500{

2. *随机模拟伯努利试验，生成临时变量`died`表示特定某一天蚊子是否死亡

3. qui gen died=rbinomial(1,mort_rate) if living==1

4. *如果蚊子死亡，则替换其存活状态变量

5. qui replace living=0 if living==1 & died == 1

6. *如果蚊子还活着，替换survived变量来表示当前是第几天

7. qui replace survived=`i' if living==1

8. *每次循环报告出当前循环是第几天，以及当天蚊子尚存活的比例

9. qui sum living

10. di "Round `i' :" r(mean)

11. *删除临时变量`died`

12. drop died

13. }

2.1 绘制图像

根据不同的死亡率分组绘制蚊子存活天数的直方图。

1. hist survived, by(mort_rate)

从图上可以看出，死亡率越高的组别，蚊子存活数量衰减越快。

2.2 不同死亡率下的平均存活率

此外，我们还可以根据不同的死亡率分组看 suivived 变量的描述性统计信息。

1. bysort mort_rate: sum survived

1. ----------------------------------------------------------------------------------

2. -> mort_rate = .01

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

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

5. survived | 2,000 100.424 96.9187 0 500

6. ----------------------------------------------------------------------------------

7. -> mort_rate = .015

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

9. -------------+---------------------------------------------------------

10. survived | 2,000 66.1455 67.45336 0 470

11. ----------------------------------------------------------------------------------

12. -> mort_rate = .02

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

14. -------------+---------------------------------------------------------

15. survived | 2,000 49.517 50.21277 0 408

16. ----------------------------------------------------------------------------------

17. -> mort_rate = .025

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

19. -------------+---------------------------------------------------------

20. survived | 2,000 40.9 40.77116 0 302

21. ----------------------------------------------------------------------------------

22. -> mort_rate = .03

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

24. -------------+---------------------------------------------------------

25. survived | 2,000 32.2755 33.19746 0 276

26. ----------------------------------------------------------------------------------

27. -> mort_rate = .035

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

29. -------------+---------------------------------------------------------

30. survived | 2,000 28.033 28.43426 0 215

从上表可得，6组观测值存活天数的数值模拟均值均与本文开篇通过负二项分布概率密度函数计算所得的期望近似，验证了结果的可靠性。

写在后面

本文通过一个有趣的例子介绍了用 Stata 做数值模拟的方法，步骤可以大致概括为：确定分布函数和相关参数，生成观测值，循环模拟随机试验，绘图呈现结果。

关于我们

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

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

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

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

联系我们

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

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

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

• 联系邮件： [e-mail]

往期精彩推文

---

欢迎加入Stata连享会(公众号: StataChina)

2019暑期“实证研究方法与经典论文”专题班