Permalink
Browse files

Initial commit

  • Loading branch information...
1 parent 9b08099 commit 0c94f92208c597c7ff5efdb54aaa0d39ef990420 @liamrevell committed Aug 18, 2015
Showing with 14,936 additions and 0 deletions.
  1. +58 −0 DESCRIPTION
  2. +84 −0 NAMESPACE
  3. +21 −0 R/add.everywhere.R
  4. +39 −0 R/add.random.R
  5. +39 −0 R/add.species.to.genus.R
  6. +48 −0 R/allFurcTrees.R
  7. +103 −0 R/anc.Bayes.R
  8. +90 −0 R/anc.ML.R
  9. +50 −0 R/anc.trend.R
  10. +334 −0 R/ancThresh.R
  11. +109 −0 R/backbonePhylo.R
  12. +121 −0 R/bmPlot.R
  13. +48 −0 R/branching.diffusion.R
  14. +113 −0 R/brownie.lite.R
  15. +41 −0 R/brownieREML.R
  16. +216 −0 R/collapseTree.R
  17. +128 −0 R/contMap.R
  18. +145 −0 R/cophylo.R
  19. +197 −0 R/densityMap.R
  20. +67 −0 R/drop.tip.simmap.R
  21. +65 −0 R/estDiversity.R
  22. +379 −0 R/evol.rate.mcmc.R
  23. +151 −0 R/evol.vcv.R
  24. +114 −0 R/evolvcv.lite.R
  25. +81 −0 R/exhaustiveMP.R
  26. +32 −0 R/export.as.xml.R
  27. +258 −0 R/fancyTree.R
  28. +40 −0 R/fastAnc.R
  29. +103 −0 R/fastBM.R
  30. +25 −0 R/fitBayes.R
  31. +109 −0 R/fitDiversityModel.R
  32. +118 −0 R/fitPagel.R
  33. +107 −0 R/locate.fossil.R
  34. +183 −0 R/locate.yeti.R
  35. +202 −0 R/ltt.R
  36. +105 −0 R/ltt95.R
  37. +53 −0 R/make.era.map.R
  38. +425 −0 R/make.simmap.R
  39. +42 −0 R/map.overlap.R
  40. +111 −0 R/map.to.singleton.R
  41. +113 −0 R/mcmcBM.R
  42. +116 −0 R/mcmcBM.full.R
  43. +137 −0 R/mcmcLambda.R
  44. +87 −0 R/mrp.supertree.R
  45. +64 −0 R/multi.mantel.R
  46. +20 −0 R/multiRF.R
  47. +148 −0 R/optim.phylo.ls.R
  48. +94 −0 R/paintSubTree.R
  49. +259 −0 R/pbtree.R
  50. +87 −0 R/pgls.Ives.R
  51. +154 −0 R/phenogram.R
  52. +35 −0 R/phyl.RMA.R
  53. +104 −0 R/phyl.cca.R
  54. +57 −0 R/phyl.pairedttest.R
  55. +153 −0 R/phyl.pca.R
  56. +57 −0 R/phyl.resid.R
  57. +68 −0 R/phylANOVA.R
  58. +128 −0 R/phylo.to.map.R
  59. +111 −0 R/phylomorphospace.R
  60. +95 −0 R/phylomorphospace3d.R
  61. +138 −0 R/phylosig.R
  62. +124 −0 R/plotBranchbyTrait.R
  63. +340 −0 R/plotSimmap.R
  64. +79 −0 R/plotTree.wBars.R
  65. +70 −0 R/ratebystate.R
  66. +121 −0 R/rateshift.R
  67. +124 −0 R/read.newick.R
  68. +347 −0 R/read.simmap.R
  69. +47 −0 R/rerootingMethod.R
  70. +67 −0 R/roundPhylogram.R
  71. +45 −0 R/sim.corrs.R
  72. +107 −0 R/sim.history.R
  73. +30 −0 R/sim.rates.R
  74. +37 −0 R/skewers.R
  75. +97 −0 R/splitplotTree.R
  76. +13 −0 R/starTree.R
  77. +49 −0 R/strahlerNumber.R
  78. +139 −0 R/threshBayes.R
  79. +49 −0 R/treeSlice.R
  80. +1,247 −0 R/utilities.R
  81. +73 −0 R/write.simmap.R
  82. +111 −0 R/writeAncestors.R
  83. +29 −0 R/writeNexus.R
  84. BIN data/anoletree.rda
  85. +15 −0 inst/CITATION
  86. +26 −0 man/add.arrow.Rd
  87. +31 −0 man/add.color.bar.Rd
  88. +25 −0 man/add.everywhere.Rd
  89. +31 −0 man/add.random.Rd
  90. +26 −0 man/add.simmap.legend.Rd
  91. +30 −0 man/add.species.to.genus.Rd
  92. +31 −0 man/allFurcTrees.Rd
  93. +34 −0 man/anc.Bayes.Rd
  94. +43 −0 man/anc.ML.Rd
  95. +40 −0 man/anc.trend.Rd
  96. +39 −0 man/ancThresh.Rd
  97. +16 −0 man/anoletree.Rd
  98. +22 −0 man/applyBranchLengths.Rd
  99. +32 −0 man/ave.rates.Rd
  100. +28 −0 man/bind.tip.Rd
  101. +44 −0 man/bmPlot.Rd
  102. +34 −0 man/branching.diffusion.Rd
  103. +49 −0 man/brownie.lite.Rd
  104. +35 −0 man/brownieREML.Rd
  105. +33 −0 man/cladelabels.Rd
  106. +25 −0 man/collapse.to.star.Rd
  107. +29 −0 man/collapseTree.Rd
  108. +46 −0 man/contMap.Rd
  109. +41 −0 man/cophylo.Rd
  110. +30 −0 man/countSimmap.Rd
  111. +47 −0 man/densityMap.Rd
  112. +28 −0 man/describe.simmap.Rd
  113. +29 −0 man/di2multi.simmap.Rd
  114. +25 −0 man/drop.clade.Rd
  115. +22 −0 man/drop.leaves.Rd
  116. +30 −0 man/drop.tip.contMap.Rd
  117. +30 −0 man/drop.tip.simmap.Rd
  118. +36 −0 man/estDiversity.Rd
  119. +35 −0 man/evol.rate.mcmc.Rd
  120. +44 −0 man/evol.vcv.Rd
  121. +30 −0 man/evolvcv.lite.Rd
  122. +34 −0 man/exhaustiveMP.Rd
  123. +20 −0 man/expm.Rd
  124. +31 −0 man/export.as.xml.Rd
  125. +59 −0 man/fancyTree.Rd
  126. +33 −0 man/fastAnc.Rd
  127. +36 −0 man/fastBM.Rd
  128. +35 −0 man/fastMRCA.Rd
  129. +33 −0 man/findMRCA.Rd
  130. +32 −0 man/fitBayes.Rd
  131. +35 −0 man/fitDiversityModel.Rd
  132. +30 −0 man/fitPagel.Rd
  133. +37 −0 man/gammatest.Rd
  134. +33 −0 man/genSeq.Rd
  135. +25 −0 man/getCladesofSize.Rd
  136. +26 −0 man/getDescendants.Rd
  137. +27 −0 man/getExtant.Rd
  138. +23 −0 man/getSisters.Rd
  139. +25 −0 man/getStates.Rd
  140. +25 −0 man/ladderize.simmap.Rd
  141. +26 −0 man/lambda.transform.Rd
  142. +28 −0 man/likMlambda.Rd
  143. +30 −0 man/locate.fossil.Rd
  144. +28 −0 man/locate.yeti.Rd
  145. +26 −0 man/ls.tree.Rd
  146. +48 −0 man/ltt.Rd
  147. +34 −0 man/ltt95.Rd
  148. +35 −0 man/make.era.map.Rd
  149. +58 −0 man/make.simmap.Rd
  150. +26 −0 man/map.overlap.Rd
  151. +29 −0 man/map.to.singleton.Rd
  152. +28 −0 man/markChanges.Rd
  153. +27 −0 man/matchNodes.Rd
  154. +29 −0 man/mergeMappedStates.Rd
  155. +31 −0 man/midpoint.root.Rd
  156. +28 −0 man/minRotate.Rd
  157. +32 −0 man/minSplit.Rd
  158. +41 −0 man/mrp.supertree.Rd
  159. +34 −0 man/multi.mantel.Rd
  160. +24 −0 man/multiC.Rd
  161. +26 −0 man/multiRF.Rd
  162. +35 −0 man/nodeHeights.Rd
  163. +43 −0 man/optim.phylo.ls.Rd
  164. +22 −0 man/orderMappedEdge.Rd
  165. +33 −0 man/paintSubTree.Rd
  166. +32 −0 man/paste.tree.Rd
  167. +40 −0 man/pbtree.Rd
  168. +45 −0 man/pgls.Ives.Rd
  169. +49 −0 man/phenogram.Rd
  170. +41 −0 man/phyl.RMA.Rd
  171. +41 −0 man/phyl.cca.Rd
  172. +43 −0 man/phyl.pairedttest.Rd
  173. +44 −0 man/phyl.pca.Rd
  174. +41 −0 man/phyl.resid.Rd
  175. +27 −0 man/phyl.vcv.Rd
  176. +43 −0 man/phylANOVA.Rd
  177. +29 −0 man/phylo.to.map.Rd
  178. +29 −0 man/phylo.toBackbone.Rd
  179. +25 −0 man/phyloDesign.Rd
  180. +37 −0 man/phylomorphospace.Rd
  181. +43 −0 man/phylomorphospace3d.Rd
  182. +57 −0 man/phylosig.Rd
  183. +23 −0 man/phytools-package.Rd
  184. +50 −0 man/plot.backbonePhylo.Rd
  185. +32 −0 man/plotBranchbyTrait.Rd
  186. +65 −0 man/plotSimmap.Rd
  187. +34 −0 man/plotThresh.Rd
  188. +29 −0 man/plotTree.Rd
  189. +33 −0 man/plotTree.wBars.Rd
  190. +29 −0 man/posterior.evolrate.Rd
  191. +26 −0 man/print.backbonePhylo.Rd
  192. +33 −0 man/ratebystate.Rd
  193. +28 −0 man/rateshift.Rd
  194. +32 −0 man/read.newick.Rd
  195. +38 −0 man/read.simmap.Rd
  196. +27 −0 man/reorder.backbonePhylo.Rd
  197. +28 −0 man/reorderSimmap.Rd
  198. +35 −0 man/rep.phylo.Rd
  199. +31 −0 man/reroot.Rd
  200. +39 −0 man/rerootingMethod.Rd
  201. +28 −0 man/rescaleSimmap.Rd
  202. +27 −0 man/rotateNodes.Rd
  203. +22 −0 man/roundBranches.Rd
  204. +38 −0 man/roundPhylogram.Rd
  205. +25 −0 man/rstate.Rd
  206. +25 −0 man/sampleFrom.Rd
  207. +28 −0 man/setMap.Rd
  208. +30 −0 man/sim.corrs.Rd
  209. +32 −0 man/sim.history.Rd
  210. +29 −0 man/sim.ratebystate.Rd
  211. +30 −0 man/sim.rates.Rd
  212. +31 −0 man/skewers.Rd
  213. +25 −0 man/splitTree.Rd
  214. +35 −0 man/splitplotTree.Rd
  215. +28 −0 man/starTree.Rd
  216. +25 −0 man/strahlerNumber.Rd
  217. +33 −0 man/threshBayes.Rd
  218. +32 −0 man/threshDIC.Rd
  219. +28 −0 man/threshState.Rd
  220. +21 −0 man/to.matrix.Rd
  221. +28 −0 man/treeSlice.Rd
  222. +22 −0 man/untangle.Rd
  223. +24 −0 man/vcvPhylo.Rd
  224. +42 −0 man/write.simmap.Rd
  225. +30 −0 man/writeAncestors.Rd
  226. +25 −0 man/writeNexus.Rd
