# this function computes a phylogenetic reduced major axis (RMA) regression # written by Liam Revell 2010,2011,2012 phyl.RMA<-function(x,y,tree,method="BM",lambda=NULL,fixed=FALSE,h0=1.0){ x<-x[tree$tip.label]; y<-y[tree$tip.label] # bind the x & y into columns X<-cbind(x,y) # function to compute C with lambda lambda.transform<-function(C,lambda){ if(lambda==1) return(C) else { V<-diag(diag(C)) C<-C-V C.lambda<-(V+lambda*C) return(C.lambda) } } # function to compute phylogenetic VCV using joint Pagel's lambda phyl.vcv<-function(X,phy,lambda){ C<-vcv(phy) C<-lambda.transform(C,lambda) one<-matrix(1,nrow(C),1) a<-solve(t(one)%*%solve(C)%*%one)%*%(t(one)%*%solve(C)%*%X) V<-t(X-one%*%a)%*%solve(C)%*%(X-one%*%a)/nrow(C) return(list(R=V,alpha=a)) } # likelihood function likelihood<-function(theta,X,phy){ C<-vcv(phy) C<-lambda.transform(C,theta) # compute R, conditioned on lambda temp<-phyl.vcv(X,phy,theta); R<-temp$R; a<-temp$alpha; rm(temp) # prep n<-nrow(X); m<-ncol(X); D<-matrix(0,n*m,m) for(i in 1:(n*m)) for(j in 1:m) if((j-1)*n