# CHARGING PACKAGES
library(survminer)
library(lme4)


# CHARGING DATA
Outbreaks <- read.csv("OutbreaksA_characteristics.csv", sep=";", header=TRUE, dec=",")
Outbreaks$Loc_ID <- factor(Outbreaks$Loc_ID)
Outbreaks$Densite_pop_1000 <- as.numeric(Outbreaks$Densite_pop_1000)
Outbreaks$Time2cities <- as.numeric(Outbreaks$Time2cities)
Outbreaks$PluieMoy <- as.numeric(Outbreaks$PluieMoy)
Outbreaks$End_1 <- Outbreaks$End + 1      

Outbreaks$CATIc7_01 <- factor(Outbreaks$CATIc7_01)
Outbreaks$CATIc7_delaiC <- factor(Outbreaks$CATIc7_delaiC, levels = c("D",0,"C","B","A"))
Outbreaks$CATIc7_delaiC2 <- factor(Outbreaks$CATIc7_delaiC2, levels = c("E",0,"D","C","B","A"))
Outbreaks$CATIc7_CasC <- factor(Outbreaks$CATIc7_CasC, levels = c("D",0,"C","B","A"))
Outbreaks$CATIc7_weekC <- factor(Outbreaks$CATIc7_weekC, levels = c("D",0,"C","B","A"))

CATIc7 <- subset(Outbreaks, Outbreaks$CATIc7_01 == 1)
CATIc7$FirstRep_1 <- CATIc7$Start + CATIc7$CATIc7_delai


# 1 CHARACTERISTICS OF RESPONSE GROUPS
    
### Semester
Rep01_Semester_stat <- glmer(CATIc7_01 ~ SemesterNb + (1 | Loc_ID), data=Outbreaks, family=binomial)
Delai_Semester_stat <- coxph(Surv(Start, FirstRep_1, Outbreak_end) ~ SemesterNb + cluster(Loc_ID), data=CATIc7)
RepCas_Semester_stat <- glmer.nb(CATIc7_Cas ~ SemesterNb + (1 | Loc_ID), data=CATIc7)
RepWeek_Semester_stat <- glmer.nb(CATIc7_week ~ SemesterNb + (1 | Loc_ID), data=CATIc7)
summary(Rep01_Semester_stat)
summary(Delai_Semester_stat)
summary(RepCas_Semester_stat)
summary(RepWeek_Semester_stat)

### Population density
Rep01_Density_stat <- glmer(CATIc7_01 ~ Densite_pop_1000 + (1 | Loc_ID), data=Outbreaks, family=binomial)
Delai_Density_stat <- coxph(Surv(Start, FirstRep_1, Outbreak_end) ~ Densite_pop_1000 + cluster(Loc_ID), data=CATIc7)
RepCas_Density_stat <- glmer.nb(CATIc7_Cas ~ Densite_pop_1000 + (1 | Loc_ID), data=CATIc7)
RepWeek_Density_stat <- glmer.nb(CATIc7_week ~ Densite_pop_1000 + (1 | Loc_ID), data=CATIc7)
summary(Rep01_Density_stat)
summary(Delai_Density_stat)
summary(RepCas_Density_stat)
summary(RepWeek_Density_stat)

### Travel time to cities
Rep01_Time2cities_stat <- glmer(CATIc7_01 ~ Time2cities + (1 | Loc_ID), data=Outbreaks, family=binomial)
Delai_Time2cities_stat <- coxph(Surv(Start, FirstRep_1, Outbreak_end) ~ Time2cities + cluster(Loc_ID), data=CATIc7)
RepCas_Time2cities_stat <- glmer.nb(CATIc7_Cas ~ Time2cities + (1 | Loc_ID), data=CATIc7)
RepWeek_Time2cities_stat <- glmer.nb(CATIc7_week ~ Time2cities + (1 | Loc_ID), data=CATIc7)
summary(Rep01_Time2cities_stat)
summary(Delai_Time2cities_stat)
summary(RepCas_Time2cities_stat)
summary(RepWeek_Time2cities_stat)

