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

发布时间：2020-10-09 阅读 2799

Stata 连享会主页 || 视频 || 推文

温馨提示： 定期清理浏览器缓存，可以获得最佳浏览体验。

✌ 课程详情： https://gitee.com/arlionn/Course | lianxh.cn

✨ 课程主页： https://gitee.com/arlionn/Course

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

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

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

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

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

\frac{r p}{1 - p}

其中， $r = 1$ .

给定若干死亡率：

死亡率=0.01时， $p = 0.99$ , $E = 0.99 / 0.01$ ；

死亡率=0.015时， $p = 0.985$ , $E = 0.985 / 0.015$ ；

死亡率=0.02时， $p = 0.98$ , $E = 0.98 / 0.02$ ；

死亡率=0.025时， $p = 0.975$ , $E = 0.975 / 0.025$ ；

死亡率=0.03时， $p = 0.97$ , $E = 0.97 / 0.03$ ；

死亡率=0.035时， $p = 0.965$ , $E = 0.99 / 0.035$ .

1. 数据生成

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

clear 
set obs 12000

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

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

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

gen living=1

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

gen survived=0

连享会计量方法专题……

2. 数值模拟

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

forvalue i=1/500{
  *随机模拟伯努利试验，生成临时变量`died`表示特定某一天蚊子是否死亡
  qui gen died=rbinomial(1,mort_rate) if living==1
  *如果蚊子死亡，则替换其存活状态变量
  qui replace living=0 if living==1 & died == 1
  *如果蚊子还活着，替换survived变量来表示当前是第几天
  qui replace survived=`i' if living==1
  *每次循环报告出当前循环是第几天，以及当天蚊子尚存活的比例
  qui sum living
   di "Round `i' :" r(mean)
  *删除临时变量`died`
   drop died
}

2.1 绘制图像

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

hist survived, by(mort_rate)

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

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

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

bysort mort_rate: sum survived

-----------------------------------------------------------------------
-> mort_rate = .01

    Variable |        Obs        Mean    Std. Dev.       Min        Max
-------------+---------------------------------------------------------
    survived |      2,000     100.424     96.9187          0        500

-----------------------------------------------------------------------
-> mort_rate = .015

    Variable |        Obs        Mean    Std. Dev.       Min        Max
-------------+---------------------------------------------------------
    survived |      2,000     66.1455    67.45336          0        470

-----------------------------------------------------------------------
-> mort_rate = .02

    Variable |        Obs        Mean    Std. Dev.       Min        Max
-------------+---------------------------------------------------------
    survived |      2,000      49.517    50.21277          0        408

-----------------------------------------------------------------------
-> mort_rate = .025

    Variable |        Obs        Mean    Std. Dev.       Min        Max
-------------+---------------------------------------------------------
    survived |      2,000        40.9    40.77116          0        302

-----------------------------------------------------------------------
-> mort_rate = .03

    Variable |        Obs        Mean    Std. Dev.       Min        Max
-------------+---------------------------------------------------------
    survived |      2,000     32.2755    33.19746          0        276

-----------------------------------------------------------------------
-> mort_rate = .035

    Variable |        Obs        Mean    Std. Dev.       Min        Max
-------------+---------------------------------------------------------
    survived |      2,000      28.033    28.43426          0        215

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

写在后面

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

相关课程

连享会-直播课 上线了！
http://lianxh.duanshu.com

免费公开课：

直击面板数据模型 - 连玉君，时长：1小时40分钟
Stata 33 讲 - 连玉君, 每讲 15 分钟.
部分直播课课程资料下载 (PPT，dofiles等)

课程一览

支持回看，所有课程可以随时购买观看。

专题	嘉宾	直播/回看视频
⭐ 最新专题 ⭐		DSGE, 因果推断, 空间计量等
⭕ Stata数据清洗	游万海	直播, 2 小时，已上线
研究设计	连玉君	我的特斯拉-实证研究设计，-幻灯片-
面板模型	连玉君	动态面板模型，-幻灯片-
面板模型	连玉君	直击面板数据模型 [免费公开课，2小时]

Note: 部分课程的资料，PPT 等可以前往连享会-直播课主页查看，下载。

关于我们

Stata连享会 由中山大学连玉君老师团队创办，定期分享实证分析经验。直播间有很多视频课程，可以随时观看。
连享会-主页和知乎专栏，300+ 推文，实证分析不再抓狂。
公众号推文分类： 计量专题 | 分类推文 | 资源工具。推文分成 内生性 | 空间计量 | 时序面板 | 结果输出 | 交乘调节 五类，主流方法介绍一目了然：DID, RDD, IV, GMM, FE, Probit 等。
公众号关键词搜索/回复 功能已经上线。大家可以在公众号左下角点击键盘图标，输入简要关键词，以便快速呈现历史推文，获取工具软件和数据下载。常见关键词：课程, 直播, 视频, 客服, 模型设定, 研究设计, stata, plus, 绘图, 编程, 面板, 论文重现, 可视化, RDD, DID, PSM, 合成控制法 等