diff --git a/Codes/Conditional_incidence_and_reconstruction.R b/Codes/Conditional_incidence_and_reconstruction.R
index d29dac7c8c14e289f771d6f674832774acff6cc7..f80944b4e68db5b39952213dd92f697194cab266 100644
--- a/Codes/Conditional_incidence_and_reconstruction.R
+++ b/Codes/Conditional_incidence_and_reconstruction.R
@@ -46,9 +46,9 @@ All_PC<<-rep(All_PC_s,each=length(dates_s))
#Change HERE depending on the variable object of study, by copy-pasting from the list below
variable_x<-"Maximum_air_temperature"
-#variable_x<-"Relative_humidity"
-variable_y<-"Mean_Precipitation"
-#variable_x<-"Mean_wind_speed"
+#variable_y<-"Relative_humidity"
+#variable_y<-"Mean_Precipitation"
+variable_y<-"Mean_wind_speed"
variable_z<-"daylength"
#variable_z<-"Relative_humidity"
@@ -60,7 +60,7 @@ variable_z<-"daylength"
delta_variable_x<-1
round_x<-0
-delta_variable_y<-1
+delta_variable_y<-5
round_y<-0
delta_variable_z<-1
round_z<-0
diff --git a/Codes/Conditional_incidence_and_reconstruction_for_constant_situation.R b/Codes/Conditional_incidence_and_reconstruction_for_constant_situation.R
index 3d78c57853546e4c2235b77d91730763ba67ea3a..7dfdcced495dd736be09503af3d6b5dc835855f7 100644
--- a/Codes/Conditional_incidence_and_reconstruction_for_constant_situation.R
+++ b/Codes/Conditional_incidence_and_reconstruction_for_constant_situation.R
@@ -81,8 +81,8 @@ by_x<-0.1
#const_variable<-"const_temp_20_hum_76"
-#const_variable<-"const_hum_76_daylength_15"
-const_variable<-"const_temp_20_daylength_15"
+const_variable<-"const_hum_76_daylength_15"
+#const_variable<-"const_temp_20_daylength_15"
if (const_variable=="const_hum_76_daylength_15"){
diff --git a/Codes/Conditional_incidence_and_reconstruction_one_factor.R b/Codes/Conditional_incidence_and_reconstruction_one_factor.R
index 5a08a5068ce0e2ef300ce56bbd9f60a50451f7d3..c87a77a0ad7e1068276742e65ee1b9e4137b20e6 100644
--- a/Codes/Conditional_incidence_and_reconstruction_one_factor.R
+++ b/Codes/Conditional_incidence_and_reconstruction_one_factor.R
@@ -41,16 +41,16 @@ 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
-variables<-c("Maximum_air_temperature","Minimum_air_temperature","Mean_wind_speed","Mean_Precipitation","Relative_humidity","daylength")
-time_lags <- c(7,14,30,60,90)
-for (i in seq(1:5)) {
- time_lag <<-time_lags[i]
+time_lag <- 14
+
time_lag_char <<-paste(time_lag)
-for (j in seq(1:6)) {
-variable<-variables[j]
+#variable<-"daylength"
+variable<-"Mean_wind_speed"
+#variable<-"Relative_humidity"
+#variable<-"Maximum_air_temperature"
+#variable<-"Mean_Precipitation"
### Set the size of the bins for Uniform conditional incidence and resolution in rounding the values of the bins.
### For example:
@@ -58,12 +58,28 @@ variable<-variables[j]
## 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
-by_z<-0.05
+by_x<-0.05
+delta_variable_x<-1
+round_x<-0
+
+######## Input of data ####
######## Input of data ####
source("Function_Input_files.R", chdir=TRUE)
+Env_Pathogen_data$Difference_temperature<-Env_Pathogen_data$Maximum_air_temperature-Env_Pathogen_data$Minimum_air_temperature
+Env_laboratory_data$Difference_temperature<-Env_laboratory_data$Maximum_air_temperature-Env_laboratory_data$Minimum_air_temperature
+
+######## 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 ######
@@ -71,9 +87,45 @@ source("Function_Input_files.R", chdir=TRUE)
# 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_Quantiles.R", chdir=TRUE)
-Pars_cond_incidence_quantiles<- c(variable,"_Quantiles",by_z)
-Conditional_incidence_quantiles_info<-Conditional_incidence_Quantiles_one_factor(Env_laboratory_data,Env_Pathogen_data,Pars_cond_incidence_quantiles)
-}}
-# use Plots_Campylobacter_environment_one_variables_quantile.R to plot
+#
+source("Function_Conditional_Incidence_Uniform.R", chdir=TRUE)
+#Conditional incidence, Uniform bins
+Pars_cond_incidence <- c(variable,"_Uniform.csv",delta_variable_x,round_x)
+Conditional_incidence_unif_info<-Conditional_incidence_Uniform_one_factor(Env_laboratory_data,Env_Pathogen_data,Pars_cond_incidence)
+
+## plot Conditional incidence for uniform, for diagnostic only
+Pars_plot_cond_prev_uniform <- c(variable,time_lag_char)
+source("Function_Plot_Conditional_incidence_Uniform.R", chdir=TRUE)
+Plot_Conditional_incidence_Uniform_one_factor(Conditional_incidence_unif_info,Pars_plot_cond_prev_uniform)
+
+
+
+
+#source("Function_Conditional_Incidence_Quantiles.R", chdir=TRUE)
+##Conditional incidence, bins divided in quantiles. As by_z=1, there is only one bin for the third variable, so practically we consider two variables
+#Pars_cond_prev_quantiles<- c(variable,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)
+
+## plot Conditional incidence for bins divided in quantiles
+#Pars_plot_cond_prev_quantiles <- c(variable_x,variable_y)
+#source("Function_Plot_Conditional_incidence_Quantiles.R", chdir=TRUE)
+#Plot_Conditional_incidence_quantiles_two_factors(Conditional_incidence_quantiles_info,Pars_plot_cond_prev_quantiles)
+
+
+###### Reconstruction ######
+
+source("Function_Reconstruction.R")
+Pars_reconstruction <- c(variable,"_reconstructed.csv")
+Reconstruction_time_series<-Reconstruction_time_series_one_factor(Conditional_incidence_unif_info,Env_laboratory_data,Pars_reconstruction)
+
+
+
+############## End reconstructions ####################
+
+###### Plot Reconstruction ######
+source("Function_Plot_Reconstruction.R")
+Pars_Plot_Reconstruction <- c(paste("../Graphs/",variable,time_lag_char,"_time_series.tiff",sep=""),paste("../Graphs/",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 a49c5edea4020a5849278840c3b673556dcacd56..6682d3f0b4b24326a349a1a0c2574775ba6154dd 100644
--- a/Codes/Conditional_incidence_and_reconstruction_two_factors.R
+++ b/Codes/Conditional_incidence_and_reconstruction_two_factors.R
@@ -22,7 +22,7 @@ library(zoo)
# Set length of the desired time-lag
-time_lag <<-14
+time_lag <<-30
time_lag_char <<-paste(time_lag)
# Set the years you are focusing on
@@ -47,12 +47,14 @@ All_PC<<-rep(All_PC_s,each=length(dates_s))
#variable_x<-"Difference_temperature"
#variable_x<-"Minimum_air_temperature"
-#variable_y<-"Mean_Precipitation"
-variable_y<-"Relative_humidity"
+#variable_x<-"Mean_Precipitation"
+#variable_y<-"Relative_humidity"
#variable_y<-"Mean_wind_speed"
-#variable_y<-"daylength"
-
variable_z<-"daylength"
+
+#variable_z<-"daylength"
+#variable_y<-"Mean_Precipitation"
+variable_y<-"Relative_humidity"
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:
@@ -62,8 +64,8 @@ variable_x<-"Maximum_air_temperature"
delta_variable_x<-1
round_x<-0
-delta_variable_y<-0.2
-round_y<-1
+delta_variable_y<-1
+round_y<-0
## partitions of the quantiles
by_z<-1
by_y<-0.25
diff --git a/Codes/Conditional_incidence_one_factor.R b/Codes/Conditional_incidence_one_factor.R
new file mode 100644
index 0000000000000000000000000000000000000000..c2fbc9268e2f44b610fb9105a10ec90ecde74585
--- /dev/null
+++ b/Codes/Conditional_incidence_one_factor.R
@@ -0,0 +1,80 @@
+# 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
+
+
+# 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
+variables<-c("Maximum_air_temperature","Minimum_air_temperature","Mean_wind_speed","Mean_Precipitation","Relative_humidity","daylength")
+time_lags <- c(7,14,30,60,90)
+
+for (i in seq(1:5)) {
+ time_lag <<-time_lags[i]
+ time_lag_char <<-paste(time_lag)
+
+for (j in seq(1:6)) {
+variable<-variables[j]
+
+### 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
+
+by_x<-0.05
+delta_variable_x<-0.2
+round_x<-1
+
+
+######## Input of data ####
+source("Function_Input_files.R", chdir=TRUE)
+
+
+###### 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_Quantiles.R", chdir=TRUE)
+Pars_cond_incidence_quantiles<- c(variable,"_Quantiles",by_x)
+Conditional_incidence_quantiles_info<-Conditional_incidence_Quantiles_one_factor(Env_laboratory_data,Env_Pathogen_data,Pars_cond_incidence_quantiles)
+}}
+# use Plots_Campylobacter_environment_one_variables_quantile.R to plot
+
diff --git a/Codes/Function_Conditional_Incidence_Uniform.R b/Codes/Function_Conditional_Incidence_Uniform.R
index f645ee2a1aff5fa605817268decb1014038e328c..a7eeefbd8c1fecd456dada5c21cc0127b106c483 100644
--- a/Codes/Function_Conditional_Incidence_Uniform.R
+++ b/Codes/Function_Conditional_Incidence_Uniform.R
@@ -238,3 +238,70 @@ Conditional_incidence_Uniform_two_factors<-function(Data_frame_1, Data_frame_2,P
})
}
+
+
+
+
+Conditional_incidence_Uniform_one_factor<-function(Data_frame_1, Data_frame_2,Pars)
+{
+
+ with(as.list(c(Data_frame_1, Data_frame_2,Pars)), {
+
+ variable_x<-Pars[1]
+ info<-Pars[2]
+ delta_x<-as.numeric(Pars[3])
+ round_digit_x<-as.numeric(Pars[4])
+
+ # Select the variables of study:
+ Env_Laboratory_data_df <- Data_frame_1[,c("PostCode","Date", variable_x,
+ "residents")]
+ Env_Pathogen_data_df <- Data_frame_2[,c("PostCode","Date", variable_x,
+ "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)
+
+
+
+ #################### Create a data frame to fill with
+ Conditional_incidence_unif <-data.frame(character(), numeric(), numeric())
+ colnames(Conditional_incidence_unif) <-c(variable_x, "counts","residents_tot")
+
+
+ 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_df[[variable_x]], breaks_x)
+ ww_x <-which(wt_x==j_x)
+ Pathogen_stratified <-Env_Pathogen_data_df[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_df[[variable_x]], breaks_x)
+ ww_x_t <-which(wt_x_t==j_x)
+ Laboratory_stratified <-Env_Laboratory_data_df[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],
+ Total_cases,
+ residents_tot
+ )
+ colnames(data_df) <-c(variable_x, "counts","residents_tot") #
+ Conditional_incidence_unif <-rbind(Conditional_incidence_unif, data_df)
+
+ print(c(j_x, variable_x)) # 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,"_",time_lag_char,info,sep=""), row.names = FALSE , sep = ",",eol = "\n")
+
+ return(list(Conditional_incidence_unif,breaks_x))
+
+ })
+}
diff --git a/Codes/Function_Detrend_analysis.R b/Codes/Function_Detrend_analysis.R
new file mode 100644
index 0000000000000000000000000000000000000000..f7eb8b2fb6d7c1d99508f423482ec7a04ca321cf
--- /dev/null
+++ b/Codes/Function_Detrend_analysis.R
@@ -0,0 +1,51 @@
+# Remove temporal trend from the time-series of salmonellosis
+
+Detrend_analysis<-function(Data_frame_1, Data_frame_2,Pars)
+{
+
+ with(as.list(c(Data_frame_1, Data_frame_2,Pars)), {
+
+
+ pathogen_data_national<-Data_frame_1
+ pathogen_enviroment_data<-Data_frame_2
+ file_name<-Pars
+
+mean_cases <- mean(pathogen_data_national$Cases)
+national_detrended <-detrend.series(y=pathogen_data_national$Cases, y.name = "", make.plot = TRUE,
+ method = c("Spline"),return.info=TRUE)
+
+pathogen_df <-data.frame(pathogen_data_national$Date, pathogen_data_national$Cases, national_detrended$series, national_detrended$curves)
+colnames(pathogen_df)<-c("Date","cases","detrend","curve")
+
+
+pathogen_cases_detrend <-merge(pathogen_enviroment_data , pathogen_df, by="Date")
+# After discussing with GLI, we concluded that applying the national curve data to a smaller PC area is similar to calculating the curve of the specific PC.
+pathogen_cases_detrend$Cases_detrend<-pathogen_cases_detrend$Cases/pathogen_cases_detrend$curve # correct to use Cases here
+pathogen_cases_detrend$Cases_detrend_adjusted <- pathogen_cases_detrend$Cases_detrend* mean_cases # mean daily cases
+
+# Plot detrend:
+pl3<-ggplot(pathogen_cases_detrend, aes(x=Date,y=Cases_detrend_adjusted))+geom_line(color="darkred",stat = "identity")
+
+pl3<- pl3 +xlab("Date")+ylab("Number of Occurences")#+scale_x_date(labels = date_format("%d/%m/%y"),breaks = pretty_breaks())
+
+pl3<- pl3+ theme(axis.title.y =element_text(face="bold", colour="#990000", size=13))
+
+pl3<- pl3+theme(axis.title.x =element_text(face="bold", colour="#990000", size=13)) + ggtitle ("Detrended cases or individual postocodes")
+
+pl3
+
+
+tiff(filename = file_name ,width = 17.35, height = 17.35, units = "cm", pointsize = 9, res = 600,compression = "lzw",antialias="default")
+
+pl3
+
+dev.off()
+
+####### Diagnostic plot ##########
+time_series <-ddply(pathogen_cases_detrend,~Date, cases=sum(Cases_detrend_adjusted))
+plot(time_series$Date,time_series$cases,type="l",main="AFTER DETREND",xlab="Date",ylab="Cases at National Level")
+
+return(pathogen_cases_detrend)
+
+})
+}
diff --git a/Codes/Function_Input_files.R b/Codes/Function_Input_files.R
index abb55b26236580e86149d9fb58789f1cb36097bb..e27c3e311553a4d565038209141cdabf77a7b138 100644
--- a/Codes/Function_Input_files.R
+++ b/Codes/Function_Input_files.R
@@ -20,8 +20,13 @@ Env_Pathogen_data_all$Temperature_amplitude <- Env_Pathogen_data_all$Maximum_air
# Select the variables of interest:
Env_Pathogen_data_all$Date<-as.Date(Env_Pathogen_data_all$Date,format = "%d/%m/%Y")
Env_Pathogen_data_all <- distinct(Env_Pathogen_data_all) # Remove possible duplicates Gianni, correct but better to do this when we create the file
+Env_Pathogen_data_rolling_mean <-Env_Pathogen_data_all %>% mutate (rolling_mean= slider::slide_mean(Cases, before=7 , after=0 )) # %>% ungroup()
+Env_Pathogen_data_all$Cases <-Env_Pathogen_data_rolling_mean$rolling_mean
Env_Pathogen_data <-(subset(Env_Pathogen_data_all, year(Env_Pathogen_data_all$Date)>=year_first & year(Env_Pathogen_data_all$Date)<=year_last))
+
+
+
# Read the weather variables for every PC, averaged for the estimated delay in the past
Env_laboratory_data_all <-read_csv(paste("../Data_Base/Laboratory_environment_",time_lag_char,".csv",sep=""))
diff --git a/Codes/Function_Plot_Conditional_Incidence_Quantiles.R b/Codes/Function_Plot_Conditional_Incidence_Quantiles.R
index d94b97c4816b9183494666b8c6b460b079d7fc06..8133f4e5fd5b39e7c688f6c2e150491be888a04c 100644
--- a/Codes/Function_Plot_Conditional_Incidence_Quantiles.R
+++ b/Codes/Function_Plot_Conditional_Incidence_Quantiles.R
@@ -94,7 +94,7 @@ Plot_Conditional_incidence_quantiles<-function(Data_frame_1, Pars)
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))
+ 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 <-
@@ -104,7 +104,7 @@ Plot_Conditional_incidence_quantiles<-function(Data_frame_1, Pars)
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))
+ 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 <-
@@ -116,7 +116,7 @@ Plot_Conditional_incidence_quantiles<-function(Data_frame_1, Pars)
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))
+ 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 <-
@@ -127,7 +127,7 @@ Plot_Conditional_incidence_quantiles<-function(Data_frame_1, Pars)
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+ scale_y_continuous(limits=c(0,5))+scale_x_continuous(limits=c(0,25))
}
@@ -218,7 +218,7 @@ Plot_Conditional_incidence_quantiles<-function(Data_frame_1, Pars)
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))
+ 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 <-
@@ -228,7 +228,7 @@ Plot_Conditional_incidence_quantiles<-function(Data_frame_1, Pars)
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))
+ 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 <-
@@ -240,7 +240,7 @@ Plot_Conditional_incidence_quantiles<-function(Data_frame_1, Pars)
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))
+ 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 <-
@@ -251,7 +251,7 @@ Plot_Conditional_incidence_quantiles<-function(Data_frame_1, Pars)
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+ scale_y_continuous(limits=c(0,5))+scale_x_continuous(limits=c(0,25))
}
if (variable_x=="daylength" & variable_y=="Maximum_air_temperature"){
Conditional_incidence_plot <-
@@ -396,25 +396,34 @@ Plot_Conditional_incidence_quantiles_two_factors<-function(Data_frame_1, Pars)
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))
+ 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(-5,30))
+ 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(-5,30))
+ 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(-5,30))
+ 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 <-
diff --git a/Codes/Function_Plot_Conditional_Incidence_Uniform.R b/Codes/Function_Plot_Conditional_Incidence_Uniform.R
index 47bfcd8dbf333a0c71a251ea8b847b79e1851797..51cc7139f70efc443d3e85a3aa538dd7de8da5fb 100644
--- a/Codes/Function_Plot_Conditional_Incidence_Uniform.R
+++ b/Codes/Function_Plot_Conditional_Incidence_Uniform.R
@@ -204,6 +204,95 @@ Plot_Conditional_incidence_Uniform_two_factors<-function(Data_frame_1, Pars)
+
+Plot_Conditional_incidence_Uniform_one_factor<-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]
+ time_lag_char<-Pars[2]
+
+ # 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_x=="Mean_wind_speed"){variable_x_char<-"Mean Wind Speed (knots)"}
+ if (variable_x=="Mean_wind_speed"){variable_x_char<-"Mean Wind Speed (knots)"}
+ if (variable_x=="Mean_Precipitation"){variable_x_char<-"Mean Precipitation (mm)"}
+ if (variable_x=="daylength"){variable_x_char<-"Day-Length (hours)"}
+ if (variable_x=="Minimum_air_temperature"){variable_x_char<-expression("Minimum Air Temperature"*~degree*C)}
+ if (variable_x=="Maximum_air_temperature"){variable_x_char<-expression("Maximum Air Temperature"*~degree*C)}
+ if (variable_x=="Difference_temperature"){variable_x_char<-expression("Difference Maximum-Minimum Air Temperature"*~degree*C)}
+
+
+
+ var_x_loc_plot <- ggplot(Conditional_incidence_unif,aes(x=.data[[variable_x]],y=preval))
+
+ var_x_loc_plot <- var_x_loc_plot+ scale_fill_continuous (type = "viridis")
+ var_x_loc_plot <- var_x_loc_plot+ scale_colour_continuous(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,"_",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)
{
diff --git a/Codes/Function_Plot_Reconstruction.R b/Codes/Function_Plot_Reconstruction.R
index a8d8a7226e1bb2de4e6efa780aedbb2033b52720..38aefd66081647ff4d369a44ba0a9ea9bfe99904 100644
--- a/Codes/Function_Plot_Reconstruction.R
+++ b/Codes/Function_Plot_Reconstruction.R
@@ -14,14 +14,16 @@ Plot_Reconstruction_time_series<-function(Data_frame_1, Data_frame_2,Pars)
colnames(time_series_Post_Codes) <-c("Date","Cases","Lambda","PostCode")
time_series<-aggregate (Cases ~ Date, time_series_Post_Codes, sum) # get the mean of the cases per day for all the PC
+
time_series_quantile <-ddply(time_series,~Date, function (x) quantile(x$Cases, c(.25,.5,.75)))
time_series_summary <-data.frame(time_series_quantile, rep("Model",times=length(time_series[,1])))
colnames(time_series_summary)<-c("Date","f_quant","Cases","s_quant","source")
-#real_cases_mean <-empirical_time_series %>% mutate (rolling_mean= slider::slide_mean(Cases, before=3 , after=3 )) # %>% ungroup() #Gianni Check this, I changed as not sure it worked
-#empirical_time_series$rolling_mean <-real_cases_mean$rolling_mean
+real_cases_mean <-empirical_time_series %>% mutate (rolling_mean= slider::slide_mean(Cases, before=7 , after=0 )) # %>% ungroup() ##rolling mean to aremove the weekend effects
+empirical_time_series$rolling_mean <-real_cases_mean$rolling_mean
+real_cases_quantile<-ddply(empirical_time_series,~Date, function (x) quantile(x$rolling_mean, c(.25,.5,.75)))
-real_cases_quantile<-ddply(empirical_time_series,~Date, function (x) quantile(x$Cases, c(.25,.5,.75)))
+#real_cases_quantile<-ddply(empirical_time_series,~Date, function (x) quantile(x$Cases, c(.25,.5,.75)))
real_cases_summary<-cbind(real_cases_quantile,rep("Cases",times=length(real_cases_quantile[,1])))
colnames(real_cases_summary)<-c("Date","f_quant","Cases","s_quant","source")
@@ -34,13 +36,14 @@ time_series_all <-rbind(time_series_summary, real_cases_summary)
pal<-wes_palette("Cavalcanti1")
time_series_plot<-ggplot(time_series_all,aes(x=Date,y=Cases,colour=source))
+time_series_plot<-time_series_plot+ theme_bw(20)
#time_series_plot<-time_series_plot+scale_color_brewer(palette=pal)
#time_series_plot<-time_series_plot+geom_line(data=real_cases_summary)
time_series_plot<-time_series_plot+geom_line(size=0.75)
#time_series_plot<-time_series_plot+scale_color_manual(values=c( "#E69F00", "#56B4E9"))
time_series_plot<-time_series_plot+scale_color_manual(values=wes_palette(n=2, name="Cavalcanti1"))
-time_series_plot<-time_series_plot+ xlab("Date")+ ylab("Campylobacter Cases")
+time_series_plot<-time_series_plot+ xlab("Date")+ ylab("Campylobacteriosis")
time_series_plot<-time_series_plot+
theme(legend.position= c(0.125,0.85),legend.title = element_blank(),
@@ -96,7 +99,7 @@ tick_pos<-yday(as.Date(as.character(
-yearly_average <-ggplot(average_data, aes(x=Day, y=Mean, colour=source))+
+yearly_average <-ggplot(average_data, aes(x=Day, y=Mean, colour=source))+ theme_bw(20)+
geom_line(size=2)+
geom_ribbon(aes(ymin=f_quant, ymax=s_quant, fill=source), alpha=0.15) +
xlab("Day of the Year") + ylab(y_lab_text) +
@@ -126,4 +129,3 @@ return()
})
}
-
diff --git a/Codes/Function_Reconstruction.R b/Codes/Function_Reconstruction.R
index cc7aab272c8ab810aa29e82df8399692e3111a57..cd78ec184aa26fe883abb0a8385db111bf710a68 100644
--- a/Codes/Function_Reconstruction.R
+++ b/Codes/Function_Reconstruction.R
@@ -1,13 +1,13 @@
-Reconstruction_time_series<-function(Conditional_prevalence_unif_info, Data_frame_2,Pars)
+Reconstruction_time_series<-function(Conditional_incidence_unif_info, Data_frame_2,Pars)
{
- with(as.list(c(Conditional_prevalence_unif_info, Data_frame_2,Pars)), {
+ with(as.list(c(Conditional_incidence_unif_info, Data_frame_2,Pars)), {
-
- Conditional_prevalence_unif<-Conditional_prevalence_unif_info[[1]]
- breaks_x<-Conditional_prevalence_unif_info[[2]]
- breaks_y<-Conditional_prevalence_unif_info[[3]]
- breaks_z<-Conditional_prevalence_unif_info[[4]]
+
+ Conditional_incidence_unif<-Conditional_incidence_unif_info[[1]]
+ breaks_x<-Conditional_incidence_unif_info[[2]]
+ breaks_y<-Conditional_incidence_unif_info[[3]]
+ breaks_z<-Conditional_incidence_unif_info[[4]]
Env_laboratory_data<-Data_frame_2
variable_x<-Pars[1]
variable_y<-Pars[2]
@@ -46,7 +46,7 @@ Reconstruction_time_series<-function(Conditional_prevalence_unif_info, Data_fram
colnames(variable_df2) <- c("PostCode","Date", "residents", variable_z, variable_y, variable_x)
dt_df <- data.table(variable_df2)
- dt_var <- data.table(Conditional_prevalence_unif)
+ dt_var <- data.table(Conditional_incidence_unif)
variable_df_dis3 <- dt_df[dt_var, on= c(variable_z,variable_y,variable_x), allow.cartesian=TRUE] # repeated rows since same conditions apply to more than one PC. Gianni this shoul be a merge
variable_df_dis3 <- na.omit(variable_df_dis3) %>% distinct()
@@ -57,7 +57,7 @@ Reconstruction_time_series<-function(Conditional_prevalence_unif_info, Data_fram
# calculation of all the potential incidence based on the weather conditions of the area. The counts used as numerator correspond to the simulation of cases per weather condition.
- variable_df_dis3$prevalence <-variable_df_dis3$counts/variable_df_dis3$residents_tot
+ variable_df_dis3$incidence <-variable_df_dis3$counts/variable_df_dis3$residents_tot
#variable_df_dis_excluded<-subset(variable_df_dis3,variable_df_dis3$residents_tot<=cutoff)
#print(variable_df_dis_excluded)
#write.csv(na.omit(variable_df_dis_excluded),"variable_df_dis_excluded.csv")
@@ -65,8 +65,8 @@ Reconstruction_time_series<-function(Conditional_prevalence_unif_info, Data_fram
# if (length(na.omit(variable_df_dis3[,1])!=0)){
- comp_cases <-variable_df_dis3$prevalence*variable_df_dis3$residents # estimation of the "real" incidence based on the possible cases due to weather and the residents number in the catchment area
- time_series <-data.frame(variable_df_dis3$Date, comp_cases, variable_df_dis3$prevalence, variable_df_dis3$PostCode)
+ comp_cases <-variable_df_dis3$incidence*variable_df_dis3$residents # estimation of the "real" incidence based on the possible cases due to weather and the residents number in the catchment area
+ time_series <-data.frame(variable_df_dis3$Date, comp_cases, variable_df_dis3$incidence, variable_df_dis3$PostCode)
colnames(time_series) <-c("Date","Cases","Lambda","Lab")
#}
@@ -74,6 +74,142 @@ Reconstruction_time_series<-function(Conditional_prevalence_unif_info, Data_fram
write.table(time_series, paste("../Data_Base/",variable_x,"_",variable_y,"_",variable_z,"_",time_lag_char,info,sep=""), col.names= NA, row.names= TRUE, sep=",",eol= "\n") # Saved on 17/11 after switching residents_total and residents from the denominator respecting the above attempt.
return(time_series)
+
+ })
+}
+
+
+Reconstruction_time_series_two_factors<-function(Conditional_incidence_unif_info, Data_frame_2,Pars)
+{
+
+ with(as.list(c(Conditional_incidence_unif_info, Data_frame_2,Pars)), {
+
+
+ Conditional_incidence_unif<-Conditional_incidence_unif_info[[1]]
+ breaks_x<-Conditional_incidence_unif_info[[2]]
+ breaks_y<-Conditional_incidence_unif_info[[3]]
+ Env_laboratory_data<-Data_frame_2
+ variable_x<-Pars[1]
+ variable_y<-Pars[2]
+ info<-Pars[3]
+
+
+ variable_df <- as.data.frame(Env_laboratory_data[,c("PostCode","Date", "residents")])
+
+ # y variable
+ dt_y = data.table(breaks_y, val = breaks_y)
+ dt_y <- na.omit (dt_y)
+ setattr(dt_y, "sorted", "breaks_y")
+ dt_y2 <- dt_y[J(Env_laboratory_data[[variable_y]]), roll = "nearest"]
+
+
+ variable_df$variable_y_break<- dt_y2$val
+
+
+
+ # x variable
+ dt_temp = data.table(breaks_x, val = breaks_x)
+ dt_temp <- na.omit (dt_temp)
+ setattr(dt_temp, "sorted", "breaks_x")
+ dt_temp2 <- dt_temp[J(Env_laboratory_data[[variable_x]]), roll = "nearest"]
+
+ variable_df$variable_x_break<- dt_temp2$val
+
+ variable_df2<-data.frame(variable_df$PostCode,variable_df$Date,variable_df$residents, variable_df$variable_y_break, variable_df$variable_x_break)
+ colnames(variable_df2) <- c("PostCode","Date", "residents", variable_y, variable_x)
+
+ dt_df <- data.table(variable_df2)
+ dt_var <- data.table(Conditional_incidence_unif)
+ variable_df_dis3 <- dt_df[dt_var, on= c(variable_y,variable_x), allow.cartesian=TRUE] # repeated rows since same conditions apply to more than one PC. Gianni this shoul be a merge
+
+ variable_df_dis3 <- na.omit(variable_df_dis3) %>% distinct()
+
+
+
+ variable_df_dis3 <-variable_df_dis3[order(as.Date(variable_df_dis3$Date)),]
+
+
+ # calculation of all the potential incidence based on the weather conditions of the area. The counts used as numerator correspond to the simulation of cases per weather condition.
+ variable_df_dis3$incidence <-variable_df_dis3$counts/variable_df_dis3$residents_tot
+ #variable_df_dis_excluded<-subset(variable_df_dis3,variable_df_dis3$residents_tot<=cutoff)
+ #print(variable_df_dis_excluded)
+ #write.csv(na.omit(variable_df_dis_excluded),"variable_df_dis_excluded.csv")
+ #variable_df_dis3<- subset(variable_df_dis3,variable_df_dis3$residents_tot>=cutoff) #Gianni discuss this
+
+ # if (length(na.omit(variable_df_dis3[,1])!=0)){
+
+ comp_cases <-variable_df_dis3$incidence*variable_df_dis3$residents # estimation of the "real" incidence based on the possible cases due to weather and the residents number in the catchment area
+ time_series <-data.frame(variable_df_dis3$Date, comp_cases, variable_df_dis3$incidence, variable_df_dis3$PostCode)
+ colnames(time_series) <-c("Date","Cases","Lambda","Lab")
+
+ #}
+
+ write.table(time_series, paste("../Data_Base/",variable_x,"_",variable_y,"_",time_lag_char,info,sep=""), col.names= NA, row.names= TRUE, sep=",",eol= "\n") # Saved on 17/11 after switching residents_total and residents from the denominator respecting the above attempt.
+ return(time_series)
+
+
+ })
+}
-})
+
+
+
+Reconstruction_time_series_one_factor<-function(Conditional_incidence_unif_info, Data_frame_2,Pars)
+{
+
+ with(as.list(c(Conditional_incidence_unif_info, Data_frame_2,Pars)), {
+
+
+ Conditional_incidence_unif<-Conditional_incidence_unif_info[[1]]
+ breaks_x<-Conditional_incidence_unif_info[[2]]
+ Env_laboratory_data<-Data_frame_2
+ variable_x<-Pars[1]
+ info<-Pars[2]
+
+
+ variable_df <- as.data.frame(Env_laboratory_data[,c("PostCode","Date", "residents")])
+
+
+ # x variable
+ dt_temp = data.table(breaks_x, val = breaks_x)
+ dt_temp <- na.omit (dt_temp)
+ setattr(dt_temp, "sorted", "breaks_x")
+ dt_temp2 <- dt_temp[J(Env_laboratory_data[[variable_x]]), roll = "nearest"]
+
+ variable_df$variable_x_break<- dt_temp2$val
+
+ variable_df2<-data.frame(variable_df$PostCode,variable_df$Date,variable_df$residents, variable_df$variable_x_break)
+ colnames(variable_df2) <- c("PostCode","Date", "residents", variable_x)
+
+ dt_df <- data.table(variable_df2)
+ dt_var <- data.table(Conditional_incidence_unif)
+ variable_df_dis3 <- dt_df[dt_var, on= c(variable_x), allow.cartesian=TRUE] # repeated rows since same conditions apply to more than one PC. Gianni this shoul be a merge
+
+ variable_df_dis3 <- na.omit(variable_df_dis3) %>% distinct()
+
+
+
+ variable_df_dis3 <-variable_df_dis3[order(as.Date(variable_df_dis3$Date)),]
+
+
+ # calculation of all the potential incidence based on the weather conditions of the area. The counts used as numerator correspond to the simulation of cases per weather condition.
+ variable_df_dis3$incidence <-variable_df_dis3$counts/variable_df_dis3$residents_tot
+ #variable_df_dis_excluded<-subset(variable_df_dis3,variable_df_dis3$residents_tot<=cutoff)
+ #print(variable_df_dis_excluded)
+ #write.csv(na.omit(variable_df_dis_excluded),"variable_df_dis_excluded.csv")
+ #variable_df_dis3<- subset(variable_df_dis3,variable_df_dis3$residents_tot>=cutoff) #Gianni discuss this
+
+ # if (length(na.omit(variable_df_dis3[,1])!=0)){
+
+ comp_cases <-variable_df_dis3$incidence*variable_df_dis3$residents # estimation of the "real" incidence based on the possible cases due to weather and the residents number in the catchment area
+ time_series <-data.frame(variable_df_dis3$Date, comp_cases, variable_df_dis3$incidence, variable_df_dis3$PostCode)
+ colnames(time_series) <-c("Date","Cases","Lambda","Lab")
+
+ #}
+
+ write.table(time_series, paste("../Data_Base/",variable_x,"_",time_lag_char,info,sep=""), col.names= NA, row.names= TRUE, sep=",",eol= "\n") # Saved on 17/11 after switching residents_total and residents from the denominator respecting the above attempt.
+ return(time_series)
+
+
+ })
}
diff --git a/Codes/Function_Wavelet_analysis.R b/Codes/Function_Wavelet_analysis.R
index df24444eca9a42882d518a36b95cc717d4e46095..d1d1b065cd6135c40d7fb0c5b14fd1b09f95a5f4 100644
--- a/Codes/Function_Wavelet_analysis.R
+++ b/Codes/Function_Wavelet_analysis.R
@@ -9,17 +9,17 @@ with(as.list(c(State, Pars)), {
filename_power_average<-Pars[2]
-Env_Pathogen_data0 <- ddply (State, ~Date, summarise, Cases=sum(Cases_computed)) # subset for all postcodes combined per day
+Env_Pathogen_data0 <- ddply (State, ~Date, summarise, Cases=sum(Cases)) # subset for all postcodes combined per day
# Wavelet spectrum plot
Env_Pathogen_data0$date<-Env_Pathogen_data0$Date
-salmonella_data_national <-Env_Pathogen_data0[order(as.Date(Env_Pathogen_data0$Date)),] # Sort from least recent to most recent
+pathogen_data_national <-Env_Pathogen_data0[order(as.Date(Env_Pathogen_data0$Date)),] # Sort from least recent to most recent
dt0<-1
-my.wt = analyze.wavelet(salmonella_data_national, "Cases",
+my.wt = analyze.wavelet(pathogen_data_national, "Cases",
loess.span=0.75,
dt=dt0, dj=1/20,
lowerPeriod=2.*dt0,
- upperPeriod=floor(nrow(salmonella_data_national)/3)*dt0,
+ upperPeriod=floor(nrow(pathogen_data_national)/3)*dt0,
make.pval=T, n.sim=100)
@@ -54,6 +54,6 @@ wt.avg(my.wt,
dev.off()
-return(salmonella_data_national[,c(1,2)])
+return(pathogen_data_national[,c(1,2)])
})
}
diff --git a/Codes/Pathogen_Linkage_fixed_time_lag.R b/Codes/Pathogen_Linkage_fixed_time_lag.R
index 53938b2b8ecddcc18c656c7bf2ca772475dcd8dd..52ab00c2ff21aaf84b1fceaf840e6f9d0673f88e 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<-90
+time_lag<-14
time_lag_char<-paste(time_lag)
## This create a database for all postocode where there was at least one case
diff --git a/Codes/Plot_reported_and_adjusted_campylobacter_cases.R b/Codes/Plot_reported_and_adjusted_campylobacter_cases.R
index 2203f6a34b151092e09ac069fb1a505fb50e2deb..1ee4ed390c53fffe5d274ed30fb8c2b6afd08b12 100644
--- a/Codes/Plot_reported_and_adjusted_campylobacter_cases.R
+++ b/Codes/Plot_reported_and_adjusted_campylobacter_cases.R
@@ -19,6 +19,8 @@ Campylobacter_adjusted_date_df<-read.csv("../Data_Base/Campylobacter_adjusted_da
Campylobacter_adjusted_date_df$Cases<-1
Campylobacter_adjusted_date_df$Reported_Cases<-1
+
+Campylobacter_adjusted_date_df$Date<-as.Date(Campylobacter_adjusted_date_df$Date,format="%d/%m/%Y")
#This line associates one case for each day and PostCode
Campylobacter_adjusted_date_df<-subset(Campylobacter_adjusted_date_df,as.Date(as.character(Campylobacter_adjusted_date_df$Date))>="1989-01-01")