/archive/source.archive/bm.and.ou/trait.rjmcmc.09132010a.R

http://github.com/eastman/auteur · R · 604 lines · 506 code · 58 blank · 40 comment · 119 complexity · 29f6a773a7d5b89b79c35b1ddb254c9f MD5 · raw file 

    

  1. ## Written by LJ HARMON & AH: 2008

    2. 

  2. ## Modified by JM EASTMAN: 09.2010

    3. 

  3. 

    4. 

  4. rjMCMC.trait<-function (phy, data, ngen=1000, sampleFreq=100, model=c("OU","BM"), probAlpha=0.05, probRegimes=0.10, probMergeSplit=0.05, probRoot=0.1, heat=1, fileBase="result") 

    5. 

  5. {

    6. 

  6. 	if(length(model)>1)stop(paste("Please specify either ", sQuote("OU"), " or ", sQuote("BM")))

    7. 

  7. 	

    8. 

  8. 	require(ouch)

    9. 

  9. 	require(geiger)

    10. 

  10. 	

    11. 

  11. ### prepare objects for rjMCMC

    12. 

  12. 	currModel		<- model								# initialize family of models explored

    13. 

  13. 

    14. 

  14. 	dataList		<- prepare.ouchdata(phy, data)			# check tree and data; convert into S4

    15. 

  15. 

    16. 

  16. 	ape.tre			<- currTree <- dataList$ape.tre			# S3 phylogeny

    17. 

  17. 	orig.dat		<- dataList$orig.dat					# S3 data	

    18. 

  18. 	ouch.tre		<- dataList$ouch.tre					# S4 OUCH phylogeny

    19. 

  19. 	ouch.dat		<- dataList$ouch.dat					# S4 OUCH data

    20. 

  20. 	

    21. 

  21. 	nn				<- Ntip(ape.tre)						# max number of rates (one per branch)

    22. 

  22. 	

    23. 

  23. 	currRates.c		<- rep(1,nn)							# initialize array for recording which branch belongs to which rate category

    24. 

  24.     currRates		<- rep(fitContinuous(ape.tre,orig.dat)[[1]]$beta,nn)							# initialize array for rate estimates branchwise

    25. 

  25. 	currRoot		<- mean(orig.dat)

    26. 

  26.     delta.rates		<- rep(0,nn)

    27. 

  27. 	currAlpha		<- 0.001								# initialize alpha

    28. 

  28. 	currTheta		<- 0									# initialize trait optima under OU

    29. 

  29. 	delta.theta		<- rep(0,nn)							# initialize local 'rates'

    30. 

  30. 	currRegimes		<- rep(1,nrow(ouch.dat))				# initial set of selective regimes

    31. 

  31. 	

    32. 

  32. 	nRateProp = 0											# proposal count: update rate parameter(s)

    33. 

  33. 	nRcatProp = 0											# proposal count: update rate categories

    34. 

  34. 	nRootProp = 0

    35. 

  35. 	nAlphaProp = 0											# proposal count: update alpha parameter

    36. 

  36. 	nRegimesProp = 0										# proposal count: update selective regimes

    37. 

  37. 	nThetaProp = 0

    38. 

  38. 	nRateOK = 0												# proposal count: successful rate changes

    39. 

  39. 	nRcatOK = 0												# proposal count: successful rate categories updates

    40. 

  40. 	nRootOK = 0

    41. 

  41. 	nAlphaOK = 0											# proposal count: successful alpha updates

    42. 

  42. 	nRegimesOK = 0											# proposal count: successful regime updates

    43. 

  43. 	nThetaOK = 0

    44. 

  44. 	l <- array(dim=c(ngen,4))								# array for storing output: size is number of MCMC generations

    45. 

  45.     l[] <- NA												# initialize array for storing lnLs

    46. 

  46. 	

    47. 

  47. 	cOK=c(													# list of proposal tallies

    48. 

  48. 		  nRateOK, 

    49. 

  49. 		  nRcatOK, 

    50. 

  50. 		  nRootOK, 

    51. 

  51. 		  nAlphaOK, 

    52. 

  52. 		  nRegimesOK, 

    53. 

  53. 		  nThetaOK

    54. 

  54. 		  )		

    55. 

  55. 	

    56. 

  56. 	errorLog=file(paste(getwd(),paste(currModel, "rjTraitErrors.log",sep="."),sep="/"),open='w+')

    57. 

  57. 	generate.error.message(i=NULL, modCurr=NULL, modNew=NULL, lnR=NULL, errorLog=errorLog, init=TRUE)

    58. 

  58. 	runLog=file(paste(getwd(),paste(currModel, "rjTrait.log",sep="."),sep="/"),open='w+')

    59. 

  59. 	generate.log(bundled.parms=NULL, currModel, file=runLog, init=TRUE)

    60. 

  60. 	

    61. 

  61. ### Begin rjMCMC

    62. 

  62.     for (i in 1:ngen) {

    63. 

  63. 		error.flag=FALSE

    64. 

  64.         lnLikelihoodRatio <- lnProposalRatio <- lnPriorRatio <- 0

    65. 

  65. 		startparms=list(rates=currRates, cats=currRates.c, root=currRoot, alph=currAlpha, reg=currRegimes, theta=currTheta)

    66. 

  66. 		

    67. 

  67. ## OU IMPLEMENTATION ## 

    68. 

  68. 		if(currModel=="OU") {

    69. 

  69. 			if (runif(1) < (2 * probMergeSplit)) {

    70. 

  70. 				if(runif(1)<0.5) { # adjust regimes

    71. 

  71. 					nr=split.or.merge(currTheta, currRegimes, decide.s.or.m(currRegimes, regimes=TRUE)->s.m, regimes=TRUE)

    72. 

  72. 					newRegimes=nr$new.categories

    73. 

  73. 					newTheta=nr$new.values

    74. 

  74. 					newRates=currRates

    75. 

  75. 					newRates.c=currRates.c

    76. 

  76. 					nRegimesProp=nRegimesProp+1

    77. 

  77. 				} else { # adjust rate categories

    78. 

  78. 					nr=split.or.merge(currRates, currRates.c, decide.s.or.m(currRates.c)->s.m)

    79. 

  79. 					newRates=nr$new.values

    80. 

  80. 					newRates.c=nr$new.categories

    81. 

  81. 					newRegimes=currRegimes

    82. 

  82. 					newTheta=currTheta

    83. 

  83. 					nRcatProp=nRcatProp+1

    84. 

  84. 				}

    85. 

  85. 				newAlpha=currAlpha

    86. 

  86. 				lnProposalRatio=nr$lnl.prop

    87. 

  87. 				lnPriorRatio=nr$lnl.prior

    88. 

  88. 			} else { # neither split nor merge

    89. 

  89. 				if(runif(1)<probAlpha) { # adjust alpha

    90. 

  90. 					newAlpha=adjust.rates(currAlpha)

    91. 

  91. 					newRates=currRates

    92. 

  92. 					newRates.c=currRates.c

    93. 

  93. 					newRegimes=currRegimes

    94. 

  94. 					newTheta=currTheta

    95. 

  95. 					nAlphaProp=nAlphaProp+1

    96. 

  96. 				} else { 

    97. 

  97. 					if(runif(1)<0.5) {# adjust rates

    98. 

  98. 						newRates=adjust.rates(currRates)

    99. 

  99. 						newRates.c=fix.categories(newRates)

    100. 

  100. 						newAlpha=currAlpha

    101. 

  101. 						newRegimes=currRegimes

    102. 

  102. 						newTheta=currTheta

    103. 

  103. 						nRateProp=nRateProp+1

    104. 

  104. 					} else { # adjust theta

    105. 

  105. 						newTheta=sapply(currTheta, adjust.value)

    106. 

  106. 						newAlpha=currAlpha

    107. 

  107. 						newRates=currRates

    108. 

  108. 						newRates.c=currRates.c

    109. 

  109. 						newRegimes=currRegimes

    110. 

  110. 						nThetaProp=nThetaProp+1

    111. 

  111. 					}

    112. 

  112. 				}

    113. 

  113. 			}

    114. 

  114. 			

    115. 

  115. 			currTree=apply.BMM.to.apetree(ape.tre, currRates)

    116. 

  116. 			newTree=apply.BMM.to.apetree(ape.tre, newRates)

    117. 

  117. 			modCurr=new.ou.lik.fn(currRates, currAlpha, currRegimes, currTheta, currTree, ouch.dat)

    118. 

  118. 			modNew=try(new.ou.lik.fn(newRates, newAlpha, newRegimes, newTheta, newTree, ouch.dat),silent=TRUE)

    119. 

  119. 		} else if(currModel=="BM"){

    120. 

  120. ## BM IMPLEMENTATION ##

    121. 

  121. 			if (runif(1) < (2 * probMergeSplit)) {

    122. 

  122. 				nr=split.or.merge(currRates, currRates.c, decide.s.or.m(currRates.c)->s.m)

    123. 

  123. 				newRates=nr$new.values

    124. 

  124. 				newRates.c=nr$new.categories

    125. 

  125. 				newRoot=currRoot

    126. 

  126. 				nRcatProp=nRcatProp+1

    127. 

  127. 				lnProposalRatio=nr$lnl.prop

    128. 

  128. 				lnPriorRatio=nr$lnl.prior

    129. 

  129. 			} else { 

    130. 

  130. 				if(runif(1)<probRoot) { # adjust root

    131. 

  131. 					newRoot=adjust.value(currRoot)

    132. 

  132. 					newRates=currRates

    133. 

  133. 					newRates.c=currRates.c

    134. 

  134. 					nRootProp=nRootProp+1

    135. 

  135. 				} else { # adjust rates

    136. 

  136. 					newRates=adjust.rates(currRates)

    137. 

  137. 					newRates.c=fix.categories(newRates)

    138. 

  138. 					newRoot=currRoot

    139. 

  139. 					nRateProp=nRateProp+1

    140. 

  140. 				}

    141. 

  141. 			}

    142. 

  142. 			

    143. 

  143. 			currTree=apply.BMM.to.apetree(ape.tre, currRates)

    144. 

  144. 			newTree=apply.BMM.to.apetree(ape.tre, newRates)

    145. 

  145. 			modCurr=new.bm.lik.fn(currRates, currRoot, currTree, ouch.dat)

    146. 

  146. 			modNew=try(new.bm.lik.fn(newRates, newRoot, newTree, ouch.dat),silent=TRUE)

    147. 

  147. 		}

    148. 

  148. 	

    149. 

  149. 		if(!inherits(modNew, "try-error")) {

    150. 

  150. 			lnLikelihoodRatio = modNew$lnL - modCurr$lnL

    151. 

  151. 		} else {

    152. 

  152. 			lnLikelihoodRatio = -Inf

    153. 

  153. 			modNew$lnL=NaN

    154. 

  154. 		}

    155. 

  155. 		r=assess.lnR((heat * lnLikelihoodRatio + heat * lnPriorRatio + lnProposalRatio)->lnR)

    156. 

  156. 		

    157. 

  157. 		if (runif(1) <= r$r) {			# adopt proposal

    158. 

  158. 			if(currModel=="OU") {

    159. 

  159. 				currAlpha <- newAlpha

    160. 

  160. 				currRegimes <- newRegimes

    161. 

  161. 				currTheta <- newTheta

    162. 

  162. 			} else {

    163. 

  163. 				currRoot <- newRoot

    164. 

  164. 			}

    165. 

  165. 			currRates <- newRates

    166. 

  166. 			currRates.c <- newRates.c

    167. 

  167. 			curr.lnL <- modNew$lnL

    168. 

  168. 			currTree <- newTree

    169. 

  169. 		} else {						# deny proposal			

    170. 

  170. 			curr.lnL <- modCurr$lnL 

    171. 

  171. 		}

    172. 

  172. 		if(r$error) generate.error.message(i, modCurr, modNew, lnR, errorLog)

    173. 

  173. 		

    174. 

  174. 		l[i,]<-c(curr.lnL, max(currRates.c), currAlpha, max(currRegimes))

    175. 

  175. 		endparms=list(rates=currRates, cats=currRates.c, root=currRoot, alph=currAlpha, reg=currRegimes, theta=currTheta)

    176. 

  176. 		cOK=tally.mcmc.parms(startparms, endparms, cOK)

    177. 

  177. 		

    178. 

  178. 		bundled.parms=list(gen=i, mrate=exp(mean(log(currRates))), cats=max(currRates.c), root=currRoot, alpha=currAlpha, reg=max(currRegimes), theta=mean(currTheta), lnL=curr.lnL)

    179. 

  179. 		if(i%%sampleFreq==0) generate.log(bundled.parms, currModel, file=runLog)

    180. 

  180. 	}

    181. 

  181. 	

    182. 

  182. # End rjMCMC

    183. 

  183. 	close(errorLog)

    184. 

  184. 	close(runLog)

    185. 

  185. 	summarize.run(cOK, (cProp<-c(nRateProp, nRcatProp, nRootProp, nAlphaProp, nRegimesProp, nThetaProp)), currModel)	

    186. 

  186. 	if(currModel=="BM") {

    187. 

  187. 		return(list(r.cat=l[,2], lnL=l[,1]))

    188. 

  188. 	} else {

    189. 

  189. 		return(list(r.cat=l[,2], r.alpha=l[,3], r.regimes=l[,4], lnL=l[,1]))

    190. 

  190. 	}

    191. 

  191. }

    192. 

  192. 

    193. 

  193. 

    194. 

  194. ## AUXILIARY FUNCTIONS

    195. 

  195. 

    196. 

  196. 

    197. 

  197. new.ou.lik.fn <- function(rates, alpha, regimes, theta, ape.tre, ouch.dat) { # from OUCH 2.7-1:::hansen()

    198. 

  198. # mod 09.13.2010

    199. 

  199. 	ouch.tre=ape2ouch(ape.tre, scale=FALSE)

    200. 

  200. 	ouch.dat$regimes=as.factor(regimes)

    201. 

  201. 	alpha.ouch=as.matrix(sqrt(alpha))

    202. 

  202. 	beta=ouch:::regime.spec(ouch.tre, ouch.dat['regimes'])

    203. 

  203. 	dat=do.call(c,lapply(ouch.dat['trait'],function(y)y[ouch.tre@term]))

    204. 

  204. 	n <- length(dat)

    205. 

  205. 	otree.df=as(ouch.tre,"data.frame")

    206. 

  206. 	ordering=otree.df[n:nrow(otree.df),'labels']

    207. 

  207. 	ev <- eigen(alpha.ouch,symmetric=TRUE)

    208. 

  208. 	w <- .Call(ouch:::ouch_weights,object=ouch.tre,lambda=ev$values,S=ev$vectors,beta=beta)

    209. 

  209. 

    210. 

  210. 	v <- geiger:::ouMatrix(vcv.phylo(ape.tre), alpha)

    211. 

  211. 	v.match <- match(row.names(v),ordering)

    212. 

  212. 	v <- v[v.match,v.match]

    213. 

  213. 		## GEIGER call does not assume ultrametricity (and thus allows different 'rates' under BM)

    214. 

  214. 		# -old-	v <- .Call(ouch:::ouch_covar,object=ouch.tre,lambda=ev$values,S=ev$vectors,sigma.sq=1)

    215. 

  215. 	

    216. 

  216. 	y <- matrix(theta,ncol=1)

    217. 

  217. 	e <- w%*%y-dat

    218. 

  218. 	dim(y) <- dim(y)[1]

    219. 

  219. 	dim(e) <- n

    220. 

  220. 	q <- e%*%solve(v,e)

    221. 

  221. 	det.v <- determinant(v,logarithm=TRUE)

    222. 

  222. 	if (det.v$sign!=1)stop("ou.lik.fn error: non-positive determinant",call.=FALSE)

    223. 

  223. 	dev <- n*log(2*pi)+as.numeric(det.v$modulus)+q[1,1]

    224. 

  224. 	list(

    225. 

  225. 		 theta=theta,

    226. 

  226. 		 weight=w,

    227. 

  227. 		 vcov=v,

    228. 

  228. 		 resids=e,

    229. 

  229. 		 lnL=-0.5*dev

    230. 

  230. 		 )

    231. 

  231. }

    232. 

  232. 

    233. 

  233. new.bm.lik.fn <- function(rates, root, ape.tre, ouch.dat) { # from OUCH 2.7-1:::brown()	

    234. 

  234. # mod 09.13.2010	

    235. 

  235. 	tree=ape2ouch(ape.tre, scale=FALSE)

    236. 

  236. 	data=ouch.dat['trait']

    237. 

  237. 	

    238. 

  238. 	if (!is(tree, "ouchtree")) {stop(sQuote("tree"), " must be an object of class ", sQuote("ouchtree"))}

    239. 

  239.     if (is.data.frame(data)) {

    240. 

  240.         nm <- rownames(data)

    241. 

  241.         data <- lapply(as.list(data), function(x) {names(x) <- nm; x})

    242. 

  242.     }

    243. 

  243.     if (is.numeric(data)) {

    244. 

  244.         nm <- deparse(substitute(data))[1]

    245. 

  245.         data <- list(data)

    246. 

  246.         names(data) <- nm

    247. 

  247.     }

    248. 

  248.     if (is.list(data)) {

    249. 

  249.         if (any(sapply(data, class) != "numeric") || any(sapply(data, length) != tree@nnodes)) {stop(sQuote("data"), " vector(s) must be numeric, with one entry per node of the tree")}

    250. 

  250.         if (any(sapply(data, function(x) (is.null(names(x))) || (!setequal(names(x), tree@nodes))))) {stop(sQuote("data"), " vector names (or data-frame row names) must match node names of ", sQuote("tree"))}

    251. 

  251.         for (xx in data) {

    252. 

  252.             no.dats <- which(is.na(xx[tree@nodes[tree@term]]))

    253. 

  253.             if (length(no.dats) > 0) {stop("missing data on terminal node(s): ", paste(tree@nodes[tree@term[no.dats]], collapse = ","))}

    254. 

  254.         }

    255. 

  255.     }

    256. 

  256.     nterm <- tree@nterm

    257. 

  257.     nchar <- length(data)

    258. 

  258.     if (is.null(names(data))) names(data) <- paste("char", 1:nchar, sep = "")

    259. 

  259.     dat.tmp <- lapply(data, function(x) x[tree@nodes])

    260. 

  260.     dat <- lapply(dat.tmp, function(y) y[tree@term])

    261. 

  261. 	dat <- dat[[1]]

    262. 

  262. 	nterm <- Ntip(ape.tre)

    263. 

  263. 	nchar <- 1

    264. 

  264. 	w <- matrix(data=1,nrow=nterm,ncol=1)

    265. 

  265. 	b <- tree@branch.times

    266. 

  266. 	y <- matrix(root)

    267. 

  267. 	e <- w%*%y-dat

    268. 

  268. 	q <- t(e)%*%solve(b,e)

    269. 

  269. 	v <- q/nterm

    270. 

  270. 	dev <- nchar*nterm*(1+log(2*pi))+nchar*log(det(b))+nterm*log(det(v))

    271. 

  271. 	list(

    272. 

  272. 		 root=root,

    273. 

  273. 		 lnL = -0.5*dev

    274. 

  274. 		 )	

    275. 

  275. }

    276. 

  276. 

    277. 

  277. rtpois<-function(N, lambda, k){

    278. 

  278. 	p=ppois(k, lambda, lower.tail=FALSE)

    279. 

  279. 	out=qpois(runif(N, min=0, max=1-p), lambda)

    280. 

  280. 	out

    281. 

  281. }

    282. 

  282. 

    283. 

  283. get.descendants.of.node <- function(node, phy) {

    284. 

  284. 	descendants.of.node <- function(node, phy) {

    285. 

  285. 		current.node=node[node%in%phy$edge[,1]]

    286. 

  286. 		if (length(current.node)>0) {

    287. 

  287. 			desc=c(sapply(current.node, function(x) {return(phy$edge[which(phy$edge[,1]==x),2])}))

    288. 

  288. 			write(desc,file="TEMP",append=TRUE)

    289. 

  289. 			descendants.of.node(desc,phy)

    290. 

  290. 		}	

    291. 

  291. 	}

    292. 

  292. 	descendants.of.node(node, phy)

    293. 

  293. 	if(file.exists("TEMP")) {

    294. 

  294. 		out=unlist(strsplit(readLines("TEMP")," "))

    295. 

  295. 		unlink("TEMP")

    296. 

  296. 		return(sort(unique(as.numeric(out))))

    297. 

  297. 	}

    298. 

  298. }	

    299. 

  299. 

    300. 

  300. get.ancestor.of.node<-function(node, phy) {

    301. 

  301. 	anc=phy$edge[which(phy$edge[,2]==node),1]

    302. 

  302. 	anc

    303. 

  303. }

    304. 

  304. 

    305. 

  305. assess.lnR <- function(lnR) {

    306. 

  306. 	if(is.na(lnR) || abs(lnR)==Inf) {

    307. 

  307. 		error=TRUE

    308. 

  308. 		r=0

    309. 

  309. 	} else {

    310. 

  310. 		if(lnR < -100) {

    311. 

  311. 			r=0 

    312. 

  312. 		} else if(lnR > 0) {

    313. 

  313. 			r=1 

    314. 

  314. 		} else {

    315. 

  315. 			r=exp(lnR)

    316. 

  316. 		}

    317. 

  317. 		error=FALSE

    318. 

  318. 	}

    319. 

  319. 	return(list(r=r, error=error))

    320. 

  320. }

    321. 

  321. 

    322. 

  322. decide.s.or.m <- function(categories, regimes=FALSE) {

    323. 

  323. 	if (max(categories) == 1) {

    324. 

  324. 		return("split")

    325. 

  325. 	}

    326. 

  326. 	else if (max(categories) == length(categories)) {

    327. 

  327. 		return("merge")

    328. 

  328. 	}

    329. 

  329. #	else if (regimes && (((length(categories)+1)/2)==max(categories))) {

    330. 

  330. # limits number of selective regimes to number of tips

    331. 

  331. #		return("merge")

    332. 

  332. #	}

    333. 

  333. 	else if (runif(1) < 0.5) {

    334. 

  334. 		return("split")

    335. 

  335. 	}

    336. 

  336. 	else return("merge")

    337. 

  337. }

    338. 

  338. 	

    339. 

  339. split.or.merge <- function(values, categories, task=c("split","merge"), regimes=FALSE) 

    340. 

  340. {	

    341. 

  341. # updated to adjust for an additional theta when regimes are split (or one fewer theta if regimes are merged)

    342. 

  342. # implement Drummon and Suchard 2010 method, adapting local relaxed clocks for local theta and BMMrates - draw changes in BMMrates from truncated Poisson

    343. 

  343. 	cc=categories

    344. 

  344. 

    345. 

  345. 	if(regimes) vv=values[cc] else vv=values

    346. 

  346. 	new.vv=vv

    347. 

  347. 	new.cc=cc

    348. 

  348. 	nn=length(vv)

    349. 

  349. 	

    350. 

  350. 	if(task=="merge") {

    351. 

  351. 		ncat <- max(cc)

    352. 

  352. 		m.cc <- sample(1:ncat)[1:2] # grab two categories

    353. 

  353. 		new.cc[new.cc == m.cc[1]] <- m.cc[2]	# merge selected categories

    354. 

  354. 		new.cc <- fix.categories(new.cc)

    355. 

  355. 		or1 <- vv[cc == m.cc[1]][1] # value affected

    356. 

  356. 		or2 <- vv[cc == m.cc[2]][1]	# value affected

    357. 

  357. 		n1 <- sum(cc == m.cc[1]) # no. entries affected

    358. 

  358. 		n2 <- sum(cc == m.cc[2]) # no. entries affected

    359. 

  359. 		nr <- (n1 * or1 + n2 * or2)/(n1 + n2) # find new value for merged category

    360. 

  360. 		new.vv[cc == m.cc[1] | cc == m.cc[2]] <- nr # new value assigned to merged classes

    361. 

  361. 		numSplittableBeforeMerge <- sum(table(cc) > 1)

    362. 

  362. 		numSplittableAfterMerge <- sum(table(new.cc) > 1)

    363. 

  363. 		if (max(cc) == nn) {

    364. 

  364. 			probMerge = 1

    365. 

  365. 		} else {

    366. 

  366. 			probMerge = 0.5

    367. 

  367. 		}

    368. 

  368. 		

    369. 

  369. 		if (max(new.cc) == 1) {

    370. 

  370. 			probSplit = 1

    371. 

  371. 		} else {

    372. 

  372. 			probSplit = 0.5

    373. 

  373. 		}

    374. 

  374. 		factor = probSplit/probMerge		

    375. 

  375. # AS WRITTEN by LJ HARMON

    376. 

  376. 		lnProposalRatio = log(factor) + log(nn) + log(ncat) - log(numSplittableAfterMerge) - log(2^(n1 + n2) - 2) - log(nr/mean(vv)) - log(n1 + n2)

    377. 

  377. 		lnPriorRatio = log(Stirling2(nn, ncat)) - log(Stirling2(nn, ncat - 1))

    378. 

  378. #	lnHastings = log((2N-2-K+1)/K) # where N is tips, K is number of local parms in tree

    379. 

  379. 	}

    380. 

  380. 	if(task=="split") {

    381. 

  381. 		ncat <- max(cc)

    382. 

  382. 		while (1) { # find non-unique category

    383. 

  383. 			rcat <- round(runif(1) * ncat + 0.5) # class to split

    384. 

  384. 			nc <- sum(cc == rcat) 

    385. 

  385. 			if (nc > 1) 

    386. 

  386. 			break

    387. 

  387. 		}

    388. 

  388. 		while (1) { # 

    389. 

  389. 			new <- round(runif(nc) * 2 + 0.5)

    390. 

  390. 			if (length(table(new)) == 2) 

    391. 

  391. 			break

    392. 

  392. 		}

    393. 

  393. 		new.cc[cc == rcat][new == 2] <- ncat + 1

    394. 

  394. 		or <- vv[cc == rcat][1] # old rate

    395. 

  395. 		n1 <- sum(new == 1)	# no. elements affected

    396. 

  396. 		n2 <- sum(new == 2)	# no. elements affected

    397. 

  397. 		lims <- c(-0.5 * n1 * or, 0.5 * n2 * or)

    398. 

  398. 		lims <- lims[order(lims)]

    399. 

  399. 		u <- runif(1, min = lims[1], max = lims[2])

    400. 

  400. 		nr1 <- or + u/n1 # new rate 1

    401. 

  401. 		nr2 <- or - u/n2 # new rate 2

    402. 

  402. 		new.vv[new.cc == rcat] <- nr1 # assign new rate 1

    403. 

  403. 		new.vv[new.cc == ncat + 1] <- nr2 # assign new rate 2

    404. 

  404. 		new.cc <- fix.categories(new.cc)

    405. 

  405. 		numSplittableBeforeSplit <- sum(table(cc) > 1)

    406. 

  406. 		numSplittableAfterSplit <- sum(table(new.cc) > 1)

    407. 

  407. 		if (max(cc) == 1) {

    408. 

  408. 			probSplit = 1

    409. 

  409. 		} else {

    410. 

  410. 			probSplit = 0.5

    411. 

  411. 		}

    412. 

  412. 		

    413. 

  413. 		if (max(new.cc) == nn) {

    414. 

  414. 			probMerge = 1

    415. 

  415. 		} else {

    416. 

  416. 			probMerge = 0.5

    417. 

  417. 		}

    418. 

  418. 		

    419. 

  419. 		factor = probMerge/probSplit

    420. 

  420. 		

    421. 

  421. # AS WRITTEN by LJ HARMON

    422. 

  422. 		lnProposalRatio = log(factor) + log(numSplittableBeforeSplit) + log(2^(n1 + n2) - 2) + log(or/mean(vv)) + log(n1 + n2) - log(nn) - log(ncat + 1)

    423. 

  423. 		lnPriorRatio = log(Stirling2(nn, ncat)) - log(Stirling2(nn, ncat + 1))

    424. 

  424. #		lnHastings=log((K+1)/(2N-2-K)) # where N is tips, K is number of local parms in tree

    425. 

  425. 

    426. 

  426. 	}

    427. 

  427. 	if(regimes) new.vv=unique(new.vv[order(new.cc)])

    428. 

  428. 	list(new.values=new.vv, new.categories=new.cc, lnl.prop=lnProposalRatio, lnl.prior=lnPriorRatio)

    429. 

  429. }

    430. 

  430. 

    431. 

  431. 

    432. 

  432. fix.categories <- function(x)

    433. 

  433. {

    434. 

  434. 	f<-factor(x)

    435. 

  435. 	n<-nlevels(f)

    436. 

  436. 	o<-numeric(n)

    437. 

  437. 	for(i in 1:n) o[i]<-which(f==levels(f)[i])[1]

    438. 

  438. 	nf<-ordered(f, levels(f)[order(o)])	

    439. 

  439. 	nx<-as.numeric(nf)

    440. 

  440. 	nx

    441. 

  441. }

    442. 

  442. 

    443. 

  443. choose.one <- function(bit.vector)

    444. 

  444. {

    445. 

  445. 	bb=bit.vector

    446. 

  446. 	s=sample(1:length(bb),1)

    447. 

  447. 	bb[s]=1-bb[s]

    448. 

  448. 	bb

    449. 

  449. }

    450. 

  450. 

    451. 

  451. fix.heritable.values<-function(bit.vector, values, phy) {

    452. 

  452. # assuming most recently changed value in correspondence with bit.vector has been changed on entry

    453. 

  453. 	bb=bit.vector

    454. 

  454. 	vv=values

    455. 

  455. 	names(vv)<-names(bb)<-phy$edge[,2]

    456. 

  456. 	new.bb=choose.one(bb)

    457. 

  457. #	new.vv=vv (do adjust value on the element that was just changed (if 0->1)

    458. 

  458. 	shifts=sort(as.numeric(names(which(bb>0))))	

    459. 

  459. 	desc=lapply(shifts, function(x)get.descendants.of.node(x,phy))

    460. 

  460. # fix if descendants of one set are included in another

    461. 

  461. 	new.vv=match()

    462. 

  462. }

    463. 

  463. 

    464. 

  464. assign.edgelengths<-function(ouchtree) {

    465. 

  465. 	ouchtree=as(ouchtree, "data.frame")

    466. 

  466. 	ouchtree$edges=NA

    467. 

  467. 	for(i in 2:nrow(ouchtree)) {

    468. 

  468. 		ouchtree$edges[i]=ouchtree$times[i]-(ouchtree[match(ouchtree[i,'ancestors'], ouchtree[,'nodes']),'times'])

    469. 

  469. 	}

    470. 

  470. 	ouchtree$edges[1]=0

    471. 

  471. 	ouchtree

    472. 

  472. }

    473. 

  473. 

    474. 

  474. apply.BMM.to.ouchtree<-function(ouchtree, rates) {

    475. 

  475. ## FIXME: ensure that the tree maintains structure before and after this function is called

    476. 

  476. ## VERY SLOW 

    477. 

  477. ## NOT CURRENTLY implemented

    478. 

  478. 	ot=as(ouchtree, "data.frame")

    479. 

  479. 	ot=assign.edgelengths(ot)

    480. 

  480. 	ot$edges.new=c(0, rates)*ot$edges

    481. 

  481. 	ot=ot[order(ot$times),]

    482. 

  482. 	ot$times.new=NA

    483. 

  483. 	ot$times.new[1]=0

    484. 

  484. 	for(i in 2:nrow(ot)) {

    485. 

  485. 		ot$times.new[i]=ot[match(ot[i,'ancestors'], ot[,'nodes']),'times.new']+ot[i,'edges.new']

    486. 

  486. 	}

    487. 

  487. 	ot$times=ot$times.new

    488. 

  488. 	return(with(ot, ouchtree(nodes=nodes, ancestors=ancestors, times=times, labels=labels)))

    489. 

  489. }

    490. 

  490. 

    491. 

  491. apply.BMM.to.apetree<-function(apetree, rates) {

    492. 

  492. ## FASTER than apply.BMM.to.ouchtree

    493. 

  493. 	apetree$edge.length=apetree$edge.length*rates

    494. 

  494. 	return(apetree)

    495. 

  495. }

    496. 

  496. 

    497. 

  497. adjust.value <- function(value) {

    498. 

  498. 	rr=value

    499. 

  499. 	rch <- runif(1, min = -0.5, max = 0.5)

    500. 

  500. 	rr <- rr + rch

    501. 

  501. 	rr

    502. 

  502. }

    503. 

  503. 

    504. 

  504. adjust.rates <- function(values) {

    505. 

  505. 	vv=values

    506. 

  506. 	ncat <- max(fix.categories(vv)->cc)

    507. 

  507. 	rcat <- round(runif(1) * (ncat) + 0.5)

    508. 

  508. 	r <- log(vv[cc == rcat][1])

    509. 

  509. 	rch <- runif(1, min = -2, max = 2)

    510. 

  510. 	nr <- r + rch

    511. 

  511. 	vv[cc == rcat] <- exp(nr)

    512. 

  512. 	vv

    513. 

  513. }

    514. 

  514. 

    515. 

  515. prepare.ouchdata <- function(phy, data) {

    516. 

  516. 	td <- treedata(phy, data, sort = T)	

    517. 

  517. 	ouch.tre <- ape2ouch(td$phy->ape.tre, scale=FALSE)		

    518. 

  518. 	ouch.dat=data.frame(as.numeric(ouch.tre@nodes), as.character(ouch.tre@nodelabels))

    519. 

  519. 	names(ouch.dat)=c("nodes","labels")

    520. 

  520. 	ouch.dat[ouch.dat[,2]=="",2]=NA

    521. 

  521. 	ouch.dat$trait=NA

    522. 

  522. 	mM=match(ouch.dat$labels,names(data))

    523. 

  523. 	for(i in 1:nrow(ouch.dat)){

    524. 

  524. 		if(!is.na(mM[i])){

    525. 

  525. 			ouch.dat$trait[i]=data[mM[i]]	

    526. 

  526. 		}

    527. 

  527. 	}

    528. 

  528. 	return(list(ouch.tre=ouch.tre, ouch.dat=ouch.dat, ape.tre=td$phy, orig.dat=td$data[,1]))

    529. 

  529. }

    530. 

  530. 

    531. 

  531. summarize.run<-function(cOK, cProp, currModel) {

    532. 

  532. 	names(cProp)<-names(cOK)<-c("rates", "rcats", "root", "alpha", "sRegimes", "theta")

    533. 

  533. 	names(cProp)=paste("prop.",names(cProp),sep="")

    534. 

  534. 	names(cOK)=paste("okay.",names(cOK),sep="")

    535. 

  535. 	cat("\n")

    536. 

  536. 	if(currModel=="BM") {

    537. 

  537. 		print(cProp[1:3])

    538. 

  538. 		print(cOK[1:3])

    539. 

  539. 	} else {

    540. 

  540. 		print(cProp[c(1,2,4,5,6)])

    541. 

  541. 		print(cOK[c(1,2,4,5,6)])

    542. 

  542. 	}

    543. 

  543. }

    544. 

  544. 

    545. 

  545. generate.log<-function(bundled.parms, currModel, file, init=FALSE) {

    546. 

  546. #	bundled.parms=list(gen=i, mrate=mean(currRates), cats=max(currRates.c), root=currRoot, alpha=currAlpha, reg=max(currRegimes), theta=mean(currTheta), lnL=curr.lnL)

    547. 

  547. 	res=bundled.parms

    548. 

  548. 

    549. 

  549. 	if(init) {

    550. 

  550. 		if(currModel=="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")

    551. 

  551. 	} else {

    552. 

  552. 		if(currModel=="BM") {

    553. 

  553. 			msg<-paste(res$gen, sQuote(currModel), sprintf("%.3f", res$mrate), res$cats, sprintf("%.3f", res$root), sprintf("%.3f", res$lnL),sep="\t")

    554. 

  554. 		} else {

    555. 

  555. 			msg<-paste(res$gen, sQuote(currModel), sprintf("%.3f", res$mrate), res$cats, sprintf("%.4f", res$alpha), res$reg, sprintf("%.3f", res$theta), sprintf("%.3f", res$lnL),sep="\t")

    556. 

  556. 		}

    557. 

  557. 	}

    558. 

  558. 	write(msg, file=file, append=TRUE)

    559. 

  559. }

    560. 

  560. 

    561. 

  561. generate.error.message<-function(i, modCurr, modNew, lnR, errorLog, init=FALSE) {

    562. 

  562. 	if(init) {

    563. 

  563. 		initmsg = write(paste("gen", "curr.lnL", "new.lnL", "lnR", sep="\t"), file=errorLog)

    564. 

  564. 	} else {

    565. 

  565. 		write(paste(i, sprintf("%.3f", modCurr$lnL), sprintf("%.3f", modNew$lnL), sprintf("%.3f", lnR), sep="\t"), file=errorLog)

    566. 

  566. 	}

    567. 

  567. }

    568. 

  568. 

    569. 

  569. tally.mcmc.parms<-function(startparms, endparms, cOK) {

    570. 

  570. 	index=array(dim=(length(startparms)-1))

    571. 

  571. 	for(i in 1:(length(startparms)-1)) {

    572. 

  572. 		index[i]=all(startparms[[i]]==endparms[[i]])

    573. 

  573. 	}

    574. 

  574. 	if(!all(index)) {

    575. 

  575. 		cOK[max(which(index==FALSE))]=cOK[max(which(index==FALSE))]+1

    576. 

  576. 	}

    577. 

  577. 	if(length(startparms$theta) == length(endparms$theta)) {

    578. 

  578. 		if(!all(sort(startparms$theta)==sort(endparms$theta))) {

    579. 

  579. 			cOK[6]=cOK[6]+1

    580. 

  580. 		}

    581. 

  581. 	}

    582. 

  582. 	cOK

    583. 

  583. }

    584. 

  584. 	

    585. 

  585. Stirling2 <- function(n,m)

    586. 

  586. {

    587. 

  587.     ## Purpose:  Stirling Numbers of the 2-nd kind

    588. 

  588.     ## 		S^{(m)}_n = number of ways of partitioning a set of

    589. 

  589.     ##                      $n$ elements into $m$ non-empty subsets

    590. 

  590.     ## Author: Martin Maechler, Date:  May 28 1992, 23:42

    591. 

  591.     ## ----------------------------------------------------------------

    592. 

  592.     ## Abramowitz/Stegun: 24,1,4 (p. 824-5 ; Table 24.4, p.835)

    593. 

  593.     ## Closed Form : p.824 "C."

    594. 

  594.     ## ----------------------------------------------------------------

    595. 

  595. 

    596. 

  596.     if (0 > m || m > n) stop("'m' must be in 0..n !")

    597. 

  597.     k <- 0:m

    598. 

  598.     sig <- rep(c(1,-1)*(-1)^m, length= m+1)# 1 for m=0; -1 1 (m=1)

    599. 

  599.     ## The following gives rounding errors for (25,5) :

    600. 

  600.     ## r <- sum( sig * k^n /(gamma(k+1)*gamma(m+1-k)) )

    601. 

  601.     ga <- gamma(k+1)

    602. 

  602.     round(sum( sig * k^n /(ga * rev(ga))))

    603. 

  603. }
    604.