/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

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