劉應興的部落格

假設 Y(i1), ..., Y(in_i), i = 1, ..., k, 是相互獨立分別來自 N(μ_i, σ^2_i) 的樣本資料，也就是說分別自 k 個常態群體獨立抽取大小為 n_i 的隨機樣本，總樣本大小為 n = Σ n_i。我們關心 k 群體之變異數 σ^2_1, ..., ο^2_n 是否相等，

H0: σ^2_1 = . . . = σ^2_k = σ^2,   Ha: some σ^2_i ≠ σ^2_j

概度比檢定 (likelihood ratio test) 統計量是

λ  = {Π_i (s'^2_i)^{n_i/2}}/(s'^2)^{n/2}

式中 s'^2_i = Σ_j (Y(ij) - Y(i.))^2/n_i, s'^2 = Σ_i Σ_j (Y(ij) - Y(i.))^2/n，而其中 Y(i.) = Σ_j Y(ij)/n_i，即樣本 i 之樣本平均數。與一般「樣本變異數」相比，

s'^2_i = (n_i - 1) s^2_i/n_i,   s'^2 = Σ_i (n_i - 1)s^2_i/n = (n-k)s^2/n

此處 s^2_i 是第 i 樣本變異數，s^2 是諸 s^2_i 以其自由度 n_1 - 1 加權的平均，即 H0 成立時共同 σ^2 的聯合不偏估計。拒絕域是 λ < λ*, 當 k = 2 時，相當於 F = s^2_1/s^2_2 偏高與偏低的部分。當 H0 成立時，前項 F 統計量服從 F(n_1-1, n_2-1)，也就是自由度 n_1 - 1, n_2 - 1 的 F 分布。這就是一般初級統計教程所謂的「兩群體變異比檢定」。

對於 k > 2, 概度比統計量一般取對數形式:

- 2 ㏑(λ) = n ㏑(s'^2) - Σ_i n_i ㏑(s'^2_i) 

習慣上是取 s^2_i 代替 s'^2_i, 結果和上列 -2 ㏑(λ) 相差只是一個與諸 n_i 有關的常數項。Bartlett 檢定則是以

T = (n - k) ㏑(s^2) - Σ (n_i - 1)㏑(s^2_i)

為基礎，實際是採用

χ^2 = T/{1 + [Σ_i 1/(n_i - 1) - 1/(n-k)]/[3(k-1)]}

與自由度 k - 1 的卡方分布右尾臨界值比較。-2 ㏑(λ) 統計量在大樣本中常用卡方近似，其自由度依概度比大樣本理論是 k - 1；而 Bartlett 檢定是修正自前項概度比轉換的 -2 ㏑(λ)，使更近似 χ^2(k-1)，原始論文為 Bartlett, M. S. (1937). "Properties of sufficiency and statistical tests". Proceedings of the Royal Statistical Society, Series A 160, 268–282.

如同兩群體變異比之 F 檢定，Bartlett 檢定對群體的非常態性很敏感，一方面是群體非常態則型Ⅰ誤機率可能遠大於名目顯著水準；另方面是如果檢定結果呈現顯著差異，可能只是因為群體非常態，而不是諸群體之變異數不等。取而代之的，Levene 檢定或 Brown-Forsythe 檢定是對群體非常態性較不敏慼的檢定方法。這兩種檢定很相似，都是把原資料轉換成離中（均）差絕對值

(Levene)  Z(ij) = |Y(ij) - Y(i.)|;   (B-F)  Z(ij) = |Y(ij) - Y*(i,.)|

式中 Y(i.) 如前述是第 i 樣本平均數，而 Y*(i,.) 是第 i 樣本中位數，或某個比例的截尾平均數 (trimmed mean，或譯：裁剪、修剪、截斷平均數) 。然後對上列轉換後資料做變異數分析，用 F 檢定。兩者何者較適當與群體真實分布有關， Brown and Forsythe (1974) 的蒙地卡羅模擬發現高狹峰（厚尾）分布群體如 Cauchy 分布，以截尾平均取代簡單樣本平均或中位數為計算離差的中心最佳；而高度偏斜如 4 個自由度的卡方分布群體，則是樣本中位數較佳；而原始的 Levene 檢定則適用於群體分布比較對稱，峰度也不太高，也就是尾巴比較普通的情形。

有個問題是 Z(ij) 並不符合 ANOVA 的各項假設，以 Levene 檢定為例，諸離均差 Zij, j = 1, ..., n_i, i = 1, ..., k, 非常態分布固不用說，相同 i 之下各 Z(ij) 相互並不獨立，不同 i 之間諸 Z(ij) 也不是等幅變異。非常態性的問題，由於 ANOVA 的 F 檢定是兩獨立樣本 t 檢定的多群體（多樣本）版本，如 t 檢定一般具有穩健性 (robustness)，也就是說對非常態性不那麼敏感。由於 Y(ij), Z(ij) 是半常態分布，即常態分布從中心折半，故尾巴不會太厚（峰度不太高；雖有偏態，其偏態系數也就相當於 8 或 9 自由度的卡方分布。這意謂因取絕對值導致的非常態性在 ANOVA 產生的問題不大。

假設 Y(ij) 是常態分布，則

D(ij) ≡ Y(ij) - Y(i.) ～ N(0, (1-1/n_i)σ^2_i)

諸 D(ij) 聯合服從多變量常態分布，且

Cov(D(ij), D(ij')) = (δ(jj') - 1/n_i)σ^2_i

其中 δ(jj') = 1 當 j = j'; 否則為 0。因此，

E[Z(ij)] = E[|D(ij)|] = √[(2/π)(1-1/n_i)σ^2_i]
Var[Z(ij)] = E[D(ij)^2] - (E[|D(ij)|])^2 = (1-2/π)(1-1/n_i)σ^2_i

上列兩式表示：

(1) 原本 Y(ij) 代表的 k 個群體間變異數相等與否的問題，幾乎轉換成 Z(ij) 所代表的群體之間平均數相等與否的問題。之所以只是「幾乎」而不是完全轉換，因為其中有 n_i 的作用，當各樣本大小 n_i 一致，或樣本均不小時，n_i 之差異造成的作用可以忽略，那麼 Z(ij) 之組間平均數差異，就代表了 Y(ij) 組間變異數的差異。

(2) 虛無假說不成立時，Z(ij) 所代表的群體之間變異數也是不等的。但檢定方法所依賴的是虛無分布，也就是虛無假說成立時檢定統計量的分布，而此時除各組樣本大小差異造成的影響外，變異數可以說是同幅的。

很明顯 Z(ij) 與 Z(ij') 不獨立，因為它們有共同的 Y(i.)；實際上可以證明 E[Z(ij)Z(ij')] ≠ E[Z(ij)]E[Z(ij')], j≠j'，不過其計算及結果太瑣碎，就不列了。注意 Z(ij) = |D(ij)|，而 Σ_j D(ij) = 0 for all i，若 n_i = 2, 則 Z(i1) = Z(i2); 但 n_i > 2 時並不能構成 Z(i1), ..., Z(in_i) 間的函數限制式，也就是說，由於取了絕對值，使得 Z(ij) 之間沒有了 D(ij) 間的線性相依現象，甚至也不存在完全函數關係。

在「線性模型：當誤差非等幅變異時」一文中，我們談到普通最小平方法和加權或一般最小平方法結果一致的條件

V = XΓX' + ZΘZ' + λI

其中 c(X) 和 c(Z) 正交。對 Z(ij) 做一因子變異數分析的問題，全模型的模型矩陣是

X = Blk Diag(J_1, . . ., J_k]

其中 J_i 是長度 n_i 元素值皆為 1 的行向量。而觀測值向量 Z = [ Z(ij) ] 的共變異矩陣正是

V = BlkDiag(V_1, . . ., V_k)
   = BlkDiag((1-ρ_1) I_1 + ρ_1 J_1 J'_1, . . ., (1-ρ_k) I_k + ρ_k J_k J'_k)
   = X Diag(τ_1, ..., τ_k)  X' + BlkDiag(λ_1  I_1, . . ., λ_k I_k)

的形式，其中 I_i 是階為 n_i 的單位矩陣。對此模型而言，各區塊之估計並不相干，也就是說

A = BlkDiak(A_1, . . ., A_k),

其中 A_i = J_i J'_i V_i^{-1}/J'_i V_i^{-1} J_i。但因

V_i = τ_i J_i J'_i +  λ_i I_i

因此普通最小平方法得到的是 BLUE。在縮減模型，也就是假設虛無假說成立，對應的模型矩陣是 J，一個元素皆為 1 的 n×1 矩陣或行向量，c(J) 上的垂直投影是 M* = JJ'/n。將 X 分解為

X = U + W,  其中  U = M* X,   W = (I - M*)X 

則 U = M* X = J(J'X/n)

V = τ XX' + λ I  = τ JΓJ' + τ WW' + λ I

其中 Γ = J'XX'J/n^2。故，普通最小平方法估計仍是 BLUE。當然如先前說過，可能 n_i 的不等會導致實際上即使 H0 成立仍存在 Z(ij) 對應的群體變異數（及平均數）略有差異，在此假設這個問題可以忽略。

在一般最小平方法的 ANOVA 中，檢定統計量是

F = [Y'(A-A*)'V^{-1}(A-A*)Y/(k-1)]／[Y'(I-A)'V^{-1}(I-A)Y/(n-k)]

其中 A 是在 c(X) 之投影而 A* 是在 c(J) 之斜交投影。此處因在全模型 (E[Z] in c(X) = c(A)) 與縮減模型 (E[Z] in c(J) = c(A*)) 中 Cov(Z) 不同，V 該對用哪個？由於決定拒絕 H0 與否是根據 H0 成立時的分布，因此

V = τ XX' + λ I;  所以  V^{-1} = α I + β XX'

而 AY = MY 且 A*Y= M*Y，結果

F = [Y'(M-M*)V^{-1}(M-M*)Y/(k-1)]／[Y'(I-M)V^{-1}(I-M)Y/(n-k)]
  = [Y'(M-M*)Y/(k-1)]／[Y'(I-M)Y/(n-k)]

總結來說，除了 Z(ij) 非常態性問題以外，普通 ANOVA 方法檢定組間平均數差異是適當的（條件是 n_i 都大於 2 並且 n_i 間差異不大，或 n_i 都夠大）。

當原資料 Y(ij) 是常態時，資料 Z(ij) 符合

Var(Z(ij)) = c(E[Z(ij)])^2

如前述，用變異數穩定化變換 (variance-stabilizing transformation) 可將資料轉換成

X(ij) = ㏒(Z(ij))

則 Z(ij) 之組間差異在 X(ij) 仍保持相同順序，但 X(ij) 將可望縮小組間變異數差距。當 Y(ij) 非常態時，Z(ij) 的資料變換及 ANOVA 的實施，基本上是認為 E[Z(ij)] 與 Var(Y(ij))，或更精確地說是與 √Var(Y(ij)) 成正向關聯，因此 Var(Z(ij)) 與 (E[Z(ij)])^2 的比例關係似可預期。換言之，前述定義的 Z(ij) 或許再取個對數變換，然後進行普通一因子變異數分析，或許更合適？不過，從 F 統計量來看，問題主要還是在於 Z(ij) 同組內的相關，而且這項相關即使對 Z(ij) 做了對數變換，其相關仍存在，但正如前面對 Z(ij) 做 ANOVA 時的討論，此項相關基本上不影響 ANOVA 的採行。

還有其他方法，例如 Jackknife 方法：

u(ij) = n_i ㏒(s^2_i) - (n_i - 1) ㏒(s^2_i(-j))

其中 s^2_i(-j) 是剔除 i 樣本中第 j 觀測值之後計算的第 i 樣本變異數，而 s^2_i 則是第 i 樣本之全樣本變異數。也是對 u(ij) 進行一因子變異數分析。