《Modern Statistics for Modern Bi

2019-03-02 本文已影响9人热衷组培的二货潜

《Modern Statistics for Modern Biology》Chapter 一：离散数据模型的预测（1.1 - 1.3）

意外收获
一个记录了很多R语言函数翻译的博客
 利用R生成隨機分佈的方法

1.4、多维正态分布：DNA

当在模拟中有四种情况时候，例如在图 1.9的box或者当我们计算A、T、C、G四种碱基的数目时候，我们需要扩展二项式模型。
回想一下，当使用二项式时，我们可以考虑两个结果的不相等的概率，方法是将概率p=P(1)=p1分配给结果1，将1−p=p(0)=p0分配给结果0。当像ATCG这样超过两种可能的输出结果时候，我们可以考虑把球扔进不同大小的盒子里，这些盒子的大小对应于不同的概率。我们可以给这些概率贴上p_A，p_C，p_G，p_T的标签。就像在二项式情形下，所有可能结果的概率之和是1：p_A+p_C+p_G+p_T=1。
尝试使用随机数生成器，该生成器通过名为runif的函数生成0到1之间的所有可能的数字。用它生成一个4个level(A，C，G，T)的随机变量，其中p_A= $1\over8$ ，p_C= $3\over8$ ，p_G= $3\over8$ ，p_T= $1\over8$ 。

数学公式. 多项式分布是计数的最重要的模型，R使用一个通式来计算计数的多项向量(x₁,...,x_m)的概率表示具有m个具有概率(p₁,...,p_m)的box的结果:

► 问题 1.7

假设我们有四个同样可能的盒子。使用这个公式，在第一个盒子中观察到4，在第二个盒子中观察到2，而在另外两个盒子中观察不到的概率是多少？

答案：

我们经常进行模拟实验，以检验我们所看到的数据是否与最简单的四箱模型一致，其中每个箱的概率是1/4。在某种意义上，它是稻草人(没有发生什么有趣的事情)。我们将在第2章中看到更多的例子。在这里，我们使用一些R命令来生成这样的计数矢量。首先，假设我们有四种不同的、同样可能的类型的8个字符：

> pvec = rep(1/4, 4)
> pvec
[1] 0.25 0.25 0.25 0.25
> t(rmultinom(1, prob = pvec, size = 8))
     [,1] [,2] [,3] [,4]
[1,]    1    0    4    3

rmultinom

## Usage（用法)
rmultinom(n, size, prob) #【产生随机样本】
dmultinom(x, size = NULL, prob, log = FALSE) # 【密度函数】

## 参数说明
x：长度为k的整数向量，从0到size。
n：要抽取的随机变量的个数
size：整数，指定多维正态试验中被放入K个盒子对象的总数量，对于dmultnom，它默认为sum(x)
prob：长度为K的非负数值向量，指定K类的概率，内部规范为1，无限和缺失值是不允许的
log：逻辑值，如果为真，计算的是对数概率

## 比如 
# https://blog.csdn.net/zhaozhn5/article/details/77933670
rmultinom(n, size, prob) 
#抛 10 次骰子为一次实验，做 1000 次实验。则 n=1000，size=10。 
#prob为每个独立结果出现的概率，其总和为1。 
#结果为 k×n 的矩阵，k即length(prob)

► 问题 1.8：t函数表示转置

► 问题 1.9：rmultinom(n = 8, prob = pvec, size = 1)和 rmultinom(n = 1, prob = pvec, size = 8)的区别? 提示：1.3.1 和 1.3.2.

> rmultinom(n = 8, prob = pvec, size = 1)
     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
[1,]    0    0    0    0    0    0    0    0
[2,]    1    0    0    0    0    0    0    0
[3,]    0    1    0    1    1    1    0    1
[4,]    0    0    1    0    0    0    1    0

> rmultinom(n = 1, prob = pvec, size = 8)
     [,1]
[1,]    3
[2,]    2
[3,]    1
[4,]    2

n = 8, size = 1表示丢1一个球为一次实验，做8次实验，得到的矩阵为length(prob) * n(即4 * 8)的矩阵，每一列代表四种情况中哪一种会成功，标记为1。四行对应的是四种情况，八列代表的是八次试验。
n = 1, size = 8表示丢八个球为，做一次实验，得到一个 4 * 1的矩阵。由于只做一次实验，所以只有一列，每一行表示每种可能出现的次数，最后总次数等于size。
总结见概率分布之rmultinom函数中当n=1和size=1傻傻分不清就好比你每次工作汇报你和你老板的身份

1.4.1 Simulating for power

