第二章：Statistical Modeling

本文参考翻译自：https://web.stanford.edu/class/bios221/book/Chap-Models.html

目录

2.2 The difference between statistical and probabilistic model

在第一章的点位的例子中，已知每个位置的误报率（假阳率false positives）为Bernoulli(0.01)，但被测患者数和确定的蛋白质长度是未知的。在这个例子中，我们假定有50个患者，每个患者在100个位置进行检测（假定是独立的），因此统计完50个病人，我们期望对于任意位置，50个患者中检测出多少个表位这个随机数服从参数为$50\times 0.01$=0.5的泊松分布。这就是我们的零假设。我们可以用数学推导的方法估计出100个位置观测到表位的最大值为7的概率约为$10^{-4}$，这是个很小的概率，然而在数据中观测到了，因此零假设很大概率是不可信的。
现在假设我们知道患者的数量和蛋白质的长度（这些由实验设计给出），但不知道分布律和假阳性率。一旦我们观测到数据，就需要从数据中估算概率模型（比如二项分布，正态分布，泊松分布），以及估计未知的模型参数。这就是本章所要讲的统计推断。

2.3 A simple example of statistical modeling

Start with the data
建模分为两个过程。首先我们需要一个合理的概率模型来对数据过程进行建模。类似我们在第一章中看到的，离散计数数据可以用简单的概率分布（二项分布，多项式分布或泊松分布）。正态分布或钟形曲线通常是连续数据的良好模型。分布也可以是这些最基本的分布的混合。
回顾第一章的表位的那个例子

load("../data/e100.RData")
e99 = e100[-which.max(e100)]

Goodness-of-fit : visual evaluation
我们第一步要从备选分布中求出拟合值；这需要从图形上和从数量上考察拟合优度。对于离散数据，我们绘制频率的条形图（barplot),对于连续变量，我们将查看直方图（histogram)。

barplot(table(e99), space = 0.8, col = "chartreuse4")

但是如果不进行比较，很难决定哪种理论分布更适合数据。一种可视化的拟合优度图称为“rootogram“（Cleveland 1988）。
根图（rootogram）将理论值的平方根显示为红点，而观测值的平方根显示为下拉矩阵。

library("vcd")
gf1 = goodfit( e99, "poisson")
rootogram(gf1, xlab = "", rect_gp = gpar(fill = "chartreuse4"))

如果观测值与理论值完全对应，则框的底部将与水平轴完全对齐。

例如：为了验证对比与一个真正服从参数为0.5的泊松分布，根图是否有差异。做模拟：

simp = rpois(100, lambda = 0.05)
gf2 = goodfit(simp, "poisson")
rootogram(gf2, xlab = "")

通过对比根图的表现，看上去用泊松分布模拟 e99数据很有道理。但别忘了，我们是先移除了outlier的。

Estimating the parameter of the Poisson distribution
泊松分布的均值取为多少的时候才使得数据看上去最有可能是那样的？首先我们来统计一下数据：

table(e100)

table(rpois(100, 3))

对比参数为3的泊松分布，看起来数据中有更多的2和3。所以e100数据不太可能是参数为3的泊松分布。

我们可以一直这样试下去，但我们还可以做得更好。下面我计算一下假如分布为参数为m的泊松分布，观察到该数据的概率有多大：

$P(58\times 0,34\times 1,7\times 2,one7|dataarePoisson(m))$
$=p(0{)}^{58}p(1{)}^{34}p(2{)}^{7}p(7)$
当m=3时，结果为：

prod(dpois(c(0, 1, 2, 7), lambda = 3) ^ (c(58, 34, 7, 1)))

得到的这个概率就是关于参数$\lambda$的似然函数(likelihood function of $\lambda$)，记作:
$L(\lambda ,x=(k_1,k_2,k_3,...))={\prod }_{i=1}^{n}f(k_i)$
比如对于泊松分布，$f(k)=e^{-\lambda }\frac{{\lambda }^k}{k!}$
下面我们通过编程来计算$\lambda$取值从0.05到0.95之间100个数时，对应的泊松分布的似然函数的值。

loglikelihood  =  function(lambda, data = e100) {
   
  sum(log(dpois(data, lambda)))
}

lambdas = seq(0.05, 0.95, length = 100)
loglik = vapply(lambdas, loglikelihood, numeric(1))
plot(lambdas, loglik, type = "l", col = "red", ylab = "", lwd = 2,
     xlab = expression(lambda))
m0 = mean(e100)
abline(v = m0, col = "blue", lwd = 2)
abline(h = loglikelihood(m0), col = "purple", lwd = 2)
m0

实际上在vcd包中，有个函数可以实现一行代码就计算出上述的所有运算。

gf  =  goodfit(e100, "poisson")
names(gf)

gf$par

结果是一样的。

2.3.1 Classical statistics for classical data

我们用数学推导得到在上面例子中$\lambda$的极大似然估计。
前面所讲的是统计推断的第一部分，从手头的数据中估计出分布的未知参数。统计中一个重要的内容是推断数据来自哪一类分布。这就是下面要讲到的内容。
在经典的统计假设框架下，我们为数据考虑一个单一的模型，称为“无效模型（null model)”。“无效模型（null model)”制定了一个“无趣”的基线，例如，所有观察值都来自相同的随机分布，而不管他们来自哪个组或接受的什么治疗。然后我们在这个假设下计算概率，检测是否有“有趣”的事情发生。通常，这就是我们所能做到的最好的事了，因为我们没有足够信息知道“有趣”的或备择模型应该是什么。而在有的时候，我们可以比较两个“竞争”模型，后文将讲到。

