前面我們一起學習了如何進行多組獨立樣本的差異比較，包括參數和非參數的檢驗方法。但備擇假設為組間的差異不全相等，如拒絕原假設的話，通常需要進行組間差異多重比較，即兩兩比較。針對不同的設計和數據特徵，我們可以選擇合適的方法，如圖1所示。
圖1 多重比較的方法
以下，我們將接著方差分析和非參數檢驗的例子進一步介紹在R中實現組間差異的多重比較。
因為我們要做的是多組獨立樣本參數檢驗後的均數間的兩兩比較，我們先簡單回顧完全隨機設計的方差分析這個例子，詳情可見《方差分析》。
某臨床實驗中心想測試五種藥物降低膽固醇的療效，將50名參與者隨機分成5組，每組分別接受一種藥物，數據資料如表1所示，變量包括使用的藥物類型(trt，分組變量)和膽固醇的降低值（response，反應變量）。前文進行方差分析，結果顯示P＜0.05，拒絕原假設，不能認為五個組的均數是相等。
表1 方差分析的數據
由此，我們需要進一步做均數間的兩兩比較。均數間的兩兩比較檢驗方法有的需要直接使用數據，有的則使用模型方差。下面為大家介紹幾種常見的檢驗方法。
Bonferroni檢驗用途最廣，幾乎可用於任何多重比較的情形，可用agricolae包中LSD.test() 或基礎包中pairwise.t.test()函數實現 。其基本思想是採用α'=α/n作為新的檢驗水準，其中，n為兩兩比較次數, α為累積I類錯誤的概率。
我們首先使用agricolae包中LSD.test()函數，注意將P值的調整方法p.adj設置為"bonferroni"。代碼和結果如下：
library(agricolae)
bon <- LSD.test(aov,"trt", p.adj="bonferroni")
bon$groups #顯示歸類結果
plot(bon)
圖2 LSD.test()函數設置Bonferroni法調整P值的結果
圖3 Groups and Range圖示
圖2和圖3顯示了比較的結果，兩組間有相同的字母表示差異不顯著，兩組間字母不同表示差異顯著。例如本例，E藥物組只有單獨字母a，與其他4組皆不同，代表E藥物組分別與其他4個藥物組的療效有顯著差異。而D藥物組與C藥物組有相同的b,與B藥物組和A藥物組無相同字母，表示D藥物組與C藥物組差異不顯著，D藥物組分別與B藥物組和A藥物組間有顯著差異。
我們還可以用基礎包中的pairwise.t.test()函數做Bonferroni檢驗。代碼和結果如下：
attach(a2)
pairwise.t.test(response,trt, p.adj = "bonferroni")
detach()
圖4 pairwise.t.test()函數運行結果（bonferroni檢驗）
結果如圖4所示，可以看出，結果與圖2、圖3一致。
Bonferroni檢驗比較保守，在組數較多時需要比較的次數較多，校正後的檢驗水準過小，有時很難達到。Holm是一種常用的修正Bonferroni過保守的方法。將p.adj設置為"holm"即可。圖5的結果顯示，相對於校正之前的Bonferroni法，holm的結果沒那麼嚴格，A藥物組與B藥物組、B藥物組與C藥物組、C藥物組與D藥物組間的P值接近檢驗水準。
圖5 pairwise.t.test()函數運行結果（holm法）
LSD-t檢驗，也叫最小顯著性差異，由agricolae包中的LSD.test()函數實現。依據t分布，是t檢驗的簡單變形。LSD檢驗適用於在專業上有特殊意義的樣本均數間的比較，是在設計之初，就已明確要比較某幾個組均數間是否有差異。LSD法側重於減小II類錯誤，但有增大I類錯誤(假陽性)的可能。代碼和結果如下：
library(agricolae)
L <- LSD.test(aov,"trt", p.adj = "none")
L$groups ##顯示歸類結果
圖6 LSD-t檢驗結果
LSD-t檢驗不像Bonferroni檢驗那麼保守，所以結果有所不同。如圖6所示，各藥物組對應的指示字母均不同，表明各藥物組的療效有顯著差異。
Dunnett-t檢驗由multcomp包中glht()函數實現，適用於k-1個試驗組與一個對照組均數差異的多重比較。這種比較也是在設計階段就確定了。
library(multcomp)
D<-glht(aov, linfct = mcp(trt = 'Dunnett'), alternative = 'two.side')
summary(D)
圖7 Dunnett-t檢驗結果
如圖7所示，模型默認第一組為對照組，並依次對比了其他組和它的差異，發現B組與A組的療效差異不顯著，而C、D、E藥物組分別與A藥物組的療效有顯著差異。
Tukey HSD檢驗由基礎包中的Tukey HSD()函數或multcomp包中的glht()函數實現，用於各組樣本均數間的全面比較，適用於組間例數相等的情況。代碼和結果如下：
TukeyHSD(aov)
plot(TukeyHSD(aov))
###或者##
library(multcomp)
T<-glht(aov, linfct = mcp(trt = "Tukey"))
summary(T)
plot(T)
圖8 Tukey檢驗結果1
圖9 效應之差的置信區間
TukeyHSD()和glht()函數的運行結果一致，圖8和圖9結果顯示，除了B-A、C-B、D-C外，其他各組間的差異均具有統計學意義。
Scheffe檢驗由agricolae包中的scheffe.test ()函數實現。各組樣本數相等或不等均可以使用，但是以各組樣本數不相等時使用較多。Scheffe也是通過指示字母的相同與否判定二者間差異是否有統計學意義，不顯示兩兩比較的P值。代碼如下：
library(agricolae)
scheffe.test(aov,"trt", console=TRUE)
SNK-q檢驗由agricolae包中的SNK.test() 函數實現，適用於多個樣本均數兩兩之間的全面比較，與LSD-t檢驗相似，可能存在假陽性。代碼如下：
library(agricolae)
S<-SNK.test(aov,"trt")#console=T表示輸出結果
在數據不滿足正態分布等條件的情況下，組間差異性比較通常會選擇非參數檢驗。同樣的，多組樣本的非參數檢驗也存在無法進行組間兩兩比較的問題。在前文《非參數檢驗》中，我們分析了這樣一個案例：
15隻家兔按體重隨機分為三組，分別注入無菌生理鹽水，0.05μg/g和0.10μg/g劑量的DON進行實驗處理，實驗期滿後測定關節衝洗液中腫瘤壞死因子（TNF-α）的水平（μg/L），數據如表2，比較這三組家兔關節衝洗液TNF-α測定結果差異是否具有統計學意義。
表2 多組獨立樣本數據(data4)
group
TNF-α(μg/L)
0(對照組)
0.218
…
…
1(低劑量組)
0.253
…
…
2(高劑量組)
0.695
…
…
我們採用多組獨立樣本數據的K-W平均秩檢驗，結果顯示拒絕原假設，可以認為這三組家兔關節衝洗液TNF-α測定結果差異有統計學意義。
多組獨立樣本非參數檢驗後組間差異的兩兩比較可以用pgirmess包中的kruskalmc ()函數或PMCMR包中的posthoc.kruskal.nemenyi.test()函數實現。我們先看kruskalmc ()函數。代碼和結果如下：
library(pgirmess)
kruskalmc(data4$TNF,data4$group)
kruskalmc(data4$TNF,data4$group, probs=0.001) ##調整檢驗水準
kruskalmc()函數以邏輯判斷的方式表示是否有顯著性差異，如圖10所示，檢驗水準為0.05時，僅高劑量組和對照組間差異有統計學意義。採用更嚴格的檢驗水準後，各組間不存在顯著性差異。
我們還可以用PMCMR包中的posthoc.kruskal.nemenyi.test()函數實現，代碼和結果如下：
library(PMCMR)
posthoc.kruskal.nemenyi.test(data4$TNF,data4$group)
如圖11所示，結果與圖10檢驗水準為0.05時的結果一致，高劑量組和對照組間差異有統計學意義。
如果上一例中的設計變為隨機區組設計，即如圖12中矩陣數據(行為區組，列為處理組)，如果在Friedman M檢驗後需要兩兩比較，可以用PMCMR包中的posthoc.friedman.nemenyi.test ()函數實現，數據和代碼如下：
library(PMCMR)
posthoc.friedman.nemenyi.test(data4.1)
小結
組間差異多重比較的方法較多，同一種方法在R中又有一到多種實現方式，因此，實際應用中我們需要根據特定的目的、設計和數據的特徵來選擇最適合的分析方法。那麼，看完上面的介紹，對於手上的數據你有一些思路了嗎？如果您需要協助分析，歡迎聯繫我們（關於我們----業務介紹）。
以上案例數據，可在公眾號輸入「多重比較」獲取~
參考來源 ：
1. 孫振球,徐勇勇.醫學統計學:第4版[M].北京:人民衛生出版社.2014.
2. 方積乾,徐勇勇,陳峰,等.衛生統計學:第7版[M].北京:人民衛生出版社.2012.
3. [美]Robert I. Kabacoff著．R語言實戰(第2版) [M].王小寧等譯．北京:人民郵電出出版社.2016.
製作：孟浩蓉、蔡敏
初審：何冠豪、胡建雄
審核：肖建鵬、劉濤
指導：馬文軍
《數據分析和應用》致力於為全國各地公共衛生與醫學工作者（機構）提供專業可靠的統計諮詢、研究設計、數據分析、高通量測序數據和序列分析、調研報告等服務（詳細可見公眾號菜單欄），歡迎有需要的人員和機構與我們聯繫。