R/chronopl.R

   1 ## chronopl.R (2007-01-17)
   2
   3 ##   Molecular Dating With Penalized Likelihood
   4
   5 ## Copyright 2005-2007 Emmanuel Paradis
   6
   7 ## This file is part of the R-package `ape'.
   8 ## See the file ../COPYING for licensing issues.
   9
  10 chronopl <- function(phy, lambda, node.age = 1, node = "root",
  11                      CV = FALSE)
  12 {
  13     n <- length(phy$tip.label)
  14     n.node <- phy$Nnode
  15     if (n != n.node + 1)
  16       stop("the tree must be rooted AND dichotomous")
  17     if (any(phy$edge.length == 0))
  18       stop("some branch lengths are equal to zero;
  19 you must remove them beforehand.")
  20     N <- dim(phy$edge)[1]
  21     ROOT <- n + 1
  22     if (node == "root") node <- ROOT
  23     ini.rate <- phy$edge.length
  24     ## `known.ages' contains the index of all nodes (internal and
  25     ## terminal) of known age:
  26     known.ages <- c(1:n, node)
  27     ## `unknown.ages' contains the index of the nodes of unknown age:
  28     unknown.ages <- ((n + 1):(n + n.node))[-(node - n)]
  29     ## `basal' contains the indices of the basal edges (ie, linked to the root):
  30     basal <- which(phy$edge[, 1] == ROOT)
  31
  32     ## `ind' contains in its 1st column the index of all nonbasal
  33     ## edges, and in its second column the index of the edges
  34     ## where these edges come from (ie, this matrix contains pairs
  35     ## of contiguous edges), eg:
  36
  37     ##         ___b___    ind:
  38     ##        |           |   |   |
  39     ## ___a___|           | b | a |
  40     ##        |           | c | a |
  41     ##        |___c___    |   |   |
  42
  43     ind <- matrix(NA, N - 2, 2)
  44     j <- 1
  45     for (i in 1:N) {
  46         if (phy$edge[i, 1] == ROOT) next
  47         ind[j, 1] <- i
  48         ind[j, 2] <- which(phy$edge[, 2] == phy$edge[i, 1])
  49         j <- j + 1
  50     }
  51
  52     age <- rep(0, 2*n - 1)
  53     age[node] <- node.age
  54
  55     tmp <- reorder(phy, "pruningwise")
  56     ini.time <- .C("node_depth", as.integer(n), as.integer(n.node),
  57                    as.integer(tmp$edge[, 1]), as.integer(tmp$edge[, 2]),
  58                    as.integer(N), double(n + n.node), DUP = FALSE,
  59                    PACKAGE = "ape")[[6]][-(1:n)] - 1
  60     ini.time <- ini.time/max(ini.time)
  61     ini.time <- ini.time*node.age/ini.time[known.ages[-(1:n)] - n]
  62     ## check that there are no negative branch lengths:
  63     ini.time[known.ages[-(1:n)] - n] <- node.age
  64     it <- c(age[1:n], ini.time)
  65     ibl <- it[phy$edge[, 1]] - it[phy$edge[, 2]]
  66     if (any(ibl < 0)) {
  67         for (i in which(ibl < 0))
  68           if (phy$edge[i, 1] %in% node)
  69             ini.time[phy$edge[i, 2]] <- ini.time[phy$edge[i, 1]] - 1e-3
  70           else
  71             ini.time[phy$edge[i, 1]] <- ini.time[phy$edge[i, 2]] + 1e-3
  72     }
  73
  74     ploglik <- function(rate, node.time) {
  75         age[unknown.ages] <- node.time
  76         real.edge.length <- age[phy$edge[, 1]] - age[phy$edge[, 2]]
  77         B <- rate*real.edge.length
  78         loglik <- sum(-B + phy$edge.length*log(B) -
  79                       lfactorial(phy$edge.length))
  80         loglik - lambda * (sum((rate[ind[, 1]] - rate[ind[, 2]])^2)
  81                            + var(rate[basal]))
  82     }
  83
  84     out <- nlm(function(p) -ploglik(p[1:N], p[-(1:N)]),
  85                p = c(ini.rate, ini.time[unknown.ages - n]),
  86                iterlim = 500)
  87
  88     attr(phy, "ploglik") <- -out$minimum
  89     attr(phy, "rates") <- out$estimate[1:N]
  90     age[unknown.ages] <- out$estimate[-(1:N)]
  91     if (CV) ophy <- phy
  92     phy$edge.length <- age[phy$edge[, 1]] - age[phy$edge[, 2]]
  93     if (CV)
  94       attr(phy, "D2") <-
  95         chronopl.cv(ophy, lambda, node.age, node, n)
  96     phy
  97 }
  98
  99 chronopl.cv <- function(ophy, lambda, node.age, nodes, n)
 100 ### ophy: the original phylogeny
 101 ### n: number of tips
 102 ### Note that we assume here that the order of the nodes
 103 ### in node.label are not modified by the drop.tip operation
 104 {
 105     cat("Doing cross-validation\n")
 106     BT <- branching.times(ophy)
 107     D2 <- numeric(n)
 108
 109     for (i in 1:n) {
 110         cat("  dropping tip", i, "\n")
 111         tr <- drop.tip(ophy, i)
 112         j <- which(ophy$edge[, 2] == i)
 113         if (ophy$edge[j, 1] %in% nodes) {
 114             k <- which(nodes == ophy$edge[j, 1])
 115             nodes <- nodes[-k]
 116             node.age <- node.age[-k]
 117         }
 118         if (length(nodes)) {
 119             chr <- chronopl(tr, lambda, node.age, nodes)
 120             ## <FIXME> à vérifier:
 121             ## tmp <- BT[as.character(ophy$edge[j, 1])]
 122             tmp <- BT[-(ophy$edge[j, 1] - n)]
 123             ## </FIXME>
 124             D2[i] <- sum((tmp - branching.times(chr))^2 / tmp)
 125         } else D2[i] <- 0
 126     }
 127     D2
 128 }