diff --git a/Codes/Conditional_incidence_and_reconstruction_factors_4_factors.R b/Codes/Conditional_incidence_and_reconstruction_factors_4_factors.R
index a7d4b2934fccc3ceec7555f6edbd40a782fd3a79..47cfd4be244dfbba18ac1ff7fbdf3e7b63049c9f 100644
--- a/Codes/Conditional_incidence_and_reconstruction_factors_4_factors.R
+++ b/Codes/Conditional_incidence_and_reconstruction_factors_4_factors.R
@@ -90,13 +90,20 @@ source("Function_Input_files.R", chdir=TRUE)
#
source("Function_Conditional_Incidence_Quantiles.R", chdir=TRUE)
#Conditional incidence, bins divided in quantiles
+Conditional_incidence_quantiles_info2<-c()
+
+
for (quarter in c(1:4)){
Pars_cond_prev_quantiles<- c(variable_x,variable_y,variable_z,"_Quantiles.csv",by_x,by_y,by_z,quarter)
Conditional_incidence_quantiles_info<-Conditional_incidence_Quantiles_4_factors(Env_laboratory_data,Env_Pathogen_data,Pars_cond_prev_quantiles)
+Conditional_incidence_quantiles_info2<-rbind(Conditional_incidence_quantiles_info2,Conditional_incidence_quantiles_info[[1]])
+
## plot Conditional incidence for bins divided in quantiles
-Pars_plot_cond_prev_quantiles <- c(variable_x,variable_y,variable_z,quarter)
-source("Function_Plot_Conditional_Incidence_Quantiles.R", chdir=TRUE)
-Plot_Conditional_incidence_quantiles_4_factors(Conditional_incidence_quantiles_info,Pars_plot_cond_prev_quantiles)
}
+
+Pars_plot_cond_prev_quantiles <- c(variable_x,variable_y,variable_z)
+source("Function_Plot_Conditional_Incidence_Quantiles.R", chdir=TRUE)
+Plot_Conditional_incidence_quantiles_4_factors(Conditional_incidence_quantiles_info2,Pars_plot_cond_prev_quantiles)
+
diff --git a/Codes/Function_Plot_Conditional_Incidence_Quantiles.R b/Codes/Function_Plot_Conditional_Incidence_Quantiles.R
index 8133f4e5fd5b39e7c688f6c2e150491be888a04c..4fb9ba6c617a55c9e11c7988a6370315349b3615 100644
--- a/Codes/Function_Plot_Conditional_Incidence_Quantiles.R
+++ b/Codes/Function_Plot_Conditional_Incidence_Quantiles.R
@@ -15,14 +15,14 @@ Plot_Conditional_incidence_quantiles<-function(Data_frame_1, Pars)
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)
+ # 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)
@@ -58,8 +58,8 @@ Plot_Conditional_incidence_quantiles<-function(Data_frame_1, Pars)
labelling_char<-c()
-
-
+
+
############# 3 variable from here
if (variable_z=="daylength"){
@@ -82,12 +82,12 @@ Plot_Conditional_incidence_quantiles<-function(Data_frame_1, Pars)
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])
@@ -95,95 +95,95 @@ Plot_Conditional_incidence_quantiles<-function(Data_frame_1, Pars)
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(0,25))
- }
+ }
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(0,25))
- }
+ }
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"){
+ 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(0,25))
- }
+ }
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(0,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)
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)
+
+
+ for (j in c(1:(n_labelling-1))) {
- dev.off()
+ 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"){
@@ -315,7 +315,7 @@ Plot_Conditional_incidence_quantiles<-function(Data_frame_1, Pars)
}
- })
+ })
}
@@ -379,7 +379,224 @@ Plot_Conditional_incidence_quantiles_two_factors<-function(Data_frame_1, Pars)
############# variables from here
- Conditional_incidence_df<-Conditional_incidence_quantiles
+ 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,25))
+
+ }
+ 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(0,25))
+ }
+ 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(0,25))
+ }
+ 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(0,25))
+ }
+ if (variable_x=="Relative_humidity" & variable_y=="daylength"){
+ Conditional_incidence_plot <-
+ ggplot(Conditional_incidence_df,aes(x=Relative_humidity,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(60,100))
+ }
+ if (variable_x=="Mean_Precipitation" & variable_y=="daylength"){
+ Conditional_incidence_plot <-
+ ggplot(Conditional_incidence_df,aes(x=Mean_Precipitation,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(0,10))
+ }
+
+ 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()
+
+
+
+
+
+ })
+}
+
+
+
+Plot_Conditional_incidence_quantiles_4_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_all <- Data_frame_1
+
+ #Conditional_incidence_quantiles<-Conditional_incidence_quantiles[,-1]
+ Conditional_incidence_quantiles_all$incidence<-Conditional_incidence_quantiles_all$counts/Conditional_incidence_quantiles_all$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 Rainfall (mm)"}
+ if (variable_y=="Mean_Precipitation"){variable_y_char<-"Mean Rainfall (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()
+ light_label_char<-c()
+ z_label_char<-c()
+ light_label<-c()
+
+ light_label<-round(as.numeric(unique(Conditional_incidence_quantiles_all$daylength)),1)
+
+
+ for (i_light in c(1:length(unique(Conditional_incidence_quantiles_all$daylength)))) {
+
+ if (i_light==1){
+ light_label_char[i_light]<-paste("Day-length < ",light_label[i_light]," hours",sep="")
+ }
+
+ if (i_light==length(light_label )){
+ light_label_char[i_light]<-paste("Day-length > ",light_label[i_light-1]," hours",sep="")
+ }
+
+ if (i_light!=1 && i_light<length(light_label )){
+ light_label_char[i_light]<-paste("Day-length ",light_label[i_light-1]," - ",light_label[i_light]," hours",sep="")
+ }
+
+ }
+
+ for (i_light in c(1:length(light_label)) ){
+
+ Conditional_incidence_quantiles<-subset(Conditional_incidence_quantiles_all, as.numeric(Conditional_incidence_quantiles_all$daylength)==light_label[i_light])
+ if (variable_z=="Maximum_air_temperature"){
+
+ z_label<-round(as.numeric(unique(Conditional_incidence_quantiles$Maximum_air_temperature)))
+ for (i_z in c(1:length(unique(Conditional_incidence_quantiles$Maximum_air_temperature)))) {
+
+ if (i_z<length(z_label)){
+ z_label_char[i_z]<-paste("Maximum Air Temperature ",z_label[i_z],"-",z_label[i_z+1]," (\u00B0C)",sep="")
+ } else {
+ z_label_char[i_z]<-paste("Maximum Air Temperature > ",z_label[i_z]," (\u00B0C)",sep="")
+ }
+
+
+
+
+ wt<-which(Conditional_incidence_quantiles$Maximum_air_temperature==unique(Conditional_incidence_quantiles$Maximum_air_temperature)[i_z])
+ Conditional_incidence_df<-Conditional_incidence_quantiles[wt,]
Conditional_incidence_df<-Conditional_incidence_df[order(Conditional_incidence_df[variable_y]),]
@@ -390,75 +607,45 @@ Plot_Conditional_incidence_quantiles_two_factors<-function(Data_frame_1, Pars)
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(1,4))+scale_x_continuous(limits=c(5,25))
-
- }
- 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(0,25))
- }
- 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(0,25))
- }
- 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(0,25))
- }
- if (variable_x=="Relative_humidity" & variable_y=="daylength"){
+
+ if (variable_x=="Relative_humidity" & variable_y=="Mean_Precipitation"){
Conditional_incidence_plot <-
- ggplot(Conditional_incidence_df,aes(x=Relative_humidity,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(60,100))
+ ggplot(Conditional_incidence_df,aes(x=Relative_humidity,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(60,100))
}
- if (variable_x=="Mean_Precipitation" & variable_y=="daylength"){
+
+ if (variable_x=="Mean_Precipitation" & variable_y=="Relative_humidity"){
Conditional_incidence_plot <-
- ggplot(Conditional_incidence_df,aes(x=Mean_Precipitation,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(0,10))
+ ggplot(Conditional_incidence_df,aes(x=Mean_Precipitation,y=incidence,color=factor(Relative_humidity),fill=factor(Relative_humidity)))
+ Conditional_incidence_plot <- Conditional_incidence_plot+ scale_y_continuous(limits=c(60,100))+scale_x_continuous(limits=c(0,5))
}
- if (variable_x=="Minimum_air_temperature" & variable_y=="Relative_humidity"){
+ if (variable_x=="Relative_humidity" & variable_y=="daylength"){
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))
+ ggplot(Conditional_incidence_df,aes(x=Relative_humidity,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(60,100))
}
- 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)
+ Conditional_incidence_plot <- Conditional_incidence_plot+ ggtitle(z_label_char[i_z])
+
+ Conditional_incidence_plot <- Conditional_incidence_plot+geom_text(aes(x = 72, y = 4.75,
+ label = light_label_char[i_light]),colour="black",size=7)
+
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]) )
@@ -483,10 +670,12 @@ Plot_Conditional_incidence_quantiles_two_factors<-function(Data_frame_1, Pars)
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
+
+ print(Conditional_incidence_plot)
- file_name<-paste("../Graphs/Campylobacter_",variable_x,"_",variable_y,"_",time_lag_char,"_quantile.tiff", sep = "")
+ file_name<-paste("../Graphs/Campylobacter_",variable_x,"_",variable_y,"_",variable_z,"_daylength_",
+ z_label[i_z],"_quarter_",i_light,"_",time_lag_char,".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)
@@ -494,12 +683,20 @@ Plot_Conditional_incidence_quantiles_two_factors<-function(Data_frame_1, Pars)
dev.off()
-
-
+
+
+
+
+ }
+
+ }
+ }
})
}
+
+