Permalink
Browse files

posthoc method for ratebytree

  • Loading branch information...
liamrevell committed Sep 18, 2017
1 parent d00c87f commit f9c2ac1f9dbf5d6acdba2ec3f6918538af4e8d5e
Showing with 131 additions and 15 deletions.
  1. +4 −4 DESCRIPTION
  2. +1 −0 NAMESPACE
  3. +120 −10 R/ratebytree.R
  4. +6 −1 man/ratebytree.Rd
View
@@ -1,6 +1,6 @@
Package: phytools
Version: 0.6-25
Date: 2017-09-11
Version: 0.6-26
Date: 2017-09-18
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-09-11 12:00:00 EST
Packaged: 2017-09-18 12:00:00 EST
Repository:
Date/Publication: 2017-09-11 12:00:00 EST
Date/Publication: 2017-09-18 12:00:00 EST
View
@@ -105,6 +105,7 @@ S3method(print, ratebytree)
S3method(print, evolvcv.lite)
S3method(print, expand.clade)
S3method(plot, expand.clade)
S3method(posthoc, 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
@@ -131,7 +131,8 @@ rbt.disc<-function(trees,x,...){
index.matrix=fit.multi[[1]]$index.matrix,
states=fit.multi[[1]]$states,
pi=fit.multi[[1]]$pi,
model=model,N=N,likelihood.ratio=LR,P.chisq=P.chisq,
model=model,N=N,n=sapply(trees,Ntip),
likelihood.ratio=LR,P.chisq=P.chisq,
type="discrete")
class(obj)<-"ratebytree"
obj
@@ -170,7 +171,7 @@ rbt.cont<-function(trees,x,...){
x<-mapply(function(x,t) x<-x[t$tip.label],x=x,t=trees,SIMPLIFY=FALSE)
se<-mapply(function(x,t) x<-x[t$tip.label],x=se,t=trees,SIMPLIFY=FALSE)
## first, fit multi-rate model
lik.multi<-function(theta,trees,x,se,model,trace=FALSE){
lik.multi<-function(theta,trees,y,se,model,trace=FALSE){
n<-sapply(trees,Ntip)
N<-length(trees)
sig<-theta[1:N]
@@ -186,7 +187,7 @@ rbt.cont<-function(trees,x,...){
V<-mapply("+",mapply("*",C,sig,SIMPLIFY=FALSE),E,SIMPLIFY=FALSE)
logL<-0
for(i in 1:N)
logL<-logL-t(x[[i]]-a[i])%*%solve(V[[i]])%*%(x[[i]]-
logL<-logL-t(y[[i]]-a[i])%*%solve(V[[i]])%*%(y[[i]]-
a[i])/2-n[i]*log(2*pi)/2-determinant(V[[i]])$modulus[1]/2
if(trace){
cat(paste(paste(round(sig,digits),collapse="\t"),
@@ -225,12 +226,16 @@ rbt.cont<-function(trees,x,...){
paste("r[",1:N,"]",sep="",collapse="\t"),"logL\n",sep="\t"))
}
}
fit.multi<-optim(p,lik.multi,trees=trees,x=x,se=se,model=model,trace=trace,
fit.multi<-optim(p,lik.multi,trees=trees,y=x,se=se,model=model,trace=trace,
method="L-BFGS-B",lower=c(rep(tol,N),rep(-Inf,N),
if(model%in%c("OU","EB")) rep(-Inf,N)),upper=c(rep(Inf,N),rep(Inf,N),
if(model%in%c("OU","EB")) rep(Inf,N)))
## compute covariance matrix
H.multi<-hessian(lik.multi,fit.multi$par,trees=trees,y=x,
se=se,model=model,trace=FALSE)
Cov.multi<-solve(H.multi)
## now fit single-rate model
lik.onerate<-function(theta,trees,x,se,model,trace=FALSE){
lik.onerate<-function(theta,trees,y,se,model,trace=FALSE){
n<-sapply(trees,Ntip)
N<-length(trees)
sig<-theta[1]
@@ -246,7 +251,7 @@ rbt.cont<-function(trees,x,...){
V<-mapply("+",lapply(C,"*",sig),E,SIMPLIFY=FALSE)
logL<-0
for(i in 1:N)
logL<-logL-t(x[[i]]-a[i])%*%solve(V[[i]])%*%(x[[i]]-
logL<-logL-t(y[[i]]-a[i])%*%solve(V[[i]])%*%(y[[i]]-
a[i])/2-n[i]*log(2*pi)/2-determinant(V[[i]])$modulus[1]/2
if(trace){
cat(paste(round(sig,digits),
@@ -257,8 +262,8 @@ rbt.cont<-function(trees,x,...){
}
-logL
}
p<-c(mean(fit.multi$par[1:N]),fit.multi$par[1:N+1])
if(model%in%c("OU","EB")) p<-c(p,mean(fit.multi$par[N+2]))
p<-c(mean(fit.multi$par[1:N]),fit.multi$par[1:N+N])
if(model%in%c("OU","EB")) p<-c(p,mean(fit.multi$par[1:N+2*N]))
if(hasArg(init)){
if(!is.null(init$sigc)) p[1]<-init$sigc
if(!is.null(init$ac)) p[1:N+1]<-init$ac
@@ -277,10 +282,14 @@ rbt.cont<-function(trees,x,...){
cat(paste("sig ","r ","logL\n",sep="\t"))
}
}
fit.onerate<-optim(p,lik.onerate,trees=trees,x=x,se=se,model=model,
fit.onerate<-optim(p,lik.onerate,trees=trees,y=x,se=se,model=model,
trace=trace,method="L-BFGS-B",
lower=c(tol,rep(-Inf,N),if(model=="OU") tol else if(model=="EB") -Inf),
upper=c(Inf,rep(Inf,N),if(model%in%c("OU","EB")) Inf))
## compute covariance matrix
H.onerate<-hessian(lik.onerate,fit.onerate$par,trees=trees,y=x,
se=se,model=model,trace=FALSE)
Cov.onerate<-solve(H.onerate)
## compare models:
LR<-2*(-fit.multi$value+fit.onerate$value)
km<-2*N+if(model=="BM") 0 else if(model%in%c("OU","EB")) N
@@ -317,23 +326,31 @@ rbt.cont<-function(trees,x,...){
}
obj<-list(
multi.rate.model=list(sig2=fit.multi$par[1:N],
SE.sig2=sqrt(diag(Cov.multi)[1:N]),
a=fit.multi$par[1:N+N],
SE.a=sqrt(diag(Cov.multi)[1:N+N]),
alpha=if(model=="OU") fit.multi$par[1:N+2*N] else NULL,
SE.alpha=if(model=="OU") sqrt(diag(Cov.multi)[1:N+2*N]) else NULL,
r=if(model=="EB") fit.multi$par[1:N+2*N] else NULL,
SE.r=if(model=="EB") sqrt(diag(Cov.multi)[1:N+2*N]) else NULL,
k=km,
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],
SE.sig2=sqrt(diag(Cov.onerate)[1]),
a=fit.onerate$par[1:N+1],
SE.a=sqrt(diag(Cov.onerate)[1:N+1]),
alpha=if(model=="OU") fit.onerate$par[N+2] else NULL,
SE.alpha=if(model=="OU") sqrt(diag(Cov.onerate)[N+2]) else NULL,
r=if(model=="EB") fit.onerate$par[N+2] else NULL,
SE.r=if(model=="EB") sqrt(diag(Cov.onerate)[N+2]) else NULL,
k=k1,
logL=-fit.onerate$value,
counts=fit.onerate$counts,convergence=fit.onerate$convergence,
message=fit.onerate$message),
model=model,
N=N,likelihood.ratio=LR,
N=N,n=sapply(trees,Ntip),likelihood.ratio=LR,
P.chisq=if(test=="chisq") P.chisq else NULL,
P.sim=if(test=="simulation") P.sim else NULL,
type="continuous")
@@ -358,6 +375,9 @@ prbt.cont<-function(x,digits=digits){
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",sep="\t"))
cat(paste("SE ",round(x$common.rate.model$SE.sig2,digits),
paste(round(x$common.rate.model$SE.a,digits),collapse="\t"),
"\n\n",sep="\t"))
cat("ML multi-rate model:\n")
cat(paste("\t",paste(paste("s^2[",1:N,"]",sep=""),collapse="\t"),"\t",
@@ -366,6 +386,9 @@ prbt.cont<-function(x,digits=digits){
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",sep="\t"))
cat(paste("SE ",paste(round(x$multi.rate.model$SE.sig2,digits),collapse="\t"),
paste(round(x$multi.rate.model$SE.a,digits),collapse="\t"),
"\n\n",sep="\t"))
} else if(x$model=="OU"){
cat("ML common-regime OU model:\n")
@@ -375,6 +398,10 @@ prbt.cont<-function(x,digits=digits){
paste(round(x$common.rate.model$a,digits),collapse="\t"),
round(x$common.rate.model$alpha,digits),
x$common.rate.model$k,round(x$common.rate.model$logL,digits),
"\n",sep="\t"))
cat(paste("SE ",round(x$common.rate.model$SE.sig2,digits),
paste(round(x$common.rate.model$SE.a,digits),collapse="\t"),
round(x$common.rate.model$SE.alpha,digits),
"\n\n",sep="\t"))
cat("ML multi-regime OU model:\n")
cat(paste("\t",paste(paste("s^2[",1:N,"]",sep=""),collapse="\t"),"\t",
@@ -385,6 +412,10 @@ prbt.cont<-function(x,digits=digits){
paste(round(x$multi.rate.model$a,digits),collapse="\t"),
paste(round(x$multi.rate.model$alpha,digits),collapse="\t"),
x$multi.rate.model$k,round(x$multi.rate.model$logL,digits),
"\n",sep="\t"))
cat(paste("SE ",paste(round(x$multi.rate.model$SE.sig2,digits),collapse="\t"),
paste(round(x$multi.rate.model$SE.a,digits),collapse="\t"),
paste(round(x$multi.rate.model$SE.alpha,digits),collapse="\t"),
"\n\n",sep="\t"))
} else if(x$model=="EB"){
cat("ML common-regime EB model:\n")
@@ -394,6 +425,10 @@ prbt.cont<-function(x,digits=digits){
paste(round(x$common.rate.model$a,digits),collapse="\t"),
round(x$common.rate.model$r,digits),
x$common.rate.model$k,round(x$common.rate.model$logL,digits),
"\n",sep="\t"))
cat(paste("SE ",round(x$common.rate.model$SE.sig2,digits),
paste(round(x$common.rate.model$SE.a,digits),collapse="\t"),
round(x$common.rate.model$SE.r,digits),
"\n\n",sep="\t"))
cat("ML multi-regime EB model:\n")
cat(paste("\t",paste(paste("s^2[",1:N,"]",sep=""),collapse="\t"),"\t",
@@ -404,6 +439,10 @@ prbt.cont<-function(x,digits=digits){
paste(round(x$multi.rate.model$a,digits),collapse="\t"),
paste(round(x$multi.rate.model$r,digits),collapse="\t"),
x$multi.rate.model$k,round(x$multi.rate.model$logL,digits),
"\n",sep="\t"))
cat(paste("SE ",paste(round(x$multi.rate.model$SE.sig2,digits),collapse="\t"),
paste(round(x$multi.rate.model$SE.a,digits),collapse="\t"),
paste(round(x$multi.rate.model$SE.r,digits),collapse="\t"),
"\n\n",sep="\t"))
}
cat(paste("Likelihood ratio:",round(x$likelihood.ratio,digits),"\n"))
@@ -443,3 +482,74 @@ prbt.disc<-function(x,digits=digits){
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"))
}
## posthoc comparison S3 method
posthoc<-function(x, ...) UseMethod("posthoc")
posthoc.ratebytree<-function(x,...){
if(hasArg(p.adjust.method)) p.adjust.method<-list(...)$p.adjust.method
else p.adjust.method<-"none"
if(x$type!="continuous"){
cat("Sorry. No posthoc method yet implemented for this data type.\n\n")
} else {
if(x$model=="BM") k<-2
else if(x$model%in%c("EB","OU")) k<-3
t<-df<-P<-matrix(0,x$N,x$N)
for(i in 1:x$N){
for(j in 1:x$N){
x1<-if(x$model=="BM") x$multi.rate.model$sig2[i]
else if(x$model=="OU") x$multi.rate.model$alpha[i]
else if(x$model=="EB") x$multi.rate.model$r[i]
x2<-if(x$model=="BM") x$multi.rate.model$sig2[j]
else if(x$model=="OU") x$multi.rate.model$alpha[j]
else if(x$model=="EB") x$multi.rate.model$r[j]
s1<-x$multi.rate.model$SE.sig2[i]^2
s2<-x$multi.rate.model$SE.sig2[j]^2
n1<-x$n[i]
n2<-x$n[j]
se<-sqrt(s1+s2)
df[i,j]<-(s1/n1+s2/n2)^2/
((s1/n1)^2/(n1-k)+(s2/n2)^2/(n2-k))
t[i,j]<-(x1-x2)/se
P[i,j]<-2*pt(abs(t[i,j]),df[i,j],lower.tail=FALSE)
p<-P[upper.tri(P)]
p<-p.adjust(p,method=p.adjust.method)
P[upper.tri(P)]<-p
P[lower.tri(P)]<-p
}
}
}
obj<-list(t=t,df=df,P=P,model=x$model,type=x$type,
p.adjust.method=p.adjust.method)
class(obj)<-"posthoc.ratebytree"
obj
}
print.posthoc.ratebytree<-function(x,...){
if(hasArg(digits)) digits<-list(...)$digits
else digits<-4
t<-x$t[upper.tri(x$t)]
df<-x$df[upper.tri(x$df)]
P<-x$P[upper.tri(x$P)]
N<-nrow(x$P)
paste("tree",1:(N-1),"vs.",2:N)
nn<-vector(mode="character",length=N*(N-1)/2)
k<-1
for(i in 1:N) for(j in i:N){
if(i!=j){
nn[k]<-paste("tree",i," vs. ",j,sep="")
k<-k+1
}
}
X<-data.frame(t=t,df=df,P=P)
rownames(X)<-nn
cat(paste("\nPost-hoc test for \"",x$model,"\" model.\n",sep=""))
cat(paste("(Comparison is of estimated values of",
if(x$model=="BM") "sigma^2.)\n\n"
else if(x$model=="OU") "alpha.)\n\n"
else if(x$model=="EB") "r.)\n\n"))
print(round(X,digits))
cat(paste("\nP-values adjusted using method=\"",x$p.adjust.method,
"\".\n\n",sep=""))
}
View
@@ -1,12 +1,15 @@
\name{ratebytree}
\alias{ratebytree}
\alias{posthoc.ratebytree}
\title{Likelihood test for rate variation among trees}
\usage{
ratebytree(trees, x, ...)
posthoc.ratebytree(x, ...)
\method{posthoc}{ratebytree}(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{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}. In the case of \code{posthoc.ratebytree}, an object of class \code{"ratebytree"}.}
\item{...}{optional arguments, including the argument \code{type} (\code{"continuous"} or \code{"discrete"}), 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}}. Other optional arguments are not yet available.}
}
\description{
@@ -19,6 +22,8 @@ ratebytree(trees, x, ...)
The latter model corresponds to an extension 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).
For discrete traits, the function instead proceeds to fit two variants of the M<i>k</i> model: 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.
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"}.
}
\details{
At present it is not possible to specify different models to fit for the different trees - although if (for instance) character evolution on tree 1 proceeded by a strong \emph{OU} process while character evolution on tree 2 was by \emph{BM}, we would probably reject a constant-process model and tree 2 should show a very low value of \code{alpha}.

0 comments on commit f9c2ac1

Please sign in to comment.