# lightweight version of ace(...,method="ML") for continuous traits # written by Liam J. Revell 2011 anc.ML<-function(tree,x,maxit=2000){ # function returns the log-likelihood likelihood<-function(par,C,invC,detC,x){ sig2<-par[1]; a<-par[2]; y<-par[3:length(par)] z<-c(x,y)-a logLik<--z%*%invC%*%z/(2*sig2)-nrow(C)*log(2*pi)/2-nrow(C)*log(sig2)/2-detC/2 return(-logLik) } # compute C C<-vcvPhylo(tree) invC<-solve(C) detC<-determinant(C,log=T)$modulus[1] x[rownames(C)[1:length(tree$tip)]]->x # assign starting values y<-runif(n=tree$Nnode-1)*(max(x)-min(x))+min(x) a<-runif(n=1)*(max(x)-min(x))+min(x) sig2<-mean(pic(x,tree)^2)*runif(n=1,max=2) # optimize res<-optim(c(sig2,a,y),fn=likelihood,C=C,invC=invC,detC=detC,x=x,method="L-BFGS-B",lower=c(1e-10,rep(-Inf,tree$Nnode)),control=list(maxit=maxit)) # return result states<-res$par[2:length(res$par)] names(states)<-c(length(tree$tip)+1,rownames(C)[(length(tree$tip)+1):nrow(C)]) return(list(sig2=res$par[1],ace=states,logLik=-res$value,convergence=res$convergence)) }