#Install package for Keplar-Meier curve analysis install.packages("survival") #Load tools from survival package library(survival) #load and install coxphf "cox proportional hazard model with firth correction" install.packages("coxphf") library(coxphf) #Load table into a new object "Awake" and attach it to the workspace Recov<- read.table("C:/Users/Challenger/Documents/Stats/TEA_flies_TEA25_recovery.txt", header = TRUE, stringsAsFactors = TRUE) attach(Recov) #Define variables time<- Recovery event<- Status sex<- Sex anesthesia<-Treatment #Checking the variables summary(time) summary(event) summary(sex) summary(anesthesia) #Kaplan-Meier non-parametric analysis kmsurvival<- survfit(Surv(time,event == 1)~1) summary(kmsurvival) plot(kmsurvival, xlab ="Time", ylab ="Survival Probability") #Kaplan-Meier non-parametric analysis by Anesthesia kmsurvivalA<- survfit(Surv(time,event == 1) ~ anesthesia) summary(kmsurvivalA) plot(kmsurvivalA, xlab ="Time", ylab ="Survival Probability") #Kaplan-Meier non-parametric analysis by Anesthesia and sex kmsurvivalB<- survfit(Surv(time,event == 1) ~ anesthesia+sex) summary(kmsurvivalB) plot(kmsurvivalA, xlab ="Time", ylab ="Survival Probability") #anesthesia and sex plot with different colors and legend; lwd = line weight and lty = line pattern ***make sure to change titles and axes plot(kmsurvivalA, col=c(1:6), lty=(1:6), lwd =4, xlab ="Time(mins)", ylab ="Anesthetised probability") legend(180,1, c("5010","5020","5030","5040","5050","5060"), lwd=3, lty=(1:6),col=(1:6), cex = 1) #anesthesia and sex plot with different colors and legend; lwd = line weight and lty = line pattern ***make sure to change titles and axes plot(kmsurvivalB, col=c(1:8), lty=(1:8), lwd =4, xlab ="Time(mins)", ylab ="Anesthetised probability") legend(200,1, c("Chilling (F)","Chilling(M)","CO2(F)","CO2(M)","Control(F)","Control(M)","TEA(F)","TEA(M)"), lwd=3, lty=(1:8),col=(1:8), cex = 1) #Modifying labels, title, and margins margin default is c(5,4,4,2) +0.1 (bottom, left, top, right) par(cex.lab = 2, cex.main = 2.5, cex.axis = 1.5, mar = c(5,5,4,2) + 0.1) #Cox proportional hazard model -coefficients and hazards rates CPsurvival <- coxph(Surv(time,event == 1) ~ anesthesia+sex, method = "breslow") summary(CPsurvival) CPsurvival2 <- coxph(Surv(time,event == 1)~anesthesia) summary(CPsurvival2) anova(CPsurvival,CPsurvival2) #Use corrected coxph with firth CPsurvivalF<- coxphf(data=Recov,Surv(time,event == 1) ~ anesthesia) summary(CPsurvivalF) #test of normality for the recovery times, if significant, data is not normal shapiro.test(time)