多組常態分佈資料之差異檢定與事後比較：R的ANOVA與Welch's anova / Parametric Tests for Comparing Many Normal Distribution Groups: ANOVA and Welch's anova in R

1月 28, 2018 Programming/R , Statistics 24 Comments Edit Copy Download

繼前一篇用的Kruskal–Wallis檢定跟Welch's anova來檢定多組非常態分佈資料之間是否有差異的無母數統計之後，這一篇則是用於常態分佈的多組資料之間差異檢定的有母數統計，也就是大家比較熟悉的ANOVA。雖然ANOVA也可以用SPSS來操作，不過這篇用R腳本來實作，不僅自動進行敘述統計、繪製直方圖與箱型圖，還會根據資料同質性檢定結果選擇使用ANOVA或Welch's ANOVA，操作起來更簡單喔。

R腳本與操作環境：RStudio / R script & Environment: RStudio

• R腳本：Parametric Tests for Comparing Many Non-normal Distribution Groups.R
• RStudio：下載RStudio Desktop網頁

R腳本必須要在R環境中執行，我推薦安裝RStudio這個R環境的IDE編輯器。在安裝RStudio的時候也會一併安裝R環境。而這個R腳本使用了許多套件，在R腳本執行過程中會自動確認並安裝。

此腳本主要改編自EV Nordheim、MK Clayton與BS Yandell所著的An Example of ANOVA using R，後面加上的Welch's anova則是來自於An R Companion for the Handbook of Biological Statistics的One-way Anova，事後多重比較Games-Howell檢定則是來自於RPubs。

因為R的版本更新的非常頻繁，你閱讀這篇的時候安裝的不一定是跟我一樣的R 3.4版，可能會有版本差異的問題。若安裝套件時發生了版本不符合的問題，請看Ubuntu中舊版R安裝套件的方法。關於R的版本問題，請看我之前學習R的感想。

那麼以下我們就來開始實作看看吧。

實作：多組常態分佈同質性資料之平均數差異檢定 / Practice: Test for difference in mean among many homoscedastic groups

我們先以最常見的三組常態分佈的資料來進行ANOVA分析。

資料格式 / Data format

• ANOVA example：Google試算表、CSV格式下載

這份資料改編自An Example of ANOVA using R中的例子。這裡面有A、B、C三組資料，我們希望檢測這三組資料之間有沒有差異。我的R腳本只能接受csv格式，第一行欄位必須為group跟value，而資料本身需以逗號(,)隔開。

請下載CSV格式，待會我們就用這個檔案來進行ANOVA分析。

執行R腳本 / Excute R script

下載好R腳本Parametric Tests for Comparing Many Non-normal Distribution Groups.R之後，請在該檔案上按右鍵，在「開啟檔案」裡面選擇RStudio。

在腳本文字的位置按Ctrl + A全選全部的腳本，然後按下上方的「Run」。

選擇你要分析的CSV檔案，例如ANOVA example - data.csv。

第一次執行的時候，因為要安裝大量的R套件，所以會花費大概十分鐘左右的時間。第二次執行的時候，大概只要三十秒左右就可以算完了。如果左下角的console出現了「Finish」，表示分析已經結束。

在CSV相同的資料夾裡面會出現三個分析結果的檔案，各別是：

• 直方圖(histogram)：檔案名稱例子為ANOVA example - data_0128-1957_histogram.png
• 箱型圖(boxplot)：檔案名稱例子為ANOVA example - data_0128-1957_boxplot.png
• 檢定結果：ANOVA example - data_0128-1957_result.txt

以下我們以檢定結果檔案為主，其他兩個統計圖表為輔，用這個來解說ANOVA的分析結果吧。

檢定結果 / Result

在result.txt為後綴的檔案中記錄著分析的過程，以下我就各個段落進行解說：

檔案名稱 / File

### File



C:\data\ANOVA example - data.csv

這是記錄此分析結果是來此於那個CSV檔案。

敘述統計 / Descriptive statistics

### Descriptive statistics


   group n     mean       sd  min   Q1 median    Q3  max

1     A 7 18.61429 1.165373 16.8 17.9   18.8 19.40 20.1

