Advertisement
Not a member of Pastebin yet?
Sign Up,
it unlocks many cool features!
- BMpruneLL <- function(traits, sig, tree){
- LLs <- 0
- for(i in 1:(length(tree$tip.label) - 1)){
- #Find two sister tips
- #preterminal nodes
- ptn <- tree$edge[sapply(1:length(tree$tip.label), function(t) which(tree$edge[,2] == t)),1]
- uniq <- unique(ptn)
- parent <- uniq[which(sapply(1:length(uniq), function (x) sum(uniq[x] == ptn)) > 1)][1]
- sisters <- tree$edge[tree$edge[,1] == parent, 2]
- sisters <- intersect(sisters, 1:(length(tree$tip.label)))
- #calculate contrast between two sister nodes
- contrast <- traits[tree$tip.label[sisters[1]],] - traits[tree$tip.label[sisters[2]],]
- #calculate BL separating two sister nodes
- BLength <- tree$edge.length[tree$edge[,2] == sisters[1]] + tree$edge.length[tree$edge[,2] == sisters[2]]
- #calculate multivariate normal density of contrast along Branch Length
- LLs[i] <- dmvnorm(x = contrast, sigma = BLength*sig, log = T)
- if(length(tree$tip.label) > 2){
- #Fix Tree and Trait Matrix
- tipName = paste(tree$tip.label[sisters], collapse = "+")
- #Fix Trait Matrix
- nodeValues <- (traits[tree$tip.label[sisters[1]],]*tree$edge.length[tree$edge[,2] == sisters[2]] + traits[tree$tip.label[sisters[2]],]*tree$edge.length[tree$edge[,2] == sisters[1]])/BLength
- traits <- rbind(traits, nodeValues)
- rownames(traits)[length(traits[,1])] <- tipName
- traits <- traits[-c(which(rownames(traits) == tree$tip.label[sisters[1]]), which(rownames(traits) == tree$tip.label[sisters[2]])),]
- #Fix Tree
- newTip <- list(edge=matrix(c(2,1),1,2),
- tip.label=tipName,
- edge.length=(prod(tree$edge.length[tree$edge[,2] == sisters[1]], tree$edge.length[tree$edge[,2] == sisters[2]])/BLength) - tree$edge.length[tree$edge[,2] == sisters[1]],
- Nnode=1)
- class(newTip)<-"phylo"
- tree <- bind.tree(x = tree, y = newTip, where = sisters[1])
- tree <- drop.tip(tree, sisters[2])
- }
- }
- return(sum(LLs))
- }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement