R/CDF.birth.death.R

   1 ## CDF.birth.death.R (2010-08-09)
   2
   3 ##    Functions to simulate and fit
   4 ##  Time-Dependent Birth-Death Models
   5
   6 ## Copyright 2010 Emmanuel Paradis
   7
   8 ## This file is part of the R-package `ape'.
   9 ## See the file ../COPYING for licensing issues.
  10
  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)
  14 {
  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
  20 }
  21
  22 ## case:
  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)
  30
  31 .getCase <- function(birth, death, BIRTH = NULL, DEATH = NULL)
  32 {
  33     if (is.numeric(birth)) {
  34         if (is.numeric(death)) 1 else {
  35             if (is.null(DEATH)) 5 else 7
  36         }
  37     } else {
  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
  41     }
  42 }
  43
  44 ###.getRHO <- function(birth, death, BIRTH = NULL, DEATH = NULL, case)
  45 ###{
  46 ###    ## build the RHO(), \rho(t), function
  47 ###    switch (case,
  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,
  51 ###            function(t1, t2)
  52 ###                DEATH(t2) - BIRTH(t2) - DEATH(t1) + BIRTH(t1), # case 3
  53 ###            function(t1, t2)
  54 ###                death - integrate(birth, t1, t2)$value, # case 4
  55 ###            function(t1, t2)
  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
  59 ###        )
  60 ###}
  61
  62 .getRHOetINT <- function(birth, death, BIRTH = NULL, DEATH = NULL, case, fast)
  63 {
  64     ## build the RHO(), \rho(t), and INT(), I(t), functions
  65     switch (case,
  66         { # case 1:
  67             RHO <- function(t1, t2) (t2 - t1)*(death - birth)
  68             INT <- function(t) {
  69                 rho <- death - birth
  70                 death*(exp(rho*(Tmax - t)) - 1)/rho
  71             }
  72         },{ # case 2:
  73             if (fast) {
  74                 RHO <- function(t1, t2)
  75                     integrateTrapeze(function(t) death(t) - birth(t), t1, t2)
  76                 INT <- function(t) {
  77                     FOO <- function(u) exp(RHO(t, u)) * death(u)
  78                     integrateTrapeze(FOO, t, Tmax)
  79                 }
  80             } else {
  81                 RHO <- function(t1, t2)
  82                     integrate(function(t) death(t) - birth(t), t1, t2)$value
  83                 INT <- function(t) {
  84                     FOO <- function(u) exp(RHO(t, u)) * death(u)
  85                     integrate(Vectorize(FOO), t, Tmax)$value # Vectorize required
  86                 }
  87             }
  88         },{ # case 3:
  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)
  93                 out <- t
  94                 for (i in 1:length(t)) {
  95                     tt <- t[i]
  96                     out[i] <- integrate(FOO, tt, Tmax)$value
  97                 }
  98                 out
  99             }
 100         },{ # case 4:
 101             if (fast) {
 102                 RHO <- function(t1, t2)
 103                     death * (t2 - t1) - integrateTrapeze(birth, t1, t2)
 104                 INT <- function(t) {
 105                     FOO <- function(u) exp(RHO(t, u)) * death
 106                     integrateTrapeze(Vectorize(FOO), t, Tmax)
 107                 }
 108             } else {
 109                 RHO <- function(t1, t2)
 110                     death * (t2 - t1) - integrate(birth, t1, t2)$value
 111                 INT <- function(t) {
 112                     FOO <- function(u) exp(RHO(t, u)) * death
 113                     integrate(Vectorize(FOO), t, Tmax)$value
 114                 }
 115             }
 116         },{ # case 5:
 117             RHO <- function(t1, t2)
 118                 integrate(death, t1, t2)$value - birth * (t2 - t1)
 119             if (fast) {
 120                 INT <- function(t) {
 121                     FOO <- function(u) exp(RHO(t, u)) * death(u)
 122                     integrateTrapeze(FOO, t, Tmax)
 123                 }
 124             } else {
 125                 INT <- function(t) {
 126                     FOO <- function(u) exp(RHO(t, u)) * death(u)
 127                     integrate(Vectorize(FOO), t, Tmax)$value
 128                 }
 129             }
 130         },{ # case 6:
 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
 134                 out <- t
 135                 for (i in 1:length(t)) {
 136                     tt <- t[i]
 137                     out[i] <- integrate(FOO, tt, Tmax)$value
 138                 }
 139                 out
 140             }
 141         },{ # case 7:
 142             RHO <- function(t1, t2) DEATH(t2) - DEATH(t1) - birth * (t2 - t1)
 143             if (fast) {
 144                 INT <- function(t) {
 145                     FOO <- function(u) exp(RHO(t, u)) * death(u)
 146                     integrateTrapeze(FOO, t, Tmax)
 147                 }
 148             } else {
 149                 INT <- function(t) {
 150                     FOO <- function(u) exp(RHO(t, u)) * death(u)
 151                     integrate(Vectorize(FOO), t, Tmax)$value
 152                 }
 153             }
 154         })
 155     list(RHO, INT)
 156 }
 157
 158 CDF.birth.death <-
 159     function(birth, death, BIRTH = NULL, DEATH = NULL, Tmax, x, case, fast = FALSE)
 160 {
 161     ff <- .getRHOetINT(birth, death, BIRTH, DEATH, case, fast)
 162     RHO <- ff[[1]]
 163     INT <- ff[[2]]
 164     environment(INT) <- environment() # so that INT() can find Tmax
 165     .CDF.birth.death2(RHO, INT, birth, death, BIRTH, DEATH,
 166                       Tmax, x, case, fast)
 167 }
 168
 169 .CDF.birth.death2 <-
 170     function(RHO, INT, birth, death, BIRTH, DEATH, Tmax, x, case, fast)
 171 {
 172     Pi <- if (case %in% c(1, 5, 7))
 173         function(t) (1/(1 + INT(t)))^2 * 2 * exp(-RHO(0, t)) * birth
 174     else
 175         function(t) (1/(1 + INT(t)))^2 * 2 * exp(-RHO(0, t)) * birth(t)
 176
 177     if (!case %in% c(1, 3, 6)) Pi <- Vectorize(Pi)
 178
 179     denom <- if (fast) integrateTrapeze(Pi, 0, Tmax) else integrate(Pi, 0, Tmax)$value
 180     n <- length(x)
 181     p <- numeric(n)
 182     if (fast) {
 183         for (i in 1:n) p[i] <- integrateTrapeze(Pi, 0, x[i])
 184     } else {
 185         for (i in 1:n) p[i] <- integrate(Pi, 0, x[i])$value
 186     }
 187     p/denom
 188 }
 189
 190 .makePhylo <- function(edge, edge.length, i)
 191 {
 192     NODES <- edge > 0
 193     edge[NODES] <- edge[NODES] + i + 1L
 194     edge[!NODES] <- 1:(i + 1L)
 195
 196     phy <- list(edge = edge, edge.length = edge.length,
 197                 tip.label = paste("t", 1:(i + 1), sep = ""), Nnode = i)
 198     class(phy) <- "phylo"
 199     phy
 200 }
 201
 202 rlineage <-
 203     function(birth, death, Tmax = 50, BIRTH = NULL, DEATH = NULL, eps = 1e-6)
 204 {
 205     case <- .getCase(birth, death, BIRTH, DEATH)
 206
 207     rTimeToEvent <- function(t)
 208     {
 209         ## CDF of the times to event (speciation or extinction):
 210         switch (case,
 211             { # case 1:
 212                 Foo <- function(t, x)
 213                     1 - exp(-(birth + death)*(x - t))
 214             },{ # case 2:
 215                 Foo <- function(t, x) {
 216                     if (t == x) return(0)
 217                     1 - exp(-integrate(function(t) birth(t) + death(t),
 218                                        t, x)$value)
 219                 }
 220             },{ # case 3:
 221                 Foo <- function(t, x) {
 222                     if (t == x) return(0)
 223                     1 - exp(-(BIRTH(x) - BIRTH(t) + DEATH(x) - DEATH(t)))
 224                 }
 225             },{ # case 4:
 226                 Foo <- function(t, x) {
 227                     if (t == x) return(0)
 228                     1 - exp(-(integrate(function(t) birth(t), t, x)$value
 229                               + death*(x - t)))
 230                 }
 231
 232             },{ # case 5:
 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))
 237                 }
 238
 239             },{ # case 6:
 240                 Foo <- function(t, x) {
 241                     if (t == x) return(0)
 242                     1 - exp(-(BIRTH(x) - BIRTH(t) + death*(x - t)))
 243                 }
 244
 245             },{ # case 7:
 246                 Foo <- function(t, x) {
 247                     if (t == x) return(0)
 248                     1 - exp(-(birth*(x - t) + DEATH(x) - DEATH(t) ))
 249                 }
 250             })
 251
 252         ## generate a random time to event by the inverse method:
 253         P <- runif(1)
 254         ## in case speciation probability is so low
 255         ## that time to speciation is infinite:
 256         if (Foo(t, Tmax) < P) return(Tmax + 1)
 257         inc <- 10
 258         x <- t + inc
 259         while (inc > eps) { # la précision influe sur le temps de calcul
 260             if (Foo(t, x) > P) {
 261                 x <- x - inc
 262                 inc <- inc/10
 263             } else x <- x + inc
 264         }
 265         x - t
 266     }
 267
 268     if (case == 1)
 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))
 276
 277     ## the recursive function implementing algorithm 1
 278     foo <- function(node) {
 279         for (k in 0:1) {
 280             X <- rTimeToEvent(t[node])
 281             tmp <- t[node] + X
 282             ## is the event a speciation or an extinction?
 283             if (tmp >= Tmax) {
 284                 Y <- 0
 285                 tmp <- Tmax
 286             } else Y <- rbinom(1, size = 1, prob = speORext(tmp))
 287             j <<- j + 1L
 288             edge.length[j] <<- tmp - t[node]
 289             if (Y) {
 290                 i <<- i + 1L
 291                 t[i] <<- tmp
 292                 ## set internal edge:
 293                 edge[j, ] <<- c(node, i)
 294                 foo(i)
 295             } else
 296                 ## set terminal edge:
 297                 edge[j, ] <<- c(node, 0L)
 298         }
 299     }
 300
 301     edge <- matrix(0L, 1e5, 2)
 302     edge.length <- numeric(1e5)
 303     j <- 0L; i <- 1; t <- 0
 304     foo(1L)
 305     .makePhylo(edge[1:j, ], edge.length[1:j], i)
 306 }
 307
 308 drop.fossil <- function(phy, tol = 0)
 309 {
 310     n <- Ntip(phy)
 311     x <- dist.nodes(phy)[n + 1, ][1:n]
 312     drop.tip(phy, which(x < max(x) - tol))
 313 }
 314
 315 rbdtree <-
 316     function(birth, death, Tmax = 50, BIRTH = NULL, DEATH = NULL, eps = 1e-6)
 317 {
 318     case <- .getCase(birth, death, BIRTH, DEATH)
 319     ff <- .getRHOetINT(birth, death, BIRTH, DEATH, case, FALSE)
 320     RHO <- ff[[1]]
 321     INT <- ff[[2]]
 322     ## so that RHO() and INT() can find Tmax:
 323     environment(RHO) <- environment(INT) <- environment()
 324
 325     rtimetospe <- function(t)
 326     {
 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))
 330         else {
 331             if (case %in% c(3, 6))
 332                 function(t, x) 1 - exp(-(BIRTH(x) - BIRTH(t)))
 333             else {
 334                 function(t, x) {
 335                     if (t == x) return(0)
 336                     1 - exp(-integrate(birth, t, x)$value)
 337                 }
 338             }
 339         }
 340         ## generate a random time to speciation by the inverse method:
 341         P <- runif(1)
 342         ## in case speciation probability is so low
 343         ## that time to speciation is infinite:
 344         if (Foo(t, Tmax) < P) return(Tmax + 1)
 345         inc <- 10
 346         x <- t + inc
 347         while (inc > eps) { # la précision influe sur le temps de calcul
 348             if (Foo(t, x) > P) {
 349                 x <- x - inc
 350                 inc <- inc/10
 351             } else x <- x + inc
 352         }
 353         x - t
 354     }
 355
 356     ## the recursive function implementing algorithm 2
 357     foo <- function(node, start) {
 358         node <- node # make a local copy
 359         for (k in 0:1) {
 360             tau <- start # because tau is changed below
 361             NoDesc <- TRUE
 362             X <- rtimetospe(tau)
 363             while (X < Tmax - tau) {
 364                 tau <- tau + X
 365                 ## does the new lineage survive until Tmax?
 366                 Y <- rbinom(1, size = 1, prob = 1/(1 + INT(tau)))
 367                 if (Y) {
 368                     i <<- i + 1L
 369                     t[i] <<- tau
 370                     ## set internal edge:
 371                     j <<- j + 1L
 372                     edge[j, ] <<- c(node, i)
 373                     edge.length[j] <<- tau - t[node]
 374                     foo(i, t[i])
 375                     NoDesc <- FALSE
 376                     break
 377                 }
 378                 X <- rtimetospe(tau)
 379             }
 380             ## set terminal edge:
 381             if (NoDesc) {
 382                 j <<- j + 1L
 383                 edge[j, 1] <<- node # the 2nd column is already set to 0
 384                 edge.length[j] <<- Tmax - t[node]
 385             }
 386         }
 387     }
 388
 389     edge <- matrix(0L, 1e5, 2)
 390     edge.length <- numeric(1e5)
 391     j <- 0L; i <- 1L; t <- 0
 392     foo(1L, 0)
 393     .makePhylo(edge[1:j, ], edge.length[1:j], i)
 394 }