version 2.1

DESCRIPTION

@@ -1,16 +1,12 @@
-Package: laser
-Version: 1.0
-Date: 2006-06-06
-Title: Likelihood Analysis of Speciation/Extinction Rates from Phylogenies
-Author: Dan Rabosky
-Maintainer: Dan Rabosky <DLR32@cornell.edu>
-Depends: R(>= 2.0), ape
-Suggests: apTreeshape
-Description: laser implements maximum likelihood methods based on the birth-death
-	process to test whether diversification rates have changed over time.  The package 
-	permits batch processing of phylogenies to generate null distributions of test statistics
-	and posterior distributions of parameter estimates.  Additional functions for manipulating 
-	branching times from molecular phylogenies and for simulating branching times under constant-rate models
-	of diversification are provided.
-License: GPL 2 or later
-Packaged: Sun Jul  9 11:29:14 2006; DRabosky
+Package: laser
+Version: 2.1
+Date: 2007-5-20
+Title: Likelihood Analysis of Speciation/Extinction Rates from Phylogenies
+Author: Dan Rabosky
+Maintainer: Dan Rabosky <DLR32@cornell.edu>
+Depends: R(>= 2.6), ape, geiger
+Suggests: 
+Description: laser implements maximum likelihood methods based on the birth-death
+	process to test whether diversification rates have changed over time and whether rates vary among lineages.
+License: GPL 2 or later
+Packaged: Sun May 25 21:16:25 2008; hornik

NAMESPACE

@@ -0,0 +1 @@
+exportPattern("^[[:alpha:]]+")

R/DDL.R

@@ -0,0 +1,21 @@
+`DDL` <-
+function(x)
+{
+  if (!is.numeric(x)) stop("object x not of class 'numeric'")
+  #calculates likelihoods under DD logistic model
+  x <- rev(sort(x))
+  N <- length(x)+1
+  b <- sort(x)
+  z <- rev(c(b[1], diff(b)))
+  ddfunc <- function(r, k)
+  {
+    -(sum(log(2:(N-1))) + (N-2)*log(r) + sum(log(1-((2:(N-1))/k)))
+    - sum((2:N)*r*z) + sum(z*r*(2:N)^2)/k)
+  }
+   res <- suppressWarnings(nlm(function(p) ddfunc(p[1], p[2]), c(.5, N*2), hessian = TRUE))
+   #may want to recode this to use 'optim' rather than 'nlm'
+   aic <- 2*res$minimum + 4
+   summ <- structure(list(LH = -res$minimum, aic = aic, r1 = res$estimate[1], kparam = res$estimate[2]))
+   return(summ)
+}
+

R/DDX.R

@@ -0,0 +1,24 @@
+`DDX` <-
+function(x)
+{
+  if (!is.numeric(x)) stop("object x not of class 'numeric'")
+  x <- rev(sort(x))
+  N <- length(x)+1
+  b <- sort(x)
+  z <- rev(c(b[1], diff(b)))
+  #initial parameterization for nlm fx: estimates under yule model
+  s1 <- yuleint2(x, x[1], 0)
+
+  ddfunc <- function(r, v)
+  {
+    -(sum(log(2:(N-1))) + (N-2)*log(r) - sum(v*log(2:(N-1)))
+    - sum(((2:N)^(1-v))*z *r) )
+  }
+  res <- suppressWarnings(nlm(function(p) ddfunc(p[1], p[2]), c(s1$smax, 0), hessian = TRUE))
+  #may want to recode this to use 'optim' rather than 'nlm'
+  aic <- 2*res$minimum + 4
+  summ <- structure(list(LH = -res$minimum, aic = aic, r1 = res$estimate[1], xparam = res$estimate[2]));
+  return(summ);
+
+}
+

R/IDDL.R

@@ -0,0 +1,21 @@
+`IDDL` <-
+function(x)
+{
+  #calculates likelihoods under DD logistic model
+  N <- length(x)+1
+  b <- sort(x)
+  z <- rev(c(b[1], diff(b)))
+  res <-list()
+  #negative log-LH function to be minimized
+  ddfunc <- function(r, k)
+  {
+    -(sum(log(2:(N-1))) + (N-2)*log(r) + sum(log(1-((2:(N-1))/k)))
+    - sum((2:N)*r*z) + sum(z*r*(2:N)^2)/k)
+  }
+   temp <- nlm(function(p) ddfunc(p[1], p[2]), c(.5, N*2), hessian = TRUE)
+   res$LH <- -temp$minimum
+   res$r1 <- temp$estimate[1]
+   res$k <- temp$estimate[2]
+   return(res)
+}
+

