PageRenderTime 44ms CodeModel.GetById 16ms app.highlight 21ms RepoModel.GetById 1ms app.codeStats 0ms

/in.progress/auteur/R/rjmcmc.R

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