拓端tecdat|R語(yǔ)言多重比較示例：Bonferroni校正法和Benjamini & Hochberg法
原文鏈接：http://tecdat.cn/?p=21825
原文出處：拓端數(shù)據(jù)部落公眾號(hào)
假設(shè)檢驗(yàn)的基本原理是小概率原理，即我們認(rèn)為小概率事件在一次試驗(yàn)中實(shí)際上不可能發(fā)生。
多重比較的問(wèn)題
當(dāng)同一研究問(wèn)題下進(jìn)行多次假設(shè)檢驗(yàn)時(shí)，不再符合小概率原理所說(shuō)的“一次試驗(yàn)”。如果在該研究問(wèn)題下只要有檢驗(yàn)是陽(yáng)性的，就對(duì)該問(wèn)題下陽(yáng)性結(jié)論的話，對(duì)該問(wèn)題的檢驗(yàn)的犯一類錯(cuò)誤的概率就會(huì)增大。如果同一問(wèn)題下進(jìn)行n次檢驗(yàn)，每次的檢驗(yàn)水準(zhǔn)為α（每次假陽(yáng)性概率為α），則n次檢驗(yàn)至少出現(xiàn)一次假陽(yáng)性的概率會(huì)比α大。假設(shè)每次檢驗(yàn)獨(dú)立的條件下該概率可增加至
常見(jiàn)的多重比較情景包括：
多組間比較
多個(gè)主要指標(biāo)
臨床試驗(yàn)中期中分析
亞組分析
控制多重比較謬誤（Familywise error rate)：Bonferroni矯正
Bonferroni法得到的矯正P值=P×n
Bonferroni法非常簡(jiǎn)單，它的缺點(diǎn)在于非常保守（大概是各種方法中最保守的了），尤其當(dāng)n很大時(shí)，經(jīng)過(guò)Bonferroni法矯正后總的一類錯(cuò)誤可能會(huì)遠(yuǎn)遠(yuǎn)小于既定α。
?
控制錯(cuò)誤發(fā)現(xiàn)率：Benjamini & Hochberg法
簡(jiǎn)稱BH法。首先將各P值從小到大排序，生成順序數(shù)
排第k的矯正P值=P×n/k
另外要保證矯正后的各檢驗(yàn)的P值大小順序不發(fā)生變化。
怎么做檢驗(yàn)
R內(nèi)置了一些方法來(lái)調(diào)整一系列p值，以控制多重比較謬誤（Familywise error rate)或控制錯(cuò)誤發(fā)現(xiàn)率。
Holm、Hochberg、Hommel和Bonferroni方法控制了多重比較謬誤（Familywise error rate)。這些方法試圖限制錯(cuò)誤發(fā)現(xiàn)的概率（I型錯(cuò)誤，在沒(méi)有實(shí)際效果時(shí)錯(cuò)誤地拒絕無(wú)效假設(shè)），因此都是相對(duì)較保守的。
方法BH（Benjamini-Hochberg，與R中的FDR相同）和BY（Benjamini & Yekutieli）控制錯(cuò)誤發(fā)現(xiàn)率，這些方法試圖控制錯(cuò)誤發(fā)現(xiàn)的期望比例。
?
請(qǐng)注意，這些方法只需要調(diào)整p值和要比較的p值的數(shù)量。這與Tukey或Dunnett等方法不同，Tukey和Dunnett也需要基礎(chǔ)數(shù)據(jù)的變異性。Tukey和Dunnett被認(rèn)為是多重比較謬誤（Familywise error rate)方法。
?
要了解這些不同調(diào)整的保守程度，請(qǐng)參閱本文下面的兩個(gè)圖。
?
關(guān)于使用哪種p值調(diào)整度量沒(méi)有明確的建議。一般來(lái)說(shuō)，你應(yīng)該選擇一種你的研究領(lǐng)域熟悉的方法。此外，可能有一些邏輯允許你選擇如何平衡犯I型錯(cuò)誤和犯II型錯(cuò)誤的概率。例如，在一項(xiàng)初步研究中，你可能希望保留盡可能多的顯著值，來(lái)避免在未來(lái)的研究中排除潛在的顯著因素。另一方面，在危及生命并且治療費(fèi)用昂貴的醫(yī)學(xué)研究中，得出一種治療方法優(yōu)于另一種治療方法的結(jié)論之前，你應(yīng)該有很高的把握。
?具有25個(gè)p值的多重比較示例
### --------------------------------------------------------------### 多重比較示例### --------------------------------------------------------------Data = read.table(Input,header=TRUE)
按p值排序數(shù)據(jù)
Data = Data[order(Data$Raw.p),]
檢查數(shù)據(jù)是否按預(yù)期的方式排序
執(zhí)行p值調(diào)整并添加到數(shù)據(jù)框
Data$Bonferroni =????? p.adjust(Data$Raw.p,?????????????? method = "bonferroni")Data$BH =????? p.adjust(Data$Raw.p,?????????????? method = "BH")Data$Holm =????? p.adjust(Data$ Raw.p,?????????????? method = "holm")Data$Hochberg =????? p.adjust(Data$ Raw.p,???????????? ??method = "hochberg")Data$Hommel =????? p.adjust(Data$ Raw.p,?????????????? method = "hommel")Data$BY =????? p.adjust(Data$ Raw.p,?????????????? method = "BY")Data
繪制圖表
plot(X, Y,xlab="原始的p值",ylab="矯正后的P值"lty=1,lwd=2
?
調(diào)整后的p值與原始的p值的圖為一系列的25個(gè)p值。虛線表示一對(duì)一的線。
5個(gè)p值的多重比較示例
### --------------------------------------------------------------### 多重比較示例，假設(shè)示例### --------------------------------------------------------------Data = read.table(Input,header=TRUE)
執(zhí)行p值調(diào)整并添加到數(shù)據(jù)幀
Data$Bonferroni =????? p.adjust(Data$Raw.p,?????????????? method = "bonferroni")Data$BH =????? signif(p.adjust(Data$Raw.p,?????????????? method = "BH"),???????????? 4)Data$Holm =????? p.adjust(Data$ Raw.p,?????????????? method = "holm")Data$Hochberg =????? p.adjust(Data$ Raw.p,?????????????? method = "hochberg")Data$Hommel =????? p.adjust(Data$ Raw.p,?????????????? method = "hommel")Data$BY =????? signif(p.adjust(Data$ Raw.p,?????????????? method = "BY"),???????????? 4)Data
?
繪制(圖表)
plot(X, Y,??????? type="l",
調(diào)整后的p值與原始p值在0到0.1之間的一系列5個(gè)p值的繪圖。請(qǐng)注意，Holm和Hochberg的值與Hommel相同，因此被Hommel隱藏。虛線表示一對(duì)一的線。
最受歡迎的見(jiàn)解
1.Matlab馬爾可夫鏈蒙特卡羅法（MCMC）估計(jì)隨機(jī)波動(dòng)率（SV，Stochastic Volatility） 模型
2.基于R語(yǔ)言的疾病制圖中自適應(yīng)核密度估計(jì)的閾值選擇方法
3.WinBUGS對(duì)多元隨機(jī)波動(dòng)率模型：貝葉斯估計(jì)與模型比較
4.R語(yǔ)言回歸中的hosmer-lemeshow擬合優(yōu)度檢驗(yàn)
5.matlab實(shí)現(xiàn)MCMC的馬爾可夫切換ARMA – GARCH模型估計(jì)
6.R語(yǔ)言區(qū)間數(shù)據(jù)回歸分析
7.R語(yǔ)言WALD檢驗(yàn) VS 似然比檢驗(yàn)
8.python用線性回歸預(yù)測(cè)股票價(jià)格
9.R語(yǔ)言如何在生存分析與Cox回歸中計(jì)算IDI，NRI指標(biāo)
