1 ## chronopl.R (2012-02-09)
3 ## Molecular Dating With Penalized Likelihood
5 ## Copyright 2005-2012 Emmanuel Paradis
7 ## This file is part of the R-package `ape'.
8 ## See the file ../COPYING for licensing issues.
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, ...)
15 n <- length(phy$tip.label)
17 if (identical(node, "root")) node <- ROOT
19 stop("node numbers should be greater than the number of tips")
20 zerobl <- which(phy$edge.length <= 0)
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.")
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)
32 if (is.null(phy$node.label)) {
34 phy$node.label <- paste("node", 1:phy$Nnode)
36 node.lab <- phy$node.label[node - n]
38 node <- match(node.lab, phy$node.label) + n
39 if (tmp) phy$node.label <- NULL
54 ## `basal' contains the indices of the basal edges
55 ## (ie, linked to the root):
56 basal <- which(e1 == ROOT)
57 Nbasal <- length(basal)
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:
70 ind <- matrix(0L, N - Nbasal, 2)
71 ind[, 1] <- EDGES[-basal]
72 ind[, 2] <- match(e1[EDGES[-basal]], e2)
76 #############################################################################
77 ### This bit sets 'ini.time' and should result in no negative branch lengths
79 seq.nod <- .Call("seq_root2tip", phy$edge, n, phy$Nnode, PACKAGE = "ape")
82 ini.time[ROOT:(n + m)] <- NA
83 ini.time[node] <- if (is.null(age.max)) age.min else (age.min + age.max) / 2
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)
89 ISnotNA.ALL <- unlist(lapply(seq.nod, function(x) sum(!is.na(ini.time[x]))))
90 o <- order(ISnotNA.ALL, decreasing = TRUE)
92 for (y in seq.nod[o]) {
93 ISNA <- is.na(ini.time[y])
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
99 while (ISNA[j]) j <- j + 1L # look for the next non-NA
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)
109 real.edge.length <- ini.time[e1] - ini.time[e2]
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 #############################################################################
116 ## because if (!is.null(age.max)), 'node' is modified,
117 ## so we copy it in case CV = TRUE:
120 ## `unknown.ages' will contain the index of the nodes of unknown age:
121 unknown.ages <- n + 1:m
123 ## define the bounds for the node ages:
124 lower <- rep(tol, length(unknown.ages))
125 upper <- rep(1/tol, length(unknown.ages))
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
138 unknown.ages <- unknown.ages[n - node] # 'n - node' is simplification for '-(node - n)'
139 lower <- lower[n - node]
140 upper <- upper[n - node]
143 ## `known.ages' contains the index of all nodes (internal and
144 ## terminal) of known age:
145 known.ages <- c(TIPS, node)
147 ## concatenate the bounds for the rates:
148 lower <- c(rep(tol, N), lower)
149 upper <- c(rep(1 - tol, N), upper)
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)) {
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
164 grad[basal[1]] + lambda*(rate[basal[1]] - rate[basal[2]])
166 grad[basal[2]] + lambda*(rate[basal[2]] - rate[basal[1]])
167 } else { # the general case
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)
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]))
180 ## gradient for the 'node times'
181 for (i in 1:length(unknown.ages)) {
182 nd <- unknown.ages[i]
183 ii <- which(e1 == nd)
185 sum(rate[ii] - el[ii]/real.edge.length[ii])#, na.rm = TRUE)
187 ii <- which(e2 == nd)
188 grad[i + N] <- grad[i + N] -
189 rate[ii] + el[ii]/real.edge.length[ii]
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)
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)
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]
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, ...)
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
230 cat("Doing cross-validation\n")
231 BT <- branching.times(ophy)
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])
241 agemin <- age.min[-k]
242 agemax <- age.max[-k]
245 chr <- chronopl(tr, lambda, age.min, age.max, node,
246 S, tol, FALSE, eval.max, iter.max, ...)
248 if (Nnode(chr) == Nnode(ophy)) BT else BT[-(ophy$edge[j, 1] - n)]
249 D2[i] <- sum((tmp - branching.times(chr))^2 / tmp)