#egg survival library(MASS) library(survival) library(survminer) library(survMisc) #input file egg<-read.csv("C:/Users/USER/Documents/Aerodramus/HY manuscript/Data/Nesting/Survival-HZ/eggsurvival.csv", header=TRUE) head(egg) summary(egg) #calculations fit1 <- survfit(Surv(days, dead) ~ species, data=egg,subset=(egg$species=="b")) fit2 <- survfit(Surv(days, dead) ~ species, data=egg,subset=(egg$species=="w")) #time range for calculation is 1 less than actual range to avoid NAs in last row sf1<-summary(fit1, times=0:41) sf1 sf2<-summary(fit2, times=0:48) sf2 cox1 <- coxph(Surv(days, dead) ~ species,data=egg) cox1 rsq(cox1) survdiff(Surv(days, dead) ~ species,data=egg) summary(cox1) res.zph1 <- cox.zph(cox1) #plot(res.zph1) print(res.zph1) pairwise_survdiff(Surv(days, dead) ~ species,data=egg) #plotting par(mar=c(5,4,2,1)) svg("C:/Users/USER/Documents/Aerodramus/HY manuscript/Data/Nesting/Survival-HZ/eggsurvival.svg") plot(fit1, xlab="Days",ylab="Egg Survivability",col="#1b1b1b",conf.int=F, lwd=2,cex.axis=0.8,cex.lab=1.2,xlim=c(0,26),ylim=c(0,1)) lines(fit2,col="#e69f00",conf.int=F,lwd=2,lty=1) polygon(c(sf1$time,rev(sf1$time)),c(sf1$lower,rev(sf1$upper)), col=adjustcolor("#1b1b1b",0.2),border = F) polygon(c(sf2$time,rev(sf2$time)),c(sf2$lower,rev(sf2$upper)), col=adjustcolor("#e69f00",0.2),border = F) #change R-square value based on 'rsq' output above text(4,0.15,bquote( R^2 == .(rs), list(rs=round(0.076,2)))) text(4,0.1,"p-value < 0.001") # text(4,0.05,"n = 933") dev.off()