Input=("
day	grass_cover	oil_pad
1	16	28
2	16	28
3	16	28
4	16	28
5	16	28
6	16	28
7	16	28
8	16	27
9	16	26
10	16	25
11	16	24
12	16	24
13	16	24
14	16	24
15	16	24
16	16	23
17	16	21
18	16	20
19	15	19
20	15	17
21	15	17
22	14	14
23	14	12
24	14	11
25	14	8
26	14	8
27	12	8
28	9	NA
")
killdeer<-read.table(textConnection(Input),header=TRUE) 
par(mfrow=c(2,2))
plot(killdeer$day,killdeer$grass_cover,ylim=c(0,28),pch=16,col="darkgreen", xlab="Survival (days)",ylab="Number of nests")
points(killdeer$day,killdeer$oil_pad,pch=16,col="brown") 
plot(killdeer$day, killdeer$grass_cover/max(killdeer$grass_cover),ylim =c(0,1),pch=16,col="darkgreen",xlab="Survival (days)",ylab="Number of nests")
points(killdeer$day,killdeer$oil_pad/max(killdeer$oil_pad,na.rm=T), pch=18,col="brown")
#####
install.packages("survival")
library(survival)
grass_survive<-killdeer[nrow(killdeer),2] 
data<-c(rep(killdeer[nrow(killdeer),1],grass_survive),rep(0,grass_survive), rep(1,grass_survive))
grass_survival<-matrix(data,ncol=3,byrow=FALSE)
z<-length(levels(as.factor(killdeer[,2])))-1 
temp<-grass_survive
for(i in 1:z){
d<-max(which(killdeer[,2]>temp))
diff<-killdeer[d,2]-temp
for(j in 1:diff) grass_survival<-rbind(grass_survival,c(d,1,1))
temp<-killdeer[d,2]}
colnames(grass_survival)<-c("time","status","habitat")
oil<-Filter(Negate(anyNA),killdeer$oil_pad)
# oil
oil_survive<-oil[length(oil)]
data<-c(rep(killdeer[length(oil),1],oil_survive),rep(0,oil_survive),rep(2,oil_survive))
oil_survival<-matrix(data,ncol=3,byrow=FALSE)
temp<-oil_survive
for(i in 1:13){
d<-max(which(killdeer[,3]>temp))
diff<-killdeer[d,3]-temp
for(j in 1:diff) oil_survival<-rbind(oil_survival,c(d,1,2))
temp<-killdeer[d,3]} # end i loop
killdeerSurv<-as.data.frame(rbind(grass_survival,oil_survival))
# killdeerSurv 
kdfit<-survfit(Surv(killdeerSurv$time, killdeerSurv$status) ~ killdeerSurv$habitat)
plot(kdfit,xlab="Days",ylab="Nest survivorship",lwd=2, lty=c(1,2),col= c("darkgreen","brown"))
summary(kdfit)
kdgrass<-Surv(killdeerSurv$time[1:16], killdeerSurv$status[1:16]) ~ killdeerSurv$habitat[1:16]
kdoil<-Surv(killdeerSurv$time[17:44], killdeerSurv$status[17:44]) ~ killdeerSurv$habitat[17:44]
g95<-survfit(kdgrass)
o95<-survfit(kdoil)
plot(g95,col="darkgreen",xlab="Days",ylab="Nest survivorship", main="Kaplan-Meier estimate with 95% confidence bounds",lwd=2)
lines(o95,col="brown",lwd=2)
kdc<-coxph(Surv(killdeerSurv$time, killdeerSurv$status) ~ killdeerSurv$habitat)
summary(kdc)