1 ## CDF.birth.death.R (2012-09-14)
3 ## Functions to Simulate and Fit
4 ## Time-Dependent Birth-Death Models
6 ## Copyright 2010-2012 Emmanuel Paradis
8 ## This file is part of the R-package `ape'.
9 ## See the file ../COPYING for licensing issues.
11 integrateTrapeze <- function(FUN, from, to, nint = 10)
12 ## compute an integral with a simple trapeze method
13 ## (apparently, Vectorize doesn't give faster calculation)
15 x <- seq(from = from, to = to, length.out = nint + 1)
16 ## reorganized to minimize the calls to FUN:
17 out <- FUN(x[1]) + FUN(x[length(x)])
18 for (i in 2:nint) out <- out + 2 * FUN(x[i])
19 (x[2] - x[1]) * out/2 # (x[2] - x[1]) is the width of the trapezes
23 ## 1: birth and death rates constant
24 ## 2: no primitive available
25 ## 3: primitives are available
26 ## 4: death rate constant, no primitive available
27 ## 5: birth rate constant, no primitive available
28 ## 6: death rate constant, primitive available for birth(t)
29 ## 7: birth rate constant, primitive available for death(t)
31 .getCase <- function(birth, death, BIRTH = NULL, DEATH = NULL)
33 if (is.numeric(birth)) {
34 if (is.numeric(death)) 1 else {
35 if (is.null(DEATH)) 5 else 7
38 if (is.numeric(death)) {
39 if (is.null(BIRTH)) 4 else 6
40 } else if (is.null(BIRTH) || is.null(DEATH)) 2 else 3
44 .getRHOetINT <- function(birth, death, BIRTH = NULL, DEATH = NULL, case, fast)
46 ## build the RHO(), \rho(t), and INT(), I(t), functions
49 RHO <- function(t1, t2) (t2 - t1)*(death - birth)
52 death*(exp(rho*(Tmax - t)) - 1)/rho
56 RHO <- function(t1, t2)
57 integrateTrapeze(function(t) death(t) - birth(t), t1, t2)
59 FOO <- function(u) exp(RHO(t, u)) * death(u)
60 integrateTrapeze(FOO, t, Tmax)
63 RHO <- function(t1, t2)
64 integrate(function(t) death(t) - birth(t), t1, t2)$value
66 FOO <- function(u) exp(RHO(t, u)) * death(u)
67 integrate(Vectorize(FOO), t, Tmax)$value # Vectorize required
71 RHO <- function(t1, t2)
72 DEATH(t2) - BIRTH(t2) - DEATH(t1) + BIRTH(t1)
73 INT <- function(t) { # vectorized
74 FOO <- function(u) exp(RHO(tt, u)) * death(u)
76 for (i in 1:length(t)) {
78 out[i] <- integrate(FOO, tt, Tmax)$value
84 RHO <- function(t1, t2)
85 death * (t2 - t1) - integrateTrapeze(birth, t1, t2)
87 FOO <- function(u) exp(RHO(t, u)) * death
88 integrateTrapeze(Vectorize(FOO), t, Tmax)
91 RHO <- function(t1, t2)
92 death * (t2 - t1) - integrate(birth, t1, t2)$value
94 FOO <- function(u) exp(RHO(t, u)) * death
95 integrate(Vectorize(FOO), t, Tmax)$value
99 RHO <- function(t1, t2)
100 integrate(death, t1, t2)$value - birth * (t2 - t1)
103 FOO <- function(u) exp(RHO(t, u)) * death(u)
104 integrateTrapeze(FOO, t, Tmax)
108 FOO <- function(u) exp(RHO(t, u)) * death(u)
109 integrate(Vectorize(FOO), t, Tmax)$value
113 RHO <- function(t1, t2) death * (t2 - t1) - BIRTH(t2) + BIRTH(t1)
114 INT <- function(t) { # vectorized
115 FOO <- function(u) exp(RHO(tt, u)) * death
117 for (i in 1:length(t)) {
119 out[i] <- integrate(FOO, tt, Tmax)$value
124 RHO <- function(t1, t2) DEATH(t2) - DEATH(t1) - birth * (t2 - t1)
127 FOO <- function(u) exp(RHO(t, u)) * death(u)
128 integrateTrapeze(FOO, t, Tmax)
132 FOO <- function(u) exp(RHO(t, u)) * death(u)
133 integrate(Vectorize(FOO), t, Tmax)$value
141 function(birth, death, BIRTH = NULL, DEATH = NULL, Tmax, x, case, fast = FALSE)
143 ff <- .getRHOetINT(birth, death, BIRTH, DEATH, case, fast)
146 environment(INT) <- environment() # so that INT() can find Tmax
147 .CDF.birth.death2(RHO, INT, birth, death, BIRTH, DEATH,
152 function(RHO, INT, birth, death, BIRTH, DEATH, Tmax, x, case, fast)
154 Pi <- if (case %in% c(1, 5, 7))
155 function(t) (1/(1 + INT(t)))^2 * 2 * exp(-RHO(0, t)) * birth
157 function(t) (1/(1 + INT(t)))^2 * 2 * exp(-RHO(0, t)) * birth(t)
159 if (!case %in% c(1, 3, 6)) Pi <- Vectorize(Pi)
162 if (fast) integrateTrapeze(Pi, 0, Tmax)
163 else integrate(Pi, 0, Tmax)$value
167 for (i in 1:n) p[i] <- integrateTrapeze(Pi, 0, x[i])
169 for (i in 1:n) p[i] <- integrate(Pi, 0, x[i])$value
174 .makePhylo <- function(edge, edge.length, i)
177 edge[NODES] <- edge[NODES] + i + 1L
178 edge[!NODES] <- 1:(i + 1L)
180 phy <- list(edge = edge, edge.length = edge.length,
181 tip.label = paste("t", 1:(i + 1), sep = ""), Nnode = i)
182 class(phy) <- "phylo"
183 attr(phy, "order") <- "cladewise"
188 function(birth, death, Tmax = 50, BIRTH = NULL, DEATH = NULL, eps = 1e-6)
190 case <- .getCase(birth, death, BIRTH, DEATH)
192 rTimeToEvent <- function(t)
194 ## CDF of the times to event (speciation or extinction):
197 Foo <- function(t, x)
198 1 - exp(-(birth + death)*(x - t))
200 Foo <- function(t, x) {
201 if (t == x) return(0)
202 1 - exp(-integrate(function(t) birth(t) + death(t),
206 Foo <- function(t, x) {
207 if (t == x) return(0)
208 1 - exp(-(BIRTH(x) - BIRTH(t) + DEATH(x) - DEATH(t)))
211 Foo <- function(t, x) {
212 if (t == x) return(0)
213 1 - exp(-(integrate(function(t) birth(t), t, x)$value
218 Foo <- function(t, x) {
219 if (t == x) return(0)
220 1 - exp(-(birth*(x - t) +
221 integrate(function(t) death(t), t, x)$value))
225 Foo <- function(t, x) {
226 if (t == x) return(0)
227 1 - exp(-(BIRTH(x) - BIRTH(t) + death*(x - t)))
231 Foo <- function(t, x) {
232 if (t == x) return(0)
233 1 - exp(-(birth*(x - t) + DEATH(x) - DEATH(t) ))
237 ## generate a random time to event by the inverse method:
239 ## in case speciation probability is so low
240 ## that time to speciation is infinite:
241 if (Foo(t, Tmax) < P) return(Tmax + 1)
244 while (inc > eps) { # la précision influe sur le temps de calcul
254 speORext <- function(t) birth/(birth + death)
255 if (case == 2 || case == 3)
256 speORext <- function(t) birth(t)/(birth(t) + death(t))
257 if (case == 4 || case == 6)
258 speORext <- function(t) birth(t)/(birth(t) + death)
259 if (case == 5 || case == 7)
260 speORext <- function(t) birth/(birth + death(t))
262 ## the recursive function implementing algorithm 1
263 foo <- function(node) {
265 X <- rTimeToEvent(t[node])
267 ## is the event a speciation or an extinction?
271 } else Y <- rbinom(1, size = 1, prob = speORext(tmp))
273 edge.length[j] <<- tmp - t[node]
277 ## set internal edge:
278 edge[j, ] <<- c(node, i)
281 ## set terminal edge:
282 edge[j, ] <<- c(node, 0L)
286 edge <- matrix(0L, 1e5, 2)
287 edge.length <- numeric(1e5)
288 j <- 0L; i <- 1; t <- 0
290 .makePhylo(edge[1:j, ], edge.length[1:j], i)
293 drop.fossil <- function(phy, tol = 1e-8)
296 x <- dist.nodes(phy)[n + 1, ][1:n]
297 drop.tip(phy, which(x < max(x) - tol))
301 function(birth, death, Tmax = 50, BIRTH = NULL, DEATH = NULL, eps = 1e-6)
303 case <- .getCase(birth, death, BIRTH, DEATH)
304 ff <- .getRHOetINT(birth, death, BIRTH, DEATH, case, FALSE)
307 ## so that RHO() and INT() can find Tmax:
308 environment(RHO) <- environment(INT) <- environment()
310 rtimetospe <- function(t)
312 ## CDF of the times to speciation:
313 Foo <- if (case %in% c(1, 5, 7))
314 function(t, x) 1 - exp(-birth*(x - t))
316 if (case %in% c(3, 6))
317 function(t, x) 1 - exp(-(BIRTH(x) - BIRTH(t)))
320 if (t == x) return(0)
321 1 - exp(-integrate(birth, t, x)$value)
325 ## generate a random time to speciation by the inverse method:
327 ## in case speciation probability is so low
328 ## that time to speciation is infinite:
329 if (Foo(t, Tmax) < P) return(Tmax + 1)
332 while (inc > eps) { # la précision influe sur le temps de calcul
341 ## the recursive function implementing algorithm 2
342 foo <- function(node, start) {
343 node <- node # make a local copy
345 tau <- start # because tau is changed below
348 while (X < Tmax - tau) {
350 ## does the new lineage survive until Tmax?
351 Y <- rbinom(1, size = 1, prob = 1/(1 + INT(tau)))
355 ## set internal edge:
357 edge[j, ] <<- c(node, i)
358 edge.length[j] <<- tau - t[node]
365 ## set terminal edge:
368 edge[j, 1] <<- node # the 2nd column is already set to 0
369 edge.length[j] <<- Tmax - t[node]
374 edge <- matrix(0L, 1e5, 2)
375 edge.length <- numeric(1e5)
376 j <- 0L; i <- 1L; t <- 0
378 .makePhylo(edge[1:j, ], edge.length[1:j], i)
381 bd.time <- function(phy, birth, death, BIRTH = NULL, DEATH = NULL,
382 ip, lower, upper, fast = FALSE, boot = 0, trace = 0)
384 guess.bounds <- if (missing(lower)) TRUE else FALSE
386 PrNt <- function(t, T, x)
388 tmp <- exp(-RHO(t, T))
389 Wt <- tmp * (1 + INT(t))
390 out <- numeric(length(x))
393 out[zero] <- 1 - tmp/Wt
394 out[!zero] <- (tmp/Wt^2)*(1 - 1/Wt)^(x[!zero] - 1)
395 } else out[] <- (tmp/Wt^2)*(1 - 1/Wt)^(x - 1)
399 case <- .getCase(birth, death, BIRTH, DEATH)
401 if (is.function(birth)) {
402 paranam <- names(formals(birth))
404 upper <- rep(BIG, length(paranam))
407 formals(birth) <- alist(t=)
408 environment(birth) <- environment()
409 if (!is.null(BIRTH)) environment(BIRTH) <- environment()
418 if (is.function(death)) {
419 tmp <- names(formals(death))
422 upper <- c(upper, rep(BIG, np2))
423 lower <- c(lower, rep(-BIG, np2))
425 paranam <- c(paranam, tmp)
426 formals(death) <- alist(t=)
427 environment(death) <- environment()
428 if (!is.null(DEATH)) environment(DEATH) <- environment()
430 paranam <- c(paranam, "death")
432 upper <- c(upper, .1)
437 np <- length(paranam)
439 ff <- .getRHOetINT(birth, death, BIRTH, DEATH, case = case, fast = fast)
442 environment(RHO) <- environment(INT) <- environment()
444 x <- branching.times(phy)
447 x <- Tmax - x # change the time scale so the root is t=0
450 foo <- function(para) {
452 assign(paranam[i], para[i], pos = sys.frame(1))
453 p <- CDF.birth.death(birth, death, BIRTH, DEATH, Tmax = Tmax,
454 x = x, case = case, fast = fast)
455 ## w is the probability of the observed tree size (= number of tips)
456 w <- PrNt(0, Tmax, Ntip(phy))
457 ## p is the expected CDF of branching times
458 ## ecdf(x)(x) is the observed CDF
459 sum((1:n/n - p)^2)/w # faster than sum((ecdf(x)(x) - p)^2)/w
462 if (missing(ip)) ip <- (upper - lower)/2
464 out <- nlminb(ip, foo, control = list(trace = trace, eval.max = 500),
465 upper = upper, lower = lower)
467 names(out$par) <- paranam
468 names(out)[2] <- "SS"
470 if (boot) { # nonparametric version
471 PAR <- matrix(NA, boot, np)
474 cat("\rDoing bootstrap no.", i, "\n")
475 x <- sort(sample(x, replace = TRUE))
476 o <- try(nlminb(ip, foo, control = list(trace = 0, eval.max = 500),
477 upper = upper, lower = lower))
478 if (class(o) == "list") {