导语：上一期跟大家介绍了假设检验的基本思想和计算方法R语言学习笔记（九），本期跟大家一起学习拟合优度检验和多重假设检验（multiple testing）。
01
拟合优度检验
➤ 1.1 定义
拟合优度检验用于检验模型对样本观测值的拟合程度，或者说检验样本是否来自一个特定的分布。比如要检验一颗骰子是否是均匀的，那么可以将该骰子抛掷若干次，记录每一面出现的次数，从这些数据出发去检验各面出现的概率是否都是1/6。
➤ 1.2 数学模型
还是以扔骰子的例子为例，我们想检验骰子是否均匀，也就是各面出现的概率是否都是1/6。现在扔60次骰子，记事件Xi为i点正面向上（i = 1, 2, 3, 4, 5, 6），我们可以提出如下的无效假设:
H0：P(Xi) = 1/6，Xi的期望为10
假如我们得到的数据为X1 = 9，X2 = 11，X3 = 12，X4 = 10，X5 = 10，X6 = 8。我们很有可能接受原假设，因为实际数据与期望数据十分接近。
假如我们得到的数据为X1 = 19，X2 = 21，X3 = 12，X4 = 2，X5 = 3，X6 = 3。我们很有可能拒绝原假设，因为实际数据与期望数据相差较大。
那么我们上述推测的数学逻辑是什么呢？其实这背后就是拟合优度假设检验的基本原则。我们需要将根据假设的分布计算出理论的观测值（expected value），然后将实际的观测值（observed value）和理论值进行比较，如果二者足够接近，则认为数据就来自这个分布。
拟合优度检验的假设统计量为：
该统计量数值越大，证明观测值和理论值的差异越大，越拒绝原假设。该统计量近似服从自由度为t-1的卡方分布（χ2），这也正是为什么通常的拟合优度检验都是卡方检验。
➤ 1.3 代码实现
我们首先用R语言简单画一下卡方分布的图像。
> x1 <- seq(0,25,0.5)
> y1 <- dchisq(x1,1)
> y2 <- dchisq(x1,5)
> y3 <- dchisq(x1,15)
> plot(x1,y1,xlab="卡方分布",ylab="Density",type="l",col="red",lwd=2,main="卡方分布图")
> lines(x1,y2,lwd=2,type="l")
> lines(x1,y3,lwd=2,type="l",col="blue")
知道了统计量之后，我们就能按照假设检验的步骤往下进行，即根据显著性水平给出拒绝域，然后做出判断。
下面我们对上面两组数据进行卡方检验：
> ob1 <- c(9,11,12,10,10,8)
> ob2 <- c(19,21,12,2,3,3)
> ex <- c(rep(10, 6))
> chisq.test(cbind(ob1, ex))
Pearson's Chi-squared test
data: cbind(ob1, ex)
X-squared = 0.50429, df = 5, p-value = 0.992
> chisq.test(cbind(ob2, ex))
Pearson's Chi-squared test
data: cbind(ob2, ex)
X-squared = 19.75, df = 5, p-value = 0.001392
可以看到两组数据检验的p值分别为0.992和0.001392，前者接近1，接受原假设，即认为骰子均匀；后者接近0，认为骰子不均匀。
再看一个群体遗传学中的卡方检验例子。假如我们获得了一个突变体材料，现在需要进行遗传分析，在一个F2分离群体有两种表型，野生型和突变体的个体数分别为150和40，现在检验二者是否符合3:1的分离比，我们用R语言计算如下：
> chisq.test(c(150,40), p = c(0.75, 0.25))
Chi-squared test for given probabilities
data: c(150, 40)
X-squared = 1.5789, df = 1, p-value = 0.2089
结果给出的p值为0.2089>0.05，接受原假设，认为符合3:1的分离比，即该突变由隐性单基因控制。
02
多重假设检验
所有的假设检验都有犯错的概率，拿最常见的t检验来说，如果犯第一类错误的概率是0.05，假如进行100次t检验，则至少犯一次错误的概率（这个概率在统计学中有一个专有名词Family Wise Error Rate，FWER）为：
1-(1-0.05)^100 = 0.994
也就是说几乎肯定会犯至少一次错误，因此为了减少犯错的概率，我们在做多重假设检验时，需要对p-value的阈值进行校正。尤其对于生物中的高通量数据，校正是必不可少的步骤。常见的校正方法有两种，下面分别说明。
➤ 2.1 Bonferroni 校正
Bonferroni校正是一种非常严厉的校正方法，方法很简单，直接用单次检验的阈值α除以检验次数n。通常在GWAS文章中，为了控制好假阳性其假设检验的p-value就是经过Bonferroni 校正过的。下面以n = 10000举例说明：
> library(ggplot2)
> library(dplyr)
> m = 10000
> ggplot(tibble(
+ alpha = seq(0, 7e-6, length.out = 100),
+ p = 1 - (1 - alpha)^m),
+ aes(x = alpha, y = p)) + geom_line() +
+xlab(expression(alpha)) +
+ylab("Prob( no false rejection )") +
+ geom_hline(yintercept = 0.05, col = "red") +
+ geom_vline(xintercept = 5.13e-06, col = "blue")
上图中的黑线表示根据不同的α控制FWER时p-value阈值，当取α为0.05时，根据Bonferroni 校正α/n = 5e-06，可以看到上述两个图像的交叉点为5.13 e-06，略大于该值。
➤ 2.2 FDR校正
FDR（False Discovery Rate）校正是一种相对比较温和的校正方法，最常见的做法是对每个p-value做校正，转换为q-value。q=p*n/rank，其中rank是指p-value从小到大排序后的次序，该算法也称为Benjamini-Hochberg（BH）算法。FDR校正在生物信息中应用最多的是通过RNA-seq数据筛选差异表达基因。其具体做法如下：
（1）将所有的p-value从小到大排序，p(1)…p(m)；
（2）对于某个φ值（即FDR值），找到最大的k值，使p(k) ≤ φk/m；
（3）拒绝1到k的假设检验。
下面用代码说明：
> #install.packages(DESeq2)
> #install.package(airway)
> library("DESeq2")
> library("airway")
> data("airway")
> aw = DESeqDataSet(se = airway, design = ~ cell + dex)
> aw = DESeq(aw)
> phi = 0.10
> awde = mutate(awde, rank = rank(pvalue))
> m = nrow(awde)
> ggplot(dplyr::filter(awde, rank <= 7000), aes(x = rank, y = pvalue)) +
+ geom_line() + geom_abline(slope = phi / m, col = "red") +
+ geom_vline(xintercept = 4099, col="blue")
黑线和红线的交叉点为4077，也就是应该拒绝rank < 4077的假设检验，即认为右侧的2923个差异基因是真实的差异表达基因。
转自：“投必得学术”微信公众号
如有侵权，请联系本站删除！
