##Appendix 6: R-code #Continental Scale then regional scale #Load Libraries library(vegan) library(boot) library(car) #Read in data Laurentia<-read.csv("Ord_occs.csv",header=TRUE,sep=",",as.is=TRUE) provremod<-read.csv("provremod.csv",header=TRUE,sep=",",as.is=TRUE) col_age<-read.csv("col_age.csv",header=TRUE,sep=",",as.is=TRUE) outgroup<-read.csv("UK-culled.csv",header=T,sep=",") bretsky<-read.csv("BretskyAppCollCull.csv",header=TRUE,sep=",",row.names=1) app2013<-read.csv("Appalachian_2013.csv",header=TRUE,sep=",",row.names=1) AbundanceMatrix<-read.csv("AbundanceMatrix.csv",header=TRUE,sep=",",row.names=1) paradist<-read.csv("paradist.csv") ###############CONTINENTAL SCALE################################ #prep data #check rows, should be 0 which(Laurentia$collection_no!=provremod$collection) #put province data in occurrence table Laurentia<-cbind(Laurentia,provremod) #Clean up the data #Nix the Lingula (but not Lingulasma) Laurentia$occurrence.genus_name<-gsub("Lingula$","Lingulid",Laurentia$occurrence.genus_name) #Nix the blank environments Laurentia<-Laurentia[which(Laurentia$environment!=""),] #Oust Dalve Laurentia<-Laurentia[which(Laurentia$collection_no!=389&Laurentia$collection_no!=391&Laurentia$collection_no!=394),] #Nix Schizambon Laurentia<-Laurentia[which(Laurentia$occurrence.genus_name!="Schizambon"),] #check collections still match (should = 0) which(Laurentia$collection_no != col_age$collection_no) #Put ages in now, this table was built off the culled Laurentia Laurentia<-cbind(Laurentia,col_age) #Correct Misspellings Laurentia$occurrence.genus_name<-gsub("Apheorthis","Apheoorthis",Laurentia$occurrence.genus_name) Laurentia$occurrence.genus_name<-gsub("Cerurinella","Ceraurinella",Laurentia$occurrence.genus_name) Laurentia$occurrence.genus_name<-gsub("Cymatonata","Cymatonota",Laurentia$occurrence.genus_name) Laurentia$occurrence.genus_name<-gsub("Drepanodos","Drepanodus",Laurentia$occurrence.genus_name) Laurentia$occurrence.genus_name<-gsub("Heterotrpya","Heterotrypa",Laurentia$occurrence.genus_name) Laurentia$occurrence.genus_name<-gsub("Hormotelus","Homotelus",Laurentia$occurrence.genus_name) Laurentia$occurrence.genus_name<-gsub("Missiquoia","Missisquoia",Laurentia$occurrence.genus_name) Laurentia$occurrence.genus_name<-gsub("Paterulus","Paterula",Laurentia$occurrence.genus_name) Laurentia$occurrence.genus_name<-gsub("Pauposira","Paupospira",Laurentia$occurrence.genus_name) Laurentia$occurrence.genus_name<-gsub("Pionodema","Pianodema",Laurentia$occurrence.genus_name) #Nix Camerates Laurentia<-Laurentia[which(Laurentia$occurrence.genus_name!="Camerate"),] #Nix Skolithos Laurentia<-Laurentia[which(Laurentia$occurrence.genus_name!="Skolithos"),] #Add outgroup Laur<-Laurentia[,c(1,3,5,11,12,17,21,27,29)] Out<-outgroup[,c(1,4,6,17,18,23,30)] Prov<-replicate(length(rownames(Out)),5) Age<-replicate(length(rownames(Out)),3) Out<-cbind(Out,Prov,Age) occdata<-rbind(Laur,Out) ####### #Generate the Occurrence Matrix #make sure province and age column are called "Prov" and "Age" #If you are trying to add data to the occurrence matrix you need to put your collections in here addnames<-sort(unique(c(names(bretsky),names(app2013)))) addrows<-c(rownames(bretsky),rownames(app2013)) GenerateOccurrenceMatrix<-function(pbdbdata,addins){ if (addins=="TRUE"){ UniqueGenera<-sort(unique(c(pbdbdata$occurrence.genus_name,addnames))) UniqueCollections<-unique(c(pbdbdata$collection_no,addrows)) } else { UniqueGenera<-sort(unique(pbdbdata$occurrence.genus_name)) UniqueCollections<-unique(pbdbdata$collection_no) } Occurrences <- matrix(data=0, nrow=length(UniqueCollections), ncol=length(UniqueGenera)+2, dimnames=list(UniqueCollections,c(UniqueGenera,"Province","Age"))) for (i in 1:nrow(pbdbdata)){ CollectionNumber <- which(UniqueCollections==pbdbdata$collection_no[i]) GenusNumber <- which(UniqueGenera==pbdbdata$occurrence.genus_name[i]) Occurrences[CollectionNumber,GenusNumber] <-1 Occurrences[CollectionNumber,length(UniqueGenera)+1]<-as.numeric(pbdbdata$Prov[i]) Occurrences[CollectionNumber,length(UniqueGenera)+2]<-as.numeric(pbdbdata$Age[i]) } Occurrences<-as.data.frame(Occurrences) Occurrences } OccurrenceMatrix<-GenerateOccurrenceMatrix(occdata,TRUE) #Make sure nothing got weird (should be 0) which(rownames(OccurrenceMatrix) != unique(c(occdata$collection_no,addrows))) #Fill in added rows OccurrenceMatrix[which(rownames(OccurrenceMatrix)%in%addrows),names(OccurrenceMatrix)=="Province"]<-1 OccurrenceMatrix[which(rownames(OccurrenceMatrix)%in%addrows),names(OccurrenceMatrix)=="Age"]<-3 for (i in 1:nrow (bretsky)){ colnumber<-which(rownames(bretsky)[i]==rownames(OccurrenceMatrix)) filllist<-names(bretsky)[which(bretsky[i,]>0)] GenusNumbers<-which(names(OccurrenceMatrix)%in%filllist) OccurrenceMatrix[colnumber,GenusNumbers]<-1 } for (i in 1:nrow (app2013)){ colnumber<-which(rownames(app2013)[i]==rownames(OccurrenceMatrix)) filllist<-names(app2013)[which(app2013[i,]>0)] GenusNumbers<-which(names(OccurrenceMatrix)%in%filllist) OccurrenceMatrix[colnumber,GenusNumbers]<-1 } #Eliminate morpho-bryozoan (unique to added app collections) and "nautiloid" OccurrenceMatrix<-OccurrenceMatrix[,names(OccurrenceMatrix)!="nautiloid"] bryos<-grep("Bryozoan.",names(OccurrenceMatrix)) OccurrenceMatrix<-OccurrenceMatrix[,names(OccurrenceMatrix)[-bryos]] #Eliminate singleton collections OccurrenceMatrix<-OccurrenceMatrix[which(rowSums(OccurrenceMatrix[,c(names(OccurrenceMatrix)!="Province" & names(OccurrenceMatrix)!="Age")])>1),] ######################################################################################### ####Aggregate, continent scale, Early/Middle/Late P11<-OccurrenceMatrix[which(OccurrenceMatrix$Province==1 & OccurrenceMatrix$Age==1),] P11<-rbind(P11,colSums(P11)) AggP11<-P11[nrow(P11),] AggP11$Province<-1 AggP11$Age<-1 P21<-OccurrenceMatrix[which(OccurrenceMatrix$Province==2 & OccurrenceMatrix$Age==1),] P21<-rbind(P21,colSums(P21)) AggP21<-P21[nrow(P21),] AggP21$Province<-2 AggP21$Age<-1 if(length(which(OccurrenceMatrix$Province==3 & OccurrenceMatrix$Age==1))==0) {AggP31<-AggP21 AggP31[,which(names(AggP31)=="Province")]<-3 AggP31[1,1:(length(names(AggP31))-2)]<-replicate(length(names(AggP31))-2,0) } if(length(which(OccurrenceMatrix$Province==3 & OccurrenceMatrix$Age==1))>0){ P31<-OccurrenceMatrix[which(OccurrenceMatrix$Province==3 & OccurrenceMatrix$Age==1),] P31<-rbind(P31,colSums(P31)) AggP31<-P31[nrow(P31),] AggP31$Province<-3 AggP31$Age<-1 } P41<-OccurrenceMatrix[which(OccurrenceMatrix$Province==4 & OccurrenceMatrix$Age==1),] P41<-rbind(P41,colSums(P41)) AggP41<-P41[nrow(P41),] AggP41$Province<-4 AggP41$Age<-1 P12<-OccurrenceMatrix[which(OccurrenceMatrix$Province==1 & OccurrenceMatrix$Age==2),] P12<-rbind(P12,colSums(P12)) AggP12<-P12[nrow(P12),] AggP12$Province<-1 AggP12$Age<-2 P22<-OccurrenceMatrix[which(OccurrenceMatrix$Province==2 & OccurrenceMatrix$Age==2),] P22<-rbind(P22,colSums(P22)) AggP22<-P22[nrow(P22),] AggP22$Province<-2 AggP22$Age<-2 P32<-OccurrenceMatrix[which(OccurrenceMatrix$Province==3 & OccurrenceMatrix$Age==2),] P32<-rbind(P32,colSums(P32)) AggP32<-P32[nrow(P32),] AggP32$Province<-3 AggP32$Age<-2 P42<-OccurrenceMatrix[which(OccurrenceMatrix$Province==4 & OccurrenceMatrix$Age==2),] P42<-rbind(P42,colSums(P42)) AggP42<-P42[nrow(P42),] AggP42$Province<-4 AggP42$Age<-2 P13<-OccurrenceMatrix[which(OccurrenceMatrix$Province==1 & OccurrenceMatrix$Age==3),] P13<-rbind(P13,colSums(P13)) AggP13<-P13[nrow(P13),] AggP13$Province<-1 AggP13$Age<-3 P23<-OccurrenceMatrix[which(OccurrenceMatrix$Province==2 & OccurrenceMatrix$Age==3),] P23<-rbind(P23,colSums(P23)) AggP23<-P23[nrow(P23),] AggP23$Province<-2 AggP23$Age<-3 P33<-OccurrenceMatrix[which(OccurrenceMatrix$Province==3 & OccurrenceMatrix$Age==3),] P33<-rbind(P33,colSums(P33)) AggP33<-P33[nrow(P33),] AggP33$Province<-3 AggP33$Age<-3 P43<-OccurrenceMatrix[which(OccurrenceMatrix$Province==4 & OccurrenceMatrix$Age==3),] P43<-rbind(P43,colSums(P43)) AggP43<-P43[nrow(P43),] AggP43$Province<-4 AggP43$Age<-3 P53<-OccurrenceMatrix[which(OccurrenceMatrix$Province==5 & OccurrenceMatrix$Age==3),] P53<-rbind(P53,colSums(P53)) AggP53<-P53[nrow(P53),] AggP53$Province<-5 AggP53$Age<-3 Aggregates<-rbind(AggP11,AggP21,AggP31,AggP41,AggP12,AggP22,AggP32,AggP42,AggP13,AggP23,AggP33,AggP43,AggP53) rownames(Aggregates)<-c(11,21,31,41,12,22,32,42,13,23,33,43,53) interprovincesimilarityAgg<-function(ProvinceAn,ProvinceBn,Agen){ provinceageA<-Aggregates[Aggregates$Province==ProvinceAn & Aggregates$Age==Agen,colnames(Aggregates)!="Province"&colnames(Aggregates)!="Age"] provinceageB<-Aggregates[Aggregates$Province==ProvinceBn & Aggregates$Age==Agen, colnames(Aggregates)!="Province" & colnames(Aggregates)!="Age"] provinceage<-rbind(provinceageA,provinceageB) similarity<-1-vegdist(provinceage, method="jaccard", binary=TRUE) similarity } #Jackknife# # Aggregate Table/Jacknife interprovincesimilaritySubset<-function(ProvinceAn,ProvinceBn,Agen){ provinceageA<-subset[subset$Province==ProvinceAn & subset$Age==Agen,colnames(subset)!="Province"&colnames(subset)!="Age"] provinceageB<-subset[subset$Province==ProvinceBn & subset$Age==Agen, colnames(subset)!="Province" & colnames(subset)!="Age"] provinceage<-rbind(provinceageA,provinceageB) similarity<-1-vegdist(provinceage, method="jaccard", binary=TRUE) similarity } calculatesubSim<-function(x,y,z){ interprovincesimilaritySubset(x,y,z) #Where x is prov1, y is prov2, z is age,interprovincesimilarityAgg calls from subset of Aggregates } #Aggregate Table (includes the actual jackknife) pseudovalues <- vector(length=length(Aggregates)-2, mode="numeric") #last two columns are province/age identifiers .: length -2. AggregateTable<-as.data.frame(matrix(data=9,nrow=22,ncol=8,dimnames=list(1:22,c("ProvinceA","ProvinceB","Age","SimIndex","Jsimilarity","JStdError","JLCL","JUCL")))) row <- 1 for (Age in 1:3){ for (ProvinceA in 1:3){ for (ProvinceB in (ProvinceA+1):4){ if (ProvinceA==3 & Age==1|ProvinceB==3 & Age==1){ #no data Jsim <- 13 SimIndex<-13 } else { SimIndex<-interprovincesimilarityAgg(ProvinceA,ProvinceB,Age) theEstimate<-interprovincesimilarityAgg(ProvinceA,ProvinceB,Age) for (i in 1:(length(Aggregates)-2)){ subset<-Aggregates[-i] subsetEstimate<-calculatesubSim(ProvinceA,ProvinceB,Age) pseudovalues[i]<-theEstimate-(length(Aggregates)-3)*(subsetEstimate-theEstimate) #lengthAggregates -3 instead of -1 because the function uses Aggregates without the last two columns (province&age) so this is (Aggregates-2)-1 } jackknifeEstimate <- mean(pseudovalues) jackknifeStdError <- sd(pseudovalues)/sqrt(length(pseudovalues)) alpha <- 0.05 lowerCL <- jackknifeEstimate + qt(p=alpha/2, df=length(pseudovalues)-1) * jackknifeStdError upperCL <- jackknifeEstimate - qt(p=alpha/2, df=length(pseudovalues)-1) * jackknifeStdError } #fill in row AggregateTable[row,1] <- ProvinceA AggregateTable[row,2] <- ProvinceB AggregateTable[row,3] <- Age AggregateTable[row,4] <- SimIndex AggregateTable[row,5] <- jackknifeEstimate AggregateTable[row,6] <- jackknifeStdError AggregateTable[row,7] <- lowerCL AggregateTable[row,8] <- upperCL # increment row number row <- row+1 } } } row<-19 for (ProvinceA in 1:4){ ProvinceB<-5 Age<-3 SimIndex<-interprovincesimilarityAgg(ProvinceA,ProvinceB,Age) theEstimate<-interprovincesimilarityAgg(ProvinceA,ProvinceB,Age) for (i in 1:(length(Aggregates)-2)){ subset<-Aggregates[-i] subsetEstimate<-calculatesubSim(ProvinceA,ProvinceB,Age) pseudovalues[i]<-theEstimate-(length(Aggregates)-3)*(subsetEstimate-theEstimate) } #lengthAggregates -3 instead of -1 because the function uses Aggregates without the last two columns (province&age) so this is (Aggregates-2)-1 jackknifeEstimate <- mean(pseudovalues) jackknifeStdError <- sd(pseudovalues)/sqrt(length(pseudovalues)) alpha <- 0.05 lowerCL <- jackknifeEstimate + qt(p=alpha/2, df=length(pseudovalues)-1) * jackknifeStdError upperCL <- jackknifeEstimate - qt(p=alpha/2, df=length(pseudovalues)-1) * jackknifeStdError AggregateTable[row,1] <- ProvinceA AggregateTable[row,2] <- ProvinceB AggregateTable[row,3] <- Age AggregateTable[row,4] <- SimIndex AggregateTable[row,5] <- jackknifeEstimate AggregateTable[row,6] <- jackknifeStdError AggregateTable[row,7] <- lowerCL AggregateTable[row,8] <- upperCL # increment row number row <- row+1 } AggregateTable write.csv(AggregateTable,"AggregateSimilarity.csv") ###########################PLOTS######################### plotagr<-AggregateTable[AggregateTable$SimIndex<9 & AggregateTable$SimIndex>0,] pdf(file="agsim_Laurentia.pdf",height=7,width=14) stripchart(plotagr$SimIndex~plotagr$Age, xlim=c(0,max(plotagr$SimIndex)+.01), ylim=c(.5,3.5), xlab="Jaccard's Similarity",ylab="Age",las=1, pch=16, col="white", main="Aggregated Similarity") points(plotagr$SimIndex[plotagr$ProvinceB==5], plotagr$Age[plotagr$ProvinceB==5], pch=16, col="grey50") points(plotagr$SimIndex[plotagr$ProvinceB!=5], plotagr$Age[plotagr$ProvinceB!=5], pch=16, col="black") points(plotagr$SimIndex[plotagr$Age==1&(plotagr$ProvinceA==2|plotagr$ProvinceB==2)],plotagr$Age[plotagr$Age==1&(plotagr$ProvinceA==2|plotagr$ProvinceB==2)],pch=16,col="white") labelsage23<-c("AM","AS","AW","MS","MW","SW","AM","AS","AW","MS","MW","SW") simsage23<-plotagr$SimIndex[plotagr$ProvinceB!=5&plotagr$Age!=1] agesage23<-plotagr$Age[plotagr$ProvinceB!=5&plotagr$Age!=1]-.15 text(c(simsage23),c(agesage23),col="black",cex=.8,label=c(labelsage23)) labelsage1<-c("AW") simsage1<-plotagr$SimIndex[plotagr$Age==1&plotagr$ProvinceA==1&plotagr$ProvinceB==4] agesage1<-plotagr$Age[plotagr$Age==1&plotagr$ProvinceA==1&plotagr$ProvinceB==4]-.15 text(c(simsage1),c(agesage1),col="black",cex=.8,label=c(labelsage1)) labelsoutgroup<-c("A","M","S","W") simsoutgroup<-plotagr$SimIndex[plotagr$ProvinceB==5] agesoutgroup<-plotagr$Age[plotagr$ProvinceB==5]+.15 text(c(simsoutgroup),c(agesoutgroup),col="grey50",cex=.8,label=c(labelsoutgroup)) for (r in 1:nrow(plotagr)){ segments(plotagr[r,7],plotagr[r,3],plotagr[r,8],plotagr[r,3],col="grey") } text(.2,1.4,col="black",cex=.8,label="A Appalachian") text(.2,1.2,col="black",cex=.8,label="M Midcontinent") text(.2,1.0,col="black",cex=.8,label="S Southern") text(.2,.8,col="black",cex=.8,label="W Western") text(.2,.6,col="grey50",cex=.8,label=" outgroup (UK)") dev.off() ############################################## #Ranked Lists AbundanceList13<-sort(P13[nrow(P13),],decreasing=TRUE) AbundanceList23<-sort(P23[nrow(P23),],decreasing=TRUE) AbundanceList33<-sort(P33[nrow(P33),],decreasing=TRUE) AbundanceList43<-sort(P43[nrow(P43),],decreasing=TRUE) AbundanceList53<-sort(P53[nrow(P53),],decreasing=TRUE) List13<-AbundanceList13[,AbundanceList13>0] List23<-AbundanceList23[,AbundanceList23>0] List33<-AbundanceList33[,AbundanceList33>0] List43<-AbundanceList43[,AbundanceList43>0] List53<-AbundanceList53[,AbundanceList53>0] sum13<-sum(List13[,names(List13)!="Province"&names(List13)!="Age"]) p13<-List13[,names(List13)!="Province"&names(List13)!="Age"]/sum13 sum23<-sum(List23[,names(List23)!="Province"&names(List23)!="Age"]) p23<-List23[,names(List23)!="Province"&names(List23)!="Age"]/sum23 sum33<-sum(List33[,names(List33)!="Province"&names(List33)!="Age"]) p33<-List33[,names(List33)!="Province"&names(List33)!="Age"]/sum33 sum43<-sum(List43[,names(List43)!="Province"&names(List43)!="Age"]) p43<-List43[,names(List43)!="Province"&names(List43)!="Age"]/sum43 sum53<-sum(List53[,names(List53)!="Province"&names(List53)!="Age"]) p53<-List53[,names(List53)!="Province"&names(List53)!="Age"]/sum53 write.csv(p13,"abundancelistAL.csv") write.csv(p23,"abundancelistML.csv") write.csv(p33,"abundancelistSL.csv") write.csv(p43,"abundancelistWL.csv") write.csv(p53,"abundancelistUL.csv") AbundanceList12<-sort(P12[nrow(P12),],decreasing=TRUE) AbundanceList22<-sort(P22[nrow(P22),],decreasing=TRUE) AbundanceList32<-sort(P32[nrow(P32),],decreasing=TRUE) AbundanceList42<-sort(P42[nrow(P42),],decreasing=TRUE) List12<-AbundanceList12[,AbundanceList12>0] List22<-AbundanceList22[,AbundanceList22>0] List32<-AbundanceList32[,AbundanceList32>0] List42<-AbundanceList42[,AbundanceList42>0] sum12<-sum(List12[,names(List12)!="Province"&names(List12)!="Age"]) p12<-List12[,names(List12)!="Province"&names(List12)!="Age"]/sum12 sum22<-sum(List22[,names(List22)!="Province"&names(List22)!="Age"]) p22<-List22[,names(List22)!="Province"&names(List22)!="Age"]/sum22 sum32<-sum(List32[,names(List32)!="Province"&names(List32)!="Age"]) p32<-List32[,names(List32)!="Province"&names(List32)!="Age"]/sum32 sum42<-sum(List42[,names(List42)!="Province"&names(List42)!="Age"]) p42<-List42[,names(List42)!="Province"&names(List42)!="Age"]/sum42 write.csv(p12,"abundancelistAM.csv") write.csv(p22,"abundancelistMM.csv") write.csv(p32,"abundancelistSM.csv") write.csv(p42,"abundancelistWM.csv") AbundanceList11<-sort(P11[nrow(P11),],decreasing=TRUE) AbundanceList21<-sort(P21[nrow(P21),],decreasing=TRUE) AbundanceList31<-sort(P31[nrow(P31),],decreasing=TRUE) AbundanceList41<-sort(P41[nrow(P41),],decreasing=TRUE) List11<-AbundanceList11[,AbundanceList11>0] List21<-AbundanceList21[,AbundanceList21>0] List31<-AbundanceList31[,AbundanceList31>0] List41<-AbundanceList41[,AbundanceList41>0] sum11<-sum(List11[,names(List11)!="Province"&names(List11)!="Age"]) p11<-List11[,names(List11)!="Province"&names(List11)!="Age"]/sum11 sum21<-sum(List21[,names(List21)!="Province"&names(List21)!="Age"]) p21<-List21[,names(List21)!="Province"&names(List21)!="Age"]/sum21 sum31<-sum(List31[,names(List31)!="Province"&names(List31)!="Age"]) p31<-List31[,names(List31)!="Province"&names(List31)!="Age"]/sum31 sum41<-sum(List41[,names(List41)!="Province"&names(List41)!="Age"]) p41<-List41[,names(List41)!="Province"&names(List41)!="Age"]/sum41 write.csv(p11,"abundancelistAE.csv") write.csv(p21,"abundancelistME.csv") write.csv(p31,"abundancelistSE.csv") write.csv(p41,"abundancelistWE.csv") ################################REGIONAL SCALE################################## #Set up Data RAmatrix1<-decostand(AbundanceMatrix[,names(AbundanceMatrix)!="Box"],"total") #Convert abundance matrix to relative abundance Box<-AbundanceMatrix$Box RAmatrix<-cbind(RAmatrix1,Box) #########Abundance lists: Most abundant taxa in each box######## #Make a table of the sums of each genera in a box, divided by the number of individuals, this is basically relative abundance of the aggregate collections of each box. At the end this will spit out a list of the most abundant taxa in each collection AggBoxTable<-as.data.frame(matrix(data=9,nrow=7,ncol=length(names(AbundanceMatrix)), dimnames=list(1:7,c(names(AbundanceMatrix))))) for (i in 1:7){ Box<-AbundanceMatrix[which(AbundanceMatrix$Box==i),] Box<-rbind(Box,colSums(Box)) aggBox<-Box[nrow(Box),] aggBox$Box<-i AggBoxTable[i,]<-aggBox } #Convert this to relative abundance AggBoxTableBox<-AggBoxTable$Box AggBoxTable<-AggBoxTable[,names(AggBoxTable)!="Box"] RAagg<-decostand(AggBoxTable,"total") Boxes<-c(1:7) RAaggbox<-cbind(RAagg,Boxes) names(RAaggbox)[which(names(RAaggbox)=="Boxes")]<-"Box" AbundanceListBox1<-round(sort(RAaggbox[RAaggbox$Box==1,],decreasing=TRUE),digits=2) AbundanceListBox2<-round(sort(RAaggbox[RAaggbox$Box==2,],decreasing=TRUE),digits=2) AbundanceListBox3<-round(sort(RAaggbox[RAaggbox$Box==3,],decreasing=TRUE),digits=2) AbundanceListBox4<-round(sort(RAaggbox[RAaggbox$Box==4,],decreasing=TRUE),digits=2) AbundanceListBox5<-round(sort(RAaggbox[RAaggbox$Box==5,],decreasing=TRUE),digits=2) AbundanceListBox6<-round(sort(RAaggbox[RAaggbox$Box==6,],decreasing=TRUE),digits=2) AbundanceListBox7<-round(sort(RAaggbox[RAaggbox$Box==7,],decreasing=TRUE),digits=2) (ListBox1<-AbundanceListBox1[,AbundanceListBox1>0]) (ListBox2<-AbundanceListBox2[,AbundanceListBox2>0]) (ListBox3<-AbundanceListBox3[,AbundanceListBox3>0]) (ListBox4<-AbundanceListBox4[,AbundanceListBox4>0]) (ListBox5<-AbundanceListBox5[,AbundanceListBox5>0]) #Nashville (ListBoxN<-AbundanceListBox6[,AbundanceListBox6>0]) #Cinci (ListBoxC<-AbundanceListBox7[,AbundanceListBox7>0]) ######Similarity###### #Combines collections from the two called regions, calculates Jaccard's similarity #binary=1 for PA, binary=0 for RA betweenBoxsimilarity<-function(BoxnA,BoxnB,binary){ boxA<-RAmatrix[RAmatrix$Box==BoxnA,colnames(RAmatrix)!="Box"] boxB<-RAmatrix[RAmatrix$Box==BoxnB,colnames(RAmatrix)!="Box"] twoboxes<-rbind(boxA,boxB) if (binary=="FALSE"|binary==0){ similarity<-1-vegdist(twoboxes, method="jaccard", binary=FALSE) } else{ similarity<-1-vegdist(twoboxes, method="jaccard", binary=TRUE) } logitsim<-car::logit(similarity,adjust=.0001) ttest<-t.test(logitsim) sd<-sd(logitsim) results<-c(ttest,sd) names(results)[10]<-"sd" intresults<-c(ttest,sd) names(intresults)[10]<-"sd" results<-c(inv.logit(intresults$estimate),inv.logit(intresults$conf.int[1]),inv.logit(intresults$conf.int[2]),inv.logit(intresults$sd)) names(results)<-c("estimate","LCL","UCL","sd") results } #Fill table:PA PAbetweenBoxtable<-as.data.frame(matrix(data=9,nrow=21, ncol=6, dimnames=list(1:21, c("BoxA","BoxB","mean","LCL","UCL","sd")))) row<-1 for (BoxA in 1:6) { for (BoxB in (BoxA+1):7){ bBs<-betweenBoxsimilarity(BoxA,BoxB,1) PAbetweenBoxtable[row,1]<-BoxA PAbetweenBoxtable[row,2]<-BoxB PAbetweenBoxtable[row,3]<-bBs[1] PAbetweenBoxtable[row,4]<-bBs[2] PAbetweenBoxtable[row,5]<-bBs[3] PAbetweenBoxtable[row,6]<-bBs[4] #increment row counter row<-row+1 } } write.csv(PAbetweenBoxtable,"PAbetweenBoxtable.csv") (PAbetweenBoxtableRound<-round(PAbetweenBoxtable,3)) #Fill table:RA RAbetweenBoxtable<-as.data.frame(matrix(data=9,nrow=21, ncol=6, dimnames=list(1:21, c("BoxA","BoxB","mean","LCL","UCL","sd")))) row<-1 for (BoxA in 1:6) { for (BoxB in (BoxA+1):7){ bBs<-betweenBoxsimilarity(BoxA,BoxB,0) RAbetweenBoxtable[row,1]<-BoxA RAbetweenBoxtable[row,2]<-BoxB RAbetweenBoxtable[row,3]<-bBs[1] RAbetweenBoxtable[row,4]<-bBs[2] RAbetweenBoxtable[row,5]<-bBs[3] RAbetweenBoxtable[row,6]<-bBs[4] #increment row counter row<-row+1 } } write.csv(RAbetweenBoxtable,"RAbetweenBoxtable.csv") (RAbetweenBoxtableRound<-round(RAbetweenBoxtable,3)) ####Plots##### #####PA Plot #Grey lines are 95% confidence intervals; plotted is mean jaccard similarity between labeled collections; pdf(file="PA_plot.pdf") withinpa<-rbind(PAbetweenBoxtable[PAbetweenBoxtable$BoxA>5 & PAbetweenBoxtable$BoxB>5,],PAbetweenBoxtable[PAbetweenBoxtable$BoxA<6 & PAbetweenBoxtable$BoxB<6,]) betweenpa<-PAbetweenBoxtable[PAbetweenBoxtable$BoxA<6 & PAbetweenBoxtable$BoxB>5,] rownames(withinpa)<-c(1:11) betweenSortpa<-betweenpa[c(1,3,5,7,9,2,4,6,8,10),] rownames(betweenSortpa)<-c(12:21) combonamespa<-c(rownames(withinpa),rownames(betweenSortpa)) combomeanspa<-c(withinpa$mean,betweenSortpa$mean) windows() plot(combonamespa~combomeanspa,ylim=c(0,21),type="n",xlim=c(0,max(c(max(withinpa$mean),max(betweenSortpa$mean)))+.05),main="Presence/Absence",ylab="",xlab="Jaccard Similarity") for (r in 1:nrow(withinpa)){ segments(withinpa[r,4],r,withinpa[r,5],r,col="grey") } for (q in 1:nrow(betweenSortpa)){ segments(betweenSortpa[q,4],11+q,betweenSortpa[q,5],11+q,col="grey") } points(withinpa$mean[withinpa$BoxA>5 & withinpa$BoxB>5],rownames(withinpa)[withinpa$BoxA>5 & withinpa$BoxB>5],pch=7,col="grey30") points(withinpa$mean[withinpa$BoxA<6 & withinpa$BoxB<6],rownames(withinpa)[withinpa$BoxA<6 & withinpa$BoxB<6],pch=16,col="black") points(betweenSortpa$mean[betweenSortpa$BoxB==6],rownames(betweenSortpa)[betweenSortpa$BoxB==6],pch=22,col="grey30") points(betweenSortpa$mean[betweenSortpa$BoxB==7],rownames(betweenSortpa)[betweenSortpa$BoxB==7],pch=15,col="grey30") text(withinpa$mean[withinpa$BoxA>5 & withinpa$BoxB>5],.5,col="grey30",cex=.6,label="C,N") withinpalabels<-c("1,2","1,3","1,4","1,5","2,3","2,4","2,5","3,4","3,5","4,5") text(c(withinpa$mean[withinpa$BoxA<6 & withinpa$BoxB<6]),c(seq(1.5,10.5,1)),cex=.6,col="black",label=c(withinpalabels)) betweenpalabels<-c(1:5) text(c(betweenSortpa$mean[betweenSortpa$BoxB==6]),c(seq(11.5,15.5,1)),cex=.6,col="grey30",label=c(betweenpalabels)) text(c(betweenSortpa$mean[betweenSortpa$BoxB==7]),c(seq(16.5,20.5,1)),cex=.6,col="grey30",label=c(betweenpalabels)) text(.18,14,cex=.6,col="grey30",label="vs. Nashville") text(.17,19,cex=.6,col="grey30",label="vs. Cincinnati") text(.09,7,cex=.6,col="black",label="Appalachian") text(.19,1,cex=.6,col="grey30",label="Southern") dev.off() #####RA Plot pdf(file="RA_plot.pdf") withinra<-rbind(RAbetweenBoxtable[RAbetweenBoxtable$BoxA>5 & RAbetweenBoxtable$BoxB>5,],RAbetweenBoxtable[RAbetweenBoxtable$BoxA<6 & RAbetweenBoxtable$BoxB<6,]) betweenra<-RAbetweenBoxtable[RAbetweenBoxtable$BoxA<6 & RAbetweenBoxtable$BoxB>5,] rownames(withinra)<-c(1:11) betweenSortra<-betweenra[c(1,3,5,7,9,2,4,6,8,10),] rownames(betweenSortra)<-c(12:21) combonamesra<-c(rownames(withinra),rownames(betweenSortra)) combomeansra<-c(withinra$mean,betweenSortra$mean) windows() plot(combonamesra~combomeansra,ylim=c(0,21),type="n",xlim=c(0,max(c(max(withinra$mean),max(betweenSortra$mean)))+.05),main="Relative Abundance",ylab="",xlab="Jaccard Similarity") for (r in 1:nrow(withinra)){ segments(withinra[r,4],r,withinra[r,5],r,col="grey") } for (q in 1:nrow(betweenSortra)){ segments(betweenSortra[q,4],11+q,betweenSortra[q,5],11+q,col="grey") } points(withinra$mean[withinra$BoxA>5 & withinra$BoxB>5],rownames(withinra)[withinra$BoxA>5 & withinra$BoxB>5],pch=7,col="grey30") points(withinra$mean[withinra$BoxA<6 & withinra$BoxB<6],rownames(withinra)[withinra$BoxA<6 & withinra$BoxB<6],pch=16,col="black") points(betweenSortra$mean[betweenSortra$BoxB==6],rownames(betweenSortra)[betweenSortra$BoxB==6],pch=22,col="grey30") points(betweenSortra$mean[betweenSortra$BoxB==7],rownames(betweenSortra)[betweenSortra$BoxB==7],pch=15,col="grey30") text(withinra$mean[withinra$BoxA>5 & withinra$BoxB>5],.5,col="grey30",cex=.6,label="C,N") withinralabels<-c("1,2","1,3","1,4","1,5","2,3","2,4","2,5","3,4","3,5","4,5") text(c(withinra$mean[withinra$BoxA<6 & withinra$BoxB<6]),c(seq(1.5,10.5,1)),cex=.6,col="black",label=c(withinralabels)) betweenralabels<-c(1:5) text(c(betweenSortra$mean[betweenSortra$BoxB==6]),c(seq(11.5,15.5,1)),cex=.6,col="grey30",label=c(betweenralabels)) text(c(betweenSortra$mean[betweenSortra$BoxB==7]),c(seq(16.5,20.5,1)),cex=.6,col="grey30",label=c(betweenralabels)) text(.10,14,cex=.6,col="grey30",label="vs. Nashville") text(.25,19,cex=.6,col="grey30",label="vs. Cincinnati") text(.16,8,cex=.6,col="black",label="Appalachian") text(.18,1,cex=.6,col="grey30",label="Southern") dev.off() ################################### ###MDS### #MDS #Run MDS mMds1<-metaMDS(RAmatrix[,names(RAmatrix)!="Box"],k=1) mMds2<-metaMDS(RAmatrix[,names(RAmatrix)!="Box"],k=2) mMds3<-metaMDS(RAmatrix[,names(RAmatrix)!="Box"],k=3) mMds4<-metaMDS(RAmatrix[,names(RAmatrix)!="Box"],k=4) mMds5<-metaMDS(RAmatrix[,names(RAmatrix)!="Box"],k=5) mMds6<-metaMDS(RAmatrix[,names(RAmatrix)!="Box"],k=6) mMds7<-metaMDS(RAmatrix[,names(RAmatrix)!="Box"],k=7) #Scree plot, pick number of axes stressvector<-c(mMds1$stress,mMds2$stress,mMds3$stress,mMds4$stress,mMds5$stress,mMds6$stress,mMds7$stress) windows() plot(stressvector,pch=16,las=1,xlab="k",ylab="stress") #deciding to use k=3 #Species effects windows() plot(mMds3, type="t", display=c("species"),cex=.7) #plot(mMds4, type="t", display=c("species"),cex=.7) #MDS sites numbered (k=4) #windows() #plot(mMds4,type="n") #text(mMds4$points[which(rownames(mMds4$points)%in%rownames(RAmatrix)[RAmatrix$Box==6]),],col="grey30",labels="N",cex=.8) #text(mMds4$points[which(rownames(mMds4$points)%in%rownames(RAmatrix)[RAmatrix$Box==7]),], col="grey30", labels="C",cex=.8) #text(mMds4$points[which(rownames(mMds4$points)%in%rownames(RAmatrix)[RAmatrix$Box==2]),], col="black",labels="2",cex=.8) #text(mMds4$points[which(rownames(mMds4$points)%in%rownames(RAmatrix)[RAmatrix$Box==3]),], col="black",labels="3",cex=.8) #text(mMds4$points[which(rownames(mMds4$points)%in%rownames(RAmatrix)[RAmatrix$Box==4]),], col="black",labels="4",cex=.8) #text(mMds4$points[which(rownames(mMds4$points)%in%rownames(RAmatrix)[RAmatrix$Box==5]),], col="black",labels="5",cex=.8) #text(mMds4$points[which(rownames(mMds4$points)%in%rownames(RAmatrix)[RAmatrix$Box==1]),], col="black",labels="1",cex=.8) #text(-.3,-1.2,label="Southern",col="grey30",cex=.8) #text(1.7,-.5, label="Appalachian",col="black",cex=.8) #MDS sites numbered (k=3) pdf(file="MDS_plain.pdf") #windows() plot(mMds3,type="n") text(mMds3$points[which(rownames(mMds3$points)%in%rownames(RAmatrix)[RAmatrix$Box==6]),],col="grey30",labels="N",cex=.8) text(mMds3$points[which(rownames(mMds3$points)%in%rownames(RAmatrix)[RAmatrix$Box==7]),], col="grey30", labels="C",cex=.8) text(mMds3$points[which(rownames(mMds3$points)%in%rownames(RAmatrix)[RAmatrix$Box==2]),], col="black",labels="2",cex=.8) text(mMds3$points[which(rownames(mMds3$points)%in%rownames(RAmatrix)[RAmatrix$Box==3]),], col="black",labels="3",cex=.8) text(mMds3$points[which(rownames(mMds3$points)%in%rownames(RAmatrix)[RAmatrix$Box==4]),], col="black",labels="4",cex=.8) text(mMds3$points[which(rownames(mMds3$points)%in%rownames(RAmatrix)[RAmatrix$Box==5]),], col="black",labels="5",cex=.8) text(mMds3$points[which(rownames(mMds3$points)%in%rownames(RAmatrix)[RAmatrix$Box==1]),], col="black",labels="1",cex=.8) text(-.3,-1.2,label="Southern",col="grey30",cex=.8) text(1.7,-.5, label="Appalachian",col="black",cex=.8) dev.off() #MDS chull pdf(file="mds_circles.pdf") plot(mMds3,type="n") mdspointsn<-mMds3$points[which(rownames(mMds3$points)%in%rownames(RAmatrix)[RAmatrix$Box==6]),] mdspointsc<-mMds3$points[which(rownames(mMds3$points)%in%rownames(RAmatrix)[RAmatrix$Box==7]),] mdspoints2<-mMds3$points[which(rownames(mMds3$points)%in%rownames(RAmatrix)[RAmatrix$Box==2]),] mdspoints3<-mMds3$points[which(rownames(mMds3$points)%in%rownames(RAmatrix)[RAmatrix$Box==3]),] mdspoints4<-mMds3$points[which(rownames(mMds3$points)%in%rownames(RAmatrix)[RAmatrix$Box==4]),] mdspoints5<-mMds3$points[which(rownames(mMds3$points)%in%rownames(RAmatrix)[RAmatrix$Box==5]),] mdspoints1<-mMds3$points[which(rownames(mMds3$points)%in%rownames(RAmatrix)[RAmatrix$Box==1]),] mdspn<-mdspointsn[,1:2] mdspc<-mdspointsc[,1:2] mdsp2<-mdspoints2[,1:2] mdsp3<-mdspoints3[,1:2] mdsp4<-mdspoints4[,1:2] mdsp5<-mdspoints5[,1:2] mdsp1<-mdspoints1[,1:2] circlesn<-chull(mdspn) circlesn<-c(circlesn,circlesn[1]) circlesc<-chull(mdspc) circlesc<-c(circlesc,circlesc[1]) circles2<-chull(mdsp2) circles2<-c(circles2,circles2[1]) circles3<-chull(mdsp3) circles3<-c(circles3,circles3[1]) circles4<-chull(mdsp4) circles4<-c(circles4,circles4[1]) circles5<-chull(mdsp5) circles5<-c(circles5,circles5[1]) circles1<-chull(mdsp1) circles1<-c(circles1,circles1[1]) text(mdspointsn,col="grey30",labels="N",cex=.8) text(mdspointsc, col="grey30", labels="C",cex=.8) text(mdspoints2, col="black",labels="2",cex=.8) text(mdspoints3, col="black",labels="3",cex=.8) text(mdspoints4, col="black",labels="4",cex=.8) text(mdspoints5, col="black",labels="5",cex=.8) text(mdspoints1, col="black",labels="1",cex=.8) text(-.3,-1.2,label="Southern",col="grey30",cex=.8) text(1.7,-.5, label="Appalachian",col="black",cex=.8) lines(mdspn[circlesn,]) lines(mdspc[circlesc,]) lines(mdsp1[circles1,]) lines(mdsp2[circles2,]) lines(mdsp3[circles3,]) lines(mdsp4[circles4,]) lines(mdsp5[circles5,]) dev.off() ########################################### #Distance paradist<-read.csv("paradist.csv") attach(paradist) pdf(file="paradist_allonone.pdf") windows() plot(RA~distr,pch=16,type="n",las=1,ylim=c(0,max(paradist$RAUCL)),ylab="Similarity (RA)", xlab="distance (km)") #Add confidence intervals, r is simply a row counter for (r in 1:nrow(paradist)){ segments(paradist[r,4],paradist[r,7],paradist[r,4],paradist[r,8],col="grey") } points(distr[box1=="C"|box2=="C"],RA[box1=="C"|box2=="C"],pch=15,col="grey30") points(distr[box1=="N"|box2=="N"],RA[box1=="N"|box2=="N"],pch=22,col="grey30") points(distr[box1=="C"&box2=="N"],RA[box1=="C"&box2=="N"],pch=7,col="white") points(distr[box2!="C"&box2!="N"],RA[box2!="C"&box2!="N"],pch=16,col="black") abline(lm(RA~distr),lty=2) averageR2<-summary(lm(RA~distr))$adj.r.squared text(800,.3,label= paste("R =",round(averageR2,3))) dev.off() ##### pdf(file="paradist_threeplots.pdf") windows() par(mfrow=c(1,3)) plot(RA~distr,pch=16,type="n",las=1,ylim=c(0,max(paradist$RAUCL)),ylab="Similarity (RA)", xlab="distance (km)", main="Cincinnati/Appalachian") #Add confidence intervals, r is simply a row counter for (r in 1:nrow(paradist)){ if (box1[r]=="C"|box2[r]=="C"){ segments(paradist[r,4],paradist[r,7],paradist[r,4],paradist[r,8],col="grey") } } points(distr[box1=="C"|box2=="C"],RA[box1=="C"|box2=="C"],pch=15,col="grey30") abline(lm(RA~distr),lty=2,col="grey") abline(lm(RA[box1=="C"|box2=="C"]~distr[box1=="C"|box2=="C"]),lty=2) cinR2<-summary(lm(RA[box1=="C"|box2=="C"]~distr[box1=="C"|box2=="C"]))$adj.r.squared text(800,.3,label= paste("R =",round(averageR2,3)),col="grey") text(800,.28,label= paste("R =",round(cinR2,3))) ##### plot(RA~distr,pch=16,type="n",las=1,ylim=c(0,max(paradist$RAUCL)),ylab="Similarity (RA)", xlab="distance (km)", main="Appalachian/Appalachian") #Add confidence intervals, r is simply a row counter for (r in 1:nrow(paradist)){ if (box2[r]!="C"&box2[r]!="N"){ segments(paradist[r,4],paradist[r,7],paradist[r,4],paradist[r,8],col="grey") } } points(distr[box2!="C"&box2!="N"],RA[box2!="C"&box2!="N"],pch=16,col="black") abline(lm(RA~distr),lty=2,col="grey") abline(lm(RA[box2!="C"&box2!="N"]~distr[box2!="C"&box2!="N"]),lty=2) appR2<-summary(lm(RA[box2!="C"&box2!="N"]~distr[box2!="C"&box2!="N"]))$adj.r.squared text(800,.3,label= paste("R =",round(averageR2,3)),col="grey") text(800,.28,label= paste("R =",round(appR2,3))) dev.off() detach(paradist)