Hi everyone, I was wondering if anyone knows a good R function or package for doing a maximum likelihood search for a phylogenetic tree using continuous characters. I'm aware of Rcontml, but I'm looking for something that I can change and modify the individual steps in the search process rather than calling out to another program. I am currently using the script below, which I modified from Paradis' ML search using molecular data. The problem that I'm running into, is that lots of topologies are getting higher likelihood scores than the tree that I'm using to simulate the characters. I'm using rtree to generate a random tree and using fastBM to simulate the characters with a Brownian Motion model. Using ace, the topology that I used to generate the data doesn't even come close to being the most likely topology. If anyone has suggestions, I am keen to hear them. Best regards, Will Gelnaw
Here's the script I've been using: tree<-rtree(n=40) dat<-fastBM(tree, nsim=20) ## calculates the sum of the log likelihoods of each of the characters sumloglik<-function(tree,dat,dattype="continuous"){ loglik<-vector() nchar<-ncol(dat) for(i in 1:nchar) loglik[i]<-ace(dat[,i],tree,type=dattype)$resloglik sumlik<-sum(loglik) return(sumlik)} ## uses the above sumloglik function to do the tree search bigtreesearch<-function(dat,dattype="continuous",logreps=200){ #search breadth is the fraction of the tree that will be changed during NNI transformations #0.05 means that 5% of the branches are moved during the NNI shift test_tree<-rtree(n=nrow(dat),tip.label=row.names(dat)) ntips<-nrow(dat) print("iteration 1") write.tree(test_tree,"tree_log.tre", append=TRUE) #appends the current tree to the tree log and creates a new log if none exists loglik<-vector(); loglik[1]<-sumloglik(test_tree, dat,dattype="continuous"); print(loglik[1]) i<-2 ## sets up the number of rounds of improvements m<-1 ## sets up the number of iterations that it tries to do a search while (i<=logreps){ iteration<-paste("iteration", i) print(iteration) TR2<-rNNI(test_tree,moves=1,n=1) d=0.1 ## sets up the amount of random noise added to the branch lengths if (m>14) d=0.2 ## as the analysis tests out alternatives, it makes progressively larger and larger changes if (m>29) d=0.3 if (m>34) d=0.4 TR2$edge.length<-TR2$edge.length+rnorm(2*ntips-2,sd=d) TR2$edge.length[TR2$edge.length<=0]<-0.001 #above 5 lines add noise to the branch lengths #negative branch lengths are eliminated by changing 0-length or negative branch lengths equal to a small size. loglik[i]<-sumloglik(TR2,dat,"continuous") print(loglik[i]) accept<- {if((loglik[i])>(loglik[i-1])) TRUE else FALSE} ##tests for improvement from previous test if(accept){ test_tree<-TR2 write.tree(tree,"tree_log.tre", append=TRUE) m<-1 ## resets m back to 1 i<-i+1 } else { ## advances to the next round of testing m<-m+1 if (m>49) break} ##if the process goes 19 iterations without improvement, it quits } best_tree<-list() best_tree$tree<-test_tree char<-backcalc(test_tree,dat) best_tree$loglik<-loglik[i-1] return(best_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/