寫在前面的話，好久沒有寫統(tǒng)計方法的選擇（1）的后續(xù)了，思路都有些打斷了，所以上文所講的當(dāng)時寫完覺得分分鐘可以把后面的寫完，結(jié)果還有些亂了。所以一鼓作氣很重要，不能拖沓，逼迫自己將想法和思維定期更新并形成產(chǎn)品（作品），好了繼續(xù)開始統(tǒng)計學(xué)學(xué)習(xí)。

對數(shù)據(jù)的分布有了了解，那么下一步就要選擇什么樣的方法了。這次先學(xué)習(xí)一下符合正態(tài)分布，并且方差齊性的數(shù)據(jù)如何檢驗。

其實對于統(tǒng)計方法的選擇，在《白話統(tǒng)計》這本書的過程中，給了我對統(tǒng)計學(xué)方法豁然開朗的認(rèn)識。一般線性模型的構(gòu)建，可以視為方差分析和線性回歸的統(tǒng)一，對于這二者的深入分析，可以在后續(xù)中詳細(xì)說明。

2. 參數(shù)檢驗方法

參數(shù)方法就是針對數(shù)據(jù)分布符合正態(tài)分布，并且是方差齊性的數(shù)據(jù)的檢驗的。

2.1 兩組數(shù)據(jù)比較

對于兩組數(shù)據(jù)的比較，那就引出了耳熟能詳?shù)?img class="math-inline" src="https://math.jianshu.com/math?formula=t" alt="t" mathimg="1">檢驗，對于很多臨床醫(yī)師或者醫(yī)學(xué)生，對于兩組之間的比較，沒有分清楚是不是正態(tài)是不是方差齊性，就用t檢驗，不管對錯與否，可以說明t檢驗的知名度還是很高的，那么t檢驗可以針對哪些數(shù)據(jù)進(jìn)行檢驗?zāi)亍?/p>

首先，t檢驗是用于兩組比較，多組比較之后的兩兩比較就不能使用t檢驗，這樣會增加犯Ⅰ類錯誤的風(fēng)險
數(shù)據(jù)如果嚴(yán)重的偏態(tài)分布，就不不能使用t檢驗，這是和我們最開始的流程呼應(yīng)的。那這里提出來，我們就要解答一個問題，為什么嚴(yán)重的偏態(tài)分布就不能使用t檢驗?zāi)兀?/strong>最重要的是，當(dāng)數(shù)據(jù)嚴(yán)重偏態(tài)分布時，均值已經(jīng)不能反應(yīng)數(shù)據(jù)的真實情況，如果仍然使用t檢驗，則會導(dǎo)致結(jié)果的偏倚。而偏態(tài)分布的數(shù)據(jù)，就是我們后面會講到的非參數(shù)檢驗方法。

需要說明的一點就是，也不能狹隘的認(rèn)為t檢驗只是使用于兩組比較，在線性回歸中，也會出現(xiàn)t檢驗，這個時候的t檢驗實質(zhì)上是檢驗回歸系數(shù)是否為0，如果回歸系數(shù)為0，就說明自變量對因變量毫無影響。因此，在線性回歸中看到t檢驗時，要明白t檢驗是一種檢驗的方式或思路。

談起t檢驗，其實最初是由數(shù)學(xué)家戈塞特（Gosset）在Guiness啤酒廠（是的你沒有看錯，就是那個現(xiàn)在還存在的健力士牌黑啤酒）工作，為了檢測啤酒質(zhì)量而發(fā)明了t分布?？墒牵静辉试S員工公開發(fā)表研究成果，于是戈塞特才被迫用筆名發(fā)表了文章?！笇W(xué)生」是發(fā)現(xiàn)這個分布的數(shù)學(xué)家戈塞特（Gosset）的筆名，他于1908年在一個叫Biometrika的雜志上，發(fā)表了關(guān)于t分布的文章，當(dāng)時就是用的這個筆名。其實，戈塞特可不是什么「諾貝爾哥」之類的民科，他在發(fā)表這篇關(guān)于t檢驗的文章之前，曾在現(xiàn)代統(tǒng)計學(xué)的開山鼻祖之一皮爾遜（KarlPearson）的實驗室訪問過一兩年。因此他很好地把基礎(chǔ)研究和實際應(yīng)用結(jié)合了起來，在統(tǒng)計學(xué)的歷史上留下了自己光輝的一頁。

2.1.1 兩組均數(shù)的檢驗

