new function to compute the probabilities of changes along edges of t…

…he tree from stochastic mapping

NAMESPACE

@@ -4,7 +4,7 @@ export(backbone.toPhylo, bmPlot, bind.tip, biplot.phyl.pca, branching.diffusion,
 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(edgelabels.cophylo, estDiversity, evol.rate.mcmc, evol.vcv, evolvcv.lite, exhaustiveMP, expm, export.as.xml)
+export(edgelabels.cophylo, edgeProbs, estDiversity, evol.rate.mcmc, evol.vcv, evolvcv.lite, exhaustiveMP, expm, export.as.xml)
 export(extract.clade.simmap, extract.strahlerNumber)
 export(fancyTree, fastAnc, fastBM, fastDist, fastHeight, fastMRCA, findMRCA, fitBayes, fitMk, fitDiversityModel, fitPagel)
 export(gammatest, genus.to.species.tree, genSeq, getCladesofSize, getDescendants, getExtant, getExtinct, getSisters, getStates)

R/utilities.R

@@ -1,6 +1,37 @@
 # some utility functions
 # written by Liam J. Revell 2011, 2012, 2013, 2014, 2015
 
+## compute the probability of states changes along edges of the tree
+## written by Liam J. Revell 2015
+edgeProbs<-function(trees){
+	SS<-sapply(trees,getStates,"tips")
+	states<-sort(unique(as.vector(SS)))
+	m<-length(states)
+	TT<-sapply(states,function(x,y) sapply(y,paste,x,sep="->"),
+		y=states)
+	nn<-c(TT[upper.tri(TT)],TT[lower.tri(TT)])
+	## this function computes for a given edge
+	fn<-function(edge,trees,states){
+		obj<-sapply(trees,function(x,e,s) 
+		if(names(x$maps[[e]])[1]==
+			s[1]&&names(x$maps[[e]])[length(x$maps[[e]])]==s[2]) TRUE
+		else FALSE,e=edge,s=states)
+		sum(obj)/length(obj)
+	}
+	edge.probs<-matrix(0,nrow(trees[[1]]$edge),m,dimnames=list(NULL,nn))
+	k<-1
+	for(i in 1:m) for(j in 1:m){
+		if(i!=j){ 
+			edge.probs[,k]<-sapply(1:nrow(trees[[1]]$edge),fn,
+				trees=trees,states=states[c(i,j)])
+			k<-k+1
+		}
+	}
+	edge.probs<-cbind(edge.probs,1-rowSums(edge.probs))
+	colnames(edge.probs)[ncol(edge.probs)]<-"no change"
+	edge.probs
+}
+
 ## get a position in the tree interactively
 ## written by Liam J. Revell 2015
 get.treepos<-function(message=TRUE){

man/edgeProbs.Rd

@@ -0,0 +1,24 @@
+\name{edgeProbs}
+\alias{edgeProbs}
+\title{Compute the relative frequencies of state changes along edges}
+\usage{
+edgeProbs(trees)
+}
+\arguments{
+	\item{trees}{an object of class \code{"multiSimmap"} containing a sample of trees that are identical in topology & branch lengths with different stochastically mapped character histories.}
+}
+\description{
+	This function computes the relative frequencies of character state changes along edges from a sample of stochastically mapped character histories. This function assumes that all trees in the sample differ only in their mapped histories & not at all in topology or branch lengths. Note that it only asks whether the starting and ending states of the edge differ in a particular way, and thus ignores multiple-hits along a single edge.
+}
+\value{
+	The object that is returned is a matrix with the state changes & the relative frequency of each state change. Rows are in the order of the matrix \code{edge} for any of the mapped trees.
+}
+\references{
+	Revell, L. J. (2012) phytools: An R package for phylogenetic comparative biology (and other things). \emph{Methods Ecol. Evol.}, \bold{3}, 217-223.
+}
+\author{Liam Revell \email{liam.revell@umb.edu}}
+\seealso{
+	\code{\link{plotSimmap}}
+}
+\keyword{phylogenetics}
+\keyword{plotting}