# The following code analyzes experimental data on extinction time in experimental Daphnia populations and reproduces figures
# from Drake, J.M. The distribution of extinction time in experimental populations. Ecology. See original paper for details.

library(survival)												#load library of functions for survival analysis
data<-read.table('data.csv', header=T, sep=","); attach(data)	#load data
data2<-Surv(dur,ext)											#create a survival object using censoring variable EXT and duration variable DUR
 

#Cox proportional hazards regression
out<-coxph(data2~treatment2, robust=F)							#perform first Cox proportional hazard regression
summary(out)

out2<-coxph(data2~treatment3, robust=F)							#perform Cox proportional hazards regression with low and high variation treatments pooled
summary(out2)

new.data<-data.frame(treatment3=c("A","C"))						#create data frame to plot estimated model
win.graph()

#Plot the Kaplan-Meier estimated survival function
plot(survfit(out2, newdata=new.data, type="kap"), conf.int=T, lty=c(1,2), lwd=2.5,  xlab="Time", ylab="estimated survival function")

#Test the proportional hazards assumption
test1<-cox.zph(out2)
win.graph()
plot(test1)
abline(h=0.5,lwd=3)

#Take differences in estimated survival function to obtain estimated propbability density function for extinction time and plot
pts<-survfit(out2, newdata=new.data)
t<-pts$time
t3<-t[2:43]
t2<-diff(t)
s<-diff(1-(pts$surv))
p<-s/t2
win.graph()
plot(t3,p[,1],main='Low and Medium Variability')
win.graph()
plot(t3,p[,2],main='High Variability')


#Create semilog plot to look for exponentiality in the tail, 
log.p.low<-log(p[,1])
win.graph()
plot(t3,log.p.low,main='Low and Medium Variability')
lines(lowess(t3,log.p.low))								#plot lowess scatterplot smoother

log.p.hi<-log(p[,2])
win.graph()
plot(t3,log.p.hi,main='High Variability')
lines(lowess(t3,log.p.hi))								#plot lowess scatterplot smoother