--- title: "Temporal, spatial and household dynamics of typhoid fever in Kasese district, Uganda " author: "ADRIAN MUWONGE" date: "12 October 2017" output: html_document --- ```{r setup, include=FALSE} knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning=FALSE) library(pander) library(knitr) library(forecast) library(ggplot2) library(dplyr) library(plyr) library(reshape) library(gridExtra) library(lubridate) library(matrixStats) library(xts) library(nlme) library(tidyverse) #install.packages("INLA", repos=c(getOption("repos"), INLA="https://inla.r-inla-download.org/R/stable"), dep=TRUE) # library(INLA) # library(inlabru) library(imputeTS) #START----- multiplot <- function(..., plotlist=NULL, file, cols=1, layout=NULL) { library(grid) # Make a list from the ... arguments and plotlist plots <- c(list(...), plotlist) numPlots = length(plots) # If layout is NULL, then use 'cols' to determine layout if (is.null(layout)) { # Make the panel # ncol: Number of columns of plots # nrow: Number of rows needed, calculated from # of cols layout <- matrix(seq(1, cols * ceiling(numPlots/cols)), ncol = cols, nrow = ceiling(numPlots/cols)) } if (numPlots==1) { print(plots[[1]]) } else { # Set up the page grid.newpage() pushViewport(viewport(layout = grid.layout(nrow(layout), ncol(layout)))) # Make each plot, in the correct location for (i in 1:numPlots) { # Get the i,j matrix positions of the regions that contain this subplot matchidx <- as.data.frame(which(layout == i, arr.ind = TRUE)) print(plots[[i]], vp = viewport(layout.pos.row = matchidx$row, layout.pos.col = matchidx$col)) } } } # END----- ``` ```{r, echo=TRUE} ##••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••• ### Data analysis: Temporal, spatial and household socio-economic dynamics of typhoid: Towards forecasting outbreaks in Kasese district, Uganda #•••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••••• # READ IN DATABASES S1 & S2 MOH_DB<- read.csv("MOH_NAT_DIS_DB.csv",sep = ",",header = T) ##*************************************************************** #### FORECASTING OUTBREAKS AT NATIONAL AND DISTRICT LEVEL ##*************************************************************** MOH_2016 <- read.csv("MOH_2016.csv",sep = ",",header = T) MOH_2016_UG<-MOH_2016[,5,] MOH_2016_KAS<-MOH_2016[,4,] x <- read.csv("MOH_Kasese.csv",sep = ",",header = T) bx <- x[,1:5] FULL_MOH<-join(bx,MOH_2016, type="full")# USE PLYR LIBRARY ``` \newpage ```{r, echo=TRUE} ##•••••••••••••••••••••••••••••••••••••• ### Case control data analysis: #••••••••••••••••••••••••••••••••••••••• ## functions for univariable analysis # function creates a table with OR, CI and p value. x required is summary(model) mmt <- function(y){ x <-summary(y) bla <-strsplit(as.character(x$call)[[2]],split='~', fixed=TRUE)[[1]][2] variable <-gsub(' ','',strsplit(bla,split='+', fixed=TRUE)[[1]][1]) #OddsRatio <- round(exp(x$coefficients[,1][2:nrow(x$coefficients)]),2) OddsRatio <- round(exp(x$coefficients[,1]),2) LL <-round(exp(x$coefficients[,1] - 1.96* x$coefficients[,3]),2) UL <-round(exp(x$coefficients[,1] + 1.96* x$coefficients[,3]),2) pv <- round(x$coefficients[,5],3) #pvalue <-c(round(an1["Pr(>Chisq)"][2,],3),rep("",nrow(x$coefficients)-1)) test <- data.frame(Variable=c(variable, 1:nrow(x$coefficients)), OR=c("",OddsRatio), LL=c("",LL), UL=c("",UL), pv=c(x$waldtest["pvalue"],pv), filter=c("yes", 1:nrow(x$coefficients))) test$Variable <- c(variable, rownames(test)[2:length(rownames(test))]) #test$NAs <- c(na_count[variable], rep("",nrow(x$coefficients)-1)) rownames(test) <- NULL return(test) } # C(Using Database S2 AND CASE CONTROL DATA SETS) # read data in #clg <-read_csv("cond_logreg_data.csv") load("cond_logreg_data.RData") # create matching variable clg$match <- substr(clg$Questionnairenumber, 1,5) # clean dataset to use in regression model clg$Sex <- as.factor(as.character(clg$Sex)) clg$Levelofeducation[is.na(clg$Levelofeducation)] <- "missing" clg$Levelofeducation <- as.factor(as.character(clg$Levelofeducation)) clg$Householdmembers_cut <- recode(as.factor(clg$Householdmembers!="1-5"), "TRUE"="More than 5", "FALSE"= "1-5") clg$Doyoupreheatleftoverfood[is.na(clg$Doyoupreheatleftoverfood)] <- "No preheated food" # change NAs to no left over food clg$Washhandsaftertoiletorlatrine<- as.factor(clg$Washhandsaftertoiletorlatrine) clg$Evenonurination<- as.factor(clg$Evenonurination) ### not likely to be independent from presvious question clg$Soapandwatertowashhands<- as.factor(clg$Soapandwatertowashhands) ### not likely to be independent from presvious question clg$Usetoiletpaper<- as.factor(clg$Usetoiletpaper) ### clg$Handlechildrenswaste<- as.factor(clg$Handlechildrenswaste) ### clg$Wherefoodisprepared<- as.factor(clg$Wherefoodisprepared) ### clg$Eatreadytoeatfoods<- as.factor(clg$Eatreadytoeatfoods) ### clg$Howdoyoupreparethevegetables[is.na(clg$Howdoyoupreparethevegetables)]<-"missing" clg$Howdoyoupreparethevegetables<- as.factor(clg$Howdoyoupreparethevegetables) ### clg$Howdoyoupreparethevegetables[clg$Howdoyoupreparethevegetables=="1"]<-"2" clg$Sharethesamebowlwhileeatingfoodesp<- as.factor(clg$Sharethesamebowlwhileeatingfoodesp) ### clg$Doyouhaveanimals<- as.factor(clg$Doyouhaveanimals) ### clg$Feverishconditions<- as.factor(clg$Feverishconditions) ### clg$Diarrhoea<- as.factor(clg$Diarrhoea) ### clg$Constipation<- as.factor(clg$Constipation) ### clg$Malaise<- as.factor(clg$Malaise) ### clg$Abdominalpains<- as.factor(clg$Abdominalpains) ### clg$Swollenpainfulabdomen<- as.factor(clg$Swollenpainfulabdomen) ### clg$Vomiting<- as.factor(clg$Vomiting) ### clg$Painfulabdomennotswollen<- as.factor(clg$Painfulabdomennotswollen) ### clg$Blood<- as.factor(clg$Blood) ### clg$Stool<- as.factor(clg$Stool) ### clg$Urine<- as.factor(clg$Urine) ### clg$Saliva<- as.factor(clg$Saliva) clg$Treatdrinkingwaterpractice<- as.factor(clg$Treatdrinkingwaterpractice) clg$Dontshareutensilswiththevictim<- as.factor(clg$Dontshareutensilswiththevictim) clg$Eatfoodwhilehot<- as.factor(clg$Eatfoodwhilehot) clg$Dontsharethesamebowloffoodwiththevictim<- as.factor(clg$Dontsharethesamebowloffoodwiththevictim)###only on positive clg$Noadvicegiven<- as.factor(clg$Noadvicegiven) clg$Watersourcesharedbyanimals<- as.factor(clg$Watersourcesharedbyanimals) clg$Humanscandirectlycontaminatewateratsource<- as.factor(clg$Humanscandirectlycontaminatewateratsource) clg$Humansindirectcontactwithatwaterwhilefetching<- as.factor(clg$Humansindirectcontactwithatwaterwhilefetching) clg$Anypitlatrinesclosetowatersource[clg$Anypitlatrinesclosetowatersource==""]<-"No" clg$Anypitlatrinesclosetowatersource<- as.factor(clg$Anypitlatrinesclosetowatersource) clg$Freeofalgae<- as.factor(clg$Freeofalgae) clg$Sharedwithanimals<- as.factor(clg$Sharedwithanimals) clg$Sharedwithanimals[clg$Sharedwithanimals==""]<-"No" clg$Hasanimaldroppingsorfaeces<- as.factor(clg$Hasanimaldroppingsorfaeces) clg$Foodisstoredonthefloor<- as.factor(clg$Foodisstoredonthefloor) clg$EvidenceofinfestationwithfliesFPA<- as.factor(clg$EvidenceofinfestationwithfliesFPA) clg$Closetothepitlatrineortoilet[clg$Closetothepitlatrineortoilet==""]<-"No" clg$Closetothepitlatrineortoilet<- as.factor(clg$Closetothepitlatrineortoilet) clg$PitlatrineUnderuse[clg$PitlatrineUnderuse==""]<-"No" clg$PitlatrineUnderuse<- as.factor(clg$PitlatrineUnderuse) clg$Contaminationintheenvironment[is.na(clg$Contaminationintheenvironment)]<-"No" clg$Contaminationintheenvironment<- as.factor(clg$Contaminationintheenvironment) clg<- droplevels.data.frame(clg) clg$DrinkingWaterTx <- rep(NA, nrow(clg)) clg$DrinkingWaterTx[clg$Treatdrinkingwaterbeforestorage=="No"] <- "No" clg$DrinkingWaterTx[clg$Boiling=="Yes"] <- "Boiling" clg$DrinkingWaterTx[clg$Chlorination=="Yes"] <- "Chlorination" clg$DrinkingWaterTx[clg$Solardisinfection=="Yes"] <- "Solardisinfection" clg$DrinkingWaterTx[clg$Filtration=="Yes"] <- "Filtration" ## new columns to add clg$Accesstosafedrinkingwater <- as.character(clg$Accesstosafedrinkingwater) which.one <- which( levels(clg$Accesstotheuntreateddrinkingwater) == "" ) levels(clg$Accesstotheuntreateddrinkingwater)[which.one] <- "Missing" clg$Accesstotheuntreateddrinkingwater <- factor(clg$Accesstotheuntreateddrinkingwater, levels = c("No", "Yes", "Missing")) clg$Fromwheredoyouwashyourhand. <- as.character(clg$Fromwheredoyouwashyourhand.) ####************************************************ #### COMMENTS ON SOME COLUMS THAT CAN BE RECORDED OR REMOVED ####************************************************ ### Treatdrinkingwaterbeforestorage = YES & NO ,if yes this is how they treat colunm (Boiling Chlorination Solardisinfection Filtration) ### Can modify this how do treatwater response, Don' treat == NO from before, Boiling, Chlorinating, Solar and Filtration as response for Yes ### Lackoffirewoodtoboil,Lackofmoneytobuychlorinetablets,Itisnotinherenttotreatdrinkingwater,Believethewaterissafe,Havenotime,Feelsickwhentheytaketreatedwater,Tastelessafterboiling,Unavailabilityofchlorinetablets,Fearforchildrensickness ## these are follow up responses for those who do not treat water.. asking why they do not treat..(Consider droping them or recode including response "treat water") ## Accesstosafedrinkingwater & Accesstotheuntreateddrinkingwater. Redundant question... if someone responded yes to Tap water ## WHICH ANIMALS TO HAVE PET=Dog,Guineapigs &Cat, LIVESTOCK=(Pigs,Goats,Cows), POULTRY=(Chicken ,Ducks) # remove useless variables! clg <- select(clg, -DrinkingWaterTx,-Dateofinterview,-Questionnairenumber, - UniqueKey, -Nameofparticipant,-Bonemarrow, -Subcounty, -Parish, -Village, -Parish.latitude, -Victimshouldnotpreparefooduntildeclaredfreeoftyphoidfever, -Parish.longitude, -Examplesofreadytoeatfoods, -Householdmembers, -Positioninthehousehold, -typeofwatersource1, -Typeofwatersource2, -Notes, -Watersourceinspected,-Watercollectionvessel, -Speciesoffaecalsample1, -Fromwheredoyouwashyourhand., -Whydontyouseekmedicalattention, -Designatedcupkept,-Wheredoyoustoreprepredorleftoverfood, -FS1collected, -SalmonellaFS1, -ShigellaFS1, -EColiFS1, -FS2collected, -Speciesoffaecalsample2, -SalmonellaFS2, -ShigellaFS2, -EColiFS2, -Contaminated.fecal.sample, -PresenceofseparatevesselforDW, - KindofseparatevesseiforDW, -Presenceofdesignatedcup, # -Accesstosafedrinkingwater, # -Accesstotheuntreateddrinkingwater, -Fromwheredoyouwashyourhand., #-Wateronlytowashhands, #-Ashandwater, #-Rubagainstthesurface, #-Usemaizecobs, -Usenaturalanalcleansingleaves, -Usewaterforanal, -Usewaterandsaopforanal, -Managementofchildrensfecalwaste, -Howsoondoyoueatthefoodafterpreparation, -Doyoupreheatleftoverfood, -Wherefoodisprepared, -Howsoondoyoueatthefoodafterpreparation, -Wheredoyoustoreprepredorleftoverfood,-Doyouusuallyhaveleftovers, -Howdoyoudealwithleftovers,-Utensilstodryafterwashing, -Utensilsinthekitchen, -Arevegetablesusuallyapartofyourdiet, -Anyfamilymemberorsomeoneclosewithtyphoid, -Noneoftheabove, -Fever, -Diarrhoea1, -Headache, -Shivering, -Nausea, -Constipation1, -Jointormusclepain, -Defaultthemedicationgiven, -Anypractisestoavoidtransmission, -Didyoufollowtheadvicegiven, -Watermoving, -Watertreatedatsource, -Cleanontheoutside, -Hasacover, -Freeofholesandleakages, -Usedforstorageofdrinkingwater,-Drinkingwaterstoragevessel, -CleanontheoutsideDW, -FreeofalgaeDW, -HasacoverDW, -FreeofholesandleakagesDW, -Designatedcuptodrawthedrinkingwaterseen, -Hasadoororbarrier,-Utensilsarekeptonthefloor, -Leftoverfoodpresent, -HandwashingfacilitypresentFPA, -HWEvidenceofuseFPA, -HWFhasassoaporashFPA, -Humanfaecalmatterclosetothearea, -PitlatrinePresent, -Evidenceofproperuseofpitlatrine, -Holeiscovered, -Evidentflyinfestationatpitlatrine, -Magnitudeoftheinfestationatpitlatrine, -Handwashingfacilitypresentatpitlatrine,-HWEvidenceofuseattoilet, -Drinkingwatersample,-SalmonellaWS1, -EColiWS1, -ShigellaWS1,-Watersourcesamplecollected2, -Typeofwatersource2, -SalmonellaDW1, -EColiDW1, -ShigellaDW1, -Seconddrinkingwatersamplecollected, -Notes1, -SalmonellaDW2, -EColiDW2, -ShigellaDW2, -OverallContaminationofDWwithSalmonella, -Watersourcesamplecollected, -typeofwatersource1, -SalmonellaWS2, -EColiWS2, -ShigellaWS2, -Contamination.of.water.at.the.source, -FS1collected, -Speciesoffaecalsample1, -SalmonellaFS1, -ShigellaFS1, -EColiFS1, -FS2collected, -Speciesoffaecalsample2, -SalmonellaFS2, -ShigellaFS2, -EColiFS2, -Contaminated.fecal.sample, -Householdmembers_cat, -Lackoffirewoodtoboil,-Lackofmoneytobuychlorinetablets,-Itisnotinherenttotreatdrinkingwater,-Believethewaterissafe,-Havenotime,-Feelsickwhentheytaketreatedwater,-Tastelessafterboiling,-Unavailabilityofchlorinetablets,-Fearforchildrensickness, -Timesofconfirmedtyphoidfever, - Swollenpainfulabdomen, -Urine, -Stool, -Blood, -TyphoidfeverconfirmedontheveryfirstvisittoHC, -Noadvicegiven, -Vomiting, -Saliva, -Eatfoodwhilehot, -Painfulabdomennotswollen,-Dontshareutensilswiththevictim,#-Evenonurination, -Feverishconditions,-Diarrhoea, -Constipation, -Malaise, -Abdominalpains, -Howoftendoyouexperiencethesesymptoms, -Seekmedicalattentionwhenwithsymptoms, -Medicalresponseafteronsetofsymptoms) clg<- droplevels.data.frame(clg) # univariable analysis library(survival) uniques<-data.frame(name=vector("character",0), levels=vector("character",0)) for( i in 1:ncol(clg)){ a<-length(unique(clg[,i])) b <- names(clg)[i] new <- data.frame(name=b, levels=a) uniques <- rbind(uniques, new) } uniques1 <- arrange(uniques, levels) %>% filter(levels==1) clg <- select(clg, -which(names(clg) %in% unique(uniques1$name))) # example model options(na.action = "na.omit") m1 <-clogit(CaseorControl ~ Levelofeducation + strata(match), data=clg) test<-clg %>% select(-CaseorControl, -match) %>% map(~clogit(clg$CaseorControl ~ .x + strata(clg$match),method="approximate", data=clg)) result <- data.frame(Variable=vector("character",0), OR=vector("character",0), LL=vector("character",0), UL=vector("character",0), pv=vector("character",0), filter=vector("character",0)) for(i in 1:length(test)){ a<-mmt(test[[i]]) a$Variable[1] <- names(test[i]) result <- rbind(result, a) } ## names of variables to consider in the model arranged smaller p-value first for forward variable selection formodel <- arrange(filter(result, filter=="yes" & pv<0.15),pv)$Variable # dataframe to use for the model df_model <- select(clg, formodel, CaseorControl, match) #df_model <- select(df_model,-KindofseparatevesseiforDW) # find columns with NAs #which(apply(df_model, 2, anyNA)) df_model <- df_model %>% drop_na() ### Forward model selection ## consider each variable in the order it appears in "formodel" ## eg the first one is formodel[1] # AIC 254.2861 m1 <-clogit(CaseorControl ~ Treatdrinkingwaterpractice + strata(match), data=df_model, method="approximate") #AIC(m1) # AIC 255.9174 m2 <-clogit(CaseorControl ~ Treatdrinkingwaterpractice + Boiling + strata(match), data=df_model, method="approximate") #AIC(m2) ## FINAL MODEL # AIC 253.517 <- final model m3 <-clogit(CaseorControl ~ Treatdrinkingwaterpractice + Soapandwatertowashhands + strata(match), data=df_model, method="approximate") #AIC(m3) ## final table select(mmt(m3), -filter) # diagnostics library(ROCR) prob<-predict(m3,type = "expected") pred <- prediction(prob, df_model$CaseorControl) perf <- performance(pred, measure = "tpr", x.measure = "fpr") plot(perf) auc <- performance(pred, measure = "auc") auc <- auc@y.values[[1]] auc # AIC 256.8196 m4 <-clogit(CaseorControl ~ Treatdrinkingwaterpractice + Soapandwatertowashhands + Levelofeducation + strata(match), data=df_model, method="approximate") #AIC(m4) # AIC 255.4725 m5 <-clogit(CaseorControl ~ Treatdrinkingwaterpractice + Soapandwatertowashhands + N3644years + strata(match), data=df_model, method="approximate") #AIC(m5) # AIC 255.0848 m6 <-clogit(CaseorControl ~ Treatdrinkingwaterpractice + Soapandwatertowashhands + Washhandsaftertoiletorlatrine + strata(match), data=df_model, method="approximate") #AIC(m6) # AIC 254.4313 m7 <-clogit(CaseorControl ~ Treatdrinkingwaterpractice + Soapandwatertowashhands + Above65years + strata(match), data=df_model, method="approximate") #AIC(m7) # AIC 254.5589 m8 <-clogit(CaseorControl ~ Treatdrinkingwaterpractice + Soapandwatertowashhands + Usepapers + strata(match), data=df_model, method="approximate") #AIC(m8) ### backward elimination # AIC 264.2622 mall <-clogit(CaseorControl ~ Treatdrinkingwaterpractice + Boiling + Soapandwatertowashhands + Levelofeducation + N3644years + Washhandsaftertoiletorlatrine + Above65years + Usepapers + strata(match), data=df_model, method="approximate") AIC(mall) # AIC 263.6646 m9 <-clogit(CaseorControl ~ Treatdrinkingwaterpractice + Boiling + Soapandwatertowashhands + Levelofeducation + N3644years + Washhandsaftertoiletorlatrine + Above65years + strata(match), data=df_model, method="approximate") AIC(m9) # AIC 262.2305 m10 <-clogit(CaseorControl ~ Treatdrinkingwaterpractice + Boiling + Soapandwatertowashhands + Levelofeducation + N3644years + Washhandsaftertoiletorlatrine + strata(match), data=df_model, method="approximate") AIC(m10) # AIC 260.477 m11 <-clogit(CaseorControl ~ Treatdrinkingwaterpractice + Boiling + Soapandwatertowashhands + Levelofeducation + N3644years + strata(match), data=df_model, method="approximate") AIC(m11) # AIC 258.4978 m12 <-clogit(CaseorControl ~ Treatdrinkingwaterpractice + Boiling + Soapandwatertowashhands + Levelofeducation + strata(match), data=df_model, method="approximate") AIC(m12) # AIC 255.273 m13 <-clogit(CaseorControl ~ Treatdrinkingwaterpractice + Boiling + Soapandwatertowashhands + strata(match), data=df_model, method="approximate") AIC(m13) # AIC 255.9174 m14 <-clogit(CaseorControl ~ Treatdrinkingwaterpractice + Boiling + strata(match), data=df_model, method="approximate") AIC(m14) # AIC 253.517 <- Final model m15 <-clogit(CaseorControl ~ Treatdrinkingwaterpractice + Soapandwatertowashhands + strata(match), data=df_model, method="approximate") AIC(m15) # AIC 312.4566 m16 <-clogit(CaseorControl ~ Soapandwatertowashhands + strata(match), data=df_model, method="approximate") AIC(m16) ## this is me cheating and ignoring Soapandwatertowashhands variable to see what happens. Model still seems upside down # AIC 313.4555 m1 <-clogit(CaseorControl ~ Boiling + strata(match), data=df_model, method="approximate") #AIC(m1) # AIC 311.5513 m2 <-clogit(CaseorControl ~ Boiling + Soapandwatertowashhands + strata(match), data=df_model, method="approximate") #AIC(m2) # AIC 309.5377 m3 <-clogit(CaseorControl ~ Boiling + Soapandwatertowashhands + Levelofeducation + strata(match), data=df_model, method="approximate") #AIC(m3) # AIC 309.2024 m4 <-clogit(CaseorControl ~ Boiling + Soapandwatertowashhands + Levelofeducation + N3644years + strata(match), data=df_model, method="approximate") #AIC(m4) # AIC 311.5176 m5 <-clogit(CaseorControl ~ Boiling + Soapandwatertowashhands + Levelofeducation + Washhandsaftertoiletorlatrine + strata(match), data=df_model, method="approximate") #AIC(m5) # AIC 309.4235 m6 <-clogit(CaseorControl ~ Boiling + Soapandwatertowashhands + Levelofeducation + Above65years + strata(match), data=df_model, method="approximate") #AIC(m6) # AIC 308.9699 <- final model if cheat m7 <-clogit(CaseorControl ~ Boiling + Soapandwatertowashhands + Levelofeducation + Usepapers + strata(match), data=df_model, method="approximate") #AIC(m7) # AIC 254.5589 m8 <-clogit(CaseorControl ~ Treatdrinkingwaterpractice + Soapandwatertowashhands + Usepapers + strata(match), data=df_model, method="approximate") #AIC(m8) ## model selection library(MuMIn) mall <-clogit(CaseorControl ~ . + strata(match), data=df_model) options(na.action = "na.fail") all_models<-dredge(mall,evaluate = F) all_models_eval<-dredge(mall,evaluate = T) #save(all_models_eval, file = "all_models_eval.RData") #all_models_get <-get.models(all_models_eval, subset = T) ma <- model.avg(get.models(all_models_eval, subset = delta <= 2)) summary(ma) save(ma, file = "ma.RData") confint(ma, full = T) confint(ma, full = F) ``` ```{r} ### Reviewer's comments ### Forward model selection ## consider each variable in the order it appears in "formodel" ## eg the first one is formodel[1] # AIC 254.2861 m1 <-clogit(CaseorControl ~ Treatdrinkingwaterpractice + strata(match), data=df_model, method="approximate") #AIC(m1) # AIC 255.9174 - remove m2 <-clogit(CaseorControl ~ Treatdrinkingwaterpractice + Boiling + strata(match), data=df_model, method="approximate") #AIC(m2) # AIC 253.517 - keep m3 <-clogit(CaseorControl ~ Treatdrinkingwaterpractice + Soapandwatertowashhands + strata(match), data=df_model, method="approximate") #AIC(m3) # AIC 256.8196 - so remove m4 <-clogit(CaseorControl ~ Treatdrinkingwaterpractice + Soapandwatertowashhands + Levelofeducation + strata(match), data=df_model, method="approximate") #AIC(m4) # AIC 255.5073 - so remove m5 <-clogit(CaseorControl ~ Treatdrinkingwaterpractice + Soapandwatertowashhands + Evenonurination + strata(match), data=df_model, method="approximate") #AIC(m5) # AIC 255.4725 - so remove m6 <-clogit(CaseorControl ~ Treatdrinkingwaterpractice + Soapandwatertowashhands + N3644years + strata(match), data=df_model, method="approximate") #AIC(m6) # AIC 255.2995 - so remove m7 <-clogit(CaseorControl ~ Treatdrinkingwaterpractice + Soapandwatertowashhands + Goats + strata(match), data=df_model, method="approximate") #AIC(m7) # AIC 255.0848 - so remove m8 <-clogit(CaseorControl ~ Treatdrinkingwaterpractice + Soapandwatertowashhands + Washhandsaftertoiletorlatrine + strata(match), data=df_model, method="approximate") #AIC(m8) # AIC 254.4313 - so remove m9 <-clogit(CaseorControl ~ Treatdrinkingwaterpractice + Soapandwatertowashhands + Above65years + strata(match), data=df_model, method="approximate") #AIC(m9) # AIC 254.5589 - so remove m10 <-clogit(CaseorControl ~ Treatdrinkingwaterpractice + Soapandwatertowashhands + Usepapers + strata(match), data=df_model, method="approximate") #AIC(m10) ### backward elimination # AIC 264.2622 mall <-clogit(CaseorControl ~ Treatdrinkingwaterpractice + Boiling + Soapandwatertowashhands + Levelofeducation + Evenonurination + N3644years + Goats + Washhandsaftertoiletorlatrine + Above65years + Usepapers + strata(match), data=df_model, method="approximate") AIC(mall) # AIC 310.0997 #add Treatdrinkingwaterpractice at the end m11 <-clogit(CaseorControl ~ Boiling + Soapandwatertowashhands + Levelofeducation + Evenonurination + N3644years + Goats + Washhandsaftertoiletorlatrine + Above65years + Usepapers + strata(match), data=df_model, method="approximate") AIC(m11) # AIC 310.0997 #removed Usepapers, AIC didn't go down so kept m12 <-clogit(CaseorControl ~ Boiling + Soapandwatertowashhands + Levelofeducation + Evenonurination + N3644years + Goats + Washhandsaftertoiletorlatrine + Above65years + strata(match), data=df_model, method="approximate") AIC(m12) # AIC 310.7291 #removed Above65years, AIC didn't go down so kept m13 <-clogit(CaseorControl ~ Boiling + Soapandwatertowashhands + Levelofeducation + Evenonurination + N3644years + Goats + Washhandsaftertoiletorlatrine + Usepapers + strata(match), data=df_model, method="approximate") AIC(m13) # AIC 308.4638 #removed Washhandsaftertoiletorlatrine, AIC went down so removed m14 <-clogit(CaseorControl ~ Boiling + Soapandwatertowashhands + Levelofeducation + Evenonurination + N3644years + Goats + Above65years + Usepapers + strata(match), data=df_model, method="approximate") AIC(m14) # AIC 308.9402 #removed Goats, AIC didn't go down so kept m15 <-clogit(CaseorControl ~ Boiling + Soapandwatertowashhands + Levelofeducation + Evenonurination + N3644years + Above65years + Usepapers + strata(match), data=df_model, method="approximate") AIC(m15) # AIC 308.9678 #removed Usepapers, AIC didn't go down so kept m16 <-clogit(CaseorControl ~ Boiling + Soapandwatertowashhands + Levelofeducation + Evenonurination + N3644years + Above65years + Goats + strata(match), data=df_model, method="approximate") AIC(m16) # AIC 308.9678 #removed Usepapers, AIC didn't go down so kept m16 <-clogit(CaseorControl ~ Boiling + Soapandwatertowashhands + Levelofeducation + Evenonurination + N3644years + Above65years + Goats + Usepapers + strata(match), data=df_model, method="approximate") AIC(m16) ``` \newpage ```{r, echo=TRUE} ######### # RC Analysis were missing values are included in the dataset and the full dataset is used ########## #Get weather data for each week between 2012-2018 weather <- read.csv("kasese_weather_20180308.csv") weather$TotalRain <- rowSums(weather[,c( "prec0", "prec3", "prec6", "prec9", "prec12", "prec15","prec18", "prec21")]) weather$MeanRain <- rowMeans(weather[,c( "prec0", "prec3", "prec6", "prec9", "prec12", "prec15","prec18", "prec21")], na.rm = TRUE) weather$meantempC <- rowMeans(weather[,c("temp0", "temp3", "temp6", "temp9", "temp12", "temp15", "temp18", "temp21")], na.rm = TRUE) weather$meanhumidity <- rowMeans(weather[,c("hum0", "hum3", "hum6", "hum9", "hum12", "hum15","hum18", "hum21")], na.rm = TRUE) weather$Week <- as.integer(week(ymd(weather$date))) weather$Year <- year(ymd(weather$date)) weather$Month <- as.character((month(ymd(weather$date),label=TRUE,abbr=FALSE))) weather$Month <- factor(weather$Month, levels=c("January","February","March","April","May","June","July","August","September","October","November","December")) weather$date <- as.Date(weather$date, "%Y-%m-%d") weekdate <- aggregate(weather$date, list(weather$Week, weather$Year), min) colnames(weekdate) <- c("Week", "Year", "StartDate") weather <- merge(weather, weekdate, by.weather = c("Week", "Year"), by.weekdate = c("Week", "Year"), all.x=TRUE) weather1 <- weather[,c("date", "Year", "Week", "Month", "maxtempC", "mintempC", "meantempC", "TotalRain", "MeanRain", "meanhumidity", "StartDate")] weekly_weather <- weather1 %>% group_by(Year, Week, Month, StartDate) %>% summarise_at(vars(-date), funs(mean(., na.rm=TRUE))) #Getting the missing weeks for 2013-2017 weekly_weather1317 <- weekly_weather[weekly_weather$Year >= 2013 & weekly_weather$Year <= 2017,] FULL_MOH_na <- merge(weekly_weather1317, FULL_MOH, by = c("Year", "Week", "Month"), all.x = TRUE) #Removing extra weeks in 2017 for which we don't have case data FULL_MOH_na <- subset(FULL_MOH_na, !(FULL_MOH_na$Year == 2017 & FULL_MOH_na$Week > 23)) dim(FULL_MOH_na) #Everthing up until 2016 FULL_MOH_na_2016 <- subset(FULL_MOH_na, FULL_MOH_na$Year < 2017) #2017 FULL_MOH_na_2017 <- subset(FULL_MOH_na, FULL_MOH_na$Year == 2017) #Adding 0 for cases when there was no outbreaks FULL_MOH_na_2016$Cases.in.Kasese[which(FULL_MOH_na_2016$Week >=1 & FULL_MOH_na_2016$Week <= 13 & FULL_MOH_na_2016$Year==2013)] <- 0 FULL_MOH_na_2016$Cases.in.Uganda[which(FULL_MOH_na_2016$Week >=1 & FULL_MOH_na_2016$Week <= 13 & FULL_MOH_na_2016$Year==2013)] <- 0 FULL_MOH_na_2016$Cases.in.Kasese[which(FULL_MOH_na_2016$Week >=39 & FULL_MOH_na_2016$Year==2014)] <- 0 FULL_MOH_na_2016$Cases.in.Kasese[which(FULL_MOH_na_2016$Week <=12 & FULL_MOH_na_2016$Year==2015)] <- 0 FULL_MOH_na_2016$Cases.in.Uganda[which(FULL_MOH_na_2016$Week >=1 & FULL_MOH_na_2016$Week <= 13 & FULL_MOH_na_2016$Year==2013)] <- 0 FULL_MOH_na_2016$Cases.in.Uganda[which(FULL_MOH_na_2016$Week >=39 & FULL_MOH_na_2016$Year==2014)] <- 0 FULL_MOH_na_2016$Cases.in.Uganda[which(FULL_MOH_na_2016$Week <=12 & FULL_MOH_na_2016$Year==2015)] <- 0 FULL_MOH_na_2016_rain <- melt(FULL_MOH_na_2016[,c("Year", "Week", "Month","StartDate","TotalRain", "MeanRain")], id=c("Year", "Week", "Month", "StartDate")) FULL_MOH_na_2016_temp <- melt(FULL_MOH_na_2016[,c("Year", "Week", "Month", "StartDate","maxtempC", "mintempC", "meantempC")], id=c("StartDate","Year", "Week", "Month")) FULL_MOH_na_2016_hum <- melt(FULL_MOH_na_2016[,c("StartDate","Year", "Week", "Month", "meanhumidity")], id=c("StartDate","Year", "Week", "Month")) #Time series using TS FULL_MOH_kasese_na_ts <- ts(FULL_MOH_na_2016$Cases.in.Kasese, start=c(2013,1), end=c(2017,1), frequency=52 ,ts.eps = getOption("ts.eps"), is.ts(x)) ts_plot_kasese <- plot(FULL_MOH_kasese_na_ts) FULL_MOH_uganda_na_ts <- ts(FULL_MOH_na_2016$Cases.in.Uganda, start=c(2013,1), end=c(2017,1), frequency=52 ,ts.eps = getOption("ts.eps"), is.ts(x)) ts_plot_uganda <- plot(FULL_MOH_uganda_na_ts) FULL_MOH_temp_na_ts <- ts(FULL_MOH_na_2016$meantempC, start=c(2013,1), end=c(2017,1), frequency=52 ,ts.eps = getOption("ts.eps"), is.ts(x)) ts_plot_temp <- plot(FULL_MOH_temp_na_ts) FULL_MOH_hum_na_ts <- ts(FULL_MOH_na_2016$meanhumidity, start=c(2013,1), end=c(2017,1), frequency=52 ,ts.eps = getOption("ts.eps"), is.ts(x)) ts_plot_hum <- plot(FULL_MOH_hum_na_ts) FULL_MOH_rain_na_ts <- ts(FULL_MOH_na_2016$TotalRain, start=c(2013,1), end=c(2017,12), frequency=52 ,ts.eps = getOption("ts.eps"), is.ts(x)) ts_plot_rain <- plot(FULL_MOH_rain_na_ts) FULL_MOH_2017_ts <- ts(FULL_MOH_na_2017$Cases.in.Kasese, start=c(2017,1), end=c(2017,12), frequency=52 ,ts.eps = getOption("ts.eps"), is.ts(x)) FULL_MOH_2017_UG_ts <- ts(FULL_MOH_na_2017$Cases.in.Uganda, start=c(2017,1), end=c(2017,12), frequency=52 ,ts.eps = getOption("ts.eps"), is.ts(x)) #=================================================== # PLOT FOR IMPUTE MISSING DATA FOR TIME SERIES time series #=================================================== imp_kalman_kasese_na <- na.kalman(FULL_MOH_kasese_na_ts) ## Kasese imp_ma_kasese_na <- na.ma(FULL_MOH_kasese_na_ts) imp_mean_kasese_na <- na.mean(FULL_MOH_kasese_na_ts) imp_random_kasese_na <- na.random(FULL_MOH_kasese_na_ts) imp_locf_kasese_na <- na.locf(FULL_MOH_kasese_na_ts) imp_interpolation_kasese_na <- na.interpolation(FULL_MOH_kasese_na_ts) imp_kalman_Uganda_na <- na.kalman(FULL_MOH_uganda_na_ts)## Uganda imp_ma_Uganda_na <- na.ma(FULL_MOH_uganda_na_ts) imp_mean__Uganda_na <- na.mean(FULL_MOH_uganda_na_ts) imp_random_Uganda_na <- na.random(FULL_MOH_uganda_na_ts) imp_locf_Uganda_na <- na.locf(FULL_MOH_uganda_na_ts) imp_interpolation_Uganda_na <- na.interpolation(FULL_MOH_uganda_na_ts) imp_kalman_kasese_tx <- na.kalman(FULL_MOH_na_2016$Cases.in.Kasese) imp_ma_kasese_na_tx <- na.ma(FULL_MOH_na_2016$Cases.in.Kasese) imp_mean_kasese_na_tx <- na.mean(FULL_MOH_na_2016$Cases.in.Kasese) imp_random_kasese_na_tx <- na.random(FULL_MOH_na_2016$Cases.in.Kasese) imp_locf_kasese_na_tx <- na.locf(FULL_MOH_na_2016$Cases.in.Kasese) imp_interpolation_kasese_na_tx <- na.interpolation(FULL_MOH_na_2016$Cases.in.Kasese) imp_dataset_kasese <- data.frame(StartDate=FULL_MOH_na_2016$StartDate, OriginalTS = FULL_MOH_na_2016$Cases.in.Kasese, Kalman = imp_kalman_kasese_tx, MA = imp_ma_kasese_na_tx, # Mean = imp_mean_kasese_na_tx, # Random = imp_random_kasese_na_tx, LOCF = imp_locf_kasese_na_tx, Interpolation = imp_interpolation_kasese_na_tx) imp_dataset_kasese1 <- melt(imp_dataset_kasese, id="StartDate") imp_dataset_kasese1$Week <- as.integer(week(ymd(imp_dataset_kasese1$StartDate))) imp_dataset_kasese1$Year <- year(ymd(imp_dataset_kasese1$StartDate)) # plotNA.imputations(x.withNA = FULL_MOH_na_2016$Cases.in.Kasese, x.withImputations = na.seadec(FULL_MOH_na_2016$Cases.in.Kasese, "kalman")) # plotNA.imputations(x.withNA = FULL_MOH_na_2016$Cases.in.Kasese, x.withImputations = na.seadec(FULL_MOH_na_2016$Cases.in.Kasese,"ma")) imp_kalman_uganda_na_tx <- na.kalman(FULL_MOH_na_2016$Cases.in.Uganda) imp_ma_uganda_na_tx <- na.ma(FULL_MOH_na_2016$Cases.in.Uganda) imp_mean_uganda_na_tx <- na.mean(FULL_MOH_na_2016$Cases.in.Uganda) imp_random_uganda_na_tx <- na.random(FULL_MOH_na_2016$Cases.in.Uganda) imp_locf_uganda_na_tx <- na.locf(FULL_MOH_na_2016$Cases.in.Uganda) imp_interpolation_uganda_na_tx <- na.interpolation(FULL_MOH_na_2016$Cases.in.Uganda) imp_dataset_uganda <- data.frame(StartDate=FULL_MOH_na_2016$StartDate, OriginalTS = FULL_MOH_na_2016$Cases.in.Uganda, Kalman = imp_kalman_uganda_na_tx, MA = imp_ma_uganda_na_tx, # Mean = imp_mean_uganda_na_tx, # Random = imp_random_uganda_na_tx, LOCF = imp_locf_uganda_na_tx, Interpolation = imp_interpolation_uganda_na_tx) imp_dataset_uganda1 <- melt(imp_dataset_uganda, id="StartDate") imp_dataset_uganda1$Week <- as.integer(week(ymd(imp_dataset_uganda1$StartDate))) imp_dataset_uganda1$Year <- year(ymd(imp_dataset_uganda1$StartDate)) # # # Kasese_2017_ts <- ts(FULL_MOH_na_2017$Cases.in.Kasese, start=c(2017, 1), end=c(2017, 12), frequency=52,ts.eps = getOption("ts.eps"), is.ts(x)) # # plot(Kasese_2017_ts) # # # #Plotting the data imput_kasese_plot <- ggplot(imp_dataset_kasese1, aes(x=Week, y=value, col=variable)) + geom_line(stat = "identity") + facet_grid(.~Year) + scale_y_continuous("Number of cases (Kasese)") + scale_color_discrete(name="", breaks=c("OriginalTS", "Kalman", "MA", "LOCF", "Interpolation"), labels=c("Original time series", "Kalman smoothing", "Weighted moving average", "Last observation carried forward", "Interpolation")) + scale_x_continuous("") + annotate("text", x = c(16,19), y = 28, label = "*") ### Rebecca please add these annotations for flooding reported for (week 16 in 2013, 19 in 2014, 17 & 27 in 2015, 15 in 2016) imput_uganda_plot <- ggplot(imp_dataset_uganda1, aes(x=Week, y=value, col=variable)) + geom_line(stat = "identity") + facet_grid(.~Year) + scale_y_continuous("Number of cases (Uganda)") + scale_color_discrete(name="", breaks=c("OriginalTS", "Kalman", "MA", "LOCF", "Interpolation"), labels=c("Original time series", "Kalman smoothing", "Weighted moving average", "Last observation carried forward", "Interpolation")) + scale_x_continuous("") rain_plot <- ggplot(FULL_MOH_na_2016_rain, aes(x = Week, y = value, col=variable)) + geom_line(stat = "identity") + scale_y_continuous("Rainfall (mm)") + facet_grid(.~Year) + scale_color_discrete(name="", breaks=c("TotalRain", "MeanRain"), labels=c("Total rainfall", "Mean rainfall")) + scale_x_continuous("") temp_plot <- ggplot(FULL_MOH_na_2016_temp, aes(x = Week, y = value, col=variable)) + geom_line(stat = "identity") + scale_y_continuous(expression(paste("Temperature ",(degree*C)))) + facet_grid(.~Year) + scale_color_discrete(name="", breaks=c("maxtempC", "mintempC", "meantempC"), labels=c("Maximum temperature", "Minimum temperature", "Mean temperature")) + scale_x_continuous("") humidity_plot <- ggplot(FULL_MOH_na_2016_hum, aes(x = Week, y = value, col=variable)) + geom_line(stat = "identity") + scale_y_continuous("Humidity (%)") + facet_grid(.~Year) + scale_color_discrete(name="", breaks=c("meanhumidity"), labels=c("Mean humidity")) # # #grid.arrange(imput_kasese_plot, imput_uganda_plot, rain_plot, temp_plot, humidity_plot, nrow = 5, common.legend = TRUE, legend="bottom")---- pdf("../../Figures/Imputed_data.pdf",width = 10,height = 14) multiplot(imput_kasese_plot, imput_uganda_plot,rain_plot,temp_plot,humidity_plot, cols = 1) dev.off() ts_meanrain <- ts(FULL_MOH_na_2016$MeanRain) ts_totalrain <- ts(FULL_MOH_na_2016$TotalRain) ts_maxtemp <- ts(FULL_MOH_na_2016$maxtempC) ts_mintemp <- ts(FULL_MOH_na_2016$mintempC) ts_meantemp <- ts(FULL_MOH_na_2016$meantempC) ts_meanhum <- ts(FULL_MOH_na_2016$meanhumidity) ####******************************************* #### TS CROSS CORRELATION WITH WEATHER IN KASESE ####******************************************* Typhoid_cases<-imp_kalman_kasese_na Precipitation_in_Kasese<-FULL_MOH_rain_na_ts Temperature_in_Kasese<-FULL_MOH_temp_na_ts Humidity_in_Kasese<-FULL_MOH_hum_na_ts ccf(Typhoid_cases,Precipitation_in_Kasese,type = "correlation",plot = TRUE) ccf(Typhoid_cases,Precipitation_in_Kasese,type = "correlation",plot = TRUE) ccf(Typhoid_cases,Temperature_in_Kasese,type = "correlation",plot = TRUE) #ccf(imp_kalman_kasese_na,ts_mintemp,type = "correlation",plot = TRUE) #ccf(imp_kalman_kasese_na,ts_meantemp,type = "correlation",plot = TRUE) ccf(Typhoid_cases,Humidity_in_Kasese,type = "correlation",plot = TRUE) ####******************************************* #### SEASONALITY ANALYSIS UG AND KASESE ####******************************************* myts_NAT12 <- ts(FULL_MOH_na_2016$Cases.in.Uganda, start=c(2013, 1), end=c(2017, 1), frequency=12,ts.eps = getOption("ts.eps"), is.ts(x)) ### NATIONAL AT myts_DIS12 <- ts(FULL_MOH_na_2016$Cases.in.Kasese, start=c(2013, 1), end=c(2017, 1), frequency=12,ts.eps = getOption("ts.eps"), is.ts(x)) ### DISTRICT AT MOTHLY FREQUENCY myts_RAIN12 <- ts(FULL_MOH_na_2016$TotalRain, start=c(2013, 1), end=c(2017, 1), frequency=12,ts.eps = getOption("ts.eps"), is.ts(x)) ### RAIN AT DISTRICT MONTHLY myts_TEMP12 <- ts(FULL_MOH_na_2016$maxtempC, start=c(2013, 1), end=c(2017, 1), frequency=12,ts.eps = getOption("ts.eps"), is.ts(x)) ### TEMPERATURE AT DISTRICT MONTHLY myts_HUM12<- ts(FULL_MOH_na_2016$meanhumidity, start=c(2013, 1), end=c(2017, 1), frequency=12,ts.eps = getOption("ts.eps"), is.ts(x)) ### HUMIDITY AT DISTRICT MONTHLY # fit_NAT_1 <- ets(myts_NAT12) # res <- cbind(Residuals = residuals(fit_NAT_1), # Response.residuals = residuals(fit_NAT_1, type='response')) # autoplot(res, facets=TRUE) # # fit_DIS_1 <- ets(myts_DIS12) # res <- cbind(Residuals = residuals(fit_DIS_1), # Response.residuals = residuals(fit_DIS_1, type='response')) # autoplot(res, facets=T) imp_kalman_myts_NAT12 <- na.kalman(myts_NAT12) ## National impute season imp_kalman_myts_DIS12 <- na.kalman(myts_DIS12) ## Kasese impute season R5<-ggseasonplot(imp_kalman_myts_NAT12, polar=TRUE) + theme_bw() + labs(title= "Typhoid fever in Uganda ") R1<-ggseasonplot(imp_kalman_myts_DIS12, polar=TRUE) + theme_bw() + labs(title= "Typhoid fever in Kasese District") R2<-ggseasonplot(myts_RAIN12, polar=TRUE) + theme_bw() + labs(title= "Rainfall in Kasese District") R3<-ggseasonplot(myts_TEMP12, polar=TRUE) + theme_bw() + labs(title= "Temperature in Kasese District") R4<-ggseasonplot(myts_HUM12, polar=TRUE) + theme_bw() + labs(title= "Humidity in Kasese District") multiplot(R1,R2,R3,R4, cols = 2) pdf("../../Figures/Seasonality_Kasese.pdf",width = 8,height = 6) multiplot(R1,R2,R3,R4, cols = 2) dev.off() ####********************** #### TS MODEL FITTING ####********************** fit_kas<-auto.arima(imp_kalman_kasese_na) ### Kasese typhoid data TS fitting using kalman imputed data checkresiduals(fit_kas) fit_UG<-auto.arima(imp_kalman_Uganda_na) ### Kasese typhoid data TS fitting using kalman imputed data checkresiduals(fit_UG) ####********************** #### FORCASTING#### ####********************** Kas_for_17_intp <- forecast(object = imp_interpolation_kasese_na, h = 16, level=c(90,95)) Kas_for_17_kal <- forecast(object = imp_kalman_kasese_na, h = 16, level=c(90,95)) Kas_for_17_ma <- forecast(object = imp_ma_kasese_na, h = 16, level=c(90,95)) Ug_for_17_intp <- forecast(object = imp_interpolation_Uganda_na, h = 16, level=c(90,95)) Ug_for_17_kal <- forecast(object = imp_kalman_Uganda_na, h = 16, level=c(90,95)) Ug_for_17_ma <- forecast(object = imp_ma_Uganda_na, h = 16, level=c(90,95)) Plot_Kas_for_17<-autoplot(Kas_for_17_intp , series="interpolation 2013-2016") + theme_bw() + labs(title=" A Forecast of Typhoid cases in Kasese district for the 1st quarter of 2017", y= "Typhoid cases") + autolayer(Kas_for_17_kal, series="Forecast Kalman") + autolayer(FULL_MOH_2017_ts, series="2017 DATA") + autolayer(fitted(Kas_for_17_ma), series="Forecast Moving Average") + autolayer(fitted(Kas_for_17_intp), series="Forecast Interpolation")+ autolayer(fitted(Kas_for_17_kal), series="fitted") Plot_Ug_for_17<-autoplot(Ug_for_17_intp , series="interpolation 2013-2016") + theme_bw() + labs(title=" A Forecast of Typhoid cases in Uganda for the 1st quarter of 2017", y= "Typhoid cases") + autolayer(Ug_for_17_kal, series="Forecast Kalman") + autolayer(FULL_MOH_2017_UG_ts, series="2017 DATA") + autolayer(fitted(Ug_for_17_ma), series="Forecast Moving Average") + autolayer(fitted(Ug_for_17_intp), series="Forecast Interpolation")+ autolayer(fitted(Ug_for_17_intp), series="fitted") multiplot(Plot_Ug_for_17,Plot_Kas_for_17, cols = 1) pdf("../../Figures/Forecasting.pdf",width = 12,height = 6) multiplot(Plot_Ug_for_17,Plot_Kas_for_17, cols = 1) dev.off() ````