# Replicate = replicate ID
# Time = time of events in hours
# Event = binary 1/dead 0/censored
# Infection and Incubation = experimental conditions
# Colony series = batch number of the colony
# Colony number = unique colony identifier
# Experimentalist = experimentalist (slightly different manipulation practice)
# Group_letter 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_socind <- fread("Fig5FG_data.csv")

data_socind$Infection <- factor(data_socind$Infection)
data_socind$Incubation <- factor(data_socind$Incubation)
data_socind$Replicate <- factor(data_socind$Replicate)

#Split dataset by pathogen
data_socind %>%
  filter(Infection %in% "Pe") -> data_socind_pe 

data_socind %>%
  filter(Infection %in% "Mb",
         Time >24) -> data_socind_mb
# death before 24h with Mb are discarded as they cannot be attributed to the fungus (removal reduces the difference between curves). 


#MB analysis

null_model <- coxme ( Surv (time = Time, event = Event) ~ 1                 + ( 1 | Experimentalist) + ( 1 | Replicate) + ( 1 | ColonyNumber), data = data_socind_mb)
full_model <- coxme ( Surv (time = Time, event = Event) ~ 1 + Incubation + ( 1 | Experimentalist) + ( 1 | Replicate) + ( 1 | ColonyNumber), data = data_socind_mb)

# anova(null_model   ,  full_model )
summary(full_model)

Anova(full_model, type="II")


#Plotting

surv_fitMb <- survfit(Surv(Time, Event) ~ Group_letter, data = data_socind_mb)

ggsurvplot(surv_fitMb,
           data = data_socind_mb, conf.int = TRUE, legend = "top", 
           palette = c("#EFB565", "#CC5B38"), 
           linetype = c("dotted", "solid"), 
           legend.title="", 
           legend.labs=c("Social", "Individual"))+ 
  guides(colour = guide_legend(nrow = 1))  -> plotSocindMb

print(plotSocindMb)

plotSocindMb <- plotSocindMb$plot

plotSocindMb + 
  ggplot2::annotate("text", x = 140, y = 0.8, 
                    label = "*", size = 8) -> plotSocindMb

print(plotSocindMb)



#Pe analysis

null_model_pe <- coxme ( Surv (time = Time, event = Event) ~ 1                 + ( 1 | Experimentalist) + ( 1 | Replicate) + ( 1 | ColonyNumber), data = data_socind_pe)
full_model_pe <- coxme ( Surv (time = Time, event = Event) ~ 1 + Incubation + ( 1 | Experimentalist) + ( 1 | Replicate) + ( 1 | ColonyNumber), data = data_socind_pe)

summary(full_model_pe)

Anova(full_model_pe, type="II")



#Plotting Pe

surv_fitPe <- survfit(Surv(Time, Event) ~ Group_letter, data = data_socind_pe)

ggsurvplot(surv_fitPe,
           data = data_socind_pe, conf.int = TRUE, legend = "top", 
           palette = c("#63AF87", "#37799E"), 
           linetype = c("dotted", "solid"), 
           legend.title="", 
           legend.labs=c("Social", "Individual"))+ 
  guides(colour = guide_legend(nrow = 1)) -> plotSocindPe

print(plotSocindPe)

plotSocindPe <- plotSocindPe$plot


ggarrange(plotSocindPe, plotSocindMb, ncol = 2)


#Median survival times: 
surv_fitPe
surv_fitMb