2     B 7 17.65714 1.187234 15.9 16.9   17.7 18.55 19.1

3     C 7 16.84286 1.180194 15.2 16.2   16.7 17.40 18.8

敘述統計會列出各個組別的組別名稱(group)、個數(n)、平均數(mean)、標準差(sd)、最小值(min)、下四分位數(Q1)、中位數(median)、上四分位數(Q3)、最大值(max)。就像我在資料的中心與離度所強調的，在說明敘述統計的時候，最好搭配圖表一起看。

以histogram.png後綴的檔案是各組的直方圖，注意Y軸為所佔比例，而非資料的個數。從這邊大概可以看得出來C組資料偏小，而A組資料偏大，B組則是介於兩者中間。

以boxplot.png後綴的檔案是各組的箱型圖，我們可以看到三組的五數綜合比較。關於箱型圖的閱讀方法，請看用R畫箱型圖。從這邊可以更清楚的看到A、B、C三組的差別，其中A組的下四分位數大於C組的上四分位數，幾乎可以說這兩組有明顯的差異了。

常態性檢定 / Normality test

### Normality test



Group: A, Shapiro-Wilk normality test W = 0.9756555 p-value = 0.9359176

Group: B, Shapiro-Wilk normality test W = 0.9483481 p-value = 0.7146411

Group: C, Shapiro-Wilk normality test W = 0.9878776 p-value = 0.9886119

這裡測試各組的資料分佈是否為常態分佈。根據柴惠敏的建議，若該組的樣本個數小於50，則是使用Shapiro-Wilk常態檢定，檢定統計量為W值；若該組的樣本數在50以上，則是使用Kolmogorov-Smirnov (K-S)常態檢定，檢定統計量為D值。STHDA推薦使用Shapiro-Wilk檢定，因為它的統計檢定力(power)高於Kolmogorov-Smirnov檢定。這邊R腳本會依照該組的樣本個數來自動判斷要使用的檢定法。

這兩個檢定的虛無假設都是該組資料為常態分佈。因此若p值(p-value)在0.05以下，則表示該組資料並非常態分佈。以這個例子來說，三組資料的p值都大於顯著水準0.05，意思是三組資料皆為常態分佈。

大部分做ANOVA分析時都不會特別去檢定資料分佈是否非常態。這個常態性檢定僅供參考即可。

變異數同質性檢定 / Test for Homogeneity of Variance

### Test for Homogeneity of Variance



Levene's Test for Homogeneity of Variance (center = mean)
       Df F value Pr(>F)

group  2  0.0162  0.984
       18               



Data are homoscedastic. Excute ANOVA.

這邊我們使用Levene變異數同質性檢定中無母數檢定的方法來檢定三組的變異數是否相同。請看「Pr(>F)」這一欄的顯著性是否在0.05以上。若在0.05以上使用的是ANOVA，否則使用Welch's anova。

在本例中，Pr(>F)為0.984，大於顯著水準0.05，表示三組資料變異數同質，因此進行ANOVA分析。

值得一提的是，這邊跟我在前一篇Kruskal–Wallis檢定中不太一樣。這裡的Levene變異數同質性檢定的資料中心(center)採用的是平均值(mean)，做法跟SPSS相同。而Kruskal–Wallis檢定是無母數檢定，Levene檢定的資料中心採用的是中位數(median)。

單因子變異數分析 / ANOVA

### ANOVA for equal variances



Analysis of Variance Table



Response: value
           Df Sum Sq Mean Sq F value  Pr(>F)  

group      2 11.007  5.5033  3.9683 0.03735 *

Residuals 18 24.963  1.3868                  

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1



Eta squared: 0.3059998

ANOVA檢定統計量為F值(F value)，此例的F值為3.9683，顯著性p值要看Pr(>F)，本例子中F值的顯著性p值為0.03735，低於顯著水準0.05，表示三組資料的平均數有顯著的差異。

最後的Eta squared是表示分組變項group對於value的影響力，也就是關聯強度。計算方式是用的R的lsr套件的etaSquared函數。關於Eta squared的說明請看我介紹ANOVA的這篇。Eta squared的判斷準則如下：

