1 ## CDF.birth.death.R (2010-08-09)
3 ## Functions to simulate and fit
4 ## Time-Dependent Birth-Death Models
6 ## Copyright 2010 Emmanuel Paradis
8 ## This file is part of the R-package `ape'.
9 ## See the file ../COPYING for licensing issues.
11 ## compute an integral with a simple trapeze method
12 ## (apparently, Vectorize doesn't give faster calculation)
13 integrateTrapeze <- function(FUN, from, to, nint = 10)
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 ###.getRHO <- function(birth, death, BIRTH = NULL, DEATH = NULL, case)
46 ### ## build the RHO(), \rho(t), function
48 ### function(t1, t2) (t2 - t1)*(death - birth), # case 1
49 ### function(t1, t2) # case 2
50 ### integrate(function(t) death(t) - birth(t), t1, t2)$value,
52 ### DEATH(t2) - BIRTH(t2) - DEATH(t1) + BIRTH(t1), # case 3
54 ### death - integrate(birth, t1, t2)$value, # case 4
56 ### integrate(death, t1, t2)$value - birth, # case 5
57 ### function(t1, t2) death - BIRTH(t2) + BIRTH(t1), # case 6
58 ### function(t1, t2) DEATH(t2) - DEATH(t1) - birth # case 7
62 .getRHOetINT <- function(birth, death, BIRTH = NULL, DEATH = NULL, case, fast)
64 ## build the RHO(), \rho(t), and INT(), I(t), functions
67 RHO <- function(t1, t2) (t2 - t1)*(death - birth)
70 death*(exp(rho*(Tmax - t)) - 1)/rho
74 RHO <- function(t1, t2)
75 integrateTrapeze(function(t) death(t) - birth(t), t1, t2)
77 FOO <- function(u) exp(RHO(t, u)) * death(u)
78 integrateTrapeze(FOO, t, Tmax)
81 RHO <- function(t1, t2)
82 integrate(function(t) death(t) - birth(t), t1, t2)$value
84 FOO <- function(u) exp(RHO(t, u)) * death(u)
85 integrate(Vectorize(FOO), t, Tmax)$value # Vectorize required
89 RHO <- function(t1, t2)
90 DEATH(t2) - BIRTH(t2) - DEATH(t1) + BIRTH(t1)
91 INT <- function(t) { # vectorized
92 FOO <- function(u) exp(RHO(tt, u)) * death(u)
94 for (i in 1:length(t)) {
96 out[i] <- integrate(FOO, tt, Tmax)$value
102 RHO <- function(t1, t2)
103 death * (t2 - t1) - integrateTrapeze(birth, t1, t2)
105 FOO <- function(u) exp(RHO(t, u)) * death
106 integrateTrapeze(Vectorize(FOO), t, Tmax)
109 RHO <- function(t1, t2)
110 death * (t2 - t1) - integrate(birth, t1, t2)$value
112 FOO <- function(u) exp(RHO(t, u)) * death
113 integrate(Vectorize(FOO), t, Tmax)$value
117 RHO <- function(t1, t2)
118 integrate(death, t1, t2)$value - birth * (t2 - t1)
121 FOO <- function(u) exp(RHO(t, u)) * death(u)
122 integrateTrapeze(FOO, t, Tmax)
126 FOO <- function(u) exp(RHO(t, u)) * death(u)
127 integrate(Vectorize(FOO), t, Tmax)$value
131 RHO <- function(t1, t2) death * (t2 - t1) - BIRTH(t2) + BIRTH(t1)
132 INT <- function(t) { # vectorized
133 FOO <- function(u) exp(RHO(tt, u)) * death
135 for (i in 1:length(t)) {
137 out[i] <- integrate(FOO, tt, Tmax)$value
142 RHO <- function(t1, t2) DEATH(t2) - DEATH(t1) - birth * (t2 - t1)
145 FOO <- function(u) exp(RHO(t, u)) * death(u)
146 integrateTrapeze(FOO, t, Tmax)
150 FOO <- function(u) exp(RHO(t, u)) * death(u)
151 integrate(Vectorize(FOO), t, Tmax)$value
159 function(birth, death, BIRTH = NULL, DEATH = NULL, Tmax, x, case, fast = FALSE)
161 ff <- .getRHOetINT(birth, death, BIRTH, DEATH, case, fast)
164 environment(INT) <- environment() # so that INT() can find Tmax
165 .CDF.birth.death2(RHO, INT, birth, death, BIRTH, DEATH,
170 function(RHO, INT, birth, death, BIRTH, DEATH, Tmax, x, case, fast)
172 Pi <- if (case %in% c(1, 5, 7))
173 function(t) (1/(1 + INT(t)))^2 * 2 * exp(-RHO(0, t)) * birth
175 function(t) (1/(1 + INT(t)))^2 * 2 * exp(-RHO(0, t)) * birth(t)
177 if (!case %in% c(1, 3, 6)) Pi <- Vectorize(Pi)
179 denom <- if (fast) integrateTrapeze(Pi, 0, Tmax) else integrate(Pi, 0, Tmax)$value
183 for (i in 1:n) p[i] <- integrateTrapeze(Pi, 0, x[i])
185 for (i in 1:n) p[i] <- integrate(Pi, 0, x[i])$value
190 .makePhylo <- function(edge, edge.length, i)
193 edge[NODES] <- edge[NODES] + i + 1L
194 edge[!NODES] <- 1:(i + 1L)
196 phy <- list(edge = edge, edge.length = edge.length,
197 tip.label = paste("t", 1:(i + 1), sep = ""), Nnode = i)
198 class(phy) <- "phylo"
203 function(birth, death, Tmax = 50, BIRTH = NULL, DEATH = NULL, eps = 1e-6)
205 case <- .getCase(birth, death, BIRTH, DEATH)
207 rTimeToEvent <- function(t)
209 ## CDF of the times to event (speciation or extinction):
212 Foo <- function(t, x)
213 1 - exp(-(birth + death)*(x - t))
215 Foo <- function(t, x) {
216 if (t == x) return(0)
217 1 - exp(-integrate(function(t) birth(t) + death(t),
221 Foo <- function(t, x) {
222 if (t == x) return(0)
223 1 - exp(-(BIRTH(x) - BIRTH(t) + DEATH(x) - DEATH(t)))
226 Foo <- function(t, x) {
227 if (t == x) return(0)
228 1 - exp(-(integrate(function(t) birth(t), t, x)$value
233 Foo <- function(t, x) {
234 if (t == x) return(0)
235 1 - exp(-(birth*(x - t) +
236 integrate(function(t) death(t), t, x)$value))
240 Foo <- function(t, x) {
241 if (t == x) return(0)
242 1 - exp(-(BIRTH(x) - BIRTH(t) + death*(x - t)))
246 Foo <- function(t, x) {
247 if (t == x) return(0)
248 1 - exp(-(birth*(x - t) + DEATH(x) - DEATH(t) ))
252 ## generate a random time to event by the inverse method:
254 ## in case speciation probability is so low
255 ## that time to speciation is infinite:
256 if (Foo(t, Tmax) < P) return(Tmax + 1)
259 while (inc > eps) { # la précision influe sur le temps de calcul
269 speORext <- function(t) birth/(birth + death)
270 if (case == 2 || case == 3)
271 speORext <- function(t) birth(t)/(birth(t) + death(t))
272 if (case == 4 || case == 6)
273 speORext <- function(t) birth(t)/(birth(t) + death)
274 if (case == 5 || case == 7)
275 speORext <- function(t) birth/(birth + death(t))
277 ## the recursive function implementing algorithm 1
278 foo <- function(node) {
280 X <- rTimeToEvent(t[node])
282 ## is the event a speciation or an extinction?
286 } else Y <- rbinom(1, size = 1, prob = speORext(tmp))
288 edge.length[j] <<- tmp - t[node]
292 ## set internal edge:
293 edge[j, ] <<- c(node, i)
296 ## set terminal edge:
297 edge[j, ] <<- c(node, 0L)
301 edge <- matrix(0L, 1e5, 2)
302 edge.length <- numeric(1e5)
303 j <- 0L; i <- 1; t <- 0
305 .makePhylo(edge[1:j, ], edge.length[1:j], i)
308 drop.fossil <- function(phy, tol = 0)
311 x <- dist.nodes(phy)[n + 1, ][1:n]
312 drop.tip(phy, which(x < max(x) - tol))
316 function(birth, death, Tmax = 50, BIRTH = NULL, DEATH = NULL, eps = 1e-6)
318 case <- .getCase(birth, death, BIRTH, DEATH)
319 ff <- .getRHOetINT(birth, death, BIRTH, DEATH, case, FALSE)
322 ## so that RHO() and INT() can find Tmax:
323 environment(RHO) <- environment(INT) <- environment()
325 rtimetospe <- function(t)
327 ## CDF of the times to speciation:
328 Foo <- if (case %in% c(1, 5, 7))
329 function(t, x) 1 - exp(-birth*(x - t))
331 if (case %in% c(3, 6))
332 function(t, x) 1 - exp(-(BIRTH(x) - BIRTH(t)))
335 if (t == x) return(0)
336 1 - exp(-integrate(birth, t, x)$value)
340 ## generate a random time to speciation by the inverse method:
342 ## in case speciation probability is so low
343 ## that time to speciation is infinite:
344 if (Foo(t, Tmax) < P) return(Tmax + 1)
347 while (inc > eps) { # la précision influe sur le temps de calcul
356 ## the recursive function implementing algorithm 2
357 foo <- function(node, start) {
358 node <- node # make a local copy
360 tau <- start # because tau is changed below
363 while (X < Tmax - tau) {
365 ## does the new lineage survive until Tmax?
366 Y <- rbinom(1, size = 1, prob = 1/(1 + INT(tau)))
370 ## set internal edge:
372 edge[j, ] <<- c(node, i)
373 edge.length[j] <<- tau - t[node]
380 ## set terminal edge:
383 edge[j, 1] <<- node # the 2nd column is already set to 0
384 edge.length[j] <<- Tmax - t[node]
389 edge <- matrix(0L, 1e5, 2)
390 edge.length <- numeric(1e5)
391 j <- 0L; i <- 1L; t <- 0
393 .makePhylo(edge[1:j, ], edge.length[1:j], i)