#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/CO2_recovery_deaths.txt", header = TRUE, stringsAsFactors = TRUE) attach(Recovdeath) #Define variables time <- AvgRecov event <- Status anesthesia <-Exposure.sec dead <- Deaths #Checking the variables summary(time) summary(event) summary(anesthesia) 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) ~ 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: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)~anesthesia) summary(CPsurvival) #Use corrected coxph with firth CPsurvivalF<- coxphf(data=Recovdeath,Surv(time,event == 1) ~ anesthesia) summary(CPsurvivalF) #Get median recovery kmsurvivalA CPsurvival #package for multiple comparisons install.packages("multcomp") library("multcomp") #run glm for deaths dead.glm <- glm(dead~anesthesia, family = poisson) summary(dead.glm) dead.glm2 <- glm(dead~anesthesia, family = quasipoisson) summary(dead.glm2) #compare models anova(dead.glm,dead.glm2, test = "Chi") #run multiple comparisons summary(glht(dead.glm2, mcp(anesthesia="Tukey")))