### Accumulated incidence
Rep01_Inc_stat <- glmer(CATIc7_01 ~ IncR_2010_14 + (1 | Loc_ID), data=Outbreaks, family=binomial)
Rep1234_Inc_stat <- coxph(Surv(Start, FirstRep_1, Outbreak_end) ~ IncR_2010_14 + cluster(Loc_ID), data=CATIc7)
RepCas_Inc_stat <- glmer.nb(CATIc7_Cas ~ IncR_2010_14 + (1 | Loc_ID), data=CATIc7)
RepWeek_Inc_stat <- glmer.nb(CATIc7_week ~ IncR_2010_14 + (1 | Loc_ID), data=CATIc7)
summary(Rep01_Inc_stat)
summary(Rep1234_Inc_stat)
summary(RepCas_Inc_stat)
summary(RepWeek_Inc_stat)

### OCV coverage
Rep01_CV_stat <- glmer(CATIc7_01 ~ CV_2010_14 + (1 | Loc_ID), data=Outbreaks, family=binomial)
Rep1234_CV_stat <- coxph(Surv(Start, FirstRep_1, Outbreak_end) ~ CV_2010_14 + cluster(Loc_ID), data=CATIc7)
RepCas_CV_stat <- glmer.nb(CATIc7_Cas ~ CV_2010_14 + (1 | Loc_ID), data=CATIc7)
RepWeek_CV_stat <- glmer.nb(CATIc7_week ~ Previous_OCV + (1 | Loc_ID), data=CATIc7)
summary(Rep01_CV_stat)
summary(Rep1234_CV_stat)
summary(RepCas_CV_stat)
summary(RepWeek_CV_stat)

### Previous cases in the same locality during the study
Rep01_PrevInc_stat <- glmer(CATIc7_01 ~ Previous_Inc + (1 | Localite_ID), data=Outbreaks, family=binomial)
Rep1234_PrevInc_stat <- coxph(Surv(Start, FirstRep_1, Outbreak_end) ~ Previous_Inc + cluster(Localite_ID), data=CATIc7)
RepCas_PrevInc_stat <- glmer.nb(CATIc7_Cas ~ Previous_Inc + (1 | Localite_ID), data=CATIc7)
RepWeek_PrevInc_stat <- glmer.nb(CATIc7_week ~ Previous_Inc + (1 | Localite_ID), data=CATIc7)
summary(Rep01_PrevInc_stat)
summary(Rep1234_PrevInc_stat)
summary(RepCas_PrevInc_stat)
summary(RepWeek_PrevInc_stat)

### Previous CATIs in the same locality during the study
Rep01_PrevCATI_stat <- glmer(CATIc7_01 ~ PreviousCATIyear + (1 | Localite_ID), data=Outbreaks, family=binomial)
Rep1234_PrevCATI_stat <- coxph(Surv(Start, FirstRep_1, Outbreak_end) ~ PreviousCATIyear + cluster(Localite_ID), data=CATIc7)
RepCas_PrevCATI_stat <- glmer.nb(CATIc7_Cas ~ PreviousCATIyear + (1 | Localite_ID), data=CATIc7)
RepWeek_PrevCATI_stat <- glmer.nb(CATIc7_week ~ PreviousCATIyear + (1 | Localite_ID), data=CATIc7)
summary(Rep01_PrevCATI_stat)
summary(Rep1234_PrevCATI_stat)
summary(RepCas_PrevCATI_stat)
summary(RepWeek_PrevCATI_stat)

### Mean daily rainfall
Rep01_Rain_stat <- glmer(CATIc7_01 ~ PluieMoy + (1 | Loc_ID), data=Outbreaks, family=binomial)
Rep1234_Rain_stat <- coxph(Surv(Start, FirstRep_1, Outbreak_end) ~ PluieMoy + cluster(Loc_ID), data=CATIc7)
RepCas_Rain_stat <- glmer.nb(CATIc7_Cas ~ PluieMoy + (1 | Loc_ID), data=CATIc7)
RepWeek_Rain_stat <- glmer.nb(CATIc7_week ~ PluieMoy + (1 | Loc_ID), data=CATIc7)
summary(Rep01_Rain_stat)
summary(Rep1234_Rain_stat)
summary(RepCas_Rain_stat)
summary(RepWeek_Rain_stat)

