Hi Everyone,
I'm working on testing the effects of a likelihood reweighting algorithm
that I think will improve the use of morphometrics in phylogenetic
analyses. To do that, I need to do a bunch of tree searches using
simulated data. I'm using a modified form of the tree search script in the
Paradis Book. I'm using ace {ape} to estimate the likelihood of each
character instead of using pmi {phangorn} because I'm interested using
morphometrics. For the sake of this discussion, I haven't included the
reweighting script.
Because all of the taxa that I'm looking at are extant, I want to
constrain my search to ultrametric trees. I've been using the chronos
function {ape} to ultrametricize the trees, but that seems really
inefficient. I was hoping someone knew a better way to only search
ultrametric trees.
The other thing I've noticed is that even a short tree search often
produces a tree over which the characters have a higher likelihood than the
true tree that I simulated using rtree {ape} and sim.char {geiger}. I was
wondering if anyone else had encountered this and had an explanation. I've
included my code below.
Best regards
- Will Gelnaw
require(ape)
require
## calculates the sum of the log likelihoods of each of the characters
sumloglik<-function(tree,dat,dattype){
loglik<-vector()
nchar<-ncol(dat)
for(i in 1:nchar)
loglik[i]<-ace(dat[,i],tree,type=dattype,method="ML")$loglik
sumlik<-sum(loglik)
return(sumlik)}
## Conducts the tree search
## maxreps is the maximum number of threes for the program to search through
mltreesearch<-function(dat,maxreps,ultrametric=TRUE){
tree<-rtree(n=nrow(dat),tip.label=row.names(dat))
if (ultrametric==TRUE) tree<-chronos(tree,lambda=1)
#if an ultrametric tree is desired, using chronos smooths out the
current tree into an ultrametric tree
#using chronos here allows for future functionality of adding
calibration points to the tree.
ntips<-length(tree$tip.label)
loglik<-vector()
loglik[1]<-sumloglik(tree, dat,dattype="continuous")
i<-2
j<-1
while (i<=maxreps){
if(j>35) break
iteration<-paste("iteration", i)
print(iteration)
TR2<-rNNI(tree,moves=1,n=1)
TR2$edge.length<-TR2$edge.length+rnorm(2*ntips-2,sd=0.1)
TR2$edge.length[TR2$edge.length<=0]<-0.001
if (ultrametric==TRUE) TR2<-chronos(TR2,lambda=1)
loglik[i]<-sumloglik(TR2,dat,"continuous")
accept<- if (loglik[i]>loglik[i-1]) TRUE
else FALSE
if (accept){
tree<-TR2
i<-i+1
print(loglik[i])}
else j<-j+1}
return(tree)}
[[alternative HTML version deleted]]
_______________________________________________
R-sig-phylo mailing list - R-sig-phylo@r-project.org
https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
Searchable archive at http://www.mail-archive.com/r-sig-phylo@r-project.org/
