/archive/source.archive/bm.and.ou/trait.rjmcmc.10022010.R

http://github.com/eastman/auteur · R · 679 lines · 583 code · 69 blank · 27 comment · 105 complexity · 079f4e060120ba67d63151052ee5e759 MD5 · raw file

  1. ## Written by LJ HARMON & AH: 2008
  2. ## Modified by JM EASTMAN: 09.2010
  3. rjMCMC.trait<-function (phy, data, ngen=1000, sampleFreq=100, model=c("OU","BM"), probAlpha=0.05, probMergeSplit=0.05, probRoot=0.1, theta.rates.ratio=1, heat=1, parallel=FALSE, fileBase="result")
  4. {
  5. if(length(model)>1)stop(paste("Please specify either ", sQuote("OU"), " or ", sQuote("BM")))
  6. ## MORE ERROR CHECKING NEEDED HERE ##
  7. if(parallel) require(multicore)
  8. require(ouch)
  9. require(geiger)
  10. ### prepare objects for rjMCMC
  11. cur.model <- model
  12. dataList <- prepare.ouchdata(phy, data)
  13. ape.tre <- cur.tre <- dataList$ape.tre
  14. orig.dat <- dataList$orig.dat
  15. ouch.tre <- dataList$ouch.tre
  16. ouch.dat <- dataList$ouch.dat
  17. nn <- length(ape.tre$edge.length)
  18. node.des <- sapply(unique(c(ape.tre$edge[1,1],ape.tre$edge[,2])), function(x) get.descendants.of.node(x, ape.tre))
  19. names(node.des) <- c(ape.tre$edge[1,1], unique(ape.tre$edge[,2]))
  20. ll <- list()
  21. # initialize parameters
  22. init.rate <- generate.starting.point(data,ape.tre,node.des,theta=FALSE)
  23. cur.rates <- init.rate$values
  24. cur.delta.rates <- init.rate$delta
  25. cur.root <- adjust.value(mean(orig.dat))
  26. cur.alpha <- adjust.rate(0.001)
  27. init.theta <- generate.starting.point(data,ape.tre,node.des,theta=TRUE)
  28. cur.theta <- init.theta$values
  29. cur.delta.theta <- init.theta$delta
  30. cur.root.theta <- assign.root.theta(cur.theta,ape.tre)
  31. cur.regimes <- as.factor(c(cur.root.theta,cur.theta))
  32. cur.tre <- apply.BMM.to.apetree(ape.tre, cur.rates)
  33. if(model=="OU") {
  34. mod.cur = new.ou.lik.fn(cur.rates, cur.alpha, cur.regimes, cur.theta, cur.tre, ouch.dat)
  35. } else {
  36. mod.cur = new.bm.lik.fn(cur.rates, cur.root, cur.tre, orig.dat)
  37. }
  38. # proposal counts
  39. nRateProp = 0
  40. nRateSwapProp = 0
  41. nRcatProp = 0
  42. nRootProp = 0
  43. nAlphaProp = 0
  44. nTcatProp = 0
  45. nThetaProp = 0
  46. nThetaSwapProp= 0
  47. nRateOK = 0
  48. nRateSwapOK = 0
  49. nRcatOK = 0
  50. nRootOK = 0
  51. nAlphaOK = 0
  52. nTcatOK = 0
  53. nThetaOK = 0
  54. nThetaSwapOK = 0
  55. cOK=c(
  56. nRateOK,
  57. nRateSwapOK,
  58. nRcatOK,
  59. nRootOK,
  60. nAlphaOK,
  61. nThetaOK,
  62. nThetaSwapOK,
  63. nTcatOK
  64. )
  65. tickerFreq=ceiling(ngen/30)
  66. # file handling
  67. errorLog=file(paste(getwd(),paste(cur.model, fileBase, "rjTraitErrors.log",sep="."),sep="/"),open='w+')
  68. generate.error.message(i=NULL, mod.cur=NULL, mod.new=NULL, lnR=NULL, errorLog=errorLog, init=TRUE)
  69. runLog=file(paste(getwd(),paste(cur.model, fileBase, "rjTrait.log",sep="."),sep="/"),open='w+')
  70. generate.log(bundled.parms=NULL, cur.model, file=runLog, init=TRUE)
  71. ### Begin rjMCMC
  72. for (i in 1:ngen) {
  73. lnLikelihoodRatio <- lnHastingsRatio <- lnPriorRatio <- 0
  74. startparms = c(nRateProp, nRateSwapProp, nRcatProp, nRootProp, nAlphaProp, nThetaProp, nThetaSwapProp, nTcatProp)
  75. ## OU IMPLEMENTATION ##
  76. if(cur.model=="OU") {
  77. if (runif(1) < (2 * probMergeSplit)) {
  78. if(runif(1)<(theta.rates.ratio/(theta.rates.ratio+1))) { # adjust theta categories
  79. nr=split.or.merge(cur.delta.theta, cur.theta, ape.tre, node.des, theta=TRUE)
  80. new.alpha=cur.alpha
  81. new.theta=nr$new.values ##
  82. new.delta.theta=nr$new.delta ##
  83. new.root.theta<-nrt<-assign.root.theta(new.theta, ape.tre) ##
  84. new.regimes=as.factor(c(nrt,new.theta)) ##
  85. new.rates=cur.rates
  86. new.delta.rates=cur.delta.rates
  87. nTcatProp=nTcatProp+1
  88. } else { # adjust rate categories
  89. nr=split.or.merge(cur.delta.rates, cur.rates, ape.tre, node.des, theta=FALSE)
  90. new.alpha=cur.alpha
  91. new.root.theta=cur.root.theta
  92. new.theta=cur.theta
  93. new.delta.theta=cur.delta.theta
  94. new.regimes=cur.regimes
  95. new.rates=nr$new.values ##
  96. new.delta.rates=nr$new.delta ##
  97. nRcatProp=nRcatProp+1
  98. }
  99. lnHastingsRatio=nr$lnHastingsRatio
  100. lnPriorRatio=nr$lnPriorRatio
  101. } else { # neither split nor merge
  102. if(runif(1)<probAlpha) { # adjust alpha
  103. new.alpha=adjust.rate(cur.alpha) ##
  104. new.root.theta=cur.root.theta
  105. new.theta=cur.theta
  106. new.delta.theta=cur.delta.theta
  107. new.regimes=cur.regimes
  108. new.rates=cur.rates
  109. new.delta.rates=cur.delta.rates
  110. nAlphaProp=nAlphaProp+1
  111. } else {
  112. if(runif(1)>(theta.rates.ratio/(theta.rates.ratio+1))) {
  113. if(runif(1)<0.2 | sum(cur.delta.rates)==0) {
  114. new.alpha=cur.alpha # adjust rates
  115. new.root.theta=cur.root.theta
  116. new.theta=cur.theta
  117. new.delta.theta=cur.delta.theta
  118. new.regimes=cur.regimes
  119. new.rates=adjust.rate(cur.rates) ##
  120. new.delta.rates=cur.delta.rates
  121. nRateProp=nRateProp+1
  122. } else { # swap rates
  123. nr=switcheroo(cur.delta.rates, cur.rates, ape.tre, node.des, theta=FALSE)
  124. new.alpha=cur.alpha
  125. new.root.theta=cur.root.theta
  126. new.theta=cur.theta
  127. new.delta.theta=cur.delta.theta
  128. new.regimes=cur.regimes
  129. new.rates=nr$new.values ##
  130. new.delta.rates=nr$new.delta ##
  131. nRateSwapProp=nRateSwapProp+1
  132. }
  133. } else {
  134. if(runif(1)<0.2 | sum(cur.delta.theta)==0) { # adjust theta
  135. new.alpha=cur.alpha
  136. new.theta=adjust.value(cur.theta) ##
  137. new.delta.theta=cur.delta.theta
  138. new.root.theta<-nrt<-assign.root.theta(new.theta, ape.tre) ##
  139. new.regimes=as.factor(c(nrt,new.theta)) ##
  140. new.rates=cur.rates
  141. new.delta.rates=cur.delta.rates
  142. nThetaProp=nThetaProp+1
  143. } else { # swap theta
  144. nr=switcheroo(cur.delta.theta, cur.theta, ape.tre, node.des, theta=TRUE)
  145. new.alpha=cur.alpha
  146. new.theta=nr$new.values ##
  147. new.delta.theta=nr$new.delta ##
  148. new.root.theta<-nrt<-assign.root.theta(new.theta, ape.tre) ##
  149. new.regimes=as.factor(c(nrt,new.theta)) ##
  150. new.rates=cur.rates
  151. new.delta.rates=cur.delta.rates
  152. nThetaSwapProp=nThetaSwapProp+1
  153. }
  154. }
  155. }
  156. }
  157. # cur.tre=apply.BMM.to.apetree(ape.tre, cur.rates)
  158. new.tre=apply.BMM.to.apetree(ape.tre, new.rates)
  159. if(!parallel) {
  160. # mod.cur=new.ou.lik.fn(cur.rates, cur.alpha, cur.regimes, cur.theta, cur.tre, ouch.dat)
  161. mod.new=new.ou.lik.fn(new.rates, new.alpha, new.regimes, new.theta, new.tre, ouch.dat)
  162. } else {
  163. mod.cur=list(cur.rates, cur.alpha, cur.regimes, cur.theta, cur.tre, ouch.dat)
  164. mod.new=list(new.rates, new.alpha, new.regimes, new.theta, new.tre, ouch.dat)
  165. res=mclapply(list(mod.cur, mod.new), function(x) new.ou.lik.fn(x[[1]], x[[2]], x[[3]], x[[4]], x[[5]], x[[6]]))
  166. res[[1]]->mod.cur
  167. res[[2]]->mod.new
  168. }
  169. } else if(cur.model=="BM"){
  170. ## BM IMPLEMENTATION ##
  171. if (runif(1) < (2 * probMergeSplit)) { # adjust rate categories
  172. nr=split.or.merge(cur.delta.rates, cur.rates, ape.tre, node.des, theta=FALSE)
  173. new.rates=nr$new.values
  174. new.delta.rates=nr$new.delta
  175. new.root=cur.root
  176. nRcatProp=nRcatProp+1
  177. lnHastingsRatio=nr$lnHastingsRatio
  178. lnPriorRatio=nr$lnPriorRatio
  179. } else {
  180. if(runif(1)<probRoot) { # adjust root
  181. new.root=adjust.value(cur.root)
  182. new.rates=cur.rates
  183. new.delta.rates=cur.delta.rates
  184. nRootProp=nRootProp+1
  185. } else { # adjust rates
  186. if(runif(1)<0.05 | sum(cur.delta.rates)==0) {
  187. new.rates=adjust.rate(cur.rates)
  188. new.delta.rates=cur.delta.rates
  189. new.root=cur.root
  190. nRateProp=nRateProp+1
  191. } else { # swap rate classes
  192. nr=switcheroo(cur.delta.rates, cur.rates, ape.tre, node.des, theta=FALSE)
  193. new.rates=nr$new.values ##
  194. new.delta.rates=nr$new.delta ##
  195. new.root=cur.root
  196. nRateSwapProp=nRateSwapProp+1
  197. }
  198. }
  199. }
  200. # cur.tre=apply.BMM.to.apetree(ape.tre, cur.rates)
  201. new.tre=apply.BMM.to.apetree(ape.tre, new.rates)
  202. if(!parallel) {
  203. # mod.cur=new.bm.lik.fn(cur.rates, cur.root, cur.tre, orig.dat)
  204. mod.new=new.bm.lik.fn(new.rates, new.root, new.tre, orig.dat)
  205. } else {
  206. mod.cur=list(cur.rates, cur.root, cur.tre, orig.dat)
  207. mod.new=list(new.rates, new.root, new.tre, orig.dat)
  208. res=mclapply(list(mod.cur, mod.new), function(x) new.bm.lik.fn(x[[1]], x[[2]], x[[3]], x[[4]]))
  209. res[[1]]->mod.cur
  210. res[[2]]->mod.new
  211. }
  212. if(is.infinite(mod.new$lnL)) {
  213. failures="./failed.model.parameters/"
  214. failed.trees=paste(failures, "trees", sep="/")
  215. lapply(list(failures, failed.trees), function(x) if(!file.exists(x)) dir.create(x))
  216. pdf(file=paste(failed.trees, paste(i,"failed.model.pdf",sep="."),sep="/"))
  217. plot(new.tre)
  218. dev.off()
  219. write.table(i,file=paste(failures, "failed.rates.txt",sep="/"),append=TRUE,col=FALSE,row=FALSE,quote=FALSE)
  220. write.table(matrix(new.rates,nrow=1),file=paste(failures,"failed.rates.txt",sep="/"),append=TRUE,col=FALSE,row=FALSE,quote=FALSE)
  221. write.table(matrix(c(i, new.root),nrow=1),file=paste(failures,"failed.root.txt",sep="/"),append=TRUE,col=FALSE,row=FALSE,quote=FALSE)
  222. }
  223. }
  224. endparms = c(nRateProp, nRateSwapProp, nRcatProp, nRootProp, nAlphaProp, nThetaProp, nThetaSwapProp, nTcatProp)
  225. if(!inherits(mod.new, "try-error")) {
  226. lnLikelihoodRatio = mod.new$lnL - mod.cur$lnL
  227. } else {
  228. lnLikelihoodRatio = -Inf
  229. mod.new$lnL=NaN
  230. }
  231. r=assess.lnR((heat * lnLikelihoodRatio + heat * lnPriorRatio + lnHastingsRatio)->lnR)
  232. if (runif(1) <= r$r) { # adopt proposal
  233. if(cur.model=="OU") {
  234. decision="adopt"
  235. cur.alpha <- new.alpha
  236. cur.root.theta <- new.root.theta
  237. cur.regimes <- new.regimes
  238. cur.theta <- new.theta
  239. cur.delta.theta <- new.delta.theta
  240. } else {
  241. decision="adopt"
  242. cur.root <- new.root
  243. }
  244. cur.rates <- new.rates
  245. cur.delta.rates <- new.delta.rates
  246. curr.lnL <- mod.new$lnL
  247. cur.tre <- new.tre
  248. cOK=determine.accepted.proposal(startparms, endparms, cOK)
  249. } else { # deny proposal
  250. decision="reject"
  251. curr.lnL <- mod.cur$lnL
  252. }
  253. if(is.infinite(curr.lnL)) stop("starting point has exceptionally poor likelihood")
  254. if(r$error) generate.error.message(i, mod.cur, mod.new, lnR, errorLog)
  255. if(i%%tickerFreq==0) {
  256. if(i==tickerFreq) cat("|",rep(" ",9),toupper("generations complete"),rep(" ",9),"|","\n")
  257. cat(". ")
  258. }
  259. if(i%%sampleFreq==0) {
  260. bundled.parms=list(gen=i, mrate=exp(mean(log(cur.rates))), cats=sum(cur.delta.rates)+1, root=cur.root, alpha=cur.alpha, reg=(sum(cur.delta.theta)+1), theta=mean(cur.theta), lnL=curr.lnL)
  261. generate.log(bundled.parms, cur.model, file=runLog)
  262. ll$lnL=rbind(ll$lnL, curr.lnL)
  263. ll$rates=rbind(ll$rates, cur.rates)
  264. if(cur.model=="OU") {
  265. ll$theta=rbind(ll$theta, cur.theta)
  266. ll$alpha=rbind(ll$alpha, cur.alpha)
  267. ll$root=rbind(ll$root, cur.root.theta)
  268. }
  269. if(cur.model=="BM") {
  270. ll$root=rbind(ll$root, cur.root)
  271. }
  272. if(i==ngen) {
  273. if(cur.model=="OU") dimnames(ll$theta)[[2]]=ape.tre$edge[,2]
  274. dimnames(ll$rates)[[2]]=ape.tre$edge[,2]
  275. }
  276. }
  277. }
  278. # End rjMCMC
  279. close(errorLog)
  280. close(runLog)
  281. summarize.run(cOK, endparms, cur.model)
  282. return(list(results=ll,data=orig.dat,tree=ape.tre))
  283. }
  284. ## AUXILIARY FUNCTIONS
  285. new.ou.lik.fn <- function(rates, alpha, regimes, theta, ape.tre, ouch.dat) { # from OUCH 2.7-1:::hansen()
  286. # mod 09.13.2010 JM Eastman
  287. ouch.tre=ape2ouch(ape.tre, scale=FALSE)
  288. ouch.dat$regimes=as.factor(regimes)
  289. theta.tmp=as.numeric(levels(as.factor(theta)))
  290. theta=theta.tmp[order(match(theta.tmp,theta))]
  291. alpha.ouch=as.matrix(sqrt(alpha))
  292. beta=ouch:::regime.spec(ouch.tre, ouch.dat['regimes'])
  293. dat=do.call(c,lapply(ouch.dat['trait'],function(y)y[ouch.tre@term]))
  294. names(dat)=ouch.dat$labels[Ntip(ape.tre):nrow(ouch.dat)]
  295. n <- length(dat)
  296. otree.df=as(ouch.tre,"data.frame")
  297. ordering=otree.df[n:nrow(otree.df),'labels']
  298. ev <- eigen(alpha.ouch,symmetric=TRUE)
  299. w <- .Call(ouch:::ouch_weights,object=ouch.tre,lambda=ev$values,S=ev$vectors,beta=beta)
  300. v <- geiger:::ouMatrix(vcv.phylo(ape.tre), alpha)
  301. v.match <- match(row.names(v),ordering)
  302. v <- v[v.match,v.match]
  303. ## GEIGER call does not assume ultrametricity (and thus allows different 'rates' under BM)
  304. # -old- v <- .Call(ouch:::ouch_covar,object=ouch.tre,lambda=ev$values,S=ev$vectors,sigma.sq=1)
  305. y <- matrix(theta,ncol=1)
  306. e <- w%*%y-dat
  307. dim(y) <- dim(y)[1]
  308. dim(e) <- n
  309. q <- e%*%solve(v,e)
  310. det.v <- determinant(v,logarithm=TRUE)
  311. if (det.v$sign!=1)stop("ou.lik.fn error: non-positive determinant",call.=FALSE)
  312. dev <- n*log(2*pi)+as.numeric(det.v$modulus)+q[1,1]
  313. list(
  314. theta=theta,
  315. weight=w,
  316. vcov=v,
  317. resids=e,
  318. lnL=-0.5*dev
  319. )
  320. }
  321. new.bm.lik.fn <- function(rates, root, ape.tre, orig.dat) { # from OUCH 2.7-1:::brown()
  322. # mod 09.16.2010 JM Eastman
  323. tree=ape.tre
  324. data=orig.dat
  325. nterm=Ntip(tree)
  326. nchar=length(data)
  327. dat=data[match(tree$tip.label,names(data))]
  328. nchar <- 1
  329. w <- matrix(data=1,nrow=nterm,ncol=1)
  330. n <- length(dat)
  331. b <- vcv.phylo(ape.tre)
  332. y <- matrix(root)
  333. e <- w%*%y-dat
  334. q <- t(e)%*%solve(b,e)
  335. v <- q/nterm
  336. dev <- nchar*nterm*(1+log(2*pi))+nchar*log(det(b))+nterm*log(det(v))
  337. list(
  338. root=root,
  339. lnL = -0.5*dev
  340. )
  341. }
  342. split.or.merge <- function(cur.delta, cur.values, phy, node.des, theta=FALSE) {
  343. # mod 09.17.2010 JM Eastman
  344. bb=cur.delta
  345. vv=cur.values
  346. names(vv)<-names(bb)<-phy$edge[,2]
  347. new.bb=choose.one(bb)
  348. new.vv=vv
  349. s=names(new.bb[bb!=new.bb])
  350. all.shifts=as.numeric(names(bb[bb>0]))
  351. all.D=node.des[[which(names(node.des)==s)]]
  352. if(length(all.D)!=0) {
  353. untouchable=unlist(lapply(all.shifts[all.shifts>s], function(x)node.des[[which(names(node.des)==x)]]))
  354. remain.unchanged=union(all.shifts, untouchable)
  355. } else {
  356. untouchable=NULL
  357. remain.unchanged=list()
  358. }
  359. marker=match(s, names(new.vv))
  360. nn=length(vv)
  361. K=sum(bb)+1
  362. N=Ntip(phy)
  363. if(sum(new.bb)>sum(bb)) { # add transition: SPLIT
  364. # print("s")
  365. if(theta) { # assign new theta to current node
  366. new.vv[marker] <- nr <- adjust.value(new.vv[marker])
  367. } else { # assign new rate to current node
  368. new.vv[marker] <- nr <- adjust.rate(new.vv[marker])
  369. }
  370. if(length(remain.unchanged)==0) { # assign new value to all descendants
  371. new.vv[match(all.D,names(new.vv))] = nr
  372. } else { # assign new value to all 'open' descendants
  373. new.vv[match(all.D[!(all.D%in%remain.unchanged)],names(new.vv))] = nr
  374. }
  375. lnHastingsRatio = log((K+1)/(2*N-2-K)) ### from Drummond and Suchard 2010: where N is tips, K is number of local parms in tree
  376. lnPriorRatio = log(ppois(K+1,log(2),lower.tail=FALSE)/ppois(K,log(2),lower.tail=FALSE))
  377. } else { # drop transition: MERGE
  378. # print("m")
  379. anc = get.ancestor.of.node(s, phy)
  380. if(!is.root(anc, phy)) { # substitute ancestral value for current
  381. new.vv[match(s, names(new.vv))] <- nr <- new.vv[match(anc, names(new.vv))]
  382. } else {
  383. sister.tmp=get.desc.of.node(anc,phy)
  384. sister=sister.tmp[sister.tmp!=s]
  385. if(bb[match(sister, names(bb))]==0) {
  386. nr=as.numeric(vv[match(sister, names(vv))])
  387. } else {
  388. if(theta) { # assign new theta to first descendant of root
  389. new.vv[marker] <- nr <- adjust.value(new.vv[marker])
  390. } else { # assign new rate to first descendant of root
  391. new.vv[marker] <- nr <- adjust.rate(new.vv[marker])
  392. }
  393. }
  394. }
  395. if(length(remain.unchanged)==0) { # assign new value to all descendants
  396. new.vv[match(all.D, names(new.vv))] = nr
  397. } else { # assign new value to all 'open' descendants
  398. new.vv[match(all.D[!(all.D%in%remain.unchanged)],names(new.vv))] = nr
  399. }
  400. lnHastingsRatio = log((2*N-2-K+1)/K) ### from Drummond and Suchard 2010: where N is tips, K is number of local parms in tree
  401. lnPriorRatio = log(ppois(K-1,log(2),lower.tail=FALSE)/ppois(K,log(2),lower.tail=FALSE))
  402. }
  403. if(((length(table(new.vv))-1)-sum(new.bb))!=0) fail=TRUE else fail=FALSE
  404. return(list(new.delta=new.bb, new.values=new.vv, lnHastingsRatio=lnHastingsRatio, lnPriorRatio=lnPriorRatio, fail=fail))
  405. }
  406. switcheroo <- function(cur.delta, cur.values, phy, node.des, theta=FALSE) {
  407. bb=cur.delta
  408. vv=cur.values
  409. new.bb=sample(bb)
  410. names(vv)<-names(bb)<-names(new.bb)<-phy$edge[,2]
  411. new.vv=vv
  412. old.shifts=as.numeric(names(bb[bb!=0]))
  413. old.values=vv[match(old.shifts, names(vv))]
  414. anc=lapply(old.shifts, function(x)get.ancestor.of.node(x,phy))
  415. if(!is.root(min(unlist(anc)),phy)) {
  416. # find most ancestral value in tree
  417. anc.value=unname(vv[match(min(unlist(anc)),names(vv))])
  418. } else {
  419. if(theta) anc.value=mean(sapply(old.values, function(x) adjust.value(x))) else anc.value=mean(sapply(old.values, function(x) adjust.rate(x)))
  420. }
  421. all.desc.anc=node.des[[which(names(node.des)==min(unlist(anc)))]]
  422. new.vv[match(all.desc.anc,names(new.vv))] = anc.value
  423. new.shifts=sort(as.numeric(names(new.bb[new.bb!=0])))
  424. new.tips=new.shifts[new.shifts<=Ntip(phy)]
  425. new.ints=new.shifts[new.shifts>Ntip(phy)]
  426. new.shifts=c(new.ints, new.tips)
  427. new.values=as.numeric(old.values)[sample(1:length(old.values))]
  428. for(z in 1:length(new.shifts)) {
  429. node=new.shifts[z]
  430. new.vv[match(node, names(new.vv))]=new.values[z]
  431. new.vv[match(node.des[[which(names(node.des)==node)]],names(new.vv))]=new.values[z]
  432. }
  433. return(list(new.delta=new.bb, new.values=new.vv))
  434. }
  435. generate.starting.point <- function(data, phy, node.des, theta=FALSE) {
  436. nn=length(phy$edge.length)
  437. ntip=Ntip(phy)
  438. nshifts=rtpois(1,log(ntip),nn)
  439. if(nshifts!=0) bb=sample(c(rep(1,nshifts),rep(0,nn-nshifts)),replace=FALSE) else bb=rep(0,nn)
  440. names(bb)=phy$edge[,2]
  441. shifts=as.numeric(names(bb)[bb==1])
  442. if(theta) {
  443. min.max=c(adjust.value(min(data)), adjust.value(max(data)))
  444. values=runif(sum(bb)+1, min=min.max[min.max==min(min.max)], max=min.max[min.max==max(min.max)])
  445. } else {
  446. init.rate=fitContinuous(phy,data)[[1]]$beta
  447. min.max=c(adjust.rate(init.rate), adjust.rate(init.rate))
  448. values=runif(sum(bb)+1, min=min.max[min.max==min(min.max)], max=min.max[min.max==max(min.max)])
  449. }
  450. internal.shifts<-tip.shifts<-numeric(0)
  451. internal.shifts=sort(shifts[shifts>ntip])
  452. tip.shifts=shifts[shifts<=ntip]
  453. if(length(internal.shifts)==0 & length(tip.shifts)==0) {
  454. vv=rep(values, nn)
  455. } else {
  456. vv=bb
  457. vv[]=values[length(values)]
  458. i=0
  459. if(length(internal.shifts)!=0) {
  460. for(i in 1:length(internal.shifts)) {
  461. d=node.des[which(names(node.des)==internal.shifts[i])]
  462. # d=get.descendants.of.node(internal.shifts[i],phy)
  463. vv[match(c(internal.shifts[i], d), names(vv))]=values[i]
  464. }
  465. }
  466. if(length(tip.shifts)!=0) {
  467. for(j in 1:length(tip.shifts)) {
  468. vv[match(tip.shifts[j], names(vv))]=values[j+i]
  469. }
  470. }
  471. }
  472. return(list(delta=unname(bb), values=unname(vv)))
  473. }
  474. assign.root.theta <- function(cur.theta, phy) {
  475. root=phy$edge[1,1]
  476. desc=phy$edge[which(phy$edge[,1]==root),2]
  477. root.theta=cur.theta[sample(which(phy$edge[,2]==desc),1)]
  478. root.theta
  479. }
  480. get.descendants.of.node <- function(node, phy) {
  481. storage=list()
  482. if(node%in%phy$edge[,1]) {
  483. desc=c(sapply(node, function(x) {return(phy$edge[which(phy$edge[,1]==x),2])}))
  484. while(any(desc%in%phy$edge[,1])) {
  485. desc.true=desc[which(desc%in%phy$edge[,1])]
  486. desc.desc=lapply(desc.true, function(x) return(get.desc.of.node(x, phy)))
  487. storage=list(storage,union(desc,desc.desc))
  488. desc=unlist(desc.desc)
  489. }
  490. storage=sort(unique(unlist(union(desc,storage))))
  491. }
  492. return(storage)
  493. }
  494. get.desc.of.node <- function(node, phy) {
  495. return(phy$edge[which(phy$edge[,1]==node),2])
  496. }
  497. get.ancestor.of.node<-function(node, phy) {
  498. anc=phy$edge[which(phy$edge[,2]==node),1]
  499. anc
  500. }
  501. choose.one <- function(cur.delta)
  502. {
  503. bb=cur.delta
  504. s=sample(1:length(bb),1)
  505. bb[s]=1-bb[s]
  506. bb
  507. }
  508. is.root<-function(node,phy) {
  509. if(node==phy$edge[1,1]) return(TRUE) else return(FALSE)
  510. }
  511. assess.lnR <- function(lnR) {
  512. if(is.na(lnR) || abs(lnR)==Inf) {
  513. error=TRUE
  514. r=0
  515. } else {
  516. if(lnR < -20) {
  517. r=0
  518. } else if(lnR > 0) {
  519. r=1
  520. } else {
  521. r=exp(lnR)
  522. }
  523. error=FALSE
  524. }
  525. return(list(r=r, error=error))
  526. }
  527. apply.BMM.to.apetree<-function(apetree, rates) {
  528. apetree$edge.length=apetree$edge.length*rates
  529. return(apetree)
  530. }
  531. adjust.value <- function(value) {
  532. vv=value
  533. vv.levels=levels(factor(vv))
  534. rch <- runif(length(vv.levels), min = -1, max = 1)
  535. for(i in 1:length(vv.levels)){
  536. vv[which(vv==vv.levels[i])]=vv[which(vv==vv.levels[i])]+rch[i]
  537. }
  538. vv
  539. }
  540. adjust.rate <- function(rate) {
  541. vv=rate
  542. v=log(vv)
  543. vv.levels=levels(factor(vv))
  544. rch <- runif(length(vv.levels), min = -1, max = 1)
  545. for(i in 1:length(vv.levels)){
  546. v[which(vv==vv.levels[i])]=v[which(vv==vv.levels[i])]+rch[i]
  547. }
  548. v=exp(v)
  549. v
  550. }
  551. prepare.ouchdata <- function(phy, data) {
  552. td <- treedata(phy, data, sort = T)
  553. ouch.tre <- ape2ouch(td$phy->ape.tre, scale=FALSE)
  554. ouch.dat=data.frame(as.numeric(ouch.tre@nodes), as.character(ouch.tre@nodelabels))
  555. names(ouch.dat)=c("nodes","labels")
  556. ouch.dat[ouch.dat[,2]=="",2]=NA
  557. ouch.dat$trait=NA
  558. mM=match(ouch.dat$labels,names(data))
  559. for(i in 1:nrow(ouch.dat)){
  560. if(!is.na(mM[i])){
  561. ouch.dat$trait[i]=data[mM[i]]
  562. }
  563. }
  564. return(list(ouch.tre=ouch.tre, ouch.dat=ouch.dat, ape.tre=td$phy, orig.dat=td$data[,1]))
  565. }
  566. summarize.run<-function(cOK, endparms, cur.model) {
  567. # end = c(nRateProp, nRateSwapProp, nRcatProp, nRootProp, nAlphaProp, nThetaProp, nThetaSwapProp, nTcatProp)
  568. names(endparms)<-names(cOK)<-c("rates", "rate.swap", "rate.catg", "root", "alpha", "theta", "theta.swap", "theta.catg")
  569. names(endparms)=paste("pr.",names(endparms),sep="")
  570. names(cOK)=paste("ok.",names(cOK),sep="")
  571. cat("\n")
  572. if(cur.model=="BM") {
  573. print(endparms[paste("pr.", c("rates", "rate.swap", "rate.catg", "root"),sep="")])
  574. print(cOK[paste("ok.", c("rates", "rate.swap", "rate.catg", "root"),sep="")])
  575. } else {
  576. print(endparms[paste("pr.", c("rates", "rate.swap", "rate.catg", "alpha", "theta", "theta.swap", "theta.catg"),sep="")])
  577. print(cOK[paste("ok.", c("rates", "rate.swap", "rate.catg", "alpha", "theta", "theta.swap", "theta.catg"),sep="")])
  578. }
  579. }
  580. generate.log<-function(bundled.parms, cur.model, file, init=FALSE) {
  581. # bundled.parms=list(gen=i, mrate=mean(cur.rates), cats=sum(cur.delta.rates)+1, root=cur.root, alpha=cur.alpha, reg=sum(cur.delta.theta)+1, theta=mean(cur.theta), lnL=curr.lnL)
  582. res=bundled.parms
  583. if(init) {
  584. if(cur.model=="BM") msg=paste("gen", "model", "mean.rate", "rates", "root", "lnL", sep="\t") else msg=paste("gen", "model", "mean.rate", "rates", "alpha", "regimes", "mean.theta", "lnL", sep="\t")
  585. } else {
  586. if(cur.model=="BM") {
  587. msg<-paste(res$gen, sQuote(cur.model), sprintf("%.3f", res$mrate), res$cats, sprintf("%.3f", res$root), sprintf("%.3f", res$lnL),sep="\t")
  588. } else {
  589. msg<-paste(res$gen, sQuote(cur.model), sprintf("%.3f", res$mrate), res$cats, sprintf("%.4f", res$alpha), res$reg, sprintf("%.3f", res$theta), sprintf("%.3f", res$lnL),sep="\t")
  590. }
  591. }
  592. write(msg, file=file, append=TRUE)
  593. }
  594. generate.error.message<-function(i, mod.cur, mod.new, lnR, errorLog, init=FALSE) {
  595. if(init) {
  596. initmsg = write(paste("gen", "curr.lnL", "new.lnL", "lnR", sep="\t"), file=errorLog)
  597. } else {
  598. write(paste(i, sprintf("%.3f", mod.cur$lnL), sprintf("%.3f", mod.new$lnL), sprintf("%.3f", lnR), sep="\t"), file=errorLog)
  599. }
  600. }
  601. determine.accepted.proposal<-function(startparms, endparms, cOK) {
  602. which(startparms!=endparms)->marker
  603. cOK[marker]=cOK[marker]+1
  604. cOK
  605. }
  606. rtpois<-function(N, lambda, k) {
  607. p=ppois(k, lambda, lower.tail=FALSE)
  608. out=qpois(runif(N, min=0, max=1-p), lambda)
  609. out
  610. }