1 ## skyline.R (2002-09-12)
3 ## Methods to construct skyline objects (data underlying skyline plot)
5 ## Copyright 2002 Korbinian Strimmer
7 ## This file is part of the R-package `ape'.
8 ## See the file ../COPYING for licensing issues.
10 skyline <- function(x, ...) UseMethod("skyline")
12 # input: phylogenetic tree
13 skyline.phylo <- function(x, ...)
15 if (class(x) != "phylo")
16 stop("object \"x\" is not of class \"phylo\"")
18 skyline(coalescent.intervals(x), ...)
21 # input: coalescent intervals and epsilon
22 skyline.coalescentIntervals <- function(x, epsilon=0, ...)
24 if (class(x) != "coalescentIntervals")
25 stop("object \"x\" is not of class \"coalescentIntervals\"")
29 eps <- find.skyline.epsilon(x, ...)
34 skyline(collapsed.intervals(x, epsilon=eps), ...)
38 # input: collapsed intervals
39 skyline.collapsedIntervals <- function(x, old.style=FALSE, ...)
41 if (class(x) != "collapsedIntervals")
42 stop("object \"x\" is not of class \"collapsedIntervals\"")
44 link <- x$collapsed.interval
45 params <- x$collapsed.interval.count
47 w <- x$interval.length
49 b <- choose(l,2) # binomial coefficients
51 sg <- rep(0,params) # sizes of collapsed intervals
52 cg <- rep(0,params) # coalescent events in interval
55 ng <- rep(0,params) # lineages at beginning of an in interval
58 ng <- rep(0,params) # sum of classic skp estimates in an interval
67 cg[[i]] <- length(sgr)
70 ng[[i]] <- l[group][[1]]
72 ng[[i]] <- sum(m.classic[group])
75 # generalized skp estimate
78 m <- sg*(ng*(ng-cg)/(2.0*cg) )
83 logL <- sum(log(b/m[link]) - b/m[link]*w)
85 # AICc corrected log-likelihood
86 K <- x$collapsed.interval.count
89 logL.AICc <- logL - K- K*(K+1)/(S-K-1)
97 parameter.count=length(t),
100 logL.AICc = logL.AICc
102 class(obj) <- "skyline"
106 # grid search for finding optimal epsilon parameter
107 find.skyline.epsilon <- function(ci, GRID=1000, MINEPS=1e-6, ...)
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.
114 cat("Searching for the optimal epsilon... ")
116 # a grid search is a naive way but still effective of doing this ...
118 size <- ci$interval.count
119 besteps <- ci$total.depth
122 cli <- collapsed.intervals(ci,eps)
123 skpk <- skyline(cli, ...)
124 bestaicc <- skpk$ logL.AICc
125 params <- skpk$parameter.count
127 delta <- besteps/GRID
132 cli <- collapsed.intervals(ci,eps)
133 skpk <- skyline(cli, ...)
134 aicc <- skpk$ logL.AICc
135 params <- skpk$parameter.count
137 if (aicc > bestaicc && params < size-1)
145 cat("epsilon =", besteps, "\n")