diff --git a/Codes/Function_Plot_Conditional_Incidence_Quantiles.R b/Codes/Function_Plot_Conditional_Incidence_Quantiles.R
new file mode 100644
index 0000000000000000000000000000000000000000..d94b97c4816b9183494666b8c6b460b079d7fc06
--- /dev/null
+++ b/Codes/Function_Plot_Conditional_Incidence_Quantiles.R
@@ -0,0 +1,496 @@
+# The code does look at how the risk of Campylobacter in humans depends on environmental variables
+
+# The code does look at how the risk of Campylobacter in humans depends on environmental variables
+# Analysis was done following an quantiles division of the range of the environmental variables, independently of the number of observations
+# This is used for the reconstruction
+
+Plot_Conditional_incidence_quantiles<-function(Data_frame_1, Pars)
+{
+
+ with(as.list(c(Data_frame_1, Pars)), {
+
+ # Pars is the list of 3 variables of interest
+ # Select the variables of study:
+
+ variable_x<-Pars[1]
+ variable_y<-Pars[2]
+ variable_z<-Pars[3]
+
+ # Select the file with relevant information:
+ Conditional_incidence_quantiles <- Data_frame_1[[1]]
+
+
+ #Conditional_incidence_quantiles<-Conditional_incidence_quantiles[,-1]
+ Conditional_incidence_quantiles$incidence<-Conditional_incidence_quantiles$counts/Conditional_incidence_quantiles$residents_tot
+ # Conditional_incidence_quantiles<-na.omit(Conditional_incidence_quantiles)
+ #For visual purposes, we are removing cases when the counts were too low resulting in exceptional large incidence
+ #Conditional_incidence_quantiles<-subset(Conditional_incidence_quantiles, Conditional_incidence_quantiles$counts>100)
+
+
+ #number_days<-as.Date(last_day)-as.Date(first_day)
+ #number_days<-as.numeric(number_days)
+ #number_year<-number_days/365
+
+ norm<-1/(1e6)
+
+
+
+ ################### Choose a specific day-length ###################
+
+
+ #Conditional_incidence_quantiles<-na.omit(Conditional_incidence_quantiles)
+ ############################# General variables ###########################
+ if (variable_x=="Relative_humidity"){variable_x_char<-"Relative humidity (%)"}
+ if (variable_y=="Relative_humidity"){variable_y_char<-"Relative humidity (%)"}
+ if (variable_x=="Mean_wind_speed"){variable_x_char<-"Mean Wind Speed (knots)"}
+ if (variable_y=="Mean_wind_speed"){variable_y_char<-"Mean Wind Speed (knots)"}
+ if (variable_x=="Mean_wind_speed"){variable_x_char<-"Mean Wind Speed (knots)"}
+ if (variable_y=="Mean_wind_speed"){variable_y_char<-"Mean Wind Speed (knots)"}
+ if (variable_x=="Mean_Precipitation"){variable_x_char<-"Mean Precipitation (mm)"}
+ if (variable_y=="Mean_Precipitation"){variable_y_char<-"Mean Precipitation (mm)"}
+ if (variable_x=="daylength"){variable_x_char<-"Day-Length (hours)"}
+ if (variable_y=="daylength"){variable_y_char<-"Day-Length (hours)"}
+ if (variable_x=="Minimum_air_temperature"){variable_x_char<-expression("Minimum Air Temperature"*~degree*C)}
+ if (variable_y=="Minimum_air_temperature"){variable_y_char<-expression("Minimum Air Temperature"*~degree*C)}
+ if (variable_x=="Maximum_air_temperature"){variable_x_char<-expression("Maximum Air Temperature"*~degree*C)}
+ if (variable_y=="Maximum_air_temperature"){variable_y_char<-expression("Maximum Air Temperature"*~degree*C)}
+
+
+ labelling_char<-c()
+
+
+
+ ############# 3 variable from here
+ if (variable_z=="daylength"){
+
+
+ for (i_light in c(1:length(unique(Conditional_incidence_quantiles$daylength)))) {
+
+ light_label<-round(as.numeric(unique(Conditional_incidence_quantiles$daylength)))
+
+ if (i_light<length(light_label)){
+ light_label_char<-paste("Day-length ",light_label[i_light],"-",light_label[i_light+1]," hours",sep="")
+ } else {
+ light_label_char<-paste("Day-length >",light_label[i_light]," hours",sep="")
+ }
+
+ wt<-which(Conditional_incidence_quantiles$daylength==unique(Conditional_incidence_quantiles$daylength)[i_light])
+ Conditional_incidence_df<-Conditional_incidence_quantiles[wt,]
+
+ Conditional_incidence_df<-Conditional_incidence_df[order(Conditional_incidence_df[variable_y]),]
+
+ conf_minus<-((Conditional_incidence_df$incidence-1.96*sqrt(Conditional_incidence_df$incidence*(1-Conditional_incidence_df$incidence)/Conditional_incidence_df$residents_tot)))/norm
+ conf_plus<- ((Conditional_incidence_df$incidence+1.96*sqrt(Conditional_incidence_df$incidence*(1-Conditional_incidence_df$incidence)/Conditional_incidence_df$residents_tot)))/norm
+ Conditional_incidence_df$incidence<-Conditional_incidence_df$incidence/norm
+
+
+ Conditional_incidence_df$conf_minus<-conf_minus
+ Conditional_incidence_df$conf_plus<-conf_plus
+ #Conditional_incidence_df<-na.omit(Conditional_incidence_df)
+
+ Conditional_incidence_df[variable_y]<-round(Conditional_incidence_df[variable_y])
+
+
+ if (variable_x=="Maximum_air_temperature" & variable_y=="Relative_humidity"){
+ Conditional_incidence_plot <-
+ ggplot(Conditional_incidence_df,aes(x=Maximum_air_temperature,y=incidence,color=factor(Relative_humidity),fill=factor(Relative_humidity)))
+ Conditional_incidence_plot <- Conditional_incidence_plot+ scale_y_continuous(limits=c(0,5))+scale_x_continuous(limits=c(-5,30))
+ }
+ if (variable_x=="Relative_humidity" & variable_y=="Maximum_air_temperature"){
+ Conditional_incidence_plot <-
+ ggplot(Conditional_incidence_df,aes(x=Relative_humidity,y=incidence,color=factor(Maximum_air_temperature),fill=factor(Maximum_air_temperature)))
+ Conditional_incidence_plot <- Conditional_incidence_plot+ scale_y_continuous(limits=c(0,5))+scale_x_continuous(limits=c(60,100))
+ }
+ if (variable_x=="Maximum_air_temperature" & variable_y=="Mean_wind_speed"){
+ Conditional_incidence_plot <-
+ ggplot(Conditional_incidence_df,aes(x=Maximum_air_temperature,y=incidence,color=factor(Mean_wind_speed),fill=factor(Mean_wind_speed)))
+ Conditional_incidence_plot <- Conditional_incidence_plot+ scale_y_continuous(limits=c(0,5))+scale_x_continuous(limits=c(-5,30))
+ }
+ if (variable_x=="Mean_wind_speed" & variable_y=="Maximum_air_temperature"){
+ Conditional_incidence_plot <-
+ ggplot(Conditional_incidence_df,aes(x=Mean_wind_speed,y=incidence,color=factor(Maximum_air_temperature),fill=factor(Maximum_air_temperature)))
+ Conditional_incidence_plot <- Conditional_incidence_plot+ scale_y_continuous(limits=c(0,5))+scale_x_continuous(limits=c(0,15))
+
+ }
+
+ if (variable_x=="Maximum_air_temperature" & variable_y=="Mean_Precipitation"){
+ Conditional_incidence_plot <-
+ ggplot(Conditional_incidence_df,aes(x=Maximum_air_temperature,y=incidence,color=factor(Mean_Precipitation),fill=factor(Mean_Precipitation)))
+ Conditional_incidence_plot <- Conditional_incidence_plot+ scale_y_continuous(limits=c(0,5))+scale_x_continuous(limits=c(-5,30))
+ }
+ if (variable_x=="Mean_Precipitation" & variable_y=="Maximum_air_temperature"){
+ Conditional_incidence_plot <-
+ ggplot(Conditional_incidence_df,aes(x=Mean_Precipitation,y=incidence,color=factor(Maximum_air_temperature),fill=factor(Maximum_air_temperature)))
+ Conditional_incidence_plot <- Conditional_incidence_plot+ scale_y_continuous(limits=c(0,5))+scale_x_continuous(limits=c(0,5))
+ }
+
+ if (variable_x=="Maximum_air_temperature" & variable_y=="daylength"){
+ Conditional_incidence_plot <-
+ ggplot(Conditional_incidence_df,aes(x=Maximum_air_temperature,y=incidence,color=factor(daylength),fill=factor(daylength)))
+ Conditional_incidence_plot <- Conditional_incidence_plot+ scale_y_continuous(limits=c(0,5))+scale_x_continuous(limits=c(-5,30))
+ }
+
+
+
+
+ Conditional_incidence_plot <- Conditional_incidence_plot+theme_bw(20)
+ Conditional_incidence_plot <- Conditional_incidence_plot+ scale_colour_discrete(name =variable_y_char)
+
+ Conditional_incidence_plot <- Conditional_incidence_plot+ylab("Conditional Incidence of Campylobacter Cases ")
+ Conditional_incidence_plot <- Conditional_incidence_plot+xlab(variable_x_char)
+ Conditional_incidence_plot <- Conditional_incidence_plot+ ggtitle(light_label_char)
+
+ labelling<-c(round(as.numeric(unique(Conditional_incidence_df[variable_y])[,1])))
+ n_labelling<-length(labelling)
+
+
+ for (j in c(1:(n_labelling-1))) {
+
+ labelling_char[j]<-paste(as.character(labelling[j]), "-", as.character(labelling[j+1]) )
+
+ }
+ labelling_char[n_labelling]<-paste(">",as.character(labelling[n_labelling]))
+
+ Conditional_incidence_plot <- Conditional_incidence_plot+ scale_fill_discrete(name =variable_y_char, labels=labelling_char)
+ Conditional_incidence_plot <- Conditional_incidence_plot+ geom_point(size=3.5)+ guides(colour="none") # guides(colour="none")
+ Conditional_incidence_plot <- Conditional_incidence_plot+ geom_line(size=1)
+ Conditional_incidence_plot <- Conditional_incidence_plot+ geom_ribbon(aes(ymin=conf_minus, ymax=conf_plus),alpha=0.35)
+
+ Conditional_incidence_plot <- Conditional_incidence_plot+
+ theme(legend.position= c(0.25,0.75),legend.title = element_text( size = 11),
+ legend.text = element_text( size = 10),legend.background = element_blank(),
+ legend.key=element_rect(colour=NA))
+
+ Conditional_incidence_plot <- Conditional_incidence_plot+ theme(axis.title.x =element_text( colour="#990000", size=15))
+ Conditional_incidence_plot <- Conditional_incidence_plot+ theme(axis.title.y =element_text( colour="#990000", size=15))
+
+
+
+ Conditional_incidence_plot <- Conditional_incidence_plot+ theme(axis.title.x =element_text(size=15))
+ Conditional_incidence_plot <- Conditional_incidence_plot+ theme(axis.title.y =element_text(size=15))
+
+ print(Conditional_incidence_plot)
+
+
+ file_name<-paste("../Graphs/Campylobacter_",variable_x,"_",variable_y,"_",variable_z,"_",
+ round(unique(Conditional_incidence_quantiles$daylength)[i_light]),"_",time_lag_char,"_quantile.tiff", sep = "")
+ tiff(filename = file_name,width = 17.35, height = 17.35, units = "cm", pointsize = 12, res = 1200,compression = "lzw",antialias="default")
+
+ print(Conditional_incidence_plot)
+
+ dev.off()
+
+ }
+
+
+ }
+
+ if (variable_z=="Relative_humidity"){
+
+
+ for (i_light in c(1:length(unique(Conditional_incidence_quantiles$Relative_humidity)))) {
+
+ light_label<-round(as.numeric(unique(Conditional_incidence_quantiles$Relative_humidity)))
+
+ if (i_light<length(light_label)){
+ light_label_char<-paste("Relative Humidity ",light_label[i_light],"-",light_label[i_light+1]," %",sep="")
+ } else {
+ light_label_char<-paste("Relative Humidity >",light_label[i_light]," %",sep="")
+ }
+
+ wt<-which(Conditional_incidence_quantiles$Relative_humidity==unique(Conditional_incidence_quantiles$Relative_humidity)[i_light])
+ Conditional_incidence_df<-Conditional_incidence_quantiles[wt,]
+
+ Conditional_incidence_df<-Conditional_incidence_df[order(Conditional_incidence_df[variable_y]),]
+
+ conf_minus<-((Conditional_incidence_df$incidence-1.96*sqrt(Conditional_incidence_df$incidence*(1-Conditional_incidence_df$incidence)/Conditional_incidence_df$residents_tot)))/norm
+ conf_plus<- ((Conditional_incidence_df$incidence+1.96*sqrt(Conditional_incidence_df$incidence*(1-Conditional_incidence_df$incidence)/Conditional_incidence_df$residents_tot)))/norm
+ Conditional_incidence_df$incidence<-Conditional_incidence_df$incidence/norm
+
+
+ Conditional_incidence_df$conf_minus<-conf_minus
+ Conditional_incidence_df$conf_plus<-conf_plus
+ #Conditional_incidence_df<-na.omit(Conditional_incidence_df)
+
+ Conditional_incidence_df[variable_y]<-round(Conditional_incidence_df[variable_y])
+
+
+ if (variable_x=="Maximum_air_temperature" & variable_y=="Relative_humidity"){
+ Conditional_incidence_plot <-
+ ggplot(Conditional_incidence_df,aes(x=Maximum_air_temperature,y=incidence,color=factor(Relative_humidity),fill=factor(Relative_humidity)))
+ Conditional_incidence_plot <- Conditional_incidence_plot+ scale_y_continuous(limits=c(0,5))+scale_x_continuous(limits=c(-5,30))
+ }
+ if (variable_x=="Relative_humidity" & variable_y=="Maximum_air_temperature"){
+ Conditional_incidence_plot <-
+ ggplot(Conditional_incidence_df,aes(x=Relative_humidity,y=incidence,color=factor(Maximum_air_temperature),fill=factor(Maximum_air_temperature)))
+ Conditional_incidence_plot <- Conditional_incidence_plot+ scale_y_continuous(limits=c(0,5))+scale_x_continuous(limits=c(60,100))
+ }
+ if (variable_x=="Maximum_air_temperature" & variable_y=="Mean_wind_speed"){
+ Conditional_incidence_plot <-
+ ggplot(Conditional_incidence_df,aes(x=Maximum_air_temperature,y=incidence,color=factor(Mean_wind_speed),fill=factor(Mean_wind_speed)))
+ Conditional_incidence_plot <- Conditional_incidence_plot+ scale_y_continuous(limits=c(0,5))+scale_x_continuous(limits=c(-5,30))
+ }
+ if (variable_x=="Mean_wind_speed" & variable_y=="Maximum_air_temperature"){
+ Conditional_incidence_plot <-
+ ggplot(Conditional_incidence_df,aes(x=Mean_wind_speed,y=incidence,color=factor(Maximum_air_temperature),fill=factor(Maximum_air_temperature)))
+ Conditional_incidence_plot <- Conditional_incidence_plot+ scale_y_continuous(limits=c(0,5))+scale_x_continuous(limits=c(0,15))
+
+ }
+
+ if (variable_x=="Maximum_air_temperature" & variable_y=="Mean_Precipitation"){
+ Conditional_incidence_plot <-
+ ggplot(Conditional_incidence_df,aes(x=Maximum_air_temperature,y=incidence,color=factor(Mean_Precipitation),fill=factor(Mean_Precipitation)))
+ Conditional_incidence_plot <- Conditional_incidence_plot+ scale_y_continuous(limits=c(0,5))+scale_x_continuous(limits=c(-5,30))
+ }
+ if (variable_x=="Mean_Precipitation" & variable_y=="Maximum_air_temperature"){
+ Conditional_incidence_plot <-
+ ggplot(Conditional_incidence_df,aes(x=Mean_Precipitation,y=incidence,color=factor(Maximum_air_temperature),fill=factor(Maximum_air_temperature)))
+ Conditional_incidence_plot <- Conditional_incidence_plot+ scale_y_continuous(limits=c(0,5))+scale_x_continuous(limits=c(0,5))
+ }
+
+ if (variable_x=="Maximum_air_temperature" & variable_y=="daylength"){
+ Conditional_incidence_plot <-
+ ggplot(Conditional_incidence_df,aes(x=Maximum_air_temperature,y=incidence,color=factor(daylength),fill=factor(daylength)))
+ Conditional_incidence_plot <- Conditional_incidence_plot+ scale_y_continuous(limits=c(0,5))+scale_x_continuous(limits=c(-5,30))
+ }
+ if (variable_x=="daylength" & variable_y=="Maximum_air_temperature"){
+ Conditional_incidence_plot <-
+ ggplot(Conditional_incidence_df,aes(x=daylength,y=incidence,color=factor(Maximum_air_temperature),fill=factor(Maximum_air_temperature)))
+ Conditional_incidence_plot <- Conditional_incidence_plot+ scale_y_continuous(limits=c(0,5))+scale_x_continuous(limits=c(6,18))
+ }
+
+
+
+
+
+ Conditional_incidence_plot <- Conditional_incidence_plot+theme_bw(20)
+ Conditional_incidence_plot <- Conditional_incidence_plot+ scale_colour_discrete(name =variable_y_char)
+
+ Conditional_incidence_plot <- Conditional_incidence_plot+ylab("Conditional Incidence of Campylobacter Cases ")
+ Conditional_incidence_plot <- Conditional_incidence_plot+xlab(variable_x_char)
+ Conditional_incidence_plot <- Conditional_incidence_plot+ ggtitle(light_label_char)
+
+ labelling<-c(round(as.numeric(unique(Conditional_incidence_df[variable_y])[,1])))
+ n_labelling<-length(labelling)
+
+
+ for (j in c(1:(n_labelling-1))) {
+
+ labelling_char[j]<-paste(as.character(labelling[j]), "-", as.character(labelling[j+1]) )
+
+ }
+ labelling_char[n_labelling]<-paste(">",as.character(labelling[n_labelling]))
+
+ Conditional_incidence_plot <- Conditional_incidence_plot+ scale_fill_discrete(name =variable_y_char, labels=labelling_char)
+ Conditional_incidence_plot <- Conditional_incidence_plot+ geom_point(size=3.5)+ guides(colour="none") # guides(colour="none")
+ Conditional_incidence_plot <- Conditional_incidence_plot+ geom_line(size=1)
+ Conditional_incidence_plot <- Conditional_incidence_plot+ geom_ribbon(aes(ymin=conf_minus, ymax=conf_plus),alpha=0.35)
+
+ Conditional_incidence_plot <- Conditional_incidence_plot+
+ theme(legend.position= c(0.25,0.75),legend.title = element_text( size = 11),
+ legend.text = element_text( size = 10),legend.background = element_blank(),
+ legend.key=element_rect(colour=NA))
+
+ Conditional_incidence_plot <- Conditional_incidence_plot+ theme(axis.title.x =element_text( colour="#990000", size=15))
+ Conditional_incidence_plot <- Conditional_incidence_plot+ theme(axis.title.y =element_text( colour="#990000", size=15))
+
+
+
+ Conditional_incidence_plot <- Conditional_incidence_plot+ theme(axis.title.x =element_text(size=15))
+ Conditional_incidence_plot <- Conditional_incidence_plot+ theme(axis.title.y =element_text(size=15))
+
+ print(Conditional_incidence_plot)
+
+
+ file_name<-paste("../Graphs/Campylobacter_",variable_x,"_",variable_y,"_",variable_z,"_",
+ round(unique(Conditional_incidence_quantiles$Relative_humidity)[i_light]),"_",time_lag_char,"_quantile.tiff", sep = "")
+ tiff(filename = file_name,width = 17.35, height = 17.35, units = "cm", pointsize = 12, res = 1200,compression = "lzw",antialias="default")
+
+ print(Conditional_incidence_plot)
+
+ dev.off()
+
+ }
+
+
+ }
+
+ })
+}
+
+
+
+
+Plot_Conditional_incidence_quantiles_two_factors<-function(Data_frame_1, Pars)
+{
+
+ with(as.list(c(Data_frame_1, Pars)), {
+
+ # Pars is the list of 3 variables of interest
+ # Select the variables of study:
+
+ variable_x<-Pars[1]
+ variable_y<-Pars[2]
+ variable_z<-Pars[3]
+
+ # Select the file with relevant information:
+ Conditional_incidence_quantiles <- Data_frame_1[[1]]
+
+
+ #Conditional_incidence_quantiles<-Conditional_incidence_quantiles[,-1]
+ Conditional_incidence_quantiles$incidence<-Conditional_incidence_quantiles$counts/Conditional_incidence_quantiles$residents_tot
+ # Conditional_incidence_quantiles<-na.omit(Conditional_incidence_quantiles)
+ #For visual purposes, we are removing cases when the counts were too low resulting in exceptional large incidence
+ #Conditional_incidence_quantiles<-subset(Conditional_incidence_quantiles, Conditional_incidence_quantiles$counts>100)
+
+
+ #number_days<-as.Date(last_day)-as.Date(first_day)
+ #number_days<-as.numeric(number_days)
+ #number_year<-number_days/365
+
+ norm<-1/(1e6)
+
+
+
+ ################### Choose a specific day-length ###################
+
+
+ #Conditional_incidence_quantiles<-na.omit(Conditional_incidence_quantiles)
+ ############################# General variables ###########################
+ if (variable_x=="Relative_humidity"){variable_x_char<-"Relative humidity (%)"}
+ if (variable_y=="Relative_humidity"){variable_y_char<-"Relative humidity (%)"}
+ if (variable_x=="Mean_wind_speed"){variable_x_char<-"Mean Wind Speed (knots)"}
+ if (variable_y=="Mean_wind_speed"){variable_y_char<-"Mean Wind Speed (knots)"}
+ if (variable_x=="Mean_wind_speed"){variable_x_char<-"Mean Wind Speed (knots)"}
+ if (variable_y=="Mean_wind_speed"){variable_y_char<-"Mean Wind Speed (knots)"}
+ if (variable_x=="Mean_Precipitation"){variable_x_char<-"Mean Precipitation (mm)"}
+ if (variable_y=="Mean_Precipitation"){variable_y_char<-"Mean Precipitation (mm)"}
+ if (variable_x=="daylength"){variable_x_char<-"Day-Length (hours)"}
+ if (variable_y=="daylength"){variable_y_char<-"Day-Length (hours)"}
+ if (variable_x=="Minimum_air_temperature"){variable_x_char<-expression("Minimum Air Temperature"*~degree*C)}
+ if (variable_y=="Minimum_air_temperature"){variable_y_char<-expression("Minimum Air Temperature"*~degree*C)}
+ if (variable_x=="Maximum_air_temperature"){variable_x_char<-expression("Maximum Air Temperature"*~degree*C)}
+ if (variable_y=="Maximum_air_temperature"){variable_y_char<-expression("Maximum Air Temperature"*~degree*C)}
+
+
+ labelling_char<-c()
+
+
+
+ ############# variables from here
+
+ Conditional_incidence_df<-Conditional_incidence_quantiles
+
+ Conditional_incidence_df<-Conditional_incidence_df[order(Conditional_incidence_df[variable_y]),]
+
+ conf_minus<-((Conditional_incidence_df$incidence-1.96*sqrt(Conditional_incidence_df$incidence*(1-Conditional_incidence_df$incidence)/Conditional_incidence_df$residents_tot)))/norm
+ conf_plus<- ((Conditional_incidence_df$incidence+1.96*sqrt(Conditional_incidence_df$incidence*(1-Conditional_incidence_df$incidence)/Conditional_incidence_df$residents_tot)))/norm
+ Conditional_incidence_df$incidence<-Conditional_incidence_df$incidence/norm
+
+
+ Conditional_incidence_df$conf_minus<-conf_minus
+ Conditional_incidence_df$conf_plus<-conf_plus
+
+ Conditional_incidence_df[variable_y]<-round(Conditional_incidence_df[variable_y])
+
+ if (variable_x=="Maximum_air_temperature" & variable_y=="Relative_humidity"){
+ Conditional_incidence_plot <-
+ ggplot(Conditional_incidence_df,aes(x=Maximum_air_temperature,y=incidence,color=factor(Relative_humidity),fill=factor(Relative_humidity)))
+ Conditional_incidence_plot <- Conditional_incidence_plot+ scale_y_continuous(limits=c(1,4))+scale_x_continuous(limits=c(-5,30))
+
+ }
+ if (variable_x=="Maximum_air_temperature" & variable_y=="Mean_wind_speed"){
+ Conditional_incidence_plot <-
+ ggplot(Conditional_incidence_df,aes(x=Maximum_air_temperature,y=incidence,color=factor(Mean_wind_speed),fill=factor(Mean_wind_speed)))
+ Conditional_incidence_plot <- Conditional_incidence_plot+ scale_y_continuous(limits=c(1,4))+scale_x_continuous(limits=c(-5,30))
+ }
+ if (variable_x=="Maximum_air_temperature" & variable_y=="Mean_Precipitation"){
+ Conditional_incidence_plot <-
+ ggplot(Conditional_incidence_df,aes(x=Maximum_air_temperature,y=incidence,color=factor(Mean_Precipitation),fill=factor(Mean_Precipitation)))
+ Conditional_incidence_plot <- Conditional_incidence_plot+ scale_y_continuous(limits=c(1,4))+scale_x_continuous(limits=c(-5,30))
+ }
+ if (variable_x=="Maximum_air_temperature" & variable_y=="daylength"){
+ Conditional_incidence_plot <-
+ ggplot(Conditional_incidence_df,aes(x=Maximum_air_temperature,y=incidence,color=factor(daylength),fill=factor(daylength)))
+ Conditional_incidence_plot <- Conditional_incidence_plot+ scale_y_continuous(limits=c(1,4))+scale_x_continuous(limits=c(-5,30))
+ }
+
+
+ if (variable_x=="Minimum_air_temperature" & variable_y=="Relative_humidity"){
+ Conditional_incidence_plot <-
+ ggplot(Conditional_incidence_df,aes(x=Minimum_air_temperature,y=incidence,color=factor(Relative_humidity),fill=factor(Relative_humidity)))
+ Conditional_incidence_plot <- Conditional_incidence_plot+ scale_y_continuous(limits=c(1,4))+scale_x_continuous(limits=c(-10,25))
+ }
+ if (variable_x=="Minimum_air_temperature" & variable_y=="Mean_wind_speed"){
+ Conditional_incidence_plot <-
+ ggplot(Conditional_incidence_df,aes(x=Minimum_air_temperature,y=incidence,color=factor(Mean_wind_speed),fill=factor(Mean_wind_speed)))
+ Conditional_incidence_plot <- Conditional_incidence_plot+ scale_y_continuous(limits=c(1,4))+scale_x_continuous(limits=c(-10,25))
+
+ }
+ if (variable_x=="Minimum_air_temperature" & variable_y=="Mean_Precipitation"){
+ Conditional_incidence_plot <-
+ ggplot(Conditional_incidence_df,aes(x=Minimum_air_temperature,y=incidence,color=factor(Mean_Precipitation),fill=factor(Mean_Precipitation)))
+ Conditional_incidence_plot <- Conditional_incidence_plot+ scale_y_continuous(limits=c(1,4))+scale_x_continuous(limits=c(-10,25))
+
+ }
+ if (variable_x=="Minimum_air_temperature" & variable_y=="daylength"){
+ Conditional_incidence_plot <-
+ ggplot(Conditional_incidence_df,aes(x=Minimum_air_temperature,y=incidence,color=factor(daylength),fill=factor(daylength)))
+ Conditional_incidence_plot <- Conditional_incidence_plot+ scale_y_continuous(limits=c(1,4))+scale_x_continuous(limits=c(-10,25))
+
+ }
+
+ Conditional_incidence_plot <- Conditional_incidence_plot+theme_bw(20)
+ Conditional_incidence_plot <- Conditional_incidence_plot+ scale_colour_discrete(name =variable_y_char)
+
+ Conditional_incidence_plot <- Conditional_incidence_plot+ylab("Conditional Incidence of Campylobacter Cases ")
+ Conditional_incidence_plot <- Conditional_incidence_plot+xlab(variable_x_char)
+
+ labelling<-c(round(as.numeric(unique(Conditional_incidence_df[variable_y])[,1])))
+ n_labelling<-length(labelling)
+
+
+ for (j in c(1:(n_labelling-1))) {
+
+ labelling_char[j]<-paste(as.character(labelling[j]), "-", as.character(labelling[j+1]) )
+
+ }
+ labelling_char[n_labelling]<-paste(">",as.character(labelling[n_labelling]))
+
+ Conditional_incidence_plot <- Conditional_incidence_plot+ scale_fill_discrete(name =variable_y_char, labels=labelling_char)
+ Conditional_incidence_plot <- Conditional_incidence_plot+ geom_point(size=3.5)+ guides(colour="none") # guides(colour="none")
+ Conditional_incidence_plot <- Conditional_incidence_plot+ geom_line(size=1)
+ Conditional_incidence_plot <- Conditional_incidence_plot+ geom_ribbon(aes(ymin=conf_minus, ymax=conf_plus),alpha=0.35)
+
+ Conditional_incidence_plot <- Conditional_incidence_plot+
+ theme(legend.position= c(0.25,0.75),legend.title = element_text( size = 11),
+ legend.text = element_text( size = 10),legend.background = element_blank(),
+ legend.key=element_rect(colour=NA))
+
+ Conditional_incidence_plot <- Conditional_incidence_plot+ theme(axis.title.x =element_text( colour="#990000", size=15))
+ Conditional_incidence_plot <- Conditional_incidence_plot+ theme(axis.title.y =element_text( colour="#990000", size=15))
+
+
+
+ Conditional_incidence_plot <- Conditional_incidence_plot+ theme(axis.title.x =element_text(size=15))
+ Conditional_incidence_plot <- Conditional_incidence_plot+ theme(axis.title.y =element_text(size=15))
+ Conditional_incidence_plot
+
+
+ file_name<-paste("../Graphs/Campylobacter_",variable_x,"_",variable_y,"_",time_lag_char,"_quantile.tiff", sep = "")
+ tiff(filename = file_name,width = 17.35, height = 17.35, units = "cm", pointsize = 12, res = 1200,compression = "lzw",antialias="default")
+
+ print(Conditional_incidence_plot)
+
+ dev.off()
+
+
+
+
+
+ })
+}
+
+
+
+