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

http://github.com/eastman/auteur · R · 355 lines · 301 code · 32 blank · 22 comment · 74 complexity · 4beb9e3517f6a4403fb5ecfed0ad67ac MD5 · raw file 

    

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

    2. 

  2. ## Updated by JM EASTMAN: 09.2010

    3. 

  3. 

    4. 

  4. rjMCMC.trait<-function (phy, data, ngen = 1000, model = c("OU","BM"), probAlpha=0.05, probRegimes=0.10, probMergeSplit = 0.05, 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. 	dataList <- prepare.ouchdata(phy, data)					# check tree and data; convert into S4

    13. 

  13. 	ape.tre <- dataList$ape.tre								# S3 phylogeny

    14. 

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

    15. 

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

    16. 

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

    17. 

  17. 	nn <- length(ape.tre$edge.length)						# max number of rates (one per branch)

    18. 

  18. 	currRates.c <- numeric(nn)								# initialize array for recording which branch belongs to which rate category

    19. 

  19.     currRates.c[] <- 1										# assign one rate to all branches initially

    20. 

  20.     currRates <- numeric(nn)								# initialize array for rate estimates branchwise

    21. 

  21.     currRates[] <- init <- 0.001							# initialize BMMrate to be small

    22. 

  22. 	currAlpha <- 0.001										# initialize alpha

    23. 

  23. 	currRegimes <- numeric(nrow(ouch.dat))					# initial set of selective regimes

    24. 

  24. 	currRegimes[] <- 1										# initialize regime to be global

    25. 

  25. 	currModel <- model										# allow user to define family of models explored

    26. 

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

    27. 

  27. 	nRcatProp = 0

    28. 

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

    29. 

  29. 	nRegimesProp = 0

    30. 

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

    31. 

  31. 	nRcatOK = 0

    32. 

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

    33. 

  33. 	nRegimesOK = 0

    34. 

  34. 	cOK=c(nRateOK, nRcatOK, nAlphaOK, nRegimesOK)

    35. 

  35. 	l <- array(dim=ngen)									# number of MCMC generations

    36. 

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

    37. 

  37. 	

    38. 

  38. ### Begin rjMCMC

    39. 

  39. 	cat("\trates\tregimes\t\talpha\t\tlnL\t\t\tmodel")

    40. 

  40.     for (i in 1:ngen) {

    41. 

  41.         lnLikelihoodRatio <- lnProposalRatio <- lnPriorRatio <- 0

    42. 

  42. 		startparms=list(currRates, currRates.c, currAlpha, currRegimes)

    43. 

  43. 

    44. 

  44. ## OU IMPLEMENTATION ## thetas currently optimized by ouch:::glssoln(); selective regimes are choosable (see Beta)

    45. 

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

    46. 

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

    47. 

  47. 				adj.regimes <- adj.rates <- FALSE

    48. 

  48. 				if(runif(1) < 0.5) adj.regimes=TRUE else adj.rates=TRUE

    49. 

  49. 				if(adj.regimes) {

    50. 

  50. 					nr=split.or.merge(currRegimes, currRegimes, decide.s.or.m(currRegimes)->s.m)

    51. 

  51. 					newRegimes=nr$new.categories

    52. 

  52. 					newRates=currRates

    53. 

  53. 					newRates.c=currRates.c

    54. 

  54. 					nRegimesProp=nRegimesProp+1

    55. 

  55. 				} else { # adjust rate categories

    56. 

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

    57. 

  57. 					newRates=nr$new.values

    58. 

  58. 					newRates.c=nr$new.categories

    59. 

  59. 					newRegimes=currRegimes

    60. 

  60. 					nRcatProp=nRcatProp+1

    61. 

  61. 				}

    62. 

  62. 				newAlpha=currAlpha

    63. 

  63. 				lnProposalRatio=nr$lnl.prop

    64. 

  64. 				lnPriorRatio=nr$lnl.prior

    65. 

  65. 			} else { # neither split nor merge

    66. 

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

    67. 

  67. 					newAlpha=adjust.rates(currAlpha)

    68. 

  68. 					newRates=currRates

    69. 

  69. 					newRates.c=currRates.c

    70. 

  70. 					nAlphaProp=nAlphaProp+1

    71. 

  71. 				} else { # adjust rates

    72. 

  72. 					newRates=adjust.rates(currRates)

    73. 

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

    74. 

  74. 					newAlpha=currAlpha

    75. 

  75. 					nRateProp=nRateProp+1

    76. 

  76. 				}

    77. 

  77. 				newRegimes=currRegimes

    78. 

  78. 			}			

    79. 

  79. 			modCurr=lnl.OU(currRates, currAlpha, currRegimes, ape.tre, ouch.dat)

    80. 

  80. 			modNew=lnl.OU(newRates, newAlpha, newRegimes, ape.tre, ouch.dat)

    81. 

  81. 			

    82. 

  82. 			lnLikelihoodRatio = modNew$lnL - modCurr$lnL

    83. 

  83. 			r=assess.lnR(heat * lnLikelihoodRatio + heat * lnPriorRatio + lnProposalRatio)

    84. 

  84. 			print(lnPriorRatio+lnProposalRatio)

    85. 

  85. 	

    86. 

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

    87. 

  87. 				currRates <- newRates

    88. 

  88. 				currRates.c <- newRates.c

    89. 

  89. 				currAlpha <- newAlpha

    90. 

  90. 				currRegimes <- newRegimes

    91. 

  91. 				best.lnL <- modNew$lnL

    92. 

  92. 			} else {						# deny proposal			

    93. 

  93. 				best.lnL <- modCurr$lnL  

    94. 

  94. 			}

    95. 

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

    96. 

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

    97. 

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

    98. 

  98. 				newRates=nr$new.values

    99. 

  99. 				newRates.c=nr$new.categories

    100. 

  100. 				nRcatProp=nRcatProp+1

    101. 

  101. 				lnProposalRatio=nr$lnl.prop

    102. 

  102. 				lnPriorRatio=nr$lnl.prior

    103. 

  103. 			} else { # neither split nor merge

    104. 

  104. 				newRates=adjust.rates(currRates)

    105. 

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

    106. 

  106. 				nRateProp=nRateProp+1

    107. 

  107. 			}

    108. 

  108. 			

    109. 

  109. 			modCurr=lnl.BM(currRates, ape.tre, ouch.dat)

    110. 

  110. 			modNew=lnl.BM(newRates, ape.tre, ouch.dat)

    111. 

  111. 			

    112. 

  112. 			lnLikelihoodRatio = modNew$lnL - modCurr$lnL

    113. 

  113. 			r=assess.lnR(heat * lnLikelihoodRatio + heat * lnPriorRatio + lnProposalRatio)

    114. 

  114. 			print(lnPriorRatio+lnProposalRatio)

    115. 

  115. 			

    116. 

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

    117. 

  117. 				currRates <- newRates

    118. 

  118. 				currRates.c <- newRates.c

    119. 

  119. 				best.lnL <- modNew$lnL

    120. 

  120. 			} else {						# deny proposal			

    121. 

  121. 				best.lnL <- modCurr$lnL  

    122. 

  122. 			}

    123. 

  123. 		}

    124. 

  124. 		

    125. 

  125. 		l[i]=best.lnL

    126. 

  126. #		if(currModel=="BM") {

    127. 

  127. #			cat(paste("\n", max(currRates.c),"\t\t\t", round(best.lnL,2), currModel, sep="\t\t"))

    128. 

  128. #		} else {

    129. 

  129. #			cat(paste("\n", max(currRates.c), max(currRegimes), round(currAlpha,5), round(best.lnL,2), currModel, sep="\t\t"))

    130. 

  130. #		}

    131. 

  131. 		

    132. 

  132. 		endparms=list(currRates, currRates.c, currAlpha, currRegimes)

    133. 

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

    134. 

  134. 	}

    135. 

  135. # End rjMCMC

    136. 

  136. 	cProp=c(nRateProp, nRcatProp, nAlphaProp, nRegimesProp)

    137. 

  137. 	names(cProp)<-names(cOK)<-c("rates", "rcats", "alpha", "sRegimes")

    138. 

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

    139. 

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

    140. 

  140. 	cat("\n")

    141. 

  141. 	print(cProp)

    142. 

  142. 	print(cOK)

    143. 

  143. }

    144. 

  144. 

    145. 

  145. 

    146. 

  146. ## AUXILLIARY FUNCTIONS

    147. 

  147. 

    148. 

  148. lnl.OU <- function(rates, alpha, regimes, apetree, ouch.dat) { # from OUCH 2.7-1:::hansen

    149. 

  149. 	ouch.dat$regimes=regimes

    150. 

  150. 	ouch.tre=apply.BMM.to.tree(rates, apetree)

    151. 

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

    152. 

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

    153. 

  153. 	res=ouch:::ou.lik.fn(ouch.tre, as.matrix(alpha), as.matrix(1), beta, dat)

    154. 

  154. 	list(theta=res$coef, lnL=-0.5*res$deviance)

    155. 

  155. }

    156. 

  156. 

    157. 

  157. lnl.BM <- function(rates, apetree, ouch.dat) { # from OUCH 2.7-1:::brown

    158. 

  158. 	ouch.tre=apply.BMM.to.tree(rates, apetree)

    159. 

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

    160. 

  160. 	nterm <- ouch.tre@nterm

    161. 

  161. 	nchar <- 1

    162. 

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

    163. 

  163. 	b <- ouch.tre@branch.times

    164. 

  164. 	sols <- ouch:::glssoln(w,dat,b)

    165. 

  165. 	e <- sols$residuals # residuals

    166. 

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

    167. 

  167. 	v <- q/nterm

    168. 

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

    169. 

  169. 	list(lnL = -0.5*dev)	

    170. 

  170. }

    171. 

  171. 

    172. 

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

    173. 

  173. 	if(lnR > 0) {

    174. 

  174. 		r=1 

    175. 

  175. 	} else if(lnR < -100) {

    176. 

  176. 		r=0 

    177. 

  177. 	} else {

    178. 

  178. 		r=exp(lnR)

    179. 

  179. 	}

    180. 

  180. 	r

    181. 

  181. }

    182. 

  182. 

    183. 

  183. decide.s.or.m <- function(categories) {

    184. 

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

    185. 

  185. 		return("split")

    186. 

  186. 	}

    187. 

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

    188. 

  188. 		return("merge")

    189. 

  189. 	}

    190. 

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

    191. 

  191. 		return("split")

    192. 

  192. 	}

    193. 

  193. 	else return("merge")

    194. 

  194. }

    195. 

  195. 	

    196. 

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

    197. 

  197. {	

    198. 

  198. 	vv=values

    199. 

  199. 	cc=categories

    200. 

  200. 	new.vv=vv

    201. 

  201. 	new.cc=cc

    202. 

  202. 	nn=length(values)

    203. 

  203. 	

    204. 

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

    205. 

  205. 		ncat <- max(cc)

    206. 

  206. 		m.cc <- sample(1:ncat)[1:2]

    207. 

  207. 		new.cc[new.cc == m.cc[1]] <- m.cc[2]

    208. 

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

    209. 

  209. 		or1 <- vv[cc == m.cc[1]][1]

    210. 

  210. 		or2 <- vv[cc == m.cc[2]][1]

    211. 

  211. 		n1 <- sum(cc == m.cc[1])

    212. 

  212. 		n2 <- sum(cc == m.cc[2])

    213. 

  213. 		nr <- (n1 * or1 + n2 * or2)/(n1 + n2)

    214. 

  214. 		new.vv[cc == m.cc[1] | cc == m.cc[2]] <- nr

    215. 

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

    216. 

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

    217. 

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

    218. 

  218. 			probMerge = 1

    219. 

  219. 		} else {

    220. 

  220. 			probMerge = 0.5

    221. 

  221. 		}

    222. 

  222. 		

    223. 

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

    224. 

  224. 			probSplit = 1

    225. 

  225. 		} else {

    226. 

  226. 			probSplit = 0.5

    227. 

  227. 		}

    228. 

  228. 		factor = probSplit/probMerge

    229. 

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

    230. 

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

    231. 

  231. 	}

    232. 

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

    233. 

  233. 		ncat <- max(cc)

    234. 

  234. 		while (1) {

    235. 

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

    236. 

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

    237. 

  237. 			if (nc > 1) 

    238. 

  238. 			break

    239. 

  239. 		}

    240. 

  240. 		while (1) {

    241. 

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

    242. 

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

    243. 

  243. 			break

    244. 

  244. 		}

    245. 

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

    246. 

  246. 		or <- vv[cc == rcat][1]

    247. 

  247. 		n1 <- sum(new == 1)

    248. 

  248. 		n2 <- sum(new == 2)

    249. 

  249. 		u <- runif(1, min = -0.5 * n1 * or, max = 0.5 * n2 * or)

    250. 

  250. 		nr1 <- or + u/n1

    251. 

  251. 		nr2 <- or - u/n2

    252. 

  252. 		new.vv[new.cc == rcat] <- nr1

    253. 

  253. 		new.vv[new.cc == ncat + 1] <- nr2

    254. 

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

    255. 

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

    256. 

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

    257. 

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

    258. 

  258. 			probSplit = 1

    259. 

  259. 		} else {

    260. 

  260. 			probSplit = 0.5

    261. 

  261. 		}

    262. 

  262. 		

    263. 

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

    264. 

  264. 			probMerge = 1

    265. 

  265. 		} else {

    266. 

  266. 			probMerge = 0.5

    267. 

  267. 		}

    268. 

  268. 		

    269. 

  269. 		factor = probMerge/probSplit

    270. 

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

    271. 

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

    272. 

  272. 	}

    273. 

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

    274. 

  274. }

    275. 

  275. 

    276. 

  276. fix.categories <- function(x)

    277. 

  277. {

    278. 

  278. 	f<-factor(x)

    279. 

  279. 	n<-nlevels(f)

    280. 

  280. 	o<-numeric(n)

    281. 

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

    282. 

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

    283. 

  283. 	nx<-as.numeric(nf)

    284. 

  284. 	nx

    285. 

  285. }

    286. 

  286. 

    287. 

  287. apply.BMM.to.tree <- function(rates, phy) { # SLOW STEP: convert -> deconvert tree

    288. 

  288. 	phy$edge.length=phy$edge.length*rates

    289. 

  289. 	ouch.tre=ape2ouch(phy)

    290. 

  290. 	ouch.tre

    291. 

  291. }

    292. 

  292. 

    293. 

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

    294. 

  294. 	rr=value

    295. 

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

    296. 

  296. 	rr <- rr + rch

    297. 

  297. 	rr

    298. 

  298. }

    299. 

  299. 

    300. 

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

    301. 

  301. 	vv=values

    302. 

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

    303. 

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

    304. 

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

    305. 

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

    306. 

  306. 	nr <- r + rch

    307. 

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

    308. 

  308. 	vv

    309. 

  309. }

    310. 

  310. 

    311. 

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

    312. 

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

    313. 

  313. 	ouch.tre <- ape2ouch(td$phy->ape.tre)		

    314. 

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

    315. 

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

    316. 

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

    317. 

  317. 	ouch.dat$trait=NA

    318. 

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

    319. 

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

    320. 

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

    321. 

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

    322. 

  322. 		}

    323. 

  323. 	}

    324. 

  324. 	return(list(ouch.tre=ouch.tre, ouch.dat=ouch.dat, ape.tre=ape.tre, orig.dat=data))

    325. 

  325. }

    326. 

  326. 

    327. 

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

    328. 

  328. 	index=array(dim=length(startparms))

    329. 

  329. 	for(i in 1:length(startparms)) {

    330. 

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

    331. 

  331. 	}

    332. 

  332. 	if(!all(index)) cOK[which(index==FALSE)]=cOK[which(index==FALSE)]+1

    333. 

  333. 	cOK

    334. 

  334. }

    335. 

  335. 	

    336. 

  336. Stirling2 <- function(n,m)

    337. 

  337. {

    338. 

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

    339. 

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

    340. 

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

    341. 

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

    342. 

  342.     ## ----------------------------------------------------------------

    343. 

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

    344. 

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

    345. 

  345.     ## ----------------------------------------------------------------

    346. 

  346. 

    347. 

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

    348. 

  348.     k <- 0:m

    349. 

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

    350. 

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

    351. 

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

    352. 

  352.     ga <- gamma(k+1)

    353. 

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

    354. 

  354. }
    355.