# six.sd <- data.frame(log.size, log.size.se, log.rw, log.rw.se, log.color, log.col.se) # six.sd$binom <- names(log.size) # write.csv(six.sd, 'lib.spp.6sd.csv') ##packages library(phytools) library(car) library(MASS) library(MuMIn) library(geiger) ## load in the data and phylogeny six.dat <- read.csv('lib.spp.6sd.csv') # load data rownames(six.dat) <- six.dat$binom phylo <- read.tree('lib.na.phylo.tre') dim(six.dat) # check to make sure the data loaded in okay. Should output 72 7 summary(phylo) # check to make sure the phylogeny loaded in okay. Should be 72 tips and 71 nodes. max should = 120.317071 ############### patterns of diversification ###### patterns of size diversification log.size <- six.dat$log.size log.size.se <- six.dat$log.size.se names(log.size) <- six.dat$binom names(log.size.se) <- six.dat$binom size.BM <- fitContinuous(phylo, log.size, SE = log.size.se, model = 'BM', control = list(hessian = TRUE)) # brownian motion size.EB <- fitContinuous(phylo, log.size, SE = log.size.se, model = 'EB', bounds = list(a = c(min = -2, max = 0.5))) # early burst size.white <- fitContinuous(phylo, log.size, SE = log.size.se, model = 'white') # white noise size.delta <- fitContinuous(phylo, log.size, SE = log.size.se, model = 'delta') # delta model size.OU <- fitContinuous(phylo, dat = log.size, log.size.se, model = 'OU') # OU model aic.size <- c(size.BM$opt$aic, size.EB$opt$aic, size.white$opt$aic, size.delta$opt$aic, size.OU$opt$aic) # save the AIC values from the three models names(aic.size) <- c('Brownian Motion', 'early burst', 'white', 'Delta', 'OU') aic.size # Brownian Motion early burst white Delta OU # -46.02120 -44.21265 20.47636 -44.22612 -44.21265 ###### patterns of rw diversification log.rw <- six.dat$log.rw log.rw.se <- six.dat$log.rw.se names(log.rw) <- six.dat$binom names(log.rw.se) <- six.dat$binom rw.BM <- fitContinuous(phylo, log.rw, SE = log.rw.se, model = 'BM') # brownian motion model rw.EB <- fitContinuous(phylo, log.rw, SE = log.rw.se, model = 'EB', control = list(hessian = TRUE)) # early burst rw.white <- fitContinuous(phylo, log.rw, SE = log.rw.se, model = 'white') # white rw.delta <- fitContinuous(phylo, log.rw, SE = log.rw.se, model = 'delta') # delta model rw.OU <- fitContinuous(phylo, log.rw, SE = log.rw.se, model = 'OU') #OU aic.rw <- c(rw.BM$opt$aic, rw.EB$opt$aic, rw.white$opt$aic, rw.delta$opt$aic, rw.OU$opt$aic) names(aic.rw) <- c('Brownian Motion', 'early burst', 'white', 'Delta', 'OU') aic.rw # Brownian Motion early burst white Delta OU # -161.3248 -162.1947 -126.2558 -161.5826 -159.3603 ###### patterns of color diversification log.color <- six.dat$log.color log.col.se <- six.dat$log.col.se names(log.color) <- six.dat$binom names(log.col.se) <- six.dat$binom color.BM <- fitContinuous(phylo, log.color, SE = log.col.se, model = 'BM', control = list(hessian = TRUE)) color.EB <- fitContinuous(phylo, log.color, SE = log.col.se, model = 'EB', bounds = list(a = c(min = -1, max = 0))) color.white <-fitContinuous(phylo, log.color, SE = log.col.se, model = 'white') color.delta <- fitContinuous(phylo, log.color, SE = log.col.se, model = 'delta', bounds = list(delta = c(min = exp(-500), max = 15))) color.OU <- fitContinuous(phy = phylo, dat = log.color, model = "OU", SE = log.col.se) aic.color <- c(color.BM$opt$aic, color.EB$opt$aic, color.white$opt$aic, color.delta$opt$aic, color.OU$opt$aic) names(aic.color) <- c('Brownian Motion', 'early burst', 'white', 'Delta', 'OU') aic.color # Brownian Motion early burst white Delta OU # 356.9669 358.9669 340.1961 335.2425 334.8998 #### Next, compare rates of phenotypic evolution using the Adams (2013) likelihood ratio approach. Need to run his custom's scripts first. ########## start here CompareRates.multTrait <- function(phy, x, TraitCov=T, ms.err=NULL, ms.cov=NULL){ #Compares LLik of R-matrix vs. LLik of R-matrix with constrained diagonal #TraitCov = TRUE assumes covariation among traits (default) #ms.err allows the incorporation of within-species measurement error. Input is a matrix of species (rows) by within-species variation for each trait (columns). #ms.cov allows the incorporation of within-species covariation between traits. Input is a matrix of species (rows) by within-species covariation for each pair of traits (columns). These must be provided in a specific order, beginning with covariation between trait 1 and the rest, then trait 2 and the rest, etc. For instance, for 4 traits, the columns are: cov_12, cov_13, cov_14, cov_23, cov_24 cov_34. #Some calculations adapted from 'evol.vcv' in phytools (Revell, 2012) library(MASS) x<-as.matrix(x) N<-nrow(x) p<-ncol(x) C<-vcv.phylo(phy) C<-C[rownames(x),rownames(x)] if (is.matrix(ms.err)){ ms.err<-as.matrix(ms.err[rownames(x),])} if (is.matrix(ms.cov)){ ms.cov<-as.matrix(ms.cov[rownames(x),])} #Cholesky decomposition function for diagonal-constrained VCV build.chol<-function(b){ c.mat<-matrix(0,nrow=p,ncol=p) c.mat[lower.tri(c.mat)] <- b[-1] c.mat[p,p]<-exp(b[1]) c.mat[1,1]<-sqrt(sum((c.mat[p,])^2)) if(p>2){ for (i in 2:(p-1)){ c.mat[i,i]<-ifelse( (c.mat[1,1]^2-sum((c.mat[i,])^2) )>0, sqrt(c.mat[1,1]^2-sum((c.mat[i,])^2)), 0) }} return(c.mat) } #Fit Rate matrix for all traits: follows code of L. Revell (evol.vcv) 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)*NF)'[1] Pvalue.pic } ######### end here ### to run this analysis, need to create matrices that includes our traits and the standard errors of the traits adult.columns <- data.frame(log.size, log.rw, log.color) # combine our three named vectors into a "data frame" adult.traits <- as.matrix(adult.columns) # then turn that dataframe into a matrix ses.frame <- data.frame(log.size.se, log.rw.se, log.col.se) # combine our three named vectors for the standard errors into a data frame adult.ses <- as.matrix(ses.frame) # then turn that dataframe into a matrix ## compare the rates with likelihood ratio test (Adams 2013) rate.test <- CompareRates.multTrait(phylo, adult.traits, TraitCov = T, ms.err = adult.ses, ms.cov = NULL) rate.test # $Robs # log.size log.rw log.color # log.size 5.957540e-04 -6.129719e-05 -0.001393426 # log.rw -6.129719e-05 1.532961e-04 0.001512303 # log.color -1.393426e-03 1.512303e-03 0.175182638 # $Rconstrained # [,1] [,2] [,3] # [1,] 1.000307 0.000000e+00 0.000000e+00 # [2,] 0.000000 1.000307e+00 1.533196e-05 # [3,] 0.000000 1.533196e-05 1.000307e+00 # $Lobs # [1] -128.8147 # $Lconstrained # [1] -626.2404 # $LRTest # [1] 994.8513 # $Prob # [1] 9.349336e-217 # $AICc.obs # [1] 269.6294 # $AICc.constrained # [1] 1260.481 # $optimmessage # [1] "Optimization has converged."