让我们来看一个例子，用蒙特卡罗来表示多项式，这与科学家们在计划他们的实验时经常要解决的问题有关：我需要多大的样本量呢？。
权重一词在统计学中有特殊的含义。它是发现某事件是否发生的概率，也被称为真正的阳性率。
按照惯例，实验者在计划实验时的目标是80%(或更多)的权重。这意味着，如果相同的实验运行多次，即使本应该发生，大约20%的试验次数将不能产生显著的结果。
零假设H₀，我们收集的DNA数据来自一个公平的过程，在这个过程中，4个碱基中的每一个都有相同的可能(p_A，p_C，p_G，p_T)=(0.25，0.25，0.25，0.25)。这意味着：没有什么有意思的事情发生。这就是我们试图反驳(用统计学家的术语来说，是“拒绝”)的稻草人，所以零假设应该是这样的，偏离它是很有趣的。。
正如您看到的，通过对8个字符和4个等概率结果(由大小相等的方框表示)运行R命令，我们并不总是在每个box中得到2个。这是不可能的，从只看8个字符，核苷酸是否来自一个公平的过程。
让我们来确定，通过观察一个长度为n=20的序列，我们是否能够检测出最初的核苷酸分布是否公平，或者它是否来自其他的(“替代”)过程。
我们使用rmulom函数从零假设中生成了1000个模拟。我们只显示前11列以节省空间。

> obsunder0 = rmultinom(1000, prob = pvec, size = 20)
> dim(obsunder0)
[1]    4 1000
> obsunder0[, 1:11]
     [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10] [,11]
[1,]    4    7    7    8    7    5    7    9    1     6     3
[2,]    3    4    7    4    8    7    4    3    6     3     6
[3,]    9    2    2    4    3    2    5    2    6     3     3
[4,]    4    7    4    4    2    6    4    6    7     8     8

矩阵obsunder0中的每一列都是一个模拟实例。您可以看到，框中的数字差异很大：有些大到8，而期望值是5 = $20\over4$ 。

创建测试

记住：我们知道这些值来自一个公平的过程。显然，仅仅知道流程的期望值是不够的。我们还需要一种可变性的度量，使我们能够描述预期的可变性有多大，以及有多大就是太大了。我们使用以下统计量作为度量，该统计量是观测值与期望值相对于期望值的差值的平方和。因此，对于每个实例，
生成的数据的前三列与我们预期的有多大不同？我们得到：

> expected0 = pvec * 20
> sum((obsunder0[, 1] - expected0)^2 / expected0)
[1] 4.4
> sum((obsunder0[, 2] - expected0)^2 / expected0)
[1] 3.6
> sum((obsunder0[, 3] - expected0)^2 / expected0)
[1] 3.6

度量的值可能会不同-您可以查看3个以上的列，我们将了解如何研究所有1000个列。为了避免重复键入，我们将stat公式(表达式(1.1))封装在一个函数中：

> stat = function(obsvd, exptd = 20 * pvec) {
+   sum((obsvd - exptd)^2 / exptd)
+ }
> stat(obsunder0[, 1])
[1] 4.4

为了更全面地了解这种变化，我们计算了所有1000个实例的度量，并将这些值存储在一个我们称为 S0 的向量中：它包含在H0下生成的值。我们可以考虑图1.10所示的S0值的直方图，这是对空分布的估计。

> S0 = apply(obsunder0, 2, stat)
> summary(S0
+ )
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  0.000   1.200   2.800   3.115   4.400  17.200 
> hist(S0, breaks = 25, col = "lavender", main = "")

summary函数向我们展示了S0具有不同值的分布。例如，从模拟数据中，我们可以近似地得到95%的分位数(这个值将较小的95%的值与5%的最大值分开)。

> q95 = quantile(S0, probs = 0.95)
> q95
95% 
7.6

因此，我们看到5%的S0值大于7.6。我们将提出这作为我们测试数据的临界值，并将拒绝这样的假设，即如果加权平方和统计数据大于7.6，数据来自一个公平的过程，具有同样可能的核苷酸。

Determining our test’s power

我们必须计算我们的检验-基于加权平方和差-将检测到数据实际上不是来自零假设的概率。通过仿真计算了拒绝率。我们从一个由pvecA参数化的可选进程生成1000个模拟实例。

> pvecA = c(3/8, 1/4, 3/12, 1/8)  ## 概率
> observed = rmultinom(1000, prob = pvecA, size = 20)  ## 实验一千次，每次20个球
> dim(observed)
[1]    4 1000
> observed[, 1:7]
     [,1] [,2] [,3] [,4] [,5] [,6] [,7]
