diff --git a/Taylor_Simulated_Campylobacter_environment_light_hum_min_for_rec_delay_original_modified.R b/Taylor_Simulated_Campylobacter_environment_light_hum_min_for_rec_delay_original_modified.R deleted file mode 100644 index 4123b107946df56320c15cfc2c65c32b7e2e6371..0000000000000000000000000000000000000000 --- a/Taylor_Simulated_Campylobacter_environment_light_hum_min_for_rec_delay_original_modified.R +++ /dev/null @@ -1,718 +0,0 @@ -# The code does look at how the risk of Campylobacter in humans depends on environmental variables -# this to calculate delay - -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(plyr) - -n_seas<-1 -width<-14 -width_char<-paste(width) -## Varaible file - - -variable_x<-"Minimum_air_temperature" -variable_y<-"Relative_humidity" -variable_z<-"daylength" - -variable_x2<-"min_air_temp" -variable_y2<-"humidity" -variable_z2<-"light" - - -## Varaible file - - - -Env_Campylobacter_data_all2<-read.csv(paste("../../Data_Base/Cases_Environment/Simulated_Campylobacter_environment_",width_char,"_original_MEDMI.csv",sep="")) -Env_Campylobacter_data_all2<-Env_Campylobacter_data_all2[,-1] -colnames(Env_Campylobacter_data_all2)<-c("PostCode","Date","Cases", - "Maximum_air_temperature", - "Minimum_air_temperature", - "Mean_wind_speed", - "Cumul_Precipitation", - "Mean_Precipitation", - "Relative_humidity", - "daylength", - "residents") - -Env_Campylobacter_data<-subset(Env_Campylobacter_data_all2,year(Env_Campylobacter_data_all2$Date)>=1990 & year(Env_Campylobacter_data_all2$Date)<=2015) - -dates_s<- dates_s<- seq(as.Date("1990-01-01"), as.Date("2015-12-31"), by = "day") -All_PC_s<-levels(as.factor(Env_Campylobacter_data$PostCode)) -All_PC<-rep(All_PC_s,each=length(dates_s)) - - - -Env_laboratory<-read.csv(paste("../../Data_Base/Cases_Environment/Simulated_Laboratory_",width_char,"_original_MEDMI.csv",sep="")) -Env_laboratory<-Env_laboratory[,-1] -colnames(Env_laboratory)<-c("PostCode","Date", - "Maximum_air_temperature", - "Minimum_air_temperature", - "Mean_wind_speed", - "Cumul_Precipitation", - "Mean_Precipitation", - "Relative_humidity", - "daylength", - "residents") - -Env_laboratory<-subset(Env_laboratory,year(Env_laboratory$Date)>=1990 & year(Env_laboratory$Date)<=2015) -Env_laboratory$Date<-as.Date(Env_laboratory$Date) - - - -#Env_laboratory_PHE<-Env_laboratory[wt[-1],] -##Env_laboratory_PHE<-merge(Env_laboratory,Env_Campylobacter_data,by='PostCode',all=T) not sure why this is not working, so loop above -Env_laboratory_PHE<-Env_laboratory -#All_residents_lab<-ddply(Env_laboratory_PHE,~PostCode,summarise,tot=mean(residents)) -#All_residents<-sum(All_residents_lab$tot) - -######################## 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") -#Env_laboratory_int2<-Env_laboratory_int2[,-c(16,17)] -#names(Env_laboratory_int2)[10]<-"lat" -#names(Env_laboratory_int2)[11]<-"long" -#names(Env_laboratory_int2)[14]<-"Date2" - -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") - - -Env_laboratory<-data.frame(Env_laboratory_int2,daylength_df) - -### repeat for the data only #### - -latitude<-Env_Campylobacter_data_int2$lat -day_of_the_year<-yday(as.Date(Env_Campylobacter_data_int2$Date)) - -daylength_int1<-mapply(daylength, latitude, day_of_the_year) -daylength_df<-data.frame(latitude, day_of_the_year,as.Date(Env_Campylobacter_data_int2$Date),daylength_int1) -colnames(daylength_df)<-c("lat","day_year","Date","daylength") - -#daylength_df$lat<-as.factor(daylength_df$lat) -#daylength_df$Date<-as.factor(daylength_df$Date) -#Env_Lyme_data_int2$lat<-as.factor(Env_Lyme_data_int2$lat) -#Env_Lyme_data_int2$Date<-as.factor(Env_Lyme_data_int2$Date) - -Env_Campylobacter_data<-data.frame(Env_Campylobacter_data_int2,daylength_df) - - - - -######################## - - - - - - - - - -var_x_loc_df_all<-read.csv(paste("../../Data_Base/Cases_Environment/",variable_z2,"_",variable_y2,"_",variable_x2,"_",width_char,"_Simulated_for_rec_original_MEDMI.csv",sep="")) -var_x_loc_df_all<-var_x_loc_df_all[,-1] -colnames(var_x_loc_df_all)<-c(variable_x,"breaks","prop","incidence","month",variable_z,variable_y,"counts","counts_tot","residents","residents_tot","Numb_Lab") - -var_x_loc_df_all2<-na.omit(var_x_loc_df_all) - - -################### - - -delta_light<-1 -delta_hum<-5 -delta_temp<-1 -delta_rain<-2 -delta_cum_rain<-2 -delta_wind<-1 - -breaks_hum<-seq(max(min(na.omit(Env_laboratory$Relative_humidity))-10,0),max(na.omit(Env_laboratory$Relative_humidity))+10,by=delta_hum) #i -breaks_min_temp<-seq(min(na.omit(Env_laboratory$Minimum_air_temperature))-2, max(na.omit(Env_laboratory$Minimum_air_temperature))+2,by=delta_temp) -breaks_max_temp<-seq(min(na.omit(Env_laboratory$Maximum_air_temperature))-2, max(na.omit(Env_laboratory$Maximum_air_temperature))+2,by=delta_temp) -breaks_rain<-seq(max(min(na.omit(Env_laboratory$Mean_Precipitation)),0), max(na.omit(Env_laboratory$Mean_Precipitation))+2,by=delta_rain) -breaks_wind<-seq(max(min(na.omit(Env_laboratory$Mean_wind_speed))-2,0),max(na.omit(Env_laboratory$Mean_wind_speed))+2,by=delta_wind) -breaks_cum_rain<-seq(max(min(na.omit(Env_laboratory$Cumul_Precipitation)),0), max(na.omit(Env_laboratory$Cumul_Precipitation))+2,by=delta_cum_rain) -breaks_mean_temp<-seq(min(na.omit(Env_laboratory$Minimum_air_temperature))-2,max(na.omit(Env_laboratory$Maximum_air_temperature))+2,by=delta_temp) -breaks_light<-seq(max(min(na.omit(Env_laboratory$daylength))-1,0),max(na.omit(Env_laboratory$daylength))+1,by=delta_light) - - -i_hum_min<-max(which(breaks_hum<=min(na.omit(Env_Campylobacter_data$Relative_humidity)))) -i_hum_max<-max(which(breaks_hum<=max(na.omit(Env_Campylobacter_data$Relative_humidity)))) - -i_min_temp_min<-max(which(breaks_min_temp<=min(na.omit(Env_Campylobacter_data$Minimum_air_temperature)))) -i_min_temp_max<-max(which(breaks_min_temp<=max(na.omit(Env_Campylobacter_data$Minimum_air_temperature)))) - -i_max_temp_min<-max(which(breaks_max_temp<=min(na.omit(Env_Campylobacter_data$Maximum_air_temperature)))) -i_max_temp_max<-max(which(breaks_max_temp<=max(na.omit(Env_Campylobacter_data$Maximum_air_temperature)))) - -i_rain_min<-max(which(breaks_rain<=min(na.omit(Env_Campylobacter_data$Mean_Precipitation)))) -i_rain_max<-max(which(breaks_rain<=max(na.omit(Env_Campylobacter_data$Mean_Precipitation)))) - -i_cum_rain_min<-max(which(breaks_cum_rain<=min(na.omit(Env_Campylobacter_data$Cumul_Precipitation)))) -i_cum_rain_max<-max(which(breaks_cum_rain<=max(na.omit(Env_Campylobacter_data$Cumul_Precipitation)))) - -i_wind_min<-max(which(breaks_wind<=min(na.omit(Env_Campylobacter_data$Mean_wind_speed)))) -i_wind_max<-max(which(breaks_wind<=max(na.omit(Env_Campylobacter_data$Mean_wind_speed)))) - - -i_light_min<-max(which(breaks_light<=min(na.omit(Env_Campylobacter_data$daylength)))) -i_light_max<-max(which(breaks_light<=max(na.omit(Env_Campylobacter_data$daylength)))) - -############## from here '################### - - -if (variable_z=="Maximum_air_temperature"){ - - breaks_var<-breaks_max_temp - i_var_min<-i_max_temp_min - i_var_max<-i_max_temp_max - Env_Campylobacter_data_var<-Env_Campylobacter_data$Maximum_air_temperature - Env_laboratory_var<-Env_laboratory$Maximum_air_temperature -} -if (variable_x=="Maximum_air_temperature"){ - - i_var_x_min<-i_max_temp_min - i_var_x_max<-i_max_temp_max - breaks_var_x<-breaks_max_temp -} - -if (variable_y=="Maximum_air_temperature"){ - - i_var_y_min<-i_max_temp_min - i_var_y_max<-i_max_temp_max - breaks_var_y<-breaks_max_temp -} - -if (variable_z=="Minimum_air_temperature"){ - - breaks_var<-breaks_min_temp - i_var_min<-i_min_temp_min - i_var_max<-i_min_temp_max - Env_Campylobacter_data_var<-Env_Campylobacter_data$min_temp - Env_laboratory_var<-Env_laboratory$min_temp -} -if (variable_x=="Minimum_air_temperature"){ - - i_var_x_min<-i_min_temp_min - i_var_x_max<-i_min_temp_max - breaks_var_x<-breaks_min_temp -} - -if (variable_y=="Minimum_air_temperature"){ - - i_var_y_min<-i_min_temp_min - i_var_y_max<-i_min_temp_max - breaks_var_y<-breaks_min_temp -} - -if (variable_z=="Relative_humidity"){ - - breaks_var<-breaks_hum - i_var_min<-i_hum_min - i_var_max<-i_hum_max - Env_Campylobacter_data_var<-Env_Campylobacter_data$Relative_humidity - Env_laboratory_var<-Env_laboratory$Relative_humidity -} -if (variable_x=="Relative_humidity"){ - - i_var_x_min<-i_hum_min - i_var_x_max<-i_hum_max - breaks_var_x<-breaks_hum -} - -if (variable_y=="Relative_humidity"){ - - i_var_y_min<-i_hum_min - i_var_y_max<-i_hum_max - breaks_var_y<-breaks_hum -} - - -if (variable_z=="mean_temp"){ - - breaks_var<-breaks_mean_temp - i_var_min<-i_mean_temp_min - i_var_max<-i_mean_temp_max - Env_Campylobacter_data_var<-Env_Campylobacter_data$mean_temp - Env_laboratory_var<-Env_laboratory$mean_temp -} -if (variable_x=="mean_temp"){ - - i_var_x_min<-i_mean_temp_min - i_var_x_max<-i_mean_temp_max - breaks_var_x<-breaks_mean_temp -} - -if (variable_y=="mean_temp"){ - - i_var_y_min<-i_mean_temp_min - i_var_y_max<-i_mean_temp_max - breaks_var_y<-breaks_mean_temp -} - -if (variable_z=="Mean_Precipitation"){ - - breaks_var<-breaks_rain - i_var_min<-i_rain_min - i_var_max<-i_rain_max - Env_Campylobacter_data_var<-Env_Campylobacter_data$Mean_Precipitation - Env_laboratory_var<-Env_laboratory$Mean_Precipitation -} -if (variable_x=="Mean_Precipitation"){ - - i_var_x_min<-i_rain_min - i_var_x_max<-i_rain_max - breaks_var_x<-breaks_rain -} - -if (variable_y=="Mean_Precipitation"){ - - i_var_y_min<-i_rain_min - i_var_y_max<-i_rain_max - breaks_var_y<-breaks_rain -} - -if (variable_z=="Cumul_Precipitation"){ - - breaks_var<-breaks_cum_rain - i_var_min<-i_cum_rain_min - i_var_max<-i_cum_rain_max - Env_Campylobacter_data_var<-Env_Campylobacter_data$Cumul_Precipitation - Env_laboratory_var<-Env_laboratory$Cumul_Precipitation -} -if (variable_x=="Cumul_Precipitation"){ - - i_var_x_min<-i_cum_rain_min - i_var_x_max<-i_cum_rain_max - breaks_var_x<-breaks_cum_rain -} - -if (variable_y=="Cumul_Precipitation"){ - - i_var_y_min<-i_cum_rain_min - i_var_y_max<-i_cum_rain_max - breaks_var_y<-breaks_cum_rain -} - - -if (variable_z=="Mean_wind_speed"){ - - breaks_var<-breaks_wind - i_var_min<-i_wind_min - i_var_max<-i_wind_max - Env_Campylobacter_data_var<-Env_Campylobacter_data$Mean_wind_speed - Env_laboratory_var<-Env_laboratory$Mean_wind_speed -} -if (variable_x=="Mean_wind_speed"){ - - i_var_x_min<-i_wind_min - i_var_x_max<-i_wind_max - breaks_var_x<-breaks_wind -} - -if (variable_y=="Mean_wind_speed"){ - - i_var_y_min<-i_wind_min - i_var_y_max<-i_wind_max - breaks_var_y<-breaks_wind -} - - - - -if (variable_z=="daylength"){ - - breaks_var<-breaks_light - i_var_min<-i_light_min - i_var_max<-i_light_max - Env_Campylobacter_data_var<-Env_Campylobacter_data$daylength - Env_laboratory_var<-Env_laboratory$daylength -} -if (variable_x=="daylength"){ - - i_var_x_min<-i_light_min - i_var_x_max<-i_light_max - breaks_var_x<-breaks_light -} - -if (variable_y=="daylength"){ - - i_var_y_min<-i_light_min - i_var_y_max<-i_light_max - breaks_var_y<-breaks_light -} - - - - -D_eta_var_x<-function(i_var,i_var_y,i_var_x) -{ - - - - -#var_x_loc_df$eta<-var_x_loc_df$prop*var_x_loc_df$Numb_Lab -var_x_loc_df$eta<-var_x_loc_df$incidence*All_residents_lab[,2]*var_x_loc_df$Numb_Lab #All_residents_lab[,2] this change by postcode only - -if(is.na(breaks_var[i_var+1])==TRUE) -{ - Yt0<-subset(var_x_loc_df,var_x_loc_df$daylength>=breaks_var[i_var]) -} else { - Yt0<-subset(var_x_loc_df,var_x_loc_df$daylength>=breaks_var[i_var] & var_x_loc_df$daylength<breaks_var[i_var+1]) -} - - -if(is.na(breaks_var_y[i_var_y+1])==TRUE) -{ - Yt1<-subset(Yt0,Yt0$Relative_humidity>=breaks_var_y[i_var_y]) -} else { - Yt1<-subset(Yt0,Yt0$Relative_humidity>=breaks_var_y[i_var_y] & Yt0$Relative_humidity<breaks_var_y[i_var_y+1]) -} - - - - if(is.na(breaks_var_x[i_var_x+1])==TRUE) - { - Yt2<-subset(Yt1,Yt1$Minimum_air_temperature>=breaks_var_x[i_var_x]) - } else { - Yt2<-subset(Yt1,Yt1$Minimum_air_temperature>=breaks_var_x[i_var_x] & Yt1$Minimum_air_temperature<breaks_var_x[i_var_x+1]) - } - - - if(is.na(breaks_var_x[i_var_x+2])==TRUE) - { - Yt3<-subset(Yt1,Yt1$Minimum_air_temperature>=breaks_var_x[i_var_x+1]) - } else { - Yt3<-subset(Yt1,Yt1$Minimum_air_temperature>=breaks_var_x[i_var_x+1] & Yt1$Minimum_air_temperature<breaks_var_x[i_var_x+2]) - } - - if(length(Yt3[,1])!=0 & length(Yt2[,1])!=0){ - D_eta_D_x<- (mean(Yt3$eta- Yt2$eta))/(breaks_var_x[i_var_x+1]-breaks_var_x[i_var_x]) - } else if (length(Yt3[,1])!=0 & length(Yt2[,1])==0){ - D_eta_D_x<- 0.5*mean(Yt3$eta)/(breaks_var_x[i_var_x+1]-breaks_var_x[i_var_x]) - } else if (length(Yt3[,1])==0 & length(Yt2[,1])!=0){ - D_eta_D_x<- 0.5*mean(Yt2$eta)/(breaks_var_x[i_var_x+1]-breaks_var_x[i_var_x]) - } else if (length(Yt3[,1])==0 & length(Yt2[,1])==0){ - D_eta_D_x<-NA - } - - return(D_eta_D_x) - } - - - - -D_eta_var_y<-function(i_var,i_var_y,i_var_x) -{ - - - #var_x_loc_df$eta<-var_x_loc_df$prop*var_x_loc_df$Numb_Lab - var_x_loc_df$eta<-var_x_loc_df$incidence*All_residents_lab[,2]*var_x_loc_df$Numb_Lab - - - if(is.na(breaks_var[i_var+1])==TRUE) - { - Yt0<-subset(var_x_loc_df,var_x_loc_df$daylength>=breaks_var[i_var]) - } else { - Yt0<-subset(var_x_loc_df,var_x_loc_df$daylength>=breaks_var[i_var] & var_x_loc_df$daylength<breaks_var[i_var+1]) - } - - - if(is.na(breaks_var_x[i_var_x+1])==TRUE) - { - Yt1<-subset(Yt0,Yt0$Minimum_air_temperature>=breaks_var_x[i_var_x]) - } else { - Yt1<-subset(Yt0,Yt0$Minimum_air_temperature>=breaks_var_x[i_var_x] & Yt0$Minimum_air_temperature<breaks_var_x[i_var_x+1]) - } - - - - if(is.na(breaks_var_y[i_var_y+1])==TRUE) - { - Yt2<-subset(Yt1,Yt1$Relative_humidity>=breaks_var_y[i_var_y]) - } else { - Yt2<-subset(Yt1,Yt1$Relative_humidity>=breaks_var_y[i_var_y] & Yt1$Relative_humidity<breaks_var_y[i_var_y+1]) - } - - - if(is.na(breaks_var_y[i_var_y+2])==TRUE) - { - Yt3<-subset(Yt1,Yt1$Relative_humidity>=breaks_var_y[i_var_y+1]) - } else { - Yt3<-subset(Yt1,Yt1$Relative_humidity>=breaks_var_y[i_var_y+1] & Yt1$Relative_humidity<breaks_var_y[i_var_y+2]) - } - - if(length(Yt3[,1])!=0 & length(Yt2[,1])!=0){ - D_eta_D_y<- (mean(Yt3$eta- Yt2$eta))/(breaks_var_y[i_var_y+1]-breaks_var_y[i_var_y]) - } else if (length(Yt3[,1])!=0 & length(Yt2[,1])==0){ - D_eta_D_y<- 0.5*mean(Yt3$eta)/(breaks_var_y[i_var_y+1]-breaks_var_y[i_var_y]) - } else if (length(Yt3[,1])==0 & length(Yt2[,1])!=0){ - D_eta_D_y<- 0.5*mean(Yt2$eta)/(breaks_var_y[i_var_y+1]-breaks_var_y[i_var_y]) - } else if (length(Yt3[,1])==0 & length(Yt2[,1])==0){ - D_eta_D_y<-NA - } - - return(D_eta_D_y) -} - - - - -D_eta_var<-function(i_var,i_var_y,i_var_x) -{ - - -#var_x_loc_df$eta<-var_x_loc_df$prop*var_x_loc_df$Numb_Lab -var_x_loc_df$eta<-var_x_loc_df$incidence*All_residents_lab[,2]*var_x_loc_df$Numb_Lab - - if(is.na(breaks_var_x[i_var_x+1])==TRUE) - { - Yt0<-subset(var_x_loc_df,var_x_loc_df$Minimum_air_temperature>=breaks_var_x[i_var_x]) - } else { - Yt0<-subset(var_x_loc_df,var_x_loc_df$Minimum_air_temperature>=breaks_var_x[i_var_x] & var_x_loc_df$Minimum_air_temperature<breaks_var_x[i_var_x+1]) - } - -if(is.na(breaks_var_y[i_var_y+1])==TRUE) -{ - Yt1<-subset(Yt0,Yt0$Relative_humidity>=breaks_var_y[i_var_y]) -} else { - Yt1<-subset(Yt0,Yt0$Relative_humidity>=breaks_var_y[i_var_y] & Yt0$Relative_humidity<breaks_var_y[i_var_y+1]) -} - - - if(is.na(breaks_var[i_var+1])==TRUE) - { - Yt2<-subset( Yt1, Yt1$daylength>=breaks_var[i_var]) - } else { - Yt2<-subset( Yt1, Yt1$daylength>=breaks_var[i_var] & Yt1$daylength<breaks_var[i_var+1]) - } - -if(is.na(breaks_var[i_var+2])==TRUE) -{ - Yt3<-subset( Yt1, Yt1$daylength>=breaks_var[i_var+1]) -} else { - Yt3<-subset( Yt1, Yt1$Relative_humiditylight>=breaks_var[i_var+1] & Yt1$daylength<breaks_var[i_var+2]) -} - - - - if(length(Yt3[,1])!=0 & length(Yt2[,1])!=0){ - D_eta_D_v<- mean(Yt3$eta- Yt2$eta)/(breaks_var[i_var+1]-breaks_var[i_var]) - } else if (length(Yt3[,1])!=0 & length(Yt2[,1])==0){ - D_eta_D_v<- 0.5*mean(Yt3$eta)/(breaks_var[i_var+1]-breaks_var[i_var]) - } else if (length(Yt3[,1])==0 & length(Yt2[,1])!=0){ - D_eta_D_v<- 0.5*mean(Yt2$eta)/(breaks_var[i_var+1]-breaks_var[i_var]) - } else if (length(Yt3[,1])==0 & length(Yt2[,1])==0){ - D_eta_D_v<-NA - } - - return(D_eta_D_v) - } - - - - - -df_all_1<-c() -#i_var<-seq(i_var_min,i_var_max) -#i_var_x<-seq(i_var_x_min,i_var_x_max) - - for (j in c(1: length(All_PC_s))){ - -var_x_loc_df<-var_x_loc_df_all2 - variable_df<-subset(Env_laboratory,Env_laboratory$PostCode==All_PC_s[j]) - -All_residents_lab<-ddply(variable_df,~PostCode,summarise,tot=mean(residents)) -for (i_var in c(seq(i_var_min,i_var_max))){ -for (i_var_y in c(seq(i_var_y_min,i_var_y_max))){ - for (i_var_x in c(seq(i_var_x_min,i_var_x_max))){ - - -#df1<-data.frame(i_var,i_var_x,t(unlist(lapply(i_var,D_eta_var_x,i_var_x=i_var_x)))) -#df1<-data.frame(i_var,i_var_x,t(unlist(lapply(i_var,D_eta_var_x,i_var_x=i_var_x)))) - -df1<-data.frame(breaks_light[i_var],breaks_hum[i_var_y],breaks_min_temp[i_var_x],D_eta_var_x(i_var,i_var_y,i_var_x),D_eta_var_y(i_var,i_var_y,i_var_x),D_eta_var(i_var,i_var_y,i_var_x),All_PC_s[j] ) -colnames(df1)<-c("daylength","Relative_humidity","Minimum_air_temperature","D_eta_D_max_temp","D_eta_D_hum","D_eta_D_light","PostCode") -df_all_1<-rbind(df_all_1,df1) -#mapply(D_eta_var_x,i_var, MoreArgs =list(i_var_x)) -} -} - -} -} - -#df_all<-ddply(df_all_1,~PostCode,summarise,tot=mean()) -df_all<-df_all_1 - -time_series<-c() - -for (i in c(1: length(All_PC_s))){ - #for (i in c(1: 15)){ - - - - #for (j in c(1: n_seas)){ - - # n_months<-12/n_seas - - variable_df<-subset(Env_laboratory,Env_laboratory$PostCode==All_PC_s[i]) - if (length(variable_df[,1])!=0){ - - - -Delta_T<-diff(variable_df$Minimum_air_temperature) -Delta_RH<-diff(variable_df$Relative_humidity) -Delta_L<-diff(variable_df$daylength) -Delta_var<-data.frame(variable_df$Date[-1],Delta_T,Delta_RH,Delta_L) -colnames(Delta_var)<-c("dates","diff_max_temp","diff_humidity","diff_light") - -x<-variable_df$Relative_humidity -y<-variable_df$Minimum_air_temperature -z<-variable_df$daylength - -variable_df_1_dis<-data.frame(variable_df$Date,variable_df$Relative_humidity) -variable_df_2_dis<-data.frame(variable_df$Date,variable_df$Minimum_air_temperature) -variable_df_3_dis<-data.frame(variable_df$Date,variable_df$daylength) - -colnames(variable_df_1_dis)<-c("dates","Relative_humidity") -colnames(variable_df_2_dis)<-c("dates","Minimum_air_temperature") -colnames(variable_df_3_dis)<-c("dates","daylength") - - variable_df_1_dis$dates<-as.Date(as.character((variable_df_1_dis$dates))) - variable_df_2_dis$dates<-as.Date(as.character((variable_df_2_dis$dates))) - variable_df_3_dis$dates<-as.Date(as.character((variable_df_3_dis$dates))) - - variable_df_1_dis<-subset(variable_df_1_dis,year(variable_df_1_dis$dates)>=1990 & year(variable_df_1_dis$dates)<=2015) - variable_df_2_dis<-subset(variable_df_2_dis,year(variable_df_2_dis$dates)>=1990 & year(variable_df_2_dis$dates)<=2015) - variable_df_3_dis<-subset(variable_df_3_dis,year(variable_df_3_dis$dates)>=1990 & year(variable_df_3_dis$dates)<=2015) - - variable_df_1_dis$Relative_humidity<-(breaks_hum[findInterval(x, breaks_hum)]) - variable_df_2_dis$Minimum_air_temperature<-(breaks_min_temp[findInterval(y, breaks_min_temp)]) - variable_df_3_dis$daylength<-(breaks_light[findInterval(z, breaks_light)]) - - -# var_x_loc_df<-subset(var_x_loc_df_all2,var_x_loc_df_all2$month==j) -# var_x_loc_df$Relative_humidity<-floor(var_x_loc_df$Relative_humidity) -# var_x_loc_df$Minimum_air_temperature<-floor(var_x_loc_df$breaks) - #df<-subset(df_all,df_all$month==j) - df<-subset(df_all,df_all$PostCode==All_PC_s[i]) - - variable_df_4_dis<-merge(variable_df_1_dis,df, by="Relative_humidity") - variable_df_5_dis<-merge(variable_df_2_dis,df, by="Minimum_air_temperature") - variable_df_6_dis<-merge(variable_df_3_dis,df, by="daylength") - - - - variable_df_4_dis$dates<-as.Date(as.character((variable_df_4_dis$dates))) - variable_df_4_dis<-variable_df_4_dis[order(variable_df_4_dis$dates),] - variable_df_5_dis$dates<-as.Date(as.character((variable_df_5_dis$dates))) - variable_df_5_dis<-variable_df_5_dis[order(variable_df_5_dis$dates),] - variable_df_6_dis$dates<-as.Date(as.character((variable_df_6_dis$dates))) - variable_df_6_dis<-variable_df_6_dis[order(variable_df_6_dis$dates),] - - - variable_df_4_dis$dates<-as.factor(variable_df_4_dis$dates) - variable_df_5_dis$dates<-as.factor(variable_df_5_dis$dates) - variable_df_6_dis$dates<-as.factor(variable_df_6_dis$dates) - - variable_df_dis0<-merge(variable_df_4_dis,variable_df_5_dis, by=c("dates","Minimum_air_temperature","Relative_humidity","daylength") ) - variable_df_dis<-merge(variable_df_dis0,variable_df_6_dis, by=c("dates","Minimum_air_temperature","Relative_humidity","daylength")) - - variable_df_dis$dates<-as.Date(as.character((variable_df_dis$dates))) - variable_df_dis<-variable_df_dis[order(variable_df_dis$dates),] - variable_df_dis$dates<-as.factor(variable_df_dis$dates) - Delta_var$dates<-as.factor(Delta_var$dates) - variable_df_Taylor<-merge(variable_df_dis,Delta_var, by=c("dates")) - - variable_df_Taylor<-variable_df_Taylor[,c(1:7,17:19)] - colnames(variable_df_Taylor)<-c("dates",variable_x, variable_y, variable_z,"D_eta_D_max_temp","D_eta_D_hum","D_eta_D_light", - "diff_max_temp","diff_humidity","diff_light") - - - lambda_temp<-variable_df_Taylor$D_eta_D_max_temp*variable_df_Taylor$diff_max_temp - lambda_hum<-variable_df_Taylor$D_eta_D_hum*variable_df_Taylor$diff_humidity - lambda_light<-variable_df_Taylor$D_eta_D_light*variable_df_Taylor$diff_light - - - - time_series_1<- - data.frame(variable_df_Taylor$dates, - lambda_temp,lambda_hum,lambda_light, - rep(All_PC_s[i],times=length(variable_df_Taylor$dates)), - variable_df_Taylor$D_eta_D_max_temp,variable_df_Taylor$diff_max_temp, - variable_df_Taylor$D_eta_D_hum,variable_df_Taylor$diff_humidity, - variable_df_Taylor$D_eta_D_light,variable_df_Taylor$diff_light) - colnames(time_series_1)<-c("Date","Contr_temp","Contr_hum","Contr_light", - "Lab","delta_eta_temp","delta_temp","delta_eta_hum","delta_hum","delta_eta_light","delta_light") - time_series<-rbind(time_series,time_series_1) - - #print(100*c(i/length(All_PC_s),comp_cases,comp_cases2,comp_cases3,lambda )) - print(100*c(i/length(All_PC_s) )) - #print(" ") - #print(lambda) - } -} - -#time_series1<-subset(time_series,time_series$Lab=="B152TG") -#time_series1$Date<-as.Date(as.character(time_series1$Date)) -#time_series1$year<-as.factor(year(as.Date(as.character(time_series1$Date)))) -#time_series1<-subset(time_series1,time_series1$year=="2005") - -#time_series1<-na.omit(time_series1) -##temp_contr<-cumsum(time_series1$Contr_temp) -#temp_hum<-cumsum(time_series1$Contr_hum) -# -#df2<-data.frame(as.Date(as.character(time_series1$Date)),temp_contr,temp_hum) -#colnames(df2)<-c("Date","contr","humidity") -#temp_contr<-ddply(time_series1,~year,summarise,Cum_Cases=cumsum(Contr_temp)) -#time_series<-na.omit(time_series) -#time_series$temp_contr<-cumsum(time_series$Contr_temp) -#time_series$temp_hum<-cumsum(time_series$Contr_hum) -#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_lab<-ddply(time_series,~Date,summarise,tot=mean(Contr_temp)) -#time_series_lab2<-ddply(time_series,~yday,summarise,tot=mean(Contr_temp)) - -#time_series_lab<-ddply(time_series,~week,summarise,tot=mean(Contr_temp)) -#time_series_lab2<-ddply(time_series,~week,summarise,tot=mean(Contr_hum)) - -#time_series_lab<-ddply(time_series,~month,summarise,tot=mean(Contr_temp)) -#time_series_lab2<-ddply(time_series,~month,summarise,tot=mean(Contr_hum)) -#time_series_lab2<-ddply(time_series,~Date,summarise,tot=mean(Contr_hum)) - -write.table(time_series,paste("../../Data_Base/Cases/Taylor_contribution_Time_series_",variable_z,"_",variable_y,"_",variable_x,width_char,"_Simulated_original_MEDMI.csv",sep=""), col.names = FALSE, sep = ",",eol = "\n")