## function to compute the marginal posterior probabilities for nodes
## using the rerooting method
## written by Liam J. Revell 2013, 2015, 2017, 2018, 2020, 2023
rerootingMethod<-function(tree,x,model=c("ER","SYM"),...){
if(hasArg(quiet)) quiet<-list(...)$quiet
else quiet<-FALSE
if(!quiet){
cat("\nNote:\n")
cat(" This function is redundant with 'phytools::ancr' in situations in\n")
cat(" which it should be used (symmetric Q matrices) & invalid for non-\n")
cat(" symmetric Q matrices (e.g., model='ARD').\n\n")
}
if(!inherits(tree,"phylo"))
stop("tree should be an object of class \"phylo\".")
if(hasArg(tips)) tips<-list(...)$tips
else tips<-NULL
if(!is.matrix(model)) model<-model[1]
n<-Ntip(tree)
# if vector convert to binary matrix
if(!is.matrix(x)){
yy<-to.matrix(x,sort(unique(x)))
if(is.null(tips)) tips<-FALSE
} else {
if(is.null(tips)) tips<-TRUE
yy<-x
}
yy<-yy[tree$tip.label,]
yy<-yy/rowSums(yy)
YY<-fitMk(tree,yy,model=model,output.liks=TRUE,...)
Q<-matrix(c(0,YY$rates)[YY$index.matrix+1],length(YY$states),
length(YY$states),dimnames=list(YY$states,YY$states))
diag(Q)<--colSums(Q,na.rm=TRUE)
nn<-if(tips) c(1:n,if(tree$Nnode>1) 2:tree$Nnode+n) else {
if(tree$Nnode>1) 2:tree$Nnode+n else vector() }
ff<-function(nn){
tt<-if(nn>Ntip(tree)) ape::root.phylo(tree,node=nn) else
reroot(tree,nn,tree$edge.length[which(tree$edge[,2]==nn)])
fitMk(tt,yy,model=model,fixedQ=Q,output.liks=TRUE)$lik.anc[1,]
}
XX<-t(sapply(nn,ff))
if(tips) XX<-rbind(XX[1:n,],YY$lik.anc[1,],if(tree$Nnode>1)
XX[(n+1):nrow(XX),])
else XX<-rbind(yy,YY$lik.anc[1,],if(tree$Nnode>1) XX)
rownames(XX)<-1:(tree$Nnode+n)
if(tips) rownames(XX)[1:n]<-tree$tip.label
XX<-if(tips) XX else XX[1:tree$Nnode+n,]
obj<-list(loglik=YY$logLik,Q=Q,marginal.anc=XX,tree=tree,x=yy)
class(obj)<-"rerootingMethod"
obj
}
print.rerootingMethod<-function(x,digits=6,printlen=NULL,...){
cat("Ancestral character estimates using re-rooting method")
cat("\nof Yang et al. (1995):\n")
if(is.null(printlen)) print(round(x$marginal.anc,digits)) else {
print(round(x$marginal.anc[1:printlen,],digits))
cat("...\n")
}
cat("\nEstimated transition matrix,\nQ =\n")
print(round(x$Q,digits))
cat("\n**Note that if Q is not symmetric the marginal")
cat("\nreconstructions may be invalid.\n")
cat(paste("\nLog-likelihood =",round(x$loglik,digits),"\n\n"))
}
plot.rerootingMethod<-function(x, ...){
args<-list(...)
if(is.null(args$lwd)) args$lwd<-1
if(is.null(args$ylim)) args$ylim<-c(-0.1*Ntip(x$tree),
Ntip(x$tree))
if(is.null(args$offset)) args$offset<-0.5
if(is.null(args$ftype)) args$ftype="i"
args$tree<-x$tree
do.call(plotTree,args)
if(hasArg(piecol)) piecol<-list(...)$piecol
else piecol<-setNames(colorRampPalette(c("blue",
"yellow"))(ncol(x$marginal.anc)),
colnames(x$marginal.anc))
if(hasArg(node.cex)) node.cex<-list(...)$node.cex
else node.cex<-0.6
nodelabels(pie=x$marginal.anc[
as.character(1:x$tree$Nnode+Ntip(x$tree)),],
piecol=piecol,cex=node.cex)
if(hasArg(tip.cex)) tip.cex<-list(...)$tip.cex
else tip.cex<-0.4
tiplabels(pie=x$x[x$tree$tip.label,],piecol=piecol,
cex=tip.cex)
legend(x=par()$usr[1],y=par()$usr[1],
legend=colnames(x$marginal.anc),
pch=21,pt.bg=piecol,pt.cex=2.2,bty="n")
}
logLik.rerootingMethod<-function(object,...) object$loglik
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.