[1,]   10   10    8    5    9   12    7
[2,]    2    3    3    5    3    3    6
[3,]    7    4    6    7    4    3    4
[4,]    1    3    3    3    4    2    3
> apply(observed, 1, mean)
[1] 7.388 5.078 5.009 2.525
> expectedA = pvecA * 20
> expectedA
[1] 7.5 5.0 5.0 2.5

与零假设的模拟一样，观测值也有很大的差异。问题是：我们的测试多久(在1000个实例中)检测数据是否偏离NULL？
测试不会拒绝第一个观察结果(10、3、5、2)，因为统计数据的值在95%以内。

> stat(observed[, 1])
[1] 10.8
> S1 = apply(observed, 2, stat)
> q95
95% 
7.6 
> sum(S1 > q95)
[1] 187
> power = mean(S1 > q95)
> power
[1] 0.187

在1000个模拟中，测试确定了187个来自alternative distribution。因此，我们计算出概率P(拒绝H₀|H_A)为0.187。
With a sequence length of n=20 we have a power of about 20% to detect the difference between the fair generating process and our alternative
在实践中，正如我们提到的，一个可接受的加权值（power）是0.8或更多。重复模拟实验，并建议一个新的序列长度n，以确保加权值（power）是可接受的。

经典数据的经典统计

我们不需要用蒙特卡罗模拟数据来计算第95个百分位数，有足够的理论可以帮助我们进行计算。
我们的统计数据实际上有一个众所周知的分布，称为卡方分布(具有3个自由度)，写作 $χ_2 ^3$ 。

这里用markdwon表示为  $χ_2 ^3$

在第2章中，我们将看到如何使用Q-Q图比较分布(参见图2.8)。我们本可以使用更标准的测试，而不是运行手工模拟。然而，我们所学到的过程扩展到许多情况下，卡方分布并不适用。例如，当某些盒子的概率极低，并且它们的计数大多为零时。

1.5、本章节总结

我们使用数学公式和R来计算各种离散事件的概率，我们可以用一些基本的分布来建模：Bernoulli分布是我们最基本的建模-它用于表示单个二进制试验(即结果只有无两种情况，用1和0表示)，如硬币翻转。我们可以将结果编码为0和1。我们称p为成功概率(成功输出1)。
在n个二元试验中，对数字1s的个数采用二项分布(The binomial distribution)，并利用R函数dbinom生成k次成功的概率。我们还看到，我们可以使用函数rbinom来模拟n次试验二项式。
泊松分布(The Poisson distribution)最适合于p较小时(数字1s很少)的情形。它只有一个参数λ，当p较小时，λ=np的泊松分布近似于(n，p)的二项分布。我们使用泊松分布来模拟在沿序列检测表位的试验中随机出现的假阳性数，假设每个位置的假阳性率很小。只要我们知道所有的参数, 我们就能看到了这样一个参数模型如何使我们能够计算极端事件的概率。
多维正态分布（The multinomial distribution）用于具有两个以上可能结果或级别的离散事件。POWER示例向我们展示了如何使用蒙特卡洛(MonteCarlo)模拟来决定我们需要收集多少数据，如果我们想要检验具有等概率的多项模型是否与数据一致的话。我们使用概率分布和概率模型来评估关于我们的数据是如何生成的假设，通过对生成模型的假设。在给定假设的情况下，我们将看到数据的概率称为p值。这和假设是真的概率是不一样的！

课后阅读

The elementary book by Freedman, Pisani, and Purves (1997) provides the best introduction to probability through the type of box models we mention here.
The book by Durbin et al. (1998) covers many useful probability distributions and provides in its appendices a more complete view of the theoretical background in probability theory and its applications to sequences in biology.
Monte Carlo methods are used extensively in modern statistics. Robert and Casella (2009) provides an introduction to these methods using R.
Chapter 6 will cover the subject of hypothesis testing. We also suggest Rice (2006) for more advanced material useful for the type of more advanced probability distributions, beta, gamma, exponentials we often use in data analyses.

参考资源

Freedman, David, Robert Pisani, and Roger Purves. 1997. Statistics. New York, NY: WW Norton.
Durbin, Richard, Sean Eddy, Anders Krogh, and Graeme Mitchison. 1998. Biological Sequence Analysis. Cambridge University Press.
Robert, Christian, and George Casella. 2009. Introducing Monte Carlo Methods with R. Springer Science & Business Media.
Rice, John. 2006. Mathematical Statistics and Data Analysis. Cengage Learning.