R/CDF.birth.death.R

   1 ## CDF.birth.death.R (2012-09-14)
   2
   3 ##    Functions to Simulate and Fit
   4 ##  Time-Dependent Birth-Death Models
   5
   6 ## Copyright 2010-2012 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 <-
 162         if (fast) integrateTrapeze(Pi, 0, Tmax)
 163         else integrate(Pi, 0, Tmax)$value
 164     n <- length(x)
 165     p <- numeric(n)
 166     if (fast) {
 167         for (i in 1:n) p[i] <- integrateTrapeze(Pi, 0, x[i])
 168     } else {
 169         for (i in 1:n) p[i] <- integrate(Pi, 0, x[i])$value
 170     }
 171     p/denom
 172 }
 173
 174 .makePhylo <- function(edge, edge.length, i)
 175 {
 176     NODES <- edge > 0
 177     edge[NODES] <- edge[NODES] + i + 1L
 178     edge[!NODES] <- 1:(i + 1L)
 179
 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"
 184     phy
 185 }
 186
 187 rlineage <-
 188     function(birth, death, Tmax = 50, BIRTH = NULL, DEATH = NULL, eps = 1e-6)
 189 {
 190     case <- .getCase(birth, death, BIRTH, DEATH)
 191
 192     rTimeToEvent <- function(t)
 193     {
 194         ## CDF of the times to event (speciation or extinction):
 195         switch (case,
 196             { # case 1:
 197                 Foo <- function(t, x)
 198                     1 - exp(-(birth + death)*(x - t))
 199             },{ # case 2:
 200                 Foo <- function(t, x) {
 201                     if (t == x) return(0)
 202                     1 - exp(-integrate(function(t) birth(t) + death(t),
 203                                        t, x)$value)
 204                 }
 205             },{ # case 3:
 206                 Foo <- function(t, x) {
 207                     if (t == x) return(0)
 208                     1 - exp(-(BIRTH(x) - BIRTH(t) + DEATH(x) - DEATH(t)))
 209                 }
 210             },{ # case 4:
 211                 Foo <- function(t, x) {
 212                     if (t == x) return(0)
 213                     1 - exp(-(integrate(function(t) birth(t), t, x)$value
 214                               + death*(x - t)))
 215                 }
 216
 217             },{ # case 5:
 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))
 222                 }
 223
 224             },{ # case 6:
 225                 Foo <- function(t, x) {
 226                     if (t == x) return(0)
 227                     1 - exp(-(BIRTH(x) - BIRTH(t) + death*(x - t)))
 228                 }
 229
 230             },{ # case 7:
 231                 Foo <- function(t, x) {
 232                     if (t == x) return(0)
 233                     1 - exp(-(birth*(x - t) + DEATH(x) - DEATH(t) ))
 234                 }
 235             })
 236
 237         ## generate a random time to event by the inverse method:
 238         P <- runif(1)
 239         ## in case speciation probability is so low
 240         ## that time to speciation is infinite:
 241         if (Foo(t, Tmax) < P) return(Tmax + 1)
 242         inc <- 10
 243         x <- t + inc
 244         while (inc > eps) { # la précision influe sur le temps de calcul
 245             if (Foo(t, x) > P) {
 246                 x <- x - inc
 247                 inc <- inc/10
 248             } else x <- x + inc
 249         }
 250         x - t
 251     }
 252
 253     if (case == 1)
 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))
 261
 262     ## the recursive function implementing algorithm 1
 263     foo <- function(node) {
 264         for (k in 0:1) {
 265             X <- rTimeToEvent(t[node])
 266             tmp <- t[node] + X
 267             ## is the event a speciation or an extinction?
 268             if (tmp >= Tmax) {
 269                 Y <- 0
 270                 tmp <- Tmax
 271             } else Y <- rbinom(1, size = 1, prob = speORext(tmp))
 272             j <<- j + 1L
 273             edge.length[j] <<- tmp - t[node]
 274             if (Y) {
 275                 i <<- i + 1L
 276                 t[i] <<- tmp
 277                 ## set internal edge:
 278                 edge[j, ] <<- c(node, i)
 279                 foo(i)
 280             } else
 281                 ## set terminal edge:
 282                 edge[j, ] <<- c(node, 0L)
 283         }
 284     }
 285
 286     edge <- matrix(0L, 1e5, 2)
 287     edge.length <- numeric(1e5)
 288     j <- 0L; i <- 1; t <- 0
 289     foo(1L)
 290     .makePhylo(edge[1:j, ], edge.length[1:j], i)
 291 }
 292
 293 drop.fossil <- function(phy, tol = 1e-8)
 294 {
 295     n <- Ntip(phy)
 296     x <- dist.nodes(phy)[n + 1, ][1:n]
 297     drop.tip(phy, which(x < max(x) - tol))
 298 }
 299
 300 rbdtree <-
 301     function(birth, death, Tmax = 50, BIRTH = NULL, DEATH = NULL, eps = 1e-6)
 302 {
 303     case <- .getCase(birth, death, BIRTH, DEATH)
 304     ff <- .getRHOetINT(birth, death, BIRTH, DEATH, case, FALSE)
 305     RHO <- ff[[1]]
 306     INT <- ff[[2]]
 307     ## so that RHO() and INT() can find Tmax:
 308     environment(RHO) <- environment(INT) <- environment()
 309
 310     rtimetospe <- function(t)
 311     {
 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))
 315         else {
 316             if (case %in% c(3, 6))
 317                 function(t, x) 1 - exp(-(BIRTH(x) - BIRTH(t)))
 318             else {
 319                 function(t, x) {
 320                     if (t == x) return(0)
 321                     1 - exp(-integrate(birth, t, x)$value)
 322                 }
 323             }
 324         }
 325         ## generate a random time to speciation by the inverse method:
 326         P <- runif(1)
 327         ## in case speciation probability is so low
 328         ## that time to speciation is infinite:
 329         if (Foo(t, Tmax) < P) return(Tmax + 1)
 330         inc <- 10
 331         x <- t + inc
 332         while (inc > eps) { # la précision influe sur le temps de calcul
 333             if (Foo(t, x) > P) {
 334                 x <- x - inc
 335                 inc <- inc/10
 336             } else x <- x + inc
 337         }
 338         x - t
 339     }
 340
 341     ## the recursive function implementing algorithm 2
 342     foo <- function(node, start) {
 343         node <- node # make a local copy
 344         for (k in 0:1) {
 345             tau <- start # because tau is changed below
 346             NoDesc <- TRUE
 347             X <- rtimetospe(tau)
 348             while (X < Tmax - tau) {
 349                 tau <- tau + X
 350                 ## does the new lineage survive until Tmax?
 351                 Y <- rbinom(1, size = 1, prob = 1/(1 + INT(tau)))
 352                 if (Y) {
 353                     i <<- i + 1L
 354                     t[i] <<- tau
 355                     ## set internal edge:
 356                     j <<- j + 1L
 357                     edge[j, ] <<- c(node, i)
 358                     edge.length[j] <<- tau - t[node]
 359                     foo(i, t[i])
 360                     NoDesc <- FALSE
 361                     break
 362                 }
 363                 X <- rtimetospe(tau)
 364             }
 365             ## set terminal edge:
 366             if (NoDesc) {
 367                 j <<- j + 1L
 368                 edge[j, 1] <<- node # the 2nd column is already set to 0
 369                 edge.length[j] <<- Tmax - t[node]
 370             }
 371         }
 372     }
 373
 374     edge <- matrix(0L, 1e5, 2)
 375     edge.length <- numeric(1e5)
 376     j <- 0L; i <- 1L; t <- 0
 377     foo(1L, 0)
 378     .makePhylo(edge[1:j, ], edge.length[1:j], i)
 379 }
 380
 381 bd.time <- function(phy, birth, death, BIRTH = NULL, DEATH = NULL,
 382                     ip, lower, upper, fast = FALSE, boot = 0, trace = 0)
 383 {
 384     guess.bounds <- if (missing(lower)) TRUE else FALSE
 385     BIG <- 1e10
 386     PrNt <- function(t, T, x)
 387     {
 388         tmp <- exp(-RHO(t, T))
 389         Wt <- tmp * (1 + INT(t))
 390         out <- numeric(length(x))
 391         zero <- x == 0
 392         if (length(zero)) {
 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)
 396         out
 397     }
 398
 399     case <- .getCase(birth, death, BIRTH, DEATH)
 400
 401     if (is.function(birth)) {
 402         paranam <- names(formals(birth))
 403         if (guess.bounds) {
 404             upper <- rep(BIG, length(paranam))
 405             lower <- -upper
 406         }
 407         formals(birth) <- alist(t=)
 408         environment(birth) <- environment()
 409         if (!is.null(BIRTH)) environment(BIRTH) <- environment()
 410     } else {
 411         paranam <- "birth"
 412         if (guess.bounds) {
 413             upper <- 1
 414             lower <- 0
 415         }
 416     }
 417
 418     if (is.function(death)) {
 419         tmp <- names(formals(death))
 420         np2 <- length(tmp)
 421         if (guess.bounds) {
 422             upper <- c(upper, rep(BIG, np2))
 423             lower <- c(lower, rep(-BIG, np2))
 424         }
 425         paranam <- c(paranam, tmp)
 426         formals(death) <- alist(t=)
 427         environment(death) <- environment()
 428         if (!is.null(DEATH)) environment(DEATH) <- environment()
 429     } else {
 430         paranam <- c(paranam, "death")
 431         if (guess.bounds) {
 432             upper <- c(upper, .1)
 433             lower <- c(lower, 0)
 434         }
 435     }
 436
 437     np <- length(paranam)
 438
 439     ff <- .getRHOetINT(birth, death, BIRTH, DEATH, case = case, fast = fast)
 440     RHO <- ff[[1]]
 441     INT <- ff[[2]]
 442     environment(RHO) <- environment(INT) <- environment()
 443
 444     x <- branching.times(phy)
 445     n <- length(x)
 446     Tmax <- x[1]
 447     x <- Tmax - x # change the time scale so the root is t=0
 448     x <- sort(x)
 449
 450     foo <- function(para) {
 451         for (i in 1:np)
 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
 460     }
 461
 462     if (missing(ip)) ip <- (upper - lower)/2
 463
 464     out <- nlminb(ip, foo, control = list(trace = trace, eval.max = 500),
 465                   upper = upper, lower = lower)
 466
 467     names(out$par) <- paranam
 468     names(out)[2] <- "SS"
 469
 470     if (boot) { # nonparametric version
 471         PAR <- matrix(NA, boot, np)
 472         i <- 1L
 473         while (i <= boot) {
 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") {
 479                 PAR[i, ] <- o$par
 480                 i <- i + 1L
 481             }
 482         }
 483         out$boot <- PAR
 484     }
 485     out
 486 }