Permalink
Browse files

new function to compare rates among trees

  • Loading branch information...
1 parent 9d03636 commit 229fdb808467fc837ae6195157ca8aba62a3f90f @liamrevell committed Feb 2, 2017
Showing with 143 additions and 6 deletions.
  1. +4 −4 DESCRIPTION
  2. +3 −2 NAMESPACE
  3. +100 −0 R/ratebytree.R
  4. +36 −0 man/ratebytree.Rd
View
@@ -1,6 +1,6 @@
Package: phytools
-Version: 0.5-71
-Date: 2017-01-31
+Version: 0.5-72
+Date: 2017-02-02
Title: Phylogenetic Tools for Comparative Biology (and Other Things)
Author: Liam J. Revell
Maintainer: Liam J. Revell <liam.revell@umb.edu>
@@ -56,6 +56,6 @@ Description: Package contains various functions for phylogenetic analysis.
research.
License: GPL (>= 2)
URL: http://github.com/liamrevell/phytools
-Packaged: 2017-01-31 12:00:00 EST
+Packaged: 2017-02-02 12:00:00 EST
Repository:
-Date/Publication: 2017-01-31 12:00:00 EST
+Date/Publication: 2017-02-02 12:00:00 EST
View
@@ -20,8 +20,8 @@ export(phyl.vcv, phylo.heatmap, phylo.toBackbone, phylo.to.map, phylANOVA, phylo
export(plotBranchbyTrait, plot.contMap, plot.cophylo, plot.densityMap, plot.fitMk, plot.fitPagel, plot.phylo.to.map, plot.gfit)
export(plotSimmap, plotThresh, plotTree, plotTree.barplot, plotTree.boxplot, plotTree.singletons, plotTree.splits, plotTree.wBars)
export(posterior.evolrate)
-export(ratebystate, rateshift, read.newick, read.simmap, reorder.backbonePhylo, reorderSimmap, rep.multiPhylo, rep.phylo, repPhylo)
-export(reroot, rerootingMethod, rescaleSimmap, resolveAllNodes, resolveNode, rotate.multi, rotateNodes, rootedge.to.singleton)
+export(ratebystate, ratebytree, rateshift, read.newick, read.simmap, reorder.backbonePhylo, reorderSimmap, rep.multiPhylo, rep.phylo)
+export(repPhylo, reroot, rerootingMethod, rescaleSimmap, resolveAllNodes, resolveNode, rotate.multi, rotateNodes, rootedge.to.singleton)
export(roundBranches, roundPhylogram, rstate)
export(sampleFrom, setMap, sim.corrs, sim.history, sim.ratebystate, sim.rates, skewers, splitEdgeColor, splitplotTree)
export(splitTree, starTree, strahlerNumber)
@@ -99,6 +99,7 @@ S3method(print, phyl.RMA)
S3method(coef, phyl.RMA)
S3method(residuals, phyl.RMA)
S3method(plot, phyl.RMA)
+S3method(print, ratebytree)
importFrom(animation, ani.options, ani.record, ani.replay, saveVideo)
importFrom(ape, .PlotPhyloEnv, .uncompressTipLabel, ace, all.equal.phylo, as.DNAbin, as.phylo, bind.tree, branching.times, collapse.singles)
View
@@ -0,0 +1,100 @@
+## method to compare the rate of evolution for a character between trees
+## closely related to 'censored' approach of O'Meara et al. (2006; Evolution)
+## written by Liam J. Revell 2017
+
+ratebytree<-function(trees,x,...){
+ if(hasArg(tol)) tol<-list(...)$tol
+ else tol<-1e-8
+ ## check trees & x
+ if(!inherits(trees,"multiPhylo"))
+ stop("trees should be object of class \"multiPhylo\".")
+ if(!is.list(x)) stop("x should be a list of vectors.")
+ N<-length(trees)
+ ## reorder the trait vectors in x
+ x<-mapply(function(x,t) x<-x[t$tip.label],x=x,t=trees)
+ ## first, fit multi-rate model
+ lik.multi<-function(theta,trees,x){
+ n<-sapply(trees,Ntip)
+ N<-length(trees)
+ sig<-theta[1:N]
+ a<-theta[(N+1):(2*N)]
+ C<-lapply(trees,vcv)
+ V<-mapply("*",C,sig)
+ logL<-0
+ for(i in 1:N)
+ logL<-logL-t(x[[i]]-a[i])%*%solve(V[[i]])%*%(x[[i]]-
+ a[i])/2-n[i]*log(2*pi)/2-determinant(V[[i]])$modulus[1]/2
+ -logL
+ }
+ foo<-function(tree,x){
+ pvcv<-phyl.vcv(as.matrix(x),vcv(tree),1)
+ c(pvcv$R[1,1],pvcv$a[1,1])
+ }
+ PP<-mapply(foo,trees,x)
+ p<-as.vector(t(PP))
+ fit.multi<-optim(p,lik.multi,trees=trees,x=x,method="L-BFGS-B",
+ lower=c(rep(tol,N),rep(-Inf,N)),upper=c(rep(Inf,N),rep(Inf,N)))
+ ## now fit single-rate model
+ lik.onerate<-function(theta,trees,x){
+ n<-sapply(trees,Ntip)
+ N<-length(trees)
+ sig<-theta[1]
+ a<-theta[1:N+1]
+ C<-lapply(trees,vcv)
+ V<-lapply(C,"*",sig)
+ logL<-0
+ for(i in 1:N)
+ logL<-logL-t(x[[i]]-a[i])%*%solve(V[[i]])%*%(x[[i]]-
+ a[i])/2-n[i]*log(2*pi)/2-determinant(V[[i]])$modulus[1]/2
+ -logL
+ }
+ p<-c(mean(fit.multi$par[1:N]),fit.multi$par[(N+1):(2*N)])
+ fit.onerate<-optim(p,lik.onerate,trees=trees,x=x,method="L-BFGS-B",
+ lower=c(tol,rep(-Inf,N)),upper=c(Inf,rep(Inf,N)))
+ ## compare models:
+ LR<-2*(-fit.multi$value+fit.onerate$value)
+ P.chisq<-pchisq(LR,df=N-1,lower.tail=FALSE)
+ obj<-list(
+ multi.rate.model=list(sig2=fit.multi$par[1:N],
+ a=fit.multi$par[(N+1):(2*N)],
+ k=2*N,
+ logL=-fit.multi$value,
+ counts=fit.multi$counts,convergence=fit.multi$convergence,
+ message=fit.multi$message),
+ common.rate.model=list(sig2=fit.onerate$par[1],
+ a=fit.onerate$par[1:N+1],
+ k=N+1,
+ logL=-fit.onerate$value,
+ counts=fit.onerate$counts,convergence=fit.onerate$convergence,
+ message=fit.onerate$message),
+ N=N,likelihood.ratio=LR,P.chisq=P.chisq)
+ class(obj)<-"ratebytree"
+ obj
+}
+
+## S3 print method for ratebytree
+print.ratebytree<-function(x,...){
+ if(hasArg(digits)) digits<-list(...)$digits
+ else digits<-4
+ N<-x$N
+ cat("ML common-rate model:\n")
+ cat(paste("\ts^2\t",paste(paste("a[",1:N,"]",sep=""),collapse="\t")),
+ "\tk\tlogL\n")
+ cat(paste("value",round(x$common.rate.model$sig2,digits),
+ paste(round(x$common.rate.model$a,digits),collapse="\t"),
+ x$common.rate.model$k,round(x$common.rate.model$logL,digits),
+ "\n\n",sep="\t"))
+ cat("ML multi-rate model:\n")
+ cat(paste("\t",paste(paste("s^2[",1:N,"]",sep=""),collapse="\t"),"\t",
+ paste(paste("a[",1:N,"]",sep=""),collapse="\t")),
+ "\tk\tlogL\n")
+ cat(paste("value",paste(round(x$multi.rate.model$sig2,digits),collapse="\t"),
+ paste(round(x$multi.rate.model$a,digits),collapse="\t"),
+ x$multi.rate.model$k,round(x$multi.rate.model$logL,digits),
+ "\n\n",sep="\t"))
+ cat(paste("Likelihood ratio:",round(x$likelihood.ratio,digits),"\n"))
+ cat(paste("P-value (based on X^2):",round(x$P.chisq,digits),"\n\n"))
+ if(x$multi.rate.model$convergence==0&&x$common.rate.model$convergence==0)
+ cat("R thinks it has found the ML solution.\n\n")
+ else cat("One or the other optimization may not have converged.\n\n")
+}
View
@@ -0,0 +1,36 @@
+\name{ratebytree}
+\alias{ratebytree}
+\title{Likelihood test for rate variation among trees}
+\usage{
+ratebytree(trees, x, ...)
+}
+\arguments{
+ \item{trees}{an object of class \code{"multiPhylo"}.}
+ \item{x}{a list of trait vectors for a continuous trait in which the names of each vectors correspond to the tip labels of \code{trees}.}
+ \item{...}{optional arguments.}
+}
+\description{
+ This function takes a object of class \code{"multiPhylo"} containing two or more phylogenies (\code{trees}), and a list of trait vectors (\code{x}), and fits two models: one in which the rate of trait evolution is equal among all trees; and a second in which the rates can differ between trees.
+
+ This model corresponds to the \emph{censored} approach of O'Meara et al. (2006; \emph{Evolution}) and should also be related to the method of Adams (2012; \emph{Systematic Biology}) for comparing rates among traits. See \code{\link{brownie.lite}} for a different implementation of the \emph{noncensored} approach of O'Meara et al. (2006).
+}
+\details{
+ The function also conducts a likelihood-ratio test to compare the two models.
+}
+\value{
+ An object of class \code{"ratebytree"}.
+}
+\references{
+ Adams, D. C. (2012) Comparing evolutionary rates for different phenotypic traits on a phylogeny using likelihood. \emph{Syst. Biol.}, \bold{62}, 181-192.
+
+ O'Meara, B. C., C. Ane, M. J. Sanderson, and P. C. Wainwright. (2006) Testing for different rates of continuous trait evolution using likelihood. \emph{Evolution}, \bold{60}, 922-933.
+
+ Revell, L. J. (2012) phytools: An R package for phylogenetic comparative biology (and other things). \emph{Methods Ecol. Evol.}, \bold{3}, 217-223.
+}
+\author{Liam Revell \email{liam.revell@umb.edu}}
+\seealso{
+ \code{\link{brownie.lite}}
+}
+\keyword{phylogenetics}
+\keyword{comparative method}
+\keyword{maximum likelihood}

0 comments on commit 229fdb8

Please sign in to comment.