# function simulates stochastic character history under some model # written by Liam J. Revell 2011 sim.history<-function(tree,Q,anc=NULL,nsim=1){ # check dependencies if(!require(ape)) warning("install ape") # reorder to cladewise tree<-reorder.phylo(tree,"cladewise") # check Q if(!all(round(colSums(Q),8)==0)){ print("some columns of Q don't sum to 0.0. fixing. . .") diag(Q)<-0 diag(Q)<--colSums(Q,na.rm=TRUE) } # does Q have names? if(is.null(dimnames(Q))) dimnames(Q)<-list(1:nrow(Q),1:ncol(Q)) # create "multiPhylo" object mtrees<-list(); class(mtrees)<-"multiPhylo" # now loop for(i in 1:nsim){ # set root state if(is.null(anc)){ temp<-rep(1/nrow(Q),nrow(Q)); names(temp)<-rownames(Q) anc<-rstate(temp) } # 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]<-anc for(j in 1:nrow(tree$edge)){ if(tree$edge.length[j]==0){ map<-vector(); map[1]<-tree$edge.length[j] names(map)[1]<-node.states[which(tree$edge[,1]==tree$edge[j,2]),1]<-node.states[j,2]<-node.states[j,1] } else { time=0; state<-node.states[j,1]; new.state<-state; dt<-vector(); map<-vector(); k<-1 while(timenode.states[j,2]->node.states[which(tree$edge[,1]==tree$edge[j,2]),1] } mtree$maps[[j]]<-map } # add a couple of elements mtree$node.states<-node.states tip.states<-node.states[tree$edge[,2]<=length(tree$tip),2] tip.states<-tip.states[order(tree$edge[tree$edge[,2]<=length(tree$tip),2])] names(tip.states)<-tree$tip.label mtree$states<-tip.states # now construct the matrix "mapped.edge" (for backward compatibility allstates<-vector() for(j in 1:nrow(mtree$edge)) allstates<-c(allstates,names(mtree$maps[[j]])) allstates<-unique(allstates) mtree$mapped.edge<-matrix(data=0,length(mtree$edge.length),length(allstates),dimnames=list(apply(mtree$edge,1,function(x) paste(x,collapse=",")),state=allstates)) for(j in 1:length(mtree$maps)) for(k in 1:length(mtree$maps[[j]])) mtree$mapped.edge[j,names(mtree$maps[[j]])[k]]<-mtree$mapped.edge[j,names(mtree$maps[[j]])[k]]+mtree$maps[[j]][k] # plotSimmap(mtree) mtrees[[i]]<-mtree } if(nsim==1) mtrees<-mtrees[[1]] return(mtrees) } # function returns random state with probability given by y # written by Liam Revell 2011 rstate<-function(y){ if(length(y)==1) return(names(y)[1]) else { cumy<-y; for(i in 2:length(y)) cumy[i]<-cumy[i-1]+y[i] r<-runif(n=1); state=names(y)[1] j=1; while(r>cumy[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) }