另一个有用的方向是回归。我们可能想知道连续的协变量（例如温度）如何影响基于计数的响应变量。你可能已经接触过回归分析，其中响应变量y通过方程y=ax+b+e取决于协变量x，并且残差服从正态分布。而计数数据，相同类型的回归模型是可能的，尽管残差的概率分布需要是非正态的。在那种情况下，我们使用广义线性回归。

2.4 Binomial distributions and maximum likelihood

二项分布有两个参数：n次独立重复的试验；每次试验事件A发生的概率p。

2.4.1 An example

抽取了120个男性数据，对蓝绿色盲进行编码，如果不是色盲，编码为0，否则为1。数据总结如下：
则p的估计量为：
与前面泊松分布的例子一样，如果我们可以计算出关于参数p的许多个不同的值对应的似然函数的值，那么我们就可以绘制出似然函数图并且查看它在哪里跌落的最大。

probs  =  seq(0, 0.3, by = 0.005)
likelihood = dbinom(sum(cb), prob = probs, size = length(cb))
plot(probs, likelihood, pch = 16, xlab = "probability of success",
       ylab = "likelihood", cex=0.6)
probs[which.max(likelihood)]

可以看到最大似然估计并不恰好=1/12。这是因为probs中不包含$1/12\approx 0.0833$这个值。所以我们取的是probs中与0.0833最接近的那个值。

Likelihood for the binomial distribution
举个例子，计算二项分布的似然函数的值，并得到参数的最大似然估计。
不妨设 n=300,观测到 y=40 次成功。
取对数为：

loglikelihood = function(theta, n = 300, k = 40) {
   
  115 + k * log(theta) + (n - k) * log(1 - theta)
}

thetas = seq(0, 1, by = 0.001)
plot(thetas, loglikelihood(thetas), xlab = expression(theta),
  ylab = expression(paste("log f(", theta, " | y)")),type = "l")

2.5 More boxes:multinomial data

2.5.1 DNA count modeling: base pairs

DNA有四个基本分子：A - 腺嘌呤, C - 胞嘧啶, G - 鸟嘌呤, T - 胸腺嘧啶。

2.5.2 Nucleotide bias

这一节在一个真实数例中给出了估计和模拟测试的例子。在staphsequence.ffn.txt中可以找到金黄色葡萄球菌细菌基因的一条DNA链中的数据，我们可以用Bioconductor 包中的 Biostrings读取。

library("Biostrings")
staph = readDNAStringSet("../data/staphsequence.ffn.txt", "fasta")
staph[1]
letterFrequency(staph[[1]], letters = "ACGT", OR = 0)

问题：测试第一个基因的核苷酸是否均匀的分布在四个核苷酸上。
由于它们不同的生理特性，进化选择可以作用于不同的核苷酸频率。因此，我们可以问，来自这一数据的前10个基因是否来自同一多项式分布。我们没有先验的参考，我们只想确定核苷酸是否在前10个基因以相同的概率出现。如果不是这样，这将为我们提供证据，改变这10个基因的选择压力（selective pressure）。

letterFrq = vapply(staph, letterFrequency, FUN.VALUE = numeric(4),
         letters = "ACGT", OR = 0)
colnames(letterFrq) = paste0("gene", seq(along = staph))
tab10 = letterFrq[, 1:10]
computeProportions = function(x) {
    x/sum(x) }
prop10 = apply(tab10, 2, computeProportions)
round(prop10, digits = 2)

p0 = rowMeans(prop10)
p0

假设p0是所有10个基因的多项式分布的概率，并使用Monte Carlo模拟来测试在此假设下观察到的字母频率与期望值之间的偏差是否在合理范围内。
我们通过对p0和每列核苷酸的计数做外积求计数的期望值。

cs = colSums(tab10)
cs

`

expectedtab10 = outer(p0, cs, FUN = "*")
round(expectedtab10)

现在我们可以用rmultinom创建随机表，该表是根据零假设生成的，即认为多项式的概率分布由p0给出。

randomtab10 = sapply(cs, function(s) {
    rmultinom(1, s, p0) } )
all(colSums(randomtab10) == cs)

现在我们将这个重复B=1000次，从每个表中计算得到测试统计量，并将它储存在向量 simulstat中。

stat = function(obsvd, exptd = 20 * pvec) {
   
   sum((obsvd - exptd)^2 / exptd)
}
B = 1000
simulstat = replicate(B, {
   
  randomtab10 = sapply(cs, function(s) {
    rmultinom(1, s, p0) })
  stat(randomtab10, expectedtab10)
})
S1 = stat(tab10, expectedtab10)
sum(simulstat >= S1)

hist(simulstat, col = "lavender", breaks = seq(0, 75, length.out=50))
abline(v = S1, col = "red")
abline(v = quantile(simulstat, probs = c(0.95, 0.99)),
       col = c("darkgreen", "blue"), lty = 2)

从直方图中看到，当零假设成立时，我们看到s1=70.1的概率很小。因此拒绝这十个基因来自同一个多项式分布的假定。

2.6 The ${\chi }^2$ distribution

simulstat统计量的理论分布为参数为30($=10\times (4-1)$）的${\chi }^2$分布，严格来说，simulstat统计量由${\chi }^2$分布近似，表中的计数越大，近似就越好。我们可以据此计算s1=70.1的概率。正如上面看到的那样，蒙特卡洛难以计算小概率：计算的粒度为1/B,因此我们无法估计任何小于该值的概率，实际上，估计的不确定性更大。因此如果有任何理论适用，那将很有用。我们可以用另一个可视化的拟合优度工具QQ图来检测理论值与观测值的匹配程度。当比较两个分布时，无论是来自两个样本的分布，还是分别来自理论上的分布，以及来自样本