R/yule.time.R

   1 ## yule.time.R (2009-02-20)
   2
   3 ##    Fits the Time-Dependent Yule Model
   4
   5 ## Copyright 2009 Emmanuel Paradis
   6
   7 ## This file is part of the R-package `ape'.
   8 ## See the file ../COPYING for licensing issues.
   9
  10 yule.time <- function(phy, birth, BIRTH = NULL, root.time = 0,
  11                       opti = "nlm", start = 0.01)
  12 {
  13     opti <- pmatch(opti, c("nlm", "nlminb", "optim"))
  14     if (is.na(opti)) stop("ambiguous argument 'opti'")
  15     LAMBDA <- function() x
  16     body(LAMBDA) <- body(birth)
  17     formals(LAMBDA) <- alist(t=)
  18     BT <- branching.times(phy)
  19     T <- BT[1]
  20     x <- BT[1] - BT + root.time
  21     m <- phy$Nnode
  22
  23     paranam <- c(names(formals(birth)))
  24     np <- length(paranam)
  25     start <- rep(start, length.out = np)
  26
  27     ## Foo is always vectorized
  28     if (is.null(BIRTH)) {
  29         Foo <- function(x) {
  30             n <- length(x)
  31             res <- numeric(n)
  32             for (i in 1:n)
  33                 res[i] <- integrate(LAMBDA, x[i], T)$value
  34             res
  35         }
  36     } else {
  37         environment(BIRTH) <- environment()
  38         Foo <- function(x) BIRTH(T) - BIRTH(x)
  39     }
  40
  41     half.dev <- function(p) {
  42         for (i in 1:np)
  43             assign(paranam[i], p[i], pos = sys.frame(1))
  44         root.term <-
  45             if (is.null(BIRTH)) integrate(LAMBDA, x[1], T)$value
  46             else BIRTH(T) - BIRTH(x[1])
  47         sum(Foo(x)) + root.term - sum(log(LAMBDA(x[2:m])))
  48     }
  49
  50     switch(opti,
  51         {
  52             out <- nlm(half.dev, start, hessian = TRUE)
  53             est <- out$estimate
  54             se <- sqrt(diag(solve(out$hessian)))
  55             loglik <- lfactorial(m) - out$minimum
  56         },{
  57             out <- nlminb(start, half.dev)
  58             est <- out$par
  59             se <- NULL
  60             loglik <- lfactorial(m) - out$objective
  61         },{
  62             out <- optim(start, half.dev, hessian = TRUE,
  63                          control = list(maxit = 1000), method = "BFGS")
  64             est <- out$par
  65             se <- sqrt(diag(solve(out$hessian)))
  66             loglik <- lfactorial(m) - out$value
  67         })
  68     names(est) <- paranam
  69     if (!is.null(se)) names(se) <- paranam
  70     structure(list(estimate = est, se = se, loglik = loglik), class = "yule")
  71 }