問題概要是關於 ANOVA 的事後比較 (post-hoc) 在多個組別 (levels) 之間依捉對比較（pair-wise）差異是否顯著的結果進行分組的問題。
但事情沒有這麼簡單，如果只是普通的 ANOVA 事後比較，無論是手動推演或是使用程式 （例如 SAS PORC GLM 項下 LSMEANS 中的 LINES 陳述） 都可以輕易解決，可這次問題的趣味點在於「不等組」(unequal sample sizes) 設計。不等組設計最大的麻煩就在於不同組之間的捉對比較，因為 n1, n2 組合的不同而造成信賴區間估計有所差異，當樣本數越懸殊，顯著性就不怎麼依從分組平均的排序了。這樣講一堆想必大家都聽不太懂，接著來給他詳細解釋一下。
但若是你對 ANOVA、post-hoc（事後比較）、甚至 multiple paired-comparison 不清楚的朋友，這篇文章大概就不這麼實用了。
至於看到 ANOVA 與事後比較（post-hoc）就驚恐萬分的朋友也不用驚慌，其實本篇文章用到 ANOVA 與事後比較的概念幾乎只有 MSE 估計。而事後比較（post-hoc）簡單來說也不過是從好幾組平均中捉對比較，看看彼此之間是否有顯著差異（在校正後的 alpha 下）而已。而文章的脈絡也聚焦在介紹一個演算法。
廢話不多說立刻來舉例，假設現在有 4 組平均數（1~4）進行事後比較，而實驗設計上是完全等組（每組樣本數相當），則在事後比較中每組將共用一個相同、由整體模型估計出來的 MSE，假如修正後 alpha 分配也相同，則理論上每組的信賴區間（confidence interval）也會等寬（如圖1）。
圖1
從圖一可以看得出來，判斷兩個組平均之間的差異是否顯著的標準若以信賴區間是否有重疊來看的話，我們可以得到以下配對結果：1-2同組、2-3同組、3-4顯著、1-3顯著、1-4顯著、2-4顯著。但事實上，前面的判斷只需要到1-3顯著就可以結束了。為什麼呢？因為在信賴區間等寬，各組平均又有“順序性”的情形下（1→2→3→4），當3-4顯著被確認後，組平均1及2當然與組平均4也會達差異顯著。所以我們可以知道，在完全等組設計的 ANOVA 下進行手動的事後比較，可以依照組平均排序後線性地一一檢查並完成分組，沒有太大的困難。
但是在圖1的四組平均數分組上，遇到一個問題：組平均2分別與組平均1、3各成一組，如此一來要怎麼表示呢？這裡引入線段表示法以及一個常用的符號區別系統，後者規則是：
「在組平均後面標示符號，任意兩個組平均之間若有共通的符號即表示彼此差異不顯著」
圖2
如圖2，我們可以看到用「線段」表示後更清楚，最精簡的情形下一共可分成三組：[1,2]、[2,3]、[4]。而這三組又可以分別以 a、b、c 標示。所以我們常在學術論文中看到類似的表達方式（如表1）。
表1
無論用線段或是符號標示，是不是都比原來的兩兩配對描述清楚得多呢。
好，回到正題。但若討論到不等組ANOVA設計下的事後比較，事情就沒有這麼簡單了。例如組1有30個受試者、組2只有5個受試者，這樣便成立不等組的條件。尤其是當組間樣本數差異越大時，事後比較就沒辦法參考前面提到的組平均順序性了。因為這時候針對各組的信賴區間估計就不再等寬。如圖4範例表示。
圖4
從圖4來看，以配對描述試說明可以得到：1-2顯著、2-3同組、3-4同組、1-3同組、1-4顯著、2-4顯著。我們明顯發現因為組平均3的信賴區間異常地寬，導致分組變得有點奇怪。除此之外，各組別的信賴區間寬度也是不相等的。
讓我們更清楚些用前面提到的線段及符號分組法表示，結果如圖5。
圖5
從圖5來看，我們可以發現除了組平均3跟大家都“有一腿”以外，其它三個組平均之間是沒有任何交集的，也就是差異達顯著。更有意思的是組平均1、3雖然可以放在同一組中，但卻跳過了理應介於他們之間的組平均2，形成標號為a的這一組線段有中斷的狀況。符號標示的結果可以整理如表2。
表2
簡單小結一下，不等組 ANOVA 在進行事後比較時，因為各組平均的信賴區間估計寬度可能不同，進而導致分組結果與組平均大小順序之間沒有直接關係。使用線段法也許可以稍微清楚地表示關係，但是當組數越來越多時，線段圖也會越複雜甚至難以繪製。如前面舉例，線段甚至會在同一組中出現中斷的情形。而剩下唯一可以保證正確描述分組關係的只有兩兩捉對比較的配對描述。符號標示法雖然清楚，但是因為兩兩配對描述的結果可能很複雜，如何手動算出“最精簡”的符號標示分組法便成了很麻煩的問題。目前也不清楚有哪套統計軟體能解決這個問題。
以 SAS 為例，當你丟入一筆不等組設計的 ANOVA 資料並進行 post-hoc 時，系統會提醒你這是一個不等組的檢定，而由於 SAS 的 LINES 指令只能繪製連續線段（SAS 會直接以符號連在一起，同時呈現“線段”與“符號”），所以在不等組的情形下跑出來的結果相當存疑。如圖6，SAS 警示提醒使用者這是不等組 ANOVA 的結果。
圖6、SAS statement LSMEANS LINES 結果，以 TUKEY 為例
所以本篇文章就是要介紹 YJJ 想出的演算法，程式化以後順利解決了這個問題：給定任意 k 組平均數的兩兩捉對比較顯著結果，計算出最精簡的符號標示分組結果。
我們可以開放式地來思考一下！如果在不等組 ANOVA 設計下，各組平均間的大小順序再也沒有參考性，那麼可能的分組關係該如何描述呢？事實上，因為組平均的分組不再適用線性排序的方式做檢查，我們可以直接將各組平均視為在平面上、甚至空間中獨立的「點」，而點與點之間的關係以線段連線表示（線段相連的點可以視為同一組、不相連的表示兩者平均差異達統計顯著）。我們可以想像這樣的分組結果可能五花八門、自由度極高，已經是一個排列組合的概念了。
讓我們馬上以六組平均間的關係來舉例，如圖7。
圖7
從圖7範例中可以看到平面（空間）中的六組平均，以線段連線表示兩兩之間的關係（線段相連表示差異不顯著，可視為同組）。這個圖可以勉強看出兩兩配對比較的結果：[1,2]、[2,5]、[3,4,5]、[6] 應該是最精簡的分組法。
有趣的是組平均 3、4、5 形成一個互不顯著的小組、而組平均 6 與其它組皆達顯著差異。這樣的表示法也許還算清楚，但是若遇到更多組別、形成交疊的網狀線段時，同樣很難人工判別。
我們繼續把圖7的關係圖化成熟悉的線段及符號標示分組法，結果如圖8。
圖8
好的，這裡要介紹一個簡單的概念，用來把“兩兩配對比較的（顯著與否）結果描述”以數學的方式表達，整合成一個「顯著矩陣」的概念，例如剛剛六組平均之間的關係可以表示如圖9。矩陣的下三角中，「1」表示兩組間差異達統計顯著、「0」表示不顯著且可視為同組。這樣的表達是否清爽且清楚多了呢？更能簡單地檢查任意兩組之間關係。值得注意的是，在不等組 ANOVA 設計的情形下，矩陣中任一元素都可以為 1 或 0，不受到組平均大小排序的影響，這也是剛剛提到的分組自由度。YJJ 我的目的就是想要達成：任意一個 k 組的顯著矩陣被輸入以後，可以得出最精簡符號分組結果的演算法。
圖9
好啦！看到這裡不知道你有沒有頭昏眼花。如果你不是剛好對統計學或基本演算法很有興趣，也不是剛好遇到像我朋友那樣有 10 組平均（共 45 個配對比較）要進行分組，通常應該已經點選關閉網頁了。但如果你堅持要繼續看下去， YJJ 在這裡提醒你一下，如果前面的“問題”部分沒有弄清楚，是不建議往下繼續看演算法的喔！而就算你已經看懂題目在問什麼，也清楚線段與符號標示分組法的內容，也建議你可以休息一下，想想看你的解法，再往下看 YJJ 的解法。這裡讓我放一張最近最夯的“香蕉花”美圖，給眼睛休息一下。
看完“香蕉花”，接著就直接進入解決方法的說明囉～
演算法邏輯概要：
1. k組會有一共 [k*(k-1)]/2 組兩兩比較的結果，YJJ 另外增加「自己與自己比」的概念，所以一共會有 [k*(k+1)]/2 組。
（eg. k=10 就有 45 種兩兩比較，再加自己比自己的 10 種等於 55 種)
***這裡定義兩兩相比的結果：1 為顯著、0 為不顯著。
2. 計算 k 組平均的任意“子分組”個數（每個子分組最少包含 1 個組平均，最多包含 k 個組平均），其所有可能性，包含單 1 組平均自己獨立為一組，一共有 2^k-1 種子分組（減1為扣除空集合）。我們把它通通列出來。
3. 將所有 2^k-1 個子分組與 [k*(k+1)]/2 個兩兩比較顯著與否的結果進行關係檢查，會發現任意一個子分組只會對應到特定一種兩兩比較的結果。將這樣的結果製成一個檢查表。這個檢查表的製作方式很簡單，待會舉例詳述。
（eg. 極端來看，若全部 k 個組平均都彼此達差異顯著，也就是最簡潔的分組法是 k 個平均數各自獨立為一組的結果，其對應到的兩兩比較結果即為：除了 k 個自己與自己比的結果為 0 以外，其它所有兩兩相比的結果都為 1）
4. 製作第二個參照用的檢查表，檢查特定子分組成立時，可以被它包含的其它子分組有哪些。數學上可以做個交叉表，可以被包含就標示 1，不行就標示 0。
（所謂的“被包含”舉個例如：當 [1,2,3] 這種子分組成立的時候，其它子分組如 [1,2]、[2,3]、[1,3]、[1]、[2]、[3] 這六個同樣也可以成立，但卻沒有存在的必要。為了符合最精簡的分組法，[1,2,3] 這個子分組成立後，被它包含的所有子分組就無需被描述了）
5. 把特定“顯著矩陣”放入步驟2產生的檢查表中，對應兩兩比較顯著性，把所有合法成立的子分組標示出來（例如合法成立標示為 1、無法成立則標示為 0）。這裡的演算法稍稍取巧了些，待會舉例詳述。
（所謂的合法成立，並不代表“最精簡分組”，而是指這樣一個特定的子分組不違背被輸入的特定顯著矩陣之組間關係。例如，有四組平均，其相關矩陣圖10暗示了最精簡分組為兩組：[1,2,3]、[4]。則可以合法成立的子分組應該有：[1]、[2]、[3]、[4]、[1,2]、[2,3]、[1,3]、[1,2,3]。你可以試試看把後面那八個子分組分別標上 a~h，並以符號標示法列出來，會發現雖然有很多多餘的符號，但結果卻完全不違反圖10的顯著矩陣，如表3）
圖10
表3
6. 把步驟5產生的所有成立之子分組與步驟4的包含關係檢查表相互比照，任何一個合法成立的子分組，只要被任一其它也成立的子分組所包含，則它就沒有存在的必要。只要被包含到一次以上，就標示為 0，而合法成立的子分組中未被其它合法子分組包含者標示為 1。
7. 透過步驟6，所有合法成立的子分組中，依然被標示為 1 的子分組就是「最精簡的分組法」。
前面是概要，後面是落落長的說明。如果你前面看不懂沒關係，搭配後面的說明以後就會看懂了。而前面的概要只是作為未來喚醒記憶的省時間線索而已。但如果後面的說明看完後還是看不懂前面概要的朋友，大概是我中文表達能力不佳，只好說聲抱歉囉。
一、步驟1＆步驟2的說明
如果這兩個步驟就卡住的話，這篇文章可能就不太適合不熟悉數學的你喔。步驟1是梯形公式計算兩兩捉對一共有幾對，只是另外加上自己與自己比較的 k 對（概念上對應“獨立自己為1組”的子分組）。步驟二則是單純的排列組合，特定子分組中每一個組平均都有存在或不存在兩種可能，所以 k 組平均進行分組，所有可能的子分組會是 2^k 種，但減去空集合 1 種。
二、步驟3，子分組與兩兩配對比較的關係檢查表
這裡請直接看圖11的表示。YJJ 以 k=10 組為例，前面文章是用1234當組別，圖11的舉例是用ABCD當組別，只是寫程式時候習慣問題，小心別混淆囉。
直列部分列出了共 1024 種子分組，橫欄部分列出了 55 個配對比較（包含自己比自己，如 A-A、B-B、C-C，千萬別忘記）。而直列旁邊簡單用函數表示出該子分組包含的組別有哪些？而要怎麼判斷某個不顯著的配對是否屬於特定子分組呢？只要直列的英文字母群包含了橫欄的兩個字母，就表示這一對配對字母的不顯著是此子分組成立的必要條件。在 excel 裡面用 if 與 countif 函數就可以完成判斷。（紅格子處可以看到強迫空集合不成立）
圖11
三、步驟4，製作子分組間包含關係檢查表
又是個看圖說故事的時間，請看圖12的範例。同樣是 k=10 情形下的配對比較，這時候直列與橫欄都各別是 1024 組。而包含與否的判斷很簡單！只要橫欄的子分組字母群有被直列的子分組所包含，那麼兩者的交叉值就會顯示為 1。這裡同樣是使用 if 與 countif 就輕易完成。同要要小心空集合強迫不成立，以及主對角線強迫為 0 的設計。主對角線如果有 1，就是子分組自己包含自己，對後續的判斷會造成影響，但看個人程式設計習慣決定。
圖12
四、步驟5，輸入特定顯著矩陣，檢查哪些子分組可合法成立
同樣以 k=10 為例，被輸入的特定顯著矩陣，有意義的兩兩比較值，包含自己比自己一共有 55 個值（0同組、1顯著）。這 55 個值剛好全數對應到圖11橫欄上的 55 對比較。若圖11每個子分組對應到的配對標示方式相同（0同組、1顯著）。那麼只要「特定子分組的 55 個值都“大於或等於”特定被輸入之顯著矩陣對應的值」，這個子分組就是合法成立的。這裡的邏輯可能需要想一下才會理解，但就留給讀者去思考了。
五、步驟6＆步驟7，把合法成立的組別套入步驟4產生的檢查表中
目的是檢查所有合法成立的子分組是否有被其它合法成立的子分組所包含。這裡的判斷其實 YJJ 已經偷偷寫在圖12之中了。而只要使用一些個人設計的巧思，這個判斷式其實很簡單就完成。重點是保留合法成立子分組中未被其它成立組包含的組別。這也是為什麼前面提到子分組自己包含自己的主對角線 YJJ 要故意保留為 0 了，避免自己被自己包含。 至於被保留的子分組就是解答囉！
如何？看到這裡是否還清楚呢？希望你也跟我一樣享受這個思考的過程。
YJJ 這裡提供一個用 excel 寫出來，適用 k 小於等於 10 的分組試算表供參考使用。因為用試算表寫的關係，檔案有點肥。以後再用 C/C++ 或 VBA 寫通用版本，現在好懶。不過看文章的你可以簡單想想看，為什麼 k=10 的試算表也適用小於 10 組的分組，而且不用修改。
適用 k=10 以下的分組試算表：點我下載
老王賣瓜後，YJJ 我出於謙遜的本性（？）來檢討一下這個演算法的優缺點。
優點大概就是簡單易懂，直接用 excel 就可以完成撰寫，連 VBA 的功能都用不到。
缺點顯而易見地是出於這個演算法的發想原初就是「索引式」的方向。利用簡單邏輯去建立中介的參考表，再用簡單的判斷式去逐步推出結果。但是建立索引的過程中，因為考量到所有可能性的關係，資料量隨著比較的組數增多呈現指數成長（罪魁禍首 → 2^k-1）。根據實測大概在 k=15 以內，近五年的筆電都有辦法在 1 分鐘內完成。但是當 k=20 時候，可能性就突破百萬筆了。所以 YJJ 在這裏不負責任建議，一般電腦如果要使用這個演算法最好實驗設計的單一變項組數不要超過 20 levels。雖然目前我也沒看過這種層級的實驗就是了。
最後是雖然 YJJ 不是數學本科系所，但深深覺得這個演算法的原初議題應該是有「數學解」的（相對於索引法）。因為這篇落落長的演算法問題，其實根本兩句話就可以講完題目：
「將 k 個單元兩兩間關係製成一個主對角線為 0 的對稱方陣（symmetric square matrix），其餘元素皆為 1 或 0。試計算此 k 個單元最精簡的分組方式，使得組內各單元兩兩間對應的矩陣元素值皆為 0。」
如果你有興趣看完這個文章，又剛好知道數學解的話，歡迎隨時告訴小的我囉，我會非常感謝滴～
沒有留言:
張貼留言
歡迎留下你對兩位作者文章的任何看法，或是寫信給我們。