PageRenderTime 28ms CodeModel.GetById 16ms app.highlight 8ms RepoModel.GetById 1ms app.codeStats 0ms

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

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