主要參考資料：來自帝國理工學院的副研究員 Fabricia Nascimento的教程網站
https://thednainus.wordpress.com/2017/03/03/tutorial-bayesian-mcmc-phylogenetics-using-r/
用到的R代碼都在：https://github.com/thednainus/Bayesian_tutorial
文章中有一些原數據來源于楊子恒老師的Molecular Evolution: A statistic approach一書旬昭，可以去經管之家下載。
MCMC已經被廣泛應用于phylogeny霜旧、divergence time和物種界定等領域萍肆。
這篇教程的主要目標是使用貝葉斯MCMC來評估兩個物種之間的分子距離（記為參數d）命斧，K80模型下核酸轉換/顛換速率的比值（transition/transversion rasio卓起； 記為參數k）尚氛。 核酸數據通過三項分布來計算寡键。對于參數d和k施加gamma分布。
不了解貝葉斯公式和MCMC的朋友請先自行百度大致了解一下
Part1：計算似然值likelihood月幌、先驗prior和后驗分布posterior distribution（un-normalized非常態(tài)化的）（為避免歧義跪腹，下面的部分術語還是用英文）
假設現(xiàn)在有一個人類和大猩猩的pairwise alignment，一共948bp飞醉，兩者之間存在 84 個 transitions 和 6個 transversions。(為了避免數字上存在問題屯阀，下面都使用對數化來計算)
n <- 948 # length of alignment in bp
ns <- 84 # total number of transitions (23+30+10+21)
nv <- 6 # total number of transversions (1+0+2+0+2+1+0+0)
首先缅帘，需要對K80模型有個基本的認識：K80模型只有v和s兩個參數，分別代表顛換率和轉換率难衰，如下表所示钦无，應該有4種轉換和8中顛換（下圖(c)中右上角打錯了，AG之間應該為s盖袭，而不是v）失暂。
log-likelihood函數可以表示為： （在已知距離d和轉換/顛換速率比值k的情況下决记，觀測到數據集D的概率函數）
那么對于數據，在K80模型中倍踪，有如下表示
其中：
別問是怎么得來的系宫，這就是K80模型的內容索昂，要看具體推導可以看1980年kimura的論文
與此對應的，寫成R的代碼應該是：
# Kimura's (1980) likelihood for two sequences
k80.lnL <- function(d, k, n=948, ns=84, nv=6){
p0 <- .25 + .25 * exp(-4*d/(k+2)) +
.5 * exp(-2*d*(k+1)/(k+2))
p1 <- .25 + .25 * exp(-4*d/(k+2)) -
.5 * exp(-2*d*(k+1)/(k+2))
p2 <- .25 - .25 * exp(-4*d/(k+2))
return ((n - ns - nv) * log(p0/4) +
ns * log(p1/4) + nv * log(p2/4))
}
# 這里return的值就是log似然值了扩借，這樣看是不是明白許多
那么進一步椒惨，根據貝葉斯公式可以推導出后驗概率（在觀測到D數據集的情況下，參數d和k的概率分布潮罪，也就是我們最終想求的）：
這里的z是常態(tài)化常數(normalized constant)康谆，
而和
分別是對參數d和k施加的的marginal priors。也就是分別定義d和k概率分布的函數错洁。在這里秉宿，我們對其施加gamma分布:
,
，這兩個參數的prior mean分別為2/20=0.1和2/0.1=20.
常數z是多重積分屯碴，基本上涉及多少個參數就是幾重積分描睦，所以很難解。其實以上各步驟都沒有涉及到MCMC导而，都是貝葉斯的內容忱叭，而正是因為z這個邊緣似然值很難解、錯誤率高今艺，所以才引入了MCMC鏈來避免計算z值韵丑。
所以我們暫時放棄計算z值，那么根據(1)式虚缎，未常態(tài)化(un-normalized)的后驗概率就應該是：
現(xiàn)在就可以在R中把likelihood撵彻、prior和posterior畫出來了，先設定一個參數網格:
dim <- 100 # dimension for the plot
d.v <- seq(from=0, to=0.3, len=dim) # vector of d values
##從0到0.3实牡，長度為100陌僵。假設這兩個物種最大分子距離不超過0.3，也就是序列差異不超過0.3创坞，0-0.3都有可能碗短。
k.v <- seq(from=0, to=100, len=dim) # vector of k values
##從0到100，長度為100题涨。同理偎谁，假設這兩個物種最大轉換/顛換速率比不超過100，0-100都有可能纲堵。
dk <- expand.grid(d=d.v, k=k.v)
##畫一個填滿的網格
par(mfrow=c(1, 3))
#分面板巡雨，1列，每列三個面板
分別繪制三個面板的參數圖：
# Prior surface, f(D)f(k)
Pri <- matrix(dgamma(dk$d, 2, 20) * dgamma(dk$k, 2, .1),
ncol=dim)
##dgamma 密度gamma函數婉支。兩個gamma分布相乘鸯隅，這兩個分布分別就是f(d)和f(k)，這個相乘的結果就是f(d)f(k)，也就是先驗了(prior)蝌以。
image(d.v, k.v, -Pri, las=1, col=heat.colors(50),
main="Prior", xlab="distance, d",
ylab="kappa, k", cex.main=2.0,
cex.lab=1.5, cex.axis=1.5)
contour(d.v, k.v, Pri, nlev=10, drawlab=FALSE, add=TRUE)
# Likelihood surface, f(D|d,k)
lnL <- matrix(k80.lnL(d=dk$d, k=dk$k), ncol=dim)
##把之前設定的k80模型拿來用炕舵，輸入擬定的d和k的向量
# for numerical reasons, we scale the likelihood to be 1
# at the maximum, i.e. we subtract max(lnL)
L <- exp(lnL - max(lnL))
##把log-likelihood轉化為likelihood。這里作者因為數字有點問題的原因跟畅，要把L控制在0-1之間咽筋。
image(d.v, k.v, -L, las=1, col=heat.colors(50),
main="Likelihood", xlab="distance, d",
ylab="kappa, k", cex.main=2.0,
cex.lab=1.5, cex.axis=1.5)
contour(d.v, k.v, L, nlev=10,
drawlab=FALSE, add=TRUE) # creates a contour plot
# Unscaled posterior surface, f(d)f(k)f(D|d,k)
Pos <- Pri * L
##f(d)f(k)就是prior，d(D|d,k)就是likelihood徊件，兩者都是之前算好了的奸攻，直接乘一起就是沒有均一化/或者說常態(tài)化的后驗概率(posterior)
image(d.v, k.v, -Pos, las=1, col=heat.colors(50),
main="Posterior", xlab="distance, d",
ylab="kappa, k", cex.main=2.0,
cex.lab=1.5, cex.axis=1.5)
contour(d.v, k.v, Pos, nlev=10,
drawlab=FALSE, add=TRUE)
一句話：unscaled posteriors=likelihood*priors。likelihood是模型的基礎虱痕，prior影響likelihood睹耐，posterior是影響后的結果。
Part2 馬爾可夫鏈蒙特卡羅方法 Markov chain Monte Carlo (MCMC)
盡管我們得到了posterior distribution部翘，但是因為沒有計算邊緣似然常數z硝训，所以要通過MCMC抽樣來獲得更加準確的posterior distribution。
MCMC算法的基本邏輯：
1新思、對參數d和k設置初始狀態(tài)窖梁。
2、提出一個新的狀態(tài) d*提案（從一個合適的提案密度（proposal density）中得出）夹囚。
3纵刘、對這個新的狀態(tài)的提案進行概率評估:
如果接受了提案（），設定d=d*荸哟，否則d=d
4假哎、現(xiàn)在d進行了一次狀態(tài)檢驗，再對k進行一次鞍历，重復2-3.
5位谋、 儲存目前狀態(tài)的d和k。
6堰燎、回到第二步繼續(xù)進行狀態(tài)的提案。
Gamma分布公式 ：
由于是個常數笋轨，因此在第3步評估狀態(tài)的時候會被消掉秆剪，所以在下面計算prior的時候也不必考慮了。
把part1中的代碼重新梳理一遍：
# function returning the logarithm of the unscaled posterior:
# f(d) * f(k) * f(D|d,k)
# by default we set the priors as:
# f(d) = Gamma(d | 2, 20) and f(k) = Gamma(k | 2, .1)
ulnPf <- function(d, k, n=948, ns=84, nv=6, a.d=2, b.d=20,
a.k=2, b.k=.1) {
# The normalizing constant in the prior densities can
# be ignored:
lnpriord <- (a.d - 1)*log(d) - b.d * d
lnpriork <- (a.k - 1)*log(k) - b.k * k
##分別計算d和k的log-prior
##這里作者寫的是log爵政，但嚴格按照上文gamma分布公式對應應該是ln仅讽。其實無所謂，就是前面多一個常數比例钾挟，最后MCMC鏈里都會消掉的洁灵。
## ulnPf：unnormalized posterior
# log-Likelihood (K80 model):
expd1 <- exp(-4*d/(k+2));
expd2 <- exp(-2*d*(k+1)/(k+2))
##這里作者解釋說宾巍，因為多了一步exp函數，所以運算會比part1快一些
p0 <- .25 + .25 * expd1 + .5 * expd2
p1 <- .25 + .25 * expd1 - .5 * expd2
p2 <- .25 - .25 * expd1
lnL <- ((n - ns - nv) * log(p0/4) +
ns * log(p1/4) + nv * log(p2/4))
# Return unnormalised posterior:
return (lnpriord + lnpriork + lnL)
## 這里就是加毅桃，不是乘珍手，因為是ln
}
終于到了寫MCMC的時候了：
mcmcf <- function(init.d, init.k, N, w.d, w.k) {
# init.d and init.k ：初始d和k的狀態(tài)
# N ：MCMC重復次數
# w.d ：d的sliding-window 提案的寬度
# w.k ：k的sliding-window 提案的寬度
# We keep the visited states (d, k) in sample.d and sample.k
# for easy plotting. In practical MCMC applications, these
# are usually written into a file.
sample.d <- sample.k <- numeric(N+1)
##這就是個計數的作用，看看d和k經歷了多少個循環(huán)双抽，先設置這么多個空百框，后面再填上狀態(tài)
# STEP 1: Initialise the MCMC chain
## 第一步，啟動MCMC鏈
d <- init.d; sample.d[1] <- init.d
k <- init.k; sample.k[1] <- init.k
ulnP <- ulnPf(d, k)
## 這個是上一個block寫好的ulnPf函數牍汹，輸入b和k铐维，輸出log-posterior
acc.d <- 0; acc.k <- 0 # number of acceptances
## 對接受狀態(tài)的次數進行計數
for (i in 1:N) {
# STEP 2: Propose new state d*
# we use a uniform sliding window of width w with reflection
# to propose new values d* and k*
# propose d* and accept/reject the proposal
## 第二步： 提出一個新的狀態(tài)d*，下面是用dprop表示的
## 使用sliding-window來進行狀態(tài)的評估
dprop <- d + runif(1, -w.d/2, w.d/2)
## 通過runif隨機生成一個數值慎菲，這個數值位于設定的滑窗之內嫁蛇，所以不會太大也不會太小。使新的d*狀態(tài)增加這個數值露该。
# reflect if dprop is negative
if (dprop < 0) dprop <- -dprop
## 如果不小心小于0了睬棚，取絕對值
ulnPprop <- ulnPf(dprop, k)
## 計算新的d*提案下的log-posterior數值
lnalpha <- ulnPprop - ulnP
## 判斷是否接受新狀態(tài)。ulnPprop > ulnP才接受狀態(tài)提案有决，所以lnalpha>0接受提案闸拿。
# STEP 3: Accept or reject the proposal
# if ru < alpha accept proposed d*:
## 第三部：判斷是否接受狀態(tài)提案
if (lnalpha > 0 || runif(1) < exp(lnalpha)){
## 這個runif什么作用我沒看懂，難道是允許存在一定的干擾和偏差书幕？
d <- dprop; ulnP <- ulnPprop;
acc.d <- acc.d + 1
}
# else reject and stay where we are (so that nothing needs
# to be done).
# STEP 4: Repeat steps 2-3 for k
# propose k* and accept/reject the proposal
## 第四步：對參數 k 重復同樣的動作
kprop <- k + runif(1, -w.k/2, w.k/2)
# reflect if kprop is negative
if (kprop < 0) kprop <- -kprop
ulnPprop <- ulnPf(d, kprop)
lnalpha <- ulnPprop - ulnP
# if ru < alpha accept proposed k*:
if (lnalpha > 0 || runif(1) < exp(lnalpha)){
k <- kprop; ulnP <- ulnPprop
acc.k <- acc.k + 1
}
# else reject and stay where we are.
# STEP 5: Save chain state
## 第五步： 儲存鏈的狀態(tài)
sample.d[i+1] <- d; sample.k[i+1] <- k
}
# print out the proportion of times
# the proposals were accepted
print("Acceptance proportions (d, k):")
print(c(acc.d/N, acc.k/N))
# return vector of d and k visited during MCMC
return (list(d=sample.d, k=sample.k))
##輸出狀態(tài)列表
}
現(xiàn)在來測試一下mcmc鏈新荤，設置init.d=0.2（初始分子距離為0.2）, init.k=20（初始轉換/顛換速率比kappa為20）, N=1e4（重復10000次）, w.d=0.12（d的window寬為0.12）, w.k（k的window寬為180）：
# Test run-time:
## 可以在自己電腦上測試一下重復一萬次用了多少時間，我這臺Mac pro A1708大概0.13秒台汇。
system.time(mcmcf(0.2, 20, 1e4, .12, 180))
# Run again and save MCMC output:
## 儲存結果
dk.mcmc <- mcmcf(0.2, 20, 1e4, .12, 180)
現(xiàn)在進行trace plot苛骨，看看參數狀態(tài)隨著MCMC鏈的變化：
par(mfrow=c(1,3))
# trace plot of d
plot(dk.mcmc$d, ty='l', xlab="Iteration", ylab="d",
main="Trace of d")
# trace plot of k
plot(dk.mcmc$k, ty='l', xlab="Iteration", ylab="k",
main="Trace of k")
# We can also plot the joint sample of d and k
# (points sampled from posterior surface)
plot(dk.mcmc$d, dk.mcmc$k, pch='.', xlab="d", ylab="k",
main="Joint of d and k")
如圖所示痒芝，MCMC鏈后期的參數混合比較好（圖中看起來變化比較大是因為比例尺和窗口大小設置的問題）。其實可以發(fā)現(xiàn)牵素，d和k的聯(lián)合作圖做出來的和之前我們畫的prior圖基本一致严衬，從公式上來說都是，觀測值比較符合我們預期的概率分布笆呆。右邊log-likelihood也可以看出最佳的后驗概率请琳。
Part3 評估MCMC鏈的效率
作者語：
MCMC鏈是自相關的(autocorrelated)，因為任何一個狀態(tài)要么是前一個狀態(tài)變化而來的赠幕，要么就是前一個狀態(tài)沒變俄精。從直覺上我們會想，一條MCMC鏈的效率與其自相關程度是密不可分的榕堰，它越是自相關竖慧，效率就越低，就要花更多時間去使其收斂以得到較好的后驗概率分布。
現(xiàn)在定義一條MCMC鏈的效率為：
為滯后
時間之后圾旨，兩者之間的自相關系數踱讨。總而言之碳胳，這是一個時間序列的自相關函數勇蝙，兩次觀測時間相隔越近（也是就是
越小）挨约，自相關應該越高味混，相隔越大（也就是
越大），自相關越低诫惭，是一個衰減的序列翁锡。
關于自相關可以來這里補課：
https://blog.csdn.net/yuting_sunshine/article/details/95317735
以下使用R語言的acf函數來計算自相關系數并繪制自相關圖：
# run a very long chain (1e6 generations take about
# 40s in my MacBook Air) to calculate efficiency
dk.mcmc2 <- mcmcf(0.2, 20, 1e7, .12, 180)
# R's acf function (for AutoCorrelation Function)
par(mfrow=c(1,2))
acf(dk.mcmc2$d)
acf(dk.mcmc2$k)
## 使用acf函數計算一組數據的自相關
# Define efficiency function
eff <- function(acf) 1 / (1 + 2 * sum(acf$acf[-1]))
## 定義efficiency函數
# the efficiencies are roughly 22% for d and
# 20% for k:
# [1] 0.2255753 # mcmcf(0.2, 20, 1e7, .12, 180)
eff(acf(dk.mcmc2$d))
# [1] 0.2015054 # mcmcf(0.2, 20, 1e7, .12, 180)
eff(acf(dk.mcmc2$k))
跑出來長這樣：
這張圖表明了隨著時間間隔的延長，兩次計算之間的自相關性下降夕土。下降得越快越好馆衔。再根據效率函數計算出k和d在重復1e7次后的鏈的效率在0.2左右。這意味著這條鏈的抽樣模擬效果幾乎相當于0.2*N個獨立樣本擬合后的計算結果（真的有這么靠譜嗎......）怨绣。作者語：20%效率的鏈就算非常高效了角溃。
這就是計算有效樣本大小(effective sample size, EFF)的公式：
另一方面，由于MCMC鏈的自相關度很高篮撑，在處理結果的時候可以間隔取值减细，以去除高自相關的iteration，這樣又可以節(jié)省磁盤空間又可以保留信息量赢笨∥打颍‘
效率低的鏈
效率低的鏈可能是由于狀態(tài)提案密度較低，或者是由于過度設置復雜參數導致后驗概率復雜茧妒。
什么情況下會出現(xiàn)效率低的鏈呢萧吠？作者聚了兩個例子，我簡單描述一下桐筏。
第一纸型，比如說和
的步長設置得不合適。比如說d作為分子距離梅忌，是在0-1之間的绊袋，而設置3就會使鏈堵住。k作為轉換/顛換速率比铸鹰，應該是比較大的，如果設置每次變動5以內皂岔，就會效率太低（作者原話叫baby-sized steps）蹋笼。如此計算得到這兩條鏈的eff大概為1.5%和0.3%。比如下圖的情況：
作者建議在調整sliding-windows步長的時候，使最后狀態(tài)的接受率在30%左右逊谋。
Part4 如何判斷MCMC鏈已經收斂
這里涉及到一個burn-in的概念擂达，其實很簡單。就是MCMC鏈在收斂的過程中胶滋，最初給定的初始值一般是不合適的板鬓，因此要把初始一定范圍內的接受狀態(tài)去除掉。
舉個例子究恤，這是分別設置較高和較低初始值的結果：
dk.mcmc.l <- mcmcf(0.01, 20, 1e4, .12, 180)
dk.mcmc.h <- mcmcf(0.4, 20, 1e4, .12, 180)
plot(dk.mcmc.l$d, xlim = c(1,200), ylim = c(0,0.4), ty = "l")
lines(dk.mcmc.h$d, col="red")
# We use the low chain to calculate the mean
# and sd of d. We could have used the high chain
# as well.
mean.d <- mean(dk.mcmc.l$d)
sd.d <- sd(dk.mcmc.l$d)
abline(h = mean.d + 2 * c(-sd.d, sd.d), lty = 2)
作者語：判斷MCMC鏈有沒有收斂的最簡單方法就是多跑幾條鏈俭令，然后看他們的95% CI、后驗均值和histogram等是不是相似部宿。
也可以計算后驗均值和標準誤：
# posterior means (similar for efficient chains):
mean(dk.mcmc$d); mean(dk.mcmc.b$d)
mean(dk.mcmc$k); mean(dk.mcmc.b$k)
# posterior means (not so similar for the inefficient chains):
mean(dk.mcmc3$d); mean(dk.mcmc3.b$d)
mean(dk.mcmc3$k); mean(dk.mcmc3.b$k)
## 計算后驗均值
# efficient chain, standard error of the means
sqrt(1/1e4 * var(dk.mcmc$d) / 0.23) # roughly 2.5e-4
sqrt(1/1e4 * var(dk.mcmc$k) / 0.20) # roughly 0.2
# inefficient chain, standard error of the means
sqrt(1/1e4 * var(dk.mcmc3$d) / 0.015) # roughly 9.7e-4
sqrt(1/1e4 * var(dk.mcmc3$k) / 0.003) # roughly 1.6
##計算標準誤，其中var()函數計算的是variance雾叭，即MCMC樣本之間的差值除以樣本量（抽樣次數）
最后：如果模型太過復雜悟耘，一定要多跑幾條鏈并且通過差別比較大的初始值來判斷是否收斂。R語言coda包中的gelman.diag可以診斷MCMC鏈的收斂拷况。
感謝Fabricia Nascimento作煌、Mario dos Reis和Ziheng Yang老師們的教程，花了一整天赚瘦，終于搞懂了
