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/