### No. of cases during the first 3 days of outbreak
Rep01_CasJ3_stat <- glmer(CATIc7_01 ~ Cas_J3 + (1 | Loc_ID), data=Outbreaks, family=binomial)
Rep1234_CasJ3_stat <- coxph(Surv(Start, FirstRep_1, Outbreak_end) ~ Cas_J3 + cluster(Loc_ID), data=CATIc7)
RepCas_CasJ3_stat <- glmer.nb(CATIc7_Cas ~ Cas_J3 + (1 | Loc_ID), data=CATIc7)
RepWeek_CasJ3_stat <- glmer.nb(CATIc7_week ~ Cas_J3 + (1 | Loc_ID), data=CATIc7)
summary(Rep01_CasJ3_stat)
summary(Rep1234_CasJ3_stat)
summary(RepCas_CasJ3_stat)
summary(RepWeek_CasJ3_stat)

### No. of positive cultures during the first 3 days
Rep01_CultJ3_stat <- glmer(CATIc7_01 ~ Cultures_pos_J3 + (1 | Loc_ID), data=Outbreaks, family=binomial)
Rep1234_CultJ3_stat <- coxph(Surv(Start, FirstRep_1, Outbreak_end) ~ Cultures_pos_J3 + cluster(Loc_ID), data=CATIc7)
RepCas_CultJ3_stat <- glmer.nb(CATIc7_Cas ~ Cultures_pos_J3 + (1 | Loc_ID), data=CATIc7)
RepWeek_CultJ3_stat <- glmer.nb(CATIc7_week ~ Cultures_pos_J3 + (1 | Loc_ID), data=CATIc7)
summary(Rep01_CultJ3_stat)
summary(Rep1234_CultJ3_stat)
summary(RepCas_CultJ3_stat)
summary(RepWeek_CultJ3_stat)


# 2 COMPARISONS BETWEEN OUTBREAKS THAT RECEIVED A RESPONSE

## 2.1 Number of cases from the 4th day of outbreak according to the response promptness

### 2.1.1 Univariate graphs
Cas4_Rep1234_plot <- ggplot(data=CATIc7, aes(x=CATIc7$CATIc7_delaiC, y=CATIc7$Cas_J4_)) +   
  theme(legend.position="none") + 
  theme_classic() +
  scale_x_discrete(name="Time to the first complete CATI", limits=c("D","C", "B", "A"), labels=c(">7 days", "3 to 7 days", "2 days", "<=1 day")) +
  scale_y_continuous(name="Number of cases from the 4th day of outbreak") +  
  geom_jitter(alpha=0.5) +
  geom_boxplot(colour=c("red", "orange", "yellow", "green"), alpha=0, outlier.shape=NA, lwd=1)

### 2.1.2 Univariate analysis with a generalized linear mixed model
Cas4_delai_lmer <- glmer.nb(Cas_J4_ ~  CATIc7_delaiC + (1 | Loc_ID) , data = CATIc7)  
summary(Cas4_delai_lmer)

### 2.1.3 Multivariate analysis with a generalized linear mixed model
Cas4_delai_glmer <- glmer.nb(Cas_J4_ ~  CATIc7_delaiC + Cultures_pos_J3 + Densite_pop_1000 + IncR_2010_14 + CV_2010_14 + SemesterNb + (1 | Localite_ID) , data = CATIc7)
summary(Cas4_delai_glmer)

## 2.2 Outbreak duration according to the response promptness

### 2.2.1 Univariate graphs
CATIc7$CATIc7_delaiC <- factor(CATIc7$CATIc7_delaiC, levels=c("D","C","B","A"))
DR_surv_delai <- survfit(Surv(Duration, Outbreak_end) ~ CATIc7_delaiC + cluster(Loc_ID), data=CATIc7)
Surv_delai_plot <- ggsurvplot(DR_surv_delai,
                              fun="event", conf.int=FALSE, pval = FALSE, risk.table=FALSE, surv.median.line = "none",     
                              xlim=c(0, 100), break.x.by=10, xlab="Duration of outbreaks (days)", ylab="Cumulative probability of outbreak end",
                              palette = c("red", "orange", "yellow", "green"),
                              legend.labs=c(">7 days","3 to 7 days","2 days","≤1 day"), legend = c(0.7, 0.3), legend.title="Time to first correct response")    

