Showing posts with label survival analysis. Show all posts
Showing posts with label survival analysis. Show all posts

Monday, June 17, 2013

Zombie Apocalypse Survival Test - R-Powered (using Concerto)

This test is the first attempt to seriously assess the ability of individuals to survive a zombie apocalypse.  This test is administered using the R powered open-source testing platform Concerto  developed at the University of Cambridge.

The test administers a number of questions which are expertly chosen in order to evaluate the likelihood of an individual surviving a zombie apocalypse.  The R code which powers the test be found on can be found on github.  The test makes extensive use of the RConcerto supplementary package that I am working.

The purpose of the package is to facilitate test creation using Concerto.  Many of the functions either create wrappers for concerto functions or build html objects such as the function mk.facebook(link) which takes a link and turns it into a facebook share button.

In order to implement the test, I run the code straight from a permalink created by dropbox which I demonstrate the use of in a previous post in the function dropbox.eval() which is also contained within the RConcerto package.

Please try out this test at:

Zombie Survival Test

Wednesday, July 11, 2012

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