R/chronopl.R

   1 ## chronopl.R (2012-02-09)
   2
   3 ##   Molecular Dating With Penalized Likelihood
   4
   5 ## Copyright 2005-2012 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 + 1L
  17     if (identical(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     e1 <- phy$edge[, 1L]
  46     e2 <- phy$edge[, 2L]
  47     N <- length(e1)
  48     TIPS <- 1:n
  49     EDGES <- 1:N
  50
  51     ini.rate <- el
  52     el <- el/S
  53
  54     ## `basal' contains the indices of the basal edges
  55     ## (ie, linked to the root):
  56     basal <- which(e1 == ROOT)
  57     Nbasal <- length(basal)
  58
  59     ## `ind' contains in its 1st column the index of all nonbasal
  60     ## edges, and in its second column the index of the edges
  61     ## where these edges come from (ie, this matrix contains pairs
  62     ## of contiguous edges), eg:
  63
  64     ##         ___b___    ind:
  65     ##        |           |   |   |
  66     ## ___a___|           | b | a |
  67     ##        |           | c | a |
  68     ##        |___c___    |   |   |
  69
  70     ind <- matrix(0L, N - Nbasal, 2)
  71     ind[, 1] <- EDGES[-basal]
  72     ind[, 2] <- match(e1[EDGES[-basal]], e2)
  73
  74     age <- numeric(n + m)
  75
  76 #############################################################################
  77 ### This bit sets 'ini.time' and should result in no negative branch lengths
  78
  79     seq.nod <- .Call("seq_root2tip", phy$edge, n, phy$Nnode, PACKAGE = "ape")
  80
  81     ini.time <- age
  82     ini.time[ROOT:(n + m)] <- NA
  83     ini.time[node] <- if (is.null(age.max)) age.min else (age.min + age.max) / 2
  84
  85     ## if no age given for the root, find one approximately:
  86     if (is.na(ini.time[ROOT]))
  87         ini.time[ROOT] <- if (is.null(age.max)) 3 * max(age.min) else 3 * max(age.max)
  88
  89     ISnotNA.ALL <- unlist(lapply(seq.nod, function(x) sum(!is.na(ini.time[x]))))
  90     o <- order(ISnotNA.ALL, decreasing = TRUE)
  91
  92     for (y in seq.nod[o]) {
  93         ISNA <- is.na(ini.time[y])
  94         if (any(ISNA)) {
  95             i <- 2L # we know the 1st value is not NA, so we start at the 2nd one
  96             while (i <= length(y)) {
  97                 if (ISNA[i]) { # we stop at the next NA
  98                     j <- i + 1L
  99                     while (ISNA[j]) j <- j + 1L # look for the next non-NA
 100                     nb.val <- j - i
 101                     by <- (ini.time[y[i - 1L]] - ini.time[y[j]]) / (nb.val + 1)
 102                     ini.time[y[i:(j - 1L)]] <- ini.time[y[i - 1L]] - by * seq_len(nb.val)
 103                     i <- j + 1L
 104                 } else i <- i + 1L
 105             }
 106         }
 107     }
 108
 109     real.edge.length <- ini.time[e1] - ini.time[e2]
 110
 111     if (any(real.edge.length <= 0))
 112         stop("some initial branch lengths are zero or negative;
 113   maybe you need to adjust the given dates -- see '?chronopl' for details")
 114 #############################################################################
 115
 116     ## because if (!is.null(age.max)), 'node' is modified,
 117     ## so we copy it in case CV = TRUE:
 118     node.bak <- node
 119
 120     ## `unknown.ages' will contain the index of the nodes of unknown age:
 121     unknown.ages <- n + 1:m
 122
 123     ## define the bounds for the node ages:
 124     lower <- rep(tol, length(unknown.ages))
 125     upper <- rep(1/tol, length(unknown.ages))
 126
 127     if (!is.null(age.max)) { # are some nodes known within some intervals?
 128         lower[node - n] <- age.min
 129         upper[node - n] <- age.max
 130         ## find nodes known within an interval:
 131         interv <- which(age.min != age.max)
 132         ## drop them from the 'node' since they will be estimated:
 133         node <- node[-interv]
 134         if (length(node)) age[node] <- age.min[-interv] # update 'age'
 135     } else age[node] <- age.min
 136
 137     if (length(node)) {
 138         unknown.ages <- unknown.ages[n - node] # 'n - node' is simplification for '-(node - n)'
 139         lower <- lower[n - node]
 140         upper <- upper[n - node]
 141     }
 142
 143     ## `known.ages' contains the index of all nodes (internal and
 144     ## terminal) of known age:
 145     known.ages <- c(TIPS, node)
 146
 147     ## concatenate the bounds for the rates:
 148     lower <- c(rep(tol, N), lower)
 149     upper <- c(rep(1 - tol, N), upper)
 150
 151     minusploglik.gr <- function(rate, node.time) {
 152         grad <- numeric(N + length(unknown.ages))
 153         age[unknown.ages] <- node.time
 154         real.edge.length <- age[e1] - age[e2]
 155         if (any(real.edge.length < 0)) {
 156             grad[] <- 0
 157             return(grad)
 158         }
 159         ## gradient for the rates:
 160         ## the parametric part can be calculated without a loop:
 161         grad[EDGES] <- real.edge.length - el/rate
 162         if (Nbasal == 2) { # the simpler formulae if there's a basal dichotomy
 163             grad[basal[1]] <-
 164                 grad[basal[1]] + lambda*(rate[basal[1]] - rate[basal[2]])
 165             grad[basal[2]] <-
 166                 grad[basal[2]] + lambda*(rate[basal[2]] - rate[basal[1]])
 167         } else { # the general case
 168             for (i in 1:Nbasal)
 169                 grad[basal[i]] <- grad[basal[i]] +
 170                     lambda*(2*rate[basal[i]]*(1 - 1/Nbasal) -
 171                             2*sum(rate[basal[-i]])/Nbasal)/(Nbasal - 1)
 172         }
 173
 174         for (i in EDGES) {
 175             ii <- c(which(e2 == e1[i]), which(e1 == e2[i]))
 176             if (!length(ii)) next
 177             grad[i] <- grad[i] + lambda*(2*length(ii)*rate[i] - 2*sum(rate[ii]))
 178         }
 179
 180         ## gradient for the 'node times'
 181         for (i in 1:length(unknown.ages)) {
 182             nd <- unknown.ages[i]
 183             ii <- which(e1 == nd)
 184             grad[i + N] <-
 185                 sum(rate[ii] - el[ii]/real.edge.length[ii])#, na.rm = TRUE)
 186             if (nd != ROOT) {
 187                 ii <- which(e2 == nd)
 188                 grad[i + N] <- grad[i + N] -
 189                     rate[ii] + el[ii]/real.edge.length[ii]
 190             }
 191         }
 192         grad
 193     }
 194
 195     minusploglik <- function(rate, node.time) {
 196         age[unknown.ages] <- node.time
 197         real.edge.length <- age[e1] - age[e2]
 198         if (any(real.edge.length < 0)) return(1e50)
 199         B <- rate*real.edge.length
 200         loglik <- sum(-B + el*log(B) - lfactorial(el))
 201         -(loglik - lambda*(sum((rate[ind[, 1]] - rate[ind[, 2]])^2)
 202                            + var(rate[basal])))
 203     }
 204
 205     out <- nlminb(c(ini.rate, ini.time[unknown.ages]),
 206                   function(p) minusploglik(p[EDGES], p[-EDGES]),
 207                   function(p) minusploglik.gr(p[EDGES], p[-EDGES]),
 208                   control = list(eval.max = eval.max, iter.max = iter.max, ...),
 209                   lower = lower, upper = upper)
 210
 211     attr(phy, "ploglik") <- -out$objective
 212     attr(phy, "rates") <- out$par[EDGES]
 213     attr(phy, "message") <- out$message
 214     age[unknown.ages] <- out$par[-EDGES]
 215     if (CV) ophy <- phy
 216     phy$edge.length <- age[e1] - age[e2]
 217     if (CV) attr(phy, "D2") <-
 218         chronopl.cv(ophy, lambda, age.min, age.max, node.bak,
 219                     n, S, tol, eval.max, iter.max, ...)
 220     phy
 221 }
 222
 223 chronopl.cv <- function(ophy, lambda, age.min, age.max, nodes,
 224                         n, S, tol, eval.max, iter.max, ...)
 225 ### ophy: the original phylogeny
 226 ### n: number of tips
 227 ### Note that we assume here that the order of the nodes
 228 ### in node.label are not modified by the drop.tip operation
 229 {
 230     cat("Doing cross-validation\n")
 231     BT <- branching.times(ophy)
 232     D2 <- numeric(n)
 233
 234     for (i in 1:n) {
 235         cat("\r  dropping tip ", i, " / ", n, sep = "")
 236         tr <- drop.tip(ophy, i)
 237         j <- which(ophy$edge[, 2] == i)
 238         if (ophy$edge[j, 1] %in% nodes) {
 239             k <- which(nodes == ophy$edge[j, 1])
 240             node <- nodes[-k]
 241             agemin <- age.min[-k]
 242             agemax <- age.max[-k]
 243         } else node <- nodes
 244         if (length(node)) {
 245             chr <- chronopl(tr, lambda, age.min, age.max, node,
 246                             S, tol, FALSE, eval.max, iter.max, ...)
 247             tmp <-
 248                 if (Nnode(chr) == Nnode(ophy)) BT else BT[-(ophy$edge[j, 1] - n)]
 249             D2[i] <- sum((tmp - branching.times(chr))^2 / tmp)
 250         } else D2[i] <- 0
 251     }
 252     cat("\n")
 253     D2
 254 }