R/IDDX.R

@@ -0,0 +1,26 @@
+`IDDX` <-
+function(x)
+{
+
+  N <- length(x)+1
+  b <- sort(x)
+  z <- rev(c(b[1], diff(b)))
+  #initial parameterization for nlm fx: estimates under yule model
+  s1 <- IpureBirth(x)
+  #print(s1$r1)
+  res <- list()
+  temp <- list()
+  ddfunc <- function(r, v)
+  {
+    -(sum(log(2:(N-1))) + (N-2)*log(r) - sum(v*log(2:(N-1)))
+    - sum(((2:N)^(1-v))*z *r) )
+  }
+  temp <- nlm(function(p) ddfunc(p[1], p[2]), c(s1$r1, 0), hessian = TRUE)
+
+  res$LH <- -temp$minimum
+  res$r1 <- temp$estimate[1]
+  res$xp <- temp$estimate[2]
+
+   return(res)
+}
+

R/Ibd.R

@@ -0,0 +1,43 @@
+`Ibd` <-
+function(x) #new 'optim' version, 6.4
+{
+  N <- length(x)+1
+  b <- sort(x)
+  z <- rev(c(b[1], diff(b)))
+  x <- c(0, x)
+  res <- list()
+  ai <- c(.1, .9)
+  mlbd <- function(v)
+  {
+    r <- v[1]
+    a <- v[2]
+    -( sum(log(1:(N-1))) + ((N-2)*log(r))
+        + (r*sum(x[3:N]))
+        +(N*log(1-a)) - 2 * sum(log(exp(r * x[2:N])-a)))
+  }
+
+  for (k in 1:length(ai))
+  {
+    temp <- suppressWarnings(optim(c(.2, ai[k]), mlbd))
+    if (temp$par[2] <= 0)
+    {
+      temp <- IpureBirth(x[2:length(x)])
+      if (k == 1 || (k > 1 && temp$LH > res$LH))
+      {
+        res$LH <- temp$LH
+        res$r1 <- temp$r1
+        res$a <- 0
+      }
+    }
+    else if (k == 1 || (k > 1 && res$LH < -temp$value))
+    {
+        res$LH <- -temp$value
+        res$r1 <- temp$par[1]
+        res$a <- temp$par[2]
+    }
+  }
+
+  return(res)
+
+}
+

R/IpureBirth.R

@@ -0,0 +1,10 @@
+`IpureBirth` <-
+function(x)
+{
+  res <- list()
+  temp <- yuleint2(x, x[1], 0)
+  res$LH <- temp$LH
+  res$r1 <- temp$smax
+  return(res)
+}
+

R/Irvbd.R

