导语:上一期跟大家介绍了假设检验的基本思想和计算方法R语言学习笔记(九),本期跟大家一起学习拟合优度检验和多重假设检验(multiple testing)。

01拟合优度检验

➤ 1.1 定义

拟合优度检验用于检验模型对样本观测值的拟合程度,或者说检验样本是否来自一个特定的分布。比如要检验一颗骰子是否是均匀的,那么可以将该骰子抛掷若干次,记录每一面出现的次数,从这些数据出发去检验各面出现的概率是否都是1/6。

➤ 1.2 数学模型

还是以扔骰子的例子为例,我们想检验骰子是否均匀,也就是各面出现的概率是否都是1/6。现在扔60次骰子,记事件Xii点正面向上(i = 1, 2, 3, 4, 5, 6),我们可以提出如下的无效假设:

H0P(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)和理论值进行比较,如果二者足够接近,则认为数据就来自这个分布。

拟合优度检验的假设统计量为:

r语言怎么做平稳性检验(R语言学习笔记十)(1)

该统计量数值越大,证明观测值和理论值的差异越大,越拒绝原假设。该统计量近似服从自由度为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")

r语言怎么做平稳性检验(R语言学习笔记十)(2)

知道了统计量之后,我们就能按照假设检验的步骤往下进行,即根据显著性水平给出拒绝域,然后做出判断。

下面我们对上面两组数据进行卡方检验:

> 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")

r语言怎么做平稳性检验(R语言学习笔记十)(3)

上图中的黑线表示根据不同的α控制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")

r语言怎么做平稳性检验(R语言学习笔记十)(4)

黑线和红线的交叉点为4077,也就是应该拒绝rank < 4077的假设检验,即认为右侧的2923个差异基因是真实的差异表达基因。

,