diff --git a/Campylobacter_environment_analysis_subset_variables_hum_rain.Rout b/Campylobacter_environment_analysis_subset_variables_hum_rain.Rout deleted file mode 100644 index 30ce4ae163e5f02f298cf0d76c0ae784e1d6b232..0000000000000000000000000000000000000000 --- a/Campylobacter_environment_analysis_subset_variables_hum_rain.Rout +++ /dev/null @@ -1,491 +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 - -> -> -> -> ## Variable 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<-14 -> width_char<-paste(width) -> -> ######################## Catchment areas ################# -> -> population_df<-read.csv(paste("../../Data_Base/Catchment_areas/LSOA2011Population_SumByLabs.csv",sep="")) -> colnames(population_df)<-c("PostCode","col2","col3","residents") -> -> -> -> merged_lab_all2<-data.frame(All_PC,dates,humidity,max_temp,min_temp,rain,wind) -> colnames(merged_lab_all2)<-c("PostCode","dates","humidity", "max_temp", "min_temp","rain","wind") -> merged_lab_all<-merge(merged_lab_all2, population_df,by="PostCode") -> # -> #merged_lab<-data.frame(dates,humidity,max_temp,min_temp) -> merged.data_all<-read.csv("../../Data_Base/OPIE_data_base/OPIE_environment_Campylobacter.csv") -> -> -> merged.data_all[,2]<-as.Date((as.character(merged.data_all[,2]))) -> merged.data_all2<-merged.data_all -> merged.data_all<-merge(merged.data_all2, population_df,by="PostCode") -> -> -> -> ########################### Weekly average ###################### -> -> week_cases<-function(lab_fac){ -+ lab<-as.character(lab_fac) -+ merged_sub<-subset(merged.data_all,merged.data_all$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(merged_lab_all,merged_lab_all$PostCode==lab) -+ merged_lab_sub2<-merged_lab_sub[order(as.Date(merged_lab_sub$dates, format="%Y-%m-%d")),] -+ -+ -+ -+ 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(merged_lab_all$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,"_test.csv",sep=""), col.names = FALSE, sep = ",",eol = "\n") -> -> -> merged.data_all$Cases<-1 -> merged.data_all_weekly<-merge(merged.data_all, merged_lab_weekly,by=c("PostCode","Date")) -> #merged.data_all_weekly<-merge(merged.data_all_weekly2, merged_cases_weekly,by=c("PostCode","Date")) -> -> -> merged_lab_all<-merged_lab_weekly -> #sort(merged_lab_all) -> merged.data_all<-data.frame(merged.data_all_weekly$PostCode, -+ merged.data_all_weekly$Date, -+ merged.data_all_weekly$Cases, -+ merged.data_all_weekly$mean_humidity, -+ merged.data_all_weekly$mean_max_temp, -+ merged.data_all_weekly$mean_min_temp, -+ merged.data_all_weekly$mean_rain, -+ merged.data_all_weekly$cum_rain, -+ merged.data_all_weekly$mean_wind, -+ merged.data_all_weekly$mean_residents) -> -> colnames(merged.data_all)<-c("PostCode","Date","Cases","humidity","max_temp","min_temp","rain","cum_rain","wind_speed","residents") -> colnames(merged_lab_all)<-c("PostCode","Date","humidity","max_temp","min_temp","rain","cum_rain","wind_speed","residents") -> -> merged.data_all<-merged.data_all[order(as.Date(merged.data_all$Date)),] -> -> ########################### END Weekly average ###################### -> -> -> #PHE_Centre<-merged.data_all$PHE_Centre_Name -> #n_Centre<-length(levels(PHE_Centre)) -> #i_centre<-6 -> #For Dorset only -> #merged.data_PHE<-subset(merged.data_all,merged.data_all$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(merged.data_all,year(merged.data_all$Date)>=1990 & year(merged.data_all$Date)<=2015) -> merged_lab<-subset(merged_lab_all,year(merged_lab_all$Date)>=1990 & year(merged_lab_all$Date)<=2015) -> ################### weekly summary -> -> -> ################### -> -> -> delta_hum<-5 -> delta_temp<-1 -> delta_rain<-1 -> delta_wind<-0.5 -> 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))-0.5,0),max(na.omit(merged_lab_weekly$mean_wind))+0.5,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(mergedab_weekly$cum_rain))+1,by=delta_rain) -Error in na.omit(mergedab_weekly$cum_rain) : - object 'mergedab_weekly' not found -Calls: seq -> seq.default -> na.omit -Execution halted