############## # Notes ############## # Code used for figures and analyses in # Stronen, A. V., Molnar, B., Ciucci, P., Darimont, C. T., Grottoli, L., Paquet, P. C., Sallows, T., Smmits, J.E.G., and Bryan, H. M. 2021. Cross-continental comparison of parasite communities in a wide-ranging carnivore suggests associations with prey diversity and host density. Ecology and Evolution. # Deposted on datadryad.org # Need to use R 3.3.3 for modeling # See README.txt for variable definitions ################################################################################################################################################# ################################################################################################################################################# # Set working directory setwd("C:/Users/") setwd("C:/Users/Heather/Wolf/wolfparasitepaper/astrid_comparison/Revisions for Ecology and Evolution/Dryad") # Read in data finaldb<-read.csv("wolfParasiteData.csv",header=T) # Load Packages library(plyr) library(dplyr) library("plotrix") library(R2admb) library(glmmADMB) library(MuMIn) ##################################################################################################################################################################### ##################################################################################################################################################################### # Figure 2: Barplot showing direct and indirect parasites by study area pdf("Fig2_Direct and indirect parasites by study area.pdf", width = 7.5, height = 3.5,pointsize = 10,bg = "white",useDingbats=F) par(mfrow=c(1,2),mgp=c(1.5,0.5,0),oma=c(0,0,0,0),mar=c(2,3,1,0.8),bty="l",las=1) IT <-table(finaldb$DirectParasites,finaldb$Study.Area) ITmat<-mat.or.vec(nrow(IT),ncol(IT)) nVec <- mat.or.vec(ncol(IT),1) for (i in 1:ncol(IT)) { n <- sum(IT[,i]) ITmat[,i] = IT[,i] / n * 100 nVec[i] <- n } #reorder by increasing wolf density ITmat2<-ITmat[,c(3,1,5,2,6,4)] nVec2<-nVec[c(3,1,5,2,6,4)] barpos<-barplot(ITmat2[2,], beside=F, horiz=F,xlab="",ylab="Percent of Faeces", main="Direct",ylim=c(0,100),col=c("white","grey","grey","grey","black","black")) text(barpos,(ITmat2[2,]+2),nVec2,cex=0.8) xlab<-c("MNP","DMPPF","RMNP","GBR","YNP","PNALM") axis(1,at=barpos,line=-0.5, labels=xlab,cex=0.8,cex.lab=0.8,cex.axis=0.8,las=1,tick=FALSE,lwd=0,padj=0) legend("topleft",c("low 11.5","medium 25-35","high 50"),fill=c("white","grey","black"),title="Wolf Density/1000 kmsq",bty="n",cex=0.8) mtext("Study Area",1,line=1,las=1,cex=1) mtext("A",2,line=2,at=103,las=1,cex=1.2) IT <-table(finaldb$IndirectParasites,finaldb$Study.Area) ITmat<-mat.or.vec(nrow(IT),ncol(IT)) nVec <- mat.or.vec(ncol(IT),1) for (i in 1:ncol(IT)) { n <- sum(IT[,i]) ITmat[,i] = IT[,i] / n * 100 nVec[i] <- n } ITmat2<-ITmat[,c(6,2,5,1,3,4)] nVec2<-nVec[c(6,2,5,1,3,4)] barpos<-barplot(ITmat2[2,], beside=F, horiz=F,xlab="",ylab="Percent of Faeces", main="Indirect",ylim=c(0,100),col=c("white","grey","grey","grey","black","black")) text(barpos,(ITmat2[2,]+2),nVec2,cex=0.8) xlab<-c("YNP","GBR","RMNP","DMPPF","MNP","PNALM") axis(1,at=barpos,line=-0.5,labels=xlab,cex=0.8,cex.lab=0.8,cex.axis=0.8,las=1,tick=FALSE,lwd=0,padj=0) legend("topright",c("low","medium","high"),fill=c("white","grey","black"),title="Diet Diversity",bty="n",cex=0.8) mtext("Study Area",1,line=1,las=1,cex=1) mtext("B",2,line=2,at=103,las=1,cex=1.2) dev.off() ########################################################################################################## ########################################################################################################## # Figure 3: Barplot showing prevalence of each of the direct parasites by study area pdf("Fig3_Direct and indirect parasites separate by study area.pdf", width = 7.5, height = 3.5,pointsize = 10,bg = "white",useDingbats=F) #5.5 inches is 140mm or 1.5 columns par(mfrow=c(1,2),bty="n",las=1,mgp=c(1.5,0.5,0),oma=c(0,0,0,0),mar=c(2,2.5,1,0.8)) directParasites<-group_by(finaldb, Study.Area) %>% summarise(number_samples=length(Study.Area),ToxascarisCount=sum(Toxascaris_binomial), CapillariaCount=sum(Capillaria_binomial),UncinariaCount=sum(Uncinaria.stenocephala_binomial), CoccidiaCount=sum(Coccidia_binomial),ToxocaraCount=sum(Toxocara.canis_binomial),TrichurisCount=sum(Trichuris.vulpis_binomial),SpiroideaCount=sum(Spiroidea_binomial)) %>% mutate(ToxascarisPct=ToxascarisCount/number_samples*100,CapillariaPct=CapillariaCount/number_samples*100,UncinariaCountPct=UncinariaCount/number_samples*100, CoccidiaPct=CoccidiaCount/number_samples*100,ToxocaraPct=ToxocaraCount/number_samples*100,TrichurisPct=TrichurisCount/number_samples*100,SpiroideaPct=SpiroideaCount/number_samples*100) directParasites$Other<-rowSums(directParasites[,6:9]) directParasites$OtherPct<-directParasites$Other/directParasites$number_samples*100 studyAreasDirect<-c("MNP","DMPPF","RMNP","GBR","YNP","PNALM") IdxAreas<-grep("^Study.Area$", colnames(directParasites)) IdxCapi<-grep("^CapillariaPct$", colnames(directParasites)) IdxToxa<-grep("^ToxascarisPct$", colnames(directParasites)) IdxUnc<-grep("^UncinariaCountPct$", colnames(directParasites)) IdxOther<-grep("^OtherPct$", colnames(directParasites)) directParasitesRed<-directParasites[c(3,1,5,2,6,4),c(IdxAreas,IdxCapi,IdxToxa,IdxUnc,IdxOther)] directParasitesMat<-t(data.matrix(directParasitesRed)) directParasitesMat2<-directParasitesMat[2:5,] # add row of NAs for spacing directParasitesMat3=rbind(NA,directParasitesMat2) directParasitesMat4<-directParasitesMat3 cols<-rep(c("white","white"," light grey","dark grey","black"),ncol(directParasitesMat3)) # create barplot and store returned value in 'a' require(plotrix) par(bty='n') a = gap.barplot(as.matrix(directParasitesMat4), gap=c(25,65), ytics=c(0,5,10,15,20,25,30,70,75,80), ylab="Percent of Faeces", las=1, xlab="", col=cols, main="Direct", ylim=c(0,42), yaxs="i", xaxt='n') # disable the default x-axis # calculate mean x-position for each group, omitting the first row # first row (NAs) is only there for spacing between groups aa = matrix(a, nrow=nrow(directParasitesMat4)) xticks = colMeans(aa[2:nrow(directParasitesMat4),]) # add axis labels at mean position axis(1, at=xticks, line=-0.6, labels=studyAreasDirect,cex=0.8,cex.lab=0.8,cex.axis=0.8,las=1,tick=FALSE,lwd=0,padj=0) legend("topleft",c("Capillaria","Toxascaris","Uncinaria","Other"),fill=c("white"," light grey","dark grey","black"),title="Parasite",bty="n",cex=0.8) mtext("A",2,line=1.75,at=43,las=1,cex=1.2) mtext("Study Area",1,line=1,las=1,cex=1) indirectParasites<-group_by(finaldb, Study.Area) %>% summarise(number_samples=length(Study.Area),TaeniidCount=sum(Taeniid_binomial), AlariaCount=sum(Alaria_binomial),SarcoCount=sum(Sarcosystis_binomial),DiphCount=sum(Diphyllobothrium_binomial),MetorchisCount=sum(Metorchis_binomial),PhysalopteraCount=sum(Physaloptera_binomial)) %>% mutate(TaeniidPct=TaeniidCount/number_samples*100,AlariaPct=AlariaCount/number_samples*100,SarcoPct=SarcoCount/number_samples*100,DiphPct=DiphCount/number_samples*100,MetorchisPct=MetorchisCount/number_samples*100,PhysalopteraPct=PhysalopteraCount/number_samples*100) studyAreasIndirect<-c("YNP","GBR","RMNP","DMPPF","MNP","PNALM") indirectParasitesRed<-indirectParasites[c(6,2,5,1,3,4),c(1,10,8,9,11)] indirectParasitesMat<-t(data.matrix(indirectParasitesRed)) indirectParasitesMat2<-indirectParasitesMat[2:5,] barpos<-barplot(indirectParasitesMat2, beside=TRUE, horiz=F,xlab="",xaxt="n",ylab="Percent of Faeces", main="Indirect",ylim=c(0,83),col=c("white","light grey","dark grey","black")) axis(1,at=c(3,8,13,18,23,28),line=-0.6, labels=studyAreasIndirect,cex=0.8,cex.lab=0.8,cex.axis=0.8,las=1,tick=FALSE,lwd=0,padj=0) legend("topright",c("Sarcocystis","Taeniid","Alaria","Diphyllobothrium"),fill=c("white"," light grey","dark grey","black"),title="Parasite",bty="n",cex=0.8) mtext("B",2,line=1.75,at=85,las=1,cex=1.2) mtext("Study Area",1,line=1,las=1,cex=1) dev.off() ############################################################################################################################################# ############################################################################################################################################# # Figure 4: Giardia prevalence by study area # Giardia was only examined in the Canadaian studies (GBR, DMNP, RMP) # Giardia should be coded as 0 or 1 but there are NAs instead of zeros in DMNP and RMNP Canada<-subset(finaldb,finaldb$Country=="Canada") Canada$Study.Area2<-factor(Canada$Study.Area2) Canada$Giardia_binomial[is.na(Canada$Giardia_binomial)]<-0 pdf("Fig4_Giardia prevalence by study area_Coast mainland island and prairies.pdf", width = 3.5, height = 3.5,pointsize = 10,bg = "white",useDingbats=F) par(mfrow=c(1,1),mgp=c(1.5,0.5,0),oma=c(0,0,0,0),mar=c(2.5,2.5,1,0.8),bty="l",las=1) title<-"Giardia" IT <-table(Canada$Giardia_binomial, Canada$Study.Area2) ITmat<-mat.or.vec(nrow(IT),ncol(IT)) nVec <- mat.or.vec(ncol(IT),1) for (i in 1:ncol(IT)) { n <- sum(IT[,i]) ITmat[,i] = IT[,i] / n * 100 nVec[i] <- n } ITmat2<-ITmat[,c(3,2,4,1)] nVec2<-nVec[c(3,2,4,1)] IT2<-IT[,c(3,2,4,1)] barpos<-barplot(ITmat2[2,], beside=F, horiz=F,xlab="Study Area",ylab="Percent of Faeces",ylim=c(0,(max(ITmat2[2,]*1.2))),main=title) text(barpos,(ITmat2[2,]+2),nVec2,cex=0.8) xlab<-c("GBRM","GBRI","RMNP","DMPPF") axis(1,at=barpos, labels=xlab,cex=0.8,cex.lab=0.8,cex.axis=0.8,las=1,tick=FALSE,lwd=0,padj=0) dev.off() ############################################################################################################################################# ############################################################################################################################################# # Models # Run using R version 3.3.3 # Add predictor variables for each study area # results were very similar when the continuous predictors were included as categorical # AIC differed by about 1, which means that it doesn't make much difference (can just use the continuous variable) finaldb$wolfDensity<-mat.or.vec(nrow(finaldb),1) finaldb$DietDiversity<-mat.or.vec(nrow(finaldb),1) # Wolf density scaled to mean of 0 and 1 vars<-c(50,11.5,50,25,25,33) scaled<-scale(vars) finaldb$wolfDensity[finaldb$Study.Area=="PNALM"]<-scaled[1] finaldb$wolfDensity[finaldb$Study.Area=="MNP"]<-scaled[2] finaldb$wolfDensity[finaldb$Study.Area=="YNP"]<-scaled[3] finaldb$wolfDensity[finaldb$Study.Area=="RMNP"]<-scaled[4] finaldb$wolfDensity[finaldb$Study.Area=="DMPF"]<-scaled[5] finaldb$wolfDensity[finaldb$Study.Area=="GBR"]<-scaled[6] # Diet diversity scaled to mean of 0 and sd of 1 vars<-c(1.78,1.71,0.14,0.57,0.57,0.33) scaled<-scale(vars) finaldb$dietDiversity[finaldb$Study.Area=="PNALM"]<-scaled[1] finaldb$dietDiversity[finaldb$Study.Area=="MNP"]<-scaled[2] finaldb$dietDiversity[finaldb$Study.Area=="YNP"]<-scaled[3] finaldb$dietDiversity[finaldb$Study.Area=="RMNP"]<-scaled[4] finaldb$dietDiversity[finaldb$Study.Area=="DMPF"]<-scaled[5] finaldb$dietDiversity[finaldb$Study.Area=="GBR"]<-scaled[6] # Recolonization hitory (NR = Not recent, R = Recent) finaldb$history[finaldb$Study.Area=="PNALM"]<-"NR" finaldb$history[finaldb$Study.Area=="MNP"]<-"R" finaldb$history[finaldb$Study.Area=="YNP"]<-"R" finaldb$history[finaldb$Study.Area=="RMNP"]<-"NR" finaldb$history[finaldb$Study.Area=="DMPF"]<-"NR" finaldb$history[finaldb$Study.Area=="GBR"]<-"NR" # Number of alternative hosts in each study area finaldb$altHosts[finaldb$Study.Area=="PNALM"]<-2 finaldb$altHosts[finaldb$Study.Area=="MNP"]<-2 finaldb$altHosts[finaldb$Study.Area=="YNP"]<-3 finaldb$altHosts[finaldb$Study.Area=="RMNP"]<-3 finaldb$altHosts[finaldb$Study.Area=="DMPF"]<-3 finaldb$altHosts[finaldb$Study.Area=="GBR"]<-0 finaldb$altHosts.cs<-scale(finaldb$altHosts) # Models for parasites with direct lifecycles direct0<-glmmadmb(DirectParasites ~ 1 + (1 | Study.Area), data = finaldb, family = "binomial", link="logit") direct1<-glmmadmb(DirectParasites ~ dietDiversity + (1 | Study.Area), data = finaldb, family = "binomial", link="logit") direct2<-glmmadmb(DirectParasites ~ wolfDensity + (1 | Study.Area), data = finaldb, family = "binomial", link="logit") direct3<-glmmadmb(DirectParasites ~ history + (1 | Study.Area), data = finaldb, family = "binomial", link="logit") direct4<-glmmadmb(DirectParasites ~ altHosts.cs + (1 | Study.Area), data = finaldb, family = "binomial", link="logit") # get confidence intervals confint(direct1, 'dietDiversity', level=0.95) confint(direct2, 'wolfDensity', level=0.95) confint(direct3, 'historyR', level=0.95) confint(direct4, 'altHosts.cs', level=0.95) Indirect0<-glmmadmb(IndirectParasites ~ 1 + (1 | Study.Area), data = finaldb, family = "binomial", link="logit") Indirect1<-glmmadmb(IndirectParasites ~ dietDiversity + (1 | Study.Area), data = finaldb, family = "binomial", link="logit") Indirect2<-glmmadmb(IndirectParasites ~ wolfDensity + (1 | Study.Area), data = finaldb, family = "binomial", link="logit") Indirect3<-glmmadmb(IndirectParasites ~ history + (1 | Study.Area), data = finaldb, family = "binomial", link="logit") Indirect4<-glmmadmb(IndirectParasites ~ altHosts.cs + (1 | Study.Area), data = finaldb, family = "binomial", link="logit") confint(Indirect1, 'dietDiversity', level=0.95) confint(Indirect2, 'wolfDensity', level=0.95) confint(Indirect3, 'historyR', level=0.95) confint(Indirect4, 'altHosts.cs', level=0.95) #Model Results for Direct Parasites model.list=list(direct0,direct1,direct2,direct3,direct4) model.names = list("null","DD","WD","H","AH") ms.mumin <- model.sel(model.list) AIC.table.full = data.frame(Model=unlist(model.names[as.numeric(rownames(ms.mumin))]),df=ms.mumin$df, logLik=round(ms.mumin$logLik,2), AICc=round(ms.mumin$AICc,2),delta.AICc=round(ms.mumin$delta,2), weight=round(ms.mumin$weight,2)) AIC.table.full # Indirect Parasites model.list=list(Indirect0,Indirect1,Indirect2,Indirect3,Indirect4) model.names = list("null","DD","WD","H","AH") ms.mumin <- model.sel(model.list) AIC.table.full = data.frame(Model=unlist(model.names[as.numeric(rownames(ms.mumin))]),df=ms.mumin$df, logLik=round(ms.mumin$logLik,2), AICc=round(ms.mumin$AICc,2),delta.AICc=round(ms.mumin$delta,2), weight=round(ms.mumin$weight,2)) AIC.table.full ############################################################################################################################################################################ ############################################################################################################################################################################ # Analyses on Corrected dataset (See Methods Section and Appendix 1) # Note: Model outputs will vary slightly on different runs as positive samples are either randomly dropped or added to the dataset # Preparing the corrected data set # Sarcosystis--subtract 13.6% from flotation study areas # Generate new columns corrected for detection problems finaldb$Sarcosystis_binomial_corrected<-finaldb$Sarcosystis_binomial ## Select the rows in the three 'flotation study areas' equal to 1 rows_with_1 = (1:length(finaldb$Sarcosystis_binomial))[(finaldb$Sarcosystis_binomial==1) & ((finaldb$Study.Area=="GBR")|(finaldb$Study.Area=="DMPF")|(finaldb$Study.Area=="RMNP"))] ## Randomly select a % of these rows and set equal to zero ## Warning: there will likely be some rounding here finaldb$Sarcosystis_binomial_corrected[sample(rows_with_1, length(rows_with_1)*0.136)] = 0 # Taeniid--add 17% to flotation study areas finaldb$Taeniid_binomial_corrected<-finaldb$Taeniid_binomial ##Select the rows in the three 'flotation study areas' equal to 0 rows_with_0 = (1:length(finaldb$Taeniid_binomial))[(finaldb$Taeniid_binomial==0) & ((finaldb$Study.Area=="GBR")|(finaldb$Study.Area=="DMPF")|(finaldb$Study.Area=="RMNP"))] # Calculate the multiplier nPostivites = nrow(finaldb[(finaldb$Taeniid_binomial==1) & ((finaldb$Study.Area=="GBR")|(finaldb$Study.Area=="DMPF")|(finaldb$Study.Area=="RMNP")),]) nSamples <- nrow(finaldb[(finaldb$Study.Area=="GBR")|(finaldb$Study.Area=="DMPF")|(finaldb$Study.Area=="RMNP"),]) pct<-nPostivites/nSamples*.17 ##Randomly select a % of these rows and set equal to 1 ##Warning: there will likely be some rounding here finaldb$Taeniid_binomial_corrected[sample(rows_with_0, length(rows_with_0)*pct)] = 1 # Diphyllobothrium--add 17% to flotation study areas # Assume only for GBR, as it wasn't detected at all in other study areas finaldb$Diphyllobothrium_binomial_corrected<-finaldb$Diphyllobothrium_binomial ##Select the rows in the three 'flotation study areas' equal to 0 rows_with_0 = (1:length(finaldb$Diphyllobothrium_binomial))[(finaldb$Diphyllobothrium_binomial==0)&(finaldb$Study.Area=="GBR")] # Calculate the multiplier nPostivites = nrow(finaldb[(finaldb$Diphyllobothrium_binomial==1) & (finaldb$Study.Area=="GBR"),]) nSamples <- nrow(finaldb[(finaldb$Study.Area=="GBR"),]) pct<-nPostivites/nSamples*.17 ##Randomly select a % of these rows and set equal to 1 ##Warning: there will likely be some rounding here finaldb$Diphyllobothrium_binomial_corrected[sample(rows_with_0, length(rows_with_0)*pct)] = 1 # Alaria--add 18.6% to flotation study areas # Assume only for DMPF and RMNP, as it wasn't detected in GBR finaldb$Alaria_binomial_corrected<-finaldb$Alaria_binomial ##Select the rows in the three 'flotation study areas' equal to 0 rows_with_0 = (1:length(finaldb$Alaria_binomial))[(finaldb$Alaria_binomial==0) & ((finaldb$Study.Area=="DMPF")|(finaldb$Study.Area=="RMNP"))] # Calculate the multiplier nPostivites = nrow(finaldb[(finaldb$Alaria_binomial==1) & ((finaldb$Study.Area=="DMPF")|(finaldb$Study.Area=="RMNP")),]) nSamples <- nrow(finaldb[(finaldb$Study.Area=="DMPF")|(finaldb$Study.Area=="RMNP"),]) pct<-nPostivites/nSamples*.186 ##Randomly select a % of these rows and set equal to 1 ##Warning: there will likely be some rounding here finaldb$Alaria_binomial_corrected[sample(rows_with_0, length(rows_with_0)*pct)] = 1 # Metorchis--add 18.6% to flotation study areas # Assume only for GBR, as it wasn't detected at all in other study areas finaldb$Metorchis_binomial_corrected<-finaldb$Metorchis_binomial ##Select the rows in the three 'flotation study areas' equal to 0 rows_with_0 = (1:length(finaldb$Metorchis_binomial))[(finaldb$Metorchis_binomial==0)] # Calculate the multiplier nPostivites = nrow(finaldb[(finaldb$Metorchis_binomial==1) & ((finaldb$Study.Area=="GBR")),]) nSamples <- nrow(finaldb[(finaldb$Study.Area=="GBR"),]) pct<-nPostivites/nSamples*.186 ##Randomly select a % of these rows and set equal to 1 ##Warning: there will likely be some rounding here finaldb$Metorchis_binomial_corrected[sample(rows_with_0, length(rows_with_0)*pct)] = 1 # Physaloptera was only detected in 2 samples from PNALM, so no adjustment made to other study areas # Direct parasite cols # Coccidia--subtract 13.6% from GBR (only found there) # Generate new columns corrected for detection problems finaldb$Coccidia_binomial_corrected<-finaldb$Coccidia_binomial ##Select the rows in the three 'flotation study areas' equal to 1 rows_with_1 = (1:length(finaldb$Coccidia_binomial))[(finaldb$Coccidia_binomial==1) & (finaldb$Study.Area=="GBR")] ##Randomly select a % of these rows and set equal to zero ##Warning: there will likely be some rounding here finaldb$Coccidia_binomial_corrected[sample(rows_with_1, length(rows_with_1)*0.136)] = 0 #Toxascaris--subtract 13.6% from flotation study areas # Generate new columns corrected for detection problems finaldb$Toxascaris_binomial_corrected<-finaldb$Toxascaris_binomial ##Select the rows in the three 'flotation study areas' equal to 1 rows_with_1 = (1:length(finaldb$Toxascaris_binomial))[(finaldb$Toxascaris_binomial==1) & ((finaldb$Study.Area=="GBR")|(finaldb$Study.Area=="DMPF")|(finaldb$Study.Area=="RMNP"))] ##Randomly select a % of these rows and set equal to zero ##Warning: there will likely be some rounding here finaldb$Toxascaris_binomial_corrected[sample(rows_with_1, length(rows_with_1)*0.136)] = 0 #Toxocara # Generate new columns corrected for detection problems finaldb$Toxocara.canis_binomial_corrected<-finaldb$Toxocara.canis_binomial ##Select the rows in the three 'flotation study areas' equal to 1 rows_with_1 = (1:length(finaldb$Toxocara.canis_binomial))[(finaldb$Toxocara.canis_binomial==1) & ((finaldb$Study.Area=="GBR")|(finaldb$Study.Area=="DMPF")|(finaldb$Study.Area=="RMNP"))] ##Randomly select a % of these rows and set equal to zero ##Warning: there will likely be some rounding here finaldb$Toxocara.canis_binomial_corrected[sample(rows_with_1, length(rows_with_1)*0.136)] = 0 # Capillaria # Generate new columns corrected for detection problems finaldb$Capillaria_binomial_corrected<-finaldb$Capillaria_binomial ##Select the rows in the three 'flotation study areas' equal to 1 rows_with_1 = (1:length(finaldb$Capillaria_binomial))[(finaldb$Capillaria_binomial==1) & ((finaldb$Study.Area=="GBR")|(finaldb$Study.Area=="DMPF")|(finaldb$Study.Area=="RMNP"))] ##Randomly select a % of these rows and set equal to zero ##Warning: there will likely be some rounding here finaldb$Capillaria_binomial_corrected[sample(rows_with_1, length(rows_with_1)*0.136)] = 0 # Spiroidea--only in 2 samples in GBR, so no change in corrected version # Generate new columns corrected for detection problems finaldb$Spiroidea_binomial_corrected<-finaldb$Spiroidea_binomial ##Select the rows in the three 'flotation study areas' equal to 1 rows_with_1 = (1:length(finaldb$Spiroidea_binomial))[(finaldb$Spiroidea_binomial==1) & ((finaldb$Study.Area=="GBR")|(finaldb$Study.Area=="DMPF")|(finaldb$Study.Area=="RMNP"))] ##Randomly select a % of these rows and set equal to zero ##Warning: there will likely be some rounding here finaldb$Spiroidea_binomial_corrected[sample(rows_with_1, length(rows_with_1)*0.136)] = 0 # Uncinaria--only in PNALM and MNP, so no correction needed # Trichuris.vulpis # Generate new columns corrected for detection problems finaldb$Trichuris.vulpis_binomial_corrected<-finaldb$Trichuris.vulpis_binomial ##Select the rows in the three 'flotation study areas' equal to 1 rows_with_1 = (1:length(finaldb$Trichuris.vulpis_binomial))[(finaldb$Trichuris.vulpis_binomial==1) & ((finaldb$Study.Area=="GBR")|(finaldb$Study.Area=="DMPF")|(finaldb$Study.Area=="RMNP"))] ##Randomly select a % of these rows and set equal to zero ##Warning: there will likely be some rounding here finaldb$Trichuris.vulpis_binomial_corrected[sample(rows_with_1, length(rows_with_1)*0.136)] = 0 # re-calculate direct and indirect parasite cols using corrected variables ToxascarisCol_corrected<-grep("Toxascaris_binomial_corrected",colnames(finaldb)) TaeniidCol_corrected<-grep("Taeniid_binomial_corrected",colnames(finaldb)) AlariaCol_corrected<-grep("Alaria_binomial_corrected",colnames(finaldb)) SarcoCol_corrected<-grep("Sarcosystis_binomial_corrected",colnames(finaldb)) CapiCol_corrected<-grep("Capillaria_binomial_corrected",colnames(finaldb)) CoccidiaCol_corrected<-grep("Coccidia_binomial_corrected",colnames(finaldb)) DiphylloCol_corrected<-grep("Diphyllobothrium_binomial_corrected",colnames(finaldb)) MetorchisCol_corrected<-grep("Metorchis_binomial_corrected",colnames(finaldb)) SpiroideaCol_corrected<-grep("Spiroidea_binomial_corrected",colnames(finaldb)) ToxocaraCol_corrected<-grep("Toxocara.canis_binomial_corrected",colnames(finaldb)) TrichurisCol_corrected<-grep("Trichuris.vulpis_binomial_corrected",colnames(finaldb)) PhysalopteraCol_corrected<-grep("Physaloptera_binomial",colnames(finaldb)) #keep this UncinariaCol_corrected<-grep("Uncinaria.stenocephala_binomial",colnames(finaldb)) paracols_corrected<-c(ToxascarisCol_corrected,TaeniidCol_corrected,AlariaCol_corrected,SarcoCol_corrected,CapiCol_corrected,CoccidiaCol_corrected,DiphylloCol_corrected,MetorchisCol_corrected,ToxocaraCol_corrected,TrichurisCol_corrected,PhysalopteraCol_corrected,UncinariaCol_corrected,SpiroideaCol_corrected) directparacols_corrected<-c(ToxascarisCol_corrected,CapiCol_corrected,CoccidiaCol_corrected,ToxocaraCol_corrected,TrichurisCol_corrected,UncinariaCol_corrected,SpiroideaCol_corrected) indirectparacols_corrected<-c(TaeniidCol_corrected,AlariaCol_corrected,SarcoCol_corrected,DiphylloCol_corrected,MetorchisCol_corrected,PhysalopteraCol_corrected) finaldb[,paracols_corrected][is.na(finaldb[,paracols_corrected])] <- 0 finaldb$DirectParasites_corrected<-mat.or.vec(nrow(finaldb),1) finaldb$IndirectParasites_corrected<-mat.or.vec(nrow(finaldb),1) # Leave in spiroidea # Calculate variable for direct parasites for (k in 1:nrow(finaldb)) { if (any(finaldb[k,directparacols_corrected]>0)) { finaldb$DirectParasites_corrected[k]<-1 } } #Calculate variable for indirect parasites for (k in 1:nrow(finaldb)) { if (any(finaldb[k,indirectparacols_corrected]>0)) { finaldb$IndirectParasites_corrected[k]<-1 } } ############################################################################################################## ############################################################################################################## # Barplot showing direct and indirect parasites by study area # Corrected for differences in methodology between study areas pdf("Fig2_Direct and indirect parasites by study area_Correction factor.pdf", width = 7.5, height = 3.5,pointsize = 10,bg = "white",useDingbats=F) par(mfrow=c(1,2),mgp=c(1.5,0.5,0),oma=c(0,0,0,0),mar=c(2,3,1,0.8),bty="l",las=1) IT <-table(finaldb$DirectParasites_corrected,finaldb$Study.Area) ITmat<-mat.or.vec(nrow(IT),ncol(IT)) nVec <- mat.or.vec(ncol(IT),1) for (i in 1:ncol(IT)) { n <- sum(IT[,i]) ITmat[,i] = IT[,i] / n * 100 nVec[i] <- n } #reorder by increasing wolf density ITmat2<-ITmat[,c(3,1,5,2,6,4)] ITmat2<-ITmat2 nVec2<-nVec[c(3,1,5,2,6,4)] barpos<-barplot(ITmat2[2,], beside=F, horiz=F,xlab="",ylab="Percent of Faeces", main="Direct",ylim=c(0,100),col=c("white","grey","grey","grey","black","black")) text(barpos,(ITmat2[2,]+2),nVec2,cex=0.8) xlab<-c("MNP","DMPPF","RMNP","GBR","YNP","PNALM") axis(1,at=barpos,line=-0.5, labels=xlab,cex=0.8,cex.lab=0.8,cex.axis=0.8,las=1,tick=FALSE,lwd=0,padj=0) legend("topleft",c("low 11.5","medium 25-35","high 50"),fill=c("white","grey","black"),title="Wolf Density/1000 kmsq",bty="n",cex=0.8) mtext("Study Area",1,line=1,las=1,cex=1) mtext("A",2,line=2,at=103,las=1,cex=1.2) IT <-table(finaldb$IndirectParasites_corrected,finaldb$Study.Area) ITmat<-mat.or.vec(nrow(IT),ncol(IT)) nVec <- mat.or.vec(ncol(IT),1) for (i in 1:ncol(IT)) { n <- sum(IT[,i]) ITmat[,i] = IT[,i] / n * 100 nVec[i] <- n } ITmat2<-ITmat[,c(6,2,5,1,3,4)] nVec2<-nVec[c(6,2,5,1,3,4)] barpos<-barplot(ITmat2[2,], beside=F, horiz=F,xlab="",ylab="Percent of Faeces", main="Indirect",ylim=c(0,100),col=c("white","grey","grey","grey","black","black")) text(barpos,(ITmat2[2,]+2),nVec2,cex=0.8) xlab<-c("YNP","GBR","RMNP","DMPPF","MNP","PNALM") axis(1,at=barpos,line=-0.5,labels=xlab,cex=0.8,cex.lab=0.8,cex.axis=0.8,las=1,tick=FALSE,lwd=0,padj=0) legend("topright",c("low","medium","high"),fill=c("white","grey","black"),title="Diet Diversity",bty="n",cex=0.8) mtext("Study Area",1,line=1,las=1,cex=1) mtext("B",2,line=2,at=103,las=1,cex=1.2) dev.off() ########################################################################################################## ########################################################################################################## # Barplot showing prevalence of each of the direct parasites by study area # with Correction Factor applied pdf("Fig3_Direct and indirect parasites separate by study area_CorrectionFactors.pdf", width = 7.5, height = 3.5,pointsize = 10,bg = "white",useDingbats=F) #5.5 inches is 140mm or 1.5 columns par(mfrow=c(1,2),bty="n",las=1,mgp=c(1.5,0.5,0),oma=c(0,0,0,0),mar=c(2,2.5,1,0.8)) directParasites<-group_by(finaldb, Study.Area) %>% summarise(number_samples=length(Study.Area),ToxascarisCount=sum(Toxascaris_binomial_corrected), CapillariaCount=sum(Capillaria_binomial_corrected),UncinariaCount=sum(Uncinaria.stenocephala_binomial), CoccidiaCount=sum(Coccidia_binomial_corrected),ToxocaraCount=sum(Toxocara.canis_binomial_corrected),TrichurisCount=sum(Trichuris.vulpis_binomial_corrected),SpiroideaCount=sum(Spiroidea_binomial_corrected)) %>% mutate(ToxascarisPct=ToxascarisCount/number_samples*100,CapillariaPct=CapillariaCount/number_samples*100,UncinariaCountPct=UncinariaCount/number_samples*100, CoccidiaPct=CoccidiaCount/number_samples*100,ToxocaraPct=ToxocaraCount/number_samples*100,TrichurisPct=TrichurisCount/number_samples*100,SpiroideaPct=SpiroideaCount/number_samples*100) directParasites$Other<-rowSums(directParasites[,6:9]) directParasites$OtherPct<-directParasites$Other/directParasites$number_samples*100 studyAreasDirect<-c("MNP","DMPPF","RMNP","GBR","YNP","PNALM") directParasitesRed<-directParasites[c(3,1,5,2,6,4),c(1,11,10,12,18)] directParasitesMat<-t(data.matrix(directParasitesRed)) directParasitesMat2<-directParasitesMat[2:5,] # add row of NAs for spacing directParasitesMat3=rbind(NA,directParasitesMat2) cols<-rep(c("white","white"," light grey","dark grey","black"),ncol(directParasitesMat3)) # create barplot and store returned value in 'a' require(plotrix) par(bty='n') a = gap.barplot(as.matrix(directParasitesMat3), gap=c(25,65), ytics=c(0,5,10,15,20,25,30,70,75,80), ylab="Percent of Faeces", las=1, xlab="", col=cols, main="Direct", ylim=c(0,42), yaxs="i", xaxt='n') # disable the default x-axis # calculate mean x-position for each group, omitting the first row # first row (NAs) is only there for spacing between groups aa = matrix(a, nrow=nrow(directParasitesMat3)) xticks = colMeans(aa[2:nrow(directParasitesMat3),]) # add axis labels at mean position axis(1, at=xticks, line=-0.6, labels=studyAreasDirect,cex=0.8,cex.lab=0.8,cex.axis=0.8,las=1,tick=FALSE,lwd=0,padj=0) legend("topleft",c("Capillaria","Toxascaris","Uncinaria","Other"),fill=c("white"," light grey","dark grey","black"),title="Parasite",bty="n",cex=0.8) mtext("A",2,line=1.75,at=43,las=1,cex=1.2) mtext("Study Area",1,line=1,las=1,cex=1) indirectParasites<-group_by(finaldb, Study.Area) %>% summarise(number_samples=length(Study.Area),TaeniidCount=sum(Taeniid_binomial), AlariaCount=sum(Alaria_binomial),SarcoCount=sum(Sarcosystis_binomial),DiphCount=sum(Diphyllobothrium_binomial),MetorchisCount=sum(Metorchis_binomial),PhysalopteraCount=sum(Physaloptera_binomial)) %>% mutate(TaeniidPct=TaeniidCount/number_samples*100,AlariaPct=AlariaCount/number_samples*100,SarcoPct=SarcoCount/number_samples*100,DiphPct=DiphCount/number_samples*100,MetorchisPct=MetorchisCount/number_samples*100,PhysalopteraPct=PhysalopteraCount/number_samples*100) studyAreasIndirect<-c("YNP","GBR","RMNP","DMPPF","MNP","PNALM") indirectParasitesRed<-indirectParasites[c(6,2,5,1,3,4),c(1,11,9,10,12)] indirectParasitesMat<-t(data.matrix(indirectParasitesRed)) indirectParasitesMat2<-indirectParasitesMat[2:5,] # Sarco is over-reported by flotation by 13.6% (Wolf et al) indirectParasitesMat2[1,2:4]<-indirectParasitesMat2[1,2:4]-.136*indirectParasitesMat2[1,2:4] # Taeniids are under-reported by 17% by flotation (Oge et al) indirectParasitesMat2[2,2:4]<-indirectParasitesMat2[2,2:4]+.17*indirectParasitesMat2[2,2:4] # Alaria (fluke egg)18.6% under-reported by flotation indirectParasitesMat2[3,2:4]<-indirectParasitesMat2[3,2:4]+.186*indirectParasitesMat2[3,2:4] # Diphyllobothrium--assume same under-reporting as Taeniids indirectParasitesMat2[4,2:4]<-indirectParasitesMat2[4,2:4]+.17*indirectParasitesMat2[4,2:4] barpos<-barplot(indirectParasitesMat2, beside=TRUE, horiz=F,xlab="",xaxt="n",ylab="Percent of Faeces", main="Indirect",ylim=c(0,83),col=c("white","light grey","dark grey","black")) axis(1,at=c(3,8,13,18,23,28),line=-0.6, labels=studyAreasIndirect,cex=0.8,cex.lab=0.8,cex.axis=0.8,las=1,tick=FALSE,lwd=0,padj=0) legend("topright",c("Sarcocystis","Taeniid","Alaria","Diphyllobothrium"),fill=c("white"," light grey","dark grey","black"),title="Parasite",bty="n",cex=0.8) mtext("B",2,line=1.75,at=85,las=1,cex=1.2) mtext("Study Area",1,line=1,las=1,cex=1) dev.off() ########################################################################################################## ########################################################################################################## # Models on Corrected data # Note that parameter estimates will vary slightly on each run due to the way positive samples were randomly added or removed finaldb$wolfDensity<-mat.or.vec(nrow(finaldb),1) finaldb$DietDiversity<-mat.or.vec(nrow(finaldb),1) # Wolf density scaled to mean of 0 and 1 vars<-c(50,11.5,50,25,25,33) scaled<-scale(vars) finaldb$wolfDensity[finaldb$Study.Area=="PNALM"]<-scaled[1] finaldb$wolfDensity[finaldb$Study.Area=="MNP"]<-scaled[2] finaldb$wolfDensity[finaldb$Study.Area=="YNP"]<-scaled[3] finaldb$wolfDensity[finaldb$Study.Area=="RMNP"]<-scaled[4] finaldb$wolfDensity[finaldb$Study.Area=="DMPF"]<-scaled[5] finaldb$wolfDensity[finaldb$Study.Area=="GBR"]<-scaled[6] # Diet diversity scaled to mean of 0 and sd of 1 vars<-c(1.78,1.71,0.14,0.57,0.57,0.33) scaled<-scale(vars) finaldb$dietDiversity[finaldb$Study.Area=="PNALM"]<-scaled[1] finaldb$dietDiversity[finaldb$Study.Area=="MNP"]<-scaled[2] finaldb$dietDiversity[finaldb$Study.Area=="YNP"]<-scaled[3] finaldb$dietDiversity[finaldb$Study.Area=="RMNP"]<-scaled[4] finaldb$dietDiversity[finaldb$Study.Area=="DMPF"]<-scaled[5] finaldb$dietDiversity[finaldb$Study.Area=="GBR"]<-scaled[6] finaldb$history[finaldb$Study.Area=="PNALM"]<-"NR" finaldb$history[finaldb$Study.Area=="MNP"]<-"R" finaldb$history[finaldb$Study.Area=="YNP"]<-"R" finaldb$history[finaldb$Study.Area=="RMNP"]<-"NR" finaldb$history[finaldb$Study.Area=="DMPF"]<-"NR" finaldb$history[finaldb$Study.Area=="GBR"]<-"NR" finaldb$altHosts[finaldb$Study.Area=="PNALM"]<-2 finaldb$altHosts[finaldb$Study.Area=="MNP"]<-2 finaldb$altHosts[finaldb$Study.Area=="YNP"]<-3 finaldb$altHosts[finaldb$Study.Area=="RMNP"]<-3 finaldb$altHosts[finaldb$Study.Area=="DMPF"]<-3 finaldb$altHosts[finaldb$Study.Area=="GBR"]<-0 finaldb$altHosts.cs<-scale(finaldb$altHosts) direct0<-glmmadmb(DirectParasites_corrected ~ 1 + (1 | Study.Area), data = finaldb, family = "binomial", link="logit") direct1<-glmmadmb(DirectParasites_corrected ~ dietDiversity + (1 | Study.Area), data = finaldb, family = "binomial", link="logit") direct2<-glmmadmb(DirectParasites_corrected ~ wolfDensity + (1 | Study.Area), data = finaldb, family = "binomial", link="logit") direct3<-glmmadmb(DirectParasites_corrected ~ history + (1 | Study.Area), data = finaldb, family = "binomial", link="logit") direct4<-glmmadmb(DirectParasites_corrected ~ altHosts.cs + (1 | Study.Area), data = finaldb, family = "binomial", link="logit") # get confidence intervals confint(direct1, 'dietDiversity', level=0.95) confint(direct2, 'wolfDensity', level=0.95) confint(direct3, 'historyR', level=0.95) confint(direct4, 'altHosts.cs', level=0.95) Indirect0<-glmmadmb(IndirectParasites_corrected ~ 1 + (1 | Study.Area), data = finaldb, family = "binomial", link="logit") Indirect1<-glmmadmb(IndirectParasites_corrected ~ dietDiversity + (1 | Study.Area), data = finaldb, family = "binomial", link="logit") Indirect2<-glmmadmb(IndirectParasites_corrected ~ wolfDensity + (1 | Study.Area), data = finaldb, family = "binomial", link="logit") Indirect3<-glmmadmb(IndirectParasites_corrected ~ history + (1 | Study.Area), data = finaldb, family = "binomial", link="logit") Indirect4<-glmmadmb(IndirectParasites_corrected ~ altHosts.cs + (1 | Study.Area), data = finaldb, family = "binomial", link="logit") # get confidence intervals confint(Indirect1, 'dietDiversity', level=0.95) confint(Indirect2, 'wolfDensity', level=0.95) confint(Indirect3, 'historyR', level=0.95) confint(Indirect4, 'altHosts.cs', level=0.95) #Model Results for Direct Parasites model.list=list(direct0,direct1,direct2,direct3,direct4) model.names = list("null","DD","WD","H","AH") ms.mumin <- model.sel(model.list) AIC.table.full = data.frame(Model=unlist(model.names[as.numeric(rownames(ms.mumin))]),df=ms.mumin$df, logLik=round(ms.mumin$logLik,2), AICc=round(ms.mumin$AICc,2),delta.AICc=round(ms.mumin$delta,2), weight=round(ms.mumin$weight,2)) AIC.table.full # Indirect Parasites model.list=list(Indirect0,Indirect1,Indirect2,Indirect3,Indirect4) model.names = list("null","DD","WD","H","AH") ms.mumin <- model.sel(model.list) AIC.table.full = data.frame(Model=unlist(model.names[as.numeric(rownames(ms.mumin))]),df=ms.mumin$df, logLik=round(ms.mumin$logLik,2), AICc=round(ms.mumin$AICc,2),delta.AICc=round(ms.mumin$delta,2), weight=round(ms.mumin$weight,2)) AIC.table.full