Delete...

Delete Taylor_Simulated_Campylobacter_environment_light_hum_min_for_rec_delay_original_modified.Rout

Taylor_Simulated_Campylobacter_environment_light_hum_min_for_rec_delay_original_modified.Rout

-
-R version 3.5.3 (2019-03-11) -- "Great Truth"
-Copyright (C) 2019 The R Foundation for Statistical Computing
-Platform: x86_64-pc-linux-gnu (64-bit)
-
-R is free software and comes with ABSOLUTELY NO WARRANTY.
-You are welcome to redistribute it under certain conditions.
-Type 'license()' or 'licence()' for distribution details.
-
-  Natural language support but running in an English locale
-
-R is a collaborative project with many contributors.
-Type 'contributors()' for more information and
-'citation()' on how to cite R or R packages in publications.
-
-Type 'demo()' for some demos, 'help()' for on-line help, or
-'help.start()' for an HTML browser interface to help.
-Type 'q()' to quit R.
-
-[Previously saved workspace restored]
-
-> # 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)
-
-Attaching package: ‘lubridate’
-
-The following object is masked from ‘package:base’:
-
-    date
-
-> library(ggplot2)
-> require(MASS)
-Loading required package: MASS
-> library(scales)
-> require(pheno)
-Loading required package: pheno
-Loading required package: nlme
-Loading required package: SparseM
-
-Attaching package: ‘SparseM’
-
-The following object is masked from ‘package:base’:
-
-    backsolve
-
-Loading required package: quantreg
-> library(timeDate)
-> library(pastecs)
-> library(stringi)
-> library(timeSeries)
-> library(wesanderson)
-> library(plyr)
-
-Attaching package: ‘plyr’
-
-The following object is masked from ‘package:lubridate’:
-
-    here
-
-> 
-> 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)
-+   }
-+ }
-[1] 0.4854369
-[1] 0.9708738
-[1] 1.456311
-[1] 1.941748
-[1] 2.427184
-[1] 2.912621
-[1] 3.398058
-[1] 3.883495
-[1] 4.368932
-[1] 4.854369
-[1] 5.339806
-[1] 5.825243
-[1] 6.31068
-[1] 6.796117
-[1] 7.281553
-[1] 7.76699
-[1] 8.252427
-[1] 8.737864
-[1] 9.223301
-[1] 9.708738
-[1] 10.19417
-[1] 10.67961
-[1] 11.16505
-[1] 11.65049
-[1] 12.13592
-[1] 12.62136
-[1] 13.1068
-[1] 13.59223
-[1] 14.07767
-[1] 14.56311
-[1] 15.04854
-[1] 15.53398
-[1] 16.01942
-[1] 16.50485
-[1] 16.99029
-[1] 17.47573
-[1] 17.96117
-[1] 18.4466
-[1] 18.93204
-[1] 19.41748
-[1] 19.90291
-[1] 20.38835
-[1] 20.87379
-[1] 21.35922
-[1] 21.84466
-[1] 22.3301
-[1] 22.81553
-[1] 23.30097
-[1] 23.78641
-[1] 24.27184
-[1] 24.75728
-[1] 25.24272
-[1] 25.72816
-[1] 26.21359
-[1] 26.69903
-[1] 27.18447
-[1] 27.6699
-[1] 28.15534
-[1] 28.64078
-[1] 29.12621
-[1] 29.61165
-[1] 30.09709
-[1] 30.58252
-[1] 31.06796
-[1] 31.5534
-[1] 32.03883
-[1] 32.52427
-[1] 33.00971
-[1] 33.49515
-[1] 33.98058
-[1] 34.46602
-[1] 34.95146
-[1] 35.43689
-[1] 35.92233
-[1] 36.40777
-[1] 36.8932
-[1] 37.37864
-[1] 37.86408
-[1] 38.34951
-[1] 38.83495
-[1] 39.32039
-[1] 39.80583
-[1] 40.29126
-[1] 40.7767
-[1] 41.26214
-[1] 41.74757
-[1] 42.23301
-[1] 42.71845
-[1] 43.20388
-[1] 43.68932
-[1] 44.17476
-[1] 44.66019
-[1] 45.14563
-[1] 45.63107
-[1] 46.1165
-[1] 46.60194
-[1] 47.08738
-[1] 47.57282
-[1] 48.05825
-[1] 48.54369
-[1] 49.02913
-[1] 49.51456
-[1] 50
-[1] 50.48544
-[1] 50.97087
-[1] 51.45631
-[1] 51.94175
-[1] 52.42718
-[1] 52.91262
-[1] 53.39806
-[1] 53.8835
-[1] 54.36893
-[1] 54.85437
-[1] 55.33981
-[1] 55.82524
-[1] 56.31068
-[1] 56.79612
-[1] 57.28155
-[1] 57.76699
-[1] 58.25243
-[1] 58.73786
-[1] 59.2233
-[1] 59.70874
-[1] 60.19417
-[1] 60.67961
-[1] 61.16505
-[1] 61.65049
-[1] 62.13592
-[1] 62.62136
-[1] 63.1068
-[1] 63.59223
-[1] 64.07767
-[1] 64.56311
-[1] 65.04854
-[1] 65.53398
-[1] 66.01942
-[1] 66.50485
-[1] 66.99029
-[1] 67.47573
-[1] 67.96117
-[1] 68.4466
-[1] 68.93204
-[1] 69.41748
-[1] 69.90291
-[1] 70.38835
-[1] 70.87379
-[1] 71.35922
-[1] 71.84466
-[1] 72.3301
-[1] 72.81553
-[1] 73.30097
-[1] 73.78641
-[1] 74.27184
-[1] 74.75728
-[1] 75.24272
-[1] 75.72816
-[1] 76.21359
-[1] 76.69903
-[1] 77.18447
-[1] 77.6699
-[1] 78.15534
-[1] 78.64078
-[1] 79.12621
-[1] 79.61165
-[1] 80.09709
-[1] 80.58252
-[1] 81.06796
-[1] 81.5534
-[1] 82.03883
-[1] 82.52427
-[1] 83.00971
-[1] 83.49515
-[1] 83.98058
-[1] 84.46602
-[1] 84.95146
-[1] 85.43689
-[1] 85.92233
-[1] 86.40777
-[1] 86.8932
-[1] 87.37864
-[1] 87.86408
-[1] 88.34951
-[1] 88.83495
-[1] 89.32039
-[1] 89.80583
-[1] 90.29126
-[1] 90.7767
-[1] 91.26214
-[1] 91.74757
-[1] 92.23301
-[1] 92.71845
-[1] 93.20388
-[1] 93.68932
-[1] 94.17476
-[1] 94.66019
-[1] 95.14563
-[1] 95.63107
-[1] 96.1165
-[1] 96.60194
-[1] 97.08738
-[1] 97.57282
-[1] 98.05825
-[1] 98.54369
-[1] 99.02913
-[1] 99.51456
-[1] 100
-> 
-> #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")
-> 
-> proc.time()
-    user   system  elapsed 
-30780.04  1294.12 32076.90