From 48b6b6fe4e38ae122c339784213732935e9b5343 Mon Sep 17 00:00:00 2001 From: paradis Date: Fri, 20 Feb 2009 09:45:11 +0000 Subject: [PATCH] new yule.time + several corrections in man pages git-svn-id: https://svn.mpl.ird.fr/ape/dev/ape@62 6e262413-ae40-0410-9e79-b911bd7a66b7 --- ChangeLog | 10 ++++++ DESCRIPTION | 4 +-- R/yule.time.R | 71 +++++++++++++++++++++++++++++++++++++++++++ man/corGrafen.Rd | 2 +- man/corPagel.Rd | 6 ++-- man/read.tree.Rd | 4 +-- man/weight.taxo.Rd | 2 +- man/yule.cov.Rd | 5 ++- man/yule.time.Rd | 76 ++++++++++++++++++++++++++++++++++++++++++++++ 9 files changed, 168 insertions(+), 12 deletions(-) create mode 100644 R/yule.time.R create mode 100644 man/yule.time.Rd diff --git a/ChangeLog b/ChangeLog index 2f05200..d3e56a5 100644 --- a/ChangeLog +++ b/ChangeLog @@ -1,3 +1,13 @@ + 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 diff --git a/DESCRIPTION b/DESCRIPTION index 538343c..d168fa6 100644 --- a/DESCRIPTION +++ b/DESCRIPTION @@ -1,6 +1,6 @@ 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, diff --git a/R/yule.time.R b/R/yule.time.R new file mode 100644 index 0000000..f82aab8 --- /dev/null +++ b/R/yule.time.R @@ -0,0 +1,71 @@ +## 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") +} diff --git a/man/corGrafen.Rd b/man/corGrafen.Rd index ec8375f..9f063da 100644 --- a/man/corGrafen.Rd +++ b/man/corGrafen.Rd @@ -10,7 +10,7 @@ corGrafen(value, phy, form=~1, fixed = FALSE) 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}} diff --git a/man/corPagel.Rd b/man/corPagel.Rd index 4509b44..096a248 100644 --- a/man/corPagel.Rd +++ b/man/corPagel.Rd @@ -11,11 +11,11 @@ corPagel(value, phy, form = ~1, fixed = FALSE) } \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 @@ -29,7 +29,7 @@ corPagel(value, phy, form = ~1, fixed = FALSE) \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{ diff --git a/man/read.tree.Rd b/man/read.tree.Rd index 75f4cb5..9428c1e 100644 --- a/man/read.tree.Rd +++ b/man/read.tree.Rd @@ -71,8 +71,8 @@ read.tree(file = "", text = NULL, tree.names = NULL, 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}} diff --git a/man/weight.taxo.Rd b/man/weight.taxo.Rd index 66b3d1e..b6bc15c 100644 --- a/man/weight.taxo.Rd +++ b/man/weight.taxo.Rd @@ -18,7 +18,7 @@ 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. diff --git a/man/yule.cov.Rd b/man/yule.cov.Rd index dbdc044..b749aa5 100644 --- a/man/yule.cov.Rd +++ b/man/yule.cov.Rd @@ -65,8 +65,7 @@ The user must obtain the values for the nodes separately. 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. @@ -92,9 +91,9 @@ yule.cov(bird.orders, ~ x) ### 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} diff --git a/man/yule.time.Rd b/man/yule.time.Rd new file mode 100644 index 0000000..a5c2689 --- /dev/null +++ b/man/yule.time.Rd @@ -0,0 +1,76 @@ +\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} -- 2.39.5