• 0.059 > eta squared >= 0.01 低度關聯強度
• 0.138 > eta squared >= 0.059 中度關聯強度
• eta squared >= 0.138 高度關聯強度

以本例的0.3059998為例，是屬於高度關聯強度。

Scheffe法事後多重比較檢定 / Scheffe Post Hot test

### Post Hoc Tests


   Posthoc multiple comparisons of means : Scheffe Test 
     95% family-wise confidence level



$group
           diff    lwr.ci      upr.ci   pval    

B-A -0.9571429 -2.635501  0.72121521 0.3370    

C-A -1.7714286 -3.449787 -0.09307051 0.0376 *  

C-B -0.8142857 -2.492644  0.86407235 0.4493    



---

Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

跟我之前介紹ANOVA一樣，我們用Scheffe法來做事後多重比較檢定。下面會列出三個組別A、B與C之間的比較。其中「C-A」組別的差異-1.7714286的顯著性pval為0.0376，低於顯著水準0.05，表示C組跟A組之間有顯著的差異。其中差異值-1.7714286表示前者「C組」小於後者「A組」。

雖然我也想跟前一篇ANOVA一樣做統計檢定力分析(power analysis)，不過用R來做的話看起來還得花一番功夫，在這裡還是先打住吧。

這樣我們就做完了ANOVA跟Scheffe事後多重比較囉。

實作：多組常態分佈異質性資料之中位數差異檢定 / Practice: Test for difference in means among many heteroscedastic groups

在我之前介紹ANOVA中，若各組之間的資料違反了同質性檢定，我只有在事後多重比較改用了Dunnett T3檢定。但後來看Salvatore S. Mangiafico的建議是可以改用Welch's anova，以及Jim Frost推薦使用適合Welch's anova的Games-Howell事後多重比較。以下我們就來看看當各組資料違反同質性檢定時是會如何分析吧。

異質性的資料 / Heteroscedastic data

• color-anova-example：Google試算表、CSV格式下載

這份資料是來自Jean and Alexander Heard Library的Conducting an ANOVA using RStudio中的CSV file with red, green, and blue ANOVA data。該教學也是用R來做ANOVA分析，但是其實它的三組資料red、green以及blue違反了變異數同質性檢定，因此照該網頁的教學來分析其實不太適合，但正好給我們作為Welch's anova的例子。

請下載CSV格式，待會我們就用這個檔案來進行Welch's anova分析。

執行R腳本 / Excute R script

操作方法跟上述相同，在此我們就不贅述。若順利分析完成，最後一樣可以獲得三個檔案：

• 直方圖：color-anova-example - data_0128-2115_histogram.png
• 箱型圖：color-anova-example - data_0128-2115_boxplot.png
• 檢定結果：color-anova-example - data_0128-2115_result.txt

以下我們依然是來從檢定結果檔案來說明分析的結果吧。

檢定結果 / Result

在result.txt為後綴的檔案中記錄著分析的過程，以下我就各個段落進行解說：

檔案名稱 / File

### File



C:\data\color-anova-example - data.csv

這是本次分析的檔案。

敘述統計 / Descriptive statistics

### Descriptive statistics


   group  n      mean       sd  min    Q1 median    Q3  max

1   red 24  2.491667 1.546350 0.80 1.405    2.0  3.40  6.7

2 green 24  8.530417 6.977992 1.00 4.000    6.4  9.85 27.2

3  blue 24 10.632083 5.976309 5.27 6.600    8.0 11.75 25.6

敘述統計中可以看到這次的資料有三組「red」、「green」跟「blue」，三組的平均數(mean)、標準差(sd)都不一樣，可以想象他們應該會有很大的差異。

在以histogram.png後綴的直方圖裡面可以看到三組的分佈。red較為集中，而green跟blue都非常分散，green甚至有較明顯的雙峰分佈。

在以boxplot.png後綴的箱型圖裡面可以清楚比較三組的分佈。可以看到red明顯的小於blue，而green分佈較廣，介於兩組之間。

常態性檢定 / Normality test

### Normality test



Group: red, Shapiro-Wilk normality test W = 0.865417 p-value = 0.004283695

