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
