Hi all,
I tried the function written by Liam, it works pretty well! As an indication for users, do not use an ultrametric tree, you'd get an error like this: Error in chol.default(x, pivot = FALSE) : the leading minor of order 140 is not positive definite Error in pd.solve(varcov) : x appears to be not positive definite I applied anc.tree to a multichotomous tree, it still works very well, judging from fitted values at nodes. Pas >----Messaggio originale---- >Da: liam.rev...@umb.edu >Data: 22/02/2011 4.58 >A: "pasquale.r...@libero.it"<pasquale.r...@libero.it> >Cc: "R Sig Phylo Listserv"<r-sig-phylo@r-project.org> >Ogg: Re: [R-sig-phylo] R: Re: R: Re: Simulation of Continuous Trait with Active Trend > >Hi Pasquale, > >To follow up on my earlier message, I have written a preliminary version >of a simple function to simultaneously estimate the evolutionary >parameters and ancestral states for Brownian evolution with a trend >using likelihood. I will paste it at the end of this email. I have >confirmed that the function returns the same trend parameter as >fitContinuous(model="trend") and the same ancestral character estimates >as ace() - when mu (the trend parameter) is forced to zero [note that >this can only be done by modifying the code]. The parameter "sig2," the >Brownian rate, is biased by a factor n/(2n-2) for n species relative to >fitContinuous(model="trend)$Trait1$beta. > >Note that this model can generally only be fit for a tree with tips that >are not contemporaneous - for instance, for a phylogeny containing >fossil species. > >Best, Liam > >Function code (also >http://anolis.oeb.harvard.edu/~liam/R-phylogenetics/anc.trend/v0.1/anc.trend. R): > ># written by Liam J. Revell 2011 >anc.trend<-function(phy,x,maxit=2000){ > > # preliminaries > # require dependencies > require(ape) > # set global > tol<-1e-8 > # compute C > D<-dist.nodes(phy) > ntips<-length(phy$tip.label) > Cii<-D[ntips+1,] > C<-D; C[,]<-0 > counts<-vector() > for(i in 1:nrow(D)) for(j in 1:ncol(D)) > C[i,j]<-(Cii[i]+Cii[j]-D[i,j])/2 > dimnames(C)[[1]][1:length(phy$tip)]<-phy$tip.label > dimnames(C)[[2]][1:length(phy$tip)]<-phy$tip.label > C<-C[c(1:ntips,(ntips+2):nrow(C)),c(1:ntips,(ntips+2):ncol(C))] > # sort x by phy$tip.label > x<-x[phy$tip.label] > > # function returns the negative log-likelihood > likelihood<-function(theta,x,C){ > a<-theta[1] > u<-theta[2] > sig2<-theta[3] > y<-theta[4:length(theta)] > > logLik<-dmnorm(x=c(x,y),mean=(a+diag(C)*u),varcov=sig2*C,log=TRUE) > return(-logLik) > } > > # get reasonable starting values for the optimizer > a<-mean(x) > sig2<-var(x)/max(C) > > # perform ML optimization > > result<-optim(par=c(a,0,sig2,rep(a,phy$Nnode-1)),likelihood,x=x,C=C,method=" L-BFGS-B",lower=c(-Inf,-Inf,tol,rep(-Inf,phy$Nnode-1)),control=list (maxit=maxit)) > > # return the result > ace<-c(result$par[c(1,4:length(result$par))]); >names(ace)<-c(as.character(phy$edge[1,1]),rownames(C)[(length(phy$tip.label) +1):nrow(C)]) > > return(list(ace=ace,mu=result$par[2],sig2=result$par[3],logL=-result$value, convergence=result$convergence,message=result$message)) > >} > >-- >Liam J. Revell >University of Massachusetts Boston >web: http://www.bio.umb.edu/facstaff/faculty_Revell.html >(new) email: liam.rev...@umb.edu >(new) blog: http://phytools.blogspot.com > >On 2/17/2011 12:45 PM, Liam J. Revell wrote: >> Hi Pasquale, >> >> I think this may be possible if all tips are not contemporaneous. >> >> As Gene points out, the values at the tips of the tree given a trend >> should be distributed as ~MVN(mu=diag(vcv(tree))*tr+ancestor, >> V=vcv(tree)*sig2). Thus, we should be able to compute the likelihood of >> any data pattern at the tips and internal nodes of the tree, given tr, >> ancestor, and sig2. If we start with only the values at the tips and the >> phylogeny, we can estimate the states at internal nodes and our model by >> maximizing the likelihood. This would be our solution. >> >> I have written down some code for this using dmnorm() (from the mnormt >> package) to compute the multivariate normal density; a home-made >> function to compute vcv(tree) but for all the nodes in the tree; and >> optim() to maximize the likelihood - but it doesn't seem to work very >> well yet. I will pass it along if I can figure it out. >> >> - Liam >> > _______________________________________________ R-sig-phylo mailing list R-sig-phylo@r-project.org https://stat.ethz.ch/mailman/listinfo/r-sig-phylo
