變異數分析 Analysis of variance
陳恩毅 陳柏勳
前言
在閱讀本節之前，我們期望讀者已習得一些先備知識(常用的分佈，假設檢定)，或是有基本獨立思考的能力。下面的介紹以實際操作為重點，並附上簡潔的說明，模型和table請參考藍色的數字連結。連結僅提供可能的幫助，任何連結內的評論不代表筆者立場。
變異數分析基本假設
介紹 1
變異數分析用以比較兩組以上的樣本平均數差異，如果使用T檢定的方式兩兩比較：
1.當\(k\)大時，要進行\(k \choose 2\)次檢定，過程相當冗長繁瑣。
2.同時進行\(k \choose 2\)次檢定會增加型一錯誤率(Type I error)。
\(\alpha_{total} = 1-(1-\alpha)^{k \choose 2} > \alpha\) ， \(\alpha_{total}\) 是整體型一誤差。
\(\bigstar\) 一律預設\(\alpha = 0.05\) 。
用變異數分析可以避免以上的情形。
若變異數分析的結果為各組樣本平均無明顯差異，則不需再作任何檢定。反之，需要作事後檢定，以便檢測哪些組平均不相同。
變異數分析之統計分析假設通常會依照各種模式型態不同而有差異，但廣義而言，變異數分析一共有三大前提假設：
1.常態性 : 各組樣本的分佈必須接近常態分佈。
2.獨立性 : 所有樣本必須互相獨立。
3.同質變異數 : 各組樣本的變異數必須相等或接近。
檢驗樣本是否為常態分佈
\(H_0\) : 樣本為常態分佈
\(H_1\) : 樣本非常態分佈
方法有很多種，我們舉了以下2個例子作為參考。
法一 : 畫Q-Q plot 2
Q-Q plot可以用來檢驗該樣本資料是否符合某種分佈(常態分佈、卡方分佈等等)，在R commander 中可以選擇。
eg 1-1 : Q-Q plot
從指數分佈抽取樣本。
機率分佈 => 連續型分佈 => 指數分佈 => 指數分佈隨機抽樣
\(\bigstar\) 取消勾選“樣本平均數”。
繪圖 => 分量比較(QQ)圖
這裡我們只選取第一組觀察值(20個樣本)。
依照原本選項的設定。
可以選擇是否要自動辨識樣本點(離群值)。
繪圖結果如下 :
圖上的sample 9 和 sample 17 看起來較為偏離常態。
由於Q-Q plot是以主觀性來判斷樣本是否偏離特定的分佈，以下介紹一種檢定方式以便我們作客觀的決策。
法二 : 用Shapiro-Wilk常態性檢定 3
eg 1-2 : Shapiro-Wilk常態性檢定
同上題樣本。
統計量 => 摘要 => 常態分佈檢定
由p值來看，p-value = 0.0003759 < 0.05，因此我們拒絕\(H_0\)，樣本不服從常態分佈。
\(\clubsuit\) 文章 : 如何檢查資料是否接近常態分佈 4
檢驗多組樣本變異數是否有差異
\(H_0\) : \(\sigma_1 = ... = \sigma_k\)
\(H_1\) : \(\sigma_i\) \(\neq\) \(\sigma_j\) , for some \(i,j \in\) \(\{1,...,k\}\)
我們舉例以下2種檢驗方法 :
1.Bartlett’s test 5
2.Levene’s test 6
wiki中的一段話:
Bartlett’s test is sensitive to departures from normality. That is, if the samples come from non-normal distributions, then Bartlett’s test may simply be testing for non-normality. Levene’s test and the Brown–Forsythe test are alternatives to the Bartlett test that are less sensitive to departures from normality.
簡單來說，Bartlett’s test較適用於接近常態分佈的樣本，若樣本偏離常態分佈，則Bartlett’s test的檢定結果可信度會下降，相對而言會比較建議使用Levene’s test或Brown–Forsythe test。
\(\clubsuit\) 實作範例 : 比較二或多組變異數Bartlett’s檢定 7
以下2個範例(常態分佈，指數分佈)分別使用Bartlett’s test和Levene’s test。
eg 2-1 : 從常態分佈抽取樣本
機率分佈 => 連續型分佈 => 常態分佈 => 常態分佈隨機抽樣
\(\bigstar\) 取消勾選“樣本平均數”。
資料 => 使用中的資料集 => 堆疊使用中資料集內的變數
在後面的完全隨機化設計，會使用到這裡的StackedData1,2。
Bartlett’s test
統計量 => 變異數檢定 => Bartlett 檢定
由p值來看，p-value = 0.2936 > 0.05，因此我們無法拒絕\(H_0\)，即各組變異數相同。
eg 2-2 : 從指數分佈抽取樣本
機率分佈 => 連續型分佈 => 指數分佈 => 指數分佈隨機抽樣
資料 => 使用中的資料 => 堆疊使用中資料集內的變數
Levene’s test
\(\bigstar\) 以下的Levene’s test 我們都設定\(\bar Y_i.\)為平均數。
統計量 => 變異數檢定 => Levene 檢定
由p值來看，p-value = 0.812 > 0.05，因此我們無法拒絕\(H_0\)，即各組變異數相同。
以上2個範例中的樣本皆符合變異數相同的假設，且假設樣本獨立，接下來我們會根據樣本分佈是否為常態來決定應該用何種檢定方式。
單因子變異數分析
完全隨機化設計 : Completely Randomized Design
介紹 8
檢驗多組樣本平均是否有差異(處理(treatment)效應是否存在)
\(H_0\) : \(\mu_1 = ... = \mu_k\)
\(H_1\) : \(\mu_i\) \(\neq\) \(\mu_j\) , for some \(i,j \in\) \(\{1,...,k\}\)
eg 3-1 : One way Anova CRD
樣本服從常態分佈且變異數相同
\(\bigstar\) 以下的StackedData1,2 分別為eg 2-1，2-2所用的資料。
A.使用 StackedData1 (常態分佈)
在進行變異數分析前，我們先將資料進行模式配適。
統計量 => 模型配適 => 線性模型
模型 => 假設檢定 => 變異數分析
下面的ANOVA列表使用預設即可。
由p值來看，p-value = 0.5042 > 0.05，因此我們無法拒絕\(H_0\)，即各組樣本平均相同。
eg 3-2 : Kruskal-Wallis rank sum test
介紹 9
檢定的是母體中位數，為單因子變異數分析CRD的無母數版本，不需同質變異數或常態分佈等假設。
把所有的n筆資料混和排序，改為等級。
\(\clubsuit\) 實作範例 10
B.使用 StackedData2 (指數分佈)
統計量 => 無母數檢定 => Kruskal-Wallis檢定
由p值來看，p-value = 0.7048 > 0.05，因此我們無法拒絕\(H_0\)，即各組樣本平均相同。
eg 3-3 例題實作(匯入csv檔資料)
\(\bigstar\) 如果輸入EXCEL檔出問題，可以轉成CSV檔。
資料 => 匯入資料 => 匯入文字檔,剪貼簿或URL檔
給定隨機採收的台灣蘋果，美國蘋果以及日本蘋果各20個樣本的重量，試檢驗各組之間的平均重量是否有差異。
資料如下圖 :
先將資料堆疊。
常態性檢定
\(H_0\) : 樣本為常態分佈
\(H_1\) : 樣本非常態分佈
法一 : 繪製Q-Q plot
結果如下 :
顯示大部分樣本都接近常態，我們再用Shapiro-Wilk常態性檢定。
法二 : 用Shapiro-Wilk常態性檢定
統計量 => 摘要 => 常態分佈檢定
結果如下 :
由p值來看，p-value = 0.4908 > 0.05，因此我們無法拒絕\(H_0\)。
Q-Q plot的結果沒有顯著偏離常態的樣本，且Shapiro-Wilk常態性檢定的結果為不拒絕常態的假設，p值為0.4908大於0.05，我們較傾向於樣本服從常態分佈。
同質變異數檢定
上面檢驗樣本大致服從常態分佈，因此我們應使用 Bartlett 檢定來檢驗各組樣本變異數是否相同。
統計量 => 變異數檢定 => Bartlett 檢定
由p值來看，p-value = 0.6976 > 0.05，因此我們無法拒絕\(H_0\)。
接下來用單因子變異數分析來檢驗這3組樣本平均是否有差異。
法一 :
統計量 => 平均數檢定 => 單因子變異數分析(One-way ANOVA)
由p值來看，p-value = 0.00561 < 0.05，因此我們拒絕\(H_0\)，即各組平均重量不相同。
法二 :
統計量 => 模型配適 => 線性模型
模型 => 假設檢定 => 變異數分析
由p值來看，p-value = 0.005612 < 0.05，因此我們拒絕\(H_0\)，即各組平均重量不相同。
事後檢定 (繪圖 )
檢定的結果顯示各組平均重量不相同，因此需要進行事後檢定，但R commander 沒有此功能，所以我們用畫圖的方式來比較各組平均。
繪圖 => 平均數圖
資料介面不須更改。
選項介面使用信賴區間，信賴水準為0.95
繪圖如下 :
由圖來看，apple3平均較高一些。
但如果只看圖不一定能正確判別，因此可以在R語法檔輸入以下的指令，執行事後檢定，以下使用Scheffe 檢定。
install.packages(“agricolae”)
library(agricolae)
scheffe.test(AnovaModel.1,trt=“factor”, group=FALSE,console=TRUE)
參數說明 :
AnovaModel.1 為Anova模型
trt 為處理
group 為分群的邏輯值
console 為print的邏輯值
由p值來看，apple1 - apple3的p-value = 0.0080 < 0.05，因此我們可以得知 Apple1 和 Apple3 的平均有顯著差異。
若group = TRUE，將不會顯示上圖的檢定區域，而是改成下圖的部分 :
顯示apple3為a群，apple1為b群。
Post Hoc Scheffe Test Usage 11
隨機區集化設計
隨機區集化設計 : Randomized Block Design
介紹 12
檢驗處理效應與區集效應是否存在
我們把可控的干擾因子透過區集化的實驗設計方法控制住，干擾因子的每個水準可視為一個區集，那麼在同一個區集內的干擾因子將會有同質性，在處理(treatment)之間比較時，可以移除此一干擾因子產生的效應。
eg 4-1 : RBD
設公司有三部機器，分別由五位作業員操作，記錄各作業員操作各部機器的產量，結果如下 :
表4.1
欲檢驗作業員的能力(處理)和機器產量(區集)是否有差異。
因此有2個假設 :
檢定處理效應
\(H_0\) : 處理效應不存在
\(H_1\) : 處理效應存在
檢定區集效應
\(H_0\) : 區集效應不存在
\(H_1\) : 區集效應存在
將表4.1轉為長表格，如下圖 :
統計量 => 模型配適 => 線性模型
點選反應變數~控制變數+控制變數
下面是英文版本的R Commander，稍微不同的地方是它有Blocked ANOVA這項功能。
影片來源 13
法一 : 因為我們的版本沒有Blocked ANOVA 這個選項，所以在R語法檔打上anova(線性模型)，然後點選執行語法。
由p值來看，p-value = 0.0121931 < 0.05，因此我們拒絕\(H_0\)，處理效應存在。
由p值來看，p-value = 0.0004389 < 0.05，因此我們拒絕\(H_0\)，區集效應存在。
法二 : 模型 => 假設檢定 => 變異數分析(ANOVA表)
使用預設即可。
由p值來看，p-value = 0.0121931 < 0.05，因此我們拒絕\(H_0\)，處理效應存在。
由p值來看，p-value = 0.0004389 < 0.05，因此我們拒絕\(H_0\)，區集效應存在。
eg 4-2 : Friedman test
介紹 14
檢定的是母體中位數，為RBD的無母數版本，不需同質變異數或常態分佈等假設。
在每一個區集內部把資料排序，改為等級。
下表為5位分析師對4種創業投資事業的偏好程度，欲檢驗分析師對4種創業投資事業的偏好程度是否有差異。
\(H_0\) : 分析師對4種創業投資事業的偏好程度無差異
\(H_1\) : 分析師對4種創業投資事業的偏好程度有差異
統計量 => 無母數檢定 => Friedman等級和檢定
選取因子。
由p值來看，p-value = 0.0211 < 0.05，因此我們拒絕\(H_0\)，亦即分析師對4種創業投資事業的偏好程度有顯著差異。
雙因子變異數分析
介紹 15
\(\clubsuit\) 文章 : SPSS操作16
檢驗各因子效應與交互作用效應是否影響反應變數(是否顯著)
考慮兩種因子的實驗設計時，若A因子有a個水準，B因子有b個水準，它們之間所有的組合(a\(\times\)b種)，都要進行實驗與研究，故會探討a\(\times\)b個處理。
eg 5 : Two way Anova
欲檢驗A因子(工廠)，B因子(天氣)，以及其交互作用效應是否顯著。
資料如下表5.1 :
已知樣本獨立且服從常態分佈，並具有同質變異數。
有3個假設 :
檢定A因子效應
\(H_0\) : A因子不影響反應變數
\(H_1\) : A因子影響反應變數
檢定B因子效應
\(H_0\) : B因子不影響反應變數
\(H_1\) : B因子影響反應變數
檢定交互作用效應
\(H_0\) : 交互作用不存在
\(H_1\) : 交互作用存在
將表5.1轉為長表格。
統計量 => 平均數檢定 => 多因子變異數分析
選取因子。
由p值來看，p-value = 0.07872 > 0.05，因此我們無法拒絕\(H_0\)，A因子不影響反應變數。
由p值來看，p-value = 0.09243 > 0.05，因此我們無法拒絕\(H_0\)，B因子不影響反應變數。
由p值來看，p-value = 0.33182 > 0.05，因此我們無法拒絕\(H_0\)，交互作用不存在。
若想看交互作用的繪圖，可以用以下的程式碼:
interaction.plot(Dataset$A,Dataset$B,Dataset$outcome)
繪圖結果如下:
若想從另一個因子的方向看，將A，B位置調換即可。
