做數(shù)據(jù)分析這些年單因素方差分析one-way ANOVA是我用得最多、也最容易被用錯的一個統(tǒng)計工具。無論是最常見的三組工藝參數(shù)對比還是評估四個地區(qū)用戶行為的差異第一步都是先做一個整體比較確認(rèn)組間是否存在系統(tǒng)性差異而不是一上來就兩兩檢驗。這篇筆記就把單因素方差分析的理論脈絡(luò)完整梳理一遍連同Matlab里的落地代碼一起講清楚。適合看這篇內(nèi)容的讀者有兩類一類正在學(xué)運籌優(yōu)化或數(shù)據(jù)分析需要把課本上的平方和、F統(tǒng)計量真正和代碼對應(yīng)起來另一類是做實驗評估的工程師手里的數(shù)據(jù)已經(jīng)跑出來了但不確定該用哪一種檢驗或者不知道anova1輸出的表格到底說明了什么。讀完之后你會搞清楚方差分析為什么是在“比方差”也會拿到一套可以直接照著改的分析流程。1. 單因素方差分析它到底在解決什么問題1.1 三組以上均值比較為什么不能直接兩兩t檢驗先說一個很多新手容易踩的坑。當(dāng)我們要比較三組數(shù)據(jù)的均值是否有差異直覺上最容易想到的做法是把三組拆成三對分別做兩個樣本的t檢驗比如A對B、A對C、B對C各算各的。這個思路本身沒錯錯在忽略了“多次比較”帶來的錯誤累積。假設(shè)我們比較的是三組樣本實際三組總體均值完全一樣我們用5%的顯著性水平做單次t檢驗時誤判為“有差異”的概率是0.05。但如果做三次兩兩t檢驗?zāi)呐旅拷M本身沒有差異至少出現(xiàn)一次“假陽性”的概率就變成了 1 - 0.953 ≈ 14.3%。組數(shù)越多這個膨脹越嚴(yán)重五組就要做十次比較假陽性概率接近40%。也就是說你會越來越容易“發(fā)現(xiàn)”根本不存在的差異這在統(tǒng)計上叫第一類錯誤膨脹。方差分析的處理思路是先做一個“總體的、全局的”檢驗把所有的組放在一起看組間均值之間的差異是否已經(jīng)大到了不太可能用隨機(jī)誤差來解釋。這一道檢驗通過之后再做兩兩比較這才是在控制整體錯誤率的前提下繼續(xù)定位差異的合理路徑。說得直白一點方差分析的作用是先判斷“到底有沒有區(qū)別”而兩兩比較回答的是“是誰和誰有區(qū)別”。1.2 一個典型的運籌優(yōu)化場景與基本判定思路用工程里的例子來說明最直觀。假設(shè)我們有三臺同型號加工設(shè)備想判斷它們在產(chǎn)出零部件的強(qiáng)度上是否存在顯著差異。每臺設(shè)備各取5個零件測量強(qiáng)度得到三組數(shù)據(jù)。表面上看三組均值一定不會完全相等肯定有大有小可問題在于這種均值差異到底是“機(jī)器決定的”還是“零件本身的隨機(jī)波動造成的”方差分析給了一個非常清晰的判定框架先把數(shù)據(jù)里的變異拆成兩個來源。一個叫組間變異反映的是各組的平均水平與總平均水平的偏離程度另一個叫組內(nèi)變異反映的是同一組內(nèi)部各個觀測值之間的波動。如果三臺機(jī)器對強(qiáng)度沒有本質(zhì)影響那么組間變異應(yīng)該和組內(nèi)變異差不多兩者的比值接近1如果機(jī)器確實決定了強(qiáng)度差異組間變異就會明顯大于組內(nèi)變異兩者的比值會遠(yuǎn)遠(yuǎn)超過1。這個比值就是F統(tǒng)計量我們用它的分布來判斷“大到了什么程度才算顯著”。之所以在這里反復(fù)強(qiáng)調(diào)“變異”而不是“均值差異”因為方差分析的本質(zhì)不是直接比較均值大小而是比較不同來源的方差是否處于同一個水平。理解了這一點后續(xù)的自由度、均方、F值這些概念就順理成章了。2. 理論拆解從變異分解到F檢驗2.1 組間變異與組內(nèi)變異方差分析的核心邏輯假設(shè)我們研究一個因素比如“工藝條件”它有k個水平這里k≥3。第i個水平下有ni個觀測值用xij表示第i組第j個數(shù)據(jù)。記所有數(shù)據(jù)的總體均值為x?第i組均值為x?i。一個觀測值和總平均值的偏差可以寫成下面的等式(xij - x?) (x?i - x?) (xij - x?i)這個式子看著簡單含義卻很深總偏差被分成了兩部分。第一部分反應(yīng)的是“組的效應(yīng)”即這個組的平均水平偏離整體有多遠(yuǎn)第二部分是“隨機(jī)誤差”即這個具體觀測值偏離自己所在組的平均水平有多遠(yuǎn)。把等式兩邊平方并累加就得到一個標(biāo)準(zhǔn)的關(guān)系式SST SSA SSE。SST是總平方和代表所有數(shù)據(jù)的整體變異SSA是組間平方和代表各組均值之間的差異SSE是組內(nèi)平方和代表各組內(nèi)部的隨機(jī)波動。這里的“拆分”是整個方差分析大廈的地基。為什么能這么拆因為交叉項在求和時會互相抵消具體到數(shù)學(xué)上就是等式右側(cè)的交叉項和為零這在所有教材里都有推導(dǎo)。理解這個分解并不需要很深的數(shù)學(xué)背景只要抓住一個核心方差分析關(guān)注的是“組間差異相對組內(nèi)噪聲是不是足夠大”。2.2 一張表格完成全部推導(dǎo)數(shù)值示例演算光講公式容易發(fā)虛我們拿一組實際數(shù)據(jù)完整演算一遍。還是用上面三臺設(shè)備的例子三組數(shù)據(jù)如下組1設(shè)備A82 86 80 84 85組2設(shè)備B90 92 88 94 91組3設(shè)備C78 80 76 82 79每組5個樣本總樣本量 n15。先算各組均值和總均值組1均值(8286808485)/5 83.4組2均值(9092889491)/5 91.0組3均值(7880768279)/5 79.0總均值所有數(shù)據(jù)相加除以15 84.47組間平方和的計算比較直接每組樣本量乘以該組均值與總均值偏差的平方再累加SSA 5×(83.4-84.47)2 5×(91.0-84.47)2 5×(79.0-84.47)2 368.53組內(nèi)平方和要把每組內(nèi)部5個數(shù)據(jù)分別減去自己組的均值再平方累加SSE (第一組內(nèi)部偏差平方和) (第二組內(nèi)部偏差平方和) (第三組內(nèi)部偏差平方和) 23.2 20.0 20.0 63.20于是總平方和 SST 368.53 63.20 431.73。自由度也要分開考慮組間自由度為 k-12組內(nèi)自由度為 n-k12總自由度為 n-114。有了平方和與自由度就可以把每個來源的變異標(biāo)準(zhǔn)化為均方組間均方 MSA SSA / (k-1) 368.53 / 2 184.27組內(nèi)均方 MSE SSE / (n-k) 63.20 / 12 5.27F統(tǒng)計量就是這兩個均方的比值F MSA / MSE 184.27 / 5.27 ≈ 34.98查F分布臨界值表在5%顯著性水平下F(2,12)對應(yīng)的臨界值約3.89。我們算出來的F值遠(yuǎn)大于臨界值說明組間變異相對于組內(nèi)隨機(jī)波動來說確實非常大三臺設(shè)備的零件強(qiáng)度存在顯著差異。最終的標(biāo)準(zhǔn)操作是把這些數(shù)值整理成一張方差分析表這也是Matlab輸出的格式來源。變異來源平方和SS自由度df均方MSF值顯著性組間368.532184.2734.98p0.05組內(nèi)63.20125.27總和431.73142.3 三個重要假設(shè)以及假設(shè)被破壞時的信號任何統(tǒng)計檢驗都有前提條件方差分析也不例外。單因素方差分析有三個基本假設(shè)各組數(shù)據(jù)相互獨立、各組數(shù)據(jù)近似服從正態(tài)分布、各組的方差齊性即方差相差不大。這三條里獨立性通常由實驗設(shè)計來保證比如樣本隨機(jī)抽取、不同設(shè)備之間互不影響正態(tài)性可以通過正態(tài)概率圖或檢驗來評估方差齊性是使用F檢驗時比較關(guān)鍵的假設(shè)。如果方差不齊F值會受到影響極端情況下會把“本不顯著”判成“顯著”或者反過來把“顯著”判成“不顯著”。判斷方差不齊最簡單的方法是用箱線圖看各組離散程度也可以用Levene檢驗或Bartlett檢驗做形式化判斷。遇到不滿足假設(shè)的數(shù)據(jù)最常見的處理方向有兩個一是改用不需要正態(tài)和方差齊性假設(shè)的非參數(shù)方法比如Kruskal-Wallis檢驗二是先做變量變換比如對數(shù)變換或平方根變換讓數(shù)據(jù)更接近正態(tài)和方差的同質(zhì)性要求。這兩條路后面都會展示具體做法。3. Matlab實現(xiàn)從數(shù)據(jù)準(zhǔn)備到完整分析鏈路3.1 數(shù)據(jù)組織和anova1的常見調(diào)用方式Matlab里做單因素方差分析的核心函數(shù)是anova1全稱是“one-way analysis of variance”。它的功能就是把前面說的那些平方和、自由度、F值全部計算出來并自動繪制箱線圖和方差分析表窗口。只要數(shù)據(jù)組織正確通常幾行代碼就能完成。假設(shè)我們要分析的就是上面那個三組設(shè)備數(shù)據(jù)在Matlab里第一種也是最簡單的方式是把三組數(shù)據(jù)拼成一個5行3列的矩陣每列代表一個水平% 三組數(shù)據(jù)分別存入變量 process1 [82 86 80 84 85]; process2 [90 92 88 94 91]; process3 [78 80 76 82 79]; % 拼成矩陣每一列是一個處理組 X [process1(:), process2(:), process3(:)]; % 執(zhí)行單因素方差分析 p anova1(X);這個調(diào)用方式要求矩陣的每一列恰好對應(yīng)一個組列數(shù)就是水平數(shù)k行數(shù)就是每組樣本量并且各組樣本量要相等。如果各組樣本量不一樣就不能用矩陣形式得用另一種“分組變量”形式把所有數(shù)據(jù)放在一個列向量里再用另一個列向量指明每個數(shù)據(jù)屬于哪一組% 長格式數(shù)據(jù)行數(shù)等于總樣本量 allData [82; 86; 80; 84; 85; 90; 92; 88; 94; 91; 78; 80; 76; 82; 79]; groupID [1; 1; 1; 1; 1; 2; 2; 2; 2; 2; 3; 3; 3; 3; 3]; % 運行方差分析返回p值、方差分析表、統(tǒng)計量結(jié)構(gòu)體 [p, tbl, stats] anova1(allData, groupID);anova1的函數(shù)簽名有好幾個可選參數(shù)完整格式是[p, tbl, stats] anova1(x, group, displayopt)。第三個參數(shù)displayopt可以設(shè)為off來關(guān)閉自動彈出的圖形窗口在批處理腳本里很實用。我建議平時保留默認(rèn)的顯示因為彈出來的方差分析表窗口能直觀看到中間計算結(jié)果對于核對數(shù)據(jù)是否組織正確很有幫助。3.2 方差分析表與p值的正確讀法調(diào)用anova1(X)之后Matlab會彈出一個圖形窗口里面包含左側(cè)的箱線圖和下方的方差分析表同時在命令窗口返回p值。初學(xué)者往往只看p值是否小于0.05這其實浪費了輸出結(jié)果。使用三輸出版本[p, tbl, stats] anova1(...)可以拿到完整的tbl把它顯示出來就能看到和手工推導(dǎo)時結(jié)構(gòu)一致的內(nèi)容[p, tbl, stats] anova1(X); disp(tbl);tbl是一個cell數(shù)組第一行是表頭第二行是組間信息第三行是組內(nèi)信息第四行是總和。具體字段包括Source變異來源、SS平方和、df自由度、MS均方、FF統(tǒng)計量、ProbFp值。這里的p值表示如果原假設(shè)為真所有組均值相同得到當(dāng)前這么大F值的概率。p值越小說明數(shù)據(jù)支持“組間存在差異”的證據(jù)越強(qiáng)。這里有一個容易忽略的點p值并不是“三組均值相同的概率”它只是一個條件概率指標(biāo)經(jīng)常被誤解。我通常這樣向同事解釋p值回答的是“如果各組真的沒有差異出現(xiàn)這么極端數(shù)據(jù)的概率有多大”而不是“各組沒有差異的概率有多大”。當(dāng)p小于0.05時我們只是認(rèn)為當(dāng)前的實驗數(shù)據(jù)與原假設(shè)不太兼容從而傾向于拒絕原假設(shè)。統(tǒng)計量結(jié)構(gòu)體stats包含各組均值、組內(nèi)方差等字段這個結(jié)構(gòu)體后面做兩兩比較時要用到先存起來。3.3 兩兩比較multcompare的進(jìn)一步定位如果方差分析整體檢驗的結(jié)果是顯著說明“至少有一對之間存在差異”但并不知道是哪幾對。這時需要調(diào)用多重比較函數(shù)multcompare它的輸入正是anova1返回的stats結(jié)構(gòu)體% 假定已經(jīng)執(zhí)行過 [p, tbl, stats] anova1(X); [c, m, h, gnames] multcompare(stats);輸出中c是一個矩陣每一行對應(yīng)一組兩兩比較格式為 [組號i, 組號j, 置信區(qū)間下限, 均值差估計, 置信區(qū)間上限, 調(diào)整后的p值]。判斷兩兩是否顯著的方法是看置信區(qū)間是否包含0如果整個區(qū)間都在0的同側(cè)即上下限同號說明這兩個組的均值差異顯著如果區(qū)間跨過了0則說明差異不顯著。這種方法好處是不僅能知道“是否顯著”還能直觀看到差異可能的大小范圍。multcompare默認(rèn)使用Tukey HSD方法適合數(shù)據(jù)量不特別大、各組樣本量接近的場景。實際操作中如果想用其他校正方法可以指定參數(shù)比如multcompare(stats, Alpha, 0.01)把置信水平改為99%或者用CType參數(shù)切換成Bonferroni、Dunn-Sidak等。默認(rèn)的Tukey在絕大多數(shù)情況下夠用不必過度糾結(jié)。3.4 非參數(shù)補(bǔ)充方案Kruskal-Wallis的適用時機(jī)上面說過如果數(shù)據(jù)不滿足正態(tài)性或方差齊性的假設(shè)替代方案之一是采用Kruskal-Wallis檢驗。Matlab里對應(yīng)的是kruskalwallis函數(shù)基本用法和anova1幾乎一致[p_kw, tbl_kw, stats_kw] kruskalwallis(X);這個檢驗的原理是把所有數(shù)據(jù)混合在一起排秩然后看各組秩和是否存在顯著差異。它不需要正態(tài)假設(shè)對方差齊性的容忍度也更高所以常用于數(shù)據(jù)分布嚴(yán)重偏斜、有離群值或者樣本量很小的場景。需要提醒的是Kruskal-Wallis檢驗比較的是各組分布位置是否一致當(dāng)各組分布形狀差異很大時即使均值沒有差異也可能得到顯著性結(jié)果這時要結(jié)合箱線圖一起解釋。還有一個選擇是Welch方差分析它通過對非整數(shù)自由度的校正來放寬方差齊性假設(shè)。Matlab基礎(chǔ)工具箱里沒有現(xiàn)成的函數(shù)但可以通過文件交換區(qū)找到第三方實現(xiàn)或者自己按公式寫一個。工程上優(yōu)先用Kruskal-Wallis因為實現(xiàn)成本更低、結(jié)果更穩(wěn)健。4. 實際數(shù)據(jù)分析中的高頻問題與排查心得4.1 方差不齊時的兩個替代路徑在真正跑數(shù)據(jù)的時候最常碰到的問題就是方差不齊。我自己的習(xí)慣是anova1跑完之前先看一眼各組箱線圖的盒子寬度如果差別太大再用vartestn做一個形式化檢驗[p_var, stats_var] vartestn(X);如果p_var小于0.05說明方差齊性假設(shè)存疑。這時候有兩條路可以走。第一條是改用Welch校正或Kruskal-Wallis檢驗這在前一小節(jié)已經(jīng)說了第二條是嘗試變量變換比如對原始數(shù)據(jù)取對數(shù)后再做方差分析。變換的意義在于把原始尺度上的異方差性壓縮掉讓數(shù)據(jù)更接近“方差穩(wěn)定”的狀態(tài)。比如強(qiáng)度數(shù)據(jù)本身可能和均值呈正相關(guān)對數(shù)變換后往往能緩解這個問題。我不建議在方差不齊時直接忽略它跑anova1因為F檢驗對方差偏離比較敏感。如果你已經(jīng)跑了結(jié)果又貼著0.05邊界解釋起來會非常被動。穩(wěn)妥的辦法是同時跑標(biāo)準(zhǔn)方差分析和兩種替代分析如果三者的結(jié)論一致就說明結(jié)論不是某個檢驗假設(shè)造成的偶然現(xiàn)象。4.2 樣本量不平衡怎么處理單因素方差分析并不要求各組樣本量完全相同但樣本量差異過大會帶來兩個問題統(tǒng)計功效下降以及多重比較結(jié)果的可信度降低。在實驗設(shè)計階段盡量讓各組樣本量接近這種“平衡設(shè)計”是最理想的狀態(tài)。如果已經(jīng)拿到了不平衡數(shù)據(jù)比如一組10個、一組5個、一組6個anova1一樣能運行前提是使用分組變量的長格式傳入數(shù)據(jù)而不是單純拼矩陣。不平衡數(shù)據(jù)做多重比較時默認(rèn)的Tukey HSD方法會做一些調(diào)整因為各組標(biāo)準(zhǔn)誤不同。此時輸出結(jié)果中的置信區(qū)間已經(jīng)考慮了樣本量的影響但解釋時要格外注意差異較小的那組估計值可能方差很大。經(jīng)驗上如果樣本量差異超過1倍以上我更傾向在報告中把分析結(jié)果和樣本量分布一起展示避免誤導(dǎo)。4.3 數(shù)據(jù)非正態(tài)與離群值的影響正態(tài)性假設(shè)被違反時anova1在樣本量足夠大的情況下往往仍然可靠這是因為中心極限定理讓組均值的抽樣分布趨向正態(tài)。但如果樣本量很小又存在嚴(yán)重偏態(tài)或離群值結(jié)果就可能失真。我的習(xí)慣是先畫箱線圖人眼觀察再用lillietest或jbtest做檢驗h1 lillietest(process1); h2 lillietest(process2); h3 lillietest(process3);lillietest的返回值h為1時表示拒絕正態(tài)性假設(shè)。如果只有一兩組偏離正態(tài)問題倒不大如果三組都明顯偏態(tài)就要認(rèn)真考慮非參數(shù)方案了。離群值更麻煩因為它對均值和方差都有拉動作用。處理離群值時我會先檢查是否為錄入錯誤而不是直接剔除。如果數(shù)據(jù)本身沒問題但確實是異常點可以做“剔除后重跑對比分析”的敏感性檢查跑一次包含離群值的結(jié)果再跑一次剔除后的結(jié)果如果結(jié)論不變就無所謂如果結(jié)論變了說明這個異常點主導(dǎo)了分析結(jié)論在報告里要如實說明。4.4 顯著性和實際意義要分開看效應(yīng)量的引入還有一個經(jīng)常被忽視的問題統(tǒng)計顯著性不等于實際重要性。樣本量很大的時候即使各組均值差異很小也可能算出很小的p值因為組內(nèi)均方的標(biāo)準(zhǔn)誤被大樣本壓得很低。反過來樣本量很小時即使差異在業(yè)務(wù)上有意義也可能因為功效不足而不顯著。為了解決這個矛盾可以補(bǔ)充計算效應(yīng)量。單因素方差分析中最常用的效應(yīng)量是η2eta squared計算公式為η2 SSA / SST以上面那組數(shù)據(jù)為例η2 368.53 / 431.73 ≈ 0.85說明總變異中有大約85%可以被組間差異解釋這是一個非常大的效應(yīng)。一般經(jīng)驗是η2在0.01附近算小效應(yīng)0.06附近算中等0.14以上算大效應(yīng)。報告p值的同時補(bǔ)上效應(yīng)量會讓分析結(jié)論更有說服力。這里也回應(yīng)了上一節(jié)提到的“p值貼線”問題如果p值恰好落在0.05附近效應(yīng)量往往能給出更穩(wěn)定的判斷依據(jù)。我不太贊成把0.05當(dāng)成一道鐵閘門0.049和0.051在實際結(jié)論上可能并沒有本質(zhì)差別。5. 做單因素方差分析的一些個人經(jīng)驗5.1 先畫圖再檢驗用boxplot快速發(fā)現(xiàn)問題任何一個正式的方差分析我都會先畫箱線圖。Matlab里最簡單的畫法是直接在anova1(X)彈出箱線圖之前先用boxplot畫一張獨立的圖boxplot(X, Labels, {設(shè)備A, 設(shè)備B, 設(shè)備C});箱線圖能一次暴露好幾個問題各組中位數(shù)是否差得很遠(yuǎn)、盒子寬度是否明顯不一致、是否存在離群點。如果在看圖階段就發(fā)現(xiàn)問題比如某組的離群點多到異常就不用急著跑方差分析了先弄清楚數(shù)據(jù)為什么長這樣。要把方差分析當(dāng)成“驗證工具”而不是“探索工具”探索階段永遠(yuǎn)是圖形優(yōu)先。5.2 報告結(jié)果時建議給出的三個信息如果要在實驗報告或論文里匯報單因素方差分析的結(jié)果我建議至少包含如下三個信息F統(tǒng)計量的具體值、對應(yīng)的自由度、p值。例如寫“單因素方差分析結(jié)果表明三種工藝條件下零件強(qiáng)度的差異顯著F(2, 12) 34.98p 0.001”。這樣寫的優(yōu)點是信息完整別人可以根據(jù)自由度驗證你確實做了三組×五樣本的分析。如果要進(jìn)一步給出兩兩比較結(jié)果可以把multcompare輸出的置信區(qū)間畫成圖或者用表格列出每對比較的均值差和置信區(qū)間。不要只寫“A組顯著高于B組”要給出數(shù)值依據(jù)和分析方法比如“Tukey HSD多重比較顯示設(shè)備B的均值比設(shè)備C高12.095%置信區(qū)間為[9.3, 14.7]”。這樣才算一個可復(fù)現(xiàn)的分析結(jié)論。5.3 寫在最后的一個小技巧擾動樣本驗證穩(wěn)定性最后分享一個我自己的土辦法。分析做完之后如果p值不是特別極端比如在0.01到0.1之間我會把原始數(shù)據(jù)稍微擾動一下比如每個觀測值加減一個小隨機(jī)數(shù)然后重新跑一遍方差分析。如果結(jié)論在擾動后依然穩(wěn)定那說明結(jié)果比較可靠如果稍加擾動p值就從顯著變成不顯著那說明數(shù)據(jù)本身的可信度要打問號結(jié)論需要更謹(jǐn)慎地表達(dá)。這個辦法算不上嚴(yán)格統(tǒng)計方法但用來快速評估分析結(jié)論的“韌性”很實用。做方差分析這幾年最大的體會是統(tǒng)計工具的作用不是替你做決定而是幫你看清數(shù)據(jù)的真實結(jié)構(gòu)。理解平方和分解的邏輯、讀懂F統(tǒng)計量比機(jī)械地記幾個Matlab函數(shù)重要得多。遇到問題多從假設(shè)條件出發(fā)排查比盲目套用各種檢驗更有用。希望這篇從理論到代碼的梳理能幫你在自己的數(shù)據(jù)上少走幾步彎路。