rm(list=ls())
setwd("~/Dropbox/GrayLab/BiosecurityExp")

library("survival")

pdf(file="CohousingSurvCurveAllFinal.pdf")

########## Process Data
# Load in raw data entered by Jenny into excel.
COHOUSING.DATA.RAW <- read.csv("~/Dropbox/GrayLab/BiosecurityExp/CohousingData.csv", 
                               nrow=355, 
                               header=TRUE)

# Divide 'treatment' variable into prevalence and time treatments, format accordingly, and create new data frame with only necessary columns (prevalence, time, and number of days post exposure animal died).
split.treatment <- strsplit(x=as.vector(COHOUSING.DATA.RAW$Treatment), split=" ")
for(i in 346:355){
  split.treatment[[i]] <- c(split.treatment[[i]][1], paste(split.treatment[[i]][3], split.treatment[[i]][4], sep=""))
}
Prevalence <- unlist(lapply(split.treatment, "[[", 1))
Time <- unlist(lapply(split.treatment, "[[", 2))
COHOUSING.DATA <- data.frame(Prevalence, Time, Day_Die=COHOUSING.DATA.RAW$Day_Die)
rm(Prevalence); rm(Time)

# Create a status variable indicating whether or not the animal actually dies on that day [0,14] or if it is to be censored [0,14+] and append to data set.
Status <- ifelse(COHOUSING.DATA$Day_Die==14, 0, 1)
COHOUSING.DATA <- data.frame(COHOUSING.DATA, Status)

# Since no mortality was observed in the control groups (and comparison with these groups are not meaningful), remove them entirely from the data set.
COHOUSING.DATA <- COHOUSING.DATA[COHOUSING.DATA$Prevalence!="Control",]

# Make sure R knows that Prevalence and Time are factors, remove control level in prevalence, and set reference levels. 
COHOUSING.DATA$Prevalence <- factor(COHOUSING.DATA$Prevalence, 
                                    levels=c("10%", "20%", "40%"))
COHOUSING.DATA$Time <- factor(COHOUSING.DATA$Time, 
                              levels=c("15min", "30min", "60min"))

########## Q1) Is there a difference among treatments?
Surv.Obj <- with(COHOUSING.DATA, Surv(Day_Die, Status))
Cox.Model.Full <- coxph(Surv.Obj ~ Prevalence + Time,
                        data=COHOUSING.DATA)
summary(Cox.Model.Full)

(Log.Rank.Full <- survdiff(Surv.Obj ~ Prevalence + Time, 
                           data=COHOUSING.DATA))

(Curve.Full <- survfit(Surv.Obj ~ Prevalence + Time,
                       data=COHOUSING.DATA))
plot(Curve.Full, 
     ylab="Survival Probability", 
     xlab="Time (Days)", 
     main="Cohousing Experiment; Survival by Prevalence x Time", 
     lty=rep(c(1,5,3),3), 
     lwd=2.5, 
     col=c(rep("gray83", 3), rep("gray50", 3), rep("black", 3)))
legend(x=0.5, 
       y=0.35, 
       legend=c("10%, 15min", "10%, 30min", "10%, 60min", "20%, 15min", "20%, 30min", "20%, 60min", "40%, 15min", "40%, 30min", "40%, 60min"), 
       lty=rep(c(1,5,3),3), 
       lwd=2.5, 
       col=c(rep("gray83", 3), rep("gray55", 3), rep("black", 3)), 
       cex=0.75, 
       title="Treatment")

########## Q2) Is there a difference among time treatments?
Cox.Model.Time <- coxph(Surv.Obj ~ Time, 
                        data=COHOUSING.DATA)
summary(Cox.Model.Time)

(Curve.Time <- survfit(Surv.Obj ~ Time,
                       data=COHOUSING.DATA))
plot(Curve.Time, 
     ylab="Survival Probability", 
     xlab="Time (Days)", 
     main="Cohousing Experiment; Survival by Time",
     lwd=2.5, 
     col=c("gray78", "gray40", "black"))
legend(x=0.5, 
       y=0.19, 
       legend=c("15min", "30min", "60min"),
       lwd=2.5, 
       col=c("gray70", "gray40", "black"), 
       cex=0.9,
       title="Time")

########## Q3) Is there a difference among prevalence treatments?
Cox.Model.Prevalence <- coxph(Surv.Obj ~ Prevalence, 
                        data=COHOUSING.DATA)
summary(Cox.Model.Prevalence)

(Curve.Prevalence <- survfit(Surv.Obj ~ Prevalence,
                             data=COHOUSING.DATA))
plot(Curve.Prevalence, 
     ylab="Survival Probability", 
     xlab="Time (Days)", 
     main="Cohousing Experiment; Survival by Prevalence",
     lwd=2.5, 
     col=c("gray78", "gray40", "black"))
legend(x=0.5, 
       y=0.19, 
       legend=c("10%", "20%", "40%"),
       lwd=2.5, 
       col=c("gray78", "gray40", "black"), 
       cex=0.9,
       title="Prevalence")

dev.off()