### 2.2.2 Univariate survival analysis using a Cox model for repeated events with Anderson-Gills correction (AGCP)
DR_Cox_Uni_delai <- coxph(Surv(Start, End_1, Outbreak_end) ~ CATIc7_delaiC + cluster(Loc_ID), data=CATIc7)  
summary(DR_Cox_Uni_delai)

### 2.2.3 Multivariate survival analysis using a Cox model for repeated events with Anderson-Gills correction (AGCP)
DR_Cox_Multi_delai <- coxph(Surv(Start, End_1, Outbreak_end) ~ CATIc7_delaiC + Cas_J3  + Cultures_pos_J3 + PluieMoy + Densite_pop_1000 + Time2cities + IncR_2010_14 + CV_2010_14 + SemesterNb + cluster(Localite_ID), data=CATIc7)  
summary(DR_Cox_Multi_delai)

## 2.3 Number of cases from the 4th day of outbreak according to the response intensity

### 2.3.1 Univariate graphs
Cas4_RepWeek_plot <- ggplot(data=CATIc7, aes(x=CATIc7$CATIc7_weekC, y=CATIc7$Cas_J4_)) +   
  theme(legend.position="none") +
  theme_classic() +
  scale_x_discrete(name="Response/week ratio", limits=c("D","C","B","A"), labels=c("<0.25","0.25 to 0.5", "0.5 to 1", ">=1")) +
  scale_y_continuous(name="Number of cases from the 4th day of outbreak") +  
  geom_jitter(alpha=0.5) +
  geom_boxplot(colour=c("red", "orange", "yellow", "green"), alpha=0, outlier.shape=NA, lwd=1)

### 2.3.3 Univariate analysis with a generalized linear mixed model
Cas4_CATIWeek_lmer <- glmer.nb(Cas_J4_ ~ CATIc7_weekC + (1 | Loc_ID) , data = CATIc7)  
summary(Cas4_CATIWeek_lmer)

### 2.1.4 Multivariate analysis with a generalized linear mixed model
Cas4_CATIWeek_glmer <- glmer.nb(Cas_J4_ ~ CATIc7_delaiC + Cultures_pos_J3 + Densite_pop_1000 + IncR_2010_14 + CV_2010_14 + SemesterNb + (1 | Localite_ID) , data = CATIc7)
summary(Cas4_CATIWeek_glmer)

## 2.4 Outbreak duration according to the response intensity

### 2.4.1 Univariate graphs 
DR_surv_CATICas <- survfit(Surv(Duration, Outbreak_end) ~ CATIc7_CasC + cluster(Loc_ID), data=CATIc7)
Surv_RepCas_plot <- ggsurvplot(DR_surv_RepCas,
fun="event", conf.int=FALSE, pval = FALSE, risk.table=FALSE, surv.median.line = "none",    
xlim=c(0, 100), break.x.by=10, xlab="Duration of outbreaks (days)", ylab="Cumulative probability of outbreak end", 
palette = c("red", "orange", "yellow", "green"),
legend.labs=c("<0.25","0.25 to 0.5","0.5 to 1",">=1"), legend = c(0.7, 0.3), legend.title="Response/case ratio")

### 2.4.2 Univariate survival analysis using a Cox model for repeated events with Anderson-Gills correction (AGCP)
DR_Cox_Uni_CATICas <- coxph(Surv(Start, End_1, Outbreak_end) ~ CATIc7_CasC + cluster(Loc_ID), data=CATIc7)
summary(DR_Cox_Uni_CATICas)

### 2.4.3 Multivariate survival analysis using a Cox model for repeated events with Anderson-Gills correction (AGCP) +++   
DR_Cox_Multi_CATICas <- glmer.nb(Cas_J4_ ~  CATIc7_weekC + IncR_2010_14 + PreviousCATIyear + SemesterNb + (1 | Localite_ID) , data = CATIc7)  
summary(DR_Cox_Multi_CATICas)


