生物統計(4)-單因素方差分析

方差分析的基本思想

在進行科學研究時，有時要按實驗設計將所研究的對象分為多個處理組進行不同的處理，其中處理因素(treatment)至少有兩個水平（level）。這類科研資料的統計分析，是通過所獲得的樣本信息來推斷各處理組均數間的差別是否有統計學意義，即處理是否有影響。常用採用的分析方法就是方差分析（ANOVA，analysis of variance），這是由英國統計學家R.A.Fisher首創，以F命名，故方差分析又稱為F檢驗。

設處理因素有g(g>= 2)個不同水平，實驗對象隨機分為g組，分別接受不同水平的幹預，第i(i=1,2,...,g)組的樣本含量為n_{i}，第i處理組的第j(j=1,2,…ni個觀測值用Xij來表示，其計算結果可能可以整理成以下面的形式，如下所示：

方差分析的目的就是在 成立的條件下，通過分析各處理組均數之間的差別大小，推斷g個總體均數之間有無差別，從面說明處理因素的效應是否存在。

記總均數為

各處理組均數為

總例數為其中，g為處理組數。

實驗數據有三個不同的變異：
1. 總變異。全部觀測值大小不同，這種變異稱為總變異。總變異的大小可能用離均差平方和(sum of squares of deviations from eman,SS)來表示，即各觀測值與總均數X差值的平方和，記為。公式略。

2. 組間變異。各處理組由於接受處理的水平不同，各組的樣本均數也大小不等，這種變異稱為組間變異，其大小用各組均數與總均數的離均差平方和表示，記為SS組間，計算公式略。

各組均數之間相關越懸殊，它們與總均數的差值越在在，就越大，反之就越小。反應了各組均數的變異，存在這種變異的原因有：①隨機誤差；②處理的不同水平可能對實驗結果的影響。

3. 組內變異。在同一處理組中，雖然每個實驗對象接受的處理相同，但觀測值仍各不相同，這種變異稱為組內變異（誤差）。組內變異用組內各觀測值與其所在組的均數的差值的平方和表示，記為，表示隨機誤差的影響。公式略。

總離均差平方和分解為組間離均差平方和與組內離均差平方和，就有了以下公式：

mark

各離均差平方和的自由度為：

mark

變異程序除與離均差平方和的大小有關外，還與其自由度有關，由於各部分自由度不相等，因此積分離均差平方和不能直接比較，須將各部分離均差平方和除以相應的自由度，其比值稱為均方差，簡稱為均方（mean square, MS)。公式為：

mark

如果各組樣本的總體均數相等()，即各處理組的樣本來自相同的總體，無處理因素的作用，則組間變異同組內變異一樣，只反映隨機誤差作用的大小，組間均方與組內無方的比值稱為F統計量，如下所示：

mark

如果F值接近於1，就沒有理由拒絕H0;反之，F值越大，拒絕H0的理由越充分，數理統計的理論證明，當H0成立時，F統計量服從F分布，方差分析是單側F檢驗。

變異是方差分析的基本思想

上面的話可能不太好理解。現在用大白話來理解一下，例如我們要研究某個化合物是否有改善肥胖的效果，我們最初肯定是要做動物實驗，動物實驗的話，例如採用C57的小鼠，分為5組，第1組，給生理鹽水，第2組，給減肥藥（相當於陽性對照），第3組，給高劑量的化合物，第4組，給中劑量的化合物，第5組，給低劑量的化合物。分別餵一段時間後，我們發現小鼠的體重有所變化，這個變化由兩部分構成，第一個就是外界的刺激因素導致的，第二種就是小鼠自身導致的。這種變化我們可以稱為變異（variance），方差分析研究的本質就是這種變異（體重的變化，不是基因變異那種變異），方差分析的英語就是Analysis of variance，如果外界的刺激的因素導致的變異程度遠遠大於小鼠自身因素導致的變異，那麼我們就可以認為，導致小鼠體重這種變異是由外界刺激引起的。

不過這樣還有一個問題，因為數據越多，變異程度就越大，為了解決這個問題，就需要用變異除以自由度（例數-1），這樣比較的就是平均的變異，因此方差分析中就出現了均方（MS）和組內均方的概念。組間均方/組內均方就是通常所說的F值，實際上代表了這樣一個含義：如果組間變異遠遠大於組內變異，那麼組間均方除以組內均方的值肯定很大，反之，這一值就會很小。但是，到底大到什麼程度才認為有統計學意義呢，那就得根據F分布來判斷。

