diff --git a/PAPER_Campylobacter_environment_analysis_subset_variables_hum_max.R b/PAPER_Campylobacter_environment_analysis_subset_variables_hum_max.R
deleted file mode 100644
index bb82aa466e0a4b3a9dc88ab9f4e9814adba6b477..0000000000000000000000000000000000000000
--- a/PAPER_Campylobacter_environment_analysis_subset_variables_hum_max.R
+++ /dev/null
@@ -1,725 +0,0 @@
-# The code does look at how the risk of Campylobacter in humans depends on environmental variables. It consider only two variables at a time: Relative humidity and maximum air temeprature
-# It uses original MEDMI data
-
-
-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)
-
-
-list.of.packages <- c("xts")
-new.packages <- list.of.packages[!(list.of.packages %in% installed.packages()[,"Package"])]
-if(length(new.packages)) install.packages(new.packages)
-library(xts)
-
-
-
-## Varaible file
-
-variable<-"humidity"
-variable_df_1<-read.csv(paste("../../Data_Base/OPIE_data_base/",variable,".csv",sep=""))
-humidity<-variable_df_1[,-c(1,2)]
-#dates<-as.Date(as.character(variable_df_1[-c(1,2),2]),format="%d/%m/%Y")
-
-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))
-
-humidity<-humidity[-c(1,2),]
-names(humidity) <- NULL
-humidity<-unlist(c(humidity))
-#wt<-which(humidity<=0)
-#humidity[wt]<-NA
-
-
-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))
-#wt<-which(max_temp>=39)
-#max_temp[wt]<-NA
-
-variable<-"min_air_temp"
-variable_df_3<-read.csv(paste("../../Data_Base/OPIE_data_base/",variable,".csv",sep=""))
-min_temp<-variable_df_3[,-c(1,2)]
-min_temp<-min_temp[-c(1,2),]
-names(min_temp) <- NULL
-min_temp<-unlist(c(min_temp))
-
-
-variable<-"rain"
-variable_df_4<-read.csv(paste("../../Data_Base/OPIE_data_base/",variable,".csv",sep=""))
-rain<-variable_df_4[,-c(1,2)]
-rain<-rain[-c(1,2),]
-names(rain) <- NULL
-rain<-unlist(c(rain))
-
-variable<-"wind_speed"
-variable_df_5<-read.csv(paste("../../Data_Base/OPIE_data_base/",variable,".csv",sep=""))
-wind<-variable_df_5[,-c(1,2)]
-wind<-wind[-c(1,2),]
-names(wind) <- NULL
-wind<-unlist(c(wind))
-
-
-width<-30
- width_char<-paste(width)
-
- ######################## 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")
- PC_df<-data.frame(All_PC,as.Date(dates))
- colnames(PC_df)<-c("PostCode","Date")
-
- Post_Codes_df<-merge(PC_df,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<-Post_Codes_df$lat
- day_of_the_year<-yday(as.Date(Post_Codes_df$Date))
-
- daylength_int1<-mapply(daylength, latitude, day_of_the_year)
- daylength_df<-data.frame(Post_Codes_df$PostCode,as.Date(Post_Codes_df$Date),daylength_int1)
- colnames(daylength_df)<-c("PostCode","Date","daylength")
-
-
-
-######################## Catchment areas #################
- catchment_population_df<-read.csv(paste("../../Data_Base/Catchment_areas/Sum_ByLab_1987_2016.csv",sep=""))
- colnames(catchment_population_df)<-c("PostCode","col2","col3",
- "residents_1987",
- "residents_1988",
- "residents_1989",
- "residents_1990",
- "residents_1991",
- "residents_1992",
- "residents_1993",
- "residents_1994",
- "residents_1995",
- "residents_1996",
- "residents_1997",
- "residents_1998",
- "residents_1999",
- "residents_2000",
- "residents_2001",
- "residents_2002",
- "residents_2003",
- "residents_2004",
- "residents_2005",
- "residents_2006",
- "residents_2007",
- "residents_2008",
- "residents_2009",
- "residents_2010",
- "residents_2011",
- "residents_2012",
- "residents_2013",
- "residents_2014",
- "residents_2015",
- "residents_2016",
- "residents_2011bis")
-
- diagnostic_laboratory_info_df<-read.csv(paste("../../Data_Base/Catchment_areas/Lab_information.csv",sep=""))
- diagnostic_laboratory_info_df<-diagnostic_laboratory_info_df[,-1]
-
-
- ## This create a database for all laboratories postcode, even if there were zero cases
- diagnostic_laboratory_df2<-data.frame(All_PC,dates,humidity,max_temp,min_temp,rain,wind)
- colnames(diagnostic_laboratory_df2)<-c("PostCode","Date","humidity", "max_temp", "min_temp","rain","wind")
-
-
- ###################################
-
-
-
- starting_year_catchment_population_df<-1987
- end_year_catchment_population_df<-2016
- numb_years<-length(seq(starting_year_catchment_population_df:end_year_catchment_population_df))
- sequence_years<-seq(starting_year_catchment_population_df,end_year_catchment_population_df)
-
- diagnostic_laboratory_df2$year<-year(diagnostic_laboratory_df2$Date)
-
- catchment_population_df2<-catchment_population_df[,-c(1:3,34)]
-
-
-
- catchment_population_df3<-data.frame(
- rep(catchment_population_df$PostCode,times=numb_years),
- rep(sequence_years,each=length(unique(catchment_population_df$PostCode))),
- stack(catchment_population_df2)
- )
- colnames(catchment_population_df3)<-c("PostCode","year","residents","ind")
-
- catchment_population_df3$year<-as.factor(catchment_population_df3$year)
- diagnostic_laboratory_df2$year<-as.factor(diagnostic_laboratory_df2$year)
-
- catchment_population_df3$PostCode<-as.factor(catchment_population_df3$PostCode)
- diagnostic_laboratory_df2$PostCode<-as.factor(diagnostic_laboratory_df2$PostCode)
- diagnostic_laboratory_df4<-merge(diagnostic_laboratory_df2, catchment_population_df3,by=c("PostCode","year"))
-
-
-
- daylength_df$PostCode<-as.factor(daylength_df$PostCode)
- daylength_df$Date<-as.factor(daylength_df$Date)
- diagnostic_laboratory_df4$PostCode<-as.factor(diagnostic_laboratory_df4$PostCode)
- diagnostic_laboratory_df4$Date<-as.factor(diagnostic_laboratory_df4$Date)
-
- diagnostic_laboratory_df<-merge(diagnostic_laboratory_df4[,c(1,3:9)], daylength_df,by=c("PostCode","Date"))
-
-
-
-
- ## This create a database for all postocode where there was at least one case
- Campylobacter_cases<-read.csv("../../Data_Base/Cases/Campylobacter_Simulated_Data_final_UK.csv")
- Campylobacter_cases<-Campylobacter_cases[,-1]
- colnames(Campylobacter_cases)<-c("Cases","Adjusted","Spec_Date", "PostCode","Date","Simulation")
-
- Campylobacter_cases$PostCode<-stri_replace_all_fixed(Campylobacter_cases$PostCode, " ", "")
- Campylobacter_cases$Date<-as.Date((as.character(Campylobacter_cases$Date)))
- Campylobacter_cases$year<-as.factor((as.character(year(Campylobacter_cases$Date))))
- Campylobacter_cases$PostCode<-as.factor(Campylobacter_cases$PostCode)
-
- Campylobacter_cases_df<-merge(Campylobacter_cases, catchment_population_df3,by=c("PostCode","year"))
-
-
-
-########################### Weekly average ######################
-
-week_cases<-function(lab_fac){
- lab<-as.character(lab_fac)
- merged_sub<-subset(Campylobacter_cases_df,Campylobacter_cases_df$PostCode==lab)
-
- if (length(merged_sub[,1])!=0){
-
- merged_sub2<-merged_sub[order(as.Date(merged_sub$Date, format="%Y-%m-%d")),]
-
- ep <- endpoints(as.xts(merged_sub2$Date),'weeks')
- cum_Cases<-period.apply(merged_sub2$Cases, INDEX=ep, FUN=function(x) sum(na.omit(x,na.rm=TRUE)))
-
-
-
- PC<-rep(as.character(unique(merged_sub$PostCode)),times=length(ep)-1)
- merged_weekly<-data.frame(PC,merged_sub$Date[ep],cum_Cases)
- return(merged_weekly)}
-}
-
-
-week_PC<-function(lab_fac){
- lab<-as.character(lab_fac)
- merged_lab_sub<-subset(diagnostic_laboratory_df,diagnostic_laboratory_df$PostCode==lab)
- merged_lab_sub2<-merged_lab_sub[order(as.Date(merged_lab_sub$Date)),]
-
-# ep <- endpoints(as.xts(merged_lab_sub2$dates),'weeks')
-
- #mean_humidity<-period.apply(merged_lab_sub2$humidity, INDEX=ep, FUN=function(x) mean(na.omit(x,na.rm=TRUE)))
- #mean_max_temp<-period.apply(merged_lab_sub2$max_temp, INDEX=ep, FUN=function(x) mean(na.omit(x,na.rm=TRUE)))
- #mean_min_temp<-period.apply(merged_lab_sub2$min_temp, INDEX=ep, FUN=function(x) mean(na.omit(x,na.rm=TRUE)))
- #mean_rain<-period.apply(merged_lab_sub2$rain, INDEX=ep, FUN=function(x) mean(na.omit(x,na.rm=TRUE)))
- #cum_rain<-period.apply(merged_lab_sub2$rain, INDEX=ep, FUN=function(x) sum(na.omit(x,na.rm=TRUE)))
- #mean_wind<-period.apply(merged_lab_sub2$wind, INDEX=ep, FUN=function(x) mean(na.omit(x,na.rm=TRUE)))
- #mean_residents<-period.apply(merged_lab_sub2$residents, INDEX=ep, FUN=function(x) mean(na.omit(x,na.rm=TRUE)))
- #PC<-rep(as.character(unique(merged_lab_sub$PostCode)),times=length(ep)-1)
-
-
- mean_humidity<-rollmean(merged_lab_sub2$humidity,width)
- mean_max_temp<-rollmean(merged_lab_sub2$max_temp, width)
- mean_min_temp<-rollmean(merged_lab_sub2$min_temp, width)
- mean_rain<-rollmean(merged_lab_sub2$rain, width)
- cum_rain<-rollsum(merged_lab_sub2$rain, width)
- mean_wind<-rollmean(merged_lab_sub2$wind,width)
- mean_residents<-rollmean(merged_lab_sub2$residents, width)
-
- PC<-rep(as.character(unique(merged_lab_sub$PostCode)),times=length(mean_residents))
- ep<-seq(width,length(mean_residents)+width-1)
-
- merged_lab_weekly<-data.frame(PC,dates[ep],mean_humidity,mean_max_temp,mean_min_temp,mean_rain,cum_rain,mean_wind,mean_residents)
- return(merged_lab_weekly)
-}
-
-merged_lab_weekly<-c()
-#merged_cases_weekly<-c()
-
-index_PC<-unique(diagnostic_laboratory_df$PostCode)
-for (i in c(1:length(index_PC))){
- merged_lab_weekly<-rbind(merged_lab_weekly,lapply(index_PC[i], week_PC)[[1]])
-
- #merged_cases_weekly<-rbind(merged_cases_weekly,lapply(index_PC[i], week_cases)[[1]])
- print(100*i/length(index_PC))
-}
-
-colnames(merged_lab_weekly)<-c("PostCode","Date","mean_humidity","mean_max_temp","mean_min_temp","mean_rain","cum_rain","mean_wind","mean_residents")
-write.table(merged_lab_weekly,paste("../../Data_Base/OPIE_data_base/Laboratory_",width_char,"_original_MEDMI.csv",sep=""), col.names = FALSE, sep = ",",eol = "\n")
-
-
-Campylobacter_cases_df$Cases<-1#WARNING CHECK THIS
-Campylobacter_cases_df_weekly<-merge(Campylobacter_cases_df, merged_lab_weekly,by=c("PostCode","Date"))
-#Campylobacter_cases_df_weekly<-merge(Campylobacter_cases_df_weekly2, merged_cases_weekly,by=c("PostCode","Date"))
-
-
-diagnostic_laboratory_df<-merged_lab_weekly
-#sort(diagnostic_laboratory_df)
-Campylobacter_cases_df<-data.frame(Campylobacter_cases_df_weekly$PostCode,
- Campylobacter_cases_df_weekly$Date,
- Campylobacter_cases_df_weekly$Cases,
- Campylobacter_cases_df_weekly$mean_humidity,
- Campylobacter_cases_df_weekly$mean_max_temp,
- Campylobacter_cases_df_weekly$mean_min_temp,
- Campylobacter_cases_df_weekly$mean_rain,
- Campylobacter_cases_df_weekly$cum_rain,
- Campylobacter_cases_df_weekly$mean_wind,
- Campylobacter_cases_df_weekly$mean_residents)
-
-colnames(Campylobacter_cases_df)<-c("PostCode","Date","Cases","humidity","max_temp","min_temp","rain","cum_rain","wind_speed","residents")
-colnames(diagnostic_laboratory_df)<-c("PostCode","Date","humidity","max_temp","min_temp","rain","cum_rain","wind_speed","residents")
-
-Campylobacter_cases_df<-Campylobacter_cases_df[order(as.Date(Campylobacter_cases_df$Date)),]
-
-########################### END Weekly average ######################
-
-
-#PHE_Centre<-Campylobacter_cases_df$PHE_Centre_Name
-#n_Centre<-length(levels(PHE_Centre))
-#i_centre<-6
-#For Dorset only
-#merged.data_PHE<-subset(Campylobacter_cases_df,Campylobacter_cases_df$PHE_Centre_Name==levels(PHE_Centre)[i_centre])
-#merged.data<-subset(merged.data_PHE,year(merged.data_PHE$Date)>=2000 & year(merged.data_PHE$Date)<2016)
-merged.data<-subset(Campylobacter_cases_df,year(Campylobacter_cases_df$Date)>=1990 & year(Campylobacter_cases_df$Date)<2015)
-merged_lab<-subset(diagnostic_laboratory_df,year(diagnostic_laboratory_df$Date)>=1990 & year(diagnostic_laboratory_df$Date)<2015)
-################### weekly summary
-
-
-###################
-
-
-delta_hum<-5
-delta_temp<-1
-delta_rain<-1
-delta_cum_rain<-20
-delta_wind<-2
-breaks_hum<-seq(max(min(na.omit(merged_lab_weekly$mean_humidity))-10,0),max(na.omit(merged_lab_weekly$mean_humidity))+10,by=delta_hum) #i
-breaks_min_temp<-seq(min(na.omit(merged_lab_weekly$mean_min_temp))-2, max(na.omit(merged_lab_weekly$mean_min_temp))+2,by=delta_temp)
-breaks_max_temp<-seq(min(na.omit(merged_lab_weekly$mean_max_temp))-2, max(na.omit(merged_lab_weekly$mean_max_temp))+2,by=delta_temp)
-breaks_rain<-seq(min(na.omit(merged_lab_weekly$mean_rain))-1, max(na.omit(merged_lab_weekly$mean_rain))+1,by=delta_rain)
-breaks_wind<-seq(max(min(na.omit(merged_lab_weekly$mean_wind))-2,0),max(na.omit(merged_lab_weekly$mean_wind))+2,by=delta_wind)
-#breaks_wind<-seq(max(min(na.omit(wind))-2,0),max(na.omit(wind))+2,by=delta_wind) #WARNING
-breaks_cum_rain<-seq(min(na.omit(merged_lab_weekly$cum_rain))-20, max(na.omit(merged_lab_weekly$cum_rain))+20,by=delta_cum_rain)
-breaks_mean_temp<-seq(min(na.omit(min_temp))-2,max(na.omit(max_temp))+2,by=delta_temp)
-
-
-
-# First find right domain where the values have no NA
-i_hum_min<-min(which(breaks_hum>=min(na.omit(merged.data$humidity))))
-i_hum_max<-max(which(breaks_hum<=max(na.omit(merged.data$humidity))))
-
-i_min_temp_min<-max(min(which(breaks_min_temp>=min(na.omit(merged.data$min_temp))))-1,1)
-i_min_temp_max<-max(max(which(breaks_min_temp<=max(na.omit(merged.data$min_temp))))+1,1)
-
-i_max_temp_min<-max(min(which(breaks_max_temp>=min(na.omit(merged.data$max_temp))))-1,1)
-i_max_temp_max<-max(max(which(breaks_max_temp<=max(na.omit(merged.data$max_temp))))+1,1)
-
-i_rain_min<-max(min(which(breaks_rain>=min(na.omit(merged.data$rain))))-1,1)
-i_rain_max<-max(which(breaks_rain<=max(na.omit(merged.data$rain))))+1
-
-i_cum_rain_min<-max(min(which(breaks_cum_rain>=min(na.omit(merged.data$cum_rain))))-1,1)
-i_cum_rain_max<-max(which(breaks_cum_rain<=max(na.omit(merged.data$cum_rain))))+1
-
-
-i_wind_min<-max(min(which(breaks_wind>=min(na.omit(merged.data$wind))))-1,1)
-i_wind_max<-max(which(breaks_wind<=max(na.omit(merged.data$wind))))+1
-
-print(c(i_hum_min,i_min_temp_min ,i_max_temp_min, i_rain_min,i_wind_min))
-#print(c(i_hum_min,i_hum_max,i_min_temp_min ,i_min_temp_max,i_max_temp_min,i_max_temp_max, i_rain_min,i_rain_max,i_wind_min,i_wind_max))
-
-
-
-
-
-
-
-############################# General variables ###########################
-
-
-
-
-
-
-variable_x<-"max_air_temp"
-#
-variable<-"cum_rain"
-#variable<-"rain"
-#variable<-"wind"
-#variable<-"humidity"
-var_x_loc_df<-c()
-
-var_x_loc_df<-unname(var_x_loc_df)
-
-if (variable=="max_air_temp"){
-
- breaks_var<-breaks_max_temp
- i_var_min<-i_max_temp_min
- i_var_max<-i_max_temp_max
- merged.data_var<-merged.data$max_temp
- merged_lab_var<-merged_lab$max_temp
-}
-if (variable_x=="max_air_temp"){
-
- i_var_x_min<-i_max_temp_min
- i_var_x_max<-i_max_temp_max
- breaks_var_x<-breaks_max_temp
-}
-
-if (variable=="min_air_temp"){
-
- breaks_var<-breaks_min_temp
- i_var_min<-i_min_temp_min
- i_var_max<-i_min_temp_max
- merged.data_var<-merged.data$min_temp
- merged_lab_var<-merged_lab$min_temp
-}
-if (variable_x=="min_air_temp"){
-
- i_var_x_min<-i_min_temp_min
- i_var_x_max<-i_min_temp_max
- breaks_var_x<-breaks_min_temp
-}
-
-
-if (variable=="humidity"){
-
- breaks_var<-breaks_hum
- i_var_min<-i_hum_min
- i_var_max<-i_hum_max
- merged.data_var<-merged.data$humidity
- merged_lab_var<-merged_lab$humidity
-}
-if (variable_x=="humidity"){
-
- i_var_x_min<-i_hum_min
- i_var_x_max<-i_hum_max
- breaks_var_x<-breaks_hum
-}
-
-
-
-if (variable=="mean_temp"){
-
- breaks_var<-breaks_mean_temp
- i_var_min<-i_mean_temp_min
- i_var_max<-i_mean_temp_max
- merged.data_var<-merged.data$mean_temp
- merged_lab_var<-merged_lab$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=="rain"){
-
- breaks_var<-breaks_rain
- i_var_min<-i_rain_min
- i_var_max<-i_rain_max
- merged.data_var<-merged.data$rain
- merged_lab_var<-merged_lab$rain
-}
-if (variable_x=="rain"){
-
- i_var_x_min<-i_rain_min
- i_var_x_max<-i_rain_max
- breaks_var_x<-breaks_rain
-}
-
-if (variable=="cum_rain"){
-
- breaks_var<-breaks_cum_rain
- i_var_min<-i_cum_rain_min
- i_var_max<-i_cum_rain_max
- merged.data_var<-merged.data$cum_rain
- merged_lab_var<-merged_lab$cum_rain
-}
-if (variable_x=="cum_rain"){
-
- i_var_x_min<-i_rain_min
- i_var_x_max<-i_rain_max
- breaks_var_x<-breaks_cum_rain
-}
-
-if (variable=="wind"){
-
- breaks_var<-breaks_wind
- i_var_min<-i_wind_min
- i_var_max<-i_wind_max
- merged.data_var<-merged.data$wind
- merged_lab_var<-merged_lab$wind
-}
-if (variable_x=="wind"){
-
- i_var_x_min<-i_wind_min
- i_var_x_max<-i_wind_max
- breaks_var_x<-breaks_wind
-}
-
-
-Yt_var_x<-function(i_var)
-{
-
- Yt1<-subset(merged.data,merged.data_var>=breaks_var[i_var] & merged.data_var<breaks_var[i_var+1])
- Yt2<-Yt1
- #Yt2<-subset(Yt1,Yt1$max_temp>=breaks_max_temp[i_max_temp] & Yt1$max_temp<=breaks_max_temp[i_max_temp+1])
- #Yt3<-subset(Yt2,Yt2$rain>=breaks_rain[i_rain] & Yt2$rain<=breaks_rain[i_rain+1])
- #Yt4<-subset(Yt3,Yt3$wind>=breaks_wind[i_wind] & Yt3$wind<=breaks_wind[i_wind+1])
- #if (length(Yt2[,1])>0){
-
- return(as.list(Yt2))
- # return(Yt4)
- #} else {
- # Yt2<-NA
- #return(Yt2)
-
- #}
-
-}
-
-
-Tot_var_x<-function(i_var)
-{
-
- Yt1<-subset(merged_lab,merged_lab_var>=breaks_var[i_var] & merged_lab_var<breaks_var[i_var+1])
- Yt2<-Yt1
- #Yt2<-subset(Yt1,Yt1$max_temp>=breaks_max_temp[i_max_temp] & Yt1$max_temp<=breaks_max_temp[i_max_temp+1])
- #Yt3<-subset(Yt2,Yt2$rain>=breaks_rain[i_rain] & Yt2$rain<=breaks_rain[i_rain+1])
- #Yt4<-subset(Yt3,Yt3$wind>=breaks_wind[i_wind] & Yt3$wind<=breaks_wind[i_wind+1])
- #if (length(Yt2[,1])>0){
-
- return(as.list(Yt2))
- # return(Yt4)
- #} else {
- # Yt2<-NA
- #return(Yt2)
-
- #}
-
-}
-var_x_loc_df<-c(0)
-n_seas<-1
-
-
-#colnames(var_x_loc_df)<-c(variable_x,"prop","prevalence","month",variable,"counts","counts_tot","residents","residents_tot")
-
-for (i_var in c(i_var_min:i_var_max))
-{
- for (i in c(1:n_seas))
- {
-
- n_months<-12/n_seas
- #if (is.na(Yt_min_temp(i_hum)$min_temp)==FALSE){
-
- wt<-which(month(Yt_var_x(i_var)$Date)>(i-1)*n_months & month(Yt_var_x(i_var)$Date)<=i*n_months)
- wt_tot<-which(month(Tot_var_x(i_var)$Date)>(i-1)*n_months & month(Tot_var_x(i_var)$Date)<=i*n_months)
-
- if (variable_x=="min_air_temp"){
- Campylobacter_var_x<-Yt_var_x(i_var)$min_temp[wt]
- var_x_tot<-Tot_var_x(i_var)$min_temp[wt_tot]
-
- }
- if (variable_x=="max_air_temp"){
- Campylobacter_var_x<-Yt_var_x(i_var)$max_temp[wt]
- var_x_tot<-Tot_var_x(i_var)$max_temp[wt_tot]
-
- }
-
- if (variable_x=="mean_temp"){
- Campylobacter_var_x<-Yt_var_x(i_var)$mean_temp[wt]
- var_x_tot<-Tot_var_x(i_var)$mean_temp[wt_tot]
-
- }
-
-
-
- if (variable_x=="humidity"){
- Campylobacter_var_x<-Yt_var_x(i_var)$humidity[wt]
- var_x_tot<-Tot_var_x(i_var)$humidity[wt_tot]
-
- }
- if (variable_x=="wind"){
- Campylobacter_var_x<-Yt_var_x(i_var)$wind[wt]
- var_x_tot<-Tot_var_x(i_var)$wind[wt_tot]
-
- }
- if (variable_x=="rain"){
- Campylobacter_var_x<-Yt_var_x(i_var)$rain[wt]
- var_x_tot<-Tot_var_x(i_var)$rain[wt_tot]
-
- }
-
- n_var_x_tot<-hist(var_x_tot,breaks=breaks_var_x)
- n_var_x<-hist((as.numeric(as.character(Campylobacter_var_x))),breaks=breaks_var_x)
-
- residents<-rep(0,times=length(n_var_x$counts))
- residents_tot<-rep(0,times=length(n_var_x$counts))
-
- wt<-which(n_var_x_tot$counts!=0)
- wt2<-which(n_var_x$counts!=0)
-
-
- if (variable_x=="min_air_temp"){
- if(length(wt)>0){
- for (j in c(1:(length(wt)))){
-
- ww<-which(Tot_var_x(i_var)$min_temp>=n_var_x_tot$breaks[wt[j]] &
- Tot_var_x(i_var)$min_temp<=n_var_x_tot$breaks[wt[j]]+delta_temp )
- residents_tot[wt[j]]<-sum(as.numeric(Tot_var_x(i_var)$residents[ww]))
- }
- if(length(wt2)>0){
- for (j in c(1:(length(wt2)))){
-
- ww<-which(Yt_var_x(i_var)$min_temp>=n_var_x$breaks[wt2[j]] & Yt_var_x(i_var)$min_temp<=n_var_x$breaks[wt2[j]]+delta_temp)
- residents[wt2[j]]<-sum(Yt_var_x(i_var)$residents[ww])
- }
- }
-
- }
-
-
- }
- if (variable_x=="max_air_temp"){
-
- if(length(wt)>0){
- for (j in c(1:(length(wt)))){
-
- ww<-which(Tot_var_x(i_var)$max_temp>=n_var_x_tot$breaks[wt[j]] &
- Tot_var_x(i_var)$max_temp<=n_var_x_tot$breaks[wt[j]]+delta_temp )
- residents_tot[wt[j]]<-sum(as.numeric(Tot_var_x(i_var)$residents[ww]))
- }
- if(length(wt2)>0){
- for (j in c(1:(length(wt2)))){
-
- ww<-which(Yt_var_x(i_var)$max_temp>=n_var_x$breaks[wt2[j]] & Yt_var_x(i_var)$max_temp<=n_var_x$breaks[wt2[j]]+delta_temp)
- residents[wt2[j]]<-sum(Yt_var_x(i_var)$residents[ww])
- }
- }
-
- }
-
- }
-
- if (variable_x=="mean_temp"){
- if(length(wt)>0){
- for (j in c(1:(length(wt)))){
-
- ww<-which(Tot_var_x(i_var)$mean_temp>=n_var_x_tot$breaks[wt[j]] &
- Tot_var_x(i_var)$mean_temp<=n_var_x_tot$breaks[wt[j]]+delta_temp )
- residents_tot[wt[j]]<-sum(as.numeric(Tot_var_x(i_var)$residents[ww]))
- }
- if(length(wt2)>0){
- for (j in c(1:(length(wt2)))){
-
- ww<-which(Yt_var_x(i_var)$mean_temp>=n_var_x$breaks[wt2[j]] & Yt_var_x(i_var)$mean_temp<=n_var_x$breaks[wt2[j]]+delta_temp)
- residents[wt2[j]]<-sum(Yt_var_x(i_var)$residents[ww])
- }
- }
-
- }
- }
-
-
- if (variable_x=="humidity"){
- if(length(wt)>0){
- for (j in c(1:(length(wt)))){
-
- ww<-which(Tot_var_x(i_var)$hum>=n_var_x_tot$breaks[wt[j]] &
- Tot_var_x(i_var)$hum<=n_var_x_tot$breaks[wt[j]]+delta_hum )
- residents_tot[wt[j]]<-sum(as.numeric(Tot_var_x(i_var)$residents[ww]))
- }
- if(length(wt2)>0){
- for (j in c(1:(length(wt2)))){
-
- ww<-which(Yt_var_x(i_var)$hum>=n_var_x$breaks[wt2[j]] & Yt_var_x(i_var)$hum<=n_var_x$breaks[wt2[j]]+delta_hum)
- residents[wt2[j]]<-sum(Yt_var_x(i_var)$residents[ww])
- }
- }
-
- }
- }
- if (variable_x=="wind"){
- if(length(wt)>0){
- for (j in c(1:(length(wt)))){
-
- ww<-which(Tot_var_x(i_var)$wind>=n_var_x_tot$breaks[wt[j]] &
- Tot_var_x(i_var)$wind<=n_var_x_tot$breaks[wt[j]]+delta_wind)
- residents_tot[wt[j]]<-sum(as.numeric(Tot_var_x(i_var)$residents[ww]))
- }
- if(length(wt2)>0){
- for (j in c(1:(length(wt2)))){
-
- ww<-which(Yt_var_x(i_var)$wind>=n_var_x$breaks[wt2[j]] & Yt_var_x(i_var)$wind<=n_var_x$breaks[wt2[j]]+delta_wind)
- residents[wt2[j]]<-sum(Yt_var_x(i_var)$residents[ww])
- }
- }
-
- }
-
- }
- if (variable_x=="rain"){
-
- if(length(wt)>0){
- for (j in c(1:(length(wt)))){
-
- ww<-which(Tot_var_x(i_var)$rain>=n_var_x_tot$breaks[wt[j]] &
- Tot_var_x(i_var)$rain<=n_var_x_tot$breaks[wt[j]]+delta_rain)
- residents_tot[wt[j]]<-sum(as.numeric(Tot_var_x(i_var)$residents[ww]))
- }
- if(length(wt2)>0){
- for (j in c(1:(length(wt2)))){
-
- ww<-which(Yt_var_x(i_var)$rain>=n_var_x$breaks[wt2[j]] & Yt_var_x(i_var)$rain<=n_var_x$breaks[wt2[j]]+delta_rain)
- residents[wt2[j]]<-sum(Yt_var_x(i_var)$residents[ww])
- }
- }
-
- }
- }
-
-
- # if(length(residents)>0){
- data_df<-data.frame(n_var_x$mids,n_var_x$counts/n_var_x_tot$counts,(n_var_x$counts)/(residents_tot),i,breaks_var[i_var],n_var_x$counts,n_var_x_tot$counts,residents,residents_tot)
-
- colnames(data_df)<-c(variable_x,"prop","incidence","month",variable,"counts","counts_tot","residents","residents_tot")
- var_x_loc_df<-rbind(var_x_loc_df,data_df)
- colnames(var_x_loc_df)<-c(variable_x,"prop","incidence","month",variable,"counts","counts_tot","residents","residents_tot")
-#}
- #}
-
-
- }
-}
-
-
-write.table(var_x_loc_df,paste("../../Data_Base/OPIE_data_base/",variable,"_",variable_x,"_",width_char,"_original_MEDMI.csv",sep=""), col.names = FALSE, sep = ",",eol = "\n")