# function is simplified version of evol.vcv # written by Liam J. Revell 2011, 2012, 2013 evolvcv.lite<-function(tree,X,maxit=2000,tol=1e-10){ # model 1: common variances & correlation lik1<-function(theta,C,D,y){ v<-theta[1:2]; r<-theta[3] R<-matrix(c(v[1],r*sqrt(v[1]*v[2]),r*sqrt(v[1]*v[2]),v[2]),2,2) V<-kronecker(R,C) a<-solve(t(D)%*%solve(V)%*%D)%*%(t(D)%*%solve(V)%*%y) logL<--t(y-D%*%a)%*%solve(V)%*%(y-D%*%a)/2-n*m*log(2*pi)/2-determinant(V)$modulus[1]/2 return(-logL) } # model 2: different variances, same correlation lik2<-function(theta,C,D,y){ p<-(length(theta)-1)/2 v<-matrix(theta[1:(2*p)],p,2,byrow=T); r<-theta[length(theta)] R<-list() for(i in 1:p) R[[i]]<-matrix(c(v[i,1],r*sqrt(v[i,1]*v[i,2]),r*sqrt(v[i,1]*v[i,2]),v[i,2]),2,2) V<-matrix(0,length(y),length(y)) for(i in 1:p) V<-V+kronecker(R[[i]],C[[i]]) a<-solve(t(D)%*%solve(V)%*%D)%*%(t(D)%*%solve(V)%*%y) logL<--t(y-D%*%a)%*%solve(V)%*%(y-D%*%a)/2-n*m*log(2*pi)/2-determinant(V)$modulus[1]/2 return(-logL) } # model 3: same variances, different correlation lik3<-function(theta,C,D,y){ p<-length(theta)-2 v<-theta[1:2]; r<-theta[3:length(theta)] R<-list() for(i in 1:p) R[[i]]<-matrix(c(v[1],r[i]*sqrt(v[1]*v[2]),r[i]*sqrt(v[1]*v[2]),v[2]),2,2) V<-matrix(0,length(y),length(y)) for(i in 1:p) V<-V+kronecker(R[[i]],C[[i]]) a<-solve(t(D)%*%solve(V)%*%D)%*%(t(D)%*%solve(V)%*%y) logL<--t(y-D%*%a)%*%solve(V)%*%(y-D%*%a)/2-n*m*log(2*pi)/2-determinant(V)$modulus[1]/2 return(-logL) } # model 4: everything different lik4<-function(theta,C,D,y){ p<-length(theta)/3 v<-matrix(theta[1:(2*p)],p,2,byrow=T); r<-theta[(2*p+1):length(theta)] R<-list() for(i in 1:p) R[[i]]<-matrix(c(v[i,1],r[i]*sqrt(v[i,1]*v[i,2]),r[i]*sqrt(v[i,1]*v[i,2]),v[i,2]),2,2) V<-matrix(0,length(y),length(y)) for(i in 1:p) V<-V+kronecker(R[[i]],C[[i]]) a<-solve(t(D)%*%solve(V)%*%D)%*%(t(D)%*%solve(V)%*%y) logL<--t(y-D%*%a)%*%solve(V)%*%(y-D%*%a)/2-n*m*log(2*pi)/2-determinant(V)$modulus[1]/2 return(-logL) } # done internal functions # bookkeeping X<-as.matrix(X) X<-X[tree$tip.label,] n<-nrow(X) # number of species m<-ncol(X) # number of traits if(m!=2) stop("number of traits must equal 2") p<-ncol(tree$mapped.edge) # number of states D<-matrix(0,n*m,m) for(i in 1:(n*m)) for(j in 1:m) if((j-1)*n