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

#Define variables
time <- AvgRecov
event <- Status
temperature <-Treatment.C
exposure <- Treatment.hrs
dead <- Deaths

#Checking the variables
summary(time)
summary(event)
summary(temperature)
summary(exposure)
summary(dead)

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

#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) ~ temperature)
summary(kmsurvivalA)
plot(kmsurvivalA, xlab ="Time", ylab ="Survival Probability")

#Kaplan-Meier non-parametric analysis by Anesthesia
kmsurvivalB<- survfit(Surv(time,event == 1) ~ temperature+exposure)
summary(kmsurvivalB)
plot(kmsurvivalB, xlab ="Time", ylab ="Survival Probability")

anova(kmsurvivalA,kmsurvivalB)

#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:7), lty=(1:7), lwd =4, xlab ="Time(mins)", ylab ="Anesthetised probability")
legend(50,0.8, c("5mins","60mins","600mins","900mins","1200mins","1500mins","1800mins"), lwd=3, lty=(1:7),col=(1:7), 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)~temperature)
summary(CPsurvival)

CPsurvival2 <- coxph(Surv(time,event == 1)~temperature+exposure)
summary(CPsurvival2)

anova(CPsurvival,CPsurvival2)

#Use corrected coxph with firth
CPsurvivalF<- coxphf(data=Recovdeath,Surv(time,event == 1) ~ temperature/exposure)
summary(CPsurvivalF)

#package for multiple comparisons
install.packages("multcomp")
library("multcomp")

#run glm for deaths
dead.glm <- glm(dead~temperature/exposure, family = poisson)
summary(dead.glm)

#compare models
anova(dead.glm,dead.glm2, test = "Chi")

#run multiple comparisons
summary(glht(dead.glm2, mcp(temperature/exposure="Tukey")))