在一因子變異數分析中，當分組，也就是子群體間變異數不等時，如果我們已知各群體間變異數的相對比值，也就是「線性模型：當誤差非等幅變異時」一文假設的 Cov(ε) = σ^2 , 其中 V 已知，並且此處假設 V 是荏對角線矩陣 V = BlkDiak(v_1 I_1, . . ., v_k I_k), 此處 BlkDiag 是區塊對角線 (block diagonal matrix)，k 個區塊表示 k 個子群體 k 組子樣本，各子樣本內是常數變異數 σ^2 v_i, i = 1, ..., k。其估計和檢定在前引文及「線性模型：誤差項共變異非滿秩問題」一文已有敘述（含推導）。在常態誤差假設之下，前二文以向量幾何方式推導的線性模型方法，或所謂 Gauss-Markov 方法，與最大概似估計 (maximum likelihood estimate) 以及概似比檢定 (likelihood ratio test) 基本是一致的，唯一的差別是 σ^2 的估計，G-M 方法以自由度除 SSE 而 MLE 是以樣本數除之。本文是將前述一般模型具體化到一因子 ANOVA，也就是 k 個群體平均數比較的具體計算式表現出來。
一因子 ANOVA 模型可以寫成
Y_ij = μ_i + ε _ij, j = 1, ..., n_i; i = 1, ..., k
或
Y_ij = μ + α_i + ε _ij, j = 1, ..., n_i; i = 1, ..., k
其差別只在後者 μ 和諸 α_i 是不可估的，除非對諸 α_i 加上一個限制如 α_1 = 0 或 Σ α_i = 0，但 μ + α_i 也就等同於前一模型的 μ_i 是可估的。在 ANOVA 中我們重點是諸 μ_i 是否相等之檢定，因此兩個模型是一樣的；以線性模型的術語來說，我們所要檢定的假說是一些關於可估函數 (estimable functions) 的假說；而兩模型的設計（或模型）矩弭 (design or model matrix) X 雖不同，但有相同的行空間 (column space)，因此兩模型無實質差別。為了方便，我們取前一個模型，其設計矩陣為
[ J_1 ]
[ J_2 ]
X = [ ． ]
[ ． ]
[ J_k ]
這個 X 是行滿秩 (full column rank) 的，使得我們以下符號有廣義反矩陣 (generalized inverse matrix) 的符號處，其實就是真正的反矩陣。式中 J_i 表階為 n_i × 1 元素皆為 1 的行向量。模型中誤差項 ε_ij 的變異數是 σ^i v_i，因此整個誤差向量 ε 的共變異矩陣是 σ^2 V, 其中
[ v_1 I_1 ]
[ v_2 I_2 ]
V = [ ． ]
[ ． ]
[ v_k I_k ]
式中 I_j 是 n_j 階單位矩陣，所以上列 V 雖以區塊對角線矩陣表示，但它實實地就是一個對角線矩陣。我們知道此模型，或說 c(X) 上，的投影矩陣 (projection matrix) 是
A = X[(X'V^{-1}X)^-]X'V^{-1}
這裡依筆者個人習慣，M' 表矩陣 M 的轉置 (transpose)；因為不涉及複數 (complex number)，因此 M* 將是用來表示一個與 M 不同的矩陣，但不一定是 M 經運算而得。另外，上標 M^{-1} 代表 M 的反矩陣，而 M^- 代表 M 的廣義反矩陣。首先計算 X'V^{-1}
[ (J_1)'/v_1 ]
[ (J_2)'/v_2 ]
X' V^{-1} = [ ． ]
[ ． ]
[ (J_k)'/v_k ]
其次是 X' V^{-1} X
[ n_1/v_1 ]
[ n_2/v_2 ]
X' V^{-1} X = [ ． ]
[ ． ]
[ n_k/v_k ]
這是一個 k × k 對角線矩陣，故其（廣義）反矩陣為主對角線取倒數：
(X' V^{-1} X)^- = diag( v_1/n_1, . . ., v_k/n_k )
最後得投影矩陣
A = X[(X'V^{-1}X)^-]X'V^{-1}
[ J(11)/n_1 ]
[ J(22)/n_2 ]
= [ ． ]
[ ． ]
[ J(kk)/n_k ]
式中 J(ij) 為 n_i × n_j 元素皆為 1 的矩陣。因此，平均值向量 Xβ 的估計為
AY = X[(X'V^{-1}X)^-]X'V^{-1}Y
= [ Ybar_1(J_1)' . . . Ybar_k(J_k)' ]'
或直接用行向量式表現：
[ J_1 Ybar_1 ]
[ ． ]
AY = [ ． ]
[ ． ]
[ J_k Ybar_k ]
故得模型平方和
Y'A' V^{-1} AY = Σ n_i(Ybar_i)^2/v_i
檢定 H0: μ_1 = . . . = μ_k 對抗 Ha: 諸 μ_i 任意，等於全模型 E[Y] = Xβ 與縮減模型 H0: E[Y] = X*γ 的比較，
X* = [ (J_1)' ． ． ． (J_k)' ]’
實際上是 n × 1 行向量，n = Σ n_I，故 X*' V^{-1} X* 是一個純量
X*' V^{-1} X* = Σ n_i/v_i
而此縮減模型或 c(X*) 上的投影矩陣為
A* = X*[(X*' V^{-1} X*)^-]X*' V^{-1}
[ J(11)/v_1 ．．． J(1k)/v_k ]
[ ． ． ]
= 1/Σ(n_i/v_i) [ ． ． ]
[ ． ． ]
[ J(k1)/v_1 ．．． J(kk)/v_k ]
結果，Y 在 c(X*) 的投影是
A*Y = Σ(n_i Ybar_i/v_i)/Σ(n_i/v_i) [ (J_1)' ．．．(J_k)' ]'
其平方和為
Y'A*' V^{-1} A*Y = Σ(n_i/v_i) [Σ(n_i Ybar_i/v_i)/Σ(n_i/v_i)]^2
兩模型平方和的差就等於兩模型平均向量 E[Y] 估計之差的平方和 (SSR)
Y'(A-A*)' V^{-1} (A-A*)Y = Y'A' V^{-1} AY - Y'A*' V^{-1} A*Y
而離差平方和等於總平方和減去模型平方和 (SSE)
Y'(I-A)Y = Σ[Σ(Y_ij - Ybar_i)^2]/v_i = Y' V^{-1} Y - Y'A' V^{-1} AY
檢定是用
F = [(SSR - SSR*)/(k-1)]/[SSE/(n-k)]
其中
SSR = Σ n_i(Ybar_i)^2/v_i
SSR* = Σ(n_i/v_i) (Ybar)^2
Ybar = Σ(n_i Ybar_i/v_i)/Σ(n_i/v_i)
SSE = ΣΣ(Y_ij-Ybar_i)^2/v_i
除了分母多個 v_i 以外，和同幅變異假設下的 ANOVA 並無不同。假設 k = 2，則
SSR - SSR* = (Ybar_1 - Ybar_2)^2/(v_1/n_1 + v_2/n_2)
= [(Ybar_1 - Ybar_2)/√(v_1/n_1 + v_2/n_2)]^2
而 Var(Ybar_1 - Ybar_2) = σ^2 (v_1/n_1 + v_2/n_2)，其中 σ^2 是用 SSE/(n-k) 估計，故前述 F 檢定統計量等於 t^2，此處 t 是自由度 n-k 的 t 統計量。
當 Var(Y_ij) 是 σ_i^2 未知，而不是 σ^2 v_i 僅一個共同的 σ^2 未知時，k = 2 的情形有一種方法是採
t = (Ybar_1 - Ybar_2)/√(s_1^2/n_1 + s_2^2/n_2)
這時 t 不是 Z/√(χ^2/r) 形式，因為
(S_1^2/n_1 + S_2^2/n_2)/(σ_1^2/n_1 + σ_2^2/n_2)
並不服從「卡方變量除以其自由度」的分布。假如我們把
S_1^2/n_1 + S_2^2/n_2 近似 c χ^2(r)
則其期望值等於 rc，變異數等於 2rc^2，由
E[S_i^2] = σ_i^2, Var[S_i^2] = 2σ_i^4/(n_i-1)
可得
[σ_1^2/n_1 + σ_2^2/n_2]^2
r' = --------------------------------------------------------
[(σ_1^2/n_1)^2/(n_1-1) + (σ_2^2/n_2)^2/(n_2-1)]
因 σ_i^2 未知，以樣本值 s_i^2 代之，此謂之 Satterthwate 近似，而以此 t 統計量配合 t(r) 數值表進行兩常態群體平均數差異之檢定，稱為 Welch t-檢定。舊的統計學教本建意先做兩群體變異數相等檢定，再視結果採行聯合 σ^2 估計之 t-檢定或 Welch t-檢定，但一些研究發現兩階段的檢定沒有必要，直接做 Welch t-檢定反而較好。我們看上述 t 統計量分母根號內的內容，它是 s_1^2/σ_1^2 和 s_2^2/σ_2^2 分別以權量 σ_1^2/n_1 和 σ_2^2/n_2 加權的平均，是兩個 χ^2/(df) 的加權平均，以單一的 χ^2/(df) 近似貌似合理；當樣本足夠大時，s_i^2 近似 σ_i^2 使 t 近似服從常態分布，即使群體不是常態也可能有不錯的近似。這現象不論 σ_i^2 與 σ_2^2 是否相等其實都成立的。反觀假設 σ_1^2 = σ_2^2 的聯合變異數估計的 t 統計量，基本上是樣本平均數差之變異數的有偏估計：
E[{(n_1-1)S_1^2 + (n_2-1)S_2^2}/(n_1+n_2-2)](1/n_1 + 1/n_2)
= {[(n_1-1)σ_1^2 + (n_2-1)σ_2^2]/(n_1+n_2-2)}(1/n_1 + 1/n_2)
≠ σ_1^2/n_1 + σ_2^2/n_2
除非 σ_1^2 = σ_2^2。而先做變異數是否相等的檢定，允許在兩群體變異數不等時用錯誤的 t 檢定進行平均數差異之檢定。因此，除非兩群體變異數「真的」一致，否則不應該用傳統的兩樣本 t 檢定，先做兩群體變異數是否相等的檢定再選擇檢定方法，也是沒必要的，不如直接做 Welch t-檢定。
在 ANOVA 問題，全模型的模型平方和和誤差平方和分別是
SSR = Σ n_i(Ybar_i)^2/s_i^2
SSE = ΣΣ(Y_ij-Ybar_i)^2/s_i^2
而縮減模型的模型平方和是
SSR* = Σ(n_i/s_i^2) (Ybar)^2
Ybar = Σ(n_i Ybar_i/s_i^2)/Σ(n_i/s_i^2)
由前面可知檢定 μ_i 間完全相等所要的組間平方和是
SSA = SSR - SSR* = Σ n_i(Ybar_i - Ybar)^2/s_i^2
Welch (1951) 基本上沒有用到 SSE，只用 SSA 做檢定，其方法是:
F = [SSA/(k-1)]/(1+B)
其中 B = [2(k-2)/(k^2-1)] Σ[1-(n_i/s_i^2)/Σ(n_j/s_j^2)]/(n_i-1)，與 F(ν_1, ν_2) 比較，
ν_1 = k - 1
ν_2 = 1/{[3/(k^2-1)] Σ[1-(n_i/s_i^2)/Σ(n_j/s_j^2)]/(n_i-1)}
首先，在 H0: μ_1 = . . . = μ_k = μ 成立之下，
SSA = SSR - SSR* = Σ n_i(Ybar_i - Ybar)^2/s_i^2
= Σ n_i(Ybar_i - μ)^2/s_i^2 - (Ybar - μ)^2Σn_i/s_i^2
第一部分是 k 個相互獨立的 F(1,n_i-1) 變量的和，其期望值為 Σ[1+2/(n_i-3)]；第二部分
(Ybar - μ)^2Σn_i/s_i^2 = [Σn_i(Ybar_i-μ)/s_i^2]^2/(Σn_i/s_i^2)
比較複雜，其期望值為
E[(Ybar - μ)^2Σn_i/s_i^2] = E[Σ(n_i/s_i^2)(σ_i^2/s_i^2)/Σ(n_i/s_i^2)]
也有點麻煩。但總體來說，以 F 近似 SSA 經適當調整的分布貌似合理，至於再以 1+B 除，並且 F 的分母自由度的計算公式，這裡就不討論了。現在的統計軟體大概都會有 Welch 方法的選項。另外，我們注意到上述 SSR 或 SSA 已經用到 s_i^2 訊息了，因此不再需要 SSE，畢竟它只是諸 s_i^2 的綜合。
以概度比檢定而言，全模型下
(μ_i)^ = Ybar_i = Σ_j Y_ij/n_i
(σ_i^2)^ = s_i^2 = Σ_j (Y_ij - Ybar_i)^2
而縮減模型下
μ^ = Ybar = Σ(n_i/s'_i^2)Ybar_i/Σ(n_i/s'_i^2)
(σ_i^2)^ = s'_i^2 = Σ_j (Y_ij - Ybar)^2
概度比統計量是
λ = Π_i (s_i^2/s;_i^2)^{n_i/2}
= Π[Σ_j (Y_ij - Ybar_i)^2/Σ_j (Y_ij - Ybar)^2]^{n_i/2}
取 T = -2 ㏑(λ)，則
T = Σ n_i ㏑[1 + (Ybar_i - Ybar)^2/s_i^2]
或者一個不同但看似有正向關係的平方和
T* = Σ n_i (Ybar_i - Ybar)^2/s_i^2 = SSA
因為 ㏑(1+x) 是 x 的嚴格增函數，所以把 T 換成 T* 看似合理，但它們並不等價；而且除非 (Ybar_i - Ybar)^2 相對於 s_i^2 甚低，否則 T 並不能用 T* 近似。不過，T* 正好等於 SSA，也給予 Welch 法一個支持。
上述 SSR 等「平方和」因為已涉及諸 s_i^2，不是常態變量的二次式，若忽略此一因素，是否 SSA 和 SSE 也能用卡方近似（Satterthwaite 近似）？Satterthwaite (1942) 有另一種近似方法，不過筆者還沒弄懂；他在 1941 和 1946 有兩篇論文大概也是這方面的研究。 Julia Volaufova 不知什麼時候發表的，一份異幅變異之下做 ANVA 的論文清單。
留言列表
統計 (8)
{{ article.title }}