diff --git a/Taylor_Simulated_Campylobacter_environment_light_hum_min_for_rec_delay_original_modified.Rout b/Taylor_Simulated_Campylobacter_environment_light_hum_min_for_rec_delay_original_modified.Rout
deleted file mode 100644
index 09a000653f324a74b750857fb44182f831d0f187..0000000000000000000000000000000000000000
--- a/Taylor_Simulated_Campylobacter_environment_light_hum_min_for_rec_delay_original_modified.Rout
+++ /dev/null
@@ -1,973 +0,0 @@
-
-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