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

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