#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/Final_anesthesia_awake_withcontrols.txt", header = TRUE, stringsAsFactors = TRUE) attach(Recov) #Define variables time <- Recovery event <- Status anesthesia <- Treatment sex <- Sex #Checking the variables summary(time) summary(event) summary(anesthesia) summary(sex) #test of normality for the recovery times, if significant, data is not normal shapiro.test(time) #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") #anesthesia 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:4), lty=(1:4), lwd =4, xlab ="Time(mins)", ylab ="Anesthetised probability") legend(210,1, c("CA","Chilling","CO2","TEA"), lwd=3, lty=(1:4),col=(1:4), cex = 2) #Kaplan-Meier non-parametric analysis by Anesthesia and sex kmsurvivalB<- survfit(Surv(time,event == 1) ~ anesthesia+sex) summary(kmsurvivalB) plot(kmsurvivalB, xlab ="Time", ylab ="Survival Probability") #anesthesia (with 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 ="Probability of a fly being anesthetised") legend(215,1, c("CA(female)","CA(male)","Chilling(female)","Chilling(male)","CO2(female)","CO2(male)","TEA(female)","TEA(male)"), lwd=3, lty=(1:8),col=(1:8), cex = 1.5) #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) summary(CPsurvival) CPsurvival2 <- coxph(Surv(time,event == 1)~anesthesia+sex) summary(CPsurvival2) anova(CPsurvival,CPsurvival2) #Use corrected coxph with firth CPsurvivalF<- coxphf(data=Recov,Surv(time,event == 1) ~ anesthesia) summary(CPsurvivalF) #Get median recovery kmsurvivalB