Group: green, Shapiro-Wilk normality test W = 0.8609989 p-value = 0.003522477

Group: blue, Shapiro-Wilk normality test W = 0.7905905 p-value = 0.0002075442

在常態性檢定中，因為三組資料的樣本個數都小於50，因此使用Shapiro-Wilk檢定。而這三組的p值都小於顯著水準0.05，表示三組資料皆非常態分佈。

一般來說，這時候應該改用Kruskal–Wallis檢定，但因為下面可以看到這三組也違反了同質性假設，因此仍會用Welch's anova來做分析。

變異數同質性檢定 / Test for Homogeneity of Variance

### Test for Homogeneity of Variance



Levene's Test for Homogeneity of Variance (center = mean)
       Df F value    Pr(>F)    

group  2  9.9731 0.0001568 ***
       69                      

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1



Data are heteroscedastic. Excute Welch's anova.

在Levene檢定中，Pr(>F)的顯著性p值0.0001568小於顯著水準0.05，表示各組資料違反了變異數同質性的假設，為異質性，因此接下來將使用Welch's anova進行分析。

Welch's anova

### Welch’s anova for unequal variances


   One-way analysis of means (not assuming equal variances)



data:  value and group

F = 27.526, num df = 2.000, denom df = 33.916, p-value = 7.898e-08

在Welch's anova中檢定的是各組的平均數是否相同，檢定統計量F值為27.526，而顯著性p值(p-value)為7.898e-08，有e表示是科學記號格式，意思是7.898乘以-10的8次方，也就是非常非常小的小數，低於顯著水準0.05，表示三組資料之間有顯著差異。但到底這三組中誰大誰小，則還要用事後多重比較來分析。

Games-Howell事後多重比較 / Games-Howell Post-Hoc Test

### Games-Howell Post-Hoc Test



### Oneway Anova for y=value and x=group (groups: red, green, blue)



Omega squared: 95% CI = [.12; .44], point estimate = .28

Eta Squared: 95% CI = [.14; .42], point estimate = .3


                                     SS Df    MS     F     p

Between groups (error + effect)  857.2  2 428.6 14.81 <.001

Within groups (error only)      1996.4 69 28.93            





### Post hoc test: games-howell


            diff ci.lo ci.hi    t    df     p

green-red  6.04  2.41  9.67 4.14 25.25  .001

blue-red   8.14  5.01 11.27 6.46 26.07 <.001

blue-green 2.10 -2.44  6.65 1.12 44.94  .507

Eta Squared的點估計值為0.3，表示高度關聯強度。然後我們看「### Post hoc test: games-howell」以下的部分，可以看到「red」、「green」、「blue」三組的比較結果。

在「green-red」的比較中，p值為0.001，小於顯著水準0.05，而差異值diff為6.04，表示前組「green」顯著大於後組「red」。同樣的，「blue-red」比較結果也顯示前組「blue」顯著大於後組「red」。但是「blue-green」比較結果沒有顯著差異。

這樣子就完成了Welch's anova跟Games-Howell事後多重比較囉。

結語 / Conclusion

本篇以R腳本實作了ANOVA跟Scheffe事後多重比較以及Welch's anova跟Games-Howell事後多重比較這兩套有母數統計方法，可以用來分析常態分佈的多組資料之間是否有顯著差異。ANOVA分析方法是相當常用的統計檢定技術，這套R腳本應該可以大幅簡化分析時的操作。

本篇的分析全部都可以在SPSS裡面實作，交互比較來看，R腳本的結果與SPSS計算結果完全相同。但是SPSS的群組必須用數字，在資料整理上比較麻煩，而各種操作也挺複雜的。現在寫成R腳本之後，未來分析就可以簡化許多囉。

這篇多組常態分佈資料的差異檢定的整理就到這裡為止。我本身並非統計專業，只是因為分析資料時會用到，所以看文獻來將它整理一下。如果有什麼統計上的建議，歡迎在下面留言提出指教。如果你覺得我整理的不錯的話，請幫我在AddThis分享工具按讚、將這篇分享到Facebook等社群媒體吧！感謝你的耐心閱讀，讓我們下一篇見。