R/chronopl.R

   1 ## chronopl.R (2011-07-04)
   2
   3 ##   Molecular Dating With Penalized Likelihood
   4
   5 ## Copyright 2005-2011 Emmanuel Paradis
   6
   7 ## This file is part of the R-package `ape'.
   8 ## See the file ../COPYING for licensing issues.
   9
  10 chronopl <-
  11     function(phy, lambda, age.min = 1, age.max = NULL,
  12              node = "root", S = 1, tol = 1e-8,
  13              CV = FALSE, eval.max = 500, iter.max = 500, ...)
  14 {
  15     n <- length(phy$tip.label)
  16     ROOT <- n + 1
  17     if (length(node) == 1 && node == "root") node <- ROOT
  18     if (any(node <= n))
  19         stop("node numbers should be greater than the number of tips")
  20     zerobl <- which(phy$edge.length <= 0)
  21     if (length(zerobl)) {
  22         if (any(phy$edge[zerobl, 2] <= n))
  23             stop("at least one terminal branch is of length zero:
  24   you should remove it to have a meaningful estimation.")
  25         else {
  26             warning("at least one internal branch is of length zero:
  27   it was collapsed and some nodes have been deleted.")
  28             if (length(node) == 1 && node == ROOT)
  29                 phy <- di2multi(phy)
  30             else {
  31                 tmp <- FALSE
  32                 if (is.null(phy$node.label)) {
  33                     tmp <- !tmp
  34                     phy$node.label <- paste("node", 1:phy$Nnode)
  35                 }
  36                 node.lab <- phy$node.label[node - n]
  37                 phy <- di2multi(phy)
  38                 node <- match(node.lab, phy$node.label) + n
  39                 if (tmp) phy$node.label <- NULL
  40             }
  41         }
  42     }
  43     m <- phy$Nnode
  44     el <- phy$edge.length
  45     e <- phy$edge
  46     N <- dim(e)[1]
  47     TIPS <- 1:n
  48     EDGES <- 1:N
  49
  50     ini.rate <- el
  51     el <- el/S
  52
  53     ## `basal' contains the indices of the basal edges
  54     ## (ie, linked to the root):
  55     basal <- which(e[, 1] == ROOT)
  56     Nbasal <- length(basal)
  57
  58     ## `ind' contains in its 1st column the index of all nonbasal
  59     ## edges, and in its second column the index of the edges
  60     ## where these edges come from (ie, this matrix contains pairs
  61     ## of contiguous edges), eg:
  62
  63     ##         ___b___    ind:
  64     ##        |           |   |   |
  65     ## ___a___|           | b | a |
  66     ##        |           | c | a |
  67     ##        |___c___    |   |   |
  68
  69     ind <- matrix(0L, N - Nbasal, 2)
  70     ind[, 1] <- EDGES[-basal]
  71     ind[, 2] <- match(e[EDGES[-basal], 1], e[, 2])
  72
  73     age <- numeric(n + m)
  74
  75     ##ini.time <- node.depth(phy)[-TIPS] - 1
  76     ini.time <- node.depth(phy) - 1
  77
  78     ## first, rescale all times with respect to
  79     ## the age of the 1st element of `node':
  80     ratio <- age.min[1]/ini.time[node[1]]
  81     ini.time <- ini.time*ratio
  82
  83     ## because if (!is.null(age.max)), 'node' is modified, so we copy it in case CV = TRUE:
  84     node.bak <- node
  85
  86     if (length(node) > 1) {
  87         ini.time[node] <- age.min
  88         real.edge.length <- ini.time[e[, 1]] - ini.time[e[, 2]]
  89         while (any(real.edge.length <= 0)) {
  90             for (i in EDGES) {
  91                 if (real.edge.length[i] > 0) next
  92                 if (e[i, 1] %in% node) {
  93                     ini.time[e[i, 2]] <- ini.time[e[1, 2]] - 2 * real.edge.length[i]
  94                     next
  95                 }
  96                 if (e[i, 2] %in% node) {
  97                     ini.time[e[i, 1]] <- ini.time[e[1, 1]] + 2 * real.edge.length[i]
  98                     next
  99                 }
 100                 ##browser()
 101                 ini.time[e[i, 2]] <- ini.time[e[1, 2]] - real.edge.length[i]
 102                 ini.time[e[i, 1]] <- ini.time[e[1, 1]] + real.edge.length[i]
 103             }
 104             real.edge.length <- ini.time[e[, 1]] - ini.time[e[, 2]]
 105             ##print(min(real.edge.length))
 106         }
 107     }
 108     ## `unknown.ages' will contain the index of the nodes of unknown age:
 109     unknown.ages <- n + 1:m
 110
 111     ## define the bounds for the node ages:
 112     lower <- rep(tol, length(unknown.ages))
 113     upper <- rep(1/tol, length(unknown.ages))
 114
 115     if (!is.null(age.max)) { # are some nodes known within some intervals?
 116         lower[node - n] <- age.min
 117         upper[node - n] <- age.max
 118         interv <- which(age.min != age.max)
 119         node <- node[-interv]
 120         if (length(node)) age[node] <- age.min[-interv]
 121     } else age[node] <- age.min
 122
 123     if (length(node)) {
 124         unknown.ages <- unknown.ages[n - node]
 125         lower <- lower[n - node]
 126         upper <- upper[n - node]
 127     }
 128
 129     ## `known.ages' contains the index of all nodes (internal and
 130     ## terminal) of known age:
 131     known.ages <- c(TIPS, node)
 132
 133     ## concatenate the bounds for the rates:
 134     lower <- c(rep(tol, N), lower)
 135     upper <- c(rep(1 - tol, N), upper)
 136
 137     minusploglik.gr <- function(rate, node.time) {
 138         grad <- numeric(N + length(unknown.ages))
 139         age[unknown.ages] <- node.time
 140         real.edge.length <- age[e[, 1]] - age[e[, 2]]
 141         if (any(real.edge.length < 0)) {
 142             grad[] <- 0
 143             return(grad)
 144         }
 145         ## gradient for the rates:
 146         ## the parametric part can be calculated without a loop:
 147         grad[EDGES] <- real.edge.length - el/rate
 148         if (Nbasal == 2) { # the simpler formulae if there's a basal dichotomy
 149             grad[basal[1]] <-
 150                 grad[basal[1]] + lambda*(rate[basal[1]] - rate[basal[2]])
 151             grad[basal[2]] <-
 152                 grad[basal[2]] + lambda*(rate[basal[2]] - rate[basal[1]])
 153         } else { # the general case
 154             for (i in 1:Nbasal)
 155                 grad[basal[i]] <- grad[basal[i]] +
 156                     lambda*(2*rate[basal[i]]*(1 - 1/Nbasal) -
 157                             2*sum(rate[basal[-i]])/Nbasal)/(Nbasal - 1)
 158         }
 159
 160         for (i in EDGES) {
 161             ii <- c(which(e[, 2] == e[i, 1]), which(e[, 1] == e[i, 2]))
 162             if (!length(ii)) next
 163             grad[i] <- grad[i] + lambda*(2*length(ii)*rate[i] - 2*sum(rate[ii]))
 164         }
 165
 166         ## gradient for the 'node times'
 167         for (i in 1:length(unknown.ages)) {
 168             nd <- unknown.ages[i]
 169             ii <- which(e[, 1] == nd)
 170             grad[i + N] <-
 171                 sum(rate[ii] - el[ii]/real.edge.length[ii])#, na.rm = TRUE)
 172             if (nd != ROOT) {
 173                 ii <- which(e[, 2] == nd)
 174                 grad[i + N] <- grad[i + N] -
 175                     rate[ii] + el[ii]/real.edge.length[ii]
 176             }
 177         }
 178         grad
 179     }
 180
 181     minusploglik <- function(rate, node.time) {
 182         age[unknown.ages] <- node.time
 183         real.edge.length <- age[e[, 1]] - age[e[, 2]]
 184         if (any(real.edge.length < 0)) return(1e50)
 185         B <- rate*real.edge.length
 186         loglik <- sum(-B + el*log(B) - lfactorial(el))
 187         -(loglik - lambda*(sum((rate[ind[, 1]] - rate[ind[, 2]])^2)
 188                            + var(rate[basal])))
 189     }
 190
 191     out <- nlminb(c(ini.rate, ini.time[unknown.ages]),
 192                   function(p) minusploglik(p[EDGES], p[-EDGES]),
 193                   function(p) minusploglik.gr(p[EDGES], p[-EDGES]),
 194                   control = list(eval.max = eval.max, iter.max = iter.max, ...),
 195                   lower = lower, upper = upper)
 196
 197     attr(phy, "ploglik") <- -out$objective
 198     attr(phy, "rates") <- out$par[EDGES]
 199     attr(phy, "message") <- out$message
 200     age[unknown.ages] <- out$par[-EDGES]
 201     if (CV) ophy <- phy
 202     phy$edge.length <- age[e[, 1]] - age[e[, 2]]
 203     if (CV) attr(phy, "D2") <-
 204         chronopl.cv(ophy, lambda, age.min, age.max, node.bak,
 205                     n, S, tol, eval.max, iter.max, ...)
 206     phy
 207 }
 208
 209 chronopl.cv <- function(ophy, lambda, age.min, age.max, nodes,
 210                         n, S, tol, eval.max, iter.max, ...)
 211 ### ophy: the original phylogeny
 212 ### n: number of tips
 213 ### Note that we assume here that the order of the nodes
 214 ### in node.label are not modified by the drop.tip operation
 215 {
 216     cat("Doing cross-validation\n")
 217     BT <- branching.times(ophy)
 218     D2 <- numeric(n)
 219
 220     for (i in 1:n) {
 221         cat("\r  dropping tip ", i, " / ", n, sep = "")
 222         tr <- drop.tip(ophy, i)
 223         j <- which(ophy$edge[, 2] == i)
 224         if (ophy$edge[j, 1] %in% nodes) {
 225             k <- which(nodes == ophy$edge[j, 1])
 226             node <- nodes[-k]
 227             agemin <- age.min[-k]
 228             agemax <- age.max[-k]
 229         } else node <- nodes
 230         if (length(node)) {
 231             chr <- chronopl(tr, lambda, age.min, age.max, node,
 232                             S, tol, FALSE, eval.max, iter.max, ...)
 233             tmp <-
 234                 if (Nnode(chr) == Nnode(ophy)) BT else BT[-(ophy$edge[j, 1] - n)]
 235             D2[i] <- sum((tmp - branching.times(chr))^2 / tmp)
 236         } else D2[i] <- 0
 237     }
 238     cat("\n")
 239     D2
 240 }