#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
Knock<- read.table("C:/Users/Challenger/Documents/Stats/TEA_flies_knockdown1.txt",
                   header = TRUE, stringsAsFactors = TRUE)
attach(Knock)

#Define variables
time<- Knockdown
event<- Status
anesthesia<-Treatment

#Checking the variables
summary(time)
summary(event)
summary(anesthesia)

#Test of normality, significant result = non-normal
shapiro.test(time)

#plotting the residuals
qqnorm(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")

#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(50,0.8, c("7510","7520","7530","7540","7550","7560"), 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=Knock,Surv(time,event == 1) ~ anesthesia)
summary(CPsurvivalF)

#test of normality for the recovery times, if significant, data is not normal
shapiro.test(time)