導語：上一期跟大家介紹了假設檢驗的基本思想和計算方法，本期跟大家一起學習拟合優度檢驗和多重假設檢驗（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個差異基因是真實的差異表達基因。