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/
