From 043351011dc1c41e3a9a5dae196de02d89ea7b26 Mon Sep 17 00:00:00 2001
From: "Lo Iacono, Giovanni Dr (School of Vet Med.)" <g.loiacono@surrey.ac.uk>
Date: Wed, 10 May 2023 14:12:52 +0000
Subject: [PATCH] Delete
PAPER_Campylobacter_environment_analysis_subset_variables_hum_max.Rout
---
...ent_analysis_subset_variables_hum_max.Rout | 1002 -----------------
1 file changed, 1002 deletions(-)
delete mode 100644 PAPER_Campylobacter_environment_analysis_subset_variables_hum_max.Rout
diff --git a/PAPER_Campylobacter_environment_analysis_subset_variables_hum_max.Rout b/PAPER_Campylobacter_environment_analysis_subset_variables_hum_max.Rout
deleted file mode 100644
index 2f5845b..0000000
--- a/PAPER_Campylobacter_environment_analysis_subset_variables_hum_max.Rout
+++ /dev/null
@@ -1,1002 +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. 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)
-
-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)
->
->
-> 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)
-Loading required package: zoo
-
-Attaching package: ‘zoo’
-
-The following object is masked from ‘package:timeSeries’:
-
- time<-
-
-The following objects are masked from ‘package:base’:
-
- as.Date, as.Date.numeric
-
-
-Attaching package: ‘xts’
-
-The following objects are masked from ‘package:pastecs’:
-
- first, last
-
->
->
->
-> ## 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))
-+ }
-[1] 0.4694836
-[1] 0.9389671
-[1] 1.408451
-[1] 1.877934
-[1] 2.347418
-[1] 2.816901
-[1] 3.286385
-[1] 3.755869
-[1] 4.225352
-[1] 4.694836
-[1] 5.164319
-[1] 5.633803
-[1] 6.103286
-[1] 6.57277
-[1] 7.042254
-[1] 7.511737
-[1] 7.981221
-[1] 8.450704
-[1] 8.920188
-[1] 9.389671
-[1] 9.859155
-[1] 10.32864
-[1] 10.79812
-[1] 11.26761
-[1] 11.73709
-[1] 12.20657
-[1] 12.67606
-[1] 13.14554
-[1] 13.61502
-[1] 14.08451
-[1] 14.55399
-[1] 15.02347
-[1] 15.49296
-[1] 15.96244
-[1] 16.43192
-[1] 16.90141
-[1] 17.37089
-[1] 17.84038
-[1] 18.30986
-[1] 18.77934
-[1] 19.24883
-[1] 19.71831
-[1] 20.18779
-[1] 20.65728
-[1] 21.12676
-[1] 21.59624
-[1] 22.06573
-[1] 22.53521
-[1] 23.00469
-[1] 23.47418
-[1] 23.94366
-[1] 24.41315
-[1] 24.88263
-[1] 25.35211
-[1] 25.8216
-[1] 26.29108
-[1] 26.76056
-[1] 27.23005
-[1] 27.69953
-[1] 28.16901
-[1] 28.6385
-[1] 29.10798
-[1] 29.57746
-[1] 30.04695
-[1] 30.51643
-[1] 30.98592
-[1] 31.4554
-[1] 31.92488
-[1] 32.39437
-[1] 32.86385
-[1] 33.33333
-[1] 33.80282
-[1] 34.2723
-[1] 34.74178
-[1] 35.21127
-[1] 35.68075
-[1] 36.15023
-[1] 36.61972
-[1] 37.0892
-[1] 37.55869
-[1] 38.02817
-[1] 38.49765
-[1] 38.96714
-[1] 39.43662
-[1] 39.9061
-[1] 40.37559
-[1] 40.84507
-[1] 41.31455
-[1] 41.78404
-[1] 42.25352
-[1] 42.723
-[1] 43.19249
-[1] 43.66197
-[1] 44.13146
-[1] 44.60094
-[1] 45.07042
-[1] 45.53991
-[1] 46.00939
-[1] 46.47887
-[1] 46.94836
-[1] 47.41784
-[1] 47.88732
-[1] 48.35681
-[1] 48.82629
-[1] 49.29577
-[1] 49.76526
-[1] 50.23474
-[1] 50.70423
-[1] 51.17371
-[1] 51.64319
-[1] 52.11268
-[1] 52.58216
-[1] 53.05164
-[1] 53.52113
-[1] 53.99061
-[1] 54.46009
-[1] 54.92958
-[1] 55.39906
-[1] 55.86854
-[1] 56.33803
-[1] 56.80751
-[1] 57.277
-[1] 57.74648
-[1] 58.21596
-[1] 58.68545
-[1] 59.15493
-[1] 59.62441
-[1] 60.0939
-[1] 60.56338
-[1] 61.03286
-[1] 61.50235
-[1] 61.97183
-[1] 62.44131
-[1] 62.9108
-[1] 63.38028
-[1] 63.84977
-[1] 64.31925
-[1] 64.78873
-[1] 65.25822
-[1] 65.7277
-[1] 66.19718
-[1] 66.66667
-[1] 67.13615
-[1] 67.60563
-[1] 68.07512
-[1] 68.5446
-[1] 69.01408
-[1] 69.48357
-[1] 69.95305
-[1] 70.42254
-[1] 70.89202
-[1] 71.3615
-[1] 71.83099
-[1] 72.30047
-[1] 72.76995
-[1] 73.23944
-[1] 73.70892
-[1] 74.1784
-[1] 74.64789
-[1] 75.11737
-[1] 75.58685
-[1] 76.05634
-[1] 76.52582
-[1] 76.99531
-[1] 77.46479
-[1] 77.93427
-[1] 78.40376
-[1] 78.87324
-[1] 79.34272
-[1] 79.81221
-[1] 80.28169
-[1] 80.75117
-[1] 81.22066
-[1] 81.69014
-[1] 82.15962
-[1] 82.62911
-[1] 83.09859
-[1] 83.56808
-[1] 84.03756
-[1] 84.50704
-[1] 84.97653
-[1] 85.44601
-[1] 85.91549
-[1] 86.38498
-[1] 86.85446
-[1] 87.32394
-[1] 87.79343
-[1] 88.26291
-[1] 88.73239
-[1] 89.20188
-[1] 89.67136
-[1] 90.14085
-[1] 90.61033
-[1] 91.07981
-[1] 91.5493
-[1] 92.01878
-[1] 92.48826
-[1] 92.95775
-[1] 93.42723
-[1] 93.89671
-[1] 94.3662
-[1] 94.83568
-[1] 95.30516
-[1] 95.77465
-[1] 96.24413
-[1] 96.71362
-[1] 97.1831
-[1] 97.65258
-[1] 98.12207
-[1] 98.59155
-[1] 99.06103
-[1] 99.53052
-[1] 100
->
-> 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))-1, max(na.omit(merged_lab_weekly$cum_rain))+1,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))
-[1] 4 5 3 2 2
-> #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")
->
-> proc.time()
- user system elapsed
-346.014 8.211 354.697
--
GitLab