《Modern Statistics for Modern Biology》Chapter 一： 离散数据模型的预测（1.1 - 1.3）
《Modern Statistics for Modern Biology》Chapter 一： 离散数据模型的预测（1.4 - 1.5）
《Modern Statistics for Modern Biology》Chapter 二: 统计建模（2.1-2.3）
《Modern Statistics for Modern Biology》Chapter 二: 统计建模（2.4-2.5）
《Modern Statistics for Modern Biology》Chapter 二 统计建模（2.6 - 2.7）
从这章开始最开始记录一些markdown等小知识
$\hat{p}=\frac{1}{12}$：
掌握R语言中的apply函数族
卡方检验
Hardy-Weinberg equilibrium( 哈迪-温伯格平衡 )
带你理解beta分布
简单介绍一下R中的几种统计分布及常用模型
- 统计分布每一种分布有四个函数：
d――density（密度函数），p――分布函数，q――分位数函数，r――随机数函数。比如，正态分布的这四个函数为dnorm，pnorm，qnorm，rnorm。下面我们列出各分布后缀，前面加前缀d、p、q或r就构成函数名：norm：正态，t：t分布，f：F分布，chisq：卡方（包括非中心）unif：均匀，exp：指数，weibull：威布尔，gamma：伽玛，beta：贝塔lnorm：对数正态，logis：逻辑分布，cauchy：柯西，binom：二项分布，geom：几何分布，hyper：超几何，nbinom：负二项，pois：泊松signrank：符号秩，wilcox：秩和，tukey：学生化极差
2.8 序列依赖关系建模：马尔可夫链(Modeling sequential dependencies: Markov chains)
如果我们想预测明天的天气，一个合理的猜测是，它最有可能与今天的天气相同，此外，我们还可以说明各种可能的变化的可能性,同样的推理也可以反向应用：我们可以从今天的天气“预测”昨天的天气。这种天气预报方法是马尔科夫假设的一个例子：对明天的预测只取决于今天的情况，而不是昨天或三个星期前的情况(我们可能使用的所有信息都已包含在今天的天气中)。天气的例子也强调，这样的假设不一定是完全正确的，但它应该是一个足够好的假设。将这一假设扩展到前k天的依赖项是相当简单的，其中k是有限的，希望不会太大。马尔科夫假设的实质是这个过程有一个有限的“记忆”，因此预测只需要回顾有限的时间。
我们也可以将其应用于生物序列，而不是时间序列。在DNA中，我们可能会看到特定的序列模式，因此，对核苷酸，称为
digram，例如，CG，CA，CC和CT并不是同样频繁的。例如，在基因组的某些部分，我们看到比我们在独立状态下预期的更频繁的CA实例：我们将序列中的这种依赖关系建模为一个
马尔可夫链
-
其中
N代表任何核苷酸，P(A|C)代表假设前面的碱基是C的情况下A的概率，。图2.14显示了图表上此类转换的示意图。
图2.14：4-状态马尔可夫链的可视化。每个可能的图(例如，CA)的概率是由相应节点之间的边的权重来给出的。例如，CA的概率是由边C→A给出的，我们将在第11章中了解如何使用R包来绘制这类网络图。
2.9 贝叶斯思维
-
到目前为止，我们遵循的是一种经典的方法，即我们分布的参数，即可能的不同结果的概率，代表了长期的频率。这些参数-至少在概念上-是确定的、可知的、固定的数字。我们可能不知道他们，所以我们估计他们从手头的数据。然而，这种方法没有考虑到我们可能已经知道的任何信息，而这些信息可能会约束我们的参数或使某些参数比其他参数更有可能出现，甚至在我们看到任何当前数据集之前。为此，我们需要一种不同的方法，在这种方法中，我们使用概率分布来表示关于参数的知识，并使用数据来更新这些知识，例如通过移动这些分布或使其变得更窄；这是由
贝叶斯范式提供的(图2.15)。
图 2.15 一路上都是海龟。分布参数不确定性的贝叶斯建模是由一个随机变量进行的，随机变量的分布可能依赖于其不确定性可以建模为一个随机变量的参数，这就是分层模型(HierarchicalModels)。 贝叶斯范式是一种实用的方法，使用先验和后验分布作为我们的知识模型，在收集一些数据和进行观察之前和之后。它对于整合或合并来自不同来源的信息特别有用。
假设我们有一个假设，叫它
H，我们想用数据来决定这个假设是不是真的。我们可以把关于H的先验知识形式化化为一个先验概率，即所谓的频率学家写的P(H)，这样的概率是不存在的。他们的观点是，虽然真理是未知的，但在现实中，假设要么是真的，要么是假的；称它为“0.7-true”是没有意义的。当我们看到数据后，我们得到了后验概率。我们把它写成P(H|D)，给出我们看到D的概率H。这可能比P(H)更高或更低，这取决于数据D是什么。
单倍型
为了尽量减少数学形式主义，我们将从一个例子开始。我们研究了一个取证的例子，使用来自Y染色体的组合标记，称为单倍型。
单倍型是一组等位基因(DNA序列变体)的集合，这些等位基因在空间上相邻于一条染色体上，通常是一起遗传的(重组往往不会将它们断开)，因此在遗传上是相互关联的。在这种情况下，我们研究的是Y染色体上的连锁变异。
-
首先我们将研究单倍型频率分析背后的动机，然后我们将重新讨论可能性的概念。在此之后，我们将解释如何将未知参数视为随机数本身，并用先验分布对其不确定性进行建模。然后，我们将看到如何将观测到的新数据合并到概率分布中，并计算有关参数的后验置信度声明。
图2.16：当两个或两个以上核苷酸重复且重复序列彼此直接相邻时，DNA中出现短串联重复序列(STR)。STR也称为微卫星。这种模式的长度可以从2到13个核苷酸，并且重复的数目在不同的个体之间是高度不同的。STR数可用作遗传标记，称为单倍型。
2.9.1 示例：单体型频率
- 我们要估计由一组不同的短串联重复序列(STR)组成的特定Y-单倍型的比例。用于DNA取证的特定STR位置的STR编号组合由特定位置的重复次数标记。美国Y染色体短串联重复序列数据库可通过url网址访问。以下是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个重复，在DXYS156Y位置有12个重复。通过使用这些Y-STR图谱创建的单倍型在同一宗法血统的男性之间共享。由于这些原因，两个不同的男人可能有相同的侧写。我们需要在感兴趣的人群中找到感兴趣的单倍型的潜在比例θ。我们将使用收集到的观测数据，将单倍型的出现视为二项分布中的“成功”.
2.9.2 二项式贝叶斯模型的模拟研究 (Simulation study of the Bayesian paradigm for the binomial)
不是假设我们的参数
θ只有一个值，贝叶斯世界视图允许我们将它看作是从统计分布中得出的结果。该分布表达了我们对参数θ的可能值的信任。原则上，我们可以使用我们喜欢的任何分布，其可能的值对于θ是允许的。当我们看一个表示比例或概率的参数，它的值在0到1之间时，使用β分布是很方便的。密度公式如下：
-
在图2.17中，我们可以看到这个函数是如何依赖于两个参数α和β的，这使得它成为一个非常灵活的发行版系列(因此它可以“适应”许多不同的情况)。它有一个很好的数学性质：如果我们在
θ上有一个β形状的先验信念，观察一个由n个二项试验组成的数据集，然后更新我们的信念，那么θ上的后验分布也会有Beta分布，尽管参数更新了。这是一个数学事实，我们不会证明它，但是我们通过模拟来证明它。
图2.17：α=10、20、50和β=30、60、150的Beta分布用作成功概率的{先驱}。这三种分布具有相同的平均值(α/(α + β))，但围绕平均值的浓度不同。
Y的分布
- 对于给定的
θ的选择，通过方程(2.5)，我们知道Y的分布。但是，如果θ本身也根据某种分布而变化，那么Y的分布是什么呢？我们称之为Y的边缘分布。让我们来模拟一下。首先，我们生成一个10000θs的随机样本。在代码块中，我们再次使用vapply在rtheta的所有元素中应用一个函数，即th的未命名(或“匿名”)函数，以获得另一个长度相同的向量y。然后，对于这些θ中的每一个，我们生成一个Y的随机样本(图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)，验证我们可以得到与上述代码块相同的结果。
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
- 现在，让我们根据
Y=40的结果，计算θ的后验分布。我们将其与理论后验的densPostTheory(我们使用上文第2.4节中定义的概念)进行比较，下面将对其进行更多的讨论。结果如图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个有效数字。
> mean(thetaPostEmp)
[1] 0.1287488
> dtheta = thetas[2]-thetas[1]
> sum(thetas * densPostTheory * dtheta)
[1] 0.1285714
- 为了近似理论密度
densPostTheory的平均值，我们从字面上用数值积分的方法计算了积分，即积分的
sum。这并不总是方便的(或可行的)，特别是如果我们的参数是高维的，即如果我们的模型不仅涉及单个标量θ参数，而且如果θ是高维对象(例如，在图像分析中通常是这样)，并且如果积分不能用解析方法计算。因此，让我们看看如何使用蒙特卡罗积分代替。这与上面的代码类似，在上面的代码中，我们通过调用R的mean函数，使用数值积分从thetaPostEmp计算后验均值。
> thetaPostMC = rbeta(n = 1e6, 90, 610)
> mean(thetaPostMC)
[1] 0.1285824
- 我们可以使用分位数图(QQ-plot，图2.20)检查
MonteCarlo示例thetaPostMC和示例thetaPostEmp之间的一致性。
> qqplot(thetaPostMC, thetaPostEmp, type = "l", asp = 1)
> abline(a = 0, b = 1, col = "blue")
► 问题
- 生成
thetaPostEmp的模拟与导致thetaPostMC的Monte Carlo模拟之间的区别是什么？
后验分布也是β分布 (Posterior distribution is also a beta)
现在我们已经看到，后验分布也是β分布。在我们的例子中，它的参数
α=90和β=610是通过将先前的参数α=50，β=350与观察到的成功y=40和观察到的失败n−y=260相加得到的，从而得到后验参数。
利用后验分布估计θ的不确定性, 我们可以把后验分布最大化的值作为我们的最佳估计，这被称为MAP估计，在这种情况下，它是
。
Suppose we had a second series of data
在看到我们以前的数据后，我们现在有了一个新的先前版本，beta(90，610)。
现在，我们收集了一组新的数据，其中
n=150次观测，y=25次成功，即125次失败。
那么我们对θ的最佳猜测是什么呢？
- 使用与以前相同的推理，新的后验将是：
beta(90+25=115，610+125=735)。这一分布的平均值是，因此对
θ的一个估计是0.135。 - 理论上的最大后验概率(MAP)估计是模型
beta(115, 735)，即。让我们用数字来检查一下。
> 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(更少的峰值)重做所有的计算，这意味着我们使用更少的先验信息。这在多大程度上改变了最终结果？通常情况下，除非后验分布达到极大的峰值，否则先验很少对后验分布进行实质性的改变。如果我们一开始就相当肯定会发生什么，情况就会是这样。有影响的另一种情况是数据很少的情况。
最好的情况是有足够的数据来覆盖先前的数据，这样它的选择就不会对最终的结果产生太大的影响。
比例参数的置信度声明
- 现在是时候总结一下，在给定数据的情况下，这一比例到底在哪里。一个总结是后验可信度区间，它是置信区间的贝叶斯模拟。我们可以用
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
又是一章头大的统计。
