R/skyline.R

   1 ## skyline.R (2002-09-12)
   2
   3 ##   Methods to construct skyline objects (data underlying skyline plot)
   4
   5 ## Copyright 2002 Korbinian Strimmer
   6
   7 ## This file is part of the R-package `ape'.
   8 ## See the file ../COPYING for licensing issues.
   9
  10 skyline <- function(x, ...) UseMethod("skyline")
  11
  12 # input: phylogenetic tree
  13 skyline.phylo <- function(x, ...)
  14 {
  15   if (class(x) != "phylo")
  16     stop("object \"x\" is not of class \"phylo\"")
  17
  18   skyline(coalescent.intervals(x), ...)
  19 }
  20
  21 # input: coalescent intervals and epsilon
  22 skyline.coalescentIntervals <- function(x, epsilon=0, ...)
  23 {
  24   if (class(x) != "coalescentIntervals")
  25     stop("object \"x\" is not of class \"coalescentIntervals\"")
  26
  27   if (epsilon < 0)
  28   {
  29     eps <- find.skyline.epsilon(x, ...)
  30   }
  31   else
  32     eps <- epsilon
  33
  34   skyline(collapsed.intervals(x, epsilon=eps), ...)
  35 }
  36
  37
  38 # input: collapsed intervals
  39 skyline.collapsedIntervals <- function(x, old.style=FALSE, ...)
  40 {
  41   if (class(x) != "collapsedIntervals")
  42     stop("object \"x\" is not of class \"collapsedIntervals\"")
  43
  44   link <- x$collapsed.interval
  45   params <- x$collapsed.interval.count
  46   l <- x$lineages
  47   w <- x$interval.length
  48
  49   b <- choose(l,2) # binomial coefficients
  50
  51   sg <- rep(0,params)   # sizes of collapsed intervals
  52   cg <- rep(0,params)   # coalescent events in interval
  53
  54   if(old.style)
  55     ng <- rep(0,params) # lineages at beginning of an in interval
  56   else
  57   {
  58     ng <- rep(0,params) # sum of classic skp estimates in an interval
  59     m.classic <- w*b
  60   }
  61
  62   for (i in 1:params)
  63   {
  64     group <- link==i
  65     sgr <- w[group]
  66     sg[[i]] <- sum(sgr)
  67     cg[[i]] <- length(sgr)
  68
  69     if(old.style)
  70       ng[[i]] <- l[group][[1]]
  71     else
  72       ng[[i]] <- sum(m.classic[group])
  73   }
  74
  75   # generalized skp estimate
  76   t <- cumsum(sg)
  77   if (old.style)
  78     m <- sg*(ng*(ng-cg)/(2.0*cg) )
  79   else
  80     m <- ng/cg
  81
  82   # log-likelihood
  83   logL <- sum(log(b/m[link]) - b/m[link]*w)
  84
  85   # AICc corrected log-likelihood
  86   K <- x$collapsed.interval.count
  87   S <- x$interval.count
  88   if (S-K > 1)
  89     logL.AICc <- logL - K- K*(K+1)/(S-K-1)
  90   else
  91     logL.AICc <- NA
  92
  93   obj <- list(
  94     time=t,
  95     interval.length=sg,
  96     population.size=m,
  97     parameter.count=length(t),
  98     epsilon = x$epsilon,
  99     logL = logL,
 100     logL.AICc = logL.AICc
 101   )
 102   class(obj) <- "skyline"
 103   return(obj)
 104 }
 105
 106 # grid search for finding optimal epsilon parameter
 107 find.skyline.epsilon <- function(ci, GRID=1000, MINEPS=1e-6, ...)
 108 {
 109   # Why MINEPS?
 110   # Because most "clock-like" trees are not properly
 111   # clock-like for a variety of reasons, i.e. the heights
 112   # of the tips are not exactly zero.
 113
 114   cat("Searching for the optimal epsilon... ")
 115
 116   # a grid search is a naive way but still effective of doing this ...
 117
 118   size <- ci$interval.count
 119   besteps <- ci$total.depth
 120   eps <- besteps
 121
 122   cli <- collapsed.intervals(ci,eps)
 123   skpk <- skyline(cli, ...)
 124   bestaicc <- skpk$ logL.AICc
 125   params <- skpk$parameter.count
 126
 127   delta <- besteps/GRID
 128
 129   eps <- eps-delta
 130   while(eps > MINEPS)
 131   {
 132     cli <- collapsed.intervals(ci,eps)
 133     skpk <- skyline(cli, ...)
 134     aicc <- skpk$ logL.AICc
 135     params <- skpk$parameter.count
 136
 137     if (aicc > bestaicc && params < size-1)
 138     {
 139       besteps <- eps
 140       bestaicc <- aicc
 141     }
 142     eps <- eps-delta
 143   }
 144
 145    cat("epsilon =", besteps, "\n")
 146
 147   besteps
 148 }
 149