/in.progress/auteurExtended/R/simulation.R
http://github.com/eastman/auteur · R · 138 lines · 104 code · 21 blank · 13 comment · 31 complexity · e3fceca9d993bef20a89adc39d06d52f MD5 · raw file
- # MAIN FUNCTION for LEVY PROCESS
- levyevolver=function(t, alpha, sigmasq.brown, sigma.jump, lambda.jump){
- y<-B<-rnorm(1, mean=alpha, sd=sqrt(sigmasq.brown*t))
- jumps=rpois(1, lambda.jump*t)
- if(jumps) {
- y=y+(L<-sum(rnorm(jumps, mean=y, sd=sigma.jump)))
- } else {
- L=alpha
- }
- return(list(y=y, Brownian.effect=B-alpha, Levy.effect=L-alpha))
- }
- phenogram<-function(phy,model=c("levy","bm"),alpha=0,rate=0.01,lambda=0.5,sigma=0.5,increments=100,plot=TRUE,node.cex=2,...){
- # written by JM Eastman 02.22.2011 as an extension of fastBM() by LJ Revell, under a Levy process (as a compound Brownian-Poisson process)
-
- # plotting
- add.transparency <-
- function (col, alpha) {
- tmp <- col2rgb(col)/255
- rgb(tmp[1, ], tmp[2, ], tmp[3, ], alpha = alpha)
- }
-
- # find segments of time within a vector
- get.lengths=function(yy) {
- yy[2:length(yy)]-yy[1:(length(yy)-1)]
- }
-
- tree=phy
- bb=branching.times(tree)
- abstimes=max(bb)-bb
- if(!is.null(increments)) increment=max(bb)/increments else increment=NULL
-
- n<-length(tree$tip.label)
-
- ## first simulate changes along each branch
- if(model=="levy") {
- ## ENSURE THAT DRAW FROM LEVY PROCESS IS ACCURATELY MODELLED ##
- # evolves trait by Brownian motion as usual
- # adds additional deviate according to Poisson-distributed jumps and a given jump size
- x<-lapply(1:length(tree$edge.length), function(x) {
- if(is.numeric(increment)) y=seq(0,tree$edge.length[x],by=increment) else y=0
- if(length(y)==1) {
- # return(sample(c(-1,1),1)*rlevy(n=1, m=0, s=sqrt(rate*tree$edge.length[x])))
- return(levyevolver(t=tree$edge.length[x], alpha=alpha, sigmasq.brown=rate, sigma.jump=sigma, lambda.jump=lambda)$y)
- } else {
- # return(sapply(get.lengths(y), function(z) return(sample(c(-1,1),1)*rlevy(n=1, m=0, s=sqrt(rate*z)))))
- return(sapply(get.lengths(y), function(z) return(levyevolver(t=z, alpha=alpha, sigmasq.brown=rate, sigma.jump=sigma, lambda.jump=lambda)$y)))
- }
- }
- )
- plot.title="Levy process"
- } else if(model=="bm") {
- x<-lapply(1:length(tree$edge.length), function(x) {
- if(is.numeric(increment)) y=seq(0,tree$edge.length[x],by=increment) else y=0
- if(length(y)==1) {
- return(rnorm(1, mean=0, sd=sqrt(rate*tree$edge.length[x])))
- } else {
- return(sapply(get.lengths(y), function(z) rnorm(n=1,mean=0,sd=sqrt(rate*z))))
- }
- }
- )
- plot.title="Brownian"
- }
-
- node.phenotypes=matrix(cbind(tree$edge,0,max(bb)),ncol=4)
- for(i in 1:length(abstimes)) {
- if(names(abstimes)[i]%in%node.phenotypes[,2]) node.phenotypes[which(node.phenotypes[,2]==names(abstimes)[i]),4]=abstimes[i]
- }
-
- ## accumulate change along each path
- traithistory=list()
- times=list()
- y<-matrix(0,nrow(tree$edge),ncol(tree$edge)+1)
- for(i in 1:length(x)){
- yy=c() # phenotypes
- ntt=c() # node times
- tt=seq(0,tree$edge.length[i],length=length(x[[i]])+1)[-1] # each time increment along a branch
- if(length(tt)==0) tt=tree$edge.length[i] # time where branch is smaller than increment
- if(node.phenotypes[i,1]==n+1) { # ancestor
- start=alpha
- start.time=0
- } else { # non-ancestor
- start=node.phenotypes[which(node.phenotypes[,2]==node.phenotypes[i,1]),3]
- start.time=node.phenotypes[which(node.phenotypes[,2]==node.phenotypes[i,1]),4]
- }
- for(j in 1:length(x[[i]])) {
- ntt[j]=start.time+tt[j]
- if(j==1) {
- yy[j]=start+x[[i]][j]
- } else {
- yy[j]=yy[j-1]+x[[i]][j]
- }
- if(j==length(x[[i]])) node.phenotypes[i,3]=yy[length(yy)]
- }
- traithistory[[i]]=c(start,yy)
- times[[i]]=c(start.time,ntt)
- }
-
- mm=max(abs(alpha-unlist(traithistory)))
-
- # deal with non-ultrametric trees
- for(i in 1:nrow(node.phenotypes)) {
- tip=node.phenotypes[i,2]
- if(!tip%in%tree$edge[,1]) {
- last.edge=node.phenotypes[match(node.phenotypes[i,1],node.phenotypes[,2]),4]
- if(is.na(last.edge)) last.edge=0
- node.phenotypes[i,4]=last.edge+tree$edge.length[which(tree$edge[,2]==tip)]
-
- }
- }
-
-
-
- ## PLOTTING ##
- if(plot) {
- plot(x=NULL, y=NULL, xlim=c(0,max(node.phenotypes[,4])), ylim=c(-2*mm+alpha, 2*mm+alpha), bty="n", xlab="time", ylab="phenotypic value",main=paste("evolution by",plot.title,sep=" "))
- for(i in 1:length(x)) {
- lines(times[[i]],traithistory[[i]],col=add.transparency("gray25",0.75),...)
- }
- for(i in 1:length(x)) {
- points(node.phenotypes[i,4],node.phenotypes[i,3],bg=add.transparency("white",0.75),pch=21,cex=node.cex)
-
- }
- points(0,alpha,bg=add.transparency("white",0.75),pch=21,cex=node.cex)
- }
- pp=data.frame(node.phenotypes)
- names(pp)=c("ancestor","descendant","phenotype","time")
-
-
- tt=which(pp$descendant<=Ntip(phy))
- tips.dat=pp$phenotype[tt]
- names(tips.dat)=phy$tip.label[pp$descendant[tt]]
-
- return(list(dat=tips.dat, phy=tree, phenotypes=pp))
- }