library(geiger) library(phytools) R.Version() ######################################################################################## ############# Ungulata analyses ######################################################## ######################################################################################## #load phylogeny phy <- read.nexus("ungulate_tree.nex") #ungulate phylogeny plot(phy, cex = 0.2) str(phy) #load data files dta <- read.delim("Bovid_cervid_data.txt",sep="\t") #load data row.names(dta)<-dta[,1] names(dta) #transform data to the same scale - here mm for length measures and g for mass measures dta$BodyM <- dta$BodyM*1000 dta$WeaponM <- dta$WeaponM*10 dta$SpHead <- dta$SpHead/1000 dta$SpMidpiece <- dta$SpMidpiece/1000 dta$SpTail <- dta$SpTail/1000 dta$Soma <- dta$BM-dta$CTM dta$MuzzleM <- dta$MuzzleM*10 #select numerical variables to log-transform them nums <- sapply(dta, is.numeric) dta[,nums==T] <- log10(dta[,nums==T]) #Remove the youngest branch in the mass (G. cuvieri and leptoceros) and the youngest in the length (C. taurinus and C. gnou) #Removing these species allows the data to fit with model assumptions (see Methods in the paper) dta <- dta[-c(52,53,61,64),] #match names between dataset and phylogeny to ensure that all species of phylogeny are listed in dataset phyl <- ladderize(phy) sp <- data.frame(Species = phyl$tip.label) #Grab species dat1 <- merge(dta, sp, by.x = "Phylo", by.y = "Species", all.y = T) #Appends species from new phylogeny to old one rownames(dat1) <- dat1$Phylo names(dat1) name.check(phyl, dat1) ### Analyses of evolutionary rates #select one of the models and follow the subsequent steps. Repeat from here when assessing another model. #Model 1 - All Head, Midpiece, Tail, Muzzle, Weapon #discard<-which(is.na(dat1$WeaponM)==T | is.na(dat1$MuzzleM)==T | is.na(dat1$SpHead)==T | is.na(dat1$SpMidpiece)==T | is.na(dat1$SpTail)==T) # choose variables to be used #Model 1a - Only Cervids #discard<-which(is.na(dat1$WeaponM)==T | is.na(dat1$MuzzleM)==T | is.na(dat1$SpHead)==T | is.na(dat1$SpMidpiece)==T | is.na(dat1$SpTail)==T | dat1$Family!="Cervidae") #Model 1b - Only Bovids #discard<-which(is.na(dat1$WeaponM)==T | is.na(dat1$MuzzleM)==T | is.na(dat1$SpHead)==T | is.na(dat1$SpMidpiece)==T | is.na(dat1$SpTail)==T | dat1$Family!="Bovidae") #Model 2 - All #discard<-which(is.na(dat1$CTM)==T | is.na(dat1$BodyM)==T ) # choose variables to be used #Model 2a - Cervids #discard<-which(is.na(dat1$CTM)==T | is.na(dat1$BodyM)==T | dat1$Family!="Cervidae") # choose variables to be used #Model 2b - Bovids #discard<-which(is.na(dat1$CTM)==T | is.na(dat1$BodyM)==T | dat1$Family!="Bovidae") # choose variables to be used discard.sp<-as.character(dat1[discard,"Phylo"]) phylo2<-drop.tip(phyl,discard.sp) # drop species with NA from phylogeny par(mfrow=c(1,1),pch=19,ps=10,cex=1.3,cex.lab=1.5,cex.axis=1) phylo2 <- ladderize(phylo2 ) plot(phylo2,cex=0.5,edge.width=0.8,use.edge.length=F,no.margin=T,label.offset=0.02) phylomat1<-vcv.phylo(phylo2) data<-dat1[-discard,] str(data) phylo2$tip.label newtree <- phylo2 plot(newtree, cex=0.5) name.check(newtree, data) ############################################ ######## Model of evolution test ########### ############################################ names(data) trait <- data$BodyM #name corresponds to column of data. Change the number to assess a new trait. #trait <- data$WeaponM #name corresponds to column of data. Change the number to assess a new trait. #trait <- data$TSL #name corresponds to column of data. Change the number to assess a new trait. #trait <- data$SpHead #name corresponds to column of data. Change the number to assess a new trait. #trait <- data$SpMidpiece #name corresponds to column of data. Change the number to assess a new trait. #trait <- data$SpTail #name corresponds to column of data. Change the number to assess a new trait. #trait <- data$CTM #name corresponds to column of data. Change the number to assess a new trait. #trait <- data$MuzzleM #name corresponds to column of data. Change the number to assess a new trait. #trait <- data$Soma #name corresponds to column of data. Change the number to assess a new trait. length(trait) names(trait)<- row.names(data) #you can look at the distribtion of the data on the tree here. plotTree.wBars(newtree,trait) plot(newtree, cex = 0.5) ##fitting different models #note, check to make sure the models fit (i.e. no error) BMfit <- fitContinuous(newtree, trait, model="BM") #fits Brownian motion model BMfit OUfit <- fitContinuous(newtree, trait, model="OU", bounds=list(alpha=c(min=exp(-40),max=exp(10)))) #fits Ornstein-Uhlenbeck model, if fitting is an issue try: , bounds=list(alpha=c(min=exp(-40),max=exp(10)))) OUfit EBfit <- fitContinuous(newtree, trait, model="EB", bounds=list(alpha=c(min=exp(-50),max=exp(10)))) #fits Early-burst model, if fitting is an issue try: , bounds=list(a=c(min=exp(-40),max=exp(7)))) EBfit #BM summary: parameter estimates, lnL, AICc BMfit$opt$sigsq #the Brownian rate parameter sigma squared BMfit$opt$lnL BMfit$opt$aicc #OU summary: parameter estimates, lnL, AICc OUfit$opt$alpha #alpha value of the model indicating the 'evolutionary pull' OUfit$opt$lnL OUfit$opt$aicc #EB summary: parameter estimates, lnL, AICc EBfit$opt$a #a is the rate change parameter EBfit$opt$lnL EBfit$opt$aicc ## Compare the models using AIC aic_bm <- BMfit$opt$aicc aic_ou <- OUfit$opt$aicc aic_eb <- EBfit$opt$aicc ## Get the AIC weights in a named vector aic_all <- c(aic_bm, aic_ou, aic_eb) names(aic_all) <- c("BM", "OU", "EB") aic_comp <- aicw(aic_all) tableBM <- as.numeric(c(BMfit$opt$sigsq, BMfit$opt$lnL, BMfit$opt$aicc, aic_comp$delta[1], aic_comp$w[1]), options(scipen = 999)) tableOU <- as.numeric(c(OUfit$opt$alpha, OUfit$opt$lnL, OUfit$opt$aicc, aic_comp$delta[2], aic_comp$w[2]), options(scipen = 999)) tableEB <- as.numeric(c(EBfit$opt$a, EBfit$opt$lnL, EBfit$opt$aicc, aic_comp$delta[3], aic_comp$w[3]), options(scipen = 0)) tableBM tableOU tableEB ############################################ ######## Rate of evolution tests ########### ############################################ names(data) #Model 1 #datmatrix <- data[,c("SpHead", "SpMidpiece", "SpTail", "MuzzleM", "WeaponM")] #datmatrix <- data[,c("SpHead", "SpMidpiece")] #datmatrix <- data[,c("SpHead", "SpTail")] #datmatrix <- data[,c("SpHead","MuzzleM")] #datmatrix <- data[,c("SpHead", "WeaponM")] #datmatrix <- data[,c("SpMidpiece", "SpTail")] #datmatrix <- data[,c("SpMidpiece", "MuzzleM")] #datmatrix <- data[,c("SpMidpiece", "WeaponM")] #datmatrix <- data[,c("SpTail", "MuzzleM")] #datmatrix <- data[,c("SpTail", "WeaponM")] #datmatrix <- data[,c("MuzzleM", "WeaponM")] #Model 2 #datmatrix <- data[,c("CTM", "BodyM")] # If significant effect. For posthoc redo analysis with only two traits, in all combinations x <- as.matrix(datmatrix) source("RATES_TEST_SOURCE_Nelder-Mead.R") #note the source script is modified from the script provided in Appendix 2 of Adams 2013, Syst. Biol. 62, 181–192. result <- CompareRates.multTrait(newtree,x,TraitCov=F) result #####plotting and calculating 95% CIs #calculating mean and 95%CIs library(numDeriv) #for hessian #x<-as.matrix(datmatrix) # is already defined above; change these values to whatever you want to compare in your analyses N<-nrow(x) p<-ncol(x) C<-vcv.phylo(newtree) C<-C[rownames(x),rownames(x)] #Fit Rate matrix for all traits a.obs<-colSums(solve(C))%*%x/sum(solve(C)) D<-matrix(0,N*p,p) for(i in 1:(N*p)) for(j in 1:p) if((j-1)*N