兩組均數(shù)的比較又稱為是成組t及檢驗，最主要的思想就是兩組均數(shù)所代表的兩總體均數(shù)是否不等。
前面談到t檢驗的理論，這里我們直接使用例子來理解t檢驗在R中的操作。

有兩組雌鼠，分別以高蛋白和低蛋白為飼料，8周后記錄各鼠的體重增加量（克），問兩組動物增重的均數(shù)差別是否顯著。
高蛋白組 134 146 104 119 124 161 107 83 113 129 97 123
低蛋白組 70 118 101 85 107 132 94

使用R對其進(jìn)行檢驗

#讀入數(shù)據(jù) a <- c(134, 146, 104, 119, 124, 161, 107, 83, 113, 129, 97, 123) b <- c( 70, 118, 101, 85, 107, 132, 94) #方差齊性檢驗 var.test(a,b) #t檢驗 t.test(a,b,var.equal = T)

查看結(jié)果

> var.test(a,b) F test to compare two variances data: a and b F = 1.0755, num df = 11, denom df = 6, p-value = 0.9788 alternative hypothesis: true ratio of variances is not equal to 1 95 percent confidence interval: 0.198811 4.173718 sample estimates: ratio of variances 1.07552

方差齊性檢驗的結(jié)果 p-value = 0.9788，大于0.05，兩組方差齊性。當(dāng)然方差齊性檢驗也可以使用之前統(tǒng)計方法選擇1所寫的方法檢驗。var.test()函數(shù)進(jìn)行方差檢驗適合于兩個分組，如果分組大于兩組，就要使用bartlett.test()或者leveneTest()函數(shù)進(jìn)行檢驗。結(jié)果不僅給出了p值，而且給出了置信區(qū)間。
再查看t檢驗的結(jié)果

> t.test(a,b,var.equal = T) Two Sample t-test data: a and b t = 1.8914, df = 17, p-value = 0.07573 alternative hypothesis: true difference in means is not equal to 0 95 percent confidence interval: -2.193679 40.193679 sample estimates: mean of x mean of y 120 101

t檢驗的p-value = 0.07573，p值大于0.05，因此不能拒絕 $H_0$ ，兩組均值之間沒有統(tǒng)計學(xué)差異，不能認(rèn)為兩組之間體重增加有差異。t檢驗的結(jié)果不僅給出了置信區(qū)間，而且計算了兩者的均值。好了有了檢驗結(jié)果，那就可視化檢驗結(jié)果。
進(jìn)行可視化，使用到的packages主要是

library(ggpubr)

因為使用的packages是基于ggplot2，所以要將數(shù)據(jù)轉(zhuǎn)換為長數(shù)據(jù)，并將其命名為example，如下

group x 1 1 134 2 1 146 3 1 104 4 1 119 5 1 124 6 1 161 7 1 107 8 1 83 9 1 113 10 1 129 11 1 97 12 1 123 13 2 70 14 2 118 15 2 101 16 2 85 17 2 107 18 2 132 19 2 94

接下來進(jìn)行作圖

