/working/auteur.extended/R/rjmcmc.bm.R

  1#function for Markov sampling from a range of model complexities under the general process of Brownian motion evolution of continuous traits
  2#author: LJ HARMON 2009, A HIPP 2009, and JM EASTMAN 2010
  3
  4rjmcmc.bm <-
  5function (phy, dat, SE=0, ngen=1000, sample.freq=100, prob.mergesplit=0.05, prob.root=0.01, lambdaK=log(2), constrainK=FALSE, jumpsize=NULL, simplestart=FALSE, internal.only=FALSE, summary=TRUE, fileBase="result") 
  6{ 
  7	model="BM"
  8	heat=1 ## heating not currently implemented
  9	require(geiger)
 10	
 11	if(is.null(jumpsize)) {
 12		cat("CALIBRATING jumpsize...\n")
 13		adjustable.jumpsize=TRUE
 14		jumpsize=calibrate.jumpsize(phy, dat, nsteps=ngen/1000, model)
 15	} else {
 16		adjustable.jumpsize=FALSE
 17	}
 18	
 19
 20### prepare data for rjMCMC
 21	cur.model		<- model
 22	dataList		<- prepare.data(phy, dat, SE)			
 23	ape.tre			<- dataList$ape.tre			
 24	orig.dat		<- dataList$orig.dat
 25	SE				<- dataList$SE
 26	node.des		<- sapply(unique(c(ape.tre$edge[1,1],ape.tre$edge[,2])), function(x) get.descendants.of.node(x, ape.tre))
 27	names(node.des) <- c(ape.tre$edge[1,1], unique(ape.tre$edge[,2]))
 28
 29# initialize parameters
 30	if(is.numeric(constrainK) & (constrainK > length(ape.tre$edge) | constrainK < 1)) stop("Constraint on rate shifts is nonsensical. Ensure that constrainK is at least 1 and less than the number of available nodes in the tree.")
 31
 32	if(simplestart | is.numeric(constrainK) | internal.only) {
 33		if(is.numeric(constrainK)) {
 34			init.rate	<- generate.starting.point(orig.dat, ape.tre, node.des, model, K=constrainK, jumpsize=jumpsize)
 35		} else {
 36			init.rate	<- list(values=rep(fit.continuous(ape.tre,orig.dat),length(ape.tre$edge.length)),delta=rep(0,length(ape.tre$edge.length)))
 37		}
 38	} else {
 39		init.rate	<- generate.starting.point(orig.dat, ape.tre, node.des, model, K=constrainK, jumpsize=jumpsize )
 40	}
 41		
 42    cur.rates		<- init.rate$values			
 43	cur.delta.rates	<- init.rate$delta
 44	cur.root		<- adjustvalue(mean(orig.dat), jumpsize)
 45	cur.vcv			<- updatevcv(ape.tre, cur.rates)
 46	
 47	mod.cur = bm.lik.fn(cur.rates, cur.root, orig.dat, cur.vcv, SE)
 48		
 49# proposal counts
 50	nRateProp =		0	
 51	nRateSwapProp = 0									
 52	nRcatProp =		0										
 53	nRootProp =		0
 54	nRateOK =		0											
 55	nRateSwapOK =	0
 56	nRcatOK =		0											
 57	nRootOK =		0
 58
 59	cOK=c(												
 60		nRateOK,
 61		nRateSwapOK,
 62		nRcatOK, 
 63		nRootOK 
 64	)
 65	prob.rates=1-sum(c(2*prob.mergesplit, prob.root))
 66	if(prob.rates<0) stop("proposal frequencies must not exceed 1; adjust 'prob.mergesplit' and (or) 'prob.root'")
 67	proposal.rates<-orig.proposal.rates<-c(2*prob.mergesplit, prob.root, prob.rates)
 68	prop.cs<-orig.prop.cs<-cumsum(proposal.rates)
 69		
 70# find jumpsize calibration generations
 71	if(adjustable.jumpsize) {
 72		tf=c()
 73		tenths=function(v)floor(0.1*v)
 74		nn=ngen
 75		while(1) {
 76			vv=tenths(nn)
 77			nn=vv
 78			if(nn<=1) break	else tf=c(nn,tf)		
 79		}
 80		tuneFreq=tf
 81	} else {
 82		tuneFreq=0
 83	}
 84	
 85	cur.acceptancerate=0
 86	
 87	tickerFreq=ceiling((ngen+max(tuneFreq))/30)
 88	
 89# file handling
 90	if(summary) parmBase=paste(model, fileBase, "parameters/",sep=".") else parmBase=paste(".", runif(1), model, fileBase, "parameters/",sep="")
 91	if(!file.exists(parmBase)) dir.create(parmBase)
 92	parlogger(model=model, init=TRUE, node.des, parmBase=parmBase)
 93	errorLog=paste(parmBase,paste(cur.model, fileBase, "rjmcmc.errors.log",sep="."),sep="/")
 94	runLog=file(paste(parmBase,paste(cur.model, fileBase, "rjmcmc.log",sep="."),sep="/"),open='w+')
 95	generate.log(bundled.parms=NULL, cur.model, file=runLog, init=TRUE)
 96
 97### Begin rjMCMC
 98    for (i in 1:(ngen+max(tuneFreq))) {
 99		
100        lnLikelihoodRatio <- lnHastingsRatio <- lnPriorRatio <- 0
101		startparms = c(nRateProp, nRateSwapProp, nRcatProp, nRootProp)
102		
103		if(internal.only) {
104				tips.tmp=which(as.numeric(names(cur.delta.rates))<=Ntip(ape.tre))
105				if(any(cur.delta.rates[tips.tmp]==1)) stop("Broken internal only sampling")
106		}
107
108		
109## BM IMPLEMENTATION ##
110		while(1) {
111			cur.proposal=min(which(runif(1)<prop.cs))
112			if (cur.proposal==1 & !constrainK) {									# adjust rate categories
113				nr=splitormerge(cur.delta.rates, cur.rates, ape.tre, node.des, lambdaK, theta=FALSE, internal.only)
114				new.rates=nr$new.values
115				new.delta.rates=nr$new.delta
116				new.root=cur.root
117				nRcatProp=nRcatProp+1
118				lnHastingsRatio=nr$lnHastingsRatio
119				lnPriorRatio=nr$lnPriorRatio
120				break()
121			} else if(cur.proposal==2) {											# adjust root
122				new.root=adjustvalue(cur.root, jumpsize)
123				new.rates=cur.rates
124				new.delta.rates=cur.delta.rates
125				nRootProp=nRootProp+1
126				break()
127			} else if(cur.proposal==3){													
128				if(runif(1)>0.05 & sum(cur.delta.rates)>1) {						# tune local rate
129					new.rates=tune.rate(cur.rates, jumpsize)
130					new.delta.rates=cur.delta.rates 
131					new.root=cur.root
132					nRateSwapProp=nRateSwapProp+1
133					break()
134				} else {															# adjust rates
135					new.rates=adjustrate(cur.rates, jumpsize)
136					new.delta.rates=cur.delta.rates
137					new.root=cur.root
138					nRateProp=nRateProp+1
139					break()
140				} 
141			}
142		}
143		if(any(new.rates!=cur.rates)) new.vcv=updatevcv(ape.tre, new.rates) else new.vcv=cur.vcv
144		mod.new=NULL
145		mod.new=try(bm.lik.fn(new.rates, new.root, orig.dat, new.vcv, SE), silent=TRUE)
146		
147		if(inherits(mod.new, "try-error")) {mod.new=as.list(mod.new); mod.new$lnL=-Inf}
148		if(!is.infinite(mod.new$lnL)) {
149			lnLikelihoodRatio = mod.new$lnL - mod.cur$lnL
150		} else {
151			mod.new$lnL=-Inf
152			lnLikelihoodRatio = -Inf
153		}
154		
155	# compare likelihoods
156		endparms = c(nRateProp=nRateProp, nRateSwapProp=nRateSwapProp, nRcatProp=nRcatProp, nRootProp=nRootProp)
157
158		r=assess.lnR((heat * lnLikelihoodRatio + heat * lnPriorRatio + lnHastingsRatio)->lnR)
159
160		# potential errors
161		if(is.infinite(mod.cur$lnL)) stop("starting point has exceptionally poor likelihood")
162		if(r$error) generate.error.message(Ntip(phy), i, mod.cur, mod.new, lnLikelihoodRatio, lnPriorRatio, lnHastingsRatio, cur.delta.rates, new.delta.rates, errorLog)
163		
164		if (runif(1) <= r$r) {			## adopt proposal ##
165			decision="adopt"
166			cur.root <- new.root
167			cur.rates <- new.rates
168			cur.delta.rates <- new.delta.rates
169			mod.cur <- mod.new
170			cur.vcv <- new.vcv
171			curr.lnL <- mod.new$lnL
172			cOK <- determine.accepted.proposal(startparms, endparms, cOK)
173		} else {						## deny proposal ##	
174			decision="reject" 
175			curr.lnL <- mod.cur$lnL 
176		}
177		
178		# iteration-specific functions
179		if(i%%tickerFreq==0 & summary) {
180			if(i==tickerFreq) cat("|",rep(" ",9),toupper("generations complete"),rep(" ",9),"|","\n")
181			cat(". ")
182		}
183		
184		if(i%%sample.freq==0 & i>max(tuneFreq)) {
185			bundled.parms=list(gen=i-max(tuneFreq), mrate=exp(mean(log(cur.rates))), cats=sum(cur.delta.rates)+1, root=cur.root, lnL=curr.lnL, rates=cur.rates)
186			generate.log(bundled.parms, cur.model, file=runLog)			
187			parlogger(model=model, init=FALSE, node.des=node.des, i=i-max(tuneFreq), curr.lnL=curr.lnL, cur.root=cur.root, cur.rates=cur.rates, cur.delta.rates=cur.delta.rates, parmBase=parmBase)
188		}
189									   
190		if(any(tuneFreq%in%i) & adjustable.jumpsize) {
191			new.acceptancerate=sum(cOK)/sum(endparms)
192			if((new.acceptancerate<cur.acceptancerate) & adjustable.jumpsize) {
193				jumpsize=calibrate.jumpsize(phy, dat, nsteps=1000, model, jumpsizes=c(exp(log(jumpsize)-log(2)), exp(log(jumpsize)+log(2))))
194			}
195			cur.acceptancerate=new.acceptancerate
196		}
197	}
198	
199# End rjMCMC
200	close(runLog)
201	
202  # clear out jumpsize calibration directories if calibration run
203	if(summary) {
204		summarize.run(cOK, endparms, cur.model)	
205		cleanup.files(parmBase, cur.model, fileBase)
206	} else {
207		unlink(parmBase, recursive=TRUE)
208	}
209	return(list(acceptance.rate=sum(cOK)/sum(endparms), jumpsize=jumpsize))
210}
211