#Libraries required ##### library(dplyr) library(ggplot2) library(tidyr) require(UpSetR) library(janitor) library(clusrank) # Import datasets # e.g. cxdata <- read.csv("...Folder location of downloaded data.../cxdata_update.csv") # data sets required for analysis outbreakdata <- read.csv("...Folder location of downloaded data.../outbreakdata.csv") epidemicdata_spacetime <- read.csv("...Folder location of downloaded data.../epidemicdata_spacetime.csv") regions <- read.csv("...Folder location of downloaded data.../regions.csv") cxdata <- read.csv("...Folder location of downloaded data.../cxdata.csv") riskdata <- read.csv("...Folder location of downloaded data.../riskdata.csv") movementriskdata <- read.csv("...Folder location of downloaded data.../movementriskdata.csv") movementriskdatasanstype <- read.csv("...Folder location of downloaded data.../movementriskdatasanstype.csv") # Analysis # Section 1 - General Descriptive #### #Number of distinct outbreaks incorporated in the total cases #### (n_outbreaks<-length(unique(outbreakdata$outbreakid))) #Number of cases #### (n_cases<-length(unique(outbreakdata$caselogid))) # Number of affected regions and total regions at risk #### length(unique(epidemicdata_spacetime$regiongid)>0) nrow(regions) # Case classification #### outbreakdata %>% group_by(casecode) %>% summarise(total = n()) %>% rowwise() %>% mutate(totalN = n_cases, prop = total/n_cases, lowerci = prop.test(total, n_cases)$conf.int[1], upperci = prop.test(total, n_cases)$conf.int[2]) # Section 2 - Testing categories and days to testing #### outbreakdata %>% group_by(casecode, positivesamples) %>% dplyr::summarise(total = n()) %>% mutate(class = if_else(positivesamples == 'vein serum','serology only','incl agent detection')) %>% ungroup() %>% mutate(totalN = sum(total)) %>% group_by(class,totalN) %>% summarise(total = sum(total)) %>% mutate( per=total/totalN) testingcategory_onlyPCR<- outbreakdata %>% group_by(caselogid , casecode, positivesamples) %>% mutate(class = if_else(grepl("vein serum", positivesamples, fixed=TRUE),'contains serology','PCR only')) %>% ungroup() %>% filter(class == 'PCR only') #Time to sampling #### samplingperioddata<-outbreakdata %>% filter(casecode != 'P3') samplingperioddata.detailed.PCRonly<-inner_join(samplingperioddata, testingcategory_onlyPCR, by = "caselogid") %>% filter(casecode.x == 'P1') samplingperioddata.detailed.PCRonly %>% group_by(casecode.x) %>% summarise( n(), mean = mean(daysdiff_samptocase.x), median = median(daysdiff_samptocase.x), iqr = IQR(daysdiff_samptocase.x), min = min(daysdiff_samptocase.x), max = max(daysdiff_samptocase.x), q0 = quantile(daysdiff_samptocase.x)[1], q25 = quantile(daysdiff_samptocase.x)[2], q50 = quantile(daysdiff_samptocase.x)[3], q75 = quantile(daysdiff_samptocase.x)[4], q100 = quantile(daysdiff_samptocase.x)[5] ) # Section 3 - IP-level descriptive data #### # no. CCs per IP #### outbreakaggdata<-outbreakdata %>% group_by(outbreakid, census) %>% dplyr::summarise(cc_on_ip = n()) median(outbreakaggdata$cc_on_ip) range(outbreakaggdata$cc_on_ip) IQR(outbreakaggdata$cc_on_ip) quantile((outbreakaggdata$cc_on_ip)) ##cumulative total of resident horses that were CCs #### outbreakaggdatacensus <- outbreakaggdata %>% filter(!is.na(census)) %>% mutate(proppos = cc_on_ip/census) outbreakaggdatacensus %>% ungroup() %>% summarise(N_IP_withcensus = n(), sum_cc = sum(outbreakaggdatacensus$cc_on_ip), sum_census = sum(outbreakaggdatacensus$census), prop = sum(outbreakaggdatacensus$cc_on_ip)/sum(outbreakaggdatacensus$census) ) # no. of resident horses on IPs #### quantile(outbreakaggdatacensus$census, na.rm=TRUE) #no. of horses with clinical signs on IPs #### outbreakclindata<-outbreakdata %>% group_by(outbreakid, census, numberclinical) %>% dplyr::summarise(cc_on_ip = n()) quantile(outbreakclindata$numberclinical, na.rm=TRUE) # cumulative total of resident horses with clinical signs #### outbreakaggdataclinical <- outbreakclindata %>% filter(!is.na(census))%>% filter(!is.na(numberclinical))%>% mutate(propclinical = numberclinical/census) outbreakaggdataclinical %>% ungroup() %>% summarise(N_IP_withcensus = n(), sum_clin = sum(numberclinical), sum_census = sum(census), prop = sum(numberclinical)/sum(census) ) # no. of vaccinated resident horses on IPs #### outbreakvaccdata<-outbreakdata %>% group_by(outbreakid, census, numbervaccinated) %>% dplyr::summarise(cc_on_ip = n()) quantile(outbreakvaccdata$numbervaccinated, na.rm=TRUE) # cumulative total of resident horses that were vaccinated #### outbreakaggdatavacc <- outbreakvaccdata %>% filter(!is.na(census))%>% filter(!is.na(numbervaccinated))%>% mutate(propvacc = numbervaccinated/census) outbreakaggdatavacc %>% ungroup() %>% summarise(N_IP_withcensus = n(), sum_numbervacc = sum(numbervaccinated), sum_census = sum(census), prop = sum(numbervaccinated)/sum(census) ) # Number of IPs that had a vaccinated CC (based on case data rather than on property level) #### outbreakvaccdata_detail<-outbreakdata %>% group_by(outbreakid, vaccstatus) %>% dplyr::summarise(totalccvaccinated = n()) %>% filter(vaccstatus == 'vaccinated') nrow(outbreakvaccdata_detail) sum(outbreakvaccdata_detail$totalccvaccinated) # Percentage of IP's with more than one CC sum(ifelse(outbreakaggdata$cc_on_ip > 1, 1, 0)) sum(ifelse(outbreakaggdata$cc_on_ip > 1, 1, 0))/nrow(outbreakaggdata) # premises type code#### epidemicdata_spacetime %>% mutate(premtype = case_when( premisestype=="event" ~ "various non-professional", premisestype=="yard - other" ~ "various non-professional", premisestype=="livery - private" ~ "private", premisestype=="yard - sales" ~ "sales preparation", premisestype=="charity" ~ "various non-professional", premisestype=="horse fair" ~ "various non-professional", premisestype=="sanctuary" ~ "various non-professional", premisestype=="yard - competition - showing" ~ "yard - competition", premisestype=="undefined" ~ "undefined", premisestype=="private" ~ "private", premisestype=="riding school" ~ "riding school", premisestype=="yard - racing" ~ "yard - racing", premisestype=="stud" ~ "stud", premisestype=="yard - competition" ~ "yard - competition", premisestype=="racing pre training" ~ "racing pre training", premisestype=="livery" ~ "livery", premisestype=="training" ~ "training")) %>% group_by(premtype) %>% summarise(total = n()) %>% rowwise() %>% mutate(totalN = n_outbreaks, prop = total/n_outbreaks, lowerci = prop.test(total, n_outbreaks)$conf.int[1], upperci = prop.test(total, n_outbreaks)$conf.int[2]) # biosecurity measures on IPs (isolation, quarantine, equipment, tack) #### riskdata.twolevel<-riskdata %>% filter(riskclass %in% c("biosecurity","management")) %>% filter(quantity !='unknown') riskdata.twolevel %>% group_by(riskclass, risksubclass, quantity) %>% summarise(n()) # basic risk categories coded#### riskdata.twolevel<-riskdata.twolevel %>% mutate(quantity2 = ifelse(quantity %in% c('none','minority','half'), 'limited',(ifelse(quantity %in% c('all','majority'), 'extensive',quantity)))) riskfactors.twolevel<-list() for (i in unique(riskdata.twolevel$riskclass)) { riskfactors.twolevel[[i]]<-riskdata.twolevel %>% filter(riskclass == i) %>% group_by(risksubclass, quantity2) %>% summarise(total = n()) %>% group_by(risksubclass) %>% mutate(totalT = sum(total)) %>% group_by(quantity2, add=TRUE) %>% mutate( per=total/totalT, lowerci = prop.test(total, totalT)$conf.int[1], upperci = prop.test(total, totalT)$conf.int[2]) %>% select(risksubclass, quantity2, totalT, total, per, lowerci, upperci) } riskfactors.twolevel # One level risk factors from table - Notable Equine Gathering (ne_gather) and movement #### riskdata.onelevel<-riskdata %>% filter(riskclass %in% c("event - ne_gather","movement")) %>% filter(quantity !='unknown') %>% mutate(quantity2 = ifelse(quantity %in% c('none','minority','half'), 'limited',(ifelse(quantity %in% c('all','majority'), 'extensive',quantity)))) riskfactorsonelevel<-list() ## those IPs with direct or indirect link with Notable Equine Gathering and IPs where a new arrival was a case #### for (i in unique(riskdata.onelevel$riskclass)) { riskfactorsonelevel[[i]]<-riskdata.onelevel %>% filter(riskclass == i) %>% group_by(risksubclass) %>% summarise(total = n()) %>% mutate(totalT = n_outbreaks) %>% group_by(risksubclass) %>% mutate( per=total/totalT, lowerci = prop.test(total, totalT)$conf.int[1], upperci = prop.test(total, totalT)$conf.int[2]) %>% select(risksubclass, totalT, total, per, lowerci, upperci) } riskfactorsonelevel # Movement Risk Factors #### movementriskdata<-movementriskdata %>% mutate(first2weekstotal = ifelse(timeperiod %in% c('<2 weeks', '< 1 week'), 'twoweekstotal','outsidetwoweeks')) movementriskfactors<-list() for (i in unique(movementriskdata$movementcategoryanimalclass)) { movementriskfactors[[i]]<-movementriskdata %>% filter(movementcategoryanimalclass == i) %>% group_by(timeperiod) %>% summarise(total = n()) %>% mutate(totalT = n_outbreaks) %>% rowwise() %>% mutate( per=total/totalT, lowerci = prop.test(total, totalT)$conf.int[1], upperci = prop.test(total, totalT)$conf.int[2]) %>% select(timeperiod, totalT, total, per, lowerci, upperci) movementriskfactors[[paste('total',i)]]<-movementriskdata %>% filter(movementcategoryanimalclass == i) %>% filter(first2weekstotal == 'twoweekstotal') %>% group_by(first2weekstotal) %>% summarise(total = n()) %>% mutate(totalT = n_outbreaks) %>% rowwise() %>% mutate( per=total/totalT, lowerci = prop.test(total, totalT)$conf.int[1], upperci = prop.test(total, totalT)$conf.int[2]) %>% select(first2weekstotal, totalT, total, per, lowerci, upperci) } movementriskfactors #Section 4 - CC level descriptive data #### # Totals per age #### outbreakdata %>% filter(!is.na(age)) %>% summarise( n(), mean(age), IQR(age), q0 = quantile(age)[1], q25 = quantile(age)[2], q50 = quantile(age)[3], q75 = quantile(age)[4], q100 = quantile(age)[5] ) # Totals per general breed #### outbreakdata %>% group_by(breedgen) %>% summarise(total = n()) %>% rowwise() %>% mutate(totalN = n_cases, prop = total/n_cases) %>% arrange(., -total) # Totals per breed - in supplementary table #### breed<-outbreakdata %>% group_by(breedgen, breed) %>% summarise(total = n()) %>% rowwise() %>% mutate(totalN = n_cases, prop = total/n_cases) %>% arrange(., -total) # Sex General ##### # Recode to general sex #### outbreakdata %>% mutate(sexgroup = case_when(is.na(sexgen) ~ "missing", sexgen=="male non-intact" ~ "male", sexgen=="male general" ~ "male", sexgen=="male intact" ~ "male", sexgen=="female" ~ "female", sexgen=="unknown" ~ "unknown" )) %>% group_by(sexgroup) %>% summarise(total = n()) %>% rowwise() %>% mutate(totalN = n_cases, prop = total/n_cases) %>% arrange(., -total) # Vaccination details #### # Proportion of vaccination status from outbreak data outbreakdata %>% group_by(vaccstatus) %>% dplyr::summarise(n = n()) %>% mutate(totalN = n_cases, prop = n / sum(n)) %>% arrange(., -n) # evaluate days since last vaccination #### # limited to an existing vacc date and where its classified as vaccinated to exclude those V1's where the vet vaccinated in the face of the outbreak dayssincelastvacc<-outbreakdata %>% filter(vaccstatus == 'vaccinated' & !is.na(vaccdate) & vaccdate != '1900-01-01') %>% arrange(dayssincevaccination) #Vacc cases - median no. per outbreak #### temp<-outbreakdata %>% group_by(outbreakid, vaccstatus) %>% dplyr::summarise(n = n()) temp2<-temp %>% pivot_wider( names_from = vaccstatus, values_from = n, values_fill = list(n = 0) ) sum(temp2$vaccinated>0) # number of outbreaks with a vaccinated case # no per outbreak with prop and CI (vacc.cc.onPI<-prop.test(sum(temp2$vaccinated>0), n_outbreaks)) #mean number of vaccinated cases per outbreak #### temp3 = janitor::clean_names(temp2) temp3$totalN<-apply(temp3[,-1], 1, sum) #sum across columns temp3$propvacc<- temp3$vaccinated/temp3$totalN #prop vacc quantile(temp3$vaccinated) #Aggregate data by holding to account for clustering tempnew<-as.data.frame(dayssincelastvacc %>% filter(!is.na(vaccclass)) %>% filter(vaccclass != "unknown")%>% filter(vaccclass != "V1") %>% filter(dayssincevaccination >7) %>% mutate(vaccgroup = case_when(is.na(vaccclass) ~ "missing", vaccclass=="V2" ~ "V2", vaccclass=="V3" ~ "V3", vaccclass=="1st booster after primary" ~ "booster (annual)", vaccclass=="booster (6mo)" ~ "booster (6mo)", vaccclass=="booster (annual)" ~ "booster (annual)" )) %>% group_by(outbreakid, vaccgroup) %>% summarise(total = n(), meanip = mean(dayssincevaccination))) length(unique(tempnew$outbreakid)) outcome<-tempnew %>% group_by(vaccgroup) %>% summarise( totalip = n(), meantotalperip = mean(total), mean = mean(meanip), median = median(meanip), IQR = IQR(meanip), min = min(meanip), max = max(meanip), sd = sd(meanip)) outcomeCombo<-tempnew %>% summarise( totalip = n(), meantotalperip = mean(total), mean = mean(meanip), median = median(meanip), IQR = IQR(meanip), min = min(meanip), max = max(meanip), sd = sd(meanip)) tempAge<-outbreakdata %>% filter(!is.na(age)) %>% group_by(outbreakid, vaccstatus) %>% summarise( IPtotal = n(), meanIP = mean(age), medianIP = median(age) ) length(unique(tempAge$outbreakid)) tempAge %>% group_by(vaccstatus) %>% summarise( ips = n(), meanTotal = mean(IPtotal), meanmeanIP = mean(meanIP), meanmedianIP = mean(medianIP), q0 = quantile(meanIP)[1], q25 = quantile(meanIP)[2], q50 = quantile(meanIP)[3], q75 = quantile(meanIP)[4], q100 = quantile(meanIP)[5] ) #vacc group and age - not normally distributed - established during EDA #### outbreakvaccdata <- outbreakdata %>% filter(!is.na(age))%>% filter(vaccstatus=="vaccinated") #unvacc group and age outbreakunvaccdata <- outbreakdata %>% filter(!is.na(age))%>% filter(vaccstatus=="unvaccinated") #compare using non-parametric test to compare median age of groups as not normally distributed - mann whitney wilcoxon test outbreakvaccdata.cluster <- outbreakdata %>% filter(!is.na(age) & vaccstatus %in% c('vaccinated','unvaccinated')) %>% select (outbreakid, vaccstatus, age) clusWilcox.test(age ~ vaccstatus + cluster(outbreakid), data = outbreakvaccdata.cluster, method = "rgl") # section 5 - clinical signs in CCs ########## cxdataall_summ <-cxdata #remove unreported cx - for upsetGraph cxdata<-cxdata %>% filter(cx != 'unreported') %>% droplevels() cxdataall_summ %>% mutate(cxsumm = ifelse(cx == 'unreported', 'unreported',(ifelse(cx == 'no clinical signs', 'subclinical','reported')))) %>% distinct(caselogid, cxsumm) %>% mutate(totalN = n()) %>% group_by(cxsumm,totalN) %>% summarise(total = n()) %>% mutate( per=total/totalN, lowerci = prop.test(total, totalN)$conf.int[1], upperci = prop.test(total, totalN)$conf.int[2]) #now overall counts for reported or not cxdataall_summ %>% mutate(cxsumm = ifelse(cx == 'unreported', 'unreported',(ifelse(cx == 'no clinical signs', 'subclinical','reported')))) %>% filter(cxsumm != 'unreported') %>% distinct(caselogid, cxsumm) %>% mutate(totalN = n()) %>% group_by(cxsumm,totalN) %>% summarise(total = n()) %>% mutate( per=total/totalN, lowerci = prop.test(total, totalN)$conf.int[1], upperci = prop.test(total, totalN)$conf.int[2]) #setval for upSetR graph setval<- as.integer(cxdata %>% group_by(cx) %>% summarise(total = n()) %>% arrange(desc(total)) %>% head(1) %>% .$total*1.5) cxdata.spread<-spread(as.data.frame(table(cxdata)), cx, Freq) upset(cxdata.spread, nsets = nrow(cxdata %>% filter(cx!='unreported') %>% group_by(cx) %>% summarise(count = n()) %>% filter(count>=10)), nintersects = nrow(aggregate(cx ~ caselogid, data = cxdata %>% filter(cx!='unreported'), c) %>% group_by(cx) %>% summarise(count = n()) %>% filter(count>=10) %>% arrange(-count)), mb.ratio = c(0.5, 0.5), order.by = c("degree","freq"), decreasing = c(FALSE,TRUE), matrix.color = "gray23", main.bar.color = "#109a49", mainbar.y.label = "Frequency of\n clincial sign combinations", mainbar.y.max = NULL, sets.bar.color = "gray23", sets.x.label = "Total cases", set_size.show = TRUE, text.scale = 2, set_size.scale_max = setval ) # For the specifics of when ND, pyrexia and Coughing cxdata.spread$total = rowSums(cxdata.spread[,-1]) nrow(cxdata.spread %>% filter(`nasal discharge` == 1 & coughing == 1 & pyrexia == 1 & total == 3))/nrow(cxdata.spread) # the subclincial cases #### horseswithnocx<-cxdata %>% filter(cx == 'no clinical signs') %>% select(caselogid) %>% distinct() outbreakdata %>% inner_join(horseswithnocx) # Section 6 - epidemic curve ####### #convert the casedate to a week of outbreak epidemicdata_spacetime$outbreakweek<-as.integer( difftime( epidemicdata_spacetime$outbreakdate, min(as.Date(epidemicdata_spacetime$outbreakdate)), units = "weeks") ) + 1 # add the classification of outbreaks as to whether they are in a new area epidemicdataclassed<-epidemicdata_spacetime %>% arrange(outbreakdate) %>% group_by(regiongid) %>% mutate(arearank = rank(outbreakdate, ties.method = 'first')) %>% mutate(outbreakclass = ifelse(arearank == 1, "newregion", "oldregion")) %>% arrange(outbreakweek) %>% # get the total number of outbreaks per week group_by(outbreakweek, outbreakclass) %>% summarise(totaloutbreaksperweek = n()) # arrange the spatial progress data by outbreak week outbreakprogressepidemicdata<-epidemicdata_spacetime %>% arrange(outbreakweek)%>% group_by(regiongid) %>% summarise(minoutbreakdate=min(as.Date(outbreakdate))) # add outbreak week to data outbreakprogressepidemicdata$outbreakweek<-as.integer( difftime( outbreakprogressepidemicdata$minoutbreakdate, min(as.Date(epidemicdata_spacetime$outbreakdate)), units = "weeks") ) + 1 # groupby week outbreakprogressepidemicdata<- outbreakprogressepidemicdata %>% group_by(outbreakweek) %>% summarise(noofnewareasperweek=n()) %>% mutate(areasprogress = cumsum(noofnewareasperweek)) # levels for y-axis yaxislevels<-epidemicdataclassed %>% group_by(outbreakweek) %>% summarise(totalsforgraph = sum(totaloutbreaksperweek)) outbreakprogressepidemicdata.month<-outbreakprogressepidemicdata %>% mutate(month = floor((outbreakweek-1)/4)+1) %>% group_by(month) %>% summarise(totalpermonth = sum(noofnewareasperweek))%>% complete(month = seq(0, floor((max(epidemicdataclassed$outbreakweek)-1)/4)+1)) %>% mutate(totalpermonth = replace_na(totalpermonth, 0)) %>% mutate(plotpoint = ifelse(month == 0, 1, (month)*4)) %>% mutate(phase = ifelse(month <= 1,'Week 1 and Phase 1',(ifelse(month >=5 & month <= 6,'Phase 2','Intermediate phases')))) #need to complete the dataset to include 0 new outbreaks between the last month of new areas and the last week of the actual outbreak outbreakprogressepidemicdata.month$phase <- factor(outbreakprogressepidemicdata.month$phase, levels=c("Week 1 and Phase 1", "Intermediate phases", "Phase 2")) ggplot(epidemicdataclassed) + geom_bar(aes(x = outbreakweek, y = totaloutbreaksperweek), stat = 'identity', fill = "#FFDB6D") + xlab(paste("Week of outbreak - starting",format(min(as.Date(epidemicdata_spacetime$outbreakdate)), '%d %b %Y'))) + ylab("Number of infected premises (bars) \n4 week cumulative new counties affected (lines)") + scale_x_continuous(breaks = seq(from= 1, to = max(epidemicdataclassed$outbreakweek)+4,by = 4), labels = seq(from= 1, to = max(epidemicdataclassed$outbreakweek)+4,by = 4)) + scale_y_continuous(expand = c(0,0), breaks=seq(0,max(yaxislevels$totalsforgraph)+2,5), limits = c(0,max(yaxislevels$totalsforgraph)+3)) + theme_classic() + theme(axis.text.x = element_text(), axis.text = element_text(size = 15, vjust=0, color = '#000000'), legend.text = element_text(size = 12), legend.title = element_text(size = 15), axis.title = element_text(size = 20, color = '#000000'), plot.title = element_text(size = 20, color = '#000000'), panel.background = element_blank() ) + geom_line(data=outbreakprogressepidemicdata.month, aes(x=plotpoint, y=totalpermonth, colour = phase), size = 1.5) + scale_colour_manual(values = c("#ff8525","#deebf7","#91522d")) + geom_point(data=outbreakprogressepidemicdata.month, aes(x=plotpoint, y=totalpermonth, fill = phase), size =2, shape=21) + scale_fill_manual(values = c("white","white","white")) + labs(color = "Phases of outbreak", fill = "Phases of outbreak") + theme(legend.position="bottom") + aes(group=NA) # Section 7 - modelling section on risk - this code to establish the dataset #### # First create large dataframe of n_outbreaks with risk column and in which phase the outbreak first occurred #get various levels: unique(riskdata.twolevel$risksubclass) # isolation risk.isolation<-riskdata.twolevel %>% filter(risksubclass == 'isolation facility') %>% dplyr::select(outbreakid, isolation = quantity2) #quarantine new arrival risk.quarantine<-riskdata.twolevel %>% filter(risksubclass == 'quarantine new arrival') %>% dplyr::select(outbreakid, quarantine = quantity2) #sharetack risk.sharetack<-riskdata.twolevel %>% filter(risksubclass == 'share tack') %>% dplyr::select(outbreakid, sharetack = quantity2) #share equipment risk.shareequipment<-riskdata.twolevel %>% filter(risksubclass == 'share equipment') %>% dplyr::select(outbreakid, shareequipment = quantity2) #notable equine gathering risk.ne_gather<-riskdata.onelevel %>% filter(riskclass == 'event - ne_gather') %>% dplyr::select(outbreakid) %>% mutate(ne_gather = 1) #newarrivals risk.newarrivalisacase<-riskdata.onelevel %>% filter(riskclass == 'movement') %>% dplyr::select(outbreakid) %>% mutate(newarrivals = 1) #recentmovement risk.recentmovement<-movementriskdatasanstype %>% mutate(first2weekstotal = ifelse(timeperiod %in% c('<2 weeks', '< 1 week'), 'twoweekstotal','outsidetwoweeks')) %>% dplyr::select(outbreakid, movementcategorydirection, first2weekstotal) %>% distinct() %>% filter(first2weekstotal == 'twoweekstotal') %>% dplyr::select(outbreakid) %>% mutate(move_in_last2weeks = 1) #phases of all outbreaks risk.phase<-epidemicdata_spacetime %>% mutate(phase = ifelse(outbreakweek <= 14,0, (ifelse(outbreakweek <=18, NA, (ifelse(outbreakweek <=34, 1, NA)))))) %>% dplyr::select(outbreakid, regiongid, premisestype, phase, outbreakdate) #One outbreak on 2 May was allocated to Phase 2 for statistical evaluation sum(risk.phase$phase==0, na.rm = TRUE) sum(risk.phase$phase==1, na.rm = TRUE) sum(is.na(risk.phase$phase), na.rm = TRUE) risk.phase[risk.phase$outbreakdate=='2019-05-02',]$phase <- 1 #now to join all tables riskdata.model<-left_join(risk.phase, risk.isolation) %>% left_join(., risk.quarantine) %>% left_join(., risk.sharetack) %>% left_join(., risk.shareequipment) %>% left_join(., risk.ne_gather) %>% left_join(., risk.newarrivalisacase) %>% left_join(., risk.recentmovement) #Recoding riskdata.model$regiongid<-as.character(riskdata.model$regiongid) riskdata.model$move_in_last2weeks[is.na(riskdata.model$move_in_last2weeks)] <- FALSE riskdata.model$ne_gather[is.na(riskdata.model$ne_gather)] <- 0 riskdata.model$newarrivals[is.na(riskdata.model$newarrivals)] <- FALSE write.csv(riskdata.model,'riskdatamodel.csv', row.names = FALSE)