Skip to content
Snippets Groups Projects
Commit 0744560c authored by Lo Iacono, Giovanni Dr (School of Vet Med.)'s avatar Lo Iacono, Giovanni Dr (School of Vet Med.)
Browse files

Delete Campylobacter_environment_analysis_subset_variables_hum_rain.Rout

parent 491aa94a
No related branches found
No related tags found
No related merge requests found
Hide whitespace changes
Inline Side-by-side
Campylobacter_environment_analysis_subset_variables_hum_rain.Rout deleted 100644 → 0
+ 0
− 491
View file @ 491aa94a
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
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment