diff --git a/paper_plot_presentation_variables.R b/paper_plot_presentation_variables.R deleted file mode 100644 index 18f6b2dabc95872cfbf0473a4a13e38f51fd830b..0000000000000000000000000000000000000000 --- a/paper_plot_presentation_variables.R +++ /dev/null @@ -1,402 +0,0 @@ -# The code does look at how the risk of Campylobacter in humans depends on environmental variables - - -rm(list=ls(all=TRUE)) -# -library(ISOweek) -library(lubridate) -library(ggplot2) -require(MASS) -library(scales) -require(pheno) -library(timeDate) -library(pastecs) -library(stringi) -library(timeSeries) -library(wesanderson) -library(zoo) -library(plyr) - - -width<-60 -width_char<-paste(width) -situation<-"_const_hum_70" -#situation<-"" - -variable<-"humidity" -variable_df_1<-read.csv(paste("../../Data_Base/OPIE_data_base/",variable,".csv",sep="")) -humidity<-variable_df_1[,-c(1,2)] - -dates_s<-as.Date(as.character(variable_df_1[-c(1,2),2]),format="%Y-%m-%d") -dates<-rep(dates_s,times=length(variable_df_1)-2) -All_PC_s<-names(variable_df_1[1,]) -All_PC_s<-All_PC_s[-c(1,2)] -All_PC<-rep(All_PC_s,each=length(dates_s)) - -time_series_1<-variable_df_1[-c(1,2),] -#names(time_series_1)<-NULL - -time_series_2<-as.data.frame(time_series_1)[,-1] -time_series<-c() -time_series$Date<-rep(as.Date(time_series_2[,1]),times=length(time_series_2)-1) -time_series$humidity<-c(unlist(time_series_2[,-1])) - -time_series$yday<-yday(time_series$Date) -time_series$week<-week(time_series$Date) -time_series$month<-month(time_series$Date) -#time_series$Lab<-as.factor(time_series$Lab) - -#time_series_mean<-apply(time_series,1,mean) -#time_series_mean<-apply(time_series,1,mean) - -hum_df<-data.frame(time_series) -colnames(hum_df)<-c("Date","humidity") -hum_df$month<-month(hum_df$Date) -hum_df$week<-month(hum_df$Date) - - - - -############## Average per day of the year ################# -rownames(hum_df)<-NULL - -hum_df$yday<-as.factor(yday(hum_df$Date)) - -#time_series$year<-as.factor(year(time_series$Date)) -#time_series_total<-ddply(time_series,~year,function (x) x$Cases/sum(x$Cases)) - -#time_series$Cases4<-time_series$Cases4/sum(time_series$Cases4) -time_series_average1<-ddply(hum_df,~yday,summarise,mean=mean(humidity)) -time_series_quantile1<-ddply(hum_df,~yday, function (x) quantile(x$humidity, c(.25,.5,.75))) -time_series_average2<-cbind(time_series_average1,time_series_quantile1[,-1]) -time_series_average<- - data.frame(time_series_average2[,1],time_series_average2[,c(2:5)]) - - -df2<-data.frame(time_series_average,rep("Model",times=length(time_series_average[,1]))) -colnames(df2)<-c("Day","Mean","f_quant","median","s_quant","source") - -average_data<-df2 -average_data$Day<-as.numeric(average_data$Day) - -yearly_average<-ggplot(average_data,aes(x=Day,y=Mean)) - -yearly_average<-yearly_average+geom_line(size=2,colour="blue") -yearly_average <-yearly_average+ geom_ribbon(aes(ymin=f_quant, ymax=s_quant),fill="blue",alpha=0.25) -yearly_average<-yearly_average+ xlab("Day of the Year")+ ylab("Humidity %") - - -#yearly_average<-yearly_average+ scale_y_continuous(limit=c(0,0.05)) - - -month_2014<-month.abb[unique(month(as.Date(average_data$Day, origin = "2014-01-01")))] - -tick_pos<-yday(as.Date(as.character( - c("2014-01-01","2014-02-01","2014-03-01","2014-04-01","2014-05-01","2014-06-01", - "2014-07-01","2014-08-01","2014-09-01","2014-10-01","2014-11-01","2014-12-01") -))) - -yearly_average <-yearly_average+scale_x_continuous(breaks = tick_pos, labels = month_2014) -yearly_average<-yearly_average+ - 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,fill=NA)) - -#scale_fill_discrete(name="Humidity ") - -yearly_average<-yearly_average+ theme(axis.title.x =element_text( colour="#990000", size=13)) -yearly_average<-yearly_average+ theme(axis.title.y =element_text( colour="#990000", size=13)) - - - -yearly_average<-yearly_average+ theme(axis.title.x =element_text(size=13)) -yearly_average<-yearly_average+ theme(axis.title.y =element_text(size=13)) -yearly_average - - - - - - -#file_name_pdf<-paste("../Graphs/yearly_average_",variable_x,"_",variable,"n_seas",n_seas,"_multiple_delays.pdf", sep = "") -file_name_pdf<-paste("../../Graphs/yearly_average_humidity.pdf", sep = "") - -pdf(file = file_name_pdf,width = 6.94, height = 6.5) -yearly_average - -dev.off() - - -#file_name<-paste("../Graphs/yearly_average_",variable_x,"_",variable,"n_seas",n_seas,"_multiple_delays.tiff", sep = "") -file_name<-paste("../../Graphs/yearly_average_humidity.tiff", sep = "") - -tiff(filename = file_name,width = 17.35, height = 17.35, units = "cm", pointsize = 9, res = 600,compression = "lzw",antialias="default") -yearly_average - -dev.off() - - - -############## temperature ################# - - - - - -variable<-"max_air_temp" -variable_df_2<-read.csv(paste("../../Data_Base/OPIE_data_base/",variable,".csv",sep="")) -max_temp<-variable_df_2[,-c(1,2)] -max_temp<-max_temp[-c(1,2),] -names(max_temp) <- NULL -max_temp<-unlist(c(max_temp)) - - - -dates_s<-as.Date(as.character(variable_df_2[-c(1,2),2]),format="%Y-%m-%d") -dates<-rep(dates_s,times=length(variable_df_2)-2) -All_PC_s<-names(variable_df_2[1,]) -All_PC_s<-All_PC_s[-c(1,2)] -All_PC<-rep(All_PC_s,each=length(dates_s)) - -time_series_1<-variable_df_2[-c(1,2),] -#names(time_series_1)<-NULL - -time_series_2<-as.data.frame(time_series_1)[,-1] -time_series<-c() -time_series$Date<-rep(as.Date(time_series_2[,1]),times=length(time_series_2)-1) -time_series$max_air_temp<-c(unlist(time_series_2[,-1])) - -time_series$yday<-yday(time_series$Date) -time_series$week<-week(time_series$Date) -time_series$month<-month(time_series$Date) -#time_series$Lab<-as.factor(time_series$Lab) - -#time_series_mean<-apply(time_series,1,mean) -#time_series_mean<-apply(time_series,1,mean) - -max_temp_df<-data.frame(time_series) -colnames(max_temp_df)<-c("Date","max_air_temp") -max_temp_df$month<-month(max_temp_df$Date) -max_temp_df$week<-month(max_temp_df$Date) - - - - -############## Average per day of the year ################# -rownames(max_temp_df)<-NULL - -max_temp_df$yday<-as.factor(yday(max_temp_df$Date)) - -#time_series$year<-as.factor(year(time_series$Date)) -#time_series_total<-ddply(time_series,~year,function (x) x$Cases/sum(x$Cases)) - -#time_series$Cases4<-time_series$Cases4/sum(time_series$Cases4) -time_series_average1<-ddply(max_temp_df,~yday,summarise,mean=mean(max_air_temp)) -time_series_quantile1<-ddply(max_temp_df,~yday, function (x) quantile(x$max_air_temp, c(.25,.5,.75))) -time_series_average2<-cbind(time_series_average1,time_series_quantile1[,-1]) -time_series_average<- - data.frame(time_series_average2[,1],time_series_average2[,c(2:5)]) - - -df2<-data.frame(time_series_average,rep("Model",times=length(time_series_average[,1]))) -colnames(df2)<-c("Day","Mean","f_quant","median","s_quant","source") - -average_data<-df2 -average_data$Day<-as.numeric(average_data$Day) - -yearly_average<-ggplot(average_data,aes(x=Day,y=Mean)) - -yearly_average<-yearly_average+geom_line(size=2,colour="red") -yearly_average <-yearly_average+ geom_ribbon(aes(ymin=f_quant, ymax=s_quant),fill="red",alpha=0.25) -yearly_average<-yearly_average+ xlab("Day of the Year")+ ylab( - expression("Maximum Air Temperature "*~degree*C)) - - -#yearly_average<-yearly_average+ scale_y_continuous(limit=c(0,0.05)) - - -month_2014<-month.abb[unique(month(as.Date(average_data$Day, origin = "2014-01-01")))] - -tick_pos<-yday(as.Date(as.character( - c("2014-01-01","2014-02-01","2014-03-01","2014-04-01","2014-05-01","2014-06-01", - "2014-07-01","2014-08-01","2014-09-01","2014-10-01","2014-11-01","2014-12-01") -))) - -yearly_average <-yearly_average+scale_x_continuous(breaks = tick_pos, labels = month_2014) -yearly_average<-yearly_average+ - 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,fill=NA)) - -#scale_fill_discrete(name="Humidity ") - -yearly_average<-yearly_average+ theme(axis.title.x =element_text( colour="#990000", size=13)) -yearly_average<-yearly_average+ theme(axis.title.y =element_text( colour="#990000", size=13)) - - - -yearly_average<-yearly_average+ theme(axis.title.x =element_text(size=13)) -yearly_average<-yearly_average+ theme(axis.title.y =element_text(size=13)) -yearly_average - - - - - - -#file_name_pdf<-paste("../Graphs/yearly_average_",variable_x,"_",variable,"n_seas",n_seas,"_multiple_delays.pdf", sep = "") -file_name_pdf<-paste("../../Graphs/yearly_average_temperature.pdf", sep = "") - -pdf(file = file_name_pdf,width = 6.94, height = 6.5) -yearly_average - -dev.off() - - -#file_name<-paste("../Graphs/yearly_average_",variable_x,"_",variable,"n_seas",n_seas,"_multiple_delays.tiff", sep = "") -file_name<-paste("../../Graphs/yearly_average_temperature.tiff", sep = "") - -tiff(filename = file_name,width = 17.35, height = 17.35, units = "cm", pointsize = 9, res = 600,compression = "lzw",antialias="default") -yearly_average - -dev.off() - - - - - - - -Env_laboratory<-read.csv(paste("../../Data_Base/Cases_Environment/Laboratory_",width_char,".csv",sep="")) -Env_laboratory<-Env_laboratory[,-1] -colnames(Env_laboratory)<-c("PostCode","Date","humidity","max_temp","min_temp","rain","cum_rain","wind","residents") -Env_laboratory<-subset(Env_laboratory,year(Env_laboratory$Date)>=1999 & year(Env_laboratory$Date)<2000) - -#Env_laboratory$PostCode<-as.character(Env_laboratory$PostCode) -#Env_Campylobacter_data$PostCode<-as.character(Env_Campylobacter_data$PostCode) - -#wt<-c(0) -#for (i in c(1:length(levels(Env_Campylobacter_data$PostCode)) )){ -# wt<-c(wt,which(Env_laboratory$PostCode==levels(Env_Campylobacter_data$PostCode)[i])) -# print(c(100*i/length(levels(Env_Campylobacter_data$PostCode)),levels(Env_Campylobacter_data$PostCode)[i]) ) -#} - - -Env_laboratory_PHE<-Env_laboratory -##Env_laboratory_PHE<-merge(Env_laboratory,Env_Campylobacter_data,by='PostCode',all=T) not sure why this is not working, so loop above - -######################## include daylength ################## - -Coord_laboratory<-read.csv(paste("../../Data_Base/Cases/Lab_PostCodes.csv",sep="")) - -lat_long_lab<-data.frame(names(Coord_laboratory),as.numeric(Coord_laboratory[1,]),as.numeric(Coord_laboratory[2,])) -colnames(lat_long_lab)<-c("PostCode","lat","long") -Env_laboratory_int2<-merge(Env_laboratory,lat_long_lab,by="PostCode") -Env_Campylobacter_data_int2<-merge(Env_Campylobacter_data,lat_long_lab,by="PostCode") - -daylength<-function(lat,day_year) -{ - #Latitude measure in degrees - P <- asin(.39795*cos(.2163108 + 2*atan(.9671396*tan(.00860*(day_year-186))))) - Denom<-cos(lat*pi/180)*cos(P) - Numer<-sin(0.8333*pi/180) + sin(lat*pi/180)*sin(P) - D<-24-(24/pi)*acos(Numer/Denom) - return(D) -} - -latitude<-Env_laboratory_int2$lat -day_of_the_year<-yday(as.Date(Env_laboratory_int2$Date)) - -daylength_int1<-mapply(daylength, latitude, day_of_the_year) -daylength_df<-data.frame(latitude, day_of_the_year,as.Date(Env_laboratory_int2$Date),daylength_int1) -colnames(daylength_df)<-c("lat","day_year","Date","daylength") - - -daylength_df$month<-month(daylength_df$Date) -daylength_df$week<-month(daylength_df$Date) - - - - -############## Average per day of the year ################# -#rownames(daylength_df)<-NULL - -daylength_df$yday<-as.factor(yday(daylength_df$Date)) - -#time_series$year<-as.factor(year(time_series$Date)) -#time_series_total<-ddply(time_series,~year,function (x) x$Cases/sum(x$Cases)) - -#time_series$Cases4<-time_series$Cases4/sum(time_series$Cases4) -time_series_average1<-ddply(daylength_df,~yday,summarise,mean=mean(daylength)) -time_series_quantile1<-ddply(daylength_df,~yday, function (x) quantile(x$mdaylength, c(.25,.5,.75))) -time_series_average2<-cbind(time_series_average1,time_series_quantile1[,-1]) -time_series_average<- - data.frame(time_series_average2[,1],time_series_average2[,c(2:5)]) - - -df2<-data.frame(time_series_average,rep("Model",times=length(time_series_average[,1]))) -colnames(df2)<-c("Day","Mean","f_quant","median","s_quant","source") - -average_data<-df2 -average_data$Day<-as.numeric(average_data$Day) - -yearly_average<-ggplot(average_data,aes(x=Day,y=Mean)) - -yearly_average<-yearly_average+geom_line(size=2,colour="yellow") -yearly_average <-yearly_average+ geom_ribbon(aes(ymin=f_quant, ymax=s_quant),fill="yellow",alpha=0.25) -yearly_average<-yearly_average+ xlab("Day of the Year")+ ylab("Day Light (hours)") - - -#yearly_average<-yearly_average+ scale_y_continuous(limit=c(0,0.05)) - - -month_2014<-month.abb[unique(month(as.Date(average_data$Day, origin = "2014-01-01")))] - -tick_pos<-yday(as.Date(as.character( - c("2014-01-01","2014-02-01","2014-03-01","2014-04-01","2014-05-01","2014-06-01", - "2014-07-01","2014-08-01","2014-09-01","2014-10-01","2014-11-01","2014-12-01") -))) - -yearly_average <-yearly_average+scale_x_continuous(breaks = tick_pos, labels = month_2014) -yearly_average<-yearly_average+ - 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,fill=NA)) - -#scale_fill_discrete(name="Humidity ") - -yearly_average<-yearly_average+ theme(axis.title.x =element_text( colour="#990000", size=13)) -yearly_average<-yearly_average+ theme(axis.title.y =element_text( colour="#990000", size=13)) - - - -yearly_average<-yearly_average+ theme(axis.title.x =element_text(size=13)) -yearly_average<-yearly_average+ theme(axis.title.y =element_text(size=13)) -yearly_average - - - - - - -#file_name_pdf<-paste("../Graphs/yearly_average_",variable_x,"_",variable,"n_seas",n_seas,"_multiple_delays.pdf", sep = "") -file_name_pdf<-paste("../../Graphs/lat_average_daylength.pdf", sep = "") - -pdf(file = file_name_pdf,width = 6.94, height = 6.5) -yearly_average - -dev.off() - - -#file_name<-paste("../Graphs/yearly_average_",variable_x,"_",variable,"n_seas",n_seas,"_multiple_delays.tiff", sep = "") -file_name<-paste("../../Graphs/lat_average_daylength.tiff", sep = "") - -tiff(filename = file_name,width = 17.35, height = 17.35, units = "cm", pointsize = 9, res = 600,compression = "lzw",antialias="default") -yearly_average - -dev.off() - - - - -