-## chronopl.R (2007-01-17)
+## chronopl.R (2008-03-26)
## Molecular Dating With Penalized Likelihood
-## Copyright 2005-2007 Emmanuel Paradis
+## Copyright 2005-2008 Emmanuel Paradis
## This file is part of the R-package `ape'.
## See the file ../COPYING for licensing issues.
-chronopl <- function(phy, lambda, node.age = 1, node = "root",
- CV = FALSE)
+chronopl <- function(phy, lambda, age.min = 1, age.max = NULL,
+ node = "root", S = 1, tol = 1e-8,
+ CV = FALSE, eval.max = 500, iter.max = 500, ...)
{
n <- length(phy$tip.label)
- n.node <- phy$Nnode
- if (n != n.node + 1)
- stop("the tree must be rooted AND dichotomous")
- if (any(phy$edge.length == 0))
- stop("some branch lengths are equal to zero;
-you must remove them beforehand.")
- N <- dim(phy$edge)[1]
ROOT <- n + 1
- if (node == "root") node <- ROOT
- ini.rate <- phy$edge.length
- ## `known.ages' contains the index of all nodes (internal and
- ## terminal) of known age:
- known.ages <- c(1:n, node)
- ## `unknown.ages' contains the index of the nodes of unknown age:
- unknown.ages <- ((n + 1):(n + n.node))[-(node - n)]
- ## `basal' contains the indices of the basal edges (ie, linked to the root):
- basal <- which(phy$edge[, 1] == ROOT)
+ if (length(node) == 1 && node == "root") node <- ROOT
+ if (any(node <= n))
+ stop("node numbers should be greater than the number of tips")
+ zerobl <- which(phy$edge.length <= 0)
+ if (length(zerobl)) {
+ if (any(phy$edge[zerobl, 2] <= n))
+ stop("at least one terminal branch is of length zero:
+ you should remove it to have a meaningful estimation.")
+ else {
+ warning("at least one internal branch is of length zero:
+ it was collapsed and some nodes have been deleted.")
+ if (length(node) == 1 && node == ROOT)
+ phy <- di2multi(phy)
+ else {
+ tmp <- FALSE
+ if (is.null(phy$node.label)) {
+ tmp <- !tmp
+ phy$node.label <- paste("node", 1:phy$Nnode)
+ }
+ node.lab <- phy$node.label[node - n]
+ phy <- di2multi(phy)
+ node <- match(node.lab, phy$node.label) + n
+ if (tmp) phy$node.label <- NULL
+ }
+ }
+ }
+ m <- phy$Nnode
+ el <- phy$edge.length
+ e <- phy$edge
+ N <- dim(e)[1]
+ TIPS <- 1:n
+ EDGES <- 1:N
+
+ ini.rate <- el
+ el <- el/S
+
+ ## `basal' contains the indices of the basal edges
+ ## (ie, linked to the root):
+ basal <- which(e[, 1] == ROOT)
+ Nbasal <- length(basal)
## `ind' contains in its 1st column the index of all nonbasal
## edges, and in its second column the index of the edges
## | | c | a |
## |___c___ | | |
- ind <- matrix(NA, N - 2, 2)
- j <- 1
- for (i in 1:N) {
- if (phy$edge[i, 1] == ROOT) next
- ind[j, 1] <- i
- ind[j, 2] <- which(phy$edge[, 2] == phy$edge[i, 1])
- j <- j + 1
+ ind <- matrix(0L, N - Nbasal, 2)
+ ind[, 1] <- EDGES[-basal]
+ ind[, 2] <- match(e[EDGES[-basal], 1], e[, 2])
+
+ age <- numeric(n + m)
+
+ ##ini.time <- node.depth(phy)[-TIPS] - 1
+ ini.time <- node.depth(phy) - 1
+
+ ## first, rescale all times with respect to
+ ## the age of the 1st element of `node':
+ ratio <- age.min[1]/ini.time[node[1]]
+ ini.time <- ini.time*ratio
+
+ if (length(node) > 1) {
+ ini.time[node] <- age.min
+ real.edge.length <- ini.time[e[, 1]] - ini.time[e[, 2]]
+ while (any(real.edge.length <= 0)) {
+ for (i in EDGES) {
+ if (real.edge.length[i] <= 0) {
+ if (e[i, 1] %in% node) {
+ ini.time[e[i, 2]] <-
+ ini.time[e[, 2]] - 2*real.edge.length[i]
+ next
+ }
+ if (e[i, 2] %in% node) {
+ ini.time[e[i, 1]] <-
+ ini.time[e[, 1]] + 2*real.edge.length[i]
+ next
+ }
+ ini.time[e[i, 2]] <-
+ ini.time[e[, 2]] - real.edge.length[i]
+ ini.time[e[i, 1]] <-
+ ini.time[e[, 1]] + real.edge.length[i]
+ }
+ }
+ real.edge.length <- ini.time[e[, 1]] - ini.time[e[, 2]]
+ }
+ }
+
+ ## `unknown.ages' will contain the index of the nodes of unknown age:
+ unknown.ages <- n + 1:m
+
+ ## define the bounds for the node ages:
+ lower <- rep(tol, length(unknown.ages))
+ upper <- rep(1/tol, length(unknown.ages))
+
+ if (!is.null(age.max)) { # are some nodes known within some intervals?
+ lower[node - n] <- age.min
+ upper[node - n] <- age.max
+ interv <- which(age.min != age.max)
+ node <- node[-interv]
+ if (length(node)) age[node] <- age.min[-interv]
+ } else age[node] <- age.min
+
+ if (length(node)) {
+ unknown.ages <- unknown.ages[n - node]
+ lower <- lower[n - node]
+ upper <- upper[n - node]
}
- age <- rep(0, 2*n - 1)
- age[node] <- node.age
-
- tmp <- reorder(phy, "pruningwise")
- ini.time <- .C("node_depth", as.integer(n), as.integer(n.node),
- as.integer(tmp$edge[, 1]), as.integer(tmp$edge[, 2]),
- as.integer(N), double(n + n.node), DUP = FALSE,
- PACKAGE = "ape")[[6]][-(1:n)] - 1
- ini.time <- ini.time/max(ini.time)
- ini.time <- ini.time*node.age/ini.time[known.ages[-(1:n)] - n]
- ## check that there are no negative branch lengths:
- ini.time[known.ages[-(1:n)] - n] <- node.age
- it <- c(age[1:n], ini.time)
- ibl <- it[phy$edge[, 1]] - it[phy$edge[, 2]]
- if (any(ibl < 0)) {
- for (i in which(ibl < 0))
- if (phy$edge[i, 1] %in% node)
- ini.time[phy$edge[i, 2]] <- ini.time[phy$edge[i, 1]] - 1e-3
- else
- ini.time[phy$edge[i, 1]] <- ini.time[phy$edge[i, 2]] + 1e-3
+ ## `known.ages' contains the index of all nodes (internal and
+ ## terminal) of known age:
+ known.ages <- c(TIPS, node)
+
+ ## concatenate the bounds for the rates:
+ lower <- c(rep(tol, N), lower)
+ upper <- c(rep(1 - tol, N), upper)
+
+ minusploglik.gr <- function(rate, node.time) {
+ grad <- numeric(N + length(unknown.ages))
+ age[unknown.ages] <- node.time
+ real.edge.length <- age[e[, 1]] - age[e[, 2]]
+ if (any(real.edge.length < 0)) {
+ grad[] <- 0
+ return(grad)
+ }
+ ## gradient for the rates:
+ ## the parametric part can be calculated without a loop:
+ grad[EDGES] <- real.edge.length - el/rate
+ if (Nbasal == 2) { # the simpler formulae if there's a basal dichotomy
+ grad[basal[1]] <-
+ grad[basal[1]] + lambda*(rate[basal[1]] - rate[basal[2]])
+ grad[basal[2]] <-
+ grad[basal[2]] + lambda*(rate[basal[2]] - rate[basal[1]])
+ } else { # the general case
+ for (i in 1:Nbasal)
+ grad[basal[i]] <- grad[basal[i]] +
+ lambda*(2*rate[basal[i]]*(1 - 1/Nbasal) -
+ 2*sum(rate[basal[-i]])/Nbasal)/(Nbasal - 1)
+ }
+
+ for (i in EDGES) {
+ ii <- c(which(e[, 2] == e[i, 1]), which(e[, 1] == e[i, 2]))
+ if (!length(ii)) next
+ grad[i] <- grad[i] + lambda*(2*length(ii)*rate[i] - 2*sum(rate[ii]))
+ }
+
+ ## gradient for the 'node times'
+ for (i in 1:length(unknown.ages)) {
+ nd <- unknown.ages[i]
+ ii <- which(e[, 1] == nd)
+ grad[i + N] <-
+ sum(rate[ii] - el[ii]/real.edge.length[ii])#, na.rm = TRUE)
+ if (nd != ROOT) {
+ ii <- which(e[, 2] == nd)
+ grad[i + N] <- grad[i + N] -
+ rate[ii] + el[ii]/real.edge.length[ii]
+ }
+ }
+ grad
}
- ploglik <- function(rate, node.time) {
+ minusploglik <- function(rate, node.time) {
age[unknown.ages] <- node.time
- real.edge.length <- age[phy$edge[, 1]] - age[phy$edge[, 2]]
+ real.edge.length <- age[e[, 1]] - age[e[, 2]]
+ if (any(real.edge.length < 0)) return(1e50)
B <- rate*real.edge.length
- loglik <- sum(-B + phy$edge.length*log(B) -
- lfactorial(phy$edge.length))
- loglik - lambda * (sum((rate[ind[, 1]] - rate[ind[, 2]])^2)
- + var(rate[basal]))
+ loglik <- sum(-B + el*log(B) - lfactorial(el))
+ -(loglik - lambda*(sum((rate[ind[, 1]] - rate[ind[, 2]])^2)
+ + var(rate[basal])))
}
- out <- nlm(function(p) -ploglik(p[1:N], p[-(1:N)]),
- p = c(ini.rate, ini.time[unknown.ages - n]),
- iterlim = 500)
+ out <- nlminb(c(ini.rate, ini.time[unknown.ages]),
+ function(p) minusploglik(p[EDGES], p[-EDGES]),
+ function(p) minusploglik.gr(p[EDGES], p[-EDGES]),
+ control = list(eval.max = eval.max, iter.max = iter.max, ...),
+ lower = lower, upper = upper)
- attr(phy, "ploglik") <- -out$minimum
- attr(phy, "rates") <- out$estimate[1:N]
- age[unknown.ages] <- out$estimate[-(1:N)]
+ attr(phy, "ploglik") <- -out$objective
+ attr(phy, "rates") <- out$par[EDGES]
+ attr(phy, "message") <- out$message
+ age[unknown.ages] <- out$par[-EDGES]
if (CV) ophy <- phy
- phy$edge.length <- age[phy$edge[, 1]] - age[phy$edge[, 2]]
- if (CV)
- attr(phy, "D2") <-
- chronopl.cv(ophy, lambda, node.age, node, n)
+ phy$edge.length <- age[e[, 1]] - age[e[, 2]]
+ if (CV) attr(phy, "D2") <-
+ chronopl.cv(ophy, lambda, age.min, age.max, node,
+ n, S, tol, eval.max, iter.max, ...)
phy
}
-chronopl.cv <- function(ophy, lambda, node.age, nodes, n)
+chronopl.cv <- function(ophy, lambda, age.min, age.max, nodes,
+ n, S, tol, eval.max, iter.max, ...)
### ophy: the original phylogeny
### n: number of tips
### Note that we assume here that the order of the nodes
BT <- branching.times(ophy)
D2 <- numeric(n)
+ cat(" dropping tip")
for (i in 1:n) {
- cat(" dropping tip", i, "\n")
+ cat(" ", i, sep = "")
tr <- drop.tip(ophy, i)
j <- which(ophy$edge[, 2] == i)
if (ophy$edge[j, 1] %in% nodes) {
k <- which(nodes == ophy$edge[j, 1])
- nodes <- nodes[-k]
- node.age <- node.age[-k]
- }
- if (length(nodes)) {
- chr <- chronopl(tr, lambda, node.age, nodes)
- ## <FIXME> à vérifier:
- ## tmp <- BT[as.character(ophy$edge[j, 1])]
- tmp <- BT[-(ophy$edge[j, 1] - n)]
- ## </FIXME>
+ node <- nodes[-k]
+ agemin <- age.min[-k]
+ agemax <- age.max[-k]
+ } else node <- nodes
+ if (length(node)) {
+ chr <- chronopl(tr, lambda, age.min, age.max, node,
+ S, tol, FALSE, eval.max, iter.max, ...)
D2[i] <- sum((tmp - branching.times(chr))^2 / tmp)
} else D2[i] <- 0
}
+ cat("\n")
D2
}
\alias{chronopl}
\title{Molecular Dating With Penalized Likelihood}
\usage{
-chronopl(phy, lambda, node.age = 1, node = "root", CV = FALSE)
+chronopl(phy, lambda, age.min = 1, age.max = NULL,
+ node = "root", S = 1, tol = 1e-8,
+ CV = FALSE, eval.max = 500, iter.max = 500, ...)
}
\arguments{
\item{phy}{an object of class \code{"phylo"}.}
- \item{lambda}{value of the smoothng parameter.}
- \item{node.age}{numeric values specifying the fixed node ages.}
+ \item{lambda}{value of the smoothing parameter.}
+ \item{age.min}{numeric values specifying the fixed node ages (if
+ \code{age.max = NULL}) or the youngest bound of the nodes known to
+ be within an interval.}
+ \item{age.max}{numeric values specifying the oldest bound of the nodes
+ known to be within an interval.}
\item{node}{the numbers of the nodes whose ages are given by
- \code{node.age}; \code{"root"} is a short-cut the number of the
- node.}
+ \code{age.min}; \code{"root"} is a short-cut for the root.}
+ \item{S}{the number of sites in the sequences; leave the default if
+ branch lengths are in mean number of substitutions.}
+ \item{tol}{the value below which branch lengths are considered
+ effectively zero.}
\item{CV}{whether to perform cross-validation.}
+ \item{eval.max}{the maximal number of evaluations of the penalized
+ likelihood function.}
+ \item{iter.max}{the maximal number of iterations of the optimization
+ algorithm.}
+ \item{\dots}{further arguments passed to control \code{nlminb}.}
}
\description{
- This function estimates the node ages of a tree using semi-parametric
- method based on penalized likelihood (Sanderson 2002). The branch
- lengths of the input tree are interpreted as (mean) numbers of
- substitutions.
+ This function estimates the node ages of a tree using a
+ semi-parametric method based on penalized likelihood (Sanderson
+ 2002). The branch lengths of the input tree are interpreted as mean
+ numbers of substitutions (i.e., per site).
}
\details{
The idea of this method is to use a trade-off between a parametric
rate for all branches).
\code{lambda} must be given. The known ages are given in
- \code{node.age}, and the correponding node numbers in \code{node}.
+ \code{age.min}, and the correponding node numbers in \code{node}.
These two arguments must obviously be of the same length. By default,
an age of 1 is assumed for the root, and the ages of the other nodes
are estimated.
+ If \code{age.max = NULL} (the default), it is assumed that
+ \code{age.min} gives exactly known ages. Otherwise, \code{age.max} and
+ \code{age.min} must be of the same length and give the intervals for
+ each node. Some node may be known exactly while the others are
+ known within some bounds: the values will be identical in both
+ arguments for the former (e.g., \code{age.min = c(10, 5), age.max =
+ c(10, 6), node = c(15, 18)} means that the age of node 15 is 10
+ units of time, and the age of node 18 is between 5 and 6).
+
+ The input tree may have multichotomies. If some internal branches are
+ of zero-length, they are collapsed (with a warning), and the returned
+ tree will have less nodes than the input one. The presence of
+ zero-lengthed terminal branches of results in an error since it makes
+ little sense to have zero-rate branches.
+
The cross-validation used here is different from the one proposed by
Sanderson (2002). Here, each tip is dropped successively and the
analysis is repeated with the reduced tree: the estimated dates for
with the full phylogeny, \eqn{t_j^{-i}}{tj[-i]} is the estimated date
for the \eqn{j}{j}th node after removing tip \eqn{i}{i} from the tree,
and \eqn{n}{n} is the number of tips.
+
+ The present version uses the \code{\link[stats]{nlminb}} to optimise
+ the penalized likelihood function: see its help page for details on
+ parameters controlling the optimisation procedure.
}
\value{
an object of class \code{"phylo"} with branch lengths as estimated by
- the function. There are two or three further attributes:
+ the function. There are three or four further attributes:
\item{ploglik}{the maximum penalized log-likelihood.}
\item{rates}{the estimated rates for each branch.}
+ \item{message}{the message returned by \code{nlminb} indicating
+ whether the optimisation converged.}
\item{D2}{the influence of each observation on overall date
estimates (if \code{CV = TRUE}).}
}
projects / ape.git / commitdiff