user supplied internal node states in contMap

R/add.species.to.genus.R

@@ -43,8 +43,8 @@ add.species.to.genus<-function(tree,species,genus=NULL,where=c("root","random"))
 genus.to.species.tree<-function(tree,species){
 	N<-Ntip(tree)
 	genera<-tree$tip.label
-	for(i in 1:N){
 	species<-gsub(" ","_",species)
+	for(i in 1:N){
 		jj<-grep(paste(genera[i],"_",sep=""),species)
 		nn<-which(tree$tip.label==genera[i])
 		if(length(jj)>1){

R/anc.ML.R

@@ -63,21 +63,24 @@ anc.OU<-function(tree,x,maxit=2000,...){
 	else tol<-1e-8
 	if(hasArg(trace)) trace<-list(...)$trace
 	else trace<-FALSE
+	if(hasArg(a.init)) a.init<-list(...)$a.init
+	else a.init<-2*log(2)/max(nodeHeights(tree))
 	likOU<-function(par,tree,x,trace){
 		sig2<-par[1]
 		alpha<-par[2]
 		a0<-par[3]
 		a<-par[1:(tree$Nnode-1)+3]
-		logLik<-sum(dmnorm(c(x,a),mean=rep(a0,length(tree$tip.label)+tree$Nnode-1),
-			varcov=sig2*vcvPhylo(tree,model="OU",alpha=alpha),log=TRUE))
+		## logLik<-sum(dmnorm(c(x,a),mean=rep(a0,length(tree$tip.label)+tree$Nnode-1),
+		##	varcov=sig2*vcvPhylo(tree,model="OU",alpha=alpha),log=TRUE))
+		logLik<-logMNORM(c(x,a),rep(a0,Ntip(tree)+tree$Nnode-1),sig2*vcvPhylo(tree,model="OU",alpha=alpha))
 		if(trace) print(c(sig2,alpha,logLik))
 		-logLik
 	}
 	x<-x[tree$tip.label]
 	pp<-rep(NA,tree$Nnode+2)
 	pp[1:tree$Nnode+2]<-fastAnc(tree,x)
 	pp[1]<-phyl.vcv(as.matrix(c(x,pp[2:tree$Nnode+2])),vcvPhylo(tree),lambda=1)$R[1,1]
-	pp[2]<-tol ## arbitrarily
+	pp[2]<-a.init ## arbitrarily
 	fit<-optim(pp,likOU,tree=tree,x=x,trace=trace,method="L-BFGS-B",
 		lower=c(tol,tol,rep(-Inf,tree$Nnode)),upper=rep(Inf,length(pp)),
 		control=list(maxit=maxit))
@@ -89,6 +92,9 @@ anc.OU<-function(tree,x,maxit=2000,...){
 	obj
 }
 
+logMNORM<-function(x,x0,vcv)
+	-t(x-x0)%*%solve(vcv)%*%(x-x0)/2-length(x)*log(2*pi)/2-determinant(vcv,logarithm=TRUE)$modulus[1]/2
+
 ## print method for "anc.ML"
 ## written by Liam J. Revell 2015
 print.anc.ML<-function(x,digits=6,printlen=NULL,...){

R/contMap.R

@@ -18,16 +18,41 @@ contMap<-function(tree,x,res=100,fsize=NULL,ftype=NULL,lwd=4,legend=NULL,
 	steps<-0:res/res*max(h)
 	H<-nodeHeights(tree)
 	if(method=="fastAnc") a<-fastAnc(tree,x) 
-	else { 
+	else if(method=="anc.ML") { 
 		fit<-anc.ML(tree,x)
 		a<-fit$ace
 		if(!is.null(fit$missing.x)) x<-c(x,fit$missing.x)
+	} else if(method=="user"){
+		if(hasArg(anc.states)) anc<-list(...)$anc.states
+		else {
+			cat("No ancestral states have been provided. Using states estimated with fastAnc.\n\n")
+			a<-fastAnc(tree,x)
+		}
+		if(length(anc.states)<tree$Nnode){
+			nodes<-as.numeric(names(anc.states))
+			tt<-tree
+			for(i in 1:length(nodes)){
+				M<-matchNodes(tt,tree,method="distances")
+				ii<-M[which(M[,2]==nodes[i]),1]
+				tt<-bind.tip(tt,nodes[i],edge.length=0,where=ii)
+			}
+			a<-fastAnc(tt,c(x,anc.states))
+			M<-matchNodes(tree,tt,method="distances")
+			a<-a[as.character(M[,2])]
+			names(a)<-M[,1]
+		} else { 
+			if(is.null(names(anc.states))) names(anc.states)<-1:tree$Nnode+Ntip(tree)
+			a<-a[as.character(1:tree$Nnode+Ntip(tree))]
+		}
 	}
-	y<-c(a,x[tree$tip.label]); names(y)[1:length(tree$tip)+tree$Nnode]<-1:length(tree$tip)
+	y<-c(a,x[tree$tip.label])
+	names(y)[1:Ntip(tree)+tree$Nnode]<-1:Ntip(tree)
 	A<-matrix(y[as.character(tree$edge)],nrow(tree$edge),ncol(tree$edge))
-	cols<-rainbow(1001,start=0,end=0.7); names(cols)<-0:1000
-	if(is.null(lims)) lims<-c(min(x),max(x))
-	trans<-0:1000/1000*(lims[2]-lims[1])+lims[1]; names(trans)<-0:1000
+	cols<-rainbow(1001,start=0,end=0.7)
+	names(cols)<-0:1000
+	if(is.null(lims)) lims<-c(min(y),max(y))
+	trans<-0:1000/1000*(lims[2]-lims[1])+lims[1]
+	names(trans)<-0:1000
 	for(i in 1:nrow(tree$edge)){
 		XX<-cbind(c(H[i,1],steps[intersect(which(steps>H[i,1]),which(steps<H[i,2]))]),
 			c(steps[intersect(which(steps>H[i,1]),which(steps<H[i,2]))],H[i,2]))-H[i,1]