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

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