#Function to compare phylogenetic evolutionary rates among traits using likelihood #Dean C. Adams #From Adams, D.C. 2013. Comparing evolutionary rates for different phenotypic traits on a phylogeny using #likelihood. Systematic Biology. 62:181-192. # 6 April, 2013: Corrected null Rate matrix for TraitCov=F to be diagonal matrix 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)*N