Nothing
## skyline.R (2002-09-12)
## Methods to construct skyline objects (data underlying skyline plot)
## Copyright 2002 Korbinian Strimmer
## This file is part of the R-package `ape'.
## See the file ../COPYING for licensing issues.
skyline <- function(x, ...) UseMethod("skyline")
# input: phylogenetic tree
skyline.phylo <- function(x, ...)
{
if (!inherits(x, "phylo"))
stop("object \"x\" is not of class \"phylo\"")
skyline(coalescent.intervals(x), ...)
}
# input: coalescent intervals and epsilon
skyline.coalescentIntervals <- function(x, epsilon=0, ...)
{
if (!inherits(x, "coalescentIntervals"))
stop("object \"x\" is not of class \"coalescentIntervals\"")
if (epsilon < 0)
{
eps <- find.skyline.epsilon(x, ...)
}
else
eps <- epsilon
skyline(collapsed.intervals(x, epsilon=eps), ...)
}
# input: collapsed intervals
skyline.collapsedIntervals <- function(x, old.style=FALSE, ...)
{
if (!inherits(x, "collapsedIntervals"))
stop("object \"x\" is not of class \"collapsedIntervals\"")
link <- x$collapsed.interval
params <- x$collapsed.interval.count
l <- x$lineages
w <- x$interval.length
b <- choose(l,2) # binomial coefficients
sg <- rep(0,params) # sizes of collapsed intervals
cg <- rep(0,params) # coalescent events in interval
if(old.style)
ng <- rep(0,params) # lineages at beginning of an in interval
else
{
ng <- rep(0,params) # sum of classic skp estimates in an interval
m.classic <- w*b
}
for (i in 1:params)
{
group <- link==i
sgr <- w[group]
sg[[i]] <- sum(sgr)
cg[[i]] <- length(sgr)
if(old.style)
ng[[i]] <- l[group][[1]]
else
ng[[i]] <- sum(m.classic[group])
}
# generalized skp estimate
t <- cumsum(sg)
if (old.style)
m <- sg*(ng*(ng-cg)/(2.0*cg) )
else
m <- ng/cg
# log-likelihood
logL <- sum(log(b/m[link]) - b/m[link]*w)
# AICc corrected log-likelihood
K <- x$collapsed.interval.count
S <- x$interval.count
if (S-K > 1)
logL.AICc <- logL - K- K*(K+1)/(S-K-1)
else
logL.AICc <- NA
obj <- list(
time=t,
interval.length=sg,
population.size=m,
parameter.count=length(t),
epsilon = x$epsilon,
logL = logL,
logL.AICc = logL.AICc
)
class(obj) <- "skyline"
return(obj)
}
# grid search for finding optimal epsilon parameter
find.skyline.epsilon <- function(ci, GRID=1000, MINEPS=1e-6, ...)
{
# Why MINEPS?
# Because most "clock-like" trees are not properly
# clock-like for a variety of reasons, i.e. the heights
# of the tips are not exactly zero.
cat("Searching for the optimal epsilon... ")
# a grid search is a naive way but still effective of doing this ...
size <- ci$interval.count
besteps <- ci$total.depth
eps <- besteps
cli <- collapsed.intervals(ci,eps)
skpk <- skyline(cli, ...)
bestaicc <- skpk$ logL.AICc
params <- skpk$parameter.count
delta <- besteps/GRID
eps <- eps-delta
while(eps > MINEPS)
{
cli <- collapsed.intervals(ci,eps)
skpk <- skyline(cli, ...)
aicc <- skpk$ logL.AICc
params <- skpk$parameter.count
if (aicc > bestaicc && params < size-1)
{
besteps <- eps
bestaicc <- aicc
}
eps <- eps-delta
}
cat("epsilon =", besteps, "\n")
besteps
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.