《Modern Statistics for Modern Biology》Chapter 一： 離散數(shù)據(jù)模型的預(yù)測（1.1 - 1.3）
《Modern Statistics for Modern Biology》Chapter 一： 離散數(shù)據(jù)模型的預(yù)測（1.4 - 1.5）
《Modern Statistics for Modern Biology》Chapter 二: 統(tǒng)計建模（2.1-2.3）
《Modern Statistics for Modern Biology》Chapter 二: 統(tǒng)計建模（2.4-2.5）
《Modern Statistics for Modern Biology》Chapter 二 統(tǒng)計建模（2.6 - 2.7）
從這章開始最開始記錄一些markdown等小知識
$\hat{p}=\frac{1}{12}$：
掌握R語言中的apply函數(shù)族
卡方檢驗
Hardy-Weinberg equilibrium( 哈迪-溫伯格平衡 )
帶你理解beta分布
簡單介紹一下R中的幾種統(tǒng)計分布及常用模型
- 統(tǒng)計分布每一種分布有四個函數(shù)：
d――density（密度函數(shù)）挎狸，p――分布函數(shù)龄毡，q――分位數(shù)函數(shù)，r――隨機(jī)數(shù)函數(shù)。比如安接，正態(tài)分布的這四個函數(shù)為dnorm禁筏，pnorm尊流，qnorm作彤，rnorm。下面我們列出各分布后綴昙读，前面加前綴d召调、p、q或r就構(gòu)成函數(shù)名：norm：正態(tài)蛮浑，t：t分布唠叛，f：F分布，chisq：卡方（包括非中心）unif：均勻沮稚，exp：指數(shù)艺沼，weibull：威布爾，gamma：伽瑪蕴掏，beta：貝塔lnorm：對數(shù)正態(tài)障般，logis：邏輯分布，cauchy：柯西盛杰，binom：二項分布挽荡，geom：幾何分布，hyper：超幾何即供，nbinom：負(fù)二項定拟，pois：泊松signrank：符號秩，wilcox：秩和逗嫡，tukey：學(xué)生化極差
2.8 序列依賴關(guān)系建模：馬爾可夫鏈(Modeling sequential dependencies: Markov chains)
如果我們想預(yù)測明天的天氣办素，一個合理的猜測是，它最有可能與今天的天氣相同祸穷，此外，我們還可以說明各種可能的變化的可能性,同樣的推理也可以反向應(yīng)用：我們可以從今天的天氣“預(yù)測”昨天的天氣勺三。這種天氣預(yù)報方法是馬爾科夫假設(shè)的一個例子：對明天的預(yù)測只取決于今天的情況雷滚，而不是昨天或三個星期前的情況(我們可能使用的所有信息都已包含在今天的天氣中)。天氣的例子也強(qiáng)調(diào)吗坚，這樣的假設(shè)不一定是完全正確的祈远，但它應(yīng)該是一個足夠好的假設(shè)呆万。將這一假設(shè)擴(kuò)展到前k天的依賴項是相當(dāng)簡單的，其中k是有限的车份，希望不會太大谋减。馬爾科夫假設(shè)的實質(zhì)是這個過程有一個有限的“記憶”，因此預(yù)測只需要回顧有限的時間扫沼。
我們也可以將其應(yīng)用于生物序列出爹，而不是時間序列。在DNA中缎除，我們可能會看到特定的序列模式严就，因此，對核苷酸器罐，稱為
digram梢为，例如，CG轰坊，CA铸董，CC和CT并不是同樣頻繁的。例如肴沫，在基因組的某些部分粟害，我們看到比我們在獨立狀態(tài)下預(yù)期的更頻繁的CA實例：我們將序列中的這種依賴關(guān)系建模為一個
馬爾可夫鏈
-
其中
N代表任何核苷酸，P(A|C)代表假設(shè)前面的堿基是C的情況下A的概率樊零，我磁。圖2.14顯示了圖表上此類轉(zhuǎn)換的示意圖。
圖2.14：4-狀態(tài)馬爾可夫鏈的可視化驻襟。每個可能的圖(例如夺艰，CA)的概率是由相應(yīng)節(jié)點之間的邊的權(quán)重來給出的。例如沉衣，CA的概率是由邊C→A給出的郁副，我們將在第11章中了解如何使用R包來繪制這類網(wǎng)絡(luò)圖。
2.9 貝葉斯思維
-
到目前為止豌习，我們遵循的是一種經(jīng)典的方法存谎，即我們分布的參數(shù)，即可能的不同結(jié)果的概率肥隆，代表了長期的頻率既荚。這些參數(shù)-至少在概念上-是確定的、可知的栋艳、固定的數(shù)字恰聘。我們可能不知道他們，所以我們估計他們從手頭的數(shù)據(jù)。然而晴叨，這種方法沒有考慮到我們可能已經(jīng)知道的任何信息凿宾，而這些信息可能會約束我們的參數(shù)或使某些參數(shù)比其他參數(shù)更有可能出現(xiàn)，甚至在我們看到任何當(dāng)前數(shù)據(jù)集之前兼蕊。為此初厚，我們需要一種不同的方法，在這種方法中孙技，我們使用概率分布來表示關(guān)于參數(shù)的知識产禾，并使用數(shù)據(jù)來更新這些知識，例如通過移動這些分布或使其變得更窄绪杏；這是由
貝葉斯范式提供的(圖2.15)下愈。
圖 2.15 一路上都是海龜。分布參數(shù)不確定性的貝葉斯建模是由一個隨機(jī)變量進(jìn)行的蕾久，隨機(jī)變量的分布可能依賴于其不確定性可以建模為一個隨機(jī)變量的參數(shù)势似，這就是分層模型(HierarchicalModels)。 貝葉斯范式是一種實用的方法僧著，使用先驗和后驗分布作為我們的知識模型履因，在收集一些數(shù)據(jù)和進(jìn)行觀察之前和之后。它對于整合或合并來自不同來源的信息特別有用盹愚。
假設(shè)我們有一個假設(shè)栅迄，叫它
H，我們想用數(shù)據(jù)來決定這個假設(shè)是不是真的皆怕。我們可以把關(guān)于H的先驗知識形式化化為一個先驗概率毅舆，即所謂的頻率學(xué)家寫的P(H)，這樣的概率是不存在的愈腾。他們的觀點是憋活，雖然真理是未知的，但在現(xiàn)實中虱黄，假設(shè)要么是真的悦即，要么是假的；稱它為“0.7-true”是沒有意義的橱乱。當(dāng)我們看到數(shù)據(jù)后辜梳，我們得到了后驗概率。我們把它寫成P(H|D)泳叠，給出我們看到D的概率H作瞄。這可能比P(H)更高或更低，這取決于數(shù)據(jù)D是什么危纫。
單倍型
為了盡量減少數(shù)學(xué)形式主義宗挥，我們將從一個例子開始节预。我們研究了一個取證的例子，使用來自Y染色體的組合標(biāo)記属韧，稱為單倍型。
單倍型是一組等位基因(DNA序列變體)的集合蛤吓，這些等位基因在空間上相鄰于一條染色體上宵喂，通常是一起遺傳的(重組往往不會將它們斷開)，因此在遺傳上是相互關(guān)聯(lián)的会傲。在這種情況下锅棕，我們研究的是Y染色體上的連鎖變異。
-
首先我們將研究單倍型頻率分析背后的動機(jī)淌山，然后我們將重新討論可能性的概念裸燎。在此之后，我們將解釋如何將未知參數(shù)視為隨機(jī)數(shù)本身泼疑，并用先驗分布對其不確定性進(jìn)行建模德绿。然后，我們將看到如何將觀測到的新數(shù)據(jù)合并到概率分布中退渗，并計算有關(guān)參數(shù)的后驗置信度聲明移稳。
圖2.16：當(dāng)兩個或兩個以上核苷酸重復(fù)且重復(fù)序列彼此直接相鄰時，DNA中出現(xiàn)短串聯(lián)重復(fù)序列(STR)会油。STR也稱為微衛(wèi)星个粱。這種模式的長度可以從2到13個核苷酸，并且重復(fù)的數(shù)目在不同的個體之間是高度不同的翻翩。STR數(shù)可用作遺傳標(biāo)記都许，稱為單倍型。
2.9.1 示例：單體型頻率
- 我們要估計由一組不同的短串聯(lián)重復(fù)序列(STR)組成的特定Y-單倍型的比例嫂冻。用于DNA取證的特定STR位置的STR編號組合由特定位置的重復(fù)次數(shù)標(biāo)記胶征。美國Y染色體短串聯(lián)重復(fù)序列數(shù)據(jù)庫可通過url網(wǎng)址訪問。以下是STR單倍體表的簡短摘錄：
> setwd("e://compute language/NGS/Modern Statistics for Modern Biology/data/")
> haplo6=read.table("haplotype6.txt",header = TRUE)
> haplo6
Individual DYS19 DXYS156Y DYS389m DYS389n DYS389p
1 H1 14 12 4 12 3
2 H3 15 13 4 13 3
3 H4 15 11 5 11 3
4 H5 17 13 4 11 3
5 H7 13 12 5 12 3
6 H8 16 11 5 12 3
這表明單倍型H1在DYS19位置有14個重復(fù)絮吵，在DXYS156Y位置有12個重復(fù)弧烤。通過使用這些Y-STR圖譜創(chuàng)建的單倍型在同一宗法血統(tǒng)的男性之間共享。由于這些原因蹬敲，兩個不同的男人可能有相同的側(cè)寫暇昂。我們需要在感興趣的人群中找到感興趣的單倍型的潛在比例θ。我們將使用收集到的觀測數(shù)據(jù)伴嗡，將單倍型的出現(xiàn)視為二項分布中的“成功”.
2.9.2 二項式貝葉斯模型的模擬研究 (Simulation study of the Bayesian paradigm for the binomial)
不是假設(shè)我們的參數(shù)
θ只有一個值急波，貝葉斯世界視圖允許我們將它看作是從統(tǒng)計分布中得出的結(jié)果。該分布表達(dá)了我們對參數(shù)θ的可能值的信任瘪校。原則上澄暮，我們可以使用我們喜歡的任何分布名段，其可能的值對于θ是允許的。當(dāng)我們看一個表示比例或概率的參數(shù)泣懊，它的值在0到1之間時伸辟，使用β分布是很方便的。密度公式如下：
-
在圖2.17中馍刮，我們可以看到這個函數(shù)是如何依賴于兩個參數(shù)α和β的信夫，這使得它成為一個非常靈活的發(fā)行版系列(因此它可以“適應(yīng)”許多不同的情況)。它有一個很好的數(shù)學(xué)性質(zhì)：如果我們在
θ上有一個β形狀的先驗信念卡啰，觀察一個由n個二項試驗組成的數(shù)據(jù)集静稻，然后更新我們的信念，那么θ上的后驗分布也會有Beta分布匈辱，盡管參數(shù)更新了振湾。這是一個數(shù)學(xué)事實，我們不會證明它亡脸，但是我們通過模擬來證明它押搪。
圖2.17：α=10、20梗掰、50和β=30嵌言、60、150的Beta分布用作成功概率的{先驅(qū)}及穗。這三種分布具有相同的平均值(α/(α + β))摧茴，但圍繞平均值的濃度不同。
Y的分布
- 對于給定的
θ的選擇埂陆，通過方程(2.5)苛白，我們知道Y的分布。但是焚虱，如果θ本身也根據(jù)某種分布而變化购裙，那么Y的分布是什么呢？我們稱之為Y的邊緣分布鹃栽。讓我們來模擬一下躏率。首先，我們生成一個10000θs的隨機(jī)樣本民鼓。在代碼塊中薇芝，我們再次使用vapply在rtheta的所有元素中應(yīng)用一個函數(shù)，即th的未命名(或“匿名”)函數(shù)丰嘉，以獲得另一個長度相同的向量y夯到。然后，對于這些θ中的每一個饮亏，我們生成一個Y的隨機(jī)樣本(圖2.18)耍贾。
方程 2.5
rbeta(n, shape1, shape2)
Random generation for the beta distribution with parameters shape1 and shape2
> rtheta = rbeta(100000, 50, 350)
> y = vapply(rtheta, function(th) {
+ rbinom(1, prob = th, size = 300)
+ }, numeric(1))
> hist(y, breaks = 50, col = "orange", main = "", xlab = "")
? 問題
- 通過使用R的向量化功能并編寫
rbinom(length(rtheta), rtheta, size = 300)阅爽，驗證我們可以得到與上述代碼塊相同的結(jié)果。
rtheta_rep <- rbinom(length(rtheta), rtheta, size = 300)
hist(rtheta_rep, breaks = 50, col = "blue", main = "rbinom", xlab = "")
Histogram of all the thetas such that Y = 40 : the posterior distribution
- 現(xiàn)在付翁，讓我們根據(jù)
Y=40的結(jié)果，計算θ的后驗分布晃听。我們將其與理論后驗的densPostTheory(我們使用上文第2.4節(jié)中定義的概念)進(jìn)行比較胆敞，下面將對其進(jìn)行更多的討論。結(jié)果如圖2.19所示杂伟。
> thetas = seq(0, 1, by = 0.001)
> thetaPostEmp = rtheta[ y == 40 ]
> hist(thetaPostEmp, breaks = 40, col = "chartreuse4", main = "",
+ probability = TRUE, xlab = expression("posterior"~theta))
> densPostTheory = dbeta(thetas, 90, 610)
> lines(thetas, densPostTheory, type="l", lwd = 3)
- 我們還可以檢查上面計算的兩種分布的平均值，并看到它們接近4個有效數(shù)字仍翰。
> mean(thetaPostEmp)
[1] 0.1287488
> dtheta = thetas[2]-thetas[1]
> sum(thetas * densPostTheory * dtheta)
[1] 0.1285714
- 為了近似理論密度
densPostTheory的平均值赫粥，我們從字面上用數(shù)值積分的方法計算了積分，即積分的
sum予借。這并不總是方便的(或可行的)越平，特別是如果我們的參數(shù)是高維的，即如果我們的模型不僅涉及單個標(biāo)量θ參數(shù)灵迫，而且如果θ是高維對象(例如秦叛，在圖像分析中通常是這樣)，并且如果積分不能用解析方法計算瀑粥。因此挣跋，讓我們看看如何使用蒙特卡羅積分代替。這與上面的代碼類似狞换，在上面的代碼中避咆，我們通過調(diào)用R的mean函數(shù)，使用數(shù)值積分從thetaPostEmp計算后驗均值修噪。
> thetaPostMC = rbeta(n = 1e6, 90, 610)
> mean(thetaPostMC)
[1] 0.1285824
- 我們可以使用分位數(shù)圖(QQ-plot查库，圖2.20)檢查
MonteCarlo示例thetaPostMC和示例thetaPostEmp之間的一致性。
> qqplot(thetaPostMC, thetaPostEmp, type = "l", asp = 1)
> abline(a = 0, b = 1, col = "blue")
? 問題
- 生成
thetaPostEmp的模擬與導(dǎo)致thetaPostMC的Monte Carlo模擬之間的區(qū)別是什么黄琼？
后驗分布也是β分布 (Posterior distribution is also a beta)
現(xiàn)在我們已經(jīng)看到樊销，后驗分布也是β分布。在我們的例子中脏款，它的參數(shù)
α=90和β=610是通過將先前的參數(shù)α=50围苫，β=350與觀察到的成功y=40和觀察到的失敗n?y=260相加得到的，從而得到后驗參數(shù)弛矛。
利用后驗分布估計θ的不確定性, 我們可以把后驗分布最大化的值作為我們的最佳估計够吩，這被稱為MAP估計，在這種情況下丈氓，它是
周循。
Suppose we had a second series of data
在看到我們以前的數(shù)據(jù)后强法，我們現(xiàn)在有了一個新的先前版本，beta(90湾笛，610)饮怯。
現(xiàn)在，我們收集了一組新的數(shù)據(jù)嚎研，其中
n=150次觀測蓖墅，y=25次成功，即125次失敗临扮。
那么我們對θ的最佳猜測是什么呢论矾？
- 使用與以前相同的推理，新的后驗將是：
beta(90+25=115杆勇，610+125=735)贪壳。這一分布的平均值是，因此對
θ的一個估計是0.135蚜退。 - 理論上的最大后驗概率(MAP)估計是模型
beta(115, 735)闰靴，即。讓我們用數(shù)字來檢查一下钻注。
> thetas = seq(0, 1, by = 0.001)
> densPost2 = dbeta(thetas, 115, 735)
> mcPost2 = rbeta(1e6, 115, 735)
> sum(thetas * densPost2 * dtheta) # mean, by numeric integration
[1] 0.1352941
> mean(mcPost2) # mean, by MC
[1] 0.1352917
> thetas[which.max(densPost2)] # MAP estimate
[1] 0.134
? 問題
用
softer prior(更少的峰值)重做所有的計算蚂且，這意味著我們使用更少的先驗信息。這在多大程度上改變了最終結(jié)果幅恋？通常情況下杏死，除非后驗分布達(dá)到極大的峰值，否則先驗很少對后驗分布進(jìn)行實質(zhì)性的改變捆交。如果我們一開始就相當(dāng)肯定會發(fā)生什么识埋，情況就會是這樣。有影響的另一種情況是數(shù)據(jù)很少的情況零渐。
最好的情況是有足夠的數(shù)據(jù)來覆蓋先前的數(shù)據(jù)窒舟，這樣它的選擇就不會對最終的結(jié)果產(chǎn)生太大的影響。
比例參數(shù)的置信度聲明
- 現(xiàn)在是時候總結(jié)一下诵盼，在給定數(shù)據(jù)的情況下惠豺，這一比例到底在哪里。一個總結(jié)是后驗可信度區(qū)間风宁，它是置信區(qū)間的貝葉斯模擬洁墙。我們可以用
R取后驗分布的2.5%和97.5%：P(L≤θ≤U)=0.95。
> quantile(mcPost2, c(0.025, 0.975))
2.5% 97.5%
0.1131668 0.1590122
又是一章頭大的統(tǒng)計戒财。
