Nothing
R/ratebytree.R
In phytools: Phylogenetic Tools for Comparative Biology (and Other Things)
Defines functions
AIC.ratebytree
print.posthoc.ratebytree
posthoc.ratebytree
posthoc
prbt.disc
prbt.cont
prbt.div
print.ratebytree
rbt.cont
rbt.disc
rbt.div
ratebytree
Documented in
posthoc
posthoc.ratebytree
ratebytree
## method to compare the rate of evolution for a character between trees
## continuous character method closely related to 'censored' approach of O'Meara et al. (2006; Evolution)
## discrete character method fits Mk model of Lewis 2001
## diversification method fits Yule or birth-death model of Nee et al. (1994) & Stadler (2012)
## written by Liam J. Revell 2017
ratebytree<-function(trees,x,...){
if(hasArg(type)) type<-list(...)$type
else if(!missing(x)&&!is.null(x)) {
if(is.factor(unlist(x))||is.character(unlist(x)))
type<-"discrete"
else type<-"continuous"
} else type<-"diversification"
if(type=="continuous") obj<-rbt.cont(trees,x,...)
else if(type=="discrete") obj<-rbt.disc(trees,x,...)
else if(type=="diversification") obj<-rbt.div(trees,...)
else {
cat(paste("type =",type,"not recognized.\n"))
obj<-NULL
}
obj
}
rbt.div<-function(trees,...){
if(hasArg(trace)) trace<-list(...)$trace
else trace<-FALSE
if(hasArg(digits)) digits<-list(...)$digits
else digits<-4
if(hasArg(test)) test<-list(...)$test
else test<-"chisq"
if(hasArg(quiet)) quiet<-list(...)$quiet
else quiet<-FALSE
if(hasArg(model)) model<-list(...)$model
else model<-"birth-death"
if(hasArg(rho)) rho<-list(...)$rho
else rho<-rep(1,length(trees))
if(hasArg(tol)) tol<-list(...)$tol
else tol<-1e-12
if(hasArg(iter)) iter<-list(...)$iter
else iter<-10
if(!inherits(trees,"multiPhylo"))
stop("trees should be object of class \"multiPhylo\".")
if(any(!sapply(trees,is.ultrametric))){
cat("One or more trees fails check is.ultrametric.\n")
cat("If you believe your tree to be ultrametric ")
cat("use force.ultrametric.\n")
stop()
}
t<-lapply(trees,function(phy) sort(branching.times(phy),
decreasing=TRUE))
if(model=="birth-death"){
fit.multi<-mapply(fit.bd,tree=trees,rho=rho,iter=iter,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))
count<-0
while(!is.finite(obj$objective)&&count<iter){
obj<-nlminb(runif(n=2,0,2)*c(init.b,0),lik.eqmu,t=t,rho=rho,
trace=trace,lower=rep(0,length(trees)+1),
upper=rep(Inf,length(trees)+1))
count<-count+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))
print(obj)
count<-0
while(!is.finite(obj$objective)&&count<iter){
obj<-nlminb(runif(n=2,0,2)*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))
count<-count+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
} else if(model=="Yule"){
interval<-if(hasArg(interval)) list(...)$interval else c(0,10*mean(sapply(trees,qb)))
fit.multi<-mapply(fit.yule,tree=trees,rho=rho,MoreArgs=list(interval=interval),
SIMPLIFY=FALSE)
logL.multi<-sum(sapply(fit.multi,logLik))
}
if(model!="Yule"){
lik.onerate<-function(theta,t,rho,model,trace=FALSE){
logL<-sum(-mapply(lik.bd,t=t,rho=rho,
MoreArgs=list(theta=theta)))
if(trace) cat(paste(theta[1],"\t",theta[2],"\t",logL,"\n"))
-logL
}
init.b<-mean(sapply(fit.multi,function(x) x$b))
init.d<-mean(sapply(fit.multi,function(x) x$d))
fit.onerate<-nlminb(c(init.b,init.d),lik.onerate,t=t,rho=rho,
model=model,trace=trace,lower=c(0,0),upper=rep(Inf,2))
count<-0
while(!is.finite(fit.onerate$objective)&&count<iter){
fit.onerate<-nlminb(runif(n=2,0,2)*c(init.b,init.d),lik.onerate,
t=t,rho=rho,model=model,trace=trace,lower=c(0,0),
upper=rep(Inf,2))
count<-count+1
}
rates.multi<-cbind(sapply(fit.multi,function(x) x$b),
sapply(fit.multi,function(x) x$d))
} else {
lik.yule<-function(theta,t,rho,trace=FALSE){
logL<-sum(-mapply(lik.bd,t=t,rho=rho,MoreArgs=list(theta=c(theta,0))))
if(trace) cat(paste(theta,logL,"\n",sep="\t"))
-logL
}
obj<-optimize(lik.yule,interval,t=t,rho=rho,trace=trace)
TOL<-diff(interval)*1e-3
convergence<-if((obj$minimum-TOL)<interval[1]||(obj$minimum+TOL)>interval[2]) 52 else 0
fit.onerate<-list(par=c(obj$minimum,0),objective=obj$objective,convergence=convergence,
message=if(convergence!=0) "Estimate at may be at limits of interval." else
"Probably converged.")
rates.multi<-cbind(sapply(fit.multi,function(x) x$b),
rep(0,length(trees)))
}
if(!is.null(names(trees))) rownames(rates.multi)<-names(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<-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=if(model=="Yule") "optimize" else "nlminb"),
common.rate.model=list(
logL=-fit.onerate$objective,
rates=setNames(fit.onerate$par,c("b","d")),
k=k1,
method=if(model=="Yule") "optimize" else "nlminb"),
model=model,N=length(trees),
n=sapply(trees,Ntip),
likelihood.ratio=LR,P.chisq=P.chisq,
type="diversification")
class(obj)<-"ratebytree"
obj
}
## discrete character ratebytree
rbt.disc<-function(trees,x,...){
if(hasArg(trace)) trace<-list(...)$trace
else trace<-FALSE
if(hasArg(digits)) digits<-list(...)$digits
else digits<-4
if(hasArg(test)) test<-list(...)$test
else test<-"chisq"
if(hasArg(quiet)) quiet<-list(...)$quiet
else quiet<-FALSE
if(hasArg(model)) model<-list(...)$model
else model<-"ER"
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)
fit.multi<-mapply(fitMk,tree=trees,x=x,MoreArgs=list(model=model),
SIMPLIFY=FALSE)
logL.multi<-sum(sapply(fit.multi,logLik))
lik.onerate<-function(theta,trees,x,model,trace=FALSE){
ss<-sort(unique(unlist(x)))
m<-length(ss)
if(is.character(model)){
rate<-matrix(NA,m,m)
if(model=="ER"){
k<-rate[]<-1
diag(rate)<-NA
}
else if(model=="ARD") {
k<-m*(m-1)
rate[col(rate)!=row(rate)]<-1:k
}
else if(model=="SYM") {
k<-m*(m-1)/2
ii<-col(rate)<row(rate)
rate[ii]<-1:k
rate<-t(rate)
rate[ii]<-1:k
}
} else {
rate<-model
k<-max(rate)
}
Q <- matrix(0, m, m)
index.matrix <- rate
tmp <- cbind(1:m, 1:m)
rate[tmp] <- 0
rate[rate == 0] <- k + 1
Q[]<-c(theta,0)[rate]
diag(Q)<--rowSums(Q)
colnames(Q)<-rownames(Q)<-ss
fit.onerate<-mapply(fitMk,tree=trees,x=x,MoreArgs=list(fixedQ=Q),
SIMPLIFY=FALSE)
-sum(sapply(fit.onerate,logLik))
}
pp<-sapply(fit.multi,function(x) x$rates)
pp<-if(is.matrix(pp)) rowMeans(pp) else mean(pp)
if(length(pp)>1){
fit.onerate<-optim(pp,lik.onerate,trees=trees,x=x,model=model)
} else {
fit.onerate<-optimize(lik.onerate,c(0,1000*pp),trees=trees,x=x,
model=model)
names(fit.onerate)<-c("par","value")
}
rates.multi<-t(sapply(fit.multi,function(x) x$rates))
if(is.character(model)){
m<-length(fit.multi[[1]]$states)
if(model=="ER"){
rates.multi<-t(rates.multi)
colnames(rates.multi)<-"q"
} else if(model=="SYM"){
if(length(fit.multi[[1]]$states)==2) rates.multi<-t(rates.multi)
colnames(rates.multi)<-
sapply(fit.multi[[1]]$states,function(x,y)
sapply(y,function(y,x) paste(x,"<->",y,sep=""),x=x),
y=fit.multi[[1]]$states)[lower.tri(matrix(0,m,m))]
} else if(model=="ARD"){
ii<-(upper.tri(matrix(0,m,m))+lower.tri(matrix(0,m,m))==TRUE)
colnames(rates.multi)<-
sapply(fit.multi[[1]]$states,function(x,y)
sapply(y,function(y,x) paste(y,"->",x,sep=""),x=x),
y=fit.multi[[1]]$states)[ii]
}
} else {
k<-max(fit.multi[[1]]$index.matrix,na.rm=TRUE)
if(k==1) rates.multi<-t(rates.multi)
foo<-function(i,index.matrix,states){
rc<-which(index.matrix==i,arr.ind=TRUE)
lab<-apply(rc,1,function(ind,ss) paste(ss[ind],collapse="->"),
ss=states)
if(length(lab)>1) paste(lab,collapse=",") else lab
}
colnames(rates.multi)<-sapply(1:k,
foo,fit.multi[[1]]$index.matrix,fit.multi[[1]]$states)
}
if(!is.null(names(trees))) rownames(rates.multi)<-names(trees)
else rownames(rates.multi)<-paste("tree",1:length(trees),sep="")
LR<-2*(logL.multi+fit.onerate$value)
km<-sum(sapply(fit.multi,function(x) max(x$index.matrix,na.rm=TRUE)))
k1<-length(pp)
P.chisq<-pchisq(LR,df=km-k1,lower.tail=FALSE)
obj<-list(
multi.rate.model=list(
logL=logL.multi,
rates=rates.multi,
method=fit.multi[[1]]$method),
common.rate.model=list(
logL=-fit.onerate$value,
rates=setNames(fit.onerate$par,colnames(rates.multi)),
method=if(length(pp)>1) "optim" else "optimize"),
index.matrix=fit.multi[[1]]$index.matrix,
states=fit.multi[[1]]$states,
pi=fit.multi[[1]]$pi,
model=model,N=N,n=sapply(trees,Ntip),
likelihood.ratio=LR,P.chisq=P.chisq,
type="discrete")
class(obj)<-"ratebytree"
obj
}
## continuous character ratebytree
rbt.cont<-function(trees,x,...){
if(hasArg(tol)) tol<-list(...)$tol
else tol<-1e-8
if(hasArg(trace)) trace<-list(...)$trace
else trace<-FALSE
if(hasArg(digits)) digits<-list(...)$digits
else digits<-4
if(hasArg(test)) test<-list(...)$test
else test<-"chisq"
if(hasArg(quiet)) quiet<-list(...)$quiet
else quiet<-FALSE
if(hasArg(maxit)) maxit<-list(...)$maxit
else maxit<-500
if(hasArg(regimes)){
regimes<-list(...)$regimes
if(!is.factor(regimes)) regimes<-as.factor(regimes)
} else regimes<-as.factor(1:length(trees))
if(hasArg(model)) model<-list(...)$model
else model<-"BM"
if(!(model%in%c("BM","OU","EB"))){
cat(paste("model =",model,"not recognized. using model = \"BM\"\n"))
model<-"BM"
}
## 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.")
if(hasArg(se)) se<-list(...)$se
else {
se<-x
for(i in 1:length(x)) se[[i]][1:length(se[[i]])]<-0
}
N<-length(trees)
m<-length(levels(regimes))
## reorder the trait vectors in x & SEs in se
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)
## likelihood functions
lik.multi<-function(theta,trees,y,se,model,regimes,trace=FALSE){
m<-length(unique(regimes))
ind<-setNames(lapply(levels(regimes),function(x,y) which(y==x),y=regimes),
levels(regimes))
THETA<-list()
for(i in 1:m){
THETA[[i]]<-c(theta[i],theta[ind[[i]]+m])
if(model%in%c("OU","EB")) THETA[[i]]<-c(THETA[[i]],theta[i+m+N])
}
logL<-0
for(i in 1:m)
logL<-logL-lik.onerate(THETA[[i]],trees[ind[[i]]],y[ind[[i]]],se[ind[[i]]],
model=model,trace=FALSE)
if(trace){
cat(paste(paste(round(theta[1:m],digits),collapse="\t"),
if(model=="OU") paste(round(theta[1:m+m+N],digits),collapse="\t"),
if(model=="EB") paste(round(theta[1:m+m+N],digits),collapse="\t"),
round(logL,digits),"\n",sep="\t"))
flush.console()
}
-logL
}
lik.onerate<-function(theta,trees,y,se,model,trace=FALSE){
n<-sapply(trees,Ntip)
N<-length(trees)
sig<-theta[1]
a<-theta[1:N+1]
if(model=="OU") alpha<-theta[N+2]
if(model=="EB") r<-theta[N+2]
if(model=="BM") C<-lapply(trees,vcv)
else if(model=="OU") C<-lapply(trees,vcvPhylo,model="OU",alpha=alpha,
anc.nodes=FALSE)
else if(model=="EB") C<-lapply(trees,vcvPhylo,model="EB",r=r,
anc.nodes=FALSE)
E<-lapply(se,diag)
V<-mapply("+",lapply(C,"*",sig),E,SIMPLIFY=FALSE)
logL<-0
for(i in 1:N)
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),
if(model=="OU") round(alpha,digits),
if(model=="EB") round(r,digits),
round(logL,digits),"\n",sep="\t"))
flush.console()
}
-logL
}
## first, fit multi-rate model
f1<-function(tree,x){
pvcv<-phyl.vcv(as.matrix(x),vcv(tree),1)
c(pvcv$R[1,1],pvcv$a[1,1])
}
P<-mapply(f1,trees,x,SIMPLIFY=FALSE)
PP<-list()
PP[[1]]<-vector()
for(i in 1:m)
PP[[1]][i]<-mean(sapply(P,function(x) x[1])[which(regimes==levels(regimes)[i])])
PP[[2]]<-sapply(P,function(x) x[2])
if(model%in%c("OU","EB")) PP[[3]]<-rep(0,m)
if(hasArg(init)){
init<-list(...)$init
if(!is.null(init$sigm)) PP[[1]]<-init$sigm
if(!is.null(init$am)) PP[[2]]<-init$am
if(model=="OU") if(!is.null(init$alpham)) PP[[3]]<-init$alpham
if(model=="EB") if(!is.null(init$rm)) PP[[3]]<-init$rm
}
p<-unlist(PP)
if(trace){
if(model=="BM"){
cat("\nOptimizing multi-rate model....\n")
cat(paste(paste("sig[",1:m,"]",sep="",collapse="\t"),"logL\n",sep="\t"))
} else if(model=="OU"){
cat("\nOptimizing multi-regime model....\n")
cat(paste(paste("sig[",1:m,"]",sep="",collapse="\t"),
paste("alpha[",1:m,"]",sep="",collapse="\t"),"logL\n",sep="\t"))
} else if(model=="EB"){
cat("\nOptimizing multi-regime model....\n")
cat(paste(paste("sig[",1:m,"]",sep="",collapse="\t"),
paste("r[",1:m,"]",sep="",collapse="\t"),"logL\n",sep="\t"))
}
}
fit.multi<-optim(p,lik.multi,trees=trees,y=x,se=se,model=model,
regimes=regimes,trace=trace,
method="L-BFGS-B",lower=c(rep(tol,m),rep(-Inf,N),
if(model%in%c("OU","EB")) rep(-Inf,m)),upper=c(rep(Inf,m),rep(Inf,N),
if(model%in%c("OU","EB")) rep(Inf,m)),control=list(maxit=maxit))
## compute covariance matrix
H.multi<-hessian(lik.multi,fit.multi$par,trees=trees,y=x,
se=se,model=model,regimes=regimes,trace=FALSE)
Cov.multi<-if(qr(H.multi)$rank!=ncol(H.multi)) ginv(H.multi) else solve(H.multi)
## now fit single-rate model
p<-c(mean(fit.multi$par[1:m]),fit.multi$par[1:N+m])
if(model%in%c("OU","EB")) p<-c(p,mean(fit.multi$par[1:m+m+N]))
if(hasArg(init)){
if(!is.null(init$sigc)) p[1]<-init$sigc
if(!is.null(init$ac)) p[1:N+1]<-init$ac
if(model=="OU") if(!is.null(init$alphac)) p[N+2]<-init$alphac
if(model=="EB") if(!is.null(init$rc)) p[N+2]<-init$rc
}
if(trace){
if(model=="BM"){
cat("\nOptimizing common-rate model....\n")
cat(paste("sig ","logL\n",sep="\t"))
} else if(model=="OU"){
cat("\nOptimizing common-regime model....\n")
cat(paste("sig ","alpha ","logL\n",sep="\t"))
} else if(model=="EB"){
cat("\nOptimizing common-regime mode.....\n")
cat(paste("sig ","r ","logL\n",sep="\t"))
}
}
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),
control=list(maxit=maxit))
## compute covariance matrix
H.onerate<-hessian(lik.onerate,fit.onerate$par,trees=trees,y=x,
se=se,model=model,trace=FALSE)
Cov.onerate<-if(qr(H.onerate)$rank!=ncol(H.onerate)) ginv(H.onerate) else solve(H.onerate)
## compare models:
LR<-2*(-fit.multi$value+fit.onerate$value)
km<-N+m+if(model=="BM") 0 else if(model%in%c("OU","EB")) m
k1<-N+if(model=="BM") 1 else if(model%in%c("OU","EB")) 2
if(test=="simulation"&&model%in%c("OU","EB")){
cat("Simulation test not yet available for chosen model. Using chi-square test.\n")
test<-"chisq"
}
if(test=="chisq") P.chisq<-pchisq(LR,df=km-k1,lower.tail=FALSE)
else if(test=="simulation"){
if(!quiet) cat("Generating null distribution via simulation -> |")
flush.console()
if(hasArg(nsim)) nsim<-list(...)$nsim
else nsim<-100
X<-mapply(fastBM,tree=trees,a=as.list(fit.onerate$par[1:N+1]),
MoreArgs=list(sig2=fit.onerate$par[1],nsim=nsim),
SIMPLIFY=FALSE)
P.sim<-1/(nsim+1)
pct<-0.1
for(i in 1:nsim){
x.sim<-lapply(X,function(x,ind) x[,ind],ind=i)
f2<-function(x,se) sampleFrom(xbar=x,xvar=se^2,n=rep(1,length(x)))
x.sim<-mapply(f2,x=x.sim,se=se,SIMPLIFY=FALSE)
fit.sim<-ratebytree(trees,x.sim,se=se)
P.sim<-P.sim+(fit.sim$likelihood.ratio>=LR)/(nsim+1)
if(i/nsim>=pct){
if(!quiet) cat(".")
flush.console()
pct<-pct+0.1
}
}
if(!quiet) cat(".|\nDone!\n")
flush.console()
}
obj<-list(
multi.rate.model=list(sig2=fit.multi$par[1:m],
SE.sig2=sqrt(diag(Cov.multi)[1:m]),
a=fit.multi$par[1:N+m],
SE.a=sqrt(diag(Cov.multi)[1:N+m]),
alpha=if(model=="OU") fit.multi$par[1:m+m+N] else NULL,
SE.alpha=if(model=="OU") sqrt(diag(Cov.multi)[1:m+m+N]) else NULL,
r=if(model=="EB") fit.multi$par[1:N+m+N] else NULL,
SE.r=if(model=="EB") sqrt(diag(Cov.multi)[1:N+m+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,n=sapply(trees,Ntip),likelihood.ratio=LR,
regimes=regimes,m=m,
P.chisq=if(test=="chisq") P.chisq else NULL,
P.sim=if(test=="simulation") P.sim else NULL,
type="continuous")
class(obj)<-"ratebytree"
obj
}
## S3 print method for ratebytree
print.ratebytree<-function(x,...){
if(hasArg(digits)) digits<-list(...)$digits
else digits<-4
if(x$type=="continuous") prbt.cont(x,digits=digits)
else if(x$type=="discrete") prbt.disc(x,digits=digits)
else if(x$type=="diversification") prbt.div(x,digits=digits)
else print(x)
}
prbt.div<-function(x,digits=digits){
cat("ML common diversification-rate model:")
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),
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"),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"))
}
prbt.cont<-function(x,digits=digits){
N<-x$N
m<-x$m
if(x$model=="BM"){
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",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[",levels(x$regimes),"]",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",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")
cat(paste("\ts^2\t",paste(paste("a[",1:N,"]",sep=""),collapse="\t")),
"\talpha\tk\tlogL\n")
cat(paste("value",round(x$common.rate.model$sig2,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[",levels(x$regimes),"]",sep=""),collapse="\t"),"\t",
paste(paste("a[",1:N,"]",sep=""),collapse="\t"),"\t",
paste(paste("alp[",levels(x$regimes),"]",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"),
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")
cat(paste("\ts^2\t",paste(paste("a[",1:N,"]",sep=""),collapse="\t")),
"\tr\tk\tlogL\n")
cat(paste("value",round(x$common.rate.model$sig2,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[",levels(x$regimes),"]",sep=""),collapse="\t"),"\t",
paste(paste("a[",1:N,"]",sep=""),collapse="\t"),"\t",
paste(paste("r[",levels(x$regimes),"]",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"),
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"))
if(!is.null(x$P.chisq))
cat(paste("P-value (based on X^2):",round(x$P.chisq,digits),"\n\n"))
else if(!is.null(x$P.sim))
cat(paste("P-value (based on simulation):",round(x$P.sim,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")
}
prbt.disc<-function(x,digits=digits){
cat("ML common-rate model:\n")
cat(paste("\t",paste(names(x$common.rate.model$rates),collapse="\t"),
"\tk\tlogL\n",sep=""))
cat(paste("value\t",paste(round(x$common.rate.model$rates,digits),
collapse="\t"),"\t",max(x$index.matrix,na.rm=T),"\t",
round(x$common.rate.model$logL,digits),"\n",sep=""))
cat(paste("\nModel fitting method was \"",x$common.rate.model$method,
"\".\n",sep=""))
cat("\nML multi-rate model:\n")
cat(paste("\t",paste(colnames(x$multi.rate.model$rates),collapse="\t"),
"\tk\tlogL\n",sep=""))
for(i in 1:nrow(x$multi.rate.model$rates)){
if(i>1) cat(paste(rownames(x$multi.rate.model$rates)[i],"\t",
paste(round(x$multi.rate.model$rates[i,],
digits),collapse="\t"),"\n",sep=""))
else if(i==1) cat(paste(rownames(x$multi.rate.model$rates)[i],"\t",
paste(round(x$multi.rate.model$rates[i,],
digits),collapse="\t"),"\t",x$N*max(x$index.matrix,na.rm=T),"\t",
round(x$multi.rate.model$logL,digits),"\n",sep=""))
}
cat(paste("\nModel 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"))
}
## 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$m,x$m)
for(i in 1:x$m){
for(j in 1:x$m){
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("reg.",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("reg.",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=""))
}
AIC.ratebytree<-function(object,...,k=2){
aic<-data.frame(AIC=c(k*object$common.rate.model$k-2*object$common.rate.model$logL,
k*object$multi.rate.model$k-2*object$multi.rate.model$logL),
df=c(object$common.rate.model$k,object$multi.rate.model$k))
model.names<-if(is.null(object$model)) 1 else object$model
addtl.obj<-list(...)
if(length(addtl.obj)>0){
for(i in 1:length(addtl.obj)){
aic<-rbind(aic,c(k*addtl.obj[[i]]$multi.rate.model$k-
2*addtl.obj[[i]]$multi.rate.model$logL,
addtl.obj[[i]]$multi.rate.model$k))
model.names<-c(model.names,
if(is.null(addtl.obj[[i]]$model)) 1+i else addtl.obj[[i]]$model)
}
rownames(aic)<-c("common-rate",paste("multi-rate:",model.names,sep=""))
} else rownames(aic)<-c("common-rate","multi-rate")
aic
}
Try the phytools package in your browser
Any scripts or data that you put into this service are public.
phytools documentation built on June 22, 2024, 10:39 a.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions AIC.ratebytree print.posthoc.ratebytree posthoc.ratebytree posthoc prbt.disc prbt.cont prbt.div print.ratebytree rbt.cont rbt.disc rbt.div ratebytree
Documented in posthoc posthoc.ratebytree ratebytree
## method to compare the rate of evolution for a character between trees
## continuous character method closely related to 'censored' approach of O'Meara et al. (2006; Evolution)
## discrete character method fits Mk model of Lewis 2001
## diversification method fits Yule or birth-death model of Nee et al. (1994) & Stadler (2012)
## written by Liam J. Revell 2017
ratebytree<-function(trees,x,...){
if(hasArg(type)) type<-list(...)$type
else if(!missing(x)&&!is.null(x)) {
if(is.factor(unlist(x))||is.character(unlist(x)))
type<-"discrete"
else type<-"continuous"
} else type<-"diversification"
if(type=="continuous") obj<-rbt.cont(trees,x,...)
else if(type=="discrete") obj<-rbt.disc(trees,x,...)
else if(type=="diversification") obj<-rbt.div(trees,...)
else {
cat(paste("type =",type,"not recognized.\n"))
obj<-NULL
}
obj
}
rbt.div<-function(trees,...){
if(hasArg(trace)) trace<-list(...)$trace
else trace<-FALSE
if(hasArg(digits)) digits<-list(...)$digits
else digits<-4
if(hasArg(test)) test<-list(...)$test
else test<-"chisq"
if(hasArg(quiet)) quiet<-list(...)$quiet
else quiet<-FALSE
if(hasArg(model)) model<-list(...)$model
else model<-"birth-death"
if(hasArg(rho)) rho<-list(...)$rho
else rho<-rep(1,length(trees))
if(hasArg(tol)) tol<-list(...)$tol
else tol<-1e-12
if(hasArg(iter)) iter<-list(...)$iter
else iter<-10
if(!inherits(trees,"multiPhylo"))
stop("trees should be object of class \"multiPhylo\".")
if(any(!sapply(trees,is.ultrametric))){
cat("One or more trees fails check is.ultrametric.\n")
cat("If you believe your tree to be ultrametric ")
cat("use force.ultrametric.\n")
stop()
}
t<-lapply(trees,function(phy) sort(branching.times(phy),
decreasing=TRUE))
if(model=="birth-death"){
fit.multi<-mapply(fit.bd,tree=trees,rho=rho,iter=iter,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))
count<-0
while(!is.finite(obj$objective)&&count<iter){
obj<-nlminb(runif(n=2,0,2)*c(init.b,0),lik.eqmu,t=t,rho=rho,
trace=trace,lower=rep(0,length(trees)+1),
upper=rep(Inf,length(trees)+1))
count<-count+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))
print(obj)
count<-0
while(!is.finite(obj$objective)&&count<iter){
obj<-nlminb(runif(n=2,0,2)*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))
count<-count+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
} else if(model=="Yule"){
interval<-if(hasArg(interval)) list(...)$interval else c(0,10*mean(sapply(trees,qb)))
fit.multi<-mapply(fit.yule,tree=trees,rho=rho,MoreArgs=list(interval=interval),
SIMPLIFY=FALSE)
logL.multi<-sum(sapply(fit.multi,logLik))
}
if(model!="Yule"){
lik.onerate<-function(theta,t,rho,model,trace=FALSE){
logL<-sum(-mapply(lik.bd,t=t,rho=rho,
MoreArgs=list(theta=theta)))
if(trace) cat(paste(theta[1],"\t",theta[2],"\t",logL,"\n"))
-logL
}
init.b<-mean(sapply(fit.multi,function(x) x$b))
init.d<-mean(sapply(fit.multi,function(x) x$d))
fit.onerate<-nlminb(c(init.b,init.d),lik.onerate,t=t,rho=rho,
model=model,trace=trace,lower=c(0,0),upper=rep(Inf,2))
count<-0
while(!is.finite(fit.onerate$objective)&&count<iter){
fit.onerate<-nlminb(runif(n=2,0,2)*c(init.b,init.d),lik.onerate,
t=t,rho=rho,model=model,trace=trace,lower=c(0,0),
upper=rep(Inf,2))
count<-count+1
}
rates.multi<-cbind(sapply(fit.multi,function(x) x$b),
sapply(fit.multi,function(x) x$d))
} else {
lik.yule<-function(theta,t,rho,trace=FALSE){
logL<-sum(-mapply(lik.bd,t=t,rho=rho,MoreArgs=list(theta=c(theta,0))))
if(trace) cat(paste(theta,logL,"\n",sep="\t"))
-logL
}
obj<-optimize(lik.yule,interval,t=t,rho=rho,trace=trace)
TOL<-diff(interval)*1e-3
convergence<-if((obj$minimum-TOL)<interval[1]||(obj$minimum+TOL)>interval[2]) 52 else 0
fit.onerate<-list(par=c(obj$minimum,0),objective=obj$objective,convergence=convergence,
message=if(convergence!=0) "Estimate at may be at limits of interval." else
"Probably converged.")
rates.multi<-cbind(sapply(fit.multi,function(x) x$b),
rep(0,length(trees)))
}
if(!is.null(names(trees))) rownames(rates.multi)<-names(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<-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=if(model=="Yule") "optimize" else "nlminb"),
common.rate.model=list(
logL=-fit.onerate$objective,
rates=setNames(fit.onerate$par,c("b","d")),
k=k1,
method=if(model=="Yule") "optimize" else "nlminb"),
model=model,N=length(trees),
n=sapply(trees,Ntip),
likelihood.ratio=LR,P.chisq=P.chisq,
type="diversification")
class(obj)<-"ratebytree"
obj
}
## discrete character ratebytree
rbt.disc<-function(trees,x,...){
if(hasArg(trace)) trace<-list(...)$trace
else trace<-FALSE
if(hasArg(digits)) digits<-list(...)$digits
else digits<-4
if(hasArg(test)) test<-list(...)$test
else test<-"chisq"
if(hasArg(quiet)) quiet<-list(...)$quiet
else quiet<-FALSE
if(hasArg(model)) model<-list(...)$model
else model<-"ER"
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)
fit.multi<-mapply(fitMk,tree=trees,x=x,MoreArgs=list(model=model),
SIMPLIFY=FALSE)
logL.multi<-sum(sapply(fit.multi,logLik))
lik.onerate<-function(theta,trees,x,model,trace=FALSE){
ss<-sort(unique(unlist(x)))
m<-length(ss)
if(is.character(model)){
rate<-matrix(NA,m,m)
if(model=="ER"){
k<-rate[]<-1
diag(rate)<-NA
}
else if(model=="ARD") {
k<-m*(m-1)
rate[col(rate)!=row(rate)]<-1:k
}
else if(model=="SYM") {
k<-m*(m-1)/2
ii<-col(rate)<row(rate)
rate[ii]<-1:k
rate<-t(rate)
rate[ii]<-1:k
}
} else {
rate<-model
k<-max(rate)
}
Q <- matrix(0, m, m)
index.matrix <- rate
tmp <- cbind(1:m, 1:m)
rate[tmp] <- 0
rate[rate == 0] <- k + 1
Q[]<-c(theta,0)[rate]
diag(Q)<--rowSums(Q)
colnames(Q)<-rownames(Q)<-ss
fit.onerate<-mapply(fitMk,tree=trees,x=x,MoreArgs=list(fixedQ=Q),
SIMPLIFY=FALSE)
-sum(sapply(fit.onerate,logLik))
}
pp<-sapply(fit.multi,function(x) x$rates)
pp<-if(is.matrix(pp)) rowMeans(pp) else mean(pp)
if(length(pp)>1){
fit.onerate<-optim(pp,lik.onerate,trees=trees,x=x,model=model)
} else {
fit.onerate<-optimize(lik.onerate,c(0,1000*pp),trees=trees,x=x,
model=model)
names(fit.onerate)<-c("par","value")
}
rates.multi<-t(sapply(fit.multi,function(x) x$rates))
if(is.character(model)){
m<-length(fit.multi[[1]]$states)
if(model=="ER"){
rates.multi<-t(rates.multi)
colnames(rates.multi)<-"q"
} else if(model=="SYM"){
if(length(fit.multi[[1]]$states)==2) rates.multi<-t(rates.multi)
colnames(rates.multi)<-
sapply(fit.multi[[1]]$states,function(x,y)
sapply(y,function(y,x) paste(x,"<->",y,sep=""),x=x),
y=fit.multi[[1]]$states)[lower.tri(matrix(0,m,m))]
} else if(model=="ARD"){
ii<-(upper.tri(matrix(0,m,m))+lower.tri(matrix(0,m,m))==TRUE)
colnames(rates.multi)<-
sapply(fit.multi[[1]]$states,function(x,y)
sapply(y,function(y,x) paste(y,"->",x,sep=""),x=x),
y=fit.multi[[1]]$states)[ii]
}
} else {
k<-max(fit.multi[[1]]$index.matrix,na.rm=TRUE)
if(k==1) rates.multi<-t(rates.multi)
foo<-function(i,index.matrix,states){
rc<-which(index.matrix==i,arr.ind=TRUE)
lab<-apply(rc,1,function(ind,ss) paste(ss[ind],collapse="->"),
ss=states)
if(length(lab)>1) paste(lab,collapse=",") else lab
}
colnames(rates.multi)<-sapply(1:k,
foo,fit.multi[[1]]$index.matrix,fit.multi[[1]]$states)
}
if(!is.null(names(trees))) rownames(rates.multi)<-names(trees)
else rownames(rates.multi)<-paste("tree",1:length(trees),sep="")
LR<-2*(logL.multi+fit.onerate$value)
km<-sum(sapply(fit.multi,function(x) max(x$index.matrix,na.rm=TRUE)))
k1<-length(pp)
P.chisq<-pchisq(LR,df=km-k1,lower.tail=FALSE)
obj<-list(
multi.rate.model=list(
logL=logL.multi,
rates=rates.multi,
method=fit.multi[[1]]$method),
common.rate.model=list(
logL=-fit.onerate$value,
rates=setNames(fit.onerate$par,colnames(rates.multi)),
method=if(length(pp)>1) "optim" else "optimize"),
index.matrix=fit.multi[[1]]$index.matrix,
states=fit.multi[[1]]$states,
pi=fit.multi[[1]]$pi,
model=model,N=N,n=sapply(trees,Ntip),
likelihood.ratio=LR,P.chisq=P.chisq,
type="discrete")
class(obj)<-"ratebytree"
obj
}
## continuous character ratebytree
rbt.cont<-function(trees,x,...){
if(hasArg(tol)) tol<-list(...)$tol
else tol<-1e-8
if(hasArg(trace)) trace<-list(...)$trace
else trace<-FALSE
if(hasArg(digits)) digits<-list(...)$digits
else digits<-4
if(hasArg(test)) test<-list(...)$test
else test<-"chisq"
if(hasArg(quiet)) quiet<-list(...)$quiet
else quiet<-FALSE
if(hasArg(maxit)) maxit<-list(...)$maxit
else maxit<-500
if(hasArg(regimes)){
regimes<-list(...)$regimes
if(!is.factor(regimes)) regimes<-as.factor(regimes)
} else regimes<-as.factor(1:length(trees))
if(hasArg(model)) model<-list(...)$model
else model<-"BM"
if(!(model%in%c("BM","OU","EB"))){
cat(paste("model =",model,"not recognized. using model = \"BM\"\n"))
model<-"BM"
}
## 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.")
if(hasArg(se)) se<-list(...)$se
else {
se<-x
for(i in 1:length(x)) se[[i]][1:length(se[[i]])]<-0
}
N<-length(trees)
m<-length(levels(regimes))
## reorder the trait vectors in x & SEs in se
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)
## likelihood functions
lik.multi<-function(theta,trees,y,se,model,regimes,trace=FALSE){
m<-length(unique(regimes))
ind<-setNames(lapply(levels(regimes),function(x,y) which(y==x),y=regimes),
levels(regimes))
THETA<-list()
for(i in 1:m){
THETA[[i]]<-c(theta[i],theta[ind[[i]]+m])
if(model%in%c("OU","EB")) THETA[[i]]<-c(THETA[[i]],theta[i+m+N])
}
logL<-0
for(i in 1:m)
logL<-logL-lik.onerate(THETA[[i]],trees[ind[[i]]],y[ind[[i]]],se[ind[[i]]],
model=model,trace=FALSE)
if(trace){
cat(paste(paste(round(theta[1:m],digits),collapse="\t"),
if(model=="OU") paste(round(theta[1:m+m+N],digits),collapse="\t"),
if(model=="EB") paste(round(theta[1:m+m+N],digits),collapse="\t"),
round(logL,digits),"\n",sep="\t"))
flush.console()
}
-logL
}
lik.onerate<-function(theta,trees,y,se,model,trace=FALSE){
n<-sapply(trees,Ntip)
N<-length(trees)
sig<-theta[1]
a<-theta[1:N+1]
if(model=="OU") alpha<-theta[N+2]
if(model=="EB") r<-theta[N+2]
if(model=="BM") C<-lapply(trees,vcv)
else if(model=="OU") C<-lapply(trees,vcvPhylo,model="OU",alpha=alpha,
anc.nodes=FALSE)
else if(model=="EB") C<-lapply(trees,vcvPhylo,model="EB",r=r,
anc.nodes=FALSE)
E<-lapply(se,diag)
V<-mapply("+",lapply(C,"*",sig),E,SIMPLIFY=FALSE)
logL<-0
for(i in 1:N)
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),
if(model=="OU") round(alpha,digits),
if(model=="EB") round(r,digits),
round(logL,digits),"\n",sep="\t"))
flush.console()
}
-logL
}
## first, fit multi-rate model
f1<-function(tree,x){
pvcv<-phyl.vcv(as.matrix(x),vcv(tree),1)
c(pvcv$R[1,1],pvcv$a[1,1])
}
P<-mapply(f1,trees,x,SIMPLIFY=FALSE)
PP<-list()
PP[[1]]<-vector()
for(i in 1:m)
PP[[1]][i]<-mean(sapply(P,function(x) x[1])[which(regimes==levels(regimes)[i])])
PP[[2]]<-sapply(P,function(x) x[2])
if(model%in%c("OU","EB")) PP[[3]]<-rep(0,m)
if(hasArg(init)){
init<-list(...)$init
if(!is.null(init$sigm)) PP[[1]]<-init$sigm
if(!is.null(init$am)) PP[[2]]<-init$am
if(model=="OU") if(!is.null(init$alpham)) PP[[3]]<-init$alpham
if(model=="EB") if(!is.null(init$rm)) PP[[3]]<-init$rm
}
p<-unlist(PP)
if(trace){
if(model=="BM"){
cat("\nOptimizing multi-rate model....\n")
cat(paste(paste("sig[",1:m,"]",sep="",collapse="\t"),"logL\n",sep="\t"))
} else if(model=="OU"){
cat("\nOptimizing multi-regime model....\n")
cat(paste(paste("sig[",1:m,"]",sep="",collapse="\t"),
paste("alpha[",1:m,"]",sep="",collapse="\t"),"logL\n",sep="\t"))
} else if(model=="EB"){
cat("\nOptimizing multi-regime model....\n")
cat(paste(paste("sig[",1:m,"]",sep="",collapse="\t"),
paste("r[",1:m,"]",sep="",collapse="\t"),"logL\n",sep="\t"))
}
}
fit.multi<-optim(p,lik.multi,trees=trees,y=x,se=se,model=model,
regimes=regimes,trace=trace,
method="L-BFGS-B",lower=c(rep(tol,m),rep(-Inf,N),
if(model%in%c("OU","EB")) rep(-Inf,m)),upper=c(rep(Inf,m),rep(Inf,N),
if(model%in%c("OU","EB")) rep(Inf,m)),control=list(maxit=maxit))
## compute covariance matrix
H.multi<-hessian(lik.multi,fit.multi$par,trees=trees,y=x,
se=se,model=model,regimes=regimes,trace=FALSE)
Cov.multi<-if(qr(H.multi)$rank!=ncol(H.multi)) ginv(H.multi) else solve(H.multi)
## now fit single-rate model
p<-c(mean(fit.multi$par[1:m]),fit.multi$par[1:N+m])
if(model%in%c("OU","EB")) p<-c(p,mean(fit.multi$par[1:m+m+N]))
if(hasArg(init)){
if(!is.null(init$sigc)) p[1]<-init$sigc
if(!is.null(init$ac)) p[1:N+1]<-init$ac
if(model=="OU") if(!is.null(init$alphac)) p[N+2]<-init$alphac
if(model=="EB") if(!is.null(init$rc)) p[N+2]<-init$rc
}
if(trace){
if(model=="BM"){
cat("\nOptimizing common-rate model....\n")
cat(paste("sig ","logL\n",sep="\t"))
} else if(model=="OU"){
cat("\nOptimizing common-regime model....\n")
cat(paste("sig ","alpha ","logL\n",sep="\t"))
} else if(model=="EB"){
cat("\nOptimizing common-regime mode.....\n")
cat(paste("sig ","r ","logL\n",sep="\t"))
}
}
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),
control=list(maxit=maxit))
## compute covariance matrix
H.onerate<-hessian(lik.onerate,fit.onerate$par,trees=trees,y=x,
se=se,model=model,trace=FALSE)
Cov.onerate<-if(qr(H.onerate)$rank!=ncol(H.onerate)) ginv(H.onerate) else solve(H.onerate)
## compare models:
LR<-2*(-fit.multi$value+fit.onerate$value)
km<-N+m+if(model=="BM") 0 else if(model%in%c("OU","EB")) m
k1<-N+if(model=="BM") 1 else if(model%in%c("OU","EB")) 2
if(test=="simulation"&&model%in%c("OU","EB")){
cat("Simulation test not yet available for chosen model. Using chi-square test.\n")
test<-"chisq"
}
if(test=="chisq") P.chisq<-pchisq(LR,df=km-k1,lower.tail=FALSE)
else if(test=="simulation"){
if(!quiet) cat("Generating null distribution via simulation -> |")
flush.console()
if(hasArg(nsim)) nsim<-list(...)$nsim
else nsim<-100
X<-mapply(fastBM,tree=trees,a=as.list(fit.onerate$par[1:N+1]),
MoreArgs=list(sig2=fit.onerate$par[1],nsim=nsim),
SIMPLIFY=FALSE)
P.sim<-1/(nsim+1)
pct<-0.1
for(i in 1:nsim){
x.sim<-lapply(X,function(x,ind) x[,ind],ind=i)
f2<-function(x,se) sampleFrom(xbar=x,xvar=se^2,n=rep(1,length(x)))
x.sim<-mapply(f2,x=x.sim,se=se,SIMPLIFY=FALSE)
fit.sim<-ratebytree(trees,x.sim,se=se)
P.sim<-P.sim+(fit.sim$likelihood.ratio>=LR)/(nsim+1)
if(i/nsim>=pct){
if(!quiet) cat(".")
flush.console()
pct<-pct+0.1
}
}
if(!quiet) cat(".|\nDone!\n")
flush.console()
}
obj<-list(
multi.rate.model=list(sig2=fit.multi$par[1:m],
SE.sig2=sqrt(diag(Cov.multi)[1:m]),
a=fit.multi$par[1:N+m],
SE.a=sqrt(diag(Cov.multi)[1:N+m]),
alpha=if(model=="OU") fit.multi$par[1:m+m+N] else NULL,
SE.alpha=if(model=="OU") sqrt(diag(Cov.multi)[1:m+m+N]) else NULL,
r=if(model=="EB") fit.multi$par[1:N+m+N] else NULL,
SE.r=if(model=="EB") sqrt(diag(Cov.multi)[1:N+m+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,n=sapply(trees,Ntip),likelihood.ratio=LR,
regimes=regimes,m=m,
P.chisq=if(test=="chisq") P.chisq else NULL,
P.sim=if(test=="simulation") P.sim else NULL,
type="continuous")
class(obj)<-"ratebytree"
obj
}
## S3 print method for ratebytree
print.ratebytree<-function(x,...){
if(hasArg(digits)) digits<-list(...)$digits
else digits<-4
if(x$type=="continuous") prbt.cont(x,digits=digits)
else if(x$type=="discrete") prbt.disc(x,digits=digits)
else if(x$type=="diversification") prbt.div(x,digits=digits)
else print(x)
}
prbt.div<-function(x,digits=digits){
cat("ML common diversification-rate model:")
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),
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"),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"))
}
prbt.cont<-function(x,digits=digits){
N<-x$N
m<-x$m
if(x$model=="BM"){
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",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[",levels(x$regimes),"]",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",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")
cat(paste("\ts^2\t",paste(paste("a[",1:N,"]",sep=""),collapse="\t")),
"\talpha\tk\tlogL\n")
cat(paste("value",round(x$common.rate.model$sig2,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[",levels(x$regimes),"]",sep=""),collapse="\t"),"\t",
paste(paste("a[",1:N,"]",sep=""),collapse="\t"),"\t",
paste(paste("alp[",levels(x$regimes),"]",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"),
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")
cat(paste("\ts^2\t",paste(paste("a[",1:N,"]",sep=""),collapse="\t")),
"\tr\tk\tlogL\n")
cat(paste("value",round(x$common.rate.model$sig2,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[",levels(x$regimes),"]",sep=""),collapse="\t"),"\t",
paste(paste("a[",1:N,"]",sep=""),collapse="\t"),"\t",
paste(paste("r[",levels(x$regimes),"]",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"),
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"))
if(!is.null(x$P.chisq))
cat(paste("P-value (based on X^2):",round(x$P.chisq,digits),"\n\n"))
else if(!is.null(x$P.sim))
cat(paste("P-value (based on simulation):",round(x$P.sim,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")
}
prbt.disc<-function(x,digits=digits){
cat("ML common-rate model:\n")
cat(paste("\t",paste(names(x$common.rate.model$rates),collapse="\t"),
"\tk\tlogL\n",sep=""))
cat(paste("value\t",paste(round(x$common.rate.model$rates,digits),
collapse="\t"),"\t",max(x$index.matrix,na.rm=T),"\t",
round(x$common.rate.model$logL,digits),"\n",sep=""))
cat(paste("\nModel fitting method was \"",x$common.rate.model$method,
"\".\n",sep=""))
cat("\nML multi-rate model:\n")
cat(paste("\t",paste(colnames(x$multi.rate.model$rates),collapse="\t"),
"\tk\tlogL\n",sep=""))
for(i in 1:nrow(x$multi.rate.model$rates)){
if(i>1) cat(paste(rownames(x$multi.rate.model$rates)[i],"\t",
paste(round(x$multi.rate.model$rates[i,],
digits),collapse="\t"),"\n",sep=""))
else if(i==1) cat(paste(rownames(x$multi.rate.model$rates)[i],"\t",
paste(round(x$multi.rate.model$rates[i,],
digits),collapse="\t"),"\t",x$N*max(x$index.matrix,na.rm=T),"\t",
round(x$multi.rate.model$logL,digits),"\n",sep=""))
}
cat(paste("\nModel 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"))
}
## 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$m,x$m)
for(i in 1:x$m){
for(j in 1:x$m){
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("reg.",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("reg.",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=""))
}
AIC.ratebytree<-function(object,...,k=2){
aic<-data.frame(AIC=c(k*object$common.rate.model$k-2*object$common.rate.model$logL,
k*object$multi.rate.model$k-2*object$multi.rate.model$logL),
df=c(object$common.rate.model$k,object$multi.rate.model$k))
model.names<-if(is.null(object$model)) 1 else object$model
addtl.obj<-list(...)
if(length(addtl.obj)>0){
for(i in 1:length(addtl.obj)){
aic<-rbind(aic,c(k*addtl.obj[[i]]$multi.rate.model$k-
2*addtl.obj[[i]]$multi.rate.model$logL,
addtl.obj[[i]]$multi.rate.model$k))
model.names<-c(model.names,
if(is.null(addtl.obj[[i]]$model)) 1+i else addtl.obj[[i]]$model)
}
rownames(aic)<-c("common-rate",paste("multi-rate:",model.names,sep=""))
} else rownames(aic)<-c("common-rate","multi-rate")
aic
}
Try the phytools package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.