# Day = day of experiment
# Time = time of events in jours
# Event = binary 1/dead 0/censored
# Infection and Injury = experimental conditions
# Colony series = batch number of the colony
# Colony number = unique colony identifier
# Experimentalist = experimentalist (slightly different manipulation practice)
# Group_full and Group_letters are dummy variables used in the script to identify experimental conditions. 

library("data.table")
library("tidyverse")
library("ggplot2")
library("survival")
library("coxme")
library("multcomp")
library("car")
library("survminer")
library("ggplotify")

#import data
data_infection_model <- fread("Fig4C_data.csv")

data_infection_model$Injury <- factor(data_infection_model$Injury, levels = c("None", "CutAntenna"))


#fitsurv_infection_model <- survfit(Surv(Time, Event) ~ Infection+Injury, data = data_infection_model)

data_infection_model$Infection <- factor(data_infection_model$Infection)
data_infection_model$Injury <- factor(data_infection_model$Injury)

null_model <- coxme ( Surv (time = Time, event = Event) ~ 1                 + ( 1 | Experimentalist) + (1 | ColonySeries)+ (1 | ColonyNumber), data = data_infection_model)
full_model <- coxme ( Surv (time = Time, event = Event) ~  1 + Infection + Injury + ( 1 | Experimentalist) + (1 | ColonySeries) + (1 | ColonyNumber), data = data_infection_model)
interaction_model_n <- coxme ( Surv (time = Time, event = Event) ~ 1 + Infection + Injury + Infection:Injury + ( 1 | Experimentalist) + (1 | ColonySeries) + (1 | ColonyNumber), data = data_infection_model)

#Test the interaction
anova(full_model, interaction_model_n)
summary(interaction_model_n)

# Anova on interaction model 
Anova(interaction_model_n, type="III")

# contrast matrix for multiple comparisons
contrast_matrix = rbind(
  "Healthy_cut - Healthy_uncut" =c(0,	1,	0),
  "Infected_uncut - Healthy_uncut" =c(	1,	0,	0),
  "Infected_cut - Healthy_uncut" =c(	1,	1,	1),
  "Infected_uncut-Healthy_cut" =c(	1,	-1,	0),
  "Infected_cut - Healthy_cut" =c(	1,	0,	1),
  "Infected_cut - Infected_uncut" =c(	0,	1,	1))
summary(glht(interaction_model_n,linfct=contrast_matrix),test=adjusted("BH"))


#Plotting

surv_fit <- survfit(Surv(Time, Event) ~ Group_letter, data = data_infection_model)

ggsurvplot(surv_fit,
           data = data_infection_model, conf.int = TRUE, legend = "top", 
           palette = c("#63AF87", "#37799E", "#EFB565", "#CC5B38"), 
           linetype = c("dotted", "solid","dotted", "solid"), 
           legend.title="", 
           legend.labs=c("Intact - Uninfected", "Intact - Infected", "Injured - Uninfected", "Injured - Infected")) + 
  guides(colour = guide_legend(nrow = 2)) -> plotInfModel

print(plotInfModel)

plotInfModel <- plotInfModel$plot