Permalink
Browse files

new function fitmultiMk

  • Loading branch information...
liamrevell committed Dec 5, 2017
1 parent 4bcfbda commit 8bd2de3fa8c129ddcf4bc7b88d7ddf683653a6d1
Showing with 180 additions and 9 deletions.
  1. +4 −4 DESCRIPTION
  2. +6 −2 NAMESPACE
  3. +160 −0 R/fitmultiMk.R
  4. +10 −3 man/fitMk.Rd
View
@@ -1,6 +1,6 @@
Package: phytools
Version: 0.6-49
Date: 2017-12-4
Version: 0.6-50
Date: 2017-12-5
Title: Phylogenetic Tools for Comparative Biology (and Other Things)
Author: Liam J. Revell
Maintainer: Liam J. Revell <liam.revell@umb.edu>
@@ -13,6 +13,6 @@ ZipData: no
Description: A wide range of functions for phylogenetic analysis. Functionality is concentrated in phylogenetic comparative biology, but also includes a diverse array of methods for visualizing, manipulating, reading or writing, and even inferring phylogenetic trees and data. Included among the functions in phylogenetic comparative biology are various for ancestral state reconstruction, model-fitting, simulation of phylogenies and data, and multivariate analysis. There are a broad range of plotting methods for phylogenies and comparative data which include, but are not restricted to, methods for mapping trait evolution on trees, for projecting trees into phenotypic space or a geographic map, and for visualizing correlated speciation between trees. Finally, there are a number of functions for reading, writing, analyzing, inferring, simulating, and manipulating phylogenetic trees and comparative data not covered by other packages. For instance, there are functions for randomly or non-randomly attaching species or clades to a phylogeny, for estimating supertrees or consensus phylogenies from a set, for simulating trees and phylogenetic data under a range of models, and for a wide variety of other manipulations and analyses that phylogenetic biologists might find useful in their research.
License: GPL (>= 2)
URL: http://github.com/liamrevell/phytools
Packaged: 2017-12-4 12:00:00 EST
Packaged: 2017-12-5 12:00:00 EST
Repository:
Date/Publication: 2017-12-4 12:00:00 EST
Date/Publication: 2017-12-5 12:00:00 EST
View
@@ -7,8 +7,8 @@ export(density.multiSimmap, densityMap, densityTree, describe.simmap, di2multi.s
export(drop.tip.contMap, drop.tip.densityMap, drop.tip.simmap, drop.tip.singleton)
export(edgelabels.cophylo, edgeProbs, errorbar.contMap, estDiversity, evol.rate.mcmc, evol.vcv, evolvcv.lite, exhaustiveMP)
export(expand.clade, export.as.xml, extract.clade.simmap, extract.strahlerNumber)
export(fancyTree, fastAnc, fastBM, fastDist, fastHeight, fastMRCA, findMRCA, fit.bd, fit.yule, fitBayes, fitMk, fitDiversityModel, fitPagel)
export(force.ultrametric)
export(fancyTree, fastAnc, fastBM, fastDist, fastHeight, fastMRCA, findMRCA, fit.bd, fit.yule, fitBayes, fitMk, fitmultiMk, fitDiversityModel)
export(fitPagel, force.ultrametric)
export(gammatest, genus.to.species.tree, genSeq, geo.legend, get.treepos, getCladesofSize, getDescendants, getExtant, getExtinct, getnode)
export(getParent, getSisters, getStates)
export(labelnodes, ladderize.simmap, lambda.transform, lik.bd, likMlambda, likSurface.rateshift, linklabels, locate.fossil, locate.yeti)
@@ -138,6 +138,10 @@ S3method(print, ctt)
S3method(print, multiCtt)
S3method(plot, ctt)
S3method(plot, multiCtt)
S3method(print, fitmultiMk)
S3method(summary, fitmultiMk)
S3method(logLik, fitmultiMk)
S3method(AIC, fitmultiMk)
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,160 @@
## new function to fit multi-regime Mk model
## written by Liam J. Revell 2017
## adapted from ape::ace by E. Paradis et al. 2013
fitmultiMk<-function(tree,x,model="ER",...){
if(!inherits(tree,"simmap")){
stop("tree should be an object of class \"simmap\". Use fitMk.\n")
} else {
regimes<-mapped.states(tree)
nregimes<-length(regimes)
}
if(hasArg(q.init)) q.init<-list(...)$q.init
else q.init<-length(unique(x))/sum(tree$edge.length)
if(hasArg(opt.method)) opt.method<-list(...)$opt.method
else opt.method<-"nlminb"
if(hasArg(min.q)) min.q<-list(...)$min.q
else min.q<-1e-12
if(is.matrix(x)){
x<-x[tree$tip.label,]
m<-ncol(x)
states<-colnames(x)
} else {
x<-to.matrix(x,sort(unique(x)))
x<-x[tree$tip.label,]
m<-ncol(x)
states<-colnames(x)
}
if(hasArg(pi)) pi<-list(...)$pi
else pi<-"equal"
if(pi[1]=="equal") pi<-setNames(rep(1/m,m),states)
else pi<-pi/sum(pi)
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 {
if(ncol(model)!=nrow(model))
stop("model is not a square matrix")
if(ncol(model)!=ncol(x))
stop("model does not have the right number of columns")
rate<-model
k<-max(rate)
}
Q<-replicate(nregimes,matrix(0,m,m),simplify=FALSE)
names(Q)<-regimes
index.matrix<-rate
tmp<-cbind(1:m,1:m)
rate[tmp]<-0
rate[rate==0]<-k+1
pw<-reorder(map.to.singleton(tree),"pruningwise")
N<-Ntip(pw)
M<-pw$Nnode
liks<-rbind(x,matrix(0,M,m,dimnames=list(1:M+N,states)))
lik<-function(pp,output.liks=FALSE,pi){
if(any(is.nan(pp))||any(is.infinite(pp))) return(1e50)
comp<-vector(length=N+M,mode="numeric")
for(i in 1:nregimes){
Q[[i]][]<-c(pp[1:k+(i-1)*k],0)[rate]
diag(Q[[i]])<--rowSums(Q[[i]])
}
parents<-unique(pw$edge[,1])
root<-min(parents)
for(i in 1:length(parents)){
anc<-parents[i]
ii<-which(pw$edge[,1]==parents[i])
desc<-pw$edge[ii,2]
el<-pw$edge.length[ii]
v<-vector(length=length(desc),mode="list")
reg<-names(pw$edge.length)[ii]
for(j in 1:length(v))
v[[j]]<-EXPM(Q[[reg[j]]]*el[j])%*%liks[desc[j],]
vv<-if(anc==root) Reduce('*',v)[,1]*pi else
Reduce('*',v)[,1]
comp[anc]<-sum(vv)
liks[anc,]<-vv/comp[anc]
}
logL<--sum(log(comp[1:M+N]))
return(if(is.na(logL)) Inf else logL)
}
if(length(q.init)!=(k*nregimes)) q.init<-rep(q.init[1],k*nregimes)
if(opt.method=="optim")
fit<-optim(q.init,function(p) lik(p,pi=pi),method="L-BFGS-B",
lower=rep(min.q,k))
else
fit<-nlminb(q.init,function(p) lik(p,pi=pi),
lower=rep(0,k*nregimes),upper=rep(1e50,k*nregimes))
obj<-list(logLik=
if(opt.method=="optim") -fit$value else -fit$objective,
rates=fit$par,
index.matrix=index.matrix,
states=states,
regimes=regimes,
pi=pi,
method=opt.method)
class(obj)<-"fitmultiMk"
return(obj)
}
## print method for objects of class "fitMk"
print.fitmultiMk<-function(x,digits=6,...){
k<-max(x$index.matrix,na.rm=TRUE)
cat("Object of class \"fitmultiMk\".\n\n")
for(i in 1:length(x$regimes)){
cat(paste("Fitted value of Q[",x$regimes[i],"]:\n",sep=""))
Q<-matrix(NA,length(x$states),length(x$states))
Q[]<-c(0,x$rates[1:k+(i-1)*k])[x$index.matrix+1]
diag(Q)<-0
diag(Q)<--rowSums(Q)
colnames(Q)<-rownames(Q)<-x$states
print(round(Q,digits))
cat("\n")
}
cat("Fitted (or set) value of pi:\n")
print(x$pi)
cat(paste("\nLog-likelihood:",round(x$logLik,digits),"\n"))
cat(paste("\nOptimization method used was \"",x$method,"\"\n\n",sep=""))
}
## summary method for objects of class "fitmultiMk"
summary.fitmultiMk<-function(object,...){
if(hasArg(digits)) digits<-list(...)$digits
else digits<-6
if(hasArg(quiet)) quiet<-list(...)$quiet
else quiet<-FALSE
Q<-list()
for(i in 1:length(object$regimes)){
if(!quiet) cat(paste("Fitted value of Q[",object$regimes[i],
"]:\n",sep=""))
Q[[i]]<-matrix(NA,length(object$states),length(object$states))
Q[[i]][]<-c(0,object$rates)[object$index.matrix+1]
diag(Q[[i]])<-0
diag(Q[[i]])<--rowSums(Q[[i]])
colnames(Q[[i]])<-rownames(Q[[i]])<-object$states
if(!quiet) print(round(Q[[i]],digits))
cat("\n")
}
names(Q)<-object$regimes
if(!quiet) cat(paste("Log-likelihood:",round(object$logLik,digits),"\n\n"))
invisible(list(Q=Q,logLik=object$logLik))
}
## logLik method for objects of class "fitmultiMk"
logLik.fitmultiMk<-function(object,...) object$logLik
## AIC method
AIC.fitmultiMk<-function(object,...,k=2){
np<-length(object$rates)
-2*logLik(object)+np*k
}
View
@@ -2,26 +2,33 @@
\alias{fitMk}
\alias{plot.fitMk}
\alias{plot.gfit}
\alias{fitmultiMk}
\title{Fits Mk model}
\usage{
fitMk(tree, x, model="SYM", fixedQ=NULL, ...)
\method{plot}{fitMk}(x, ...)
\method{plot}{gfit}(x, ...)
fitmultiMk(tree, x, model="ER", ...)
}
\arguments{
\item{tree}{an object of class \code{"phylo"}.}
\item{tree}{an object of class \code{"phylo"}. In the case of \code{fitmultiMk} an object of class \code{"simmap"} with a mapped discrete character.}
\item{x}{a vector of tip values for species; \code{names(x)} should be the species names. In the case of \code{plot.fitMk}, an object of class \code{"fitMk"}.}
\item{model}{model. See \code{make.simmap} or \code{ace} for details.}
\item{fixedQ}{fixed value of transition matrix \code{Q}, if one is desired.}
\item{...}{optional arguments, including \code{pi}, the prior distribution at the root node (defaults to \code{pi="equal"}). For \code{plot} method optional arguments include (but may not be limited to): \code{signif}, the number of digits for the rates to be plotted; \code{main}, a character vector of length two with the headings for each subplot; \code{cex.main}, \code{cex.traits}, and \code{cex.rates}, font sizes for the various text elements of the plot; and \code{show.zeros}, a logical argument specifying whether or not to plot arrows with the ML estimated transition rate is not different from zero (with tolerance specified by the optional argument \code{tol}).}
}
\description{
This function fits a so-called M\emph{k} model for discrete character evolution. It is primarily designed to be used inside of \code{make.simmap}.
The function \code{fitMk} fits a so-called M\emph{k} model for discrete character evolution. It is primarily designed to be used inside of \code{make.simmap}.
Two \code{plot} methods are available. \code{plot.fitMk} plots an object of class \code{"fitMk"} returned by \code{fitMk}. \code{plot.gfit} plots an object of class \code{"gfit"} from geiger's \code{fitDiscrete} function. Both plots portray the fitted model using a graph of arrows connecting states.
The newer function \code{fitmultiMk} fits an M\emph{k} model in which the transition rates between character states are allowed to vary depending on the mapped state of a discrete character on the tree. It can be combined with, for example, \code{\link{paintSubTree}} to test hypotheses about how the process of discrete character evolution for \code{x} varies between different parts of the tree.
}
\details{
Note that both \code{fitMk} and \code{fitmultiMk} recycle code from \code{\link{ace}} in the \emph{ape} package for computing the likelihood.
}
\value{
An object of class \code{"fitMk"}. In the case of \code{plot.fitMk}, a plotted M\emph{k} model.
An object of class \code{"fitMk"} or \code{"fitmultiMk"}. In the case of \code{plot.fitMk}, a plotted M\emph{k} model.
}
\references{
Revell, L. J. (2012) phytools: An R package for phylogenetic comparative biology (and other things). \emph{Methods Ecol. Evol.}, \bold{3}, 217-223.

0 comments on commit 8bd2de3

Please sign in to comment.