@@ -0,0 +1,73 @@
+`Irvbd` <-
+function(x, ai = c(.1, .5, .9), ints = NULL)
+{
+    res <- data.frame()
+    temp <- data.frame()
+
+    N <- length(x)+1
+    Nvec <- 2:N
+     x <- c(0, x)
+
+    rlist <- list()
+
+    if (is.null(ints))
+      stvec <- x[4:(length(x)-2)]
+    else
+      stvec <- seq(x[4], x[length(x)-2], length.out = ints)
+
+   LHra  <- function(v)   #v[1], v[2], v[3] v[4]
+   {
+    r1 <- v[1]
+    r2 <- v[2]
+    a <- v[3]
+   -(sum(log(1:(N-1))) + (i-2)*log(r1) + (N - i)*log(r2)
+    + sum((x[3:i]-st)*r1 + st*r2) + sum(x[(i+1):N]*r2)
+    + N*log(1-a)
+    - 2 * sum(log((exp((x[2:i]-st)*r1 + (st*r2)) - a)))
+    - 2 * sum(log((exp(x[(i+1):N]*r2) - a))))
+   }
+
+  for (j in 1:length(stvec))
+  {
+    for (z in 2:(length(x)-1))
+    {
+      if (x[z] >= stvec[j] && x[z+1] < stvec[j])
+      {
+        i = z
+      }
+    }
+    for (k in 1:length(ai))
+    {
+      st <- stvec[j]
+
+      temp <- suppressWarnings(optim(c(.3, .3, ai[k]), LHra))
+      if (temp$par[3] < 0)
+      {
+        t1 <- yuleint2(x[2:length(x)], x[2], st)
+        t2 <- yuleint2(x[2:length(x)], st, 0)
+        temp$par[1] <- t1$smax
+        temp$par[2] <- t2$smax
+        temp$par[3] <- 0
+        temp$value <- -(t1$LH + t2$LH)
+      }
+      if (k == 1)
+        res <- temp
+      else if ((k > 1) && (temp$value < res$value))
+        res <- temp
+    }
+   rlist$r1[j] <- res$par[1]
+   rlist$r2[j] <- res$par[2]
+   rlist$LH[j] <- -res$value
+   rlist$st[j] <- st
+   rlist$a[j] <- res$par[3]
+  }
+  res2 <- list()
+  res2$LH <- max(rlist$LH)
+  res2$r1 <- rlist$r1[rlist$LH == max(rlist$LH)]
+  res2$r2 <- rlist$r2[rlist$LH == max(rlist$LH)]
+  res2$a <- rlist$a[rlist$LH == max(rlist$LH)]
+  res2$st <- rlist$st[rlist$LH == max(rlist$LH)]
+  return(res2)
+
+}
+

R/Iyule2rate.R

@@ -0,0 +1,27 @@
+`Iyule2rate` <-
+function(x, ints = NULL)
+{
+ res <- list()
+ temp <- list()
+
+ if (is.null(ints))
+      stvec <- x[3:(length(x)-2)]
+ else
+      stvec <- seq(x[3], x[length(x)-2], length.out = ints)
+
+ for (i in 1:length(stvec))
+ {
+    v1 <- yuleint2(x, x[1], stvec[i])
+    v2 <- yuleint2(x, stvec[i], 0)
+    temp$LH[i] <- v1$LH + v2$LH
+    temp$st[i] <- stvec[i]
+    temp$r1[i] <- v1$smax
+    temp$r2[i] <- v2$smax
+ }
+    res$LH <- max(temp$LH)
+    res$st <- temp$st[temp$LH == max(temp$LH)]
+    res$r1 <- temp$r1[temp$LH == max(temp$LH)]
+    res$r2 <- temp$r2[temp$LH == max(temp$LH)]
+    return (res)
+}
+

R/Iyule3rate.R

@@ -0,0 +1,47 @@
+`Iyule3rate` <-
+function(x, ints = NULL)
+{
+  x <- rev(sort(x))
+  N <- length(x)+1
+
+  if (is.null(ints))
+  {
+    tv <- x[2:(N-2)]
+    tv <- unique(tv)
+    stvec <- combinations(length(tv), 2, tv)
+  }
+  else
+  {
+    inc <- (x[2] - x[length(x)])/ints
+    iv <- seq(x[2], length.out = ints, by = -inc)
+    stvec <- combinations(length(iv), 2, iv)
+   }
+  for (i in 1:(length(stvec[,1])))
+  {
+    stvec[i,] <- rev(sort(stvec[i,]))
+  }
+
+  res <- list()
+  for (i in 1:(length(stvec[,1])))
+  {
+     v1 <- 0
+     v2 <- 0
+     v3 <- 0
+      v1 <- yuleint2(x, x[1], stvec[i, 1])
+      v2 <- yuleint2(x, stvec[i, 1], stvec[i, 2])
+      v3 <- yuleint2(x, stvec[i,2], 0)
+      res$LH[i] <- v1$LH + v2$LH + v3$LH
+      res$r1[i] <- v1$smax
+      res$r2[i] <- v2$smax
+      res$r3[i] <- v3$smax
+      res$st[i] <- stvec[i, 1]
+      res$st2[i] <- stvec[i, 2]
+
+  }
+
+     res <- as.data.frame(res)
+     res <- na.omit(res)
+     summ <- res[res$LH == max(res$LH), ]
+     return(summ)
+}
+