方差分析的應用條件

多個樣本均數比較的方差分析其應用條件為：①各樣本是相互獨立的隨機樣本，均來自**正態分布總體；②相互比較的各樣本的總體方差相等，即具有等方差齊性。

R中的方差分析函數

所用的函數為aov()：
語法為：aov（formula,data=dataframe）
其中formula可使用的特殊符號如下，其中y為因變量，A、B、C為自變量：

mark單因素方差分析

單因素方差分析（one-way ANOVA）是指對單因素試驗結果進行分析，檢驗因素對試驗結果有無顯著性影響的方法。單因素方差分析是用來檢驗多個平均數之間的差異，從而確定因素對試驗結果有無顯著性影響的一種統計方法。對於完全隨機設計試驗且處理數大於2時可以用單因素方差分析（等於2 時用t檢驗）。離差平方和的分解公式為：SST(總和)=SSR(組間)+SSE(組內)，F統計量為MSR/MSE，MSR=SSR/k-1,MSE=SSE/n-k。其中SST為總離差、SSR為組間平方和、SSE為組內平方和或殘差平方和、MSR為組間均方差、MSE為組內均方差。

案例分析

單因素方差分析例4-2：某醫生為了研究一種降血脂新藥的臨床療效，按統一納入標準選擇120名高血脂患者，採用完全隨機設計方法將患者等分為4組（具體分組方法見例4-1），進行雙盲試驗。6周後測得低密度脂蛋白作為試驗結果，見表4-3。問4個處理組患者的低密度脂蛋白含量總體均數有無差別?
數據如下所示：

計算過程如下所示：

1. 導入數據

anova1 <- read.csv("https://raw.githubusercontent.com/20170505a/raw_data/master/data_szq_402.csv",sep=",")
library(reshape2)
anova1 <- melt(anova1)

2.正態檢驗

方差分析需要一定的假設，即數據集應該符合正態和各組的方差相等，可以分別用shapiro.test和bartlett.test檢驗從P值觀察到這兩個假設是符合的。對於不符合假設的情況，我們就要用到非參數方法，例如Kruskal-Wallis秩和檢驗

shapiro.test(anova1$value)

結果如下所示：

P值大於0.05說明數據正態P值大於0.05說明數據正態

3. 方差齊性檢驗:

方差齊性檢驗就是檢驗各組樣本所代表的總體方差是否一致的檢驗，兩樣本方差齊性檢驗使用Bartlett法，同樣，它也適用於多樣本的方差齊性檢驗，它是它要求所檢驗的樣本總體符合正態分頁，當不符合正態分布的時候，就不能使用，則要用Levene檢驗。Levene檢驗不受數據頒販限制，是一種穩健的檢驗，因而被廣泛地認為是一種標準的檢驗方差齊性的檢驗。

方差齊性通常用bartlett檢驗

bartlett.test(anova1$value~anova1$variable)

結果如下所示：

結果顯示p值大於0.05，可認為方差齊性。

或者用levene檢驗：

library(car)
leveneTest(anova1$value~anova1$variable)

結果如下所示：

結果發現sig值大於0.05,表明符合方差齊性假設,可以進行進一步的參數檢驗。

4. 檢驗整體均值是否有差異

result <- aov(value~variable,data=anova1)
summary(result)

結果如下所示：

其中p值小於0.001，因此各組之間存在顯著性差異。另外，R給出的F值是24.88，而書中的例子是24.93，書中的值是由查F表得來的，是個範圍，R中的是具體值。

另外也可以採用oneway.test()函數進行檢驗，如下所示：

> result2 <- oneway.test(value~variable,data=anova1,var.equal=TRUE)
> result2

    One-way analysis of means

data:  value and variable
F = 24.884, num df = 3, denom df = 116, p-value = 1.674e-12

5. 均數間的多重比較

方差分析得出總體之間有差異，要進一步知道哪兩組之間有差異，就要使用均數間的多重比較，常用的比較方法有SNK檢驗（q檢驗），LSD檢驗，Bonferroni檢驗，Dunnett檢驗，TurkeyHSD檢驗。

現在計算各組之間的均數與SD