View
@@ -0,0 +1,58 @@
+Package: phytools
+Version: 0.4-99
+Date: 2015-8-18
+Title: Phylogenetic Tools for Comparative Biology (and Other Things)
+Author: Liam J. Revell
+Maintainer: Liam J. Revell <liam.revell@umb.edu>
+Depends: R (>= 2.10), ape (>= 3.0-10), maps
+Imports: animation, clusterGeneration, graphics, grDevices, methods,
+ mnormt, msm, numDeriv, phangorn (>= 1.6-3), plotrix,
+ scatterplot3d, stats, utils
+Suggests: geiger, rgl
+ZipData: no
+Description: Package contains various functions for phylogenetic analysis.
+ This functionality is concentrated in the phylogenetic analysis of
+ comparative data from species. For example, the package includes
+ functions for Bayesian and ML ancestral state estimation; visual
+ simulation of trait evolution; fitting models of trait evolution
+ with multiple Brownian rates and correlations; visualizing
+ discrete and continuous character evolution using colors or
+ projections into trait space; identifying the location of a change
+ in the rate of character evolution on the tree; fast Brownian motion
+ simulation and simulation under several other models of
+ continuous trait evolution; fitting a model of correlated binary
+ trait evolution; locating the position of a fossil or an recently
+ extinct lineage on a tree using continuous character data with ML;
+ plotting lineage accumulation through time, including across
+ multiple trees (such as a Bayesian posterior sample); conducting
+ an analysis called stochastic character mapping, in which character
+ histories for a discrete trait are sampled from their posterior
+ probability distribution under a model; conducting a multiple
+ (i.e., partial) Mantel test; fitting a phylogenetic regression model
+ with error in predictor and response variables; conducting a
+ phylogenetic principal components analysis, a phylogenetic
+ regression, a reduced major axis regression, a phylogenetic
+ canonical correlation analysis, and a phylogenetic ANOVA; projecting
+ a tree onto a geographic map; simulating discrete character
+ histories on the tree; and fitting a model in which a discrete
+ character evolves under the threshold model. In addition to this
+ phylogenetic comparative method functionality, the package also
+ contains functions for a wide range of other purposes in
+ phylogenetic biology. For instance, functionality in this package
+ includes (but is not restricted to): adding taxa to a tree
+ (including randomly, everywhere, or automatically to genera);
+ generating all bi- and multi-furcating trees for a set of taxa;
+ reducing a phylogeny to its backbone tree; dropping tips or adding
+ tips to special types of phylogenetic trees; exporting a tree as an
+ XML file; converting a tree with a mapped character to a tree with
+ singleton nodes and one character state per edge; estimating a
+ phylogeny using the least squares method; simulating birth-death
+ trees under a range of conditions; rerooting trees; a wide range
+ of visualizations of trees; and a variety of other manipulations
+ and analyses that phylogenetic biologists may find useful for
+ their research.
+License: GPL (>= 2)
+URL: http://www.phytools.org
+Packaged: 2015-8-18 12:00:00 EDT
+Repository:
+Date/Publication: 2015-8-18 12:00:00 EDT
View
@@ -0,0 +1,84 @@
+export(add.arrow, add.color.bar, add.everywhere, add.random, add.simmap.legend, add.species.to.genus, allFurcTrees, anc.Bayes)
+export(anc.ML, anc.trend, ancThresh, applyBranchLengths, ave.rates)
+export(backbone.toPhylo, bmPlot, bind.tip, biplot.phyl.pca, branching.diffusion, brownie.lite, brownieREML)
+export(cladelabels, collapse.to.star, collapseTree, contMap, cophylo, countSimmap)
+export(densityMap, describe.simmap, di2multi.simmap, drop.clade, drop.leaves, drop.tip.contMap, drop.tip.densityMap)
+export(drop.tip.simmap, drop.tip.singleton)
+export(estDiversity, evol.rate.mcmc, evol.vcv, evolvcv.lite, exhaustiveMP, expm, export.as.xml, extract.clade.simmap)
+export(extract.strahlerNumber)
+export(fancyTree, fastAnc, fastBM, fastHeight, fastMRCA, findMRCA, fitBayes, fitDiversityModel, fitPagel)
+export(gammatest, genSeq, getCladesofSize, getDescendants, getExtant, getExtinct, getSisters, getStates)
+export(ladderize.simmap, lambda.transform, likMlambda, locate.fossil, locate.yeti, ls.tree, ltt, ltt95)
+export(make.era.map, make.simmap, map.overlap, map.to.singleton, markChanges, matchNodes, mergeMappedStates)
+export(midpoint.root, minRotate, minSplit, mrp.supertree, multi.mantel, multiC, multiRF)
+export(nodeheight, nodeHeights)
+export(optim.phylo.ls, orderMappedEdge)
+export(paintBranches, paintSubTree, paste.tree, pbtree, pgls.Ives,phenogram, phyl.cca, phyl.pairedttest, phyl.pca, phyl.resid, phyl.RMA)
+export(phyl.vcv, phylo.toBackbone,phylo.to.map, phylANOVA, phyloDesign, phylomorphospace, phylomorphospace3d, phylosig)
+export(plotBranchbyTrait, plot.contMap, plot.cophylo, plot.densityMap, plot.phylo.to.map, plotSimmap, plotThresh, plotTree)
+export(plotTree.singletons, plotTree.splits, plotTree.wBars, posterior.evolrate)
+export(ratebystate, rateshift, read.newick, read.simmap, reorder.backbonePhylo, reorderSimmap, rep.multiPhylo, rep.phylo, repPhylo)
+export(reroot, rerootingMethod, rescaleSimmap, rotateNodes, roundBranches, roundPhylogram, rstate)
+export(sampleFrom, setMap, sim.corrs, sim.history, sim.ratebystate, sim.rates, skewers, splitplotTree)
+export(splitTree, starTree, strahlerNumber)
+export(threshBayes, threshDIC, threshState, tipRotate, to.matrix, treeSlice)
+export(untangle)
+export(vcvPhylo)
+export(write.simmap, writeAncestors, writeNexus)
+
+S3method(plot, densityMap)
+S3method(plot, contMap)
+S3method(plot, phylo.to.map)
+S3method(plot, backbonePhylo)
+S3method(reorder, backbonePhylo)
+S3method(print, backbonePhylo)
+S3method(print, phyl.pca)
+S3method(summary, phyl.pca)
+S3method(biplot, phyl.pca)
+S3method(print, brownie.lite)
+S3method(print, densityMap)
+S3method(print, contMap)
+S3method(print, rateshift)
+S3method(logLik, rateshift)
+S3method(print, evol.vcv)
+S3method(plot, ltt95)
+S3method(print, ltt95)
+S3method(plot, describe.simmap)
+S3method(print, describe.simmap)
+S3method(rep, phylo)
+S3method(rep, multiPhylo)
+S3method(print, fitPagel)
+S3method(print, ltt)
+S3method(print, multiLtt)
+S3method(plot, ltt)
+S3method(plot, multiLtt)
+S3method(plot, cophylo)
+S3method(print, cophylo)
+S3method(plot, simmap)
+S3method(print, simmap)
+S3method(summary, simmap)
+S3method(plot, multiSimmap)
+S3method(print, multiSimmap)
+S3method(summary, multiSimmap)
+S3method(reorder, simmap)
+
+importFrom(animation, ani.options, ani.record, ani.replay, saveVideo)
+importFrom(ape, .PlotPhyloEnv, .uncompressTipLabel, ace, all.equal.phylo, as.DNAbin, bind.tree, branching.times, collapse.singles)
+importFrom(ape, compute.brlen, cophenetic.phylo, di2multi, dist.dna, dist.nodes, drop.tip, edgelabels, extract.clade, is.binary.tree)
+importFrom(ape, is.monophyletic, is.rooted, is.ultrametric, ladderize, matexpo, mrca, multi2di, neworder_phylo, neworder_pruningwise)
+importFrom(ape, nodelabels, Ntip, pic, plot.phylo, prop.part, read.tree, reorder.phylo, root, rotate, rtree, stree, tiplabels)
+importFrom(ape, unroot, vcv, vcv.phylo, write.tree)
+importFrom(clusterGeneration, genPositiveDefMat)
+importFrom(maps, map)
+importFrom(mnormt, dmnorm, pd.solve)
+importFrom(msm, MatrixExp)
+importFrom(numDeriv, hessian)
+importFrom(phangorn, allTrees, Ancestors, as.phyDat, dist.hamming, NJ, nni, optim.parsimony, parsimony, phyDat, pratchet)
+importFrom(plotrix, draw.arc, draw.circle, textbox)
+importFrom(scatterplot3d, scatterplot3d)
+importFrom(stats, cophenetic, reorder, runif, setNames, optim, dexp, dnorm, rnorm, var, nlminb, biplot, pchisq, rexp, optimize, logLik, as.dist, pnorm, median, cor, aggregate, rgamma, dgamma, lm, rbinom, rgeom, pt, anova, p.adjust, screeplot, rchisq, optimHess, rmultinom, sd, dunif)
+importFrom(methods, hasArg)
+importFrom(graphics, strwidth, par, segments, locator, lines, text, strheight, symbols, plot, layout, plot.new, title, axis, points, plot.window, polygon, rect, curve)
+importFrom(utils, flush.console, installed.packages, str, capture.output)
+importFrom(grDevices, palette, colors, dev.hold, dev.flush, rainbow, heat.colors, gray, colorRampPalette, dev.new, rgb)
+
View
@@ -0,0 +1,21 @@
+# function takes a tree and adds a tip in all possible places
+# returns the set of unrooted trees without branch lengths
+# written by Liam J. Revell 2011, 2013, 2015
+
+add.everywhere<-function(tree,tip.name){
+ if(!inherits(tree,"phylo")) stop("tree should be an object of class \"phylo\".")
+ tree<-unroot(tree) # unroot tree
+ tree$edge.length<-rep(1,nrow(tree$edge)) # set all edge lengths to 1.0
+ # create new tip
+ new.tip<-list(edge=matrix(c(2L,1L),1,2),tip.label=tip.name,edge.length=1,Nnode=1L)
+ class(new.tip)<-"phylo"
+ # add the new tip to all edges of the tree
+ trees<-list()
+ class(trees)<-"multiPhylo"
+ for(i in 1:nrow(tree$edge)){
+ trees[[i]]<-bind.tree(tree,new.tip,where=tree$edge[i,2],position=0.5)
+ trees[[i]]$edge.length<-NULL
+ }
+ return(trees)
+}
+
View
@@ -0,0 +1,39 @@
+# function adds a set of tips to random positions in the tree
+# written by Liam J. Revell 2013, 2015
+
+add.random<-function(tree,n=NULL,tips=NULL,edge.length=NULL,order=c("random","input")){
+ if(!inherits(tree,"phylo")) stop("tree should be an object of class \"phylo\".")
+ # internal function
+ randomPosn<-function(tree){
+ # sum edges cumulatively
+ cum.edge<-cumsum(tree$edge.length)
+ index<-tree$edge[,2]
+ # pick random position
+ pos<-runif(1)*sum(tree$edge.length)
+ edge<-1; while(pos>cum.edge[edge]) edge<-edge+1
+ return(list(node=index[edge],posn=cum.edge[edge]-pos))
+ }
+ # check if tree is ultrametric (required)
+ if(is.ultrametric(tree)) um<-TRUE
+ else um<-FALSE
+ # set n and/or name tips (if not provided)
+ if(is.null(tips)){
+ if(is.null(n)) n<-1
+ tips<-paste("t",length(tree$tip)+1:n,sep="")
+ } else n<-length(tips)
+ if(is.null(edge.length)) if(!um) edge.length<-runif(n=n,min=min(tree$edge.length),max=max(tree$edge.length))
+ # set order
+ if(order[1]=="random"){
+ o<-sample(1:n)
+ tips<-tips[o]
+ if(!is.null(edge.length)) edge.length<-edge.length[o]
+ }
+ # add tips
+ for(i in 1:n){
+ where<-randomPosn(tree)
+ if(is.null(edge.length)) tree<-bind.tip(tree,tips[i],where=where$node,position=where$posn)
+ else tree<-bind.tip(tree,tips[i],where=where$node,position=where$posn,edge.length=edge.length[i])
+ }
+ return(tree)
+}
+
@@ -0,0 +1,39 @@
+## function identifies the genus from a species name & attempts to add it to the tree
+## written by Liam J. Revell 2013
+
+add.species.to.genus<-function(tree,species,genus=NULL,where=c("root","random")){
+ if(!inherits(tree,"phylo")) stop("tree should be an object of class \"phylo\".")
+ if(!is.ultrametric(tree))
+ warning("this code has only been tested with ultrametric tree\n your tree may be returned without edge lengths")
+ where<-where[1]
+ if(is.null(genus)){
+ ## get genus from species name
+ x<-strsplit(species,"")[[1]]
+ i<-1
+ while(x[i]!="_"&&x[i]!=" ") i<-i+1
+ genus<-paste(x[2:i-1],collapse="")
+ }
+ ii<-grep(paste(genus,"_",sep=""),tree$tip.label)
+ if(length(ii)>1){
+ if(!is.monophyletic(tree,tree$tip.label[ii]))
+ warning(paste(genus,"may not be monophyletic\n attaching to the most inclusive group containing members of this genus"))
+ nn<-findMRCA(tree,tree$tip.label[ii])
+ if(where=="root") tree<-bind.tip(tree,gsub(" ","_",species),where=nn)
+ else if(where=="random"){
+ tt<-splitTree(tree,list(node=nn,bp=tree$edge.length[which(tree$edge[,2]==nn)]))
+ tt[[2]]<-add.random(tt[[2]],tips=gsub(" ","_",species))
+ tree<-paste.tree(tt[[1]],tt[[2]])
+ } else stop("option 'where' not recognized")
+ } else if(length(ii)==1){
+ nn<-ii
+ if(where=="root")
+ tree<-bind.tip(tree,gsub(" ","_",species),where=nn,position=0.5*tree$edge.length[which(tree$edge[,2]==nn)])
+ else if(where=="random")
+ tree<-bind.tip(tree,gsub(" ","_",species),where=nn,position=runif(n=1)*tree$edge.length[which(tree$edge[,2]==nn)])
+ else
+ stop("option 'where' not recognized")
+ } else
+ warning("could not match your species to a genus\n check spelling, including case")
+ tree
+}
+
View
@@ -0,0 +1,48 @@
+# function returns all unrooted multi & bifurcating tree topologies
+# written by Liam Revell 2011, 2013
+
+allFurcTrees<-function(n,tip.label=NULL,to.plot=TRUE){
+ # check to see if tip labels have been provided
+ if(is.null(tip.label)) tip.label=as.character(1:n)
+ # now pick three species at random to start
+ new<-list(stree(n=3,tip.label=sample(tip.label,3)))
+ class(new)<-"multiPhylo"
+ added<-new[[1]]$tip.label; remaining<-setdiff(tip.label,added)
+ # loop
+ while(length(remaining)>0){
+ old<-new; new<-list()
+ new.tip<-sample(remaining,1)
+ for(i in 1:length(old)){
+ temp<-add.to.branches.and.nodes(old[[i]],new.tip)
+ new<-unlist(list(new,temp),recursive=FALSE); class(new)<-"multiPhylo"
+ }
+ added<-c(added,new.tip)
+ remaining<-setdiff(tip.label,added)
+ }
+ if(to.plot) new<-read.tree(text=write.tree(new))
+ return(new) # return all trees
+}
+
+# function adds a tip to all branches and nodes
+# written by Liam Revell 2011
+
+add.to.branches.and.nodes<-function(tree,tip.name){
+ n.edge<-nrow(tree$edge)
+ tree$edge.length<-rep(1,n.edge) # set all edge lengths to 1.0
+ # create new tip
+ new.tip<-list(edge=matrix(c(2L,1L),1,2),tip.label=tip.name,edge.length=1,Nnode=1L); class(new.tip)<-"phylo"
+ # first add the tip to all edges
+ trees<-list(); class(trees)<-"multiPhylo"
+ for(i in 1:n.edge){
+ trees[[i]]<-bind.tree(tree,new.tip,where=tree$edge[i,2],position=0.5) # add tip to edge
+ trees[[i]]$edge.length<-NULL
+ }
+ # ok, now add the tip to all internal nodes
+ intnodes<-unique(tree$edge[,1])
+ for(i in 1:tree$Nnode){
+ trees[[i+n.edge]]<-bind.tree(tree,new.tip,where=intnodes[i]) # add tip to node
+ trees[[i+n.edge]]$edge.length<-NULL
+ }
+ return(trees)
+}
+
View
@@ -0,0 +1,103 @@
+# function does Bayes ancestral character estimation
+# written by Liam J. Revell 2011, 2013, 2015
+
+anc.Bayes<-function(tree,x,ngen=10000,control=list()){
+ if(!inherits(tree,"phylo")) stop("tree should be an object of class \"phylo\".")
+ # give the function some defaults (in case none are provided)
+ temp<-phyl.vcv(as.matrix(x),vcv(tree),1)
+ sig2<-temp$R[1,1]
+ a<-temp$alpha
+ y<-rep(a,tree$Nnode-1)
+ pr.mean<-c(1000,rep(0,tree$Nnode))
+ pr.var<-c(pr.mean[1]^2,rep(1000,tree$Nnode))
+ prop<-rep(0.01*max(temp$C)*sig2,tree$Nnode+1)
+
+ # populate control list
+ con=list(sig2=sig2,a=a,y=y,pr.mean=pr.mean,pr.var=pr.var,prop=prop,sample=100)
+ con[(namc<-names(control))]<-control
+ con<-con[!sapply(con,is.null)]
+ # print control parameters to screen
+ message("Control parameters (set by user or default):"); str(con)
+
+ # function returns the log-likelihood
+ likelihood<-function(C,invC,detC,x,sig2,a,y){
+ z<-c(x,y)-a
+ logLik<--z%*%invC%*%z/(2*sig2)-nrow(C)*log(2*pi)/2-nrow(C)*log(sig2)/2-detC/2
+ return(logLik)
+ }
+
+ # function returns the prior
+ log.prior<-function(pr.mean,pr.var,sig2,a,y)
+ pp<-dexp(sig2,rate=1/pr.mean[1],log=TRUE)+sum(dnorm(c(a,y),mean=pr.mean[2:length(pr.mean)],sd=sqrt(pr.var[1+1:tree$Nnode]),log=TRUE))
+
+ # compute C
+ C<-vcvPhylo(tree)
+ # check to make sure that C will be non-singular
+ if(any(tree$edge.length<=(10*.Machine$double.eps)))
+ stop("some branch lengths are 0 or nearly zero")
+ invC<-solve(C)
+ detC<-determinant(C,logarithm=TRUE)$modulus[1]
+
+ # now set starting values for MCMC
+ sig2<-con$sig2; a<-con$a; y<-con$y
+ x<-x[tree$tip.label]
+ if(is.null(names(y))) names(y)<-length(tree$tip)+2:tree$Nnode
+ else y[as.character(length(tree$tip)+2:tree$Nnode)]
+ L<-likelihood(C,invC,detC,x,sig2,a,y)
+ Pr<-log.prior(con$pr.mean,con$pr.var,sig2,a,y)
+
+ # store
+ X<-matrix(NA,ngen/con$sample+1,tree$Nnode+3,dimnames=list(NULL,c("gen","sig2",length(tree$tip)+1:tree$Nnode,"logLik")))
+ X[1,]<-c(0,sig2,a,y,L)
+
+ message("Starting MCMC...")
+
+ # start MCMC
+ for(i in 1:ngen){
+ j<-(i-1)%%(tree$Nnode+1)
+ if(j==0){
+ # update sig2
+ sig2.prime<-sig2+rnorm(n=1,sd=sqrt(con$prop[j+1]))
+ if(sig2.prime<0) sig2.prime<--sig2.prime
+ L.prime<-likelihood(C,invC,detC,x,sig2.prime,a,y)
+ Pr.prime<-log.prior(con$pr.mean,con$pr.var,sig2.prime,a,y)
+ post.odds<-min(1,exp(Pr.prime+L.prime-Pr-L))
+ if(post.odds>runif(n=1)){
+ if(i%%con$sample==0) X[i/con$sample+1,]<-c(i,sig2.prime,a,y,L.prime)
+ sig2<-sig2.prime
+ L<-L.prime
+ Pr<-Pr.prime
+ } else if(i%%con$sample==0) X[i/con$sample+1,]<-c(i,sig2,a,y,L)
+ } else if(j==1){
+ # update a
+ a.prime<-a+rnorm(n=1,sd=sqrt(con$prop[j+1]))
+ L.prime<-likelihood(C,invC,detC,x,sig2,a.prime,y)
+ Pr.prime<-log.prior(con$pr.mean,con$pr.var,sig2,a.prime,y)
+ post.odds<-min(1,exp(Pr.prime+L.prime-Pr-L))
+ if(post.odds>runif(n=1)){
+ if(i%%con$sample==0) X[i/con$sample+1,]<-c(i,sig2,a.prime,y,L.prime)
+ a<-a.prime
+ L<-L.prime
+ Pr<-Pr.prime
+ } else if(i%%con$sample==0) X[i/con$sample+1,]<-c(i,sig2,a,y,L)
+ } else {
+ k<-j-1 # update node k
+ y.prime<-y
+ y.prime[k]<-y[k]+rnorm(n=1,sd=sqrt(con$prop[j+1]))
+ L.prime<-likelihood(C,invC,detC,x,sig2,a,y.prime)
+ Pr.prime<-log.prior(con$pr.mean,con$pr.var,sig2,a,y.prime)
+ post.odds<-min(1,exp(Pr.prime+L.prime-Pr-L))
+ if(post.odds>runif(n=1)){
+ if(i%%con$sample==0) X[i/con$sample+1,]<-c(i,sig2,a,y.prime,L.prime)
+ y<-y.prime
+ L<-L.prime
+ Pr<-Pr.prime
+ } else if(i%%con$sample==0) X[i/con$sample+1,]<-c(i,sig2,a,y,L)
+ }
+ }
+
+ # done MCMC
+ message("Done MCMC.")
+ return(X)
+}
+
Oops, something went wrong.

0 comments on commit 0c94f92

Please sign in to comment.