Permalink
Browse files

model EB for anc.ML

  • Loading branch information...
liamrevell committed Aug 26, 2017
1 parent 762158a commit 69426297144ba4f4c9238f423c7d1e606939b32d
Showing with 93 additions and 8 deletions.
  1. +70 −2 R/anc.ML.R
  2. +19 −2 R/utilities.R
  3. +4 −4 man/anc.ML.Rd
View
@@ -2,10 +2,11 @@
## also allows missing data in x, in which case missing data are also estimated
## written by Liam J. Revell 2011, 2013, 2014, 2015
anc.ML<-function(tree,x,maxit=2000,model=c("BM","OU"),...){
anc.ML<-function(tree,x,maxit=2000,model=c("BM","OU","EB"),...){
if(!inherits(tree,"phylo")) stop("tree should be an object of class \"phylo\".")
if(model[1]=="BM") obj<-anc.BM(tree,x,maxit,...)
else if(model[1]=="OU") obj<-anc.OU(tree,x,maxit,...)
else if(model[1]=="EB") obj<-anc.EB(tree,x,maxit,...)
else stop(paste("Do not recognize method",model))
obj
}
@@ -119,13 +120,74 @@ anc.OU<-function(tree,x,maxit=2000,...){
obj
}
## EB is the Early-burst model (Harmon et al. 2010) and also called the ACDC model
## (accelerating-decelerating; Blomberg et al. 2003). Set by the a rate parameter, EB fits a model where
## the rate of evolution increases or decreases exponentially through time, under the model
## r[t] = r[0] * exp(a * t), where r[0] is the initial rate, a is the rate change parameter, and t is
## time. The maximum bound is set to -0.000001, representing a decelerating rate of evolution. The minimum
## bound is set to log(10^-5)/depth of the tree.
## internal to estimate ancestral states under an EB model
## written by Liam J. Revell 2017
anc.EB<-function(tree,x,maxit=2000,...){
## check to see if any tips are missing data
xx<-setdiff(tree$tip.label,names(x))
if(length(xx)>0) stop("Some tips of the tree do not have data. Try model=\"BM\".")
if(hasArg(tol)) tol<-list(...)$tol
else tol<-1e-8
if(hasArg(trace)) trace<-list(...)$trace
else trace<-FALSE
if(hasArg(r.init)){
r.init<-list(...)$r.init
obj<-phyl.vcv(as.matrix(x[tree$tip.label]),
vcv(ebTree(tree,r.init)),1)
s2.init<-obj$R[1,1]
a0.init<-obj$alpha[1,1]
} else {
## optimize r.init
lik<-function(p,tree,x)
logLik<--logMNORM(x,rep(p[3],Ntip(tree)),
p[1]*vcvPhylo(tree,model="EB",r=p[2],anc.nodes=F))
obj<-phyl.vcv(as.matrix(x[tree$tip.label]),vcv(tree),1)
fit.init<-optim(c(obj$R[1,1],0,obj$alpha[1,1]),
lik,tree=tree,x=x,method="L-BFGS-B",lower=c(tol,-Inf,-Inf),
upper=rep(Inf,3))
r.init<-fit.init$par[2]
s2.init<-fit.init$par[1]
a0.init<-fit.init$par[3]
}
likEB<-function(par,tree,x,trace){
sig2<-par[1]
r<-par[2]
obj<-fastAnc(ebTree(tree,r),x)
a0<-obj[1]
a<-obj[2:length(obj)]
logLik<-logMNORM(c(x,a),rep(a0,Ntip(tree)+tree$Nnode-1),
sig2*vcvPhylo(tree,model="EB",r=r))
if(trace) print(c(sig2,r,logLik))
-logLik
}
x<-x[tree$tip.label]
pp<-rep(NA,2)
pp[1]<-s2.init
pp[2]<-r.init
fit<-optim(pp,likEB,tree=tree,x=x,trace=trace,method="L-BFGS-B",
lower=c(tol,-Inf),upper=rep(Inf,2),control=list(maxit=maxit))
obj<-list(sig2=fit$par[1],r=fit$par[2],
ace=unclass(fastAnc(ebTree(tree,fit$par[2]),x)),
logLik=-fit$value,counts=fit$counts,convergence=fit$convergence,
message=fit$message,model="EB")
class(obj)<-"anc.ML"
obj
}
logMNORM<-function(x,x0,vcv)
-t(x-x0)%*%solve(vcv)%*%(x-x0)/2-length(x)*log(2*pi)/2-determinant(vcv,logarithm=TRUE)$modulus[1]/2
## print method for "anc.ML"
## written by Liam J. Revell 2015, 2016
print.anc.ML<-function(x,digits=6,printlen=NULL,...){
cat(paste("Ancestral character estimates using anc.ML under a",
cat(paste("Ancestral character estimates using anc.ML under a(n)",
x$model,"model:\n"))
Nnode<-length(x$ace)
if(is.null(printlen)||printlen>=Nnode) print(round(x$ace,digits))
@@ -153,6 +215,12 @@ print.anc.ML<-function(x,digits=6,printlen=NULL,...){
colnames(obj)<-c("sigma^2","alpha","logLik")
rownames(obj)<-""
print(obj)
} else if(x$model=="EB"){
obj<-data.frame(round(x$sig2,digits),round(x$r,digits),
round(x$logLik,digits))
colnames(obj)<-c("sigma^2","r","logLik")
rownames(obj)<-""
print(obj)
}
if(x$convergence==0) cat("\nR thinks it has found the ML solution.\n\n")
else cat("\nOptimization may not have converged.\n\n")
View
@@ -1,6 +1,18 @@
## some utility functions
## written by Liam J. Revell 2011, 2012, 2013, 2014, 2015, 2016, 2017
## function to rescale a tree according to an EB model
## written by Liam J. Revell 2017
ebTree<-function(tree,r){
if(r!=0){
H<-nodeHeights(tree)
e<-(exp(r*H[,2])-exp(r*H[,1]))/r
tree$edge.length<-e
}
tree
}
## function to expand clades in a plot by a given factor
## written by Liam J. Revell 2017
expand.clade<-function(tree,node,factor=5){
@@ -1222,7 +1234,6 @@ vcvPhylo<-function(tree,anc.nodes=TRUE,...){
}
if(hasArg(model)) model<-list(...)$model
else model<-"BM"
if(hasArg(tol)) tol<-list(...)$tol
else tol<-1e-12
if(model=="OU"){
@@ -1236,7 +1247,13 @@ vcvPhylo<-function(tree,anc.nodes=TRUE,...){
tree<-lambdaTree(tree,lambda)
} else model<-"BM"
model<-"BM"
}
}
if(model=="EB"){
if(hasArg(r)){
r<-list(...)$r
tree<-ebTree(tree,r)
} else model<-"BM"
}
# done settings
n<-length(tree$tip.label)
h<-nodeHeights(tree)[order(tree$edge[,2]),2]
View
@@ -2,23 +2,23 @@
\alias{anc.ML}
\title{Ancestral character estimation using likelihood}
\usage{
anc.ML(tree, x, maxit=2000, model=c("BM","OU"), ...)
anc.ML(tree, x, maxit=2000, model=c("BM","OU","EB"), ...)
}
\arguments{
\item{tree}{an object of class \code{"phylo"}.}
\item{x}{a vector of tip values for species; \code{names(x)} should be the species names.}
\item{maxit}{an optional integer value indicating the maximum number of iterations for optimization.}
\item{model}{model of continuous character evolution ont he tree. It's possible that only \code{model="BM"} works in the present version as \code{model="OU"} has not be thoroughly tested & some bugs were reported for an earlier version.}
\item{model}{model of continuous character evolution ont he tree. It's possible that only \code{model="BM"} & \code{model="EB"} work in the present version as \code{model="OU"} has not be thoroughly tested & some bugs were reported for an earlier version.}
\item{...}{other arguments.}
}
\description{
This function estimates the evolutionary parameters and ancestral states for Brownian evolution using likelihood. It is also possible (for \code{model="BM"}) to allow for missing data for some tip taxa.
}
\details{
Because this function relies on a high dimensional numerical optimization of the likelihood function, \code{\link{fastAnc}} should probably be preferred for most purposes. If using \code{\link{anc.ML}}, users should be cautious to ensure convergence. This has been ameliorated in phytools>=0.2-48 by seeding the ML optimization with the result from \code{\link{fastAnc}}.
Because this function relies on a high dimensional numerical optimization of the likelihood function, \code{\link{fastAnc}} should probably be preferred for most purposes. If using \code{\link{anc.ML}}, users should be cautious to ensure convergence. This has been ameliorated in phytools>=0.2-48 by seeding the ML optimization with the result from \code{\link{fastAnc}}. For \code{model="EB"} this should also not be a problem as the numerical optimization is performed for only \code{sig2} and \code{r}, while the ML values of the ancestral states are obtained during every iteration of the optimization algorithmically using the re-rooting method.
}
\value{
A list with the following components:
An object of class \code{"anc.ML"} with at least the following four elements (if not more, depending on \code{model}):
\item{sig2}{the variance of the BM process.}
\item{ace}{a vector with the ancestral states.}
\item{logLik}{the log-likelihood.}

0 comments on commit 6942629

Please sign in to comment.