R/CDF.birth.death.R

   1 ## CDF.birth.death.R (2010-09-27)
   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 integrateTrapeze <- function(FUN, from, to, nint = 10)
  12 ## compute an integral with a simple trapeze method
  13 ## (apparently, Vectorize doesn't give faster calculation)
  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 .getRHOetINT <- function(birth, death, BIRTH = NULL, DEATH = NULL, case, fast)
  45 {
  46     ## build the RHO(), \rho(t), and INT(), I(t), functions
  47     switch (case,
  48         { # case 1:
  49             RHO <- function(t1, t2) (t2 - t1)*(death - birth)
  50             INT <- function(t) {
  51                 rho <- death - birth
  52                 death*(exp(rho*(Tmax - t)) - 1)/rho
  53             }
  54         },{ # case 2:
  55             if (fast) {
  56                 RHO <- function(t1, t2)
  57                     integrateTrapeze(function(t) death(t) - birth(t), t1, t2)
  58                 INT <- function(t) {
  59                     FOO <- function(u) exp(RHO(t, u)) * death(u)
  60                     integrateTrapeze(FOO, t, Tmax)
  61                 }
  62             } else {
  63                 RHO <- function(t1, t2)
  64                     integrate(function(t) death(t) - birth(t), t1, t2)$value
  65                 INT <- function(t) {
  66                     FOO <- function(u) exp(RHO(t, u)) * death(u)
  67                     integrate(Vectorize(FOO), t, Tmax)$value # Vectorize required
  68                 }
  69             }
  70         },{ # case 3:
  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)
  75                 out <- t
  76                 for (i in 1:length(t)) {
  77                     tt <- t[i]
  78                     out[i] <- integrate(FOO, tt, Tmax)$value
  79                 }
  80                 out
  81             }
  82         },{ # case 4:
  83             if (fast) {
  84                 RHO <- function(t1, t2)
  85                     death * (t2 - t1) - integrateTrapeze(birth, t1, t2)
  86                 INT <- function(t) {
  87                     FOO <- function(u) exp(RHO(t, u)) * death
  88                     integrateTrapeze(Vectorize(FOO), t, Tmax)
  89                 }
  90             } else {
  91                 RHO <- function(t1, t2)
  92                     death * (t2 - t1) - integrate(birth, t1, t2)$value
  93                 INT <- function(t) {
  94                     FOO <- function(u) exp(RHO(t, u)) * death
  95                     integrate(Vectorize(FOO), t, Tmax)$value
  96                 }
  97             }
  98         },{ # case 5:
  99             RHO <- function(t1, t2)
 100                 integrate(death, t1, t2)$value - birth * (t2 - t1)
 101             if (fast) {
 102                 INT <- function(t) {
 103                     FOO <- function(u) exp(RHO(t, u)) * death(u)
 104                     integrateTrapeze(FOO, t, Tmax)
 105                 }
 106             } else {
 107                 INT <- function(t) {
 108                     FOO <- function(u) exp(RHO(t, u)) * death(u)
 109                     integrate(Vectorize(FOO), t, Tmax)$value
 110                 }
 111             }
 112         },{ # case 6:
 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
 116                 out <- t
 117                 for (i in 1:length(t)) {
 118                     tt <- t[i]
 119                     out[i] <- integrate(FOO, tt, Tmax)$value
 120                 }
 121                 out
 122             }
 123         },{ # case 7:
 124             RHO <- function(t1, t2) DEATH(t2) - DEATH(t1) - birth * (t2 - t1)
 125             if (fast) {
 126                 INT <- function(t) {
 127                     FOO <- function(u) exp(RHO(t, u)) * death(u)
 128                     integrateTrapeze(FOO, t, Tmax)
 129                 }
 130             } else {
 131                 INT <- function(t) {
 132                     FOO <- function(u) exp(RHO(t, u)) * death(u)
 133                     integrate(Vectorize(FOO), t, Tmax)$value
 134                 }
 135             }
 136         })
 137     list(RHO, INT)
 138 }
 139
 140 CDF.birth.death <-
 141     function(birth, death, BIRTH = NULL, DEATH = NULL, Tmax, x, case, fast = FALSE)
 142 {
 143     ff <- .getRHOetINT(birth, death, BIRTH, DEATH, case, fast)
 144     RHO <- ff[[1]]
 145     INT <- ff[[2]]
 146     environment(INT) <- environment() # so that INT() can find Tmax
 147     .CDF.birth.death2(RHO, INT, birth, death, BIRTH, DEATH,
 148                       Tmax, x, case, fast)
 149 }
 150
 151 .CDF.birth.death2 <-
 152     function(RHO, INT, birth, death, BIRTH, DEATH, Tmax, x, case, fast)
 153 {
 154     Pi <- if (case %in% c(1, 5, 7))
 155         function(t) (1/(1 + INT(t)))^2 * 2 * exp(-RHO(0, t)) * birth
 156     else
 157         function(t) (1/(1 + INT(t)))^2 * 2 * exp(-RHO(0, t)) * birth(t)
 158
 159     if (!case %in% c(1, 3, 6)) Pi <- Vectorize(Pi)
 160
 161     denom <- if (fast) integrateTrapeze(Pi, 0, Tmax) else integrate(Pi, 0, Tmax)$value
 162     n <- length(x)
 163     p <- numeric(n)
 164     if (fast) {
 165         for (i in 1:n) p[i] <- integrateTrapeze(Pi, 0, x[i])
 166     } else {
 167         for (i in 1:n) p[i] <- integrate(Pi, 0, x[i])$value
 168     }
 169     p/denom
 170 }
 171
 172 .makePhylo <- function(edge, edge.length, i)
 173 {
 174     NODES <- edge > 0
 175     edge[NODES] <- edge[NODES] + i + 1L
 176     edge[!NODES] <- 1:(i + 1L)
 177
 178     phy <- list(edge = edge, edge.length = edge.length,
 179                 tip.label = paste("t", 1:(i + 1), sep = ""), Nnode = i)
 180     class(phy) <- "phylo"
 181     phy
 182 }
 183
 184 rlineage <-
 185     function(birth, death, Tmax = 50, BIRTH = NULL, DEATH = NULL, eps = 1e-6)
 186 {
 187     case <- .getCase(birth, death, BIRTH, DEATH)
 188
 189     rTimeToEvent <- function(t)
 190     {
 191         ## CDF of the times to event (speciation or extinction):
 192         switch (case,
 193             { # case 1:
 194                 Foo <- function(t, x)
 195                     1 - exp(-(birth + death)*(x - t))
 196             },{ # case 2:
 197                 Foo <- function(t, x) {
 198                     if (t == x) return(0)
 199                     1 - exp(-integrate(function(t) birth(t) + death(t),
 200                                        t, x)$value)
 201                 }
 202             },{ # case 3:
 203                 Foo <- function(t, x) {
 204                     if (t == x) return(0)
 205                     1 - exp(-(BIRTH(x) - BIRTH(t) + DEATH(x) - DEATH(t)))
 206                 }
 207             },{ # case 4:
 208                 Foo <- function(t, x) {
 209                     if (t == x) return(0)
 210                     1 - exp(-(integrate(function(t) birth(t), t, x)$value
 211                               + death*(x - t)))
 212                 }
 213
 214             },{ # case 5:
 215                 Foo <- function(t, x) {
 216                     if (t == x) return(0)
 217                     1 - exp(-(birth*(x - t) +
 218                               integrate(function(t) death(t), t, x)$value))
 219                 }
 220
 221             },{ # case 6:
 222                 Foo <- function(t, x) {
 223                     if (t == x) return(0)
 224                     1 - exp(-(BIRTH(x) - BIRTH(t) + death*(x - t)))
 225                 }
 226
 227             },{ # case 7:
 228                 Foo <- function(t, x) {
 229                     if (t == x) return(0)
 230                     1 - exp(-(birth*(x - t) + DEATH(x) - DEATH(t) ))
 231                 }
 232             })
 233
 234         ## generate a random time to event by the inverse method:
 235         P <- runif(1)
 236         ## in case speciation probability is so low
 237         ## that time to speciation is infinite:
 238         if (Foo(t, Tmax) < P) return(Tmax + 1)
 239         inc <- 10
 240         x <- t + inc
 241         while (inc > eps) { # la précision influe sur le temps de calcul
 242             if (Foo(t, x) > P) {
 243                 x <- x - inc
 244                 inc <- inc/10
 245             } else x <- x + inc
 246         }
 247         x - t
 248     }
 249
 250     if (case == 1)
 251         speORext <- function(t) birth/(birth + death)
 252     if (case == 2 || case == 3)
 253         speORext <- function(t) birth(t)/(birth(t) + death(t))
 254     if (case == 4 || case == 6)
 255         speORext <- function(t) birth(t)/(birth(t) + death)
 256     if (case == 5 || case == 7)
 257         speORext <- function(t) birth/(birth + death(t))
 258
 259     ## the recursive function implementing algorithm 1
 260     foo <- function(node) {
 261         for (k in 0:1) {
 262             X <- rTimeToEvent(t[node])
 263             tmp <- t[node] + X
 264             ## is the event a speciation or an extinction?
 265             if (tmp >= Tmax) {
 266                 Y <- 0
 267                 tmp <- Tmax
 268             } else Y <- rbinom(1, size = 1, prob = speORext(tmp))
 269             j <<- j + 1L
 270             edge.length[j] <<- tmp - t[node]
 271             if (Y) {
 272                 i <<- i + 1L
 273                 t[i] <<- tmp
 274                 ## set internal edge:
 275                 edge[j, ] <<- c(node, i)
 276                 foo(i)
 277             } else
 278                 ## set terminal edge:
 279                 edge[j, ] <<- c(node, 0L)
 280         }
 281     }
 282
 283     edge <- matrix(0L, 1e5, 2)
 284     edge.length <- numeric(1e5)
 285     j <- 0L; i <- 1; t <- 0
 286     foo(1L)
 287     .makePhylo(edge[1:j, ], edge.length[1:j], i)
 288 }
 289
 290 drop.fossil <- function(phy, tol = 1e-8)
 291 {
 292     n <- Ntip(phy)
 293     x <- dist.nodes(phy)[n + 1, ][1:n]
 294     drop.tip(phy, which(x < max(x) - tol))
 295 }
 296
 297 rbdtree <-
 298     function(birth, death, Tmax = 50, BIRTH = NULL, DEATH = NULL, eps = 1e-6)
 299 {
 300     case <- .getCase(birth, death, BIRTH, DEATH)
 301     ff <- .getRHOetINT(birth, death, BIRTH, DEATH, case, FALSE)
 302     RHO <- ff[[1]]
 303     INT <- ff[[2]]
 304     ## so that RHO() and INT() can find Tmax:
 305     environment(RHO) <- environment(INT) <- environment()
 306
 307     rtimetospe <- function(t)
 308     {
 309         ## CDF of the times to speciation:
 310         Foo <- if (case %in% c(1, 5, 7))
 311             function(t, x) 1 - exp(-birth*(x - t))
 312         else {
 313             if (case %in% c(3, 6))
 314                 function(t, x) 1 - exp(-(BIRTH(x) - BIRTH(t)))
 315             else {
 316                 function(t, x) {
 317                     if (t == x) return(0)
 318                     1 - exp(-integrate(birth, t, x)$value)
 319                 }
 320             }
 321         }
 322         ## generate a random time to speciation by the inverse method:
 323         P <- runif(1)
 324         ## in case speciation probability is so low
 325         ## that time to speciation is infinite:
 326         if (Foo(t, Tmax) < P) return(Tmax + 1)
 327         inc <- 10
 328         x <- t + inc
 329         while (inc > eps) { # la précision influe sur le temps de calcul
 330             if (Foo(t, x) > P) {
 331                 x <- x - inc
 332                 inc <- inc/10
 333             } else x <- x + inc
 334         }
 335         x - t
 336     }
 337
 338     ## the recursive function implementing algorithm 2
 339     foo <- function(node, start) {
 340         node <- node # make a local copy
 341         for (k in 0:1) {
 342             tau <- start # because tau is changed below
 343             NoDesc <- TRUE
 344             X <- rtimetospe(tau)
 345             while (X < Tmax - tau) {
 346                 tau <- tau + X
 347                 ## does the new lineage survive until Tmax?
 348                 Y <- rbinom(1, size = 1, prob = 1/(1 + INT(tau)))
 349                 if (Y) {
 350                     i <<- i + 1L
 351                     t[i] <<- tau
 352                     ## set internal edge:
 353                     j <<- j + 1L
 354                     edge[j, ] <<- c(node, i)
 355                     edge.length[j] <<- tau - t[node]
 356                     foo(i, t[i])
 357                     NoDesc <- FALSE
 358                     break
 359                 }
 360                 X <- rtimetospe(tau)
 361             }
 362             ## set terminal edge:
 363             if (NoDesc) {
 364                 j <<- j + 1L
 365                 edge[j, 1] <<- node # the 2nd column is already set to 0
 366                 edge.length[j] <<- Tmax - t[node]
 367             }
 368         }
 369     }
 370
 371     edge <- matrix(0L, 1e5, 2)
 372     edge.length <- numeric(1e5)
 373     j <- 0L; i <- 1L; t <- 0
 374     foo(1L, 0)
 375     .makePhylo(edge[1:j, ], edge.length[1:j], i)
 376 }
 377
 378 bd.time <- function(phy, birth, death, BIRTH = NULL, DEATH = NULL,
 379                     ip, lower, upper, fast = FALSE, boot = 0, trace = 0)
 380 {
 381     guess.bounds <- if (missing(lower)) TRUE else FALSE
 382     BIG <- 1e10
 383     PrNt <- function(t, T, x)
 384     {
 385         tmp <- exp(-RHO(t, T))
 386         Wt <- tmp * (1 + INT(t))
 387         out <- numeric(length(x))
 388         zero <- x == 0
 389         if (length(zero)) {
 390             out[zero] <- 1 - tmp/Wt
 391             out[!zero] <- (tmp/Wt^2)*(1 - 1/Wt)^(x[!zero] - 1)
 392         } else out[] <- (tmp/Wt^2)*(1 - 1/Wt)^(x - 1)
 393         out
 394     }
 395
 396     case <- .getCase(birth, death, BIRTH, DEATH)
 397
 398     if (is.function(birth)) {
 399         paranam <- names(formals(birth))
 400         if (guess.bounds) {
 401             upper <- rep(BIG, length(paranam))
 402             lower <- -upper
 403         }
 404         formals(birth) <- alist(t=)
 405         environment(birth) <- environment()
 406         if (!is.null(BIRTH)) environment(BIRTH) <- environment()
 407     } else {
 408         paranam <- "birth"
 409         if (guess.bounds) {
 410             upper <- 1
 411             lower <- 0
 412         }
 413     }
 414
 415     if (is.function(death)) {
 416         tmp <- names(formals(death))
 417         np2 <- length(tmp)
 418         if (guess.bounds) {
 419             upper <- c(upper, rep(BIG, np2))
 420             lower <- c(lower, rep(-BIG, np2))
 421         }
 422         paranam <- c(paranam, tmp)
 423         formals(death) <- alist(t=)
 424         environment(death) <- environment()
 425         if (!is.null(DEATH)) environment(DEATH) <- environment()
 426     } else {
 427         paranam <- c(paranam, "death")
 428         if (guess.bounds) {
 429             upper <- c(upper, .1)
 430             lower <- c(lower, 0)
 431         }
 432     }
 433
 434     np <- length(paranam)
 435
 436     ff <- .getRHOetINT(birth, death, BIRTH, DEATH, case = case, fast = fast)
 437     RHO <- ff[[1]]
 438     INT <- ff[[2]]
 439     environment(RHO) <- environment(INT) <- environment()
 440
 441     x <- branching.times(phy)
 442     n <- length(x)
 443     Tmax <- x[1]
 444     x <- Tmax - x # change the time scale so the root is t=0
 445     x <- sort(x)
 446
 447     foo <- function(para) {
 448         for (i in 1:np)
 449             assign(paranam[i], para[i], pos = sys.frame(1))
 450         p <- CDF.birth.death(birth, death, BIRTH, DEATH, Tmax = Tmax,
 451                              x = x, case = case, fast = fast)
 452         ## w is the probability of the observed tree size (= number of tips)
 453         w <- PrNt(0, Tmax, Ntip(phy))
 454         ## p is the expected CDF of branching times
 455         ## ecdf(x)(x) is the observed CDF
 456         sum((1:n/n - p)^2)/w # faster than sum((ecdf(x)(x) - p)^2)/w
 457     }
 458
 459     if (missing(ip)) ip <- (upper - lower)/2
 460
 461     out <- nlminb(ip, foo, control = list(trace = trace, eval.max = 500),
 462                   upper = upper, lower = lower)
 463
 464     names(out$par) <- paranam
 465     names(out)[2] <- "SS"
 466
 467     if (boot) { # nonparametric version
 468         PAR <- matrix(NA, boot, np)
 469         i <- 1L
 470         while (i <= boot) {
 471             cat("\rDoing bootstrap no.", i, "\n")
 472             x <- sort(sample(x, replace = TRUE))
 473             o <- try(nlminb(ip, foo, control = list(trace = 0, eval.max = 500),
 474                             upper = upper, lower = lower))
 475             if (class(o) == "list") {
 476                 PAR[i, ] <- o$par
 477                 i <- i + 1L
 478             }
 479         }
 480         out$boot <- PAR
 481     }
 482     out
 483 }