Skip to content
Snippets Groups Projects
Commit b7c3cddb authored by Lo Iacono, Giovanni Dr (School of Vet Med.)'s avatar Lo Iacono, Giovanni Dr (School of Vet Med.)
Browse files

minor corrections

parent e061b535
No related branches found
No related tags found
No related merge requests found
...@@ -118,7 +118,7 @@ Pars_cond_prev_quantiles<- c(variable_x,variable_y,variable_z,"_Quantiles.csv",b ...@@ -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) 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 ## 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) source("Function_Plot_Conditional_Incidence_Quantiles.R", chdir=TRUE)
Plot_Conditional_incidence_quantiles(Conditional_incidence_quantiles_info,Pars_plot_cond_prev_quantiles) Plot_Conditional_incidence_quantiles(Conditional_incidence_quantiles_info,Pars_plot_cond_prev_quantiles)
......
# 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 ####################
...@@ -45,7 +45,7 @@ All_PC<<-rep(All_PC_s,each=length(dates_s)) ...@@ -45,7 +45,7 @@ All_PC<<-rep(All_PC_s,each=length(dates_s))
# Selection of the variables of interest. # Selection of the variables of interest.
#Change HERE depending on the variable object of study, by copy-pasting from the list below #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_x<-"Minimum_air_temperature"
#variable_y<-"Mean_Precipitation" #variable_y<-"Mean_Precipitation"
variable_y<-"Relative_humidity" variable_y<-"Relative_humidity"
...@@ -53,7 +53,7 @@ variable_y<-"Relative_humidity" ...@@ -53,7 +53,7 @@ variable_y<-"Relative_humidity"
#variable_y<-"daylength" #variable_y<-"daylength"
variable_z<-"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. ### Set the size of the bins for Uniform conditional incidence and resolution in rounding the values of the bins.
### For example: ### 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 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 ...@@ -91,7 +91,7 @@ Env_laboratory_data$Difference_temperature<-Env_laboratory_data$Maximum_air_temp
###### helpful but not necessary. So this is commented ###### helpful but not necessary. So this is commented
source("Function_Wavelet_analysis.R", chdir=TRUE) 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 Env_Pathogen_data$Cases_computed<-Env_Pathogen_data$Cases
campylobacter_data_national<-wavelet_analysis(Env_Pathogen_data, Pars_wavelet_analysis) campylobacter_data_national<-wavelet_analysis(Env_Pathogen_data, Pars_wavelet_analysis)
......
...@@ -2,7 +2,7 @@ ...@@ -2,7 +2,7 @@
# Analysis was done following an uniform division of the range of the environmental variables, independently of the number of observations # 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 # 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)), { 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 ...@@ -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 #################### Create a data frame to fill with
Conditional_prevalence_unif <-data.frame(character(), character(), character(), numeric(), numeric()) Conditional_incidence_unif <-data.frame(character(), character(), character(), numeric(), numeric())
colnames(Conditional_prevalence_unif) <-c(variable_x, variable_y, variable_z,"counts","residents_tot") 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)){ ...@@ -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") # 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)){ ...@@ -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") 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_prevalence_unif,breaks_x,breaks_y,breaks_z)) 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))
})
}
...@@ -28,10 +28,8 @@ Plot_Conditional_incidence_Uniform<-function(Data_frame_1, Pars) ...@@ -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) 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) norm<-1/(1e6)
number_year<-number_days/365
norm<-(number_year)/(1e6)
...@@ -105,11 +103,200 @@ Plot_Conditional_incidence_Uniform<-function(Data_frame_1, Pars) ...@@ -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()
})
}
...@@ -174,7 +174,7 @@ write.table(diagnostic_laboratory_df[-c(length(diagnostic_laboratory_df))],paste ...@@ -174,7 +174,7 @@ write.table(diagnostic_laboratory_df[-c(length(diagnostic_laboratory_df))],paste
################ Define the time lag here ############################### ################ Define the time lag here ###############################
time_lag<-7 time_lag<-90
time_lag_char<-paste(time_lag) time_lag_char<-paste(time_lag)
## This create a database for all postocode where there was at least one case ## This create a database for all postocode where there was at least one case
......
0% Loading or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment