Skip to content

Instantly share code, notes, and snippets.

@nmatzke
Last active May 21, 2019 00:50
  • Star 0 You must be signed in to star a gist
  • Fork 0 You must be signed in to fork a gist
Star You must be signed in to star a gist
Save nmatzke/8bb4f9944f582614deecf8d523e5964b to your computer and use it in GitHub Desktop.
# Laying out the likelihood equation in ape::birthdeath
# Calculating LnL of a tree under a birth-death process
#######################################################
# Located in:
# /drives/GDrive/z_help/ClaSSE_LnLs/grok_classe_setup/
#######################################################
# SEE: EQUATION 21, p. 308, of:
# Nee, Sean; May, Robert M.; Harvey, Paul H. (1994). "The reconstructed evolutionary
# process." Philosophical Transactions of the Royal Society B, Biological Sciences. 344,
# 305-311. doi: 10.1098/rstb.1994.0068
# Model
birthRate = 0.2222222
deathRate = 0
a = deathRate / birthRate
r = birthRate - deathRate # =growthRate or diversification rate
# Tree branchlengths (data)
phy = tr
N = length(phy$tip.label)
x <- c(NA, branching.times(phy))
# Deviance:
-2 * (lfactorial(N - 1) + (N - 2) * log(r) + r * sum(x[3:N]) +
N * log(1 - a) - 2 * sum(log(exp(r * x[2:N]) - a)))
# LnL:
# (lfactorial means log(factorial))
# (factorial means gamma(x+1))
(
lfactorial(N - 1)
+ (N - 2) * log(r)
+ r * sum(x[3:N])
+ N * log(1 - a) - 2 * sum(log(exp(r * x[2:N]) - a))
)
# Log (factorial(numtips-1))
lfactorial(N - 1)
log(factorial(N-1))
log(gamma(N))
log(3*2*1)
# log(diversification rate (r=b-d)) added (numtips-2) times
# log(growthRate) = expected waiting time
(N - 2) * log(r)
# N-2 expected waiting times for 2 internal branches
# diversification rate times age of internal nodes (except root)
# (expected amount of diversification on these 2 branches)
r * sum(x[3:N])
# number of tips (number of branches that DIDN'T speciate)
# times 1-p(extinction)
# (waiting times to non-events)
N * log(1 - a)
# Divide by
#
- 2 * sum(log(exp(r * x[2:N]) - a))
# Internal branch indexes
2:N
# ages of internal branches
x[2:N]
# diversification rate times age of internal nodes (including root)
# total amount of diversification rate from each node
r * x[2:N]
# expected amount of diversification from each node
exp(r * x[2:N])
# expected amount of diversification from each node,
# minus number of lineages expected to go extinct from each node
exp(r * x[2:N]) - a
# amount of waiting time from each internal node
sum(log(exp(r * x[2:N]) - a))
# divide by 2 lineages times amount of waiting time from each time point
# - 2 * sum(log(exp(r * x[2:N]) - a))
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment