Survival Rates and Negative Binomial

* This estimates average survival rates (days)
* Given different daily mortality rates (predation/accident etc.)

* Given a constant failure rate and a defined number of eliminations then the expected survival rate follows the negative binomial distribution.

* Negative binomial results are as follows:
* mean = pr/(1-p)
* r is number of failures:
* r = 1

* mortality rate = .01
* if p=.99 ->
di "E= " .99/.01

* mortality rate = .015
* if p=.985 ->
di "E= " .985/.015


* mortality rate = .02
* if p=.98 ->
di "E= " .98/.02

* mortality rate = .025
* if p=.975 ->
di "E= " .975/.025

* mortality rate = .03
* if p=.97 ->
di "E= " .97/.03

* mortality rate = .035
* if p=.965 ->
di "E= " .965/.035


* Let's see if the simulation gives us the expected results
clear
set obs 12000
gen mort_rate = mod(_n,6)/200+.01
* This will split the data into 6 groups of 200 with each group having mortality rates from .01 to .035

tab mort_rate

* An indicator that the mosquito is alive
gen living=1

* An indicator that the mosquito survived
gen survived=0

* Sample 500 days, see how many mosquitos survive to live for each day
forv i = 1/500 {
  * Draw a single temporary binomial draw to see if this day the mosquitos die
  qui gen died = rbinomial(1,mort_rate) if living==1
  * If the mosquitos died then change their "living" status
  qui replace living=0 if living==1 & died == 1
  * Replace survived to indicate the current day if they are still alive
  qui replace survived=`i' if living==1

  * Displays user feedback telling what round the days currently are on and what percentage of mosquitoes live.
  qui sum living
  di "Round `i' :" r(mean)

  * Drop the temporary variable "died"
  drop died
}

hist survived, by(mort_rate)

bysort mort_rate: sum survived