aggregate(anova1$value,by=list(anova1$variable),FUN=mean)
aggregate(anova1$value,by=list(anova1$variable),FUN=sd)

結果如下所示：

繪製箱形圖可能觀察到不同因素對於因變量的影響

plot.design(value ~ variable,data =anova1, main = "Group means") 
# plot.design was inclued in packages "graphics"

繪製有置信區間的組均值圖

library(gplots)
plotmeans(value ~ variable,data =anova1 )

圖片如下所示：

6. 組間均值的兩兩比較

通過方差分析後，如果整體有差異，則進一步進行兩兩比較，常用的方法有LSD，TukeyHSD，Scheffe檢驗，如下所示：

LSD檢驗

library(agricolae) 
result <- LSD.test(result,"variable",p.adj="bonferroni")
result

結果如下所示：

註：R給出的F值是24.88，而書中的例子是24.93，書中的值是由查F表得來的，是個範圍，R中的是具體值。
Bonferroni校正("bonferroni")，用於多重比較的p值校正，次數不多時適用。如果在同一數據集上同時檢驗n個獨立的假設，那麼用於每一假設的統計顯著水平，應為僅檢驗一個假設時的顯著水平的1/n。舉個例子：如要在同一數據集上檢驗兩個獨立的假設，顯著水平設為常見的0.05。此時用於檢驗該兩個假設應使用更嚴格的0.025。即0.05* (1/2)。該方法是由Carlo Emilio Bonferroni發展的，因此稱Bonferroni校正。在藥理實驗中，動物數通常不會太多，用Bonferroni校正的居多。

TukeyHSD檢驗

result2 <- aov(value~variable,data=anova1)
result.tukey <- TukeyHSD(result2)
result.tukey
plot(result.tukey)

結果如下所示：

圖片結果：

glht()函數做Tukey檢驗

此外，multcomp包中glht()函數提供了多重均值更全面的方法，適用於線性模型和廣義線性模型，下面的代碼重現Tukey HSD檢驗，如下所示：

library(multcomp)
result4 <- aov(value ~ variable, data = anova1)
tukey4 <- glht(result4, linfct=mcp(variable="Tukey"))
summary(tukey4)

計算結果如下所示：

> summary(tukey4)

     Simultaneous Tests for General Linear Hypotheses

Multiple Comparisons of Means: Tukey Contrasts


Fit: aov(formula = value ~ variable, data = anova1)

Linear Hypotheses:
                      Estimate Std. Error t value Pr(>|t|)    
low - control == 0    -0.71500    0.16946  -4.219 0.000291 ***
middle - control == 0 -0.73233    0.16946  -4.322 0.000171 ***
high - control == 0   -1.46400    0.16946  -8.639  < 1e-04 ***
middle - low == 0     -0.01733    0.16946  -0.102 0.999615    
high - low == 0       -0.74900    0.16946  -4.420 0.000108 ***
high - middle == 0    -0.73167    0.16946  -4.318 0.000170 ***
---
Signif. codes:  0 『***』 0.001 『**』 0.01 『*』 0.05 『.』 0.1 『 』 1
(Adjusted p values reported -- single-step method)

繪製不同組的箱線圖，如下所示：

plot(cld(tukey4,level = .05),col="lightgrey")

mark7.殘差分析

這一步做殘差分析就是為了再次驗證原始數據是否服從正態分布，如下所示：

residual <- rstudent(result4)
qqnorm(residual, pch=20, cex=2)
qqline(residual, col="gray60", lwd=2)

殘差的QQ圖如下所示（如果不了解QQ圖，可以參考這篇文章《分位數及其應用》：

mark

shapiro檢驗，如下所示：

> shapiro.test(residual)

    Shapiro-Wilk normality test

data:  residual
W = 0.9884, p-value = 0.4026

p值的為0.4026，可以認為殘差滿足正態性。

繪製殘差圖，如下所示：

plot(residual ~ anova1$variable, main="各組的殘差圖")

mark

對殘差進行方差齊性檢驗，如下所示：

> library(car)
> leveneTest(result4)
Levene's Test for Homogeneity of Variance (center = median)
       Df F value Pr(>F)
group   3   1.493 0.2201
      116

P值大於0.05，可以認為殘差滿足方差齊性。

參考資料

白話統計.馮國雙

醫學統計學.第四版.孫振球