##Primary code for "Conith et al. #Evolution of a soft-tissue foraging adaptation in African cichlids: #Roles for novelty, convergence, and constraint ####Data set-up#### #Set your working directory setwd("/Users/Home/Desktop/Supp") #Use this source script to install/import packages and load in a wrapper to import your data and tree source("Scripts/Ectodini_SOURCE.R") ####Data Means#### setwd("/Users/Home/Desktop/Supp") DatMean<-read.csv("OriginalData/MalawiTangNasutaMaleOnly.csv") #PC Scores# d1<-unique(ddply(DatMean, .(PhyloName), summarize, mean=Lake)) d2<-ddply(DatMean, .(PhyloName), summarize, mean=mean(BD)) d3<-ddply(DatMean, .(PhyloName), summarize, mean=mean(Flap.Length)) d4<-ddply(DatMean, .(PhyloName), summarize, mean=mean(Flap.Depth)) Dat<-cbind(d1, d2[,2], d3[,2], d4[,2]) colnames(Dat)<-c("Species", "Lake", "BD", "Flap_Length", "Flap_Depth") write.csv(Dat, "CichlidMaleMean_20180618.csv", row.names = F) #Tree and data setup Files<-MakeFiles() #Note that MakeeFiles will only work with the nexus file (LTandLM_con.nex) #LTandLM_Comb.trees is the full posterior distribution of trees from 4 BEAST runs #You can generate this file using the attached BEAST XML file #Your tree and data is stored as 'Tree' and 'Data' respectively. Tree<-Files[[1]] Data<-Files[[2]] #For Cichlid numeric data listed as factors for(i in 2:ncol(Data)){ Data[,i]<-as.numeric(as.character(Data[,i])) } ###Phylogenetically corrected residuals ##Assess model fit and extract residuals based on it #Depth pglsModelOUDepth <- gls(Flap_Depth ~ BD, correlation = corMartins(1, phy = Tree), data = Data, method = "ML") pglsModelBMDepth <- gls(Flap_Depth ~ BD, correlation = corBrownian(1, phy = Tree), data = Data, method = "ML") pglsModelPaDepth <- gls(Flap_Depth ~ BD, correlation = corPagel(1, phy = Tree, fixed = T), data = Data, method = "ML") pglsModelBlomDepth <- gls(Flap_Depth ~ BD, correlation = corBlomberg(1, phy = Tree, fixed = T), data = Data, method = "ML") summary(pglsModelOUDepth) summary(pglsModelBMDepth) summary(pglsModelPaDepth) summary(pglsModelBlomDepth) DepthOU<-pglsModelOUDepth$residuals[1:length(pglsModelOUDepth$residuals)] Data$FlapDepthBDCorr<-as.numeric(DepthOU) #Length pglsModelOULength <- gls(Flap_Length ~ BD, correlation = corMartins(1, phy = Tree), data = Data, method = "ML") pglsModelBMLength <- gls(Flap_Length ~ BD, correlation = corBrownian(1, phy = Tree), data = Data, method = "ML") pglsModelPaLength <- gls(Flap_Length ~ BD, correlation = corPagel(1, phy = Tree, fixed = T), data = Data, method = "ML") pglsModelBlomLength <- gls(Flap_Length ~ BD, correlation = corBlomberg(1, phy = Tree, fixed = T), data = Data, method = "ML") summary(pglsModelOULength) summary(pglsModelBMLength) summary(pglsModelPaLength) summary(pglsModelBlomLength) LengthOU<-pglsModelOULength$residuals[1:length(pglsModelOULength$residuals)] Data$FlapLengthBDCorr<-as.numeric(LengthOU) write.csv(Data, "FlapFemaleResiduals_20180618.csv") write.csv(Data, "FlapMaleResiduals_20180618.csv") ####Convergence#### #Male setwd("/Users/Home/Desktop/Supp") Files<-MakeFiles() #Your tree and data is stored as 'Tree' and 'Data' respectively. Tree<-Files[[1]] Data<-Files[[2]] #For Cichlids as listed as factors for(i in 6:ncol(Data)){ Data[,i]<-as.numeric(as.character(Data[,i])) } Flaps<-c("Asprotilapia_leptura", "Labeotropheus_fuelleborni", "Labeotropheus_trewavasae", "Ophthalmotilapia_nasuta") Flaps<-c("Asprotilapia_leptura", "Labeotropheus_fuelleborni", "Labeotropheus_trewavasae") Flaps<-c("Asprotilapia_leptura", "Ophthalmotilapia_nasuta") FlapConvData<-Data[,c(10,9)] par(pty='s') Flaps_ConNumSig<-convnumsig(Tree, FlapConvData, Flaps, 1000) ####Correlations#### ##Correlation w/o extremes #Tree and data setup setwd("/Users/Home/Desktop/Supp") Tr<-read.nexus("LTandLM_con.nex") Tre<-drop.tip(Tr, tip=c("Asprotilapia_leptura", "Labeotropheus_trewavasae", "Labeotropheus_fuelleborni")) EctoData<-read.csv("CichlidFemaleMean_20181126.csv", row.names = 1) Pruning<-treedata(Tre, EctoData) EctoTreePrune<-Pruning$phy EctoTreePrune<-ladderize(EctoTreePrune) EctoDataPrune<-Pruning$data EctoDataPrune<- EctoDataPrune[match(EctoTreePrune$tip.label, rownames(EctoDataPrune)),] EctoData<-as.data.frame(EctoDataPrune) EctoTree<-ladderize(EctoTreePrune,T) #For Cichlids as listed as factors for(i in 6:ncol(EctoData)){ EctoData[,i]<-as.numeric(as.character(EctoData[,i])) } #Extreme removed EctoData$Species<-rownames(EctoData) ComData<-comparative.data(EctoTree, EctoData, 'Species') LengDepFit<-pgls(FlapDepthBDCorr~FlapLengthBDCorr, data=ComData, lambda='ML') summary(LengDepFit) plot(ComData$data$FlapDepthBDCorr~ComData$data$FlapLengthBDCorr, pch=19, xlab="Relative Snout Length", ylab="Relative Snout Depth") abline(LengDepFit, lwd=2) #Now with all species Data$Species<-rownames(Data) ComData<-comparative.data(Tree, Data, 'Species') FullFit<-pgls(FlapDepthBDCorr~FlapLengthBDCorr, data=ComData, lambda='ML') summary(FullFit) plot(ComData$data$FlapDepthBDCorr~ComData$data$FlapLengthBDCorr, pch=19, xlab="Relative Snout Length", ylab="Relative Snout Depth") abline(LengDepFit, lwd=2) ####Evolutionary Models#### setwd("/Users/Home/Desktop/Supp") source("Scripts/Ectodini_SOURCE.R") #Analysis with residual corrections #Tree and data setup Files<-MakeFiles() #Your tree and data is stored as 'Tree' and 'Data' respectively. Tree<-Files[[1]] Data<-Files[[2]] #For Cichlids as listed as factors for(i in 6:ncol(Data)){ Data[,i]<-as.numeric(as.character(Data[,i])) } #Generate trees with stocastically mapped characters over a posterior distribution. ManyTrees<-read.nexus("LTandLM_Comb.trees") #Randomly subsample from a posterior distribution of trees SampleTrees<-sample(ManyTrees[10004:32004],100) #Drop additional taxa from the trees SampleDropTrees<-lapply(SampleTrees,drop.tip,tip=c("Aulonocranus_dewindti", "Callochromis_pleurospilus", "Callochromis_stappersii", "Ectodus_descampsii", "Enantiopus_melanogenys", "Microdontochromis_rotundiventralis", "Oreochromis_tanganicae", "Xenotilapia_cf_bathyphila", "Xenotilapia_flavipinnis", "Xenotilapia_papilio_sunflower")) class(SampleDropTrees)<-"multiPhylo" ###Diet### FlapDiet<-Data$Diet SimmapTreesDiet<-make.simmap(SampleDropTrees, FlapDiet, nsim=10, message=FALSE) save(SimmapTreesDiet, file = "Female_Diet_SIMMAP_OUwie.rda") MyOUwieModelsDepthDiet<-OUwieModels(SimmapTreesDiet, Data, "FlapDepthBDCorr", "Diet", "OUM") save(MyOUwieModelsDepthDiet, file = "Female_Diet_Depth_OUwie.rda") MyOUwieModelsLengthDiet<-OUwieModels(SimmapTreesDiet, Data, "FlapLengthBDCorr", "Diet", "OUM") save(MyOUwieModelsLengthDiet, file = "Female_Diet_Length_OUwie.rda") ###Habitat### FlapHabitat<-Data$Habitat SimmapTreesHabitat<-make.simmap(SampleDropTrees, FlapHabitat, nsim=10, message=FALSE) save(SimmapTreesHabitat, file = "Female_Habitat_SIMMAP_OUwie.rda") MyOUwieModelsDepthHabitat<-OUwieModels(SimmapTreesHabitat, Data, "FlapDepthBDCorr", "Habitat", "OUM") save(MyOUwieModelsDepthHabitat, file = "Female_Habitat_Depth_OUwie.rda") MyOUwieModelsLengthHabitat<-OUwieModels(SimmapTreesHabitat, Data, "FlapLengthBDCorr", "Habitat", "OUM") save(MyOUwieModelsLengthHabitat, file = "Female_Habitat_Length_OUwie.rda") ###Feeding### FlapFeeding<-Data$Feeding SimmapTreesFeeding<-make.simmap(SampleDropTrees, FlapFeeding, nsim=10, message=FALSE) save(SimmapTreesFeeding, file = "Female_Feeding_SIMMAP_OUwie.rda") MyOUwieModelsDepthFeeding<-OUwieModels(SimmapTreesFeeding, Data, "FlapDepthBDCorr", "Feeding", "OUM") save(MyOUwieModelsDepthFeeding, file = "Female_Feeding_Depth_OUwie.rda") MyOUwieModelsLengthFeeding<-OUwieModels(SimmapTreesFeeding, Data, "FlapLengthBDCorr", "Feeding", "OUM") save(MyOUwieModelsLengthFeeding, file = "Female_Feeding_Length_OUwie.rda") ###Saving OUwie### #Diet depth OUwieResults<-matrix(NA,length(MyOUwieModelsDepthDiet),2) for(i in 1:length(MyOUwieModelsDepthDiet)){ OUwieResults[i,1] <-MyOUwieModelsDepthDiet[[i]]$loglik OUwieResults[i,2] <-MyOUwieModelsDepthDiet[[i]]$AICc } colnames(OUwieResults)<-c("log likelihood", "AICc") write.csv(OUwieResults, "Female_Diet_Depth_OUwieResults.csv", row.names=F) #Diet length OUwieResults<-matrix(NA,length(MyOUwieModelsLengthDiet),2) for(i in 1:length(MyOUwieModelsLengthDiet)){ OUwieResults[i,1] <-MyOUwieModelsLengthDiet[[i]]$loglik OUwieResults[i,2] <-MyOUwieModelsLengthDiet[[i]]$AICc } colnames(OUwieResults)<-c("log likelihood", "AICc") write.csv(OUwieResults, "Female_Diet_Length_OUwieResults.csv", row.names=F) #Habitat depth OUwieResults<-matrix(NA,length(MyOUwieModelsDepthHabitat),2) for(i in 1:length(MyOUwieModelsDepthHabitat)){ OUwieResults[i,1] <-MyOUwieModelsDepthHabitat[[i]]$loglik OUwieResults[i,2] <-MyOUwieModelsDepthHabitat[[i]]$AICc } colnames(OUwieResults)<-c("log likelihood", "AICc") write.csv(OUwieResults, "Female_Habitat_Depth_OUwieResults.csv", row.names=F) #Habitat length OUwieResults<-matrix(NA,length(MyOUwieModelsLengthHabitat),2) for(i in 1:length(MyOUwieModelsLengthHabitat)){ OUwieResults[i,1] <-MyOUwieModelsLengthHabitat[[i]]$loglik OUwieResults[i,2] <-MyOUwieModelsLengthHabitat[[i]]$AICc } colnames(OUwieResults)<-c("log likelihood", "AICc") write.csv(OUwieResults, "Female_Habitat_Length_OUwieResults.csv", row.names=F) #Feeding depth OUwieResults<-matrix(NA,length(MyOUwieModelsDepthFeeding),2) for(i in 1:length(MyOUwieModelsDepthFeeding)){ OUwieResults[i,1] <-MyOUwieModelsDepthFeeding[[i]]$loglik OUwieResults[i,2] <-MyOUwieModelsDepthFeeding[[i]]$AICc } colnames(OUwieResults)<-c("log likelihood", "AICc") write.csv(OUwieResults, "Female_Feeding_Depth_OUwieResults.csv", row.names=F) #Feeding length OUwieResults<-matrix(NA,length(MyOUwieModelsLengthFeeding),2) for(i in 1:length(MyOUwieModelsLengthFeeding)){ OUwieResults[i,1] <-MyOUwieModelsLengthFeeding[[i]]$loglik OUwieResults[i,2] <-MyOUwieModelsLengthFeeding[[i]]$AICc } colnames(OUwieResults)<-c("log likelihood", "AICc") write.csv(OUwieResults, "Female_Feeding_Length_OUwieResults.csv", row.names=F) ###OUCH Evolutionary Models### #Note - only intrested in Brownian and single-peak OU SampleTrees<-sample(ManyTrees[10004:32004],1000) #Drop additional taxa from the trees SampleDropTrees<-lapply(SampleTrees,drop.tip,tip=c("Aulonocranus_dewindti", "Callochromis_pleurospilus", "Callochromis_stappersii", "Ectodus_descampsii", "Enantiopus_melanogenys", "Microdontochromis_rotundiventralis", "Oreochromis_tanganicae", "Xenotilapia_cf_bathyphila", "Xenotilapia_flavipinnis", "Xenotilapia_papilio_sunflower")) class(SampleDropTrees)<-"multiPhylo" #Single Variable OuchOutputDepth = lapply(SampleDropTrees, function(x) { OuchModels(x, Data, "FlapDepthBDCorr", "Feeding") }) OuchOutputLength = lapply(SampleDropTrees, function(x) { OuchModels(x, Data, "FlapLengthBDCorr", "Feeding") }) save(OuchOutputDepth, file = "Female_Depth_OUCH.rda") save(OuchOutputLength, file = "Female_Length_OUCH.rda") #OUCH saving #Depth OUCHResults<-matrix(NA,length(OuchOutputDepth),4) for(i in 1:length(OuchOutputDepth)){ OUCHResults[i,1] <-OuchOutputDepth[[i]][1] OUCHResults[i,2] <-OuchOutputDepth[[i]][2] OUCHResults[i,3] <-OuchOutputDepth[[i]][4] OUCHResults[i,4] <-OuchOutputDepth[[i]][5] } colnames(OUCHResults)<-c("log likelihood BM", "log likelihood SP OU", "AICc BM", "AICc SP OU") write.csv(OUCHResults, "Female_Depth_OUCHresults.csv", row.names=F) #Length OUCHResults<-matrix(NA,length(OuchOutputLength),4) for(i in 1:length(OuchOutputLength)){ OUCHResults[i,1] <-OuchOutputLength[[i]][1] OUCHResults[i,2] <-OuchOutputLength[[i]][2] OUCHResults[i,3] <-OuchOutputLength[[i]][4] OUCHResults[i,4] <-OuchOutputLength[[i]][5] } colnames(OUCHResults)<-c("log likelihood BM", "log likelihood SP OU", "AICc BM", "AICc SP OU") write.csv(OUCHResults, "Female_Length_OUCHresults.csv", row.names=F) ####Evolvability Test#### setwd("/Users/Home/Desktop/Supp") source("Scripts/Ectodini_SOURCE.R") ##Remove Extreme depth individuals #Tree and data setup Tr<-read.nexus("LTandLM_con.nex") Tre<-drop.tip(Tr, tip=c("Asprotilapia_leptura", "Labeotropheus_trewavasae", "Labeotropheus_fuelleborni")) EctoData<-read.csv("CichlidFemaleMean_20181126.csv", row.names = 1) Pruning<-treedata(Tre, EctoData) EctoTreePrune<-Pruning$phy EctoTreePrune<-ladderize(EctoTreePrune) EctoDataPrune<-Pruning$data EctoDataPrune<- EctoDataPrune[match(EctoTreePrune$tip.label, rownames(EctoDataPrune)),] EctoData<-as.data.frame(EctoDataPrune) EctoTree<-ladderize(EctoTreePrune,T) #For Cichlids as listed as factors for(i in 6:ncol(EctoData)){ EctoData[,i]<-as.numeric(as.character(EctoData[,i])) } #Evolutionary Models W/O Extreme Individuals# #Generate trees with stocastically mapped characters over a posterior distribution. ManyTrees<-read.nexus("~/Documents/LTandLM_Comb.trees") #Randomly subsample from a posterior distribution of trees SampleTrees<-sample(ManyTrees[10004:32004],1000) #Drop additional taxa from the trees - inc. depth exception species (w/o nasuta) SampleDropTrees<-lapply(SampleTrees,drop.tip,tip=c("Aulonocranus_dewindti", "Callochromis_pleurospilus", "Callochromis_stappersii", "Ectodus_descampsii", "Enantiopus_melanogenys", "Microdontochromis_rotundiventralis", "Oreochromis_tanganicae", "Xenotilapia_cf_bathyphila", "Xenotilapia_flavipinnis", "Xenotilapia_papilio_sunflower", "Asprotilapia_leptura", "Labeotropheus_trewavasae", "Labeotropheus_fuelleborni")) class(SampleDropTrees)<-"multiPhylo" save(SampleDropTrees, file = "SampleDropTrees_Alpha.rda") ###Gain depth and length evolutionary parameters### ###Diet### FlapDiet<-EctoData$Diet SimmapAlphaTreesDiet<-make.simmap(SampleDropTrees, FlapDiet, nsim=1, message=FALSE) save(SimmapAlphaTreesDiet, file = "Female_Diet_SIMMAPAlpha_OUwie.rda") MyOUwieModelsDepthAlphaDiet<-OUwieModels(SimmapAlphaTreesDiet, EctoData, "FlapDepthBDCorr", "Diet", "OUM") save(MyOUwieModelsDepthAlphaDiet, file = "Female_Diet_DepthAlpha_OUwie.rda") #MyOUwieModelsLengthAlphaDiet<-OUwieModels(SimmapAlphaTreesDiet, EctoData, "FlapLengthBDCorr", "Diet", "OUM") #save(MyOUwieModelsLengthAlphaDiet, file = "Female_Diet_LengthAlpha_OUwie.rda") ###Habitat### FlapHabitat<-EctoData$Habitat SimmapAlphaTreesHabitat<-make.simmap(SampleDropTrees, FlapHabitat, nsim=1, message=FALSE) save(SimmapAlphaTreesHabitat, file = "Female_Habitat_SIMMAPAlpha_OUwie.rda") MyOUwieModelsDepthAlphaHabitat<-OUwieModels(SimmapAlphaTreesHabitat, EctoData, "FlapDepthBDCorr", "Habitat", "OUM") save(MyOUwieModelsDepthAlphaHabitat, file = "Female_Habitat_DepthAlpha_OUwie.rda") #MyOUwieModelsLengthAlphaHabitat<-OUwieModels(SimmapAlphaTreesHabitat, EctoData, "FlapLengthBDCorr", "Habitat", "OUM") #save(MyOUwieModelsLengthAlphaHabitat, file = "Female_Habitat_LengthAlpha_OUwie.rda") ###Feeding### FlapFeeding<-EctoData$Feeding SimmapAlphaTreesFeeding<-make.simmap(SampleDropTrees, FlapFeeding, nsim=1, message=FALSE) save(SimmapAlphaTreesFeeding, file = "Female_Feeding_SIMMAPAlpha_OUwie.rda") MyOUwieModelsDepthAlphaFeeding<-OUwieModels(SimmapAlphaTreesFeeding, EctoData, "FlapDepthBDCorr", "Feeding", "OUM") save(MyOUwieModelsDepthAlphaFeeding, file = "Female_Feeding_DepthAlpha_OUwie.rda") #MyOUwieModelsLengthAlphaFeeding<-OUwieModels(SimmapAlphaTreesFeeding, EctoData, "FlapLengthBDCorr", "Feeding", "OUM") #save(MyOUwieModelsLengthAlphaFeeding, file = "Female_Feeding_LengthAlpha_OUwie.rda") ###OUwie Single Peak/Single Rate set up### EctoDatas<-EctoData EctoDatas$Species<-rownames(EctoDatas) MyData<-EctoDatas[,c("Species", "Feeding", "FlapDepthBDCorr")] #Need to add a null two regime column for it to work MyData$Feeding<-as.factor(rep(x = c("SP","SR"), nrow(MyData)/2)) #Single peak OU OUwieWOextremeOU = mclapply(SampleDropTrees, function(x) { OUwie(x, MyData, "OU1") }) save(OUwieWOextremeOU, file = "SingleOU_Results_SimAlphaOUwie.rda") #Now extract parameters and simulate on the full tree setwd("/Users/Home/Desktop/Supp") source("Scripts/Ectodini_SOURCE.R") #Analysis with residual corrections #Tree and data setup Files<-MakeFiles() #Your tree and data is stored as 'Tree' and 'Data' respectively. Tree<-Files[[1]] Data<-Files[[2]] #For Cichlids as listed as factors for(i in 6:ncol(Data)){ Data[,i]<-as.numeric(as.character(Data[,i])) } ###Simulate with Depth### ###W/ Labeotropheus_fuelleborni, Labeotropheus trewavasae, Asprotilapia leptura ###Parameters from evolultionary models without extreme individuals ##Diet #Extract diet data and perform simulation on full tree DietSimList<-rep(list(list()), 1000) for(i in 1:1000){ ThetaP<-setNames(MyOUwieModelsDepthAlphaDiet[[i]]$theta[,1],c("algae","large evasive","mollusks","small evasive")) DietSimList[[i]]<-multiOU(tree = SimmapTreesDiet[[i]], alpha = MyOUwieModelsDepthAlphaDiet[[i]]$solution[1,], sig2 = MyOUwieModelsDepthAlphaDiet[[i]]$solution[2,], theta = ThetaP) } ##Feeding #Extract feeding data and perform simulation FeedingSimList<-rep(list(list()), 1000) for(i in 1:1000){ ThetaP<-setNames(MyOUwieModelsDepthAlphaFeeding[[i]]$theta[,1],c("crushing","scraper","sifter","suction")) FeedingSimList[[i]]<-multiOU(tree = SimmapTreesFeeding[[i]], alpha = MyOUwieModelsDepthAlphaFeeding[[i]]$solution[1,], sig2 = MyOUwieModelsDepthAlphaFeeding[[i]]$solution[2,], theta = ThetaP) } ##Habitat #Extract habitat data and perform simulation HabitatSimList<-rep(list(list()), 1000) for(i in 1:1000){ ThetaP<-setNames(MyOUwieModelsDepthAlphaHabitat[[i]]$theta[,1],c("intermediate","open water","rocky","sandy")) HabitatSimList[[i]]<-multiOU(tree = SimmapTreesHabitat[[i]], alpha = MyOUwieModelsDepthAlphaHabitat[[i]]$solution[1,], sig2 = MyOUwieModelsDepthAlphaHabitat[[i]]$solution[2,], theta = ThetaP) } ##Single Peak #Extract single optima and perform simulation SPSimList<-rep(list(list()), 1000) for(i in 1:1000){ SPSimList[[i]]<-fastBM(tree = SampleDropTrees[[i]], sig2 = OUwieWOextremeOU[[i]]$solution[2], alpha = OUwieWOextremeOU[[i]]$solution[1], theta = OUwieWOextremeOU[[i]]$theta[1]) } ##Plots DepthData<-setNames(Data$FlapDepthBDCorr,rownames(Data)) DepthData[c("Labeotropheus_fuelleborni","Labeotropheus_trewavasae","Asprotilapia_leptura")] par(mfrow=c(2,2), pty='s') #Single Peak plot(density(DepthData[c("Labeotropheus_fuelleborni","Labeotropheus_trewavasae","Asprotilapia_leptura")]), xlim=c(-2,2), ylim=c(0,5), main="Single Peak") for(i in 1:1000){ lines(density(as.numeric(SPSimList[[i]])), col="dark green") } #Diet plot(density(DepthData[c("Labeotropheus_fuelleborni","Labeotropheus_trewavasae","Asprotilapia_leptura")]), xlim=c(-2,2), ylim=c(0,5), main="Diet") for(i in 1:1000){ lines(density(as.numeric(DietSimList[[i]])), col="red") } #Feeding plot(density(DepthData[c("Labeotropheus_fuelleborni","Labeotropheus_trewavasae","Asprotilapia_leptura")]), xlim=c(-2,2), ylim=c(0,5), main="Feeding") for(i in 1:1000){ lines(density(as.numeric(FeedingSimList[[i]])), col="blue") } #Habitat plot(density(DepthData[c("Labeotropheus_fuelleborni","Labeotropheus_trewavasae","Asprotilapia_leptura")]), xlim=c(-2,2), ylim=c(0,5), main="Habitat") for(i in 1:1000){ lines(density(as.numeric(HabitatSimList[[i]])), col="purple") } ##How often does real data overlap with the extreme snout depths #Single Peak CollateSPData<-matrix(unlist(lapply(SPSimList,range)),nrow = length(SPSimList), ncol = 2, byrow = T) sum(range(DepthData[c("Labeotropheus_fuelleborni","Labeotropheus_trewavasae","Asprotilapia_leptura")])[1]>CollateSPData[,2]) #Diet CollateDietData<-matrix(unlist(lapply(DietSimList,range)),nrow = length(DietSimList), ncol = 2, byrow = T) sum(range(DepthData[c("Labeotropheus_fuelleborni","Labeotropheus_trewavasae","Asprotilapia_leptura")])[1]>CollateDietData[,2]) #Feeding CollateFeedingData<-matrix(unlist(lapply(FeedingSimList,range)),nrow = length(FeedingSimList), ncol = 2, byrow = T) sum(range(DepthData[c("Labeotropheus_fuelleborni","Labeotropheus_trewavasae","Asprotilapia_leptura")])[1]>CollateFeedingData[,2]) #Habitat CollateHabitatData<-matrix(unlist(lapply(HabitatSimList,range)),nrow = length(HabitatSimList), ncol = 2, byrow = T) sum(range(DepthData[c("Labeotropheus_fuelleborni","Labeotropheus_trewavasae","Asprotilapia_leptura")])[1]>CollateHabitatData[,2]) #Summary min(DepthData[c("Labeotropheus_fuelleborni","Labeotropheus_trewavasae","Asprotilapia_leptura")]) max(CollateSPData[,2]) max(CollateDietData[,2]) max(CollateFeedingData[,2]) max(CollateHabitatData[,2]) max(CollateFeedingData[,2]) ###Simulate with Length A### ###W/ Labeotropheus_fuelleborni, Labeotropheus trewavasae, Asprotilapia leptura ###Parameters from evolultionary models without extreme individuals ##Diet #Extract diet data and perform simulation on full tree DietSimList<-rep(list(list()), 1000) for(i in 1:1000){ ThetaP<-setNames(MyOUwieModelsLengthAlphaDiet[[i]]$theta[,1],c("algae","large evasive","mollusks","small evasive")) DietSimList[[i]]<-multiOU(tree = SimmapTreesDiet[[i]], alpha = MyOUwieModelsLengthAlphaDiet[[i]]$solution[1,], sig2 = MyOUwieModelsLengthAlphaDiet[[i]]$solution[2,], theta = ThetaP) } ##Feeding #Extract feeding data and perform simulation FeedingSimList<-rep(list(list()), 1000) for(i in 1:1000){ ThetaP<-setNames(MyOUwieModelsLengthAlphaFeeding[[i]]$theta[,1],c("crushing","scraper","sifter","suction")) FeedingSimList[[i]]<-multiOU(tree = SimmapTreesFeeding[[i]], alpha = MyOUwieModelsLengthAlphaFeeding[[i]]$solution[1,], sig2 = MyOUwieModelsLengthAlphaFeeding[[i]]$solution[2,], theta = ThetaP) } ##Habitat #Extract habitat data and perform simulation HabitatSimList<-rep(list(list()), 1000) for(i in 1:1000){ ThetaP<-setNames(MyOUwieModelsLengthAlphaHabitat[[i]]$theta[,1],c("intermediate","open water","rocky","sandy")) HabitatSimList[[i]]<-multiOU(tree = SimmapTreesHabitat[[i]], alpha = MyOUwieModelsLengthAlphaHabitat[[i]]$solution[1,], sig2 = MyOUwieModelsLengthAlphaHabitat[[i]]$solution[2,], theta = ThetaP) } ##Single Peak #Extract single optima and perform simulation SPSimList<-rep(list(list()), 1000) for(i in 1:1000){ SPSimList[[i]]<-fastBM(tree = SampleDropTrees[[i]], sig2 = OUwieWOextremeOU[[i]]$solution[2], alpha = OUwieWOextremeOU[[i]]$solution[1], theta = OUwieWOextremeOU[[i]]$theta[1]) } #Plot LengthData<-setNames(Data$FlapLengthBDCorr,rownames(Data)) LengthData[c("Labeotropheus_fuelleborni","Labeotropheus_trewavasae","Asprotilapia_leptura")] par(mfrow=c(2,2), pty='s') #Single Peak plot(density(LengthData[c("Labeotropheus_fuelleborni","Labeotropheus_trewavasae","Asprotilapia_leptura")]), xlim=c(-5,5), ylim=c(0,2), main="Single Peak") for(i in 1:1000){ lines(density(as.numeric(SPSimList[[i]])), col="dark green") } #Diet plot(density(LengthData[c("Labeotropheus_fuelleborni","Labeotropheus_trewavasae","Asprotilapia_leptura")]), xlim=c(-5,5), ylim=c(0,2), main="Diet") for(i in 1:1000){ lines(density(as.numeric(DietSimList[[i]])), col="red") } #Feeding plot(density(LengthData[c("Labeotropheus_fuelleborni","Labeotropheus_trewavasae","Asprotilapia_leptura")]), xlim=c(-5,5), ylim=c(0,2), main="Feeding") for(i in 1:1000){ lines(density(as.numeric(FeedingSimList[[i]])), col="blue") } #Habitat plot(density(LengthData[c("Labeotropheus_fuelleborni","Labeotropheus_trewavasae","Asprotilapia_leptura")]), xlim=c(-5,5), ylim=c(0,2), main="Habitat") for(i in 1:1000){ lines(density(as.numeric(HabitatSimList[[i]])), col="purple") } ##How often does real data overlap with the extreme snout lengths #Single Peak CollateSPData<-matrix(unlist(lapply(SPSimList,range)),nrow = length(SPSimList), ncol = 2, byrow = T) sum(range(LengthData[c("Labeotropheus_fuelleborni","Labeotropheus_trewavasae","Asprotilapia_leptura")])[1]>CollateSPData[,2]) #Diet CollateDietData<-matrix(unlist(lapply(DietSimList,range)),nrow = length(DietSimList), ncol = 2, byrow = T) sum(range(LengthData[c("Labeotropheus_fuelleborni","Labeotropheus_trewavasae","Asprotilapia_leptura")])[1]>CollateDietData[,2]) #Feeding CollateFeedingData<-matrix(unlist(lapply(FeedingSimList,range)),nrow = length(FeedingSimList), ncol = 2, byrow = T) sum(range(LengthData[c("Labeotropheus_fuelleborni","Labeotropheus_trewavasae","Asprotilapia_leptura")])[1]>CollateFeedingData[,2]) #Habitat CollateHabitatData<-matrix(unlist(lapply(HabitatSimList,range)),nrow = length(HabitatSimList), ncol = 2, byrow = T) sum(range(LengthData[c("Labeotropheus_fuelleborni","Labeotropheus_trewavasae","Asprotilapia_leptura")])[1]>CollateHabitatData[,2]) #Summary min(LengthData[c("Labeotropheus_fuelleborni","Labeotropheus_trewavasae","Asprotilapia_leptura")]) max(CollateSPData[,2]) max(CollateDietData[,2]) max(CollateFeedingData[,2]) max(CollateHabitatData[,2]) max(CollateFeedingData[,2]) ###Length Simulations### setwd("/Users/Home/Desktop/Supp") source("Scripts/Ectodini_SOURCE.R") ##Remove Extreme length individuals #Tree and data setup Tr<-read.nexus("LTandLM_con.nex") Tre<-drop.tip(Tr, tip=c("Xenotilapia_caudafasciata", "Xenotilapia_longispinis")) EctoData<-read.csv("CichlidFemaleMean_20181126.csv", row.names = 1) Pruning<-treedata(Tre, EctoData) EctoTreePrune<-Pruning$phy EctoTreePrune<-ladderize(EctoTreePrune) EctoDataPrune<-Pruning$data EctoDataPrune<- EctoDataPrune[match(EctoTreePrune$tip.label, rownames(EctoDataPrune)),] EctoData<-as.data.frame(EctoDataPrune) EctoTree<-ladderize(EctoTreePrune,T) #For Cichlids as listed as factors for(i in 6:ncol(EctoData)){ EctoData[,i]<-as.numeric(as.character(EctoData[,i])) } #Evolutionary Models W/O Extreme Individuals# #Generate trees with stocastically mapped characters over a posterior distribution. ManyTrees<-read.nexus("~/Documents/LTandLM_Comb.trees") #Randomly subsample from a posterior distribution of trees SampleTrees<-sample(ManyTrees[10004:32004],1000) #Drop additional taxa from the trees - inc. depth exception species (w/o nasuta) SampleDropTrees<-lapply(SampleTrees,drop.tip,tip=c("Aulonocranus_dewindti", "Callochromis_pleurospilus", "Callochromis_stappersii", "Ectodus_descampsii", "Enantiopus_melanogenys", "Microdontochromis_rotundiventralis", "Oreochromis_tanganicae", "Xenotilapia_cf_bathyphila", "Xenotilapia_flavipinnis", "Xenotilapia_papilio_sunflower", "Xenotilapia_caudafasciata", "Xenotilapia_longispinis")) class(SampleDropTrees)<-"multiPhylo" save(SampleDropTrees, file = "SampleDropTrees_Alpha.rda") ###Diet### FlapDiet<-EctoData$Diet SimmapAlphaTreesDiet<-make.simmap(SampleDropTrees, FlapDiet, nsim=1, message=FALSE) save(SimmapAlphaTreesDiet, file = "Female_Diet_SIMMAPAlpha_OUwie.rda") MyOUwieModelsLengthAlphaDiet<-OUwieModels(SimmapAlphaTreesDiet, EctoData, "FlapLengthBDCorr", "Diet", "OUM") save(MyOUwieModelsLengthAlphaDiet, file = "Female_Diet_LengthAlpha_OUwie.rda") ###Habitat### FlapHabitat<-EctoData$Habitat SimmapAlphaTreesHabitat<-make.simmap(SampleDropTrees, FlapHabitat, nsim=1, message=FALSE) save(SimmapAlphaTreesHabitat, file = "Female_Habitat_SIMMAPAlpha_OUwie.rda") MyOUwieModelsLengthAlphaHabitat<-OUwieModels(SimmapAlphaTreesHabitat, EctoData, "FlapLengthBDCorr", "Habitat", "OUM") save(MyOUwieModelsLengthAlphaHabitat, file = "Female_Habitat_LengthAlpha_OUwie.rda") ###Feeding### FlapFeeding<-EctoData$Feeding SimmapAlphaTreesFeeding<-make.simmap(SampleDropTrees, FlapFeeding, nsim=1, message=FALSE) save(SimmapAlphaTreesFeeding, file = "Female_Feeding_SIMMAPAlpha_OUwie.rda") MyOUwieModelsLengthAlphaFeeding<-OUwieModels(SimmapAlphaTreesFeeding, EctoData, "FlapLengthBDCorr", "Feeding", "OUM") save(MyOUwieModelsLengthAlphaFeeding, file = "Female_Feeding_LengthAlpha_OUwie.rda") ###OUwie Single Peak/Single Rate set up### EctoDatas<-EctoData EctoDatas$Species<-rownames(EctoDatas) MyData<-EctoDatas[,c("Species", "Feeding", "FlapDepthBDCorr")] #Need to add a null two regime column for it to work MyData$Feeding<-as.factor(rep(x = c("SP","SR"), nrow(MyData)/2)) #Single peak OU OUwieWOextremeOU = mclapply(SampleDropTrees, function(x) { OUwie(x, MyData, "OU1") }) save(OUwieWOextremeOU, file = "SingleOU_Results_SimAlphaOUwie.rda") ###Simulate with Length B### ###W/ Xenotilapia_caudafasciata and Xenotilapia_longispinis ###Parameters from evolutionary models without extreme individuals ##Diet #Extract diet data and perform simulation on full tree DietSimList<-rep(list(list()), 1000) for(i in 1:1000){ ThetaP<-setNames(MyOUwieModelsLengthAlphaDiet[[i]]$theta[,1],c("algae","large evasive","mollusks","small evasive")) DietSimList[[i]]<-multiOU(tree = SimmapAlphaTreesDiet[[i]], alpha = MyOUwieModelsLengthAlphaDiet[[i]]$solution[1,], sig2 = MyOUwieModelsLengthAlphaDiet[[i]]$solution[2,], theta = ThetaP) } ##Feeding #Extract feeding data and perform simulation FeedingSimList<-rep(list(list()), 1000) for(i in 1:1000){ ThetaP<-setNames(MyOUwieModelsLengthAlphaFeeding[[i]]$theta[,1],c("crushing","scraper","sifter","suction")) FeedingSimList[[i]]<-multiOU(tree = SimmapAlphaTreesFeeding[[i]], alpha = MyOUwieModelsLengthAlphaFeeding[[i]]$solution[1,], sig2 = MyOUwieModelsLengthAlphaFeeding[[i]]$solution[2,], theta = ThetaP) } ##Habitat #Extract habitat data and perform simulation HabitatSimList<-rep(list(list()), 1000) for(i in 1:1000){ ThetaP<-setNames(MyOUwieModelsLengthAlphaHabitat[[i]]$theta[,1],c("intermediate","open water","rocky","sandy")) HabitatSimList[[i]]<-multiOU(tree = SimmapAlphaTreesHabitat[[i]], alpha = MyOUwieModelsLengthAlphaHabitat[[i]]$solution[1,], sig2 = MyOUwieModelsLengthAlphaHabitat[[i]]$solution[2,], theta = ThetaP) } ##Single Peak #Extract single optima and perform simulation SPSimList<-rep(list(list()), 1000) for(i in 1:1000){ SPSimList[[i]]<-fastBM(tree = SampleDropTrees[[i]], sig2 = OUwieWOextremeOU[[i]]$solution[2], alpha = OUwieWOextremeOU[[i]]$solution[1], theta = OUwieWOextremeOU[[i]]$theta[1]) } #Plot LengthData<-setNames(Data$FlapLengthBDCorr,rownames(Data)) LengthData[c("Xenotilapia_caudafasciata","Xenotilapia_longispinis")] par(mfrow=c(2,2), pty='s') #Single Peak plot(density(LengthData[c("Xenotilapia_caudafasciata","Xenotilapia_longispinis")]), xlim=c(-5,5), ylim=c(0,2), main="Single Peak") for(i in 1:1000){ lines(density(as.numeric(SPSimList[[i]])), col="dark green") } #Diet plot(density(LengthData[c("Xenotilapia_caudafasciata","Xenotilapia_longispinis")]), xlim=c(-5,5), ylim=c(0,2), main="Diet") for(i in 1:1000){ lines(density(as.numeric(DietSimList[[i]])), col="red") } #Feeding plot(density(LengthData[c("Xenotilapia_caudafasciata","Xenotilapia_longispinis")]), xlim=c(-5,5), ylim=c(0,2), main="Feeding") for(i in 1:1000){ lines(density(as.numeric(FeedingSimList[[i]])), col="blue") } #Habitat plot(density(LengthData[c("Xenotilapia_caudafasciata","Xenotilapia_longispinis")]), xlim=c(-5,5), ylim=c(0,2), main="Habitat") for(i in 1:1000){ lines(density(as.numeric(HabitatSimList[[i]])), col="purple") } ##How often does real data overlap with the extreme snout lengths #Single Peak CollateSPData<-matrix(unlist(lapply(SPSimList,range)),nrow = length(SPSimList), ncol = 2, byrow = T) sum(range(LengthData[c("Xenotilapia_caudafasciata","Xenotilapia_longispinis")])[1]>CollateSPData[,2]) #Diet CollateDietData<-matrix(unlist(lapply(DietSimList,range)),nrow = length(DietSimList), ncol = 2, byrow = T) sum(range(LengthData[c("Xenotilapia_caudafasciata","Xenotilapia_longispinis")])[1]>CollateDietData[,2]) #Feeding CollateFeedingData<-matrix(unlist(lapply(FeedingSimList,range)),nrow = length(FeedingSimList), ncol = 2, byrow = T) sum(range(LengthData[c("Xenotilapia_caudafasciata","Xenotilapia_longispinis")])[1]>CollateFeedingData[,2]) #Habitat CollateHabitatData<-matrix(unlist(lapply(HabitatSimList,range)),nrow = length(HabitatSimList), ncol = 2, byrow = T) sum(range(LengthData[c("Xenotilapia_caudafasciata","Xenotilapia_longispinis")])[1]>CollateHabitatData[,2]) #Summary min(LengthData[c("Xenotilapia_caudafasciata","Xenotilapia_longispinis")]) max(CollateSPData[,2]) max(CollateDietData[,2]) max(CollateFeedingData[,2]) max(CollateHabitatData[,2]) max(CollateFeedingData[,2]) ####Bias in Evolutionary Models#### setwd("/Users/Home/Desktop/Supp") source("Scripts/Ectodini_SOURCE.R") #Analysis with residual corrections #Tree and data setup Files<-MakeFiles() #Your tree and data is stored as 'Tree' and 'Data' respectively. Tree<-Files[[1]] Data<-Files[[2]] #For Cichlids as listed as factors for(i in 6:ncol(Data)){ Data[,i]<-as.numeric(as.character(Data[,i])) } #1000 BM sims run over 1000*1 post-distro SIMMAP trees #Generate trees with stocastically mapped characters over a posterior distribution. ManyTrees<-read.nexus("LTandLM_Comb.trees") #Randomly subsample from a posterior distribution of trees SampleTrees<-sample(ManyTrees[10004:32004],1000) #Drop additional taxa from the trees SampleDropTrees<-lapply(SampleTrees,drop.tip,tip=c("Aulonocranus_dewindti", "Callochromis_pleurospilus", "Callochromis_stappersii", "Ectodus_descampsii", "Enantiopus_melanogenys", "Microdontochromis_rotundiventralis", "Oreochromis_tanganicae", "Xenotilapia_cf_bathyphila", "Xenotilapia_flavipinnis", "Xenotilapia_papilio_sunflower")) class(SampleDropTrees)<-"multiPhylo" save(SampleDropTrees, file = "Female_SampleDropTrees_SimOUwie.rda") ##Diet## FlapDiet<-Data$Diet SimmapTreesDiet<-make.simmap(SampleDropTrees, FlapDiet, nsim=1, message=FALSE) save(SimmapTreesDiet, file = "Female_Diet_SIMMAP_SimOUwie.rda") ##Habitat## FlapHabitat<-Data$Habitat SimmapTreesHabitat<-make.simmap(SampleDropTrees, FlapHabitat, nsim=1, message=FALSE) save(SimmapTreesHabitat, file = "Female_Habitat_SIMMAP_SimOUwie.rda") ##Feeding## FlapFeeding<-Data$Feeding SimmapTreesFeeding<-make.simmap(SampleDropTrees, FlapFeeding, nsim=1, message=FALSE) save(SimmapTreesFeeding, file = "Female_Feeding_SIMMAP_SimOUwie.rda") ###Simulate Trait Data### BMSimList<-rep(list(list()), 1000) for(i in 1:1000){ BMSims<-fastBM(Tree) BMSims<-cbind.data.frame(rownames(Data), as.character(Data$Diet), as.character(Data$Habitat), as.character(Data$Feeding), BMSims) colnames(BMSims)<-c("PhyloName", "Diet", "Habitat", "Feeding", "BMSims") BMSimList[[i]]<-BMSims } save(BMSimList, file = "Female_SimData_SimOUwie.rda") ###Diet### OUwieSimDiet<-NULL OUwieSimDiet<-as.list(OUwieSimDiet) for(i in 1:1000){ OUwieSimDiet[[i]]<-OUwie(SimmapTreesDiet[[i]], BMSimList[[i]][,c(1,2,5)], model="OUM", simmap.tree=T) } save(OUwieSimDiet, file = "Female_Diet_Results_SimOUwie.rda") ###Habitat### OUwieSimHabitat<-NULL OUwieSimHabitat <-as.list(OUwieSimHabitat) for(i in 1:1000){ OUwieSimHabitat[[i]]<-OUwie(SimmapTreesHabitat[[i]], BMSimList[[i]][,c(1,3,5)], model="OUM", simmap.tree=T) } save(OUwieSimHabitat, file = "Female_Habitat_Results_SimOUwie.rda") ###Feeding### OUwieSimFeeding<-NULL OUwieSimFeeding<-as.list(OUwieSimFeeding) for(i in 1:1000){ OUwieSimFeeding[[i]]<-OUwie(SimmapTreesFeeding[[i]], BMSimList[[i]][,c(1,4,5)], model="OUM", simmap.tree=T) } save(OUwieSimFeeding, file = "Female_Feeding_Results_SimOUwie.rda") ###Saving Sim### #Diet OUwieResults<-matrix(NA,length(OUwieSimDiet),2) for(i in 1:length(OUwieSimDiet)){ OUwieResults[i,1] <-OUwieSimDiet[[i]]$loglik OUwieResults[i,2] <-OUwieSimDiet[[i]]$AICc } colnames(OUwieResults)<-c("log likelihood", "AICc") write.csv(OUwieResults, "Diet_Sim_OUwieResults.csv", row.names=F) #Habitat OUwieResults<-matrix(NA,length(OUwieSimHabitat),2) for(i in 1:length(OUwieSimHabitat)){ OUwieResults[i,1] <-OUwieSimHabitat[[i]]$loglik OUwieResults[i,2] <-OUwieSimHabitat[[i]]$AICc } colnames(OUwieResults)<-c("log likelihood", "AICc") write.csv(OUwieResults, "Habitat_Sim_OUwieResults.csv", row.names=F) #Feeding OUwieResults<-matrix(NA,length(OUwieSimFeeding),2) for(i in 1:length(OUwieSimFeeding)){ OUwieResults[i,1] <-OUwieSimFeeding[[i]]$loglik OUwieResults[i,2] <-OUwieSimFeeding[[i]]$AICc } colnames(OUwieResults)<-c("log likelihood", "AICc") write.csv(OUwieResults, "Feeding_Sim_OUwieResults.csv", row.names=F) ###OUCH sims### #Sim traits OUCHSimData<-NULL OUCHSimData <-as.list(OUCHSimData) for(i in 1:1000){ OUCHSimData[[i]]<-OuchModels(SampleDropTrees[[i]], BMSimList[[i]], "BMSims", "Diet") } save(OUCHSimData, file = "SinglePeakOU_Results_SimOUCH.rda") #Saving OUCHSimSave<-matrix(NA,length(OUCHSimData),4) for(i in 1:length(OUCHSimData)){ OUCHSimSave[i,1] <-OUCHSimData[[i]][1] OUCHSimSave[i,2] <-OUCHSimData[[i]][4] OUCHSimSave[i,3] <-OUCHSimData[[i]][2] OUCHSimSave[i,4] <-OUCHSimData[[i]][5] } colnames(OUCHSimSave)<-c("log likelihood BM", "AICc BM", "log likelihood SP OU", "AICc SP OU") write.csv(OUCHSimSave, "Single_Sim_OUCHResults.csv", row.names=F) ###Calcualte dAICc### setwd("/Users/Home/Desktop/Supp") #Combine all the AICc data into a single spreadsheet SimComb<-read.csv("AllAICcSimCombine.csv") SimCombSave<-matrix(NA,length(OUCHSimData), 5) for(i in 1:nrow(SimComb)){ AICcDataStore <- akaike.weights(SimComb[i,]) SimCombSave[i,] <- AICcDataStore$deltaAIC } write.csv(SimCombSave, "AlldAICcSimCombine.csv", row.names=F) ##Is there bias in the selection of complex models?## SimCombSave<-read.csv("AlldAICcSimCombine.csv") ##Frequency of selecting different models when traits generated under BM #When BM is the best model colSums(SimCombSave==0) #BM SPOU Diet Feeding Habitat #675 158 56 46 65 #When BM is <2/4 units from any other model #BM is selected, or, other model is not significantly favored over BM (sum(SimCombSave$BMdAICc<2)/1000)*100 #81.2% (sum(SimCombSave$BMdAICc<4)/1000)*100 #90.5% #Of the 46 instances where feeding is selected, how often is it <2 units from BM sum(SimCombSave$BMdAICc[SimCombSave$FeedingdAICc==0]<2) #16, 46-16=30 instances where feeding is outright selected (sum(SimCombSave$BMdAICc[SimCombSave$FeedingdAICc==0]<2)/46)*100 #34.78% ####Ancestral State Reconstruction with Rphylopars#### setwd("/Users/Home/Desktop/Supp") source("Scripts/Ectodini_SOURCE.R") #Analysis with residual corrections #Tree and data setup Files<-MakeFiles() #Your tree and data is stored as 'Tree' and 'Data' respectively. Tree<-Files[[1]] Data<-Files[[2]] #For Cichlids as listed as factors for(i in 6:ncol(Data)){ Data[,i]<-as.numeric(as.character(Data[,i])) } #Extract required data columns Dataz<-Data Dataz$Species<-rownames(Data) AncCalc<-Dataz[,c(1,9,10)] colnames(AncCalc)<-c("species", "FlapDepthBDCorr", "FlapLengthBDCorr") #Make tree ultrametric (it is already very close) NewTree<-force.ultrametric(Tree) #Calculate ancestral states based on a branch-variable OU parameter p_mvOU <- phylopars(trait_data = AncCalc, tree = NewTree, model = "mvOU") #Set-up data for contour map plotting #Depth DepthBD<-Data$FlapDepthBDCorr names(DepthBD)<-rownames(Data) DepthObj<-contMap(tree = Tree, x = DepthBD, method="user", anc.states=p_mvOU$anc_recon[34:65,][,1]) plot(setMap(DepthObj, invert=TRUE)) #Length LengBD<-Data$FlapLengthBDCorr names(LengBD)<-rownames(Data) LengthObj<-contMap(tree = Tree, x = LengBD, method="user", anc.states=p_mvOU$anc_recon[34:65,][,2]) plot(setMap(LengthObj, invert=TRUE)) ####Phylomorphospace#### par(pty='s') tip.cols<-rep("black",33) names(tip.cols)<-Tree$tip.label cols<-c(tip.cols[Tree$tip.label],rep("gray",Tree$Nnode)) names(cols)<-1:(length(Tree$tip)+Tree$Nnode) phylomorphospaceOUAnc(tree = Tree, X = Data[,c(10,9)], AncData = p_mvOU$anc_recon[34:65,c(2,1)], control=list(col.node=cols)) ####Plotting Extras#### #Plot convex hull, length x depth plot par(pty='s') MySym<-c(19, 17) plot(Data$FlapLengthBDCorr, Data$FlapDepthBDCorr, col=unclass(Data$Lake), pch=MySym[unclass(Data$Lake)]) LM<-Data[which(Data$Lake == "Malawi"),] LT<-Data[which(Data$Lake == "Tanganyika"),] ChullA<-chull(LM[,12:13]) ChullA<-c(ChullA, ChullA[1]) lines(LM[,12:13][ChullA, ], lty=1, col="black", lwd=2) ChullB<-chull(LT[,12:13]) ChullB<-c(ChullB, ChullB[1]) lines(LT[,12:13][ChullB, ], lty=2, col="red", lwd=2) #Plotting diet, feeding, and habitat plot(Data$FlapLengthBDCorr, Data$FlapDepthBDCorr, col=unclass(Data$Diet), pch=MySym[unclass(Data$Lake)]) plot(Data$FlapLengthBDCorr, Data$FlapDepthBDCorr, col=unclass(Data$Feeding), pch=MySym[unclass(Data$Lake)]) plot(Data$FlapLengthBDCorr, Data$FlapDepthBDCorr, col=unclass(Data$Habitat), pch=MySym[unclass(Data$Lake)])