# Load packages
library(tidyverse)
library(survival)
library(survminer)
library(lmerTest)
library(paletteer)
library(cowplot)

#the line below removes everything from the R environment.
#This ensures there is no other data etc that could interfere with the analysis. 
rm(list=ls())

# Survival dataset, analysis and plotting
survival<-read.csv("survival_analysis.csv", 
                   header = T, 
                   stringsAsFactors = T)

#data exploration
survival
str(survival)
levels(as.factor(survival$treatment))
## Reordering survival$treatment
survival$treatment <- factor(survival$treatment,
                             levels = c(
                               "Control", "Deltamethrin", "S10", "S26", "S31", "Deltamethrin_S10",
                               "Deltamethrin_S26", "Deltamethrin_S31", "S10_Deltamethrin", "S26_Deltamethrin",
                               "S31_Deltamethrin", "Deltamethrin+S10", "Deltamethrin+S26", "Deltamethrin+S31"))
survival$treatment
unique(survival$treatment)

## Subseting the dataset for each level of treatment
survival1 <- filter(survival, 
                    treatment %in% c("Control", 
                                     "Deltamethrin"))

plotfit1<- survfit(Surv(day, death) ~ treatment,data=survival1)

ggsurv1 <- ggsurvplot(plotfit1, conf.int = F, palette='treatment',
                      legend.title = "Treatment",
                      surv.scale = "percent",
                      xlab="Days post exposure", size=0.8)
p1 <- ggsurv1$plot +theme_bw()+scale_color_paletteer_d("ggthemes::stata_s1color")+ 
  theme(strip.background = element_rect(fill="white"),
        strip.text = element_text(size=13),
        axis.text = element_text(size = 10),
        axis.title = element_text(size = 13),
        legend.text = element_text(size = 10),
        legend.title = element_text(size = 12))
p1


survival2 <- filter(survival,
                    treatment %in%
                      c("Control", 
                        "S10",
                        "S26",
                        "S31"))
plotfit2<- survfit(Surv(day, death) ~ treatment,data=survival2)
ggsurv2 <- ggsurvplot(plotfit2, conf.int = F, palette='treatment',
                      legend.title = "Treatment",
                      surv.scale = "percent",
                      xlab="Days post exposure", size=0.8)
p2 <- ggsurv2$plot +theme_bw()+scale_color_paletteer_d("ggthemes::stata_s1color")+ 
  theme(strip.background = element_rect(fill="white"),
        strip.text = element_text(size=13),
        axis.text = element_text(size = 10),
        axis.title = element_text(size = 13),
        legend.text = element_text(size = 10),
        legend.title = element_text(size = 12))
p2


survival3 <- filter(survival, 
                    treatment %in%
                      c("Control",
                        "Deltamethrin_S10",
                        "Deltamethrin_S26",
                        "Deltamethrin_S31"))
plotfit3<- survfit(Surv(day, death) ~ treatment,data=survival3)
ggsurv3 <- ggsurvplot(plotfit3, conf.int = F, palette='treatment',
                      legend.title = "Treatment",
                      surv.scale = "percent",
                      xlab="Days post exposure", size=0.8)
p3 <- ggsurv3$plot +theme_bw()+scale_color_paletteer_d("ggthemes::stata_s1color")+ 
  theme(strip.background = element_rect(fill="white"),
        strip.text = element_text(size=13),
        axis.text = element_text(size = 10),
        axis.title = element_text(size = 13),
        legend.text = element_text(size = 10),
        legend.title = element_text(size = 12))
p3


survival4 <- filter(survival,
                    treatment %in%
                      c("Control", 
                        "S10_Deltamethrin",
                        "S26_Deltamethrin",
                        "S31_Deltamethrin"))
plotfit4<- survfit(Surv(day, death) ~ treatment,data=survival4)
ggsurv4 <- ggsurvplot(plotfit4, conf.int = F, palette='treatment',
                      legend.title = "Treatment",
                      surv.scale = "percent",
                      xlab="Days post exposure", size=0.8)
p4 <- ggsurv4$plot +theme_bw()+scale_color_paletteer_d("ggthemes::stata_s1color")+ 
  theme(strip.background = element_rect(fill="white"),
        strip.text = element_text(size=13),
        axis.text = element_text(size = 10),
        axis.title = element_text(size = 13),
        legend.text = element_text(size = 10),
        legend.title = element_text(size = 12))
p4

survival5 <- filter(survival, 
                    treatment %in%
                      c("Control",
                        "Deltamethrin+S10",
                        "Deltamethrin+S26",
                        "Deltamethrin+S31"))
plotfit5<- survfit(Surv(day, death) ~ treatment,data=survival5)
ggsurv5 <- ggsurvplot(plotfit5, conf.int = F, palette='treatment',
                      legend.title = "Treatment",
                      surv.scale = "percent",
                      xlab="Days post exposure", size=0.8)
p5 <- ggsurv5$plot +theme_bw()+scale_color_paletteer_d("ggthemes::stata_s1color")+ 
  theme(strip.background = element_rect(fill="white"),
        strip.text = element_text(size=13),
        axis.text = element_text(size = 10),
        axis.title = element_text(size = 13),
        legend.text = element_text(size = 10),
        legend.title = element_text(size = 12))
p5

## For plot grouping
plot_grid(p2, p3, p4, p5, labels = c("A", "B", "C", "D"), align = "hv")

#Test effect of treatment on mosquitoes mortality
Model.1<- coxph(Surv(day, death) ~ treatment +  
                  frailty(replicate,
                          dist='gamma'), 
                data=survival)

Model.2 <- coxph(Surv(day, death) ~ frailty(replicate, dist='gamma'),
                 data=survival)
anova(Model.2, Model.1)
### treatment is significant Chisq = 386.11 Df= 13 p< 2.2e-16
summary(Model.1)