p <- ggboxplot(example, x="group", y="x", color = "group", palette ="jco", #用jco配圖，也可以改為“l(fā)ancet”，柳葉刀 add = "jitter") p p+stat_compare_means(method = 't.test')

上述的比較函數(shù)中method = 't.test'，默認(rèn)函數(shù)是method = 'wilcoxon',默認(rèn)方法是wilcoxon秩和檢驗，也就是非參數(shù)檢驗，代碼如下，只是把method去掉，

p <- ggboxplot(example7_10, x="group", y="x", color = "group", palette ="jco", #用jco配圖，也可以改為“l(fā)ancet”，柳葉刀 add = "jitter") p p+stat_compare_means()

做出結(jié)果如下

從兩幅圖看到，t檢驗和wilcoxon檢驗得到的p值是不同的，使用ggpubr package做的圖形，t檢驗得到的p值和使用t.test()函數(shù)得到的結(jié)果接近，但是不是完全一樣，這一點我也稍微有點疑惑，要弄清究竟，需要查看函數(shù)的源代碼，暫且擱置，至少現(xiàn)在兩者都能用。

2.2 多組數(shù)據(jù)之間的比較

多組數(shù)據(jù)之間的比較可以直接使用ANOVA方法。

什么是ANOVA？
ANOVA即就是方差分析，analysis of variance。

2.2.1 ANOVA

方差分析的最基本思想是什么呢，為什么能比較兩組或者多組數(shù)據(jù)之間的差異呢？

方差分析的基本思想就是根據(jù)研究目的和設(shè)計，將總變異中的離均差平方和及其自由度分解為相應(yīng)的若干部分，然后求解相應(yīng)部分的變異，再用各部分的變異與組內(nèi)變異進(jìn)行比較，得出統(tǒng)計量F，最后根據(jù)F值來計算出p值，做出統(tǒng)計推斷。

完全隨機設(shè)計中的方差分詞就是將總變異分解為組間變異和組合兩部分。

還是直接用實例說明問題，并加入操作

為了解燙傷后不同時期切痂對肝臟ATP含量的的影響，將30大鼠隨機分為3組，每組10只，A組為對照組，B組為燙傷后24小時切痂，C組為96小時切痂，全部動物再燙傷后168小時處死并測量肝臟中ATP含量。

數(shù)據(jù)如下

要進(jìn)行數(shù)據(jù)分析需要將以上數(shù)據(jù)寬數(shù)據(jù)轉(zhuǎn)換為長數(shù)據(jù)，具體轉(zhuǎn)換方式可見R 數(shù)據(jù)長寬轉(zhuǎn)換，得到的數(shù)據(jù)為如下形式

group atp 1 1 7.76 2 2 11.14 3 3 10.85 4 1 7.71 5 2 11.60 6 3 8.58

轉(zhuǎn)換之后的數(shù)據(jù)為example
進(jìn)行如下分析

attach(example) fit <- aov(atp~group) summary(fit)

結(jié)果如下

> fit <- aov(atp~group) > summary(fit) Df Sum Sq Mean Sq F value Pr(>F) group 2 119.8 59.92 14.32 5.76e-05 *** Residuals 27 113.0 4.18 --- Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

p值小于0.05 , 拒絕 $H_0$ ，所以認(rèn)為三組之間有差異。
然后進(jìn)行可視化

library(ggpubr) ggboxplot(Example, x="group", y="atp", color = "group", palette = "lancet", add = "jitter")+ stat_compare_means(method = "anova") #設(shè)置method為anova

好的，和之前結(jié)果一致，那么下一步要知道那兩組有差異

2.2.2 查看組間比較

采用ANOVA統(tǒng)計檢驗后，發(fā)現(xiàn)具有統(tǒng)計學(xué)差異，那么具體那兩組之間的差異是怎么樣的，可以采用以下的函數(shù)進(jìn)行查看

> compare_means(atp~group, data=Example8_1) # A tibble: 3 x 8 .y. group1 group2 p p.adj p.format p.signif method <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr> 1 atp 1 2 0.00150 0.0045 0.0015 ** Wilcoxon 2 atp 1 3 0.0524 0.052 0.0524 ns Wilcoxon 3 atp 2 3 0.00209 0.0045 0.0021 ** Wilcoxon

根據(jù)上面的結(jié)果，可以看出組間比較，默認(rèn)方法是Wilcoxon,可以看到第2組和第1組有統(tǒng)計學(xué)差異，第2組和第3組比較有統(tǒng)計學(xué)差異，接下來進(jìn)行組間比較的可視化作圖

#先用非參方法試試 my_comparisons <- list(c("2", "1"), c("1", "3"), c("3", "2")) ggboxplot(Example8_1, x="group", y="atp", color = "group",palette = "lancet")+ #ylim(0,25)+ # Add pairwise comparisons p-value 即是各組兩兩比較的P值，默認(rèn)wilcox.test stat_compare_means(comparisons=my_comparisons)+ # Add global p-value 即是全局比較，這里是3組比較的p值，默認(rèn)Kruskal-Wallis stat_compare_means(label.y = 22)

全局比較采用Kruskal-Wallis，組間比較采用默認(rèn)wilcox.test

默認(rèn)的方法中全局檢驗使用的是Kruskal-Wallis方法，組間比較默認(rèn)使用的是wilcox.test方法，這些都是可以通過改變參數(shù)實現(xiàn)的。
首先查看組間比較的結(jié)果

> compare_means(atp~group, data=Example8_1,method = 't.test') # A tibble: 3 x 8 .y. group1 group2 p p.adj p.format p.signif method <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr> 1 atp 1 2 0.000411 0.00120 0.00041 *** T-test 2 atp 1 3 0.100 0.1 0.10010 ns T-test 3 atp 2 3 0.00325 0.0065 0.00325 ** T-test

并進(jìn)行可視化，全局比較使用ANOVA方法，組間比較使用t檢驗

ggboxplot(Example8_1, x="group", y="atp", color = "group",palette = "lancet")+ #ylim(0,25)+ # Add pairwise comparisons p-value 即是各組兩兩比較的P值 stat_compare_means(comparisons=my_comparisons, method = 't.test')+ # Add global p-value 即是全局比較，這里是3組比較的p值 stat_compare_means(label.y = 22,method = 'anova')

全局比較采用ANOVA，組間比較采用t檢驗

注意：這里使用了t檢驗進(jìn)行組間比較，只是為了說明代碼的作用，以及如何調(diào)整代碼。正如前面提到的組與組之間在全局比較之后，仍然使用t檢驗進(jìn)行兩兩比較，實際上Ⅰ類錯誤的概率會增大。就是統(tǒng)計方法的選擇（1）第一幅圖涉及到的事后比較。

以某一組為基點，進(jìn)行組間比較，并可視化

compare_means(atp~group, data=Example8_1, ref.group = "1", #以1組為參考組 method = "t.test" ) 采用t檢驗的方法

結(jié)果如下

# A tibble: 2 x 8 .y. group1 group2 p p.adj p.format p.signif method <chr> <chr> <chr> <dbl> <dbl> <chr> <chr> <chr> 1 atp 1 2 0.000411 0.00082 0.00041 *** T-test 2 atp 1 3 0.100 0.1 0.10010 ns T-test

作圖展示

#可視化 ggboxplot(Example8_1, x="group", y="atp", color = "group", palette = "lancet", add = 'jitter')+ # Add global p-value，全局比較，設(shè)定比較方法為anova,這個數(shù)據(jù)應(yīng)用anova stat_compare_means(method = "anova", label.y = 30)+ # Pairwise comparison against reference,增加兩兩比較的，這個數(shù)據(jù)應(yīng)用t.test stat_compare_means(label = "p.signif", method = "t.test", ref.group = "1",label.y = 23)

作圖展示中需要注意的是基點的調(diào)節(jié)使用的ref.group這個參數(shù)，也就是reference group。label = "p.signif"參數(shù)設(shè)置主要就是選擇差異結(jié)果的展示。

以第1組為基點進(jìn)行比較

參考組也可以設(shè)置為.all.，即所有的平均值
代碼如下

## 參考組也可以設(shè)置為.all.即所有的平均值 compare_means(atp~group, data=Example8_1, ref.group = ".all.", method = "t.test") #可視化 ggboxplot(Example8_1, x="group", y="atp", color = "group", palette = "lancet")+ # Add global p-value，全局比較，設(shè)定比較方法為anova,這個數(shù)據(jù)應(yīng)用anova stat_compare_means(method = "anova", label.y = 25 )+ #Pairwise comparison against all,這個數(shù)據(jù)應(yīng)用t.test stat_compare_means(label = "p.signif", method = "t.test", ref.group = ".all.",label.y = 18)

參考組為全組平均值

注意：這個例子采用全組平均值，實際上意義不大，只是為了展示這種方法。如果是多組，并且看基因的表達(dá)情況，還是有意義的。

2.2.3 組間比較和事后比較概念了解

compare_means（）的幫助文檔中有這么一句話：

If the grouping variable contains more than two levels, then a pairwise comparison is performed.

什么意思，就是說使用這個函數(shù)，不論多組比較之后，函數(shù)都會自動進(jìn)行兩兩比較，通過上面的示例，大概也看出，函數(shù)自動兩兩比較，參數(shù)方法為t檢驗，非參數(shù)檢驗是wilcox.test，默認(rèn)的方法是wilcox.test，我們可以根據(jù)數(shù)據(jù)的分布狀況更換檢驗方法。但這個與這幅圖后面的方法是一樣嗎，顯然是不一樣的。

事后比較或者事后兩兩比較是指檢驗得知多組比較后，各組均數(shù)不全相等，然后采取相應(yīng)的方法進(jìn)行兩兩比較。而以上我們使用的方法其實是不論多組比較之后，函數(shù)都會自動進(jìn)行兩兩比較。一定程度也是可以的。

在接下來的分析中，參數(shù)檢驗和非參數(shù)檢驗的事后比較可以一起學(xué)習(xí)

參考書籍及文章

R語言統(tǒng)計分析與應(yīng)用-汪海波
白話統(tǒng)計-馮國雙
一張圖說明統(tǒng)計方法的選擇
 R語言中方差齊性檢驗丨數(shù)析學(xué)院
 方差分析與R實現(xiàn)
統(tǒng)計方法如何選以及全代碼作圖實現(xiàn)
想玩轉(zhuǎn)t檢驗？你得從這一篇看起