13 比較多組平均值(單因子變異數分析)
譯者註 20240414完成實例演練以外的編修。
這個單元要介紹心理學研究最常使用的一種統計方法~“變異數分析”,通常簡稱ANOVA。羅蘭.費雪爵士在20世紀初奠定了今天這套方法的基本運算原則,同時他給的命名也困擾著今天的學習者,變異數分析這個名稱通常會造成兩種誤會。其一是「變異數」,實際上ANOVA是比較平均數之間的差異。其二是有好幾套統計方法都是奠基於變異數分析,然而有些方法與變異數的關係非常微弱。後面的單元裡,讀者會學到各式各樣的變異數分析方法,分別有各自適用的條件。這個單元要學的只是最簡單的單因子變異數分析,適用研究設計只有幾個實驗組,研究者想要分析每個實驗組在各獨變項條件之間的測量結果差異。
這個單元的學習順序是:首先介紹這個單元用來解說及示範jamovi操作的虛擬資料集。接著說明單因子變異數分析的運算原理,然後說明如何使用 jamovi的變異數分析模組執行變異數分析程序。這兩小節是這個單元的重點。
接下來分別討論在執行變異數分析時必須考慮的一系列重要課題,像是如何計算效果量大小、事後檢定和多重比較的校正,以及變異數分析的適用條件。我們還會討論如何檢查這些條件,以及適用條件不成立時有什麼樣的補救措施。最後一節,我們會學習[單因子重覆量數變數分析]。
13.1 獨立樣本變異數分析示範資料
想像你正在協助進行一項臨床試驗,測試一種名為Joyzepam的新型抗憂鬱藥物的藥效。為了能公平地測試這種新藥的效果,需要分別測試包括新藥的三種藥物,另外兩種藥物之一是安慰劑,還有已經上市的抗憂鬱/抗焦慮藥物,名為Anxifree。研究一開始招募18位患有中度至重度抑鬱症的參與者。其中有一半參與者不只是服藥,同時進行認知行為治療(CBT),另一半參與者未同時進行任何心理治療。藥物以雙盲隨機方法分派給參與者,因此每種藥物分派給3位有進行CBT的參與者及3位未進行心理治療的參與者。每位參與者各自使用藥物3個月後,研究者再評估參與者的情緒改善狀況,以\(-5\)到\(+5\)的數值代表每位參與者的情緒改善狀況。讀者可以載入資料集 Clinical Trial,瀏覽範例資料的內容,其中的變項分別是藥物、治療和情緒改善分數。
為了學習如何使用單因子變異數分析,這裡的目的是要評估各種藥物改善情緒狀況的效果。首先要進行描述統計及繪製統計圖表,我們從 單元 4 已經學到如何使用jamovi完成描述統計,報表會如同 圖 13.1 。
從 圖 13.1 可以看出,服用Joyzepam的參與者,情緒的改善程度優於服用Anxifree及安慰劑。Anxifree的情緒提升程度優於安慰劑,但是沒有像Joyzepam那麼明顯。這裡要回答的問題是,這些藥物的效果是否“真正有效”,還是只是一次偶然的發現?
13.2 單因子變異數分析的運算原理
為了運用臨床試驗資料回答以上的問題,我們要學習使用單因素變異數分析(one-way ANOVA)。如果讀者不知如何操作jamovi的ANOVA模組選單裡眼花撩亂的選項,先仔細閱讀這一節說明的基本原理,了解ANOVA程序每個步驟的運算概念,跟著實例演練操作一兩次,掌握概念後,後續與ANOVA有關的統計方法就不必如此學習了。
獨立樣本變異數分析示範資料的說明提到,研究人員有興趣的是三種藥物改善參與者憂鬱情緒的效果,這個研究設計的分析問題類似 單元 11 介紹的t檢定範例,不過要比較的不只兩組。在此先定義\(\mu_P\)代表安慰劑的情緒變化母群平均值,\(\mu_A\)和\(\mu_J\)分別對對應Anxifree和Joyzepam兩種藥物效果的平均值,所以要檢定的虛無假設就是:三組的母群平均值是相等的。也就是有點悲觀的預測,兩種藥物的效果都沒有比安慰劑好。這樣的虛無假設可以寫成:
\[H_0: \text{結果顯示 } \mu_P=\mu_A=\mu_J\]
因此對立假設就是:三種藥物之中,至少有一種的效果不同於其他兩種。用數學式表示的話,可能會讓有些同學困惑,由於有很多方式能表達虛無假設是錯誤的,我們先將對立假設寫成:
\[H_1: \text{結果}\underline{沒有}\text{顯示 } \mu_P=\mu_A=\mu_J\]
這道虛無假設比之前單元遇到的要棘手得多,應該要如何進行檢定呢?因為這個單元的標題是變異數分析,聰明的讀者應該猜到就是用「變異數分析」,但是如初學的讀者至此可能還不太了解為何這個方法其實是用來處理平均值。這是許多學生第一次上到變異數分析時,最大的挑戰。為了說明運算原理,我們要從變異數的組成談起,請先參考 圖 13.2 了解什麼是組間變異(Between-group variation)及組內變異(Within-group variation)。
13.2.1 計算依變項變異數的兩套公式
我們先定義幾個運算用的符號:G代表分組的數量,因為資料集有三種藥物,所以有 \(G = 3\) 個分組。然後定義\(N\)表示總樣本量,這個資料集有 \(N = 18\) 位參與者,任可一組的人數同樣用 \(N_k\) 表示。這份資料集的三組樣本量都是\(N_k = 6\)。1 最後是定義代表結果變項的Y,也就是每位參與者的情緒狀況改善數值,這裡用\(Y_{ik}\) 代表第 k 組的第 i 位參與者改善數值。因此, \(\bar{Y}\)是所有18位參與者的平均改善數值,\(\bar{Y}_k\)就是第 k 組的第6位參與者的改善狀況。
至此會用的符號都已經就定位,可以開始寫公式了。先回想一下談描述統計的 小單元 4.2 提過的變異數公式,這裡的結果變項Y的樣本變異數公式是 \[Var(Y)=\frac{1}{N}\sum_{k=1}^{G}\sum_{i=1}^{N_k}(Y_{ik}-\bar{Y})^2\] 這個公式看起來和 小單元 4.2 提到的變異數公式長得幾乎一樣。唯一的區別是這個公式有兩個連加記號:各組\(k\)的平均值總和及組內所有參與者 \(i\) 個人數值的總和。請留意一下符號表示的差別,若結果變項符號是 \(Y_p\),代表資料裡第p位參與者的數值,這樣子只會有所有參與者個人數值的總和。這裡要寫兩個連加符號,原因是要先將每筆數值歸到其中一組,再指定各組內所代表的個人數值。
這裡用具體的例子來理解應該會有用。來看看 表 13.1 的例子,一供有 \(N = 5\) 個人分為\(G = 2\) 組。我們可以武斷地指定「酷」的人是第 1 組,「不酷」的人是第 2 組。最後列出其中有三個人很酷(\(N_1 = 3\))和兩個人不算酷(\(N_2 = 2\))。
name | person P | group | group num. k | index in group | grumpiness \( Y_{ik} \) or \( Y_p \) |
---|---|---|---|---|---|
Ann | 1 | cool | 1 | 1 | 20 |
Ben | 2 | cool | 1 | 2 | 55 |
Cat | 3 | cool | 1 | 3 | 21 |
Tim | 4 | uncool | 2 | 1 | 91 |
Egg | 5 | uncool | 2 | 2 | 22 |
這個表格結合兩種標記方式。變項 p代表個人,所以用\(Y_p\)代表第 p 人的沮喪指數。像是第四位是Tim,就用 \(p = 4\) 代表他。我們若要用數字討論「Tim」這個人的沮喪程度,可以用\(Y_4 = 91\)來溝通。而這不是唯一可以描述Tim的方式,另一種方式是根據Tim的分組。因為Tim 是「不酷」組(\(k = 2\))的第一人(\(i = 1\)),也可以用\(Y_{12} = 91\)代表Tim的沮喪程度。
也就是說,每個人 p 都對應一個獨一無二的 ik 組合,這也就能解釋為他我說上面的變異數公式,與以下更早學到的變異數公式相同的 \[Var(Y)=\frac{1}{N}\sum_{p=1}^{N}(Y_p-\bar{Y})^2\]
這兩個公式都是求樣本資料裡所有觀察值的總和,因為\(Y_p\)的公式較簡單,做運算練習的功課大都是用第二個公式。但是變異數分析必須要區別那位參與組是屬於那一組,就需要用\(Y_{ij}\)的公式來做運算。
13.2.2 變異數與離均差平方和
好啦,更進一步認識變異數的計算方法後,就能討論什麼是總離均差平方和(total sum of squares),簡記為\(SS_{tot}\)。變異數是離均差平方和的平均結果,計算\(SS_{tot}\)只要算總和就好。2
在ANOVA的單元討論如何分析變異數,其實是談如何處理離均差平方和,而不是講如何處理變異數。3
接著來談較難講的組間離均差是怎麼回事,這裡就要來拆解分組平均值 \(\bar{Y}_k\) 和總平均值 \(\bar{Y}\)是怎麼構成了。4
如此看來,所有實驗參與者之間的總離均差平方和(\(SS_{tot}\)),等於組間離均差平方和(\(SS_b\)),加上組內離均差平方和(\(SS_w\))。其實不怎麼需要證明
\[SS_w+SS_b=SS_{tot}\]
太棒了!
所以我們學到什麼?我們已經懂得與結果變項相關的總離均差平方和(\(SS_{tot}\)),可以被劃分為“分組的樣本平均值之間的差異所造成的變異”的總和(\(SS_b\)),加上“除此之外的任何變異”的總和(\(SS_w\))5。
那要怎麼幫研究人員確認各組的母群平均值相不相等呢?嗯,稍等一下,有沒有看到每種離均差平方和都是各種平均值之間的差異?若是虛無假設的推測是成立的,代表各組樣本平均值\(\bar{Y_k}\) 應該是非常相等的,沒錯吧?也就是說\(SS_b\)要非常小,至少要比“除此之外的任何變異”(\(SS_w\))明顯的小。有沒有聞到假設檢定的味道了?
13.2.3 離均差平方和與F檢定
前一節我們已經了解ANOVA運算的核心思想是比較\(SS_b\) 和 \(SS_w\),若是組間離均差平方和\(SS_b\)是相對大於組內離均差平方和\(SS_w\),就有理由懷疑各組的母群平均值並不相等。為了轉化為可操作的假設檢定程序,我們要進行一些“微調”。我們先來認識如何計算檢驗統計值——F值(F ratio),然後再來討論要計算F值的原因。
為了將離均差平方和轉換為F值,首先要計算\(SS_b\) 和 \(SS_w\)的自由度。通常自由度是指為了計算估計值,需要貢獻的“資料點”數量,減去滿足“要估計”的參數數量。組內離均差平方和的計算元素是\(N\)筆個別觀察值與\(G\) 個分組平均值之間的變異;而組間離均差平方和的計算元素是\(G\) 個分組平均值與唯一一個總體平均值之間的變異。所以有兩個自由度:
\[df_b=G-1\] \[df_w=N-G\]
這看起來很簡單,是吧。接著是將平方和轉換為“平方和的平均”,也就是各自除以自由度:
\[MS_b=\frac{SS_b}{df_b}\] \[MS_w=\frac{SS_w}{df_w}\]
最後,將組間 MS 除以組內 MS 就能算出 F 值:
\[F=\frac{MS_b}{MS_w}\]
表面上來看,F統計值的運算相當直覺又好懂。F 值越大,就表示組間變異與組內變異之間的差異越大。所以F值起大,能反駁虛無假設的證據力就越強。但是 \(F\) 要有多大才能確實拒絕 \(H_0\)?為了清楚理解,我們要更深入了解 ANOVA 是什麼,以及什麼是離均差平方和的平均?
這些問題在實例演練將有詳細討論,不過如果讀者對於詳細的回答不感興趣,這裡提供一個簡短的說明。為了完成假設檢定,分析人員需要知道虛無假設符合推測時的F值取樣分佈。這一點都不奇怪,因為由虛無假設生成的F統計值取樣分佈,必定是一個F分佈。回想 單元 7 有關F 分佈的說明,\(F\) 分佈的兩個參數各對應兩個自由度。第一個 \(df_1\) 對應組間自由度 \(df_b\),第二個 \(df_2\) 對應組內自由度 \(df_w\)。
between
groups | within
groups | |
---|---|---|
df | \( df_b=G-1 \) | \( df_w=N-G \) |
sum of squares | \( SS_b=\sum_{k=1}^{G} N_k (\bar{Y}_k-\bar{Y})^2 \) | \( SS_w=\sum_{k=1}^{G} \sum_{i=1}^{N_k} (Y_{ik}-\bar{Y}_k)^2 \) |
mean squares | \( MS_b=\frac{SS_b}{df_b} \) | \( MS_w=\frac{SS_w}{df_w} \) |
F-statistic | \( F=\frac{MS_b}{MS_w} \) | - |
p-value | [complicated] | - |
表 13.2 摘要所有單因子變異數計算過程所產生的關鍵計數,包括各項計數的計算公式。
[額外的技術細節 6]
13.2.4 實例演練
前面的步驟說明包括很多抽象符號,而且只讀字面不能充分理解,在此我們先做個小結。回到一開始的獨立樣本變異數分析示範資料,根據這份資料的描述統計表,我們知道各種藥劑的情緒改善平均分數:安慰劑是 \(0.45\),Anxifree 是 \(0.72\),Joyzepam 是 \(1.48\)。有了這些資訊和計算公式,回到1899年就能用鉛筆和紙辦個一人玩的統計桌遊7。不過今天不是1899年,人人手上都有計算機,就讓我們偷個懶,只要做前5筆資料計算就好。我們先從計算組內離均差平方和\(SS_w\) 開始,要畫張像 表 13.3 的表格來輔助計算。
group k | outcome \( Y_{ik} \) |
---|---|
placebo | 0.5 |
placebo | 0.3 |
placebo | 0.1 |
anxifree | 0.6 |
anxifree | 0.4 |
第一步只有列出原始資料,也就是每位參與者服用的藥物(分組變項),以及情緒改善狀況(結果變項)。請注意,結果變項列出的數值,對應前面提到的\(\bar{Y}_{ik}\)。接下來就是計算每位參與者所屬分組平均值,\(\bar{Y}_k\)。這個步驟並不困難,因為描述統計表裡都有了,加上後就變成 表 13.4 。
group k | outcome \( Y_{ik} \) | group mean \( \bar{Y}_k \) |
---|---|---|
placebo | 0.5 | 0.45 |
placebo | 0.3 | 0.45 |
placebo | 0.1 | 0.45 |
anxifree | 0.6 | 0.72 |
anxifree | 0.4 | 0.72 |
算到這一步,接著就能計算每位參與者的離均差,也就是\(Y_{ik} - \bar{Y}_k\),順便計算每個離均差的平方值\((Y_{ik} - \bar{Y}_k)^2\),就得到 表 13.5
group k | outcome \( Y_{ik} \) | group mean \( \bar{Y}_k \) | dev. from group mean \( Y_{ik} - \bar{Y}_k \) | squared deviation \( (Y_{ik}-\bar{Y}_k)^2 \) |
---|---|---|---|---|
placebo | 0.5 | 0.45 | 0.05 | 0.0025 |
placebo | 0.3 | 0.45 | -0.15 | 0.0225 |
placebo | 0.1 | 0.45 | -0.35 | 0.1225 |
anxifree | 0.6 | 0.72 | -0.12 | 0.0136 |
anxifree | 0.4 | 0.72 | -0.32 | 0.1003 |
計算組內離均差平方和就簡單多了,將所有離均差平方加起來就好:
\[ \begin{split} SS_w & = 0.0025 + 0.0225 + 0.1225 + 0.0136 + 0.1003 \\ & = 0.2614 \end{split} \]
當然,真正完整的計算是要將全部18筆觀察值進行以上四個計算步驟,而不是只算五筆資料。有心的讀者可以繼續用筆計算,只是這麼做會相當繁瑣。其實讀者也可以用LibreOffice或Excel試算表軟體演練,計算難度會大幅減低。讀者可以用自已的Excel開一個檔案,命名為clinicaltrial_anova.xls。完成這裡示範的四個步驟,就能得到組內離均差平方和的值是\(1.39\)。
到這裡已經算出組內離均差平方和\(SS_w\)的值,接下來就是處理組間離均差平方和\(SS_b\) 。計算步驟相當類似,不同的地方是要改成計算每個分組的平均值 \(\bar{Y}_k\) 與總平均值 \(\bar{Y}\)的差異,計算表格就變成 表 13.6 ,得到 \(0.88\)。
group k | group mean \( \bar{Y}_k \) | grand mean \( \bar{Y} \) | deviation \( \bar{Y}_k - \bar{Y} \) | squared deviation \( ( \bar{Y}_k-\bar{Y})^2 \) |
---|---|---|---|---|
placebo | 0.45 | 0.88 | -0.43 | 0.19 |
anxifree | 0.72 | 0.88 | -0.16 | 0.03 |
joyzepam | 1.48 | 0.88 | 0.60 | 0.36 |
不過因為是計算組間的平方和,還要將每個離均差平方先乘以 \(N_k\),也就是各組的觀察值數量。理由是因為同組都有\(N_k\)個觀察值,計算組間差異需要加權。安慰劑組有六筆資料,而安慰劑組的平均值與總平均值相差 \(0.19\),因此六位參與者的組間離均差平方加權總和是\(6 \times 0.19 = 1.14\)。所以計算表格會變成 表 13.7 。
group k | ... | squared deviations \( (\bar{Y}_k-\bar{Y})^2 \) | sample size \( N_k \) | weighted squared dev \( N_k (\bar{Y}_k-\bar{Y})^2 \) |
---|---|---|---|---|
placebo | ... | 0.19 | 6 | 1.14 |
anxifree | ... | 0.03 | 6 | 0.18 |
joyzepam | ... | 0.36 | 6 | 2.16 |
接著就可以將各分組的“加權的離均差平方”加起來,得到組間離均差平方和:
\[\begin{aligned} SS_b & = 1.14 + 0.18 + 2.16 \\ &= 3.48 \end{aligned}\]
讀者肯定會發現,計算分組離均差平方和的步驟比較少 8。即然己經算好\(SS_b\) 和 \(SS_w\)的值,剩下的工作就相當簡單了。接下來是計算自由度。由於知道 \(G = 3\) 個分組以及 \(N = 18\) 個觀察值,自由度用簡單的減法就可以計算:
\[ \begin{split} df_b & = G-1 = 2 \\ df_w & = N-G = 15 \end{split} \]
有了離均差平方和和自由度,用前者除以後者,就能算出組內與組間變異數,也就是離均差平方和的平均:
\[ \begin{split} MS_b & = \frac{SS_b}{df_b} = \frac{3.48}{2} = 1.74 \\ MS_w & = \frac{SS_w}{df_w} = \frac{1.39}{15} = 0.09 \end{split} \]
快要完成ANOVA的計算了。有了兩個平均值,就能計算最想知道的F值。就是將\(MS_b\)除以\(MS_w\)。
\[ \begin{split} F & = \frac{MS_b}{MS_w} = \frac{1.74}{0.09} \\ & = 19.3 \end{split} \]
我們終於算完了,值得慶祝!有了統計值,最後就是檢定這個統計值是否代表結果是顯著的。如同在 單元 9 討論過的手作方法,有心的讀者可以翻開任何一本統計教科書最後的附錄,查找其中有關F檢定的表,只要找到指定 \(\alpha\)值所對應的閾值,例如 \(0.05\),\(0.01\) 或 \(0.001\),搭配 自由度2 和15。設定\(\alpha\)為 \(0.001\) 的話,會查到F的閾值是\(11.34\)。因為小於最後算出的F值,可以宣稱\(p < 0.001\)。但這是時代的眼淚,今天各式各樣統計軟體都能為我們算出精確的p值。這個範例的p值是 \(0.000071\)。除非我們對型一錯誤率採取非常謹慎的立場,我們幾乎 可以結論這樣的研究結果能拒絕虛無假設。
至此己經完成了ANOVA的基本計算,將以上步驟算出的數值,整理成 表 13.1 的報表,這是傳統的變異數分析報告規範。完整的ANOVA報表請見 表 13.8 。
df | sum of squares | mean squares | F-statistic | p-value | |
---|---|---|---|---|---|
between groups | 2 | 3.48 | 1.74 | 19.3 | 0.000071 |
within groups | 15 | 1.39 | 0.09 | - | - |
演練到這裡,讀者應該不大想純靠手工整理這樣的報表,幾乎所有專業統計軟體,包括 jamovi,都能將ANOVA 的結果統整成像 表 13.8 的表格,所以大家最好習慣一下。儘管軟體能輸出完整的 ANOVA 報表,真正的報告並不需要置入整份表格,現代學術寫作規範所建議的報告格式是像以下的句型:
單因子變異數分析顯示,測試的藥物對情緒改善有顯著影響,F(2,15) = 19.3,p < .001。
天哪!我們做得要死要活,只為了寫出這短短的句子。
13.3 jamovi的變異數分析模組
若是讀者認真地跟著前一節的示範步驟,在計算紙上演算,肯定能了解變異數分析程序有多麼複雜,而且一不小心就會出錯。接下來還有更多計算,如果還是用紙筆演練,各位會覺得幹麼要浪費時間做這些計算?其他要學的事都會被ANOVA的計算耽擱了!
13.3.1 使用jamovi進行變異數分析
為了讓大家從複雜的變異數分析計算中解脫,就要靠jamovi了!請開啟ANOVA
模組的ANOVA設定選單,將mood.gain
搬到”Dependent Variable”對話框,再將drug
搬到”Fixed Factors”對話框,就能在報表區看到如同 圖 13.3 的結果。9 請留意設定選單列出的選項,要勾選 \(\eta^2\) ,讓報表顯示效果量的計算結果,稍後會進一步說明這項效果量的意義。
jamovi 的ANOVA結果表包括離均差平方和、自由度以及目前還沒有真正談到的其他指標。不過,請各位留意,印在jamovi報表標籤並不是「組間」(Between)或「組內」(Within),反而是資料檔的變項名詞(drug
)和計算過程的產物(Residual
)。因為這個範例的組間離均差平方和,對應藥物對結果變項的影響,組內離均差平方和對應的是“殘餘”變異,所以習慣用“殘差”(Residual)稱呼。比較jamovi計算的數值與 小單元 13.2.4 的手算結果,會看到彼此或多或少是相同的。除了四捨五入的差異,組間平方和為 \(SS_b = 3.45\),組內平方和為 \(SS_w = 1.39\),自由度分別是 \(2\) 和 \(15\)。還有報告 F 值和 p 值,也是和手動計算的數字差不多,只是四捨五入的差異。
13.4 效果量
衡量 ANOVA結果的效果量有好幾種指標,最常用的是 \(\eta^2\)( eta sqaured )和淨 \(\eta^2\)( partial eta sqaured )。就單因子變異數分析來說,這兩種指標的解讀方式是相同的,因此這裡只解釋 \(\eta^2\)。\(\eta^2\) 定義相當簡單明瞭:
\[\eta^2=\frac{SS_b}{SS_{tot}}\]
前面有認真演練計算過程的同學一看就會懂。拿 圖 13.3 的ANOVA 報表裡的 \(\eta^2\) 驗算一下,會發現\(SS_b = 3.45\) 還有 \(SS_tot = 3.45 + 1.39 = 4.84\),所以可以算出:
\[\eta^2=\frac{3.45}{4.84}=0.71\]
\(\eta^2\)的解讀同樣直接了當,這項效果量公式就是直指結果變項的變異,可以被預測變項解釋的比例。\(\eta^2=0\) 代表兩種變項之間完全無關,\(\eta^2=1\)則表示兩者關係密切。值得一提的是,\(\eta^2\)與 小單元 12.6.1 提過的 \(R^2\) 是屬於同一種效果量家族的成員,有效力相等的解釋力。許多統計教科書建議使用 \(\eta^2\) 作為ANOVA的預設效果量指標,不過荷蘭心理學者Daniel Lakens在一篇有趣的blog文章提到,以統計實務處理的真實資料來看,\(\eta^2\)可能不是最好的效果量指標,因為這是一個有偏誤的估計值。還好jamovi可以選擇另一個估計偏誤較少的指標 \(\omega^2\) (omega squared) 。
13.5 多重比較與事後檢定
ANOVA報表顯示至少三組樣本之間存在顯著差異時,研究者會想知道是那幾組之間存在差異。以藥物臨床試驗的範例來說,由於虛無假設是三種藥物的情緒改善狀況完全一樣。不過仔細一想的話,這個範例的虛無假設包括三項不一樣的預測:
- Anxifree這款舊藥的效果沒有比安慰劑更好: \(\mu_A = \mu_P\)
- 新藥Joyzepam的效果沒有比安慰劑更好: \(\mu_J = \mu_P\)
- Anxifree 和 Joyzepam 兩款新藥的效果沒什麼差別: \(\mu_A = \mu_J\)
如果以上三項預測有一項不成立,整個虛無假設也是不成立。因為ANOVA的分析結果已經確認可以拒絕虛無假設,所以至少有一項預測是不成立的。儘管任何一項預測被拒絕都是一個有趣的結論,但是是那個預測被拒絕呢?即然關注的是新藥 Joyzepam的效果有沒有比安慰劑好,新藥與舊藥Anxifree的比較就非常重要。想清楚真正關注的目標,真正有意義的對比更有可能是Anxifree與安慰劑的效果差異。雖然已經有許多研究人員比較過Anxifree與安慰劑的效果,這次分析能再次確認早期的研究成果可以重現
將虛無假設分解為三項預測後,就能將八種可能的分析結果,表列如 表 13.9 。
possibility: | is \( \mu_P = \mu_A \)? | is \( \mu_P = \mu_J \)? | is \( \mu_A = \mu_J \)? | which hypothesis? |
---|---|---|---|---|
1 | \( \checkmark \) | \( \checkmark \) | \( \checkmark \) | null |
2 | \( \checkmark \) | \( \checkmark \) | alternative | |
3 | \( \checkmark \) | \( \checkmark \) | alternative | |
4 | \( \checkmark \) | alternative | ||
5 | \( \checkmark \) | \( \checkmark \) | \( \checkmark \) | alternative |
6 | \( \checkmark \) | alternative | ||
7 | \( \checkmark \) | alternative | ||
8 | alternative |
既然已經拒絕虛無假設,第一種分析結果已經出局。接著就要問,其餘七種可能的分析結果,那一個才是正確的結論?面對這樣的局面,要再看一眼描述統計報告,像是 圖 13.1 的統計圖,我們會發現 Joyzepam 比安慰劑和 Anxifree好,但是Anxifree 和安慰劑的效果沒有差別是最有可能的結論。不過,若是要提出精確的報告,就需要做進一步的檢定。
13.5.1 成對t檢定
要如何進行多重檢定?這裡要比較三對分組平均值:安慰劑 vs. Anxifree,安慰劑 vs. Joyzepam,和 Anxifree vs. Joyzepam,最簡單的方法是執行三次獨立樣本t檢定,逐一確認檢定結果。用 jamovi 很容易做到。打開 ANOVA 模組設定選單的 ‘Post Hoc Tests’(事後檢驗)選項,將 ‘drug’(藥物)變項移到右側的活動框中,然後勾選 ‘No correction’(無校正)核取方塊。就會產生一張整齊的表格,顯示藥物變項的三對比較 t 檢定結果,如同 圖 13.4 的展示。
13.5.2 多重檢定的校正
前面提到一次做多組t檢定可能會有問題,最大的問題是沒有明確的理論指引時,同時執行多組t檢定,就只是看有那幾組會出現顯著結果來做判斷。這樣一來,科學分析很容易就會變成一場競爭「誰是最後贏家」的遊戲。這種缺乏理論基礎的前提,直接比較多組差異的分析,稱為事後分析(“post hoc” 是拉丁語,意思是 “在真正要緊的事務之後”)。10
事後分析不是不可以做,而是要非常小心。例如成對t檢定介紹的方法應該儘量不要做,因為每組t檢定的型一錯誤率都設定為.05的話,執行了三組t檢定,會導致型一錯誤率過度膨脹。想像有筆資料要用ANOVA分析10組平均值,然後挑出其中45筆組間差異進行「事後」t檢定,找出那些組別之間有顯著差異。如同 單元 9 討論的例子,會那幾組出現顯著差異純粹是偶然的。虛無假設檢定有效的基本原則是,分析者要嚴格控制型一錯誤率,為了確定造成ANOVA顯著結果的原因,不管三七二十一執行多組t檢定,會導致全部45組差異的型一錯誤率膨脹到天邊了。
常用的解決方法是校正p值,控制總體型一錯誤率(Shaffer, 1995)。這種用於事後分析的校正方法,通常稱為多重比較校正,有時也叫「即時推論」(simultaneous inference)。其實統計學界已經開發出許多校正方法,這個單元還有 小單元 14.8 會介紹其中幾種,本書介紹的只是最常見的校正方法(有與趣的讀者請參考 Hsu, 1996 )。
13.5.3 Bonferroni校正
Bonferroni校正是其中最簡單的一種方法(Dunn, 1961)。如果現在要處理m個獨立的計檢定,這個方法能確保即使出現任何顯著結果,維持總體型一錯誤機率上限在\(\alpha\)。11 如同字面的定義,Bonferroni校正只是「將所有原始 p 值乘以 m」。做法是設定校正前的p值為\(p\),校正後的p值為\(p_j^{'}\),根據Bonferroni校正,校正前後的p值關係是:
\[p_j^{'}=m \times p\]
只要確認\(p_j^{'} < \alpha\),就能正確拒絕虛無假設。Bonferroni校正的邏輯非常簡單:現在要進行m組獨立的統計檢定,如果設定每個檢定的型一錯誤率上限是\(\frac{\alpha}{m}\),所有檢定的總型一錯誤率就不會大於 \(\alpha\)。這種方法簡單到Bonferroni在他的論文裡寫道:
因為這個方法如此簡單,又容易使用,我肯定以前一定有人用過。但是我翻遍任何文獻都找不到,只能推測:正是因為這種方法簡單到讓每個自認聰明絕頂的統計學家,不值得寫論文聲張這種方法的一點點好處(Dunn, 1961,第52-53頁)。
要用 jamovi 執行Bonferroni校正,只需要點選「校正」選項中的「Bonferroni」核選項目, ANOVA 結果報表將會出現另一列以Bonferroni校正調整的 p 值( 表 13.8 )。比較一下這三個 p 值與未校正的成對 t 檢驗的 p 值,會發現 jamovi 只是將這些數值乘以 \(3\)而已。
13.5.4 Holm校正
雖然Bonferroni校正非常簡單,但並不是最好的校正方法。另一種常用的方法是Holm校正(Holm, 1979)。這種方法背後的思路是假設按照p 值的大小,從最小的比較開始,調整到第j個比較的p值為止。校正p值的公式是:
\[p_j^{'}=j \times p_j\]
這就是說最大的 p 值保持不變,第二大的 p 值增加兩倍,第三大的 p 值增加三倍,依此類推。或者用以下等式挑出較大者。
\[p_j^{'}=p_{j+1}^{'}\]
這樣解釋不大好懂,讓我慢慢說明Holm校正的原理。首先將\(m\)組檢定的p值按照數值大小排序,從最小的排到最大的。最小的 p 值校正程序只有乘以 \(m\),其他組p值的校正程序則有兩個步驟。像是第二小的 p 值首先要乘以 \(m - 1\),若是這個校正後的p值大於前一次校正後的p值( \((m-1)\times p_{m-1} > m \times p_m\) ),就可以停止校正;若是小於前一次校正後的p值,就要將原來的p值,改換為前一次校正後的p值。 表 13.10 展示五組多重比較檢定的p值,如何用Holm校正來調整每組檢定的p值。
raw p | rank j | p \( \times \) j | Holm p |
---|---|---|---|
.001 | 5 | .005 | .005 |
.005 | 4 | .020 | .020 |
.019 | 3 | .057 | .057 |
.022 | 2 | .044 | .057 |
.103 | 1 | .103 | .103 |
雖然校正程序稍微麻煩,Holm校正有一些優點,像是保證比Bonferroni校正有更低的型二錯誤率,進而提高統計考驗力。此外,Holm校正保障每組檢定的型一錯誤率一致,雖然這套方法執行下來看起來不大像。要跑統計的研究者,搞清楚兩種校正方法的原理後,多數人會支持Holm校正。所以對於初學統計的讀者來說,Holm校正應該是首選的校正方法。 圖 13.4 展示Holm校正的 p 值,最大的 p 值( Anxifree vs. 安慰劑)和未校正的p值一樣,都是.15。最小的 p 值(Joyzepam vs. 安慰劑)則校正為原來的三倍。
13.5.5 事後檢定的報告格式
根據事後比較,確定是那幾組有真正的顯著差異,我們可以來寫正式報告:
以Holm校正調整後 p 值的事後檢定顯示,與 Anxifree(p = .001)和安慰劑(\((p = 9.0 \times{10^{-5}}\))相比,Joyzepam 有更顯著的情緒改善效果。此次研究沒有發現 Anxifree 的效果比安慰劑好的證據(\(p = .15\))。
若是不想報告精確的 p 值,可以將數值改為 \(p < .01\)、\(p < .001\) 和 \(p > .05\)。無論哪種表達方式,關鍵是要說明p 值是用Holm校正調整 。當然,這段報告之外的其他部分,應該要包括組平均值和標準差的描述統計資訊,因為 p 值並不是什麼充分的資訊。
13.6 單因子變異數分析的適用條件
變異數分析也和其他統計方法一樣,需要確認資料性質符合幾個適用條件,尤其是殘差。這裡特別談三個關鍵條件:常態性、變異同質性和獨立性。那要如何檢查殘差的適用條件呢?以下分別討論這三個條件的檢測方法。
[額外的技術細節 12]
- 變異同質性。因為會用單因子變異數分析只會有一個母群標準差\(\sigma\),各分組資料不會有各自的母群標準差\(\sigma_k\)。設定各分組的標準差都與總標準差相等,就是變異同質性的真義,因此又被稱為等變異性。檢核變異同質性會詳細說明如何檢測。
- 常態性。殘差必定符合常態分佈。如同 小單元 11.9 介紹過的,可以察看Q-Q圖,或執行Shapiro-Wilk檢定來做評估。進一步說明可見檢核常態性。
- 獨立性。解釋這個條件有點棘手。基本意思是,任何一筆觀察值的殘差與其他觀察值的殘差彼此之間毫無關聯。符合獨立性的同一組資料,所有資料點的 \(\epsilon_{ik}\) 來源不存在任何交集。這個條件沒有現成或簡易的檢核方法,只能從研究設計檢視有無違反獨立性的狀況(參考 單元 2 )。不存在獨立性的最經典狀況是重複測量設計,每個研究參與者會在多個條件中被測量,如此一來就違反獨立性的條件,不可能以單因子變異數分析處理資料。如果一項研究具備這樣的設計,就需要使用單因子重覆量數變異數分析。
13.6.1 檢核變異同質性
檢核變異數分析的第一步,就像坐大型遊輪出海,要仔細觀察海面是不是平靜無波,能讓遊輪安全離港!
– 喬治·博克斯 (Box, 1953)
據說殺貓有千百種方法,檢核變異同質性也有千百種方法。當然在動物保護意識極高的現代社會,這種話不能亂說。依照原作者曾閱讀的文獻,最常見到的檢核方法是Levene檢定法(Levene, 1960),以及Brown-Forsythe檢定法(Brown & Forsythe, 1974)。
無論是Levene還是Brown-Forsythe,檢定統計值\(F\)或\(W\)的計算程序都和ANOVA的F統計值一樣,只是結果變項的符號是用\(Z_{ik}\)代表而不是\(Y_{ik}\)。知道了這些,我們就可以了解如何使用jamovi做檢定。
[額外的技術細節13]
13.6.2 jamovi的Levene檢定
如何進行Levene檢驗呢?其實程序很簡單 - 在ANOVA模組設定選單的”條件檢定”(Assumptions Check)選項裡,點選”變異數同質性檢定”。接著查看報表區的輸出,可以看到檢定結果,如同 圖 13.5 ,顯示無顯著差異(\(F_{2,15} = 1.45, p = .266\)),所以看起來符合變異數同質性沒有問題。然而,只看外表可能會讓人上當!如果資料樣本量相當大,即使違反變異數同質性的程度沒有到影響ANOVA的穩健性,Levene檢定也可能顯示顯著結果(p < .05)。這正是George Box的論文中所指出的重點。若是資料的樣本量相當小,那麼可能無法滿足變異數同質性,使得Levene檢定結果不顯著(p > .05)。也就是說,使用任何統計方法檢定適合條件是否滿足時,也要檢查分組統計圖的標準差視覺化範圖,親自判斷各組變異範圍是否是相等或不相等。範圍相等才算是符合變異數同質性。
13.6.3 校正異質性的分析結果
以這個單元的範例分析來說,符合變異數同質性是沒有問題的:Levene檢定顯示無顯著結果,雖然更謹慎的做法是還要看看分組統計圖的變異分佈。不過現實的統計實務不會總是如此幸運,若是分組殘差的變異違反變異數同質性,要如何拯救ANOVA呢?回想一下 小單元 11.4.3 ,t檢定也會遇到同樣的問題。Student t檢定要預設分組變異相等,解決方法是使用不需要依賴這項條件的Welch t檢定。其實, Welch (1951) 也提到如何單因子變異數分析的版本:Welch單因素檢定。jamovi的ANOVA模組包括的One-Way ANOVA有執行此方法的功能,開啟設定選單後,勾選Welch檢定,如同 圖 13.6 的展示。圖 13.6 的ANOVA報表標示Fisher的檢定結果,與使用jamovi進行變異數分析這一節展示的結果相同:\(F(2, 15) = 18.611, p = .00009\) 。
之前的ANOVA分析結果是\(F(2, 15) = 18.6\),Welch單因素檢定的計算結果是\(F(2, 9.49) = 26.32\)。換句話說,Welch檢定將組內自由度從15降低到9.49,F值從18.6上升到26.32。
13.6.4 檢核常態性
常態性的檢核方法相對簡單, 小單元 11.9 已經介紹需要執行的工作。最簡單的做法是畫張Q-Q圖,可以的話就跑個Shapiro-Wilk檢定。以 圖 13.7 的Q-Q圖來看,這份資料相當符合常態性。若是Shapiro-Wilk檢定結果也是不顯著(\(p > .05\)),就能確認這筆資料符合常態性。不過要注意的是,若是樣本量很大,那麼顯著的Shapiro-Wilk檢定結果很可能是偽陽性,也就是說這筆資料即使違反常態性,並不會對分析結果造成任何問題。同樣地,樣本量太小可能會造成偽陰的檢定結果。這就是為什麼需要Q-Q圖輔助檢核。
這份範例資料殘差用Q-Q圖來看,並沒有明顯偏離常態的資料,Shapiro-Wilk檢定結果並未顯示顯著,p = 0.6053,報表詳見 圖 13.6 。因此,兩種檢核方法都沒有發現明顯違反常態性的證據。
13.6.5 排除非常態性的分析結果
除了檢核常態性,也要學會如何處理違反常態性的資料。以單因子ANOVA的來說,最簡單的解決方法是改用不必在意任何適用條件的無母數統計方法。 單元 11 已經介紹分析兩組平均值的無母數統計方法:Mann-Whitney以及Wilcoxon檢定,若要分析三組以上的平均值,可以使用Kruskal-Wallis秩和檢定法(Kruskal & Wallis, 1952)。
13.6.6 Kruskal-Wallis檢定的運算原理
Kruskal-Wallis檢定與單因子ANOVA有許多相似之處。ANOVA的計算是從第k組的第i個觀察對象,標定在結果變項的數值\(Y_{ik}\)開始。Kruskal-Wallis檢定的第一步是排序所有\(Y_{ik}\),以資料的序位進行分析。14
13.6.7 更多分析細節
前一節說明Kruskal-Wallis檢定的運算原理,我認為這是一種思考統計檢定如何運作的正確方法。15
不過,這裡要跟讀者抱歉一下,還要交待一個使用限制!什麼限制?因為前面的範例是原始資料沒有任何兩個觀察值的數值相等,如此才能使用Kruskal-Wallis檢定。也就是說,如果資料內存在相同的數值,就必須引進校正因子才能計算。由於非常有耐心的讀者看到這裡,也許再也不想讀下去了,若是真的遇到問題,再去查什麼是tie-correction factor(TFC)就好了。這裡只用範例資料做個簡單的說明,以及解釋為何可以暫視無視這個限制。以下根據clinicaltrails的原始資料變項mood.gain建立一個對應觀察值出現次數的資料表單,以\(f_j\)表示每個觀察值數值在資料裡的出現次數,下標\(j\)代表資料數值的排序。全部例出後就如 表 13.11 的第二列。
0.1 | 0.2 | 0.3 | 0.4 | 0.5 | 0.6 | 0.8 | 0.9 | 1.1 | 1.2 | 1.3 | 1.4 | 1.7 | 1.8 |
---|---|---|---|---|---|---|---|---|---|---|---|---|---|
1 | 1 | 2 | 1 | 1 | 2 | 1 | 1 | 1 | 1 | 2 | 2 | 1 | 1 |
次數表排序第三的情緒改善分數是0.3,出現次數是2,表示所有參與者裡有兩位的情緒改善是0.3分。16
所以,我們不用太擔心使用限制,jamovi輸出的Kruskall-Wallis統計值是校正處理後的。我們終於走完Kruskal-Wallis檢定的學習之路了,各位也了解為什麼不需要在意手上的資料有沒有使用限制的理由了吧?
13.6.8 使用jamovi完成Kruskal-Wallis檢定
儘管學會Kruskal Wallis檢定法是讓許多學生望而生畏的任務,實際操作卻相當輕鬆,因為jamovi的ANOVA模組提供“Non-parametic”(無母數)-“One-way ANOVA(Kruskall-Wallis)”的分析功能。類似clinicaltrial.csv的大多數資料,會有結果變項mood.gain
和分組變項drug
。只要符合這樣的結構,就能用jamovi執行Kruskall-Wallis檢定,得到 圖 13.8 展示的結果,寫成符合學術格式的報告\(\chi^2 =12.076, df = 2, p = 0.00239\)。
13.7 單因子重覆量數變異數分析
單因子重覆量數變數分析是一種用於每位參與者有參與三種以上實驗條件,或者各組參與者有密切相依的測量條件,有三個以上組平均差異要比較的狀況。所以各實驗條件應該有數量相等的觀察值。搭配這種研究設計的分析也可叫做「相依樣本變異數分析」或「參與者內變異數分析」。
重覆量數變數分析的運算原理與獨立樣本ANOVA非常像,後者有時也叫做「參與者間變異數分析」。之前討論參與者間變異數分析時提到組間離均差平方和(\(SS_b\)),以及組內離均差平方和(\(SS_w\)),各自除以對應的自由度就會得到\(MS_b\)和\(MS_w\),如同 表 13.1 。F統計值的計算公式是:
\[F=\frac{MS_b}{MS_w}\]
重覆量數變異數分析的F統計值的計算方法也是差不多,只是\(SS_w\)要再分成兩部分,不像獨立樣本ANOVA直接用\(SS_w\)做為\(MS_w\)的分母。因為每一組的參與者都是同樣的人,可以從組內離均差平方和移出參與者間個別差異\(SS_{subjects}\)。在此不深談如何移出的細節,只要想像每位參與者都是名為「參與者」這個變項的其中一個變項水準。只要將\(SS_w\)減掉\(SS_{subjects}\),就能得到一個較小的\(SS_{error}\):
\[\text{獨立樣本變異數分析: } SS_{error} = SS_w\] \[\text{重覆量數變異數分析: } SS_{error} = SS_w - SS_{subjects}\] 減少的\(SS_{error}\)可以增加統計檢定的考驗力,不過還要看減少的\(SS_{error}\)有沒有比減少的自由度更多,因為獨立組的參與者較多,重覆量數的自由度從\((n - k)\)17變為\((n - 1)(k - 1)\)
13.7.1 重覆量數變異數分析的jamovi實作
這裡用 Geschwind (1972) 的研究資料做示範。這個研究的對象是中風後出現語言缺陷的患者,研究者想了解導致語言缺陷的大腦損傷區域,找來六位確診Broca失語者的患者,研究資料列在 表 13.12 。
Participant | Speech | Conceptual | Syntax |
---|---|---|---|
1 | 8 | 7 | 6 |
2 | 7 | 8 | 6 |
3 | 9 | 5 | 3 |
4 | 5 | 4 | 5 |
5 | 6 | 6 | 2 |
6 | 8 | 7 | 4 |
患者需要完成三項語言能力檢核任務。在第一項言語生成任務,患者需要重複朗讀研究者給他的單詞。第二項概念理解任務要測試理解意義的能力,患者需要將一系列圖片與正確名稱配對。第三項句子修正任務測試文法知識,患者要對語法不正確的句子重新排列字詞順序。每位患者都完成了所有任務,各項任務的順序在參與者之間進行平衡隨機。每個任務有10次嘗試。每位患者成功完成的嘗試次數列在 表 13.11 。我們將這些資料輸入jamovi進行分析(或者載入Broca.csv資料檔)。
使用jamovi執行單因子重覆量數ANOVA,只要點選ANOVA - Repeated Measures ANOVA,就能開啟如 圖 13.9 的分析設定選單。
依序進行以下設定:
- 更改重複測量的獨變項名稱(RM Factor 1)。獨變項是涵括您所選擇的資料變項標籤,用於描述所有參與者被重複測量的條件。例如,要描述所有參與者完成的言語、概念和語法任務,一個合適的名稱是“任務”(Task)。請注意,這個新的名稱代表變異數分析的獨變項。
- 在重複測量因子對話框(Repeated Measures Factors)裡填入三個級別名稱,因為三個級別分別代表三個任務:言語、概念和語法。有需要的話可再更改級別的名稱。
- 然後將代表每個級別的資料變項,移動到重複測量因子欄位框(Repated Measures Cells)。
- 最後,在“適用條件檢核”(Assumptions Check)選項裡,勾選“球形性檢查”(Sphericity Test)。
報表區會輸出如 圖 13.10 到 圖 13.13 的單因子重覆量數變數分析報表。建議各位先查看有沒有符合Mauchly球形檢定( 圖 13.10 ),這是檢核各測量條件之間的變異是否相等(也就是研究條件之間的得分差異值分佈大致相同)。 圖 13.10 顯示Mauchly檢定的\(p = .720\),未低於一般的顯著水準(p > .05),我們有理由相信各測量條件之間的差異值分佈無顯著差異,符合球性分佈的適用條件。
另一方面,如果Mauchly檢定結果是顯著的話(p < .05),表示各測量條件之間的差異值分佈有顯著差異,不符合球性分佈的適用條件。這麼一來,我們必須校正單因子重覆量數變異數分析的F值:
- 如果”球形檢定”輸出的Greenhouse-Geisser值大於 .75,應該使用Huynh-Feldt法校正。
- 如果輸出的Greenhouse-Geisser值小於.75,應該使用Greenhouse-Geisser法校正。
這兩個校正的F值都可以用勾選“適用條件檢核”(Assumptions Check)選單裡的球形修正選項來指定,修正的F值將顯示在變異數分析結果表,如同 圖 13.11 。
根據這個範例資料的分析報表,可以看到Mauchly的球形檢定並不顯著(p = .720,即p > 0.05)。因此,這表示我們可以接受已滿足資料差異的分佈符合球形假設,不需要校正F值。報告可以用原始分析的結果說明”任務”的主要效果:\(F = 6.93\),\(df = 2\),\(p = .013\),所以三種任務正確率的差異,確實是因為各種任務的測試內容所造成的(\(F(2, 10) = 6.93\),\(p = .013\))。
與獨立樣本變異數分析一樣,jamovi也可以執行重覆量數變數分析的事後檢定。示範結果如同 圖 13.12 。這顯示言語生成任務與句子修正任務的正確率有顯著差異,但是其他各項任務之間的正確率沒有顯著差異。
運用描述統計可以幫助解釋結果,只要勾選邊際平均值選項,生成如同 圖 13.13 的表格。比較不同任務的平均正確次數,可以看出布洛卡失語症患者的言語生成任務(平均= 7.17)和概念理解任務(平均= 6.17)表現相對較好。不過,句子修正任務的表現明顯較差(平均= 4.33),事後檢定顯示言語生成和句子修正的正確率有明顯差異。
13.8 Friedman無母數重覆量數變異數分析
Friedman檢定是單因子重覆量數變異數分析的無母數版本,用來分析三組或更多組排序資料之間的差異,只要每組參與者相同,或者各條件的參與者彼此有密切相依的測量條件。如果依變項是序列數值,或者資料未滿足常態性條件,就可以使用這種方法進行假設檢定。
和Kruskall-Wallis檢定相同,徹底解釋的話需要先知道複雜的數學知識,所以這裡不會詳細介紹。對於本書的目標讀者來說,只要了解jamovi輸出的Friedman檢定報表是已經校正的結果, 圖 13.14 展示的報表是用失語症患者資料做Friedman檢定的範例。
使用jamovi執行Friedman檢定,只要從功能模組選單,選擇”ANOVA” - Repeated Measures ANOVA (Non-parametric),就能開啟如同 圖 13.14 的選單。然後將要比較的重複測量資料變項名稱(言語、概念、語法)放到“測量”(Measures)對話框。想要生成描述統計表格,只要勾選描述(Descriptives)。
jamovi報表區顯示描述統計、卡方統計值、自由度和p值。由於p值小於通常用於一般的顯著水準(p < .05),根據報表結果可以結論:布洛卡失語症患者在言語生成(中位數= 7.5)和概念理解(中位數= 6.5)的表現較好。然而,他們的句子修正任務表現明顯較差(中位數= 4.5),事後檢定顯示言語生成和句子修正任務的表現有顯著差異。
13.9 變異數分析與t檢定的關係
這個單元結束之前,還有一件事需要知道。初學者通常會對此感到訝異,不過了解後的收穫是值得的:以ANOVA和學生t檢定分析兩組組間差異的資料,結果是一樣的。不只輸出的結果相似,各種資訊的統計意義也是等價的。這裡用實際的示範取代公式推導讓各位明白。以clinical trials資料集來當例子,這次不納入drug
,只用therapy
做唯一獨變項。執行ANOVA會得到一個F統計值 \(F(1,16) = 1.71\),和p值 = \(0.21\)。因為只有兩組,也可以做學生t檢定。t檢定的結果得到t統計值 \(t(16) = -1.3068\) 和 \(p = 0.21\)。你應該發現,p值是相同的,都是\(p = .21\)。那麼統計值有什麼關係?因為t統計值是負的,顯然和ANOVA的F值不一樣。其實其中有相當直接的轉換關係。只要將t統計值平方,就會得到F統計值:\(-1.3068^{2} = 1.7077\)。
13.10 單元小結
這一章份量不少,但是有一些細節我並未提到18。最明顯的是在此並未討論處理不只一個分組變項的資料,下一個 單元 14 會討論其中一部分。這個單元的學習重點有:
理解單因子變異數分析的運算原理 以及使用jamovi的變異數分析模組
學習如何計算變異數分析的效果量
若是各組的的觀察值數目相等,這樣的研究設計就是“平衡設計”。以這個單位介紹的單因子變異數分析來說,設計是否平衡並不重要。不過要進行較複雜的變異數分析運算,設計是否平衡大有關係。↩︎
所以總離均差平方和的公式與變異數的公式長得幾乎一樣 \[SS_{tot}=\sum_{k=1}^{G} \sum_{i=1}^{N_k} (Y_{ik} - \bar{Y})^2\]↩︎
總離均差平方和能被分解成兩種離均差,帶來計算的便利性。第一種離均差是組內離均差平方和,是每個資料與所屬分組平均值之間差異\[SS_{w}= \sum_{k=1}^{G} \sum_{i=1}^{N_k} (Y_{ik} - \bar{Y}_k)^2\] \(\bar{Y}_k\) 是某一分組平均值,在此代表服用第 k 種藥物的參與者情緒狀態改善狀況,也就是並非比較每個人的數值與所有參與者的總平均值,而是各分組之內的參與者彼此比較。所以\(SS_w\)的值一定會小於總離均差平方和,因為完全不包括分組之間的差異,也就是各種藥物的效果。↩︎
這一步要知道的是如何計算組間離均差平方和 \[ \begin{aligned} SS_{b} &= \sum_{k=1}^{G} \sum_{i=1}^{N_k} ( \bar{Y}_{k} - \bar{Y} )^2 \\ &= \sum_{k=1}^{G} N_k ( \bar{Y}_{k} - \bar{Y} )^2 \end{aligned} \]↩︎
在獨立樣本變異數分的報告裡, \(SS_w\) 也被稱為誤差總和 \(SS_{error}\)。↩︎
追根究底來說,ANOVA 就是兩種不同的統計模型,\(H_0\)和\(H_1\),的對決過程。這個單元一開始提到虛無假設和對立假設所用描述方式,其實不夠精確。在此要做點補救,儘管這樣會讓一些讀者感到厭煩。回憶一下,這個單元設定的虛無假設是各分組平均值彼此相等。若是如此,拆解結果變項 \(Y_{ik}\) 的方式應該是將個別數值當成唯一的母群平均值\(\mu\),加上因為測量方式,實際數值與母群平均值之間的偏差,通常會用\(\epsilon_{ik}\)表示,學術界習慣稱為誤差或殘差。不過請小心用詞,如同“顯著”這個詞義的歷史演變,“誤差”的詞義在統計學的場域不完全等於日常生活的模糊含義。日常用語的“誤差”暗示某種錯誤,但是統計學所指的不完全是任何一種錯誤。因此在統計學場域用“殘差”,會比“誤差”不易讓人誤會。這兩個詞都表示“剩餘變異”,也就是模型無法解釋的“部分”。不論是那種稱呼,將虛無假設寫成統計模型的話,看起來就是 \[Y_{ik}=\mu+\epsilon_{ik}\] 本書稍後會討論到,通常會設定各分組的殘差值 \(\epsilon_{ik}\)皆符合平均值為 \(0\),標準差為 \(\sigma\)的常態分佈。運用從機率入門學到的符號,殘差的數學模型可以寫成 \[\epsilon_{ik} \sim Normal(0,\sigma^2)\] 那要如何設定對立假設 \(H_1\) 呢?虛無假設和對立假設的唯一區別是,研究者認為各組的母群平均值並不相等。所以只要設定\(\mu_k\) 表示第 k 組的母群平均值,\(H_1\) 的統計模型就是 \[Y_{ik}=\mu_k+\epsilon_{ik}\] 同樣地,其中的分組殘差也都是平均值為 \(0\),標準差為 \(\sigma\)的常態分佈。也就是說,對立假設還帶上 \(\epsilon_{ik} \sim Normal(0,\sigma^2)\) 這個條件。
現在\(H_0\) 和 \(H_1\)的統計模型都正式登場,可以清楚說明如何評估離均差平方和平均,以及如何解釋\(F\)的意思。這裡不需要做任何數學證明,我們只要知道,組內離均差平方和平均 \(MS_w\) 可以當成各組殘差的變異數估計值。組間離均差平方和平均\(MS_b\)也是一組估計值,除了各組殘差的變異數估計值,還有各組平均值之間的差異值。以\(Q\)代表這項差異值,就能將F統計數的公式改寫為 \(^a\) \[F=\frac{\hat{Q}+\hat{\sigma}^2}{\hat{\sigma}^2}\] 若是虛無假設的推測成立,\(Q = 0\);若是對立假設的推測成立,\(Q < 0\)(請參考 Hays (1994) ,ch. 10)
這就是為什麼在變異數分析的狀況,\(F\) 值必須大於 1,是拒絕虛無假設的基本條件。讀者要注意的是,這不是說F 值不可能小於 1。嚴格地說,如果虛無假設的推測成立,估計值的取樣分佈會逼近期望值是1的F分佈 \(^b\),所以F 值必須要大於 1,才能安全地拒絕虛無假設。
從上一段的F值公式可知,\(MS_b\)和\(MS_w\),都有殘差 \(\epsilon_{ik}\) 變異數的估計值,如果虛無假設的推測成立,兩者就只有殘差變異數的估計值。若殘差符合常態分佈,根據 小單元 7.6 學過的卡方分佈,讀者會想到\(\epsilon_{ik}\)的變異數估計值其實應該符合卡方分佈。這就是卡方分佈的真正面貌:無數個符合常態分佈的隨機變數平方後,累加形成的機率分佈。F分佈則是兩件符合卡方分佈的隨機變數之比值,集合形成的機率分佈。當然,以上說明省略許多細節,不過這是至此介紹過的各種統計檢定方法,所依據的取樣分佈來源。
—
\(^a\) 如果讀者已經讀過 單元 14 ,學到如何用 \(\alpha_k\) 定義因子內水準 k 的「處理效應」(參考 小單元 14.1 及 小單元 14.2 ),會知道 \(Q\) 是處理效應平方的加權平均值,\(Q = \frac{(\sum_{k=1}^{G}N_k \alpha_k^2)}{(G-1)}\)
\(^b\) 或者是更精確的F分佈期望值 \(1+\frac{2}{df_2-2}\)。↩︎更真實的描述應該是 “1899年的統計學家著迷於計算的主要原因,是因為現實生活沒有什麼朋友,也沒有比算術更有趣的事情可做。而且直到 1920 年,ANOVA 才正式降臨人類世界。”↩︎
使用 Excel 計算完整的\(SS_b\),會因為四捨五入得到不大一樣的數值,有可能是 \(3.45\)。↩︎
與手動計算的數值相比,jamovi 的計算結果更為準確,這是因為四捨五入的差異。↩︎
若是研究開始前已經有些理論,預期某幾組比較是確實有意義的,那就是另一回事了。這種情況所執行多重檢定,其實不是「事後分析」,而是「預先計劃的比較」。 小單元 14.9 有更進一步介紹。↩︎
順帶一提,並不是所有校正方法都是這樣搞,因為Bonferroni校正是著眼於控制「整體的型一錯誤率」(family wise Type I error rate)。有些校正方法則是控制「偽發現率」(False Discovery Rate)。↩︎
若是讀者沒有真正親手演練單因子變異數分析的運算原理的範例,最好至少要讀過一遍。這個範例裡提到支撐ANOVA運算流程的統計模型:\[H_0:Y_{ik}=\mu + \epsilon_{ik}\] \[H_1:Y_{ik}=\mu_k + \epsilon_{ik}\]其中\(\mu\)是指處於所有分組平均值中心的總平均之母群期望值,\(\mu_k\)是指第k組平均值的母群期望值。至此一直反覆討論的問題是,研究資料是最適配虛無假設代表的單一總平均,還是對應假設代表的分組平均值。這些討論是有必要的,因為統計模型的設定就是研究問題的符號化!但是在檢定過程的有效性其實是建立在殘差 \(\epsilon_{ik}\) 符合常態分佈的程度,也就是\[\epsilon_{ik} \sim Normal(0,\sigma^2)\]。如果沒有符合這個條件,所有ANOVA的計算步驟都無法執行。也可以說,就算進行所有計算過程,最後得到了F統計值,但是違反適用條件造成這個F統計值並不能代表研究結果支持任何一個模型的指標,最後無法提出有科學意義的結論。↩︎
Levene檢定法非常簡單,先設定結果變項\(Y_{ik}\)。第一步是定義一個新變項\(Z_{ik}\),表示各組觀察值與該組平均值的偏誤 \[Z_{ik}=Y_{ik}-\bar{Y}_{k}\] 製造新變項有什麼好處?讓我們花一點時間來思考一下\(Z_{ik}\)到底是什麼,以及我們要檢檢什麼。\(Z_{ik}\)代表量測第\(k\)組的第\(i\)次觀察得到的數值,與其組平均值的誤差。變異數分析的虛無假設設定各分組的平均值相等,因此各組平均值的總變異應該也是相等!所以Levene檢定法的虛無假設是設定各組的\(Z\)之母群平均值相等。嗯,那麼有什麼統計方法能用來檢核這個變項?沒錯,就是ANOVA,所以Levene檢定法的流程就是以\(Z_{ik}\)為依變項做分析的ANOVA。那Brown-Forsythe檢定法是怎麼做?有什麼特別之處嗎?其實與Levene檢定法的唯一不同是建立新變項\(Z\)的方式,是用分組中位數計算所有觀察值的誤差分數而不是用分組平均值。所以說,Brown-Forsythe檢定法建立新變項的公式是 \[Z_{ik}=Y_{ik}-median_k(Y)\] \(median_k(Y)\)代表第k組的中位數。↩︎
實作方法是先創建新的結果變項\(R_{ik}\),表示第k組的第i位成員排序。接著計算第k組觀察值的平均排序 \(\bar{R}_k\): \[\bar{R}_k=\frac{1}{N_k}\sum_i R_{ik}\],以及總平均排序\(\bar{R}\): \[\bar{R}=\frac{1}{N}\sum_i\sum_k R_{ik}\] 完成後,就可以計算與總平均排名\(\bar{R}\)的離均差平方。\((R_{ik} - \bar{R})^2\)代表第ik個觀察值排序與總平均排序的離均差,是一種無母數的計量;還有\((R_{ik} - \bar{R_k})^2\)代表第k組平均排序與總平均排序的離均差,也是一種無母數的計量。因此,離均差平方\((R_{ik} - \bar{R})^2\)就是總體誤差的無母數計量。如此一來,就能按照ANOVA的運算原理,定義排序變異的估計值。首先是“總排序平方和”\[RSS_{tot}=\sum_k\sum_i (R_{ik}-\bar{R})^2\],以及“組間排序平方和” \[\begin{aligned} RSS_{b}& =\sum{k}\sum_{i}(\bar{R}_{k}-\bar{R})^2 \\ &= \sum_{k} N_k (\bar{R}_{k}-\bar{R})^2 \end{aligned}\] 若是虛無假設成立,且根本沒有真正的組間差異,組間排序和\(RSS_b\)應該非常微小,遠小於總排序和\(RSS_{tot}\)。就計算公式來看,這和計算ANOVA 的F統計值非常像,不過因為背後的取樣分佈不同,Kruskal-Wallis檢定統計值記號是K,而K的計算方法是\[K=(N-1) \times \frac{RSS_b}{RSS_{tot}}\] 若是虛無假設成立,K的取樣分佈會逼近自由度為\(G-1\)的卡方分佈。 K的值越大,資料就越不適配虛無假設,因此這是一個單側檢定。K值夠大的話,就有機會拒絕\(H_0\)。↩︎
不過,從純數學的觀點來看,這樣的複雜解說是不必要的。雖然這本書不會談如何推導公式,讀者可以使用一些線性代數\(^b\)來推想K的公式,為什麼可以轉換為K統計值\[K=\frac{12}{N(N-1)}\sum_k N_k \bar{R}_k^2 -3(N+1)\] K的運算比在前面描述的F要容易得多,但是對於實際操作統計方法的分析人員來說完全沒有意義,將K類比為分析排序資料的ANOVA統計值可能是最好的解讀方法。但是請記住,檢定統計值K與前面的ANOVA統計值F是很不一樣的公式。
—
\(b\)就只是數學運算符號。↩︎根據上一段介紹的數學方法,我們可以知道\(f_3 = 2\)。想通這一點,就可以了解TCF的公式是用來算什麼了:\[TCF=1-\frac{\sum_j f_j^3 - f_j}{N^3 - N}\]jamovi計算出來的K值,正是除以TCF之後,輸出在Kruskal-Wallis檢定報表。↩︎
(n-k)=(參與者數量-組別數量)↩︎
就像其他章節,這個單元有許多參考來源,其中原作者參考最多的專書是 Sahai & Ageel (2000) 。這本書對初學者來說偏難,不過如果學到這裡,想知道更多變異數分析的數學原理,這本書是不錯的參考資源。↩︎