# function creates a stochastic character mapped tree as a modified "phylo" object # written by Liam Revell 2011 make.simmap<-function(tree,x,model="SYM"){ # check dependencies if(!require(ape)) warning("install ape") # fit discrete character model XX<-ace(x,tree,type="discrete",model=model) # get Q matrix Q<-matrix(XX$rates[XX$index.matrix],nrow(XX$index.matrix),ncol(XX$index.matrix),dimnames=list(colnames(XX$lik.anc),colnames(XX$lik.anc))) L<-XX$lik.anc; rownames(L)<-length(tree$tip)+1:nrow(L) diag(Q)<--colSums(Q,na.rm=TRUE) # create the map tree object mtree<-tree; mtree$maps<-list() # now we want to simulate the node states on the tree node.states<-matrix(NA,nrow(tree$edge),ncol(tree$edge)) node.states[which(tree$edge[,1]==(length(tree$tip)+1)),1]<-rstate(L[1,]) for(i in 1:nrow(tree$edge)){ if(tree$edge[i,2]>length(tree$tip)){ p<-expm(Q*tree$edge.length[i])[node.states[i,1],]*L[as.character(tree$edge[i,2]),] node.states[i,2]<-rstate(p/max(p)) node.states[which(tree$edge[,1]==tree$edge[i,2]),1]<-node.states[i,2] } else { node.states[i,2]<-x[tree$tip.label[tree$edge[i,2]]] } # now simulate on the branches accept=FALSE while(!accept){ time=0; state<-node.states[i,1]; new.state<-state; dt<-vector(); map<-vector(); j<-1 while(timecumy[j]){ state=names(y)[j+1]; j<-j+1 } return(state) } } # wraps around MatrixExp # written by Liam Revell 2011 expm<-function(Y){ Z<-MatrixExp(Y); dimnames(Z)<-dimnames(Y) return(Z) }