From b7c3cddb32cd0cb77f91edbcc627ab4fc8d347be Mon Sep 17 00:00:00 2001
From: gl0020 <g.loiacono@surrey.ac.uk>
Date: Fri, 5 May 2023 13:20:10 +0100
Subject: [PATCH] minor corrections
---
...Conditional_incidence_and_reconstruction.R | 2 +-
...nd_reconstruction_for_constant_situation.R | 223 ++++++++++++++++++
...incidence_and_reconstruction_two_factors.R | 6 +-
.../Function_Conditional_Incidence_Uniform.R | 121 +++++++++-
...ction_Plot_Conditional_Incidence_Uniform.R | 197 +++++++++++++++-
Codes/Pathogen_Linkage_fixed_time_lag.R | 2 +-
6 files changed, 535 insertions(+), 16 deletions(-)
create mode 100644 Codes/Conditional_incidence_and_reconstruction_for_constant_situation.R
diff --git a/Codes/Conditional_incidence_and_reconstruction.R b/Codes/Conditional_incidence_and_reconstruction.R
index f3f0f9d..c2699e2 100644
--- a/Codes/Conditional_incidence_and_reconstruction.R
+++ b/Codes/Conditional_incidence_and_reconstruction.R
@@ -118,7 +118,7 @@ Pars_cond_prev_quantiles<- c(variable_x,variable_y,variable_z,"_Quantiles.csv",b
Conditional_incidence_quantiles_info<-Conditional_incidence_Quantiles(Env_laboratory_data,Env_Pathogen_data,Pars_cond_prev_quantiles)
## plot Conditional incidence for bins divided in quantiles
-Pars_plot_cond_prev_quantiles <- c(variable_x,variable_y,variable_z)
+Pars_plot_cond_prev_quantiles <- c(variable_x,variable_y,variable_z,"_quantile")
source("Function_Plot_Conditional_Incidence_Quantiles.R", chdir=TRUE)
Plot_Conditional_incidence_quantiles(Conditional_incidence_quantiles_info,Pars_plot_cond_prev_quantiles)
diff --git a/Codes/Conditional_incidence_and_reconstruction_for_constant_situation.R b/Codes/Conditional_incidence_and_reconstruction_for_constant_situation.R
new file mode 100644
index 0000000..5e77509
--- /dev/null
+++ b/Codes/Conditional_incidence_and_reconstruction_for_constant_situation.R
@@ -0,0 +1,223 @@
+# The code estimate the incidence of the diseases conditional to local environmental variables. It perform a series of diagnostic and
+# then reconstruct the time series of cases given the local environmental factors in England and Wales
+
+### Require Relevant R Packages ###########
+rm(list=ls(all=TRUE))
+
+library(data.table)
+library(readr)
+library(dplyr)
+library(lubridate)
+library(plyr)
+library(dplR)
+library(MALDIquant)
+library(ggplot2)
+library(slider)
+library(WaveletComp)
+library(gdata)
+library(wesanderson)
+library(zoo)
+
+# Global Parameters
+
+# Set length of the desired time-lag
+
+time_lag <<-14
+time_lag_char <<-paste(time_lag)
+
+# Set the years you are focusing on
+
+year_first<<-1990
+year_last<<-2009
+first_day<<-as.Date(paste(as.character(year_first),"-01-01",sep=""))
+last_day<<-as.Date(paste(as.character(year_last),"-12-31",sep=""))
+dates_s<<-seq(first_day,last_day,by="1 day")
+
+## Gather Information on laboratory Postcodes
+
+variable_df_1<-read.csv("../Data_Base/Laboratories_Post_Codes.csv")
+
+All_PC_s<-variable_df_1[,1]
+dates<-rep(dates_s,times=length(All_PC_s))
+All_PC<<-rep(All_PC_s,each=length(dates_s))
+
+
+# Selection of the variables of interest.
+#Change HERE depending on the variable object of study, by copy-pasting from the list below
+
+variable_x<-"Maximum_air_temperature"
+variable_y<-"Relative_humidity"
+#variable_y<-"Mean_Precipitation"
+#variable_x<-"Mean_wind_speed"
+variable_z<-"daylength"
+#variable_z<-"Relative_humidity"
+
+### Set the size of the bins for Uniform conditional incidence and resolution in rounding the values of the bins.
+### For example:
+## if the size of the bin for humidity is 5 and we want the interval of he bins like ...70,75,80... we choose delta_variable<-5, round<-0
+## if the size of the bin for precipitation is 0.1 and we want the interval of he bins like ...0,0.1.0.2 ..., we choose delta_variable<-0.1, round_x<-1
+## to distinguishes one decimal point
+
+delta_variable_x<-1
+round_x<-0
+delta_variable_y<-3
+round_y<-0
+delta_variable_z<-0.5
+round_z<-1
+## partitions of the quantiles
+by_z<-0.25
+by_y<-0.25
+by_x<-0.1
+
+# Potential Variables to use:
+#"Maximum_air_temperature",
+#"Minimum_air_temperature",
+#"Mean_wind_speed",
+#"Mean_Precipitation"
+#"Relative_humidity"
+#"daylength"
+
+
+
+
+#const_variable<-"const_temp_20_hum_76"
+#const_variable<-"const_hum_76_daylength_15"
+const_variable<-"const_temp_20_daylength_15"
+
+
+if (const_variable=="const_hum_76_daylength_15"){
+ Relative_humidity_value<-76
+ daylength_value <-15
+}
+if (const_variable=="const_temp_20_hum_76"){
+ Maximum_air_temperature_value <-20
+ Relative_humidity_value<-76
+}
+if (const_variable=="const_temp_20_daylength_15"){
+ Maximum_air_temperature_value <-20
+ daylength_value <-15
+}
+
+
+
+
+######## Input of data ####
+source("Function_Input_files.R", chdir=TRUE)
+
+
+######## Perform Wavelet analysis for the selected 17 years to check for periodicity in the data ####
+###### helpful but not necessary. So this is commented
+source("Function_Wavelet_analysis.R", chdir=TRUE)
+
+Pars_wavelet_analysis<-c("../Graphs/Power_Spectrum.tiff","../Graphs/Power_Average.tiff","No")
+Env_Pathogen_data$Cases_computed<-Env_Pathogen_data$Cases
+campylobacter_data_national<-wavelet_analysis(Env_Pathogen_data, Pars_wavelet_analysis)
+
+
+###### Calculate Conditional incidence Uniform ######
+# The code does look at how the risk of the disease in humans depends on environmental variables
+# Analysis was done following an uniform division of the range of the environmental variables, independently of the number of observations
+# This is used for the reconstruction
+
+#
+source("Function_Conditional_Incidence_Uniform.R", chdir=TRUE)
+#Conditional incidence, Uniform bins
+Pars_cond_incidence <- c(variable_x,variable_y,variable_z,"_Uniform.csv",delta_variable_x,round_x,delta_variable_y,round_y,delta_variable_z,round_z)
+Conditional_incidence_unif_info<-Conditional_incidence_Uniform(Env_laboratory_data,Env_Pathogen_data,Pars_cond_incidence)
+
+## plot Conditional incidence for uniform, for diagnostic only
+#Pars_plot_cond_prev_uniform <- c(variable_x,variable_y,variable_z,time_lag_char)
+#source("Function_Plot_Conditional_incidence_Uniform.R", chdir=TRUE)
+#Plot_Conditional_incidence_Uniform(Conditional_incidence_unif_info,Pars_plot_cond_prev_uniform)
+
+#
+#source("Function_Conditional_Incidence_Quantiles.R", chdir=TRUE)
+#Conditional incidence, bins divided in quantiles
+#Pars_cond_prev_quantiles<- c(variable_x,variable_y,variable_z,"_Quantiles.csv",by_x,by_y,by_z)
+#Conditional_incidence_quantiles_info<-Conditional_incidence_Quantiles(Env_laboratory_data,Env_Pathogen_data,Pars_cond_prev_quantiles)
+
+if (const_variable=="const_hum_76_daylength_15"){
+
+
+
+
+ ###########
+ wt<-which(abs((Conditional_incidence_unif_info[[1]]$daylength) - daylength_value) ==
+ min(abs((Conditional_incidence_unif_info[[1]]$daylength) - daylength_value)))
+ Conditional_incidence_unif_info_subset1<-Conditional_incidence_unif_info[[1]][wt,]
+
+ ww<-which(abs((Conditional_incidence_unif_info_subset1$Relative_humidity) - Relative_humidity_value) ==
+ min(abs((Conditional_incidence_unif_info_subset1$Relative_humidity) - Relative_humidity_value)))
+ Conditional_incidence_unif_info_subset<-Conditional_incidence_unif_info_subset1[ww,]
+
+ ## plot Conditional incidence for bins divided in unif
+ Pars_plot_cond_prev_unif <- c(variable_x,paste("_",const_variable))
+ source("Function_Plot_Conditional_Incidence_Uniform.R", chdir=TRUE)
+Plot_Conditional_incidence_unif_constant(Conditional_incidence_unif_info_subset,Pars_plot_cond_prev_unif)
+
+}
+
+if (const_variable=="const_temp_20_hum_76"){
+
+ wt<-which(abs((Conditional_incidence_unif_info[[1]]$Maximum_air_temperature) - Maximum_air_temperature_value) == min(abs((Conditional_incidence_unif_info[[1]]$Maximum_air_temperature) - Maximum_air_temperature_value)))
+ Conditional_incidence_unif_info_subset1<-Conditional_incidence_unif_info[[1]][wt,]
+
+ ww<-which(abs((Conditional_incidence_unif_info_subset1$Relative_humidity) - Relative_humidity_value) ==
+ min(abs((Conditional_incidence_unif_info_subset1$Relative_humidity) - Relative_humidity_value)))
+ Conditional_incidence_unif_info_subset<-Conditional_incidence_unif_info_subset1[ww,]
+
+ ## plot Conditional incidence for bins divided in quantiles
+ Pars_plot_cond_prev_unif <- c(variable_z,paste("_",const_variable))
+ source("Function_Plot_Conditional_Incidence_Uniform.R", chdir=TRUE)
+ Plot_Conditional_incidence_unif_constant(Conditional_incidence_unif_info_subset,Pars_plot_cond_prev_unif)
+}
+if (const_variable=="const_temp_20_daylength_15"){
+ wt<-which(abs((Conditional_incidence_unif_info[[1]]$Maximum_air_temperature) - Maximum_air_temperature_value) == min(abs((Conditional_incidence_unif_info[[1]]$Maximum_air_temperature) - Maximum_air_temperature_value)))
+ Conditional_incidence_unif_info_subset1<-Conditional_incidence_unif_info[[1]][wt,]
+
+ ww<-which(abs((Conditional_incidence_unif_info_subset1$daylength) - daylength_value) == min(abs((Conditional_incidence_unif_info_subset1$daylength) - daylength_value)))
+ Conditional_incidence_unif_info_subset<-Conditional_incidence_unif_info_subset1[ww,]
+
+ ## plot Conditional incidence for bins divided in quantiles
+ Pars_plot_cond_prev_unif <- c(variable_y,paste("_",const_variable))
+ source("Function_Plot_Conditional_Incidence_Uniform.R", chdir=TRUE)
+ Plot_Conditional_incidence_unif_constant(Conditional_incidence_unif_info_subset,Pars_plot_cond_prev_unif)
+
+}
+
+
+
+
+
+
+###### Reconstruction ######
+
+
+
+if (const_variable=="const_hum_76_daylength_15"){
+ Env_laboratory_data$Relative_humidity<-Relative_humidity_value
+ Env_laboratory_data$daylength <-daylength_value
+ }
+if (const_variable=="const_temp_20_hum_76"){
+ Env_laboratory_data$Maximum_air_temperature <-Maximum_air_temperature_value
+ Env_laboratory_data$Relative_humidity<-Relative_humidity_value
+}
+if (const_variable=="const_temp_20_daylength_15"){
+ Env_laboratory_data$Maximum_air_temperature <-Maximum_air_temperature_value
+ Env_laboratory_data$daylength <-daylength_value
+}
+
+
+source("Function_Reconstruction.R")
+Pars_reconstruction <- c(variable_x,variable_y,variable_z,"_reconstructed.csv")
+Reconstruction_time_series<-Reconstruction_time_series(Conditional_incidence_unif_info,Env_laboratory_data,Pars_reconstruction)
+
+############## End reconstructions ####################
+
+###### Plot Reconstruction ######
+source("Function_Plot_Reconstruction.R")
+Pars_Plot_Reconstruction <- c(paste("../Graphs/",const_variable,"_",time_lag_char,"_time_series.tiff",sep=""),paste("../Graphs/",const_variable,"_",time_lag_char,"_yearly_average.tiff",sep=""),"Campylobacteriosis")
+Plot_Reconstruction_time_series(Reconstruction_time_series,campylobacter_data_national,Pars_Plot_Reconstruction)
+#
+
+############## End Plot Reconstruction ####################
diff --git a/Codes/Conditional_incidence_and_reconstruction_two_factors.R b/Codes/Conditional_incidence_and_reconstruction_two_factors.R
index 1705017..127e7fb 100644
--- a/Codes/Conditional_incidence_and_reconstruction_two_factors.R
+++ b/Codes/Conditional_incidence_and_reconstruction_two_factors.R
@@ -45,7 +45,7 @@ All_PC<<-rep(All_PC_s,each=length(dates_s))
# Selection of the variables of interest.
#Change HERE depending on the variable object of study, by copy-pasting from the list below
-variable_x<-"Difference_temperature"
+#variable_x<-"Difference_temperature"
#variable_x<-"Minimum_air_temperature"
#variable_y<-"Mean_Precipitation"
variable_y<-"Relative_humidity"
@@ -53,7 +53,7 @@ variable_y<-"Relative_humidity"
#variable_y<-"daylength"
variable_z<-"daylength"
-#variable_x<-"Maximum_air_temperature"
+variable_x<-"Maximum_air_temperature"
### Set the size of the bins for Uniform conditional incidence and resolution in rounding the values of the bins.
### For example:
## if the size of the bin for humidity is 5 and we want the interval of he bins like ...70,75,80... we choose delta_variable<-5, round<-0
@@ -91,7 +91,7 @@ Env_laboratory_data$Difference_temperature<-Env_laboratory_data$Maximum_air_temp
###### helpful but not necessary. So this is commented
source("Function_Wavelet_analysis.R", chdir=TRUE)
-Pars_wavelet_analysis<-c("../Graphs/Power_Spectrum.tiff","../Graphs/Power_Average.tiff")
+Pars_wavelet_analysis<-c("../Graphs/Power_Spectrum.tiff","../Graphs/Power_Average.tiff","No")
Env_Pathogen_data$Cases_computed<-Env_Pathogen_data$Cases
campylobacter_data_national<-wavelet_analysis(Env_Pathogen_data, Pars_wavelet_analysis)
diff --git a/Codes/Function_Conditional_Incidence_Uniform.R b/Codes/Function_Conditional_Incidence_Uniform.R
index bc5240a..f645ee2 100644
--- a/Codes/Function_Conditional_Incidence_Uniform.R
+++ b/Codes/Function_Conditional_Incidence_Uniform.R
@@ -2,7 +2,7 @@
# Analysis was done following an uniform division of the range of the environmental variables, independently of the number of observations
# This is used for the reconstruction
-Conditional_prevalence_Uniform<-function(Data_frame_1, Data_frame_2,Pars)
+Conditional_incidence_Uniform<-function(Data_frame_1, Data_frame_2,Pars)
{
with(as.list(c(Data_frame_1, Data_frame_2,Pars)), {
@@ -47,8 +47,8 @@ Env_Pathogen_data_df <- Data_frame_2[,c("PostCode","Date", variable_x,variable_y
#################### Create a data frame to fill with
-Conditional_prevalence_unif <-data.frame(character(), character(), character(), numeric(), numeric())
-colnames(Conditional_prevalence_unif) <-c(variable_x, variable_y, variable_z,"counts","residents_tot")
+Conditional_incidence_unif <-data.frame(character(), character(), character(), numeric(), numeric())
+colnames(Conditional_incidence_unif) <-c(variable_x, variable_y, variable_z,"counts","residents_tot")
@@ -114,7 +114,7 @@ for (j_z in c(j_z_min:j_z_max)){
)
colnames(data_df) <-c(variable_x, variable_y, variable_z,"counts","residents_tot") #
- Conditional_prevalence_unif <-rbind(Conditional_prevalence_unif, data_df)
+ Conditional_incidence_unif <-rbind(Conditional_incidence_unif, data_df)
}
}
}
@@ -124,8 +124,117 @@ for (j_z in c(j_z_min:j_z_max)){
}
-write.table(Conditional_prevalence_unif, paste("../Data_Base/Conditional_prevalence_",variable_x,"_",variable_y,"_",variable_z,"_",time_lag_char,info,sep=""), row.names = FALSE , sep = ",",eol = "\n")
-return(list(Conditional_prevalence_unif,breaks_x,breaks_y,breaks_z))
+write.table(Conditional_incidence_unif, paste("../Data_Base/Conditional_incidence_",variable_x,"_",variable_y,"_",variable_z,"_",time_lag_char,info,sep=""), row.names = FALSE , sep = ",",eol = "\n")
+return(list(Conditional_incidence_unif,breaks_x,breaks_y,breaks_z))
})
}
+
+
+
+
+Conditional_incidence_Uniform_two_factors<-function(Data_frame_1, Data_frame_2,Pars)
+{
+
+ with(as.list(c(Data_frame_1, Data_frame_2,Pars)), {
+
+ # Reshape the residents information to make a data frame with the postcode, year and number of residents
+
+
+ # Pars is the list of 3 variables of interest
+
+
+
+ variable_x<-Pars[1]
+ variable_y<-Pars[2]
+ info<-Pars[3]
+ delta_x<-as.numeric(Pars[4])
+ round_digit_x<-as.numeric(Pars[5])
+ delta_y<-as.numeric(Pars[6])
+ round_digit_y<-as.numeric(Pars[7])
+
+ # Select the variables of study:
+ Env_Laboratory_data_df <- Data_frame_1[,c("PostCode","Date", variable_x,variable_y,variable_z,
+ "residents")]
+ Env_Pathogen_data_df <- Data_frame_2[,c("PostCode","Date", variable_x,variable_y,variable_z,
+ "residents","Cases")]
+
+
+ ################### Divide the domain of weather variable in bins of delta size for the uniform distribution.
+
+ breaks_x <-round(seq(min(na.omit(Env_Laboratory_data_df[[variable_x]])-delta_x),
+ max(na.omit(Env_Laboratory_data_df[[variable_x]])+delta_x) , by=delta_x), round_digit_x)
+
+ breaks_y <-round(seq(min(na.omit(Env_Laboratory_data_df[[variable_y]])-delta_y),
+ max(na.omit(Env_Laboratory_data_df[[variable_y]])+delta_y) , by=delta_y), round_digit_y)
+
+
+
+
+ #################### Create a data frame to fill with
+ Conditional_incidence_unif <-data.frame(character(), character(), numeric(), numeric())
+ colnames(Conditional_incidence_unif) <-c(variable_x, variable_y,"counts","residents_tot")
+
+
+
+ Env_Pathogen_data_z <-Env_Pathogen_data_df
+
+ Env_Laboratory_data_z <-Env_Laboratory_data_df
+
+ j_y_min <- 1
+ j_y_max <-length(breaks_y)
+
+ for (j_y in c(j_y_min:j_y_max)){ # from the bins resulting above, find another variable stratification
+
+ wt_y <- match.closest (Env_Pathogen_data_z[[ variable_y]], breaks_y)
+ ww_y <-which(wt_y==j_y)
+ Env_Pathogen_data_y <-Env_Pathogen_data_z[ww_y, ]
+
+
+ wt_y_t <- match.closest (Env_Laboratory_data_z[[ variable_y]], breaks_y)
+ ww_y_t <-which(wt_y_t==j_y)
+ Env_Laboratory_data_y <-Env_Laboratory_data_z[ww_y_t, ]
+
+
+ #if (length(breaks(variable_x, delta_temp))!=0)
+ {
+ j_x_min <-1
+ j_x_max <- length( breaks_x)
+
+ for (j_x in c(j_x_min:j_x_max)){
+
+ wt_x <- match.closest (Env_Pathogen_data_y[[variable_x]], breaks_x)
+ ww_x <-which(wt_x==j_x)
+ Pathogen_stratified <-Env_Pathogen_data_y[ww_x,] #where there are cases for specific values of the 3 weather variables. residents is the number of people exposed to these situations during the entire duration of the dataset. It is the number of times the people is exposed to this caracteristics, only when a case is found. It can count a same catchment area twice. Residents total is the same irrespectively of the presence of a case
+
+
+ wt_x_t <- match.closest (Env_Laboratory_data_y[[variable_x]], breaks_x)
+ ww_x_t <-which(wt_x_t==j_x)
+ Laboratory_stratified <-Env_Laboratory_data_y[ww_x_t, ] #NA in residents at Env_Laboratory_data. we do not care if there is a case reported
+
+ Total_cases <-sum((as.numeric(na.omit(Pathogen_stratified$Cases))))
+ residents_tot <-sum((as.numeric(na.omit(Laboratory_stratified$residents)))) #number of people exposed to these weather conditions. Will be used to calculate the incidence
+
+ data_df <-data.frame (
+ breaks_x[j_x],
+ breaks_y[j_y],
+ Total_cases,
+ residents_tot
+ )
+
+ colnames(data_df) <-c(variable_x, variable_y,"counts","residents_tot") #
+ Conditional_incidence_unif <-rbind(Conditional_incidence_unif, data_df)
+ }
+ }
+ }
+
+
+ print(c(j_x, j_y,variable_x,variable_y)) # check in the output log file if x y and z cover all the conditions established in breaks()
+
+
+
+ write.table(Conditional_incidence_unif, paste("../Data_Base/Conditional_incidence_",variable_x,"_",variable_y,"_",time_lag_char,info,sep=""), row.names = FALSE , sep = ",",eol = "\n")
+ return(list(Conditional_incidence_unif,breaks_x,breaks_y))
+
+ })
+}
diff --git a/Codes/Function_Plot_Conditional_Incidence_Uniform.R b/Codes/Function_Plot_Conditional_Incidence_Uniform.R
index 22a3251..47bfcd8 100644
--- a/Codes/Function_Plot_Conditional_Incidence_Uniform.R
+++ b/Codes/Function_Plot_Conditional_Incidence_Uniform.R
@@ -28,10 +28,8 @@ Plot_Conditional_incidence_Uniform<-function(Data_frame_1, Pars)
Conditional_incidence_unif<-subset(Conditional_incidence_unif, Conditional_incidence_unif$counts>100)
- number_days<-as.Date(last_day)-as.Date(first_day)
- number_days<-as.numeric(number_days)
- number_year<-number_days/365
- norm<-(number_year)/(1e6)
+
+ norm<-1/(1e6)
@@ -105,11 +103,200 @@ Plot_Conditional_incidence_Uniform<-function(Data_frame_1, Pars)
+Plot_Conditional_incidence_Uniform_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]
+ time_lag_char<-Pars[3]
+
+ # Select the file with relevant information:
+ Conditional_incidence_unif <- Data_frame_1[[1]]
+
+
+ #Conditional_incidence_unif<-Conditional_incidence_unif[,-1]
+ Conditional_incidence_unif$incidence<-Conditional_incidence_unif$counts/Conditional_incidence_unif$residents_tot
+ Conditional_incidence_unif<-na.omit(Conditional_incidence_unif)
+ #For visual purposes, we are removing cases when the counts were too low resulting in exceptional large incidence
+ Conditional_incidence_unif<-subset(Conditional_incidence_unif, Conditional_incidence_unif$counts>100)
+
+
+
+ norm<-1/(1e6)
+
+
+
+ Conditional_incidence_unif$preval<-Conditional_incidence_unif$incidence/norm
+
+
+
+ ################### Choose a specific day-length ###################
+
+
+ #Conditional_incidence_unif<-na.omit(Conditional_incidence_unif)
+ ############################# 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)}
+ if (variable_x=="Difference_temperature"){variable_x_char<-expression("Difference Maximum-Minimum Air Temperature"*~degree*C)}
+ if (variable_y=="Difference_temperature"){variable_y_char<-expression("Difference Maximum-Minimum Air Temperature"*~degree*C)}
+
+
+ Conditional_incidence_unif[,c(variable_y)]<-round(Conditional_incidence_unif[,c(variable_y)])
+
+
+ var_x_loc_plot <- ggplot(Conditional_incidence_unif,aes(x=.data[[variable_x]],y=preval,color=(.data[[variable_y]]), fill=(.data[[variable_y]])))
+
+ var_x_loc_plot <- var_x_loc_plot+ scale_fill_continuous (name =variable_y_char,type = "viridis")
+ var_x_loc_plot <- var_x_loc_plot+ scale_colour_continuous(name =variable_y_char,type = "viridis")
+
+
+ var_x_loc_plot <- var_x_loc_plot+ylab("Conditional incidence of Campylobacter Cases (per year, per million)")
+ var_x_loc_plot <- var_x_loc_plot+xlab(variable_x_char)
+
+ var_x_loc_plot <- var_x_loc_plot+ geom_point(size=3.5)
+ var_x_loc_plot <- var_x_loc_plot+ geom_line(size=1)
+
+ var_x_loc_plot <- var_x_loc_plot+ scale_x_continuous(limits=c(min(Conditional_incidence_unif[,c(variable_x)]),max(Conditional_incidence_unif[,c(variable_x)])))
+
+ var_x_loc_plot <- var_x_loc_plot+
+ theme(legend.position= c(0.125,0.75),legend.title = element_text( size = 10),
+ legend.text = element_text( size = 10),legend.background = element_blank(),
+ legend.key=element_rect(colour=NA))
+
+ var_x_loc_plot <- var_x_loc_plot+ theme(axis.title.x =element_text( colour="#990000", size=13))
+ var_x_loc_plot <- var_x_loc_plot+ theme(axis.title.y =element_text( colour="#990000", size=13))
+
+ var_x_loc_plot <- var_x_loc_plot+ theme(axis.title.x =element_text(size=13))
+ var_x_loc_plot <- var_x_loc_plot+ theme(axis.title.y =element_text(size=13))
+
+ print(var_x_loc_plot)
+
+
+
+ file_name<-paste("../Graphs/Conditional_probability_",variable_x,"_",variable_y,"_",time_lag_char,"_uniform.tiff", sep = "")
+ tiff(filename = file_name,width = 17.35, height = 17.35, units = "cm", pointsize = 9, res = 1200,compression = "lzw",antialias="default")
+
+ print(var_x_loc_plot)
+
+ dev.off()
+
+
+ })
+}
-
+Plot_Conditional_incidence_unif_constant<-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<-Pars[1]
+ info<- Pars[2]
+ # Select the file with relevant information:
+ Conditional_incidence_quantiles <- Data_frame_1
+
+ norm<-1/(1e6)
+ #Conditional_incidence_quantiles<-Conditional_incidence_quantiles[,-1]
+ Conditional_incidence_quantiles$incidence<-(Conditional_incidence_quantiles$counts/Conditional_incidence_quantiles$residents_tot)/norm
+
+
+
+
+ Conditional_incidence_df<-Conditional_incidence_quantiles
+
+ ################### Choose a specific day-length ###################
+
+
+ #Conditional_incidence_quantiles<-na.omit(Conditional_incidence_quantiles)
+ ############################# General variables ###########################
+ if (variable=="Relative_humidity"){variable_x_char<-"Relative humidity (%)"}
+ if (variable=="daylength"){variable_x_char<-"Day-Length (hours)"}
+ if (variable=="Maximum_air_temperature"){variable_x_char<-expression("Maximum Air Temperature"*~degree*C)}
+
+
+
+
+ if (variable=="Maximum_air_temperature" ){
+ Conditional_incidence_plot <-
+ ggplot(Conditional_incidence_df,aes(x=Maximum_air_temperature,y=incidence))
+ Conditional_incidence_plot <- Conditional_incidence_plot+ scale_x_continuous(limits=c(10,30))
+ Conditional_incidence_plot <- Conditional_incidence_plot+ geom_point(size=3.5, color="red")+geom_line(size=1, color="red")
+
+ }
+ if (variable=="Relative_humidity" ){
+ Conditional_incidence_plot <-
+ ggplot(Conditional_incidence_df,aes(x=Relative_humidity,y=incidence))
+ Conditional_incidence_plot <- Conditional_incidence_plot+ scale_x_continuous(limits=c(60,95))
+ Conditional_incidence_plot <- Conditional_incidence_plot+ geom_point(size=3.5, color="darkblue")+geom_line(size=1, color="darkblue")
+ }
+ if (variable=="daylength"){
+ Conditional_incidence_plot <-
+ ggplot(Conditional_incidence_df,aes(x=daylength,y=incidence))
+ Conditional_incidence_plot <- Conditional_incidence_plot+ scale_x_continuous(limits=c(10,18))
+ Conditional_incidence_plot <- Conditional_incidence_plot+ geom_point(size=3.5, color="darkgreen")+geom_line(size=1, color="darkgreen")
+ }
+
+
+
+
+ Conditional_incidence_plot <- Conditional_incidence_plot+theme_bw(20)
+ 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+ guides(colour="none") # guides(colour="none")
+ Conditional_incidence_plot <- Conditional_incidence_plot
+ Conditional_incidence_plot <- Conditional_incidence_plot
+
+ 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_",info,".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()
+
+
+
+
+ })
+}
diff --git a/Codes/Pathogen_Linkage_fixed_time_lag.R b/Codes/Pathogen_Linkage_fixed_time_lag.R
index 0d147e8..53938b2 100644
--- a/Codes/Pathogen_Linkage_fixed_time_lag.R
+++ b/Codes/Pathogen_Linkage_fixed_time_lag.R
@@ -174,7 +174,7 @@ write.table(diagnostic_laboratory_df[-c(length(diagnostic_laboratory_df))],paste
################ Define the time lag here ###############################
-time_lag<-7
+time_lag<-90
time_lag_char<-paste(time_lag)
## This create a database for all postocode where there was at least one case
--
GitLab