/in.progress/auteurExtended/R/simulation.R

  1
  2# MAIN FUNCTION for LEVY PROCESS
  3levyevolver=function(t, alpha, sigmasq.brown, sigma.jump, lambda.jump){
  4	y<-B<-rnorm(1, mean=alpha, sd=sqrt(sigmasq.brown*t))
  5	jumps=rpois(1, lambda.jump*t)
  6	if(jumps) {
  7		y=y+(L<-sum(rnorm(jumps, mean=y, sd=sigma.jump)))
  8	} else {
  9		L=alpha
 10	}
 11	return(list(y=y, Brownian.effect=B-alpha, Levy.effect=L-alpha))
 12}
 13
 14phenogram<-function(phy,model=c("levy","bm"),alpha=0,rate=0.01,lambda=0.5,sigma=0.5,increments=100,plot=TRUE,node.cex=2,...){
 15# 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)
 16	
 17# plotting
 18	add.transparency <-
 19	function (col, alpha) {
 20		tmp <- col2rgb(col)/255
 21		rgb(tmp[1, ], tmp[2, ], tmp[3, ], alpha = alpha)
 22	}
 23	
 24# find segments of time within a vector 
 25	get.lengths=function(yy) {
 26		yy[2:length(yy)]-yy[1:(length(yy)-1)]
 27	}
 28
 29	
 30	tree=phy 
 31	bb=branching.times(tree)
 32	abstimes=max(bb)-bb
 33	if(!is.null(increments)) increment=max(bb)/increments else increment=NULL
 34	
 35	n<-length(tree$tip.label)
 36	
 37## first simulate changes along each branch
 38	if(model=="levy") {
 39		## ENSURE THAT DRAW FROM LEVY PROCESS IS ACCURATELY MODELLED ##
 40				# evolves trait by Brownian motion as usual
 41				# adds additional deviate according to Poisson-distributed jumps and a given jump size
 42		x<-lapply(1:length(tree$edge.length), function(x) {
 43				  if(is.numeric(increment)) y=seq(0,tree$edge.length[x],by=increment) else y=0
 44				  if(length(y)==1) {
 45#					return(sample(c(-1,1),1)*rlevy(n=1, m=0, s=sqrt(rate*tree$edge.length[x])))
 46					return(levyevolver(t=tree$edge.length[x], alpha=alpha, sigmasq.brown=rate, sigma.jump=sigma, lambda.jump=lambda)$y)
 47				  } else {
 48#					return(sapply(get.lengths(y), function(z) return(sample(c(-1,1),1)*rlevy(n=1, m=0, s=sqrt(rate*z)))))
 49					return(sapply(get.lengths(y), function(z) return(levyevolver(t=z, alpha=alpha, sigmasq.brown=rate, sigma.jump=sigma, lambda.jump=lambda)$y)))
 50
 51				  }
 52			}
 53		)
 54		plot.title="Levy process"
 55	} else if(model=="bm") {
 56		x<-lapply(1:length(tree$edge.length), function(x) {
 57				  if(is.numeric(increment)) y=seq(0,tree$edge.length[x],by=increment) else y=0
 58				  if(length(y)==1) {
 59					return(rnorm(1, mean=0, sd=sqrt(rate*tree$edge.length[x])))
 60				  } else {
 61					return(sapply(get.lengths(y), function(z) rnorm(n=1,mean=0,sd=sqrt(rate*z))))
 62				  }
 63			}
 64		)		
 65		plot.title="Brownian"
 66	}
 67	
 68	node.phenotypes=matrix(cbind(tree$edge,0,max(bb)),ncol=4)
 69	for(i in 1:length(abstimes)) {
 70		if(names(abstimes)[i]%in%node.phenotypes[,2]) node.phenotypes[which(node.phenotypes[,2]==names(abstimes)[i]),4]=abstimes[i]
 71	}
 72	
 73## accumulate change along each path
 74	traithistory=list()
 75	times=list()
 76	y<-matrix(0,nrow(tree$edge),ncol(tree$edge)+1)
 77	for(i in 1:length(x)){
 78		yy=c() # phenotypes
 79		ntt=c() # node times
 80		tt=seq(0,tree$edge.length[i],length=length(x[[i]])+1)[-1] # each time increment along a branch
 81		if(length(tt)==0) tt=tree$edge.length[i] # time where branch is smaller than increment
 82		if(node.phenotypes[i,1]==n+1) { # ancestor
 83			start=alpha 
 84			start.time=0
 85		} else {						# non-ancestor
 86			start=node.phenotypes[which(node.phenotypes[,2]==node.phenotypes[i,1]),3]
 87			start.time=node.phenotypes[which(node.phenotypes[,2]==node.phenotypes[i,1]),4]
 88		}
 89		for(j in 1:length(x[[i]])) {
 90			ntt[j]=start.time+tt[j]
 91			if(j==1) {
 92				yy[j]=start+x[[i]][j] 
 93			} else {
 94				yy[j]=yy[j-1]+x[[i]][j]
 95			}
 96			if(j==length(x[[i]])) node.phenotypes[i,3]=yy[length(yy)]
 97		}
 98		traithistory[[i]]=c(start,yy)
 99		times[[i]]=c(start.time,ntt)
100	}
101	
102	mm=max(abs(alpha-unlist(traithistory)))
103	
104	# deal with non-ultrametric trees
105	for(i in 1:nrow(node.phenotypes)) {
106		tip=node.phenotypes[i,2]
107		if(!tip%in%tree$edge[,1]) {
108			last.edge=node.phenotypes[match(node.phenotypes[i,1],node.phenotypes[,2]),4]
109			if(is.na(last.edge)) last.edge=0
110			node.phenotypes[i,4]=last.edge+tree$edge.length[which(tree$edge[,2]==tip)]
111			
112		}
113	}
114	
115	
116	
117	## PLOTTING ##
118	if(plot) {
119		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=" "))
120		for(i in 1:length(x)) {
121			lines(times[[i]],traithistory[[i]],col=add.transparency("gray25",0.75),...)		
122		}
123		for(i in 1:length(x)) {
124			points(node.phenotypes[i,4],node.phenotypes[i,3],bg=add.transparency("white",0.75),pch=21,cex=node.cex)
125			
126		}
127		points(0,alpha,bg=add.transparency("white",0.75),pch=21,cex=node.cex)		
128	}
129	pp=data.frame(node.phenotypes)
130	names(pp)=c("ancestor","descendant","phenotype","time")
131	
132	
133	tt=which(pp$descendant<=Ntip(phy))
134	tips.dat=pp$phenotype[tt]
135	names(tips.dat)=phy$tip.label[pp$descendant[tt]]
136	
137	return(list(dat=tips.dat, phy=tree, phenotypes=pp))
138}