# This software has been approved for release by the U.S. Geological Survey (USGS). # Although the software has been subjected to rigorous review, the USGS reserves # the right to update the software as needed pursuant to further analysis and # review. No warranty, expressed or implied, is made by the USGS or the U.S. # Government as to the functionality of the software and related material nor # shall the fact of release constitute any such warranty. Furthermore, the # software is released on condition that neither the USGS nor the U.S. Government # shall be held liable for any damages resulting from its authorized or # unauthorized use. ### Last updated: 24Aug2017 ### Gavin Cotterill, gavin.cotterill@aggiemail.usu.edu ### Figure 5. Dell Creek and Greys River plots illustrating the correlation between seroprevalence and ### previous 8-year average end date ### clear environment rm(list=ls()) dev.off() ### set working directory # if data are stored in same folder as code and you're using RStudio, # these lines obviate the need to change any filepaths. Requires package "rstudioapi". # Works on most recent versions of R and RStudio for Windows or OSx as of summer 2017. # Otherwise, specify the correct filepath for read.csv() and source() lines. rstudioapi::getActiveDocumentContext setwd(dirname(rstudioapi::getActiveDocumentContext()$path)) getwd() ### read in end dates ### 'End' is the day, for a given site and year, that feeding ceased, where Nov 1 of the previous year = 1. end<-read.csv("figure 5 end dates.csv",header=T) ### calculate average end date for a site in the previous 8 years region.list<-unique(end$Feedground) year.list<-unique(end$Year) outd<-data.frame(NA) for(i in 1:length(region.list)){ temp<-subset(end, Feedground==region.list[i]) temp$P8YE<-NA for(j in 9:nrow(temp)){ x8<-temp[j-1:8,"End"] temp$P8YE[j]=mean(x8,na.rm=T) } outd<-merge(outd,temp,all=T) } outd<-outd[,-which(names(outd)=="NA.")] end<-outd rm(list=c("outd","temp")) ### in order to plot using calendar days, create new variable: ### first, round end$P8YEdate<-round(end$P8YE,0) ### subtract 61 days to account for Nov1=1 end$P8YEdate<-end$P8YEdate-61 # convert to appropriate format for plotting, year is now meaningless, these are average julian dates end$P8YEdate<-as.Date(end$P8YEdate, format="%m/%d/%Y", origin="01/01/1990") ### read in serology data sdata<-read.csv("figure 5 serology.csv",header=T) ### calculate prevalence estimates and confidence intervals by site source("agg_sero_fxn.R") sero.agg.loc <- agg.sero(year = sdata$Capture.Year, state = sdata$State, disease = sdata$Interp, region = sdata$Capture.Location) source("pred_sero_fxn_edit.R") ### define gam terms knumber<--1 splinetype<-"tp" fam<-binomial pred.loc.gam<- pred.sero(year = sdata$Capture.Year, state = sdata$State, region = sdata$Capture.Location, disease = sdata$Interp, line = F, gam.smooth = T) sero.agg.loc<-merge(sero.agg.loc, pred.loc.gam, all=T) ### prep data for plotting: ### dell dell.smooth<-subset(sero.agg.loc, region=="Dell Creek") dell.smooth<-subset(dell.smooth, year>1997) dell.ends<-subset(end, Feedground=="Dell Creek") dell.ends<-subset(dell.ends, Year>1997) names(dell.ends)[1:2]=c("region","year") dell<-merge(dell.ends, dell.smooth, all=T) rm(list=c("dell.ends", "dell.smooth")) ### greys grey.smooth<-subset(sero.agg.loc, region=="Greys River") grey.ends<-subset(end, Feedground=="Greys River") grey.ends<-subset(grey.ends, Year>1992) names(grey.ends)[1:2]=c("region","year") greys<-merge(grey.ends, grey.smooth, all=T) rm(list=c("grey.ends", "grey.smooth")) dev.off() ###################### ### generate plot: ### ###################### # plot fed sites: par(mar = c(5,5,2,5), mfrow=c(1,2)) ### Dell Creek ### ### raw prevalence estimate with(dell, plot(year, Prev, col="red3", pch=1, xlab=NA, ylab="Seroprevalence Estimate",cex=N/25, xlim=c(1998,2017))) mtext(side = 3, 'Dell Creek') ### smoother with(dell, lines(year, prev, col="red3")) ### end dates par(new = T) with(dell, plot(year, as.Date(P8YEdate, format="%d-%b"), pch=16, axes=F, xlab=NA, ylab=NA, cex=1.2)) axis(side = 4, end$P8YEdate, format(end$P8YEdate, "%d-%b")) mtext(side = 4, line = 3, 'Average End Date') ### Greys River ### ### raw prevalence estimate with(greys, plot(year, Prev, col="red3", pch=1, xlab=NA, ylab="Seroprevalence Estimate",cex=N/25, xlim=c(1993,2017))) mtext(side = 3, 'Greys River') ### smoother with(greys, lines(year, prev, col="red3")) ### end dates par(new = T) with(greys, plot(year, as.Date(P8YEdate, format="%d-%b"), pch=16, axes=F, xlab=NA, ylab=NA, cex=1.2)) axis(side = 4, end$P8YEdate, format(end$P8YEdate, "%d-%b")) mtext(side = 4, line = 3, 'Average End Date') ###################### ### to write tiff: ### ###################### # tiff("figure 5 end dates.tiff",width=3000, height = 1200,units="px",res=275) # # # plot fed sites: # par(mar = c(5,5,2,5), mfrow=c(1,2)) # # ### Dell Creek ### # ### raw prevalence estimate # with(dell, plot(year, Prev, col="red3", pch=1, xlab=NA, # ylab="Seroprevalence Estimate",cex=N/25, # xlim=c(1998,2017))) # mtext(side = 3, 'Dell Creek') # ### smoother # with(dell, lines(year, prev, col="red3")) # # ### end dates # par(new = T) # with(dell, plot(year, as.Date(P8YEdate, format="%d-%b"), pch=16, axes=F, xlab=NA, ylab=NA, cex=1.2)) # axis(side = 4, end$P8YEdate, format(end$P8YEdate, "%d-%b")) # mtext(side = 4, line = 3, 'Average End Date') # # # ### Greys River ### # ### raw prevalence estimate # with(greys, plot(year, Prev, col="red3", pch=1, xlab=NA, # ylab="Seroprevalence Estimate",cex=N/25, # xlim=c(1993,2017))) # mtext(side = 3, 'Greys River') # ### smoother # with(greys, lines(year, prev, col="red3")) # # ### end dates # par(new = T) # with(greys, plot(year, as.Date(P8YEdate, format="%d-%b"), pch=16, axes=F, xlab=NA, ylab=NA, cex=1.2)) # axis(side = 4, end$P8YEdate, format(end$P8YEdate, "%d-%b")) # mtext(side = 4, line = 3, 'Average End Date') # # dev.off() ### end