Permalink
Browse files

more models for type=diversification in ratebytree

  • Loading branch information...
liamrevell committed Oct 29, 2017
1 parent 242b1cd commit f8f6b898b7c8d740b1f67d033750e99282f73d3a
Showing with 63 additions and 16 deletions.
  1. +4 −4 DESCRIPTION
  2. +57 −10 R/ratebytree.R
  3. +2 −2 man/ratebytree.Rd
View
@@ -1,6 +1,6 @@
Package: phytools
Version: 0.6-35
Date: 2017-10-26
Version: 0.6-36
Date: 2017-10-28
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-10-26 12:00:00 EST
Packaged: 2017-10-28 12:00:00 EST
Repository:
Date/Publication: 2017-10-26 12:00:00 EST
Date/Publication: 2017-10-28 12:00:00 EST
View
@@ -36,12 +36,49 @@ rbt.div<-function(trees,...){
else tol<-1e-12
if(!inherits(trees,"multiPhylo"))
stop("trees should be object of class \"multiPhylo\".")
foo<-if(model=="birth-death") fit.bd else
if(model=="Yule"||model=="yule") fit.yule
fit.multi<-mapply(foo,tree=trees,rho=rho,SIMPLIFY=FALSE)
logL.multi<-sum(sapply(fit.multi,logLik))
t<-lapply(trees,function(phy) sort(branching.times(phy),
decreasing=TRUE))
if(model=="birth-death"){
fit.multi<-mapply(fit.bd,tree=trees,rho=rho,SIMPLIFY=FALSE)
logL.multi<-sum(sapply(fit.multi,logLik))
} else if(model=="equal-extinction"){
lik.eqmu<-function(theta,t,rho,trace=FALSE){
lam<-theta[1:length(t)]
mu<-theta[length(t)+1]
logL<-0
for(i in 1:length(t)) logL<-logL-lik.bd(c(lam[i],mu),t[[i]],rho[i])
if(trace) cat(paste(paste(c(lam,mu,logL),collapse="\t"),"\n",sep=""))
-logL
}
init.b<-sapply(trees,qb)
obj<-nlminb(c(init.b,0),lik.eqmu,t=t,rho=rho,trace=trace,
lower=rep(0,length(trees)+1),upper=rep(Inf,length(trees)+1))
fit.multi<-vector(mode="list",length=length(t))
for(i in 1:length(fit.multi)){
fit.multi[[i]]$b<-obj$par[i]
fit.multi[[i]]$d<-obj$par[length(obj$par)]
}
logL.multi<--obj$objective
} else if(model=="equal-speciation"){
lik.eqlam<-function(theta,t,rho,trace=FALSE){
lam<-theta[1]
mu<-theta[1:length(t)+1]
logL<-0
for(i in 1:length(t)) logL<-logL-lik.bd(c(lam,mu[i]),t[[i]],rho[i])
if(trace) cat(paste(paste(c(lam,mu,logL),collapse="\t"),"\n",sep=""))
-logL
}
init.b<-mean(sapply(trees,qb))
obj<-nlminb(c(init.b,rep(0,length(trees))),lik.eqlam,t=t,rho=rho,
trace=trace,lower=rep(0,length(trees)+1),
upper=rep(Inf,length(trees)+1))
fit.multi<-vector(mode="list",length=length(t))
for(i in 1:length(fit.multi)){
fit.multi[[i]]$b<-obj$par[1]
fit.multi[[i]]$d<-obj$par[i+1]
}
logL.multi<--obj$objective
}
lik.onerate<-function(theta,t,rho,model,trace=FALSE){
logL<-sum(-mapply(lik.bd,t=t,rho=rho,
MoreArgs=list(theta=theta)))
@@ -58,17 +95,22 @@ rbt.div<-function(trees,...){
else rownames(rates.multi)<-paste("tree",1:length(trees),sep="")
colnames(rates.multi)<-c("b","d")
LR<-2*(logL.multi+fit.onerate$objective)
km<-2*length(trees)
k1<-2
km<-if(model=="birth-death") 2*length(trees) else
if(model=="equal-extinction") length(trees)+1 else
if(model=="equal-speciation") length(trees)+1 else
if(model=="Yule") length(trees)
k1<-if(model=="Yule") 1 else 2
P.chisq<-pchisq(LR,df=km-k1,lower.tail=FALSE)
obj<-list(
multi.rate.model=list(
logL=logL.multi,
rates=rates.multi,
k=km,
method="nlminb"),
common.rate.model=list(
logL=-fit.onerate$objective,
rates=setNames(fit.onerate$par,c("b","d")),
k=k1,
method="nlminb"),
model=model,N=length(trees),
n=sapply(trees,Ntip),
@@ -78,6 +120,9 @@ rbt.div<-function(trees,...){
obj
}
## used (for now) to get a starting value for optimization in type="diversification"
qb<-function(tree) (log(Ntip(tree))-log(2))/max(nodeHeights(tree))
## discrete character ratebytree
rbt.disc<-function(trees,x,...){
if(hasArg(trace)) trace<-list(...)$trace
@@ -443,16 +488,18 @@ prbt.div<-function(x,digits=digits){
cat("\n\tb\td\tk\tlog(L)")
cat(paste("\nvalue",round(x$common.rate.model$rates[1],digits),
round(x$common.rate.model$rates[2],digits),
2,round(x$common.rate.model$logL,digits),sep="\t"))
x$common.rate.model$k,round(x$common.rate.model$logL,digits),
sep="\t"))
cat("\n\nML multi diversification-rate model:")
cat(paste("\n",paste(paste("b[",1:x$N,"]",sep=""),collapse="\t"),
paste(paste("d[",1:x$N,"]",sep=""),collapse="\t"),"k\tlog(L)",
sep="\t"))
cat(paste("\nvalue",paste(round(x$multi.rate.model$rates[,"b"],digits),
collapse="\t"),paste(round(x$multi.rate.model$rates[,"d"],digits),
collapse="\t"),2*x$N,round(x$multi.rate.model$logL,digits),
sep="\t"))
cat(paste("\n\nModel fitting method was \"",x$multi.rate.model$method,
collapse="\t"),x$multi.rate.model$k,round(x$multi.rate.model$logL,
digits),sep="\t"))
cat(paste("\n\nDiversification model was \"",x$model,"\".\n",sep=""))
cat(paste("Model fitting method was \"",x$multi.rate.model$method,
"\".\n",sep=""))
cat(paste("\nLikelihood ratio:",round(x$likelihood.ratio,digits),"\n"))
cat(paste("P-value (based on X^2):",round(x$P.chisq,digits),"\n\n"))
View
@@ -9,7 +9,7 @@ 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}. This is not used if \code{type="diversification"}. In the case of \code{posthoc.ratebytree}, an object of class \code{"ratebytree"}.}
\item{...}{optional arguments, including the argument \code{type} (\code{"continuous"}, \code{"discrete"}, or \code{"diversification"}), which, if not specified, the function will attempt to ascertain. For continuous character, optional arguments presently include the following: \code{model}, the model of continuous trait evolution (options are \code{"BM"}, the default, \code{"OU"}, and \code{"EB"}); \code{tol}, used as a minimum value for the fitting rates, to prevent problems in optimization; \code{trace}, a logical value indicating whether or not to report progress in the optimization; \code{test}, the method for hypothesis testing (options are \code{"chisq"} and \code{"simulation"}); \code{quiet}, a logical value indicating whether or not to run perfectly quietly; and \code{se}, a list of vectors containing the standard errors for each value of \code{x}. For \code{type="discrete"} the optional arguments are slightly different. The argument \code{model} can be used, but it must assume the values \code{"ER"}, \code{"SYM"}, \code{"ARD"}, or a numeric matrix following \code{\link{ace}}. Finally, for \code{type="diversification"} only \code{model="birth-death"} is presently available, although I intend to add \code{model="Yule"} at a later date. For \code{type="diversification"} it is also important to consider supplying the sampling fractions, \code{rho}, which is a vector of values between 0 and 1 of the same length as \code{trees}. If not provided the method will assume a sampling fraction of 1.0 for all trees - which is seldom true of empirical studies. Other optional arguments are not yet available for these two types.}
\item{...}{optional arguments, including the argument \code{type} (\code{"continuous"}, \code{"discrete"}, or \code{"diversification"}), which, if not specified, the function will attempt to ascertain. For continuous character, optional arguments presently include the following: \code{model}, the model of continuous trait evolution (options are \code{"BM"}, the default, \code{"OU"}, and \code{"EB"}); \code{tol}, used as a minimum value for the fitting rates, to prevent problems in optimization; \code{trace}, a logical value indicating whether or not to report progress in the optimization; \code{test}, the method for hypothesis testing (options are \code{"chisq"} and \code{"simulation"}); \code{quiet}, a logical value indicating whether or not to run perfectly quietly; and \code{se}, a list of vectors containing the standard errors for each value of \code{x}. For \code{type="discrete"} the optional arguments are slightly different. The argument \code{model} can be used, but it must assume the values \code{"ER"}, \code{"SYM"}, \code{"ARD"}, or a numeric matrix following \code{\link{ace}}. Finally, for \code{type="diversification"} models are so far \code{model="birth-death"}, \code{"equal-extinction"}, and \code{"equal-specation"}, and I intend to add \code{model="Yule"} soon. For \code{type="diversification"} it is also important to consider supplying the sampling fractions, \code{rho}, which is a vector of values between 0 and 1 of the same length as \code{trees}. If not provided the method will assume a sampling fraction of 1.0 for all trees - which is seldom true of empirical studies. Other optional arguments are not yet available for these two types.}
}
\description{
This function essentially implements three different methods for comparing the rate or process of evolution between trees: one for continuously-valued traits, a second for discrete characters, and a third for the rate of diversification (speciation & extinction).
@@ -22,7 +22,7 @@ ratebytree(trees, x, ...)
For discrete traits, the function instead proceeds to fit two variants of the M<i>k</i> model (Lewis 2001): one in which the parameters values (transition rates) of the process are free to vary between trees, and a second in which they are fixed to be the same.
For diversification alone, the function fits two different diversification (speciation & extinction) models (Nee et al. 1994; Stadler 2012): one in which the birth (speciation) and death (extinction) rates are identical between the trees, and a second in which they are permitted to differ.
For diversification alone, the function fits two different diversification (speciation & extinction) models (Nee et al. 1994; Stadler 2012): one in which the birth (speciation) and death (extinction) rates are identical between the trees, and a second in which they are permitted to differ in various ways depending on the value of \code{"model"}.
The method \code{posthoc} conducts a post-hoc comparison of parameter estimates between trees in the multi-rate or multi-process model. The parameter that is compared depends on the fitted model. For instance, in \code{model="BM"} posthoc comparison is made of \code{sig2}; if \code{model="OU"} fitted values of \code{alpha} are compared; and so on. The argument \code{p.adjust.method} can be used to specify a method for adjusting P-values for multiple tests following \code{p.adjust} (defaults to \code{p.adjust.method="none"}.
}

0 comments on commit f8f6b89

Please sign in to comment.