+ CHANGES IN APE VERSION 2.2-5
+
+
+NEW FEATURES
+
+ o The new function yule.time fits a user-defined time-dependent
+ Yule model by maximum likelihood.
+
+
+
CHANGES IN APE VERSION 2.2-4
BUG FIXES
Package: ape
-Version: 2.2-4
-Date: 2009-01-19
+Version: 2.2-5
+Date: 2009-02-20
Title: Analyses of Phylogenetics and Evolution
Author: Emmanuel Paradis, Ben Bolker, Julien Claude, Hoa Sien Cuong,
Richard Desper, Benoit Durand, Julien Dutheil, Olivier Gascuel,
--- /dev/null
+## yule.time.R (2009-02-20)
+
+## Fits the Time-Dependent Yule Model
+
+## Copyright 2009 Emmanuel Paradis
+
+## This file is part of the R-package `ape'.
+## See the file ../COPYING for licensing issues.
+
+yule.time <- function(phy, birth, BIRTH = NULL, root.time = 0,
+ opti = "nlm", start = 0.01)
+{
+ opti <- pmatch(opti, c("nlm", "nlminb", "optim"))
+ if (is.na(opti)) stop("ambiguous argument 'opti'")
+ LAMBDA <- function() x
+ body(LAMBDA) <- body(birth)
+ formals(LAMBDA) <- alist(t=)
+ BT <- branching.times(phy)
+ T <- BT[1]
+ x <- BT[1] - BT + root.time
+ m <- phy$Nnode
+
+ paranam <- c(names(formals(birth)))
+ np <- length(paranam)
+ start <- rep(start, length.out = np)
+
+ ## Foo is always vectorized
+ if (is.null(BIRTH)) {
+ Foo <- function(x) {
+ n <- length(x)
+ res <- numeric(n)
+ for (i in 1:n)
+ res[i] <- integrate(LAMBDA, x[i], T)$value
+ res
+ }
+ } else {
+ environment(BIRTH) <- environment()
+ Foo <- function(x) BIRTH(T) - BIRTH(x)
+ }
+
+ half.dev <- function(p) {
+ for (i in 1:np)
+ assign(paranam[i], p[i], pos = sys.frame(1))
+ root.term <-
+ if (is.null(BIRTH)) integrate(LAMBDA, x[1], T)$value
+ else BIRTH(T) - BIRTH(x[1])
+ sum(Foo(x)) + root.term - sum(log(LAMBDA(x[2:m])))
+ }
+
+ switch(opti,
+ {
+ out <- nlm(half.dev, start, hessian = TRUE)
+ est <- out$estimate
+ se <- sqrt(diag(solve(out$hessian)))
+ loglik <- lfactorial(m) - out$minimum
+ },{
+ out <- nlminb(start, half.dev)
+ est <- out$par
+ se <- NULL
+ loglik <- lfactorial(m) - out$objective
+ },{
+ out <- optim(start, half.dev, hessian = TRUE,
+ control = list(maxit = 1000), method = "BFGS")
+ est <- out$par
+ se <- sqrt(diag(solve(out$hessian)))
+ loglik <- lfactorial(m) - out$value
+ })
+ names(est) <- paranam
+ if (!is.null(se)) names(se) <- paranam
+ structure(list(estimate = est, se = se, loglik = loglik), class = "yule")
+}
covariate = getCovariate(object), corr = TRUE, ...)
}
\arguments{
- \item{value}{The \eqn{\alpha}{rho} parameter}
+ \item{value}{The \eqn{\rho}{rho} parameter}
\item{phy}{An object of class \code{phylo} representing the phylogeny
(branch lengths are ignored) to consider}
\item{object}{An (initialized) object of class \code{corGrafen}}
}
\arguments{
\item{value}{the (initial) value of the parameter
- \eqn{\gamma}{gamma}.}
+ \eqn{\lambda}{lambda}.}
\item{phy}{an object of class \code{"phylo"}.}
\item{form}{(ignored).}
\item{fixed}{a logical specifying whether \code{gls} should
- estimate \eqn{\gamma}{gamma} (the default) or keep it fixed.}
+ estimate \eqn{\lambda}{lambda} (the default) or keep it fixed.}
\item{object}{an (initialized) object of class \code{"corPagel"}.}
\item{covariate}{(ignored).}
\item{corr}{a logical value specifying whether to return the
\description{
The correlation structure from the present model is derived from the
Brownian motion model by multiplying the off-diagonal elements (i.e.,
- the covariances) by \eqn{\gamma}{gamma}. The variances are thus the
+ the covariances) by \eqn{\lambda}{lambda}. The variances are thus the
same than for a Brownian motion model.
}
\value{
Olsen, G. Interpretation of the "Newick's 8:45" tree format standard.
\url{http://evolution.genetics.washington.edu/phylip/newick_doc.html}
- Paradis, E. (2006) Definition of Formats for Coding Phylogenetic Trees
- in R. \url{http://ape.mpl.ird.fr/misc/FormatTreeR_4Dec2006.pdf}
+ Paradis, E. (2008) Definition of Formats for Coding Phylogenetic Trees
+ in R. \url{http://ape.mpl.ird.fr/misc/FormatTreeR_28July2008.pdf}
}
\author{Emmanuel Paradis \email{Emmanuel.Paradis@mpl.ird.fr}}
The diagonal [i, i] is always set to 0.
- This returned matrix can be used as a weight matrix in
+ The returned matrix can be used as a weight matrix in
\code{\link{Moran.I}}. \code{x} and \code{y} may be vectors of
factors.
change in a species trait is more or less continuous between two nodes
or between a node and a tip. Thus reconstructing the ancestral values
with a Brownian motion model may be consistent with the present
- method. This can be done with the function \code{\link{pic}} but
- currently needs some hacking!
+ method. This can be done with the function \code{\link{ace}}.
}
\value{
A NULL value is returned, the results are simply printed.
### with a likelihood ratio test
yule(bird.orders)
### another example with a tree that has a multichotomy
-### but we cannot run yule() because of this!
data(bird.families)
y <- rnorm(272) # 137 tips + 135 nodes
yule.cov(bird.families, ~ y)
+yule(multi2di(bird.families))
}
\keyword{models}
--- /dev/null
+\name{yule.time}
+\alias{yule.time}
+\title{Fits the Time-Dependent Yule Model}
+\usage{
+yule.time(phy, birth, BIRTH = NULL, root.time = 0, opti = "nlm", start = 0.01)
+}
+\arguments{
+ \item{phy}{an object of class \code{"phylo"}.}
+ \item{birth}{a (vectorized) function specifying how the birth
+ (speciation) probability changes through time (see details).}
+ \item{BIRTH}{a (vectorized) function giving the primitive of
+ \code{birth}.}
+ \item{root.time}{a numeric value giving the time of the root node
+ (time is measured from the past towards the present).}
+ \item{opti}{a character string giving the function used for
+ optimisation of the likelihood function. Three choices are possible:
+ \code{"nlm"}, \code{"nlminb"}, or \code{"optim"}, or any unambiguous
+ abbreviation of these.}
+ \item{start}{the initial values used in the optimisation.}
+}
+\description{
+ This function fits by maximum likelihood the time-dependent Yule
+ model. The time is measured from the past (\code{root.time}) to the
+ present.
+}
+\details{
+ The model fitted is a straightforward extension of the Yule model with
+ covariates (see \code{\link{yule.cov}}). Rather than having
+ heterogeneity among lineages, the speciation probability is the same
+ for all lineages at a given time, but can change through time.
+
+ The function \code{birth} \emph{must} meet these two requirements: (i)
+ the parameters to be estimated are the formal arguments; (ii) time is
+ named \code{t} in the body of the function. However, this is the
+ opposite for the primitive \code{BIRTH}: \code{t} is the formal
+ argument, and the parameters are used in its body. See the examples.
+
+ It is recommended to use \code{BIRTH} if possible, and required if
+ speciation probability is constant on some time interval. If this
+ primitive cannot be provided, a numerical integration is done with
+ \code{\link[base]{integrate}}.
+
+ The standard-errors of the parameters are computed with the Hessian of
+ the log-likelihood function.
+}
+\value{
+ An object of class "yule" (see \code{\link{yule}}).
+}
+\author{Emmanuel Paradis \email{Emmanuel.Paradis@mpl.ird.fr}}
+\seealso{
+ \code{\link{branching.times}}, \code{\link{ltt.plot}},
+ \code{\link{birthdeath}}, \code{\link{yule}}, \code{\link{yule.cov}}
+}
+\examples{
+### define two models...
+birth.logis <- function(a, b) 1/(1 + exp(-a*t - b)) # logistic
+birth.step <- function(l1, l2, Tcl) { # 2 rates with one break-point
+ ans <- rep(l1, length(t))
+ ans[t > Tcl] <- l2
+ ans
+}
+### ... and their primitives:
+BIRTH.logis <- function(a, b) log(exp(-a*t) + exp(b))/a + t
+BIRTH.step <- function(t)
+{
+ if (t <= Tcl) return(t*l1)
+ Tcl*l1 + (t - Tcl)*l2
+}
+data(bird.families)
+### fit both models:
+yule.time(bird.families, birth.logis)
+yule.time(bird.families, birth.logis, BIRTH.logis) # same but faster
+yule.time(bird.families, birth.step) # fails
+yule.time(bird.families, birth.step, BIRTH.step)
+}
+\keyword{models}
projects / ape.git / commitdiff