# 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 4. Seroprevalence at test-and-slaughter feedgrounds ### 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() library(plyr) library(arm) library(mgcv) library(ggplot2) ### data includes females only, aged >=1.5 years old at capture sdata<-read.csv("figure 4 test and slaughter.csv",header=T) ### generate estimates and confidence intervals by feedground source("agg_sero_fxn.R") sero.agg.loc.sub <- agg.sero(year = sdata$Capture.Year, state = sdata$State, disease = sdata$Interp, region = sdata$Capture.Location) ### The main purpose of these gams is to "connect the dots" for interpretability of general seroprevalence # trends at a site over time, nothing more. We invoke parsimony just as a general practice. ### Fall ### fal.all<-subset(sdata, Capture.Location=="Fall Creek") f3=gam(Interp~s(Capture.Year, k=3),family=binomial, data=fal.all) f4=gam(Interp~s(Capture.Year, k=4),family=binomial, data=fal.all) f5=gam(Interp~s(Capture.Year, k=5),family=binomial, data=fal.all) f6=gam(Interp~s(Capture.Year, k=6),family=binomial, data=fal.all) f7=gam(Interp~s(Capture.Year, k=7),family=binomial, data=fal.all) f8=gam(Interp~s(Capture.Year, k=8),family=binomial, data=fal.all) f9=gam(Interp~s(Capture.Year, k=9),family=binomial, data=fal.all) anova(f3,f4,f5,f6,f7,f8,f9, test="Chisq") ### the model with 3 knots (model 1 in the output) is as good as any other. plot(f3) new.data=data.frame(Capture.Year=seq(1993,2017,1)) # just do this once new.data$fall=predict.gam(f3, newdata=new.data, type="response") plot(new.data$fall) ### Scab ### sca.all<-subset(sdata, Capture.Location=="Scab Creek") s3=gam(Interp~s(Capture.Year, k=3),family=binomial, data=sca.all) s4=gam(Interp~s(Capture.Year, k=4),family=binomial, data=sca.all) s5=gam(Interp~s(Capture.Year, k=5),family=binomial, data=sca.all) s6=gam(Interp~s(Capture.Year, k=6),family=binomial, data=sca.all) s7=gam(Interp~s(Capture.Year, k=7),family=binomial, data=sca.all) s8=gam(Interp~s(Capture.Year, k=8),family=binomial, data=sca.all) s9=gam(Interp~s(Capture.Year, k=9),family=binomial, data=sca.all) s10=gam(Interp~s(Capture.Year, k=10),family=binomial, data=sca.all) anova(s3,s4,s5,s6,s7,s8,s9,s10, test="Chisq") plot(s3) ### The model with 3 knots (model 1 in the output) is as good as any other new.data$scab=predict.gam(s3, newdata=new.data, type="response") plot(new.data$scab) ### muddy ### mud.all<-subset(sdata, Capture.Location=="Muddy Creek") m3=gam(Interp~s(Capture.Year, k=3),family=binomial, data=mud.all) m4=gam(Interp~s(Capture.Year, k=4),family=binomial, data=mud.all) m5=gam(Interp~s(Capture.Year, k=5),family=binomial, data=mud.all) m6=gam(Interp~s(Capture.Year, k=6),family=binomial, data=mud.all) m7=gam(Interp~s(Capture.Year, k=7),family=binomial, data=mud.all) m8=gam(Interp~s(Capture.Year, k=8),family=binomial, data=mud.all) m9=gam(Interp~s(Capture.Year, k=9),family=binomial, data=mud.all) m10=gam(Interp~s(Capture.Year, k=10),family=binomial, data=mud.all) anova(m3,m4,m5,m6,m7,m8,m9,m10, test="Chisq") ### The model with 6 knots (model 4 in the output) achieves the best balance of # low res. dev. and fewest knots plot(m6) new.data$mud=predict.gam(m6, newdata=new.data, type="response") plot(new.data$mud) ### save estimates names(new.data)[2:4]=c("f.fall","f.scab","f.mud") Sero.Estimate=reshape(new.data,varying = 2:4,direction="long",idvar = "Capture.Year") names(Sero.Estimate)[1:3]=c("year","region","Estimate") Sero.Estimate[which(Sero.Estimate$region=="fall"),2]="Fall Creek" Sero.Estimate[which(Sero.Estimate$region=="scab"),2]="Scab Creek" Sero.Estimate[which(Sero.Estimate$region=="mud"),2]="Muddy Creek" ### check: test<-merge(Sero.Estimate, sero.agg.loc.sub, all=T) sero.agg.loc.sub<-test ### could plot all data, but looks better to be able to zoom in on the relevant years. Omitting earlier years ### doesn't affect interpretation. Subset further: sero.agg.loc.sub<-subset(sero.agg.loc.sub, year>2004) ### create data frame with treatment years: m1_line_data<-data.frame(xint=c(2006:2010), region="Muddy Creek") s1_line_data<-data.frame(xint=c(2009:2010), region="Scab Creek") f1_line_data<-data.frame(xint=c(2008:2009), region="Fall Creek") line_data <- do.call("rbind",list(m1_line_data, s1_line_data, f1_line_data)) ### specify order of panels by resetting factor levels: neworder <- c("Muddy Creek","Scab Creek","Fall Creek") sero.agg.loc.sub <- arrange(transform(sero.agg.loc.sub, region=factor(region,levels=neworder)),region) ### plot titlesize<-20 axistitlesize<-titlesize prevsize<-15 tas3<-ggplot(data=sero.agg.loc.sub, aes(x=year, y=Prev))+ geom_point(size=1.5)+ geom_errorbar(data = sero.agg.loc.sub, mapping=aes(x=year, ymin=CI95lo, ymax=CI95hi), width=0, size=1.1, color="black")+ geom_line(data=sero.agg.loc.sub, mapping=aes(x=year,y=Estimate),size=0.8, color="red")+ geom_vline(data=line_data, aes(xintercept=xint), linetype="dotted",color="grey50",size=0.9) + facet_wrap(~region)+ scale_x_continuous("",seq(2000,2015,5),limits = c(2005,2018))+ scale_y_continuous(limits = c(0,1))+ ylab("Seroprevalence")+ theme_bw()+ theme(axis.title.y = element_text(size=axistitlesize, margin = margin(r=10,l=10)), axis.text.x = element_text(size = axistitlesize, angle = 45, margin=margin(t=30)), #year values axis.text.y = element_text(size = prevsize), #prevalence strip.text = element_text(size=titlesize, hjust=-0), strip.background = element_blank(), panel.border = element_blank(), panel.grid.major = element_blank(), panel.grid.minor = element_blank())+ annotate("segment", x=-Inf, xend=Inf, y=-Inf, yend=-Inf)+ annotate("segment", x=-Inf, xend=-Inf, y=-Inf, yend=Inf) tas3 ### to write TIFF: # tiff("figure 4 test and slaughter2.tiff",width=3000, height = 1000,units="px",res=275) # tas3 # dev.off() ### end