library(ggplot2) # For plotting functions library(gridExtra) # For making compound graph SelGrad=T # Set to true to run analysis on selection gradients instead of differentials ## Import data SData<-read.csv("Selection_Data.csv",header=T) SData$AuthorYear<-paste(SData$Author,SData$Year) ## Data Processing SData$s<-abs(SData$Selection.differential..s.) # Make selction differentials absolute value SData$B<-abs(SData$Selection.gradient..B.) # Make selction gradients absolute value SData$Region<-gsub("yes","N",SData$Native) SData$Region<-gsub("no","I",SData$Region) SData$Species<-gsub(" \\(.*","",SData$Species) unique(SData$Author) # Unique authors unique(SData$AuthorYear) sum(!is.na(SData$s[SData$Region=="N"])) # Number of selection differentials sum(!is.na(SData$s[SData$Region=="I"])) # Number of selection differentials sum(!is.na(SData$B[SData$Region=="N"])) # Number of selection gradients sum(!is.na(SData$B[SData$Region=="I"])) # Number of selection gradients # Check for typos & inconsistencies in key variables: unique(SData$Region) hist(SData$s) SelData<-data.frame(SData[,c("Species","Region","Trait","s","B")]) # TRUE to run script for Selection Gradients; FALSE to run script for Selection Differentials if (SelGrad==TRUE){SelData$s<-SelData$B} ## Bootstrap estimates # Input data # Selection Differentials NatDat<-SelData[SelData$Region=="N" & complete.cases(SelData$s),] # Data for Native species, excluding NA IntDat<-SelData[SelData$Region=="I" & complete.cases(SelData$s),] # Data for Introduced species, excluding NA #ptm <- proc.time() # Check processing time - useful to figure out run-time for # Bootstrap Iterations # Bootstrap setup IterN<-10000 # Number of bootstrap iterations NatSims<-{} IntSims<-{} # Begin bootstrap for (i in 1:IterN){ # Native Species NatBootData<-{} for(NSp in unique(NatDat$Species)){ NatTrList<-unique(paste(NatDat$Trait[NatDat$Species==NSp])) NatTr<-sample(NatTrList,length(NatTrList),replace=T) for(NTr in NatTr){ NatSel<-NatDat$s[NatDat$Species==NSp & NatDat$Trait==NTr] NatSel<-mean(sample(NatSel,length(NatSel))) # IF > 1 measurement, bootstrap and take mean NatBootData<-rbind(NatBootData,NatSel) } } # Calculate mean s for bootstrap iteration NatSims<-c(NatSims,mean(NatBootData)) # Introduced Species IntBootData<-{} for(NSp in unique(IntDat$Species)){ IntTrList<-unique(paste(IntDat$Trait[IntDat$Species==NSp])) IntTr<-sample(IntTrList,length(IntTrList),replace=T) for(NTr in IntTr){ IntSel<-IntDat$s[IntDat$Species==NSp & IntDat$Trait==NTr] IntSel<-mean(sample(IntSel,length(IntSel))) # IF > 1 measurement, bootstrap and take mean IntBootData<-rbind(IntBootData,IntSel) } } # Calculate mean s for bootstrap iteration IntSims<-c(IntSims,mean(IntBootData)) } # proc.time() - ptm # Sort bootstrap outputs from low to high NatSims<-sort(NatSims) IntSims<-sort(IntSims) # Make into data frame for graphing HistData<-data.frame(s=SelData$s,Region=SelData$Region) # Extract 95% CIs data from SimData Nat95<-sort(NatSims)[round(IterN*0.025,0):round(IterN*0.975,0)] Int95<-sort(IntSims)[round(IterN*0.025,0):round(IterN*0.975,0)] # Combine SimData SimData95<-rbind(data.frame(s=Nat95,Region="N"),data.frame(s=Int95,Region="I")) # Calculate 95% CIs from sims CIs<-c(sort(NatSims)[round(IterN*0.025,0)], # Native, lower 2.5% sort(NatSims)[round(IterN*0.975,0)], # Native, upper 97.5% sort(IntSims)[round(IterN*0.025,0)], # Intro, lower 2.5% sort(IntSims)[round(IterN*0.975,0)]) # Intro, upper 97.5% ###### GRAPHING theme_format<-function(base_size = 24, base_family = ""){ theme_grey(base_size = base_size, base_family = base_family) %+replace% theme(axis.text = element_text(size = rel(0.8),colour="black"), axis.ticks = element_line(colour = "black"), legend.key = element_rect(colour = NA), legend.position = "none", panel.border = element_rect(fill = NA, colour = NA), panel.grid.major = element_line(colour = NA,size = 0), panel.grid.minor = element_line(colour = NA,size = 0), panel.background = element_rect(fill=NA), axis.line = element_line(colour ="black") ) } ## Plot Histograms + Bootstrap mean + CI on same graph ciH<-0.01 # Height (i.e. thickness) of SE bars p <- ggplot() p <- p + geom_freqpoly(data=HistData[HistData$Region=="N",], aes(s,y=(..count..)/sum(..count..)),alpha = 0.6,colour="#1fcedd",size=2) p <- p + geom_freqpoly(data=HistData[HistData$Region=="I",], aes(s,y=(..count..)/sum(..count..)),alpha = 0.5,colour="#f53751",size=2) p <- p + geom_line(aes(x=mean(NatSims),y=c(0,ciH)),colour="#15897d",size=1) p <- p + geom_rect(aes(xmin=CIs[1],xmax=CIs[2],ymin=0,ymax=ciH),colour="white",fill="#1fcedd88") p <- p + geom_line(aes(x=mean(IntSims),y=c(0,ciH)),colour="#f53751",size=1) p <- p + geom_rect(aes(xmin=CIs[3],xmax=CIs[4],ymin=0,ymax=ciH),colour="white",fill="#f5375188") p <- p + theme_format() + ylab("Frequency") + scale_x_continuous(limits = c(0, 1.5)) p # Open .svg file to save figure if(SelGrad==F){svg("SelectionDiffWithinSp.svg")} else{ svg("SelectionGradWithinSp.svg")} p dev.off() # close file # Some summary info paste ("Means: Native =",round(mean(NatSims),3)," Introduced =",round(mean(IntSims),3)) paste ("95%CI: Native =",round(CIs[1],3),"-",round(CIs[2],3)) paste ("95%CI: Introduced =",round(CIs[3],3),"-",round(CIs[4],3)) paste ("(I-N)/N = ",round((mean(IntSims)-mean(NatSims))/mean(NatSims)*100,1),"%") paste ("N: Native s =",length(!is.na(NatDat$s)),"Intro s = ",length(!is.na(IntDat$s)))