Advertisement
Guest User

Untitled

a guest
Feb 20th, 2018
80
0
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
text 2.03 KB | None | 0 0
  1.  
  2. BMpruneLL <- function(traits, sig, tree){
  3. LLs <- 0
  4. for(i in 1:(length(tree$tip.label) - 1)){
  5. #Find two sister tips
  6. #preterminal nodes
  7. ptn <- tree$edge[sapply(1:length(tree$tip.label), function(t) which(tree$edge[,2] == t)),1]
  8. uniq <- unique(ptn)
  9. parent <- uniq[which(sapply(1:length(uniq), function (x) sum(uniq[x] == ptn)) > 1)][1]
  10. sisters <- tree$edge[tree$edge[,1] == parent, 2]
  11. sisters <- intersect(sisters, 1:(length(tree$tip.label)))
  12.  
  13. #calculate contrast between two sister nodes
  14. contrast <- traits[tree$tip.label[sisters[1]],] - traits[tree$tip.label[sisters[2]],]
  15.  
  16. #calculate BL separating two sister nodes
  17. BLength <- tree$edge.length[tree$edge[,2] == sisters[1]] + tree$edge.length[tree$edge[,2] == sisters[2]]
  18.  
  19. #calculate multivariate normal density of contrast along Branch Length
  20. LLs[i] <- dmvnorm(x = contrast, sigma = BLength*sig, log = T)
  21.  
  22. if(length(tree$tip.label) > 2){
  23. #Fix Tree and Trait Matrix
  24. tipName = paste(tree$tip.label[sisters], collapse = "+")
  25. #Fix Trait Matrix
  26. 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
  27. traits <- rbind(traits, nodeValues)
  28. rownames(traits)[length(traits[,1])] <- tipName
  29. traits <- traits[-c(which(rownames(traits) == tree$tip.label[sisters[1]]), which(rownames(traits) == tree$tip.label[sisters[2]])),]
  30.  
  31. #Fix Tree
  32. newTip <- list(edge=matrix(c(2,1),1,2),
  33. tip.label=tipName,
  34. 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]],
  35. Nnode=1)
  36. class(newTip)<-"phylo"
  37. tree <- bind.tree(x = tree, y = newTip, where = sisters[1])
  38. tree <- drop.tip(tree, sisters[2])
  39. }
  40. }
  41. return(sum(LLs))
  42. }
Advertisement
Add Comment
Please, Sign In to add comment
Advertisement