birthdeath_treemodel.R

# 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))