本文是一般線性模型的書籍主要談的模型，算是補足「線性模型：誤差項共變異非滿秩問題」。為什麼前引文只談一般少見人談的誤差項共變異矩陣 σ^2 V 非滿秩的情況？因為更早曾談過一般線性模型，其中 V 假設是滿秩的 (full rank)。不過，該文未特別討論 V = I 的特殊情況，即最基本情況；而且重點是模型的參數估計，或所謂模型配適 (model fitting, 今多譯為「模型校估」, 此處譯「配適」算是先入為主及個人偏好)，未考慮參數檢定或模型比較問題。
線性模型 (Linear models) 指（矩陣表示）：
Y = Xβ + ε, E[ε] = 0, Cov(ε) = σ^2 V, V 已知
式中 Y 是 n×1 資料向量，X 是 n×k 模型矩陣或設計矩陣，β 是 k×1 未知的模型參數向量，ε 是 n×1 的不可觀測誤差向量。
此線性模型可分 4 個類別（或： 4 個等級）：
(1) V = I; (1') V = I;
(2) V 可逆； (2') V ≠ I 但可逆；
(3) C(X) 是 C(V) 的子空間； (3') V 不可逆，C(X) 是 C(V) 的子空間；
(4) 任意 V。 (4') V 不可逆，C(X) 非 C(V) 的子空間。
上面左邊的列法是等級式的列法；而右邊等於將左邊各等級排除前一等級，是分類式的列法。本文只談 (1), (2) 兩種情形；「線性模型：誤差項共變異非滿秩問題」談的是 (3), (4) 兩種情形。
模型 Y = Xβ + ε 用普通代數式表示，可以表示為
Y_i = β_0 X_{0i} + β_1 X_{1i} + . . . + β_p X_{pi} + ε_i, i = 1, ..., n
此處 p = k-1; X_{0i} 常見的是 X_{0i} ≡ 1，所謂「常數項」, 不一定要有。資料項 Y_i 是反應變數 (response variable) 或依變數 (dependent variable) Y 的第 i 個觀測值，是伴隨著 (X_{0i}, X_{1i}, ..., X_{pi}) 被觀測的，X_{ji} 是變數 X_j 的第 i 「觀測值」。(X_{0i}, X_{1i}, ..., X_{pi}) 可能是事先被設計的，如設計好的統計實驗，事先定好諸 X_{ji}，而後在設計的 X_{ji} 值組合之下觀測 Y_i，這就是矩陣 X 被稱為設計矩陣的由來。當然諸 X_{ji} 也可能完全伴隨 Y_i 被同時觀測到的（此處「同時」不一定是時間上的相同，可以先有 X_{ji} 的出現，後來才出現 Y_i），如調查資料。諸變數 X_j 稱解釋變數 (explanatory variable)、預測變數 (predictor variable) 或自變數 (independent variable)，如前述它們可以是預先設計（定值）的或被與反應變數同時觀測的，所以矩陣 X 稱模型矩陣可能更適當，因為固然「觀測」也是在設計之下，但設計一詞卻隱有預先定值的意思。
線性模型最基本的有線性迴歸模型、實驗設計模型及其他線性模型。線性迴歸模型通常有常數項，即 X_{0i} ≡ 0；但也非必然，在某些特定問題，無常數項的迴歸模型可能更適當。當然，迴歸模型的常數項可能只是個調整項，特別是當我們用來分析資料的模型只是一個局部近似而非全域模型時，即使理論上不應有常數項，也可能加個常數項使得模型可以更切合資料。迴歸模型的 X 矩陣常是行滿秩 (full column rank)，這使得 X'X 可逆，或 X'V^{-1}X 可逆，若 X'X 不可逆，或只是數字上的近乎不可逆，就產生「線性重合問題」或稱多重共線性 (multicollinearity) 問題，使得模型中的迴歸係數參數向量 β 不可估計（完全線性重合）或誤差膨脹而使得 β 的估計或計算不穩定（高度線性重合）。線性重合問題及其解決可詳參線性迴歸專書，此處不討論。完全線性重合意謂有多餘的解釋變數項，例如 X_3 是 X_1 與 X_2 的線性組合，則有一個是多餘的變數；從模型參數來看就是過度參數化 (over-parameterization), 模型中有超過必要的參數，例如變異數分析模型
Y_{ji} = μ_j + ε_{ji}, j = 1, ..., k; i = 1, ..., n_j
這是以組平均數 μ_j 為參數的平均數模型‵ (cell means model)，其設計矩陣 X 是線性獨立的 k 行；但在 ANOVA 中常用因子效應模型 (factor effects model)
Y_{ji} = μ + α_j + ε_{ji}, j = 1, ..., k; i = 1, ..., n_j
表示，其 X 矩陣有 k+1 行，第一行是常數 1，是後面 k 行的和，因此 μ 及諸 α_j 不是唯一定義的，或說是不可辨認的 (unidentifiable)，在線性模型它們是不可估計的 (non-estimable)。在 ANOVA 中，通常會給予一個不可估的限制式如
Σ_j α_j = 0 或 Σ_j n_j α_j = 0
使得參數 μ 與諸 α_j 可以估計，但不同限制條件使得 μ 與諸 α_j 的估計結果不同。不可估限制式使模型參數可估，其實就是藉由限制式縮減模型中的多餘參數，消除過度參數化；若限制條件太多，則相當於對一個適度參數化的模型附加限制，用可估函數建立限制條件也一樣，造成模型的縮減。雖然不同的不可估限制條件會產生不同參數估計，但其實這些不同估計結果就是對應了不同限制條件下的參數真值，以上述 ANOVA 的因子效應模型為例，因為多出一個常數項，使得參數過多個別參數不可估計。若加了 Σ_j n_j α_j = 0 的限制，令
Y_j. = Σ_i Y_{ji}/n_j, Y.. = (Σ_j n_j Y_j.)/(Σ_j n_j)
則 Y.. 是 μ 的不偏估計而 Y_j. - Y.. 是 α_j 的不偏估計。而若限制式改成 Σ_j α_j = 0，令 Y'.. = Σ_j Y_j./k，則 Y'.. 是 μ 的不偏估計而 Y_j. - Y'.. 是 α_j 的不偏估計。注意兩種估計不同，其期望值也不同，分別對應到在限制條件下的參數值。例如在 Σ_j α_j = 0 限制下，
μ = Σ_j (μ + α_j)/k = E[ Σ_j Y_j./k ] = E[ Y'.. ]
注意在上述兩種情形，μ 的不偏估計加上 α_j 的不偏估計都是 Y_j.，其期望值都是 μ + α_j，這結果在無上述附加條件限制式時也成立，像這種情形，參數的線型函數式 μ + α_j 具有不偏估計 Y_j.，稱此種參數函數為可估的 (estimable)。
回到一般模型，第一級的線性模型 V = I 最基本的，我們考慮普通最小平方法，極小化
Q(β) = ||Y - Xβ||^2 = (Y-Xβ)'(Y-Xβ)
用微積分的方法，得「標準方程式」 X'Xβ = X'Y，稱之為「標準方程式 (normal equations)」,或譯「法線方程式」，因該方程式亦可表示為
X'(Y - Xβ) = 0
其中符合方程式的 β，以 β^ 表示，是參數 β 的最小平方估計式 (LSE, least squared estimator)。而 X β^ 即為 E[Y] = μ 向量的不偏估計，是「配適值 (fitted value)」向量，Y - X β^ 是反應向量與配適值向量的差，稱為殘差 (residuals) 向量。法線方程式或標準方程式的意思是：最小平方估計使殘差向量與 X 矩陣的各行向量垂直，或在線性代數中稱之為正交 (orthogonal)。
標準方程式除了如上述可以用微分法得到以外，也可以由配方法得到：
Q(β) = (Y-Xβ)'(Y-Xβ)
= β'X'Xβ - Y'Xβ - β'X'Y + Y'Y
= (β-b)'X'X(β-b) + b'X'Xβ + β'X'Xb - b'X'Xb
- Y'Xβ - β'X'Y + Y'Y
= (β-b)'X'X(β-b) - 2β'X'(Y-Xb) + (Y'Y - b'X'Xb)
如果 b 滿足 X'(Y-Xb) = 0，則 β = b 時 Q(β) 達最小值；不過，如果 X'X 不是可逆的，其實只需要 X(β-b) = 0，或即 Xβ = Xb，其中 b 滿足標準方程式。
另一種方式，以幾何方法來看，不過這裡的「幾何」不是中學平面和三度空間的幾何，而是 n 度空間的幾何。反應資料 Y 是 R^n 的一個點或一個 R^n 向量，依我們把 R^n 看成是 n 度歐氏點集或 n 度向量空間而定。E[Y] = Xβ 的意思是 μ 在 X 的行空間中，那麼對 μ 的估計也應在 X 的行空間，μ^ = Xb，最小平方法就是要極小化殘差向量 Y - Xb 的長度，也就是 Y 到 Xb 的距離。如同平面上一點到直線的距離，空間中一點到一直線或一平面的距離，都以垂直距離為最短。現在 Y 在 n 度空間 R^n，與 C(X) 即 X 的行空間之最短距離也應是垂直距離，因為在 R^n 中的距離，除了維度可能較多以外，和 R^2（平面）上一點到一直線或 R^3（空間）中一點到一平面或一直線的情形是一樣的，所以垂足 Xb 滿足
Y - Xb ⊥ C(X), 等價地， X'(Y-Xb) = 0
跟先前得到的結果一致：b 是 β 的最負平方解，若且唯若標準方程式 X'(Y-Xb) = 0。
對於最小平方法的標準方程式 X'Xβ = X'Y，如果 X'X 是可逆的，則得解
β^ = (X'X)^{-1} X'Y
但如果 X'X 不可逆呢？也就是 X 不是行滿秩的情形，前面說過此時模型是過度參數化的，β 是不可估的。
［定理］線型函數 λ'β 可估若且唯若存在 ρ 使 λ = X'ρ.
［證］ λ'β 可估意謂存在實數 a 及向量 ρ 使 E[a+ρ'Y] = λ'β 對所有 β 成立。
故若 λ'β 可估，則上列 a, ρ' 存在，使 a + ρ'Xβ = λ'β，則 a = 0 且 ρ'X = λ'。
反之，若存在 ρ 使 λ = X'ρ 則 E[ρ'Y] = ρ'Xβ = λ'β。
上面是考慮實數值線型函數式，但把 λ'β 換成向量值線型函數 Λ'β，相似的推論可知其為可估，充分且必要條件是存在矩陣 P 使 Λ = X'P。一個特例是 Xβ 可估，因為 Y 本身就是其不偏估計。但是，前面說過，既然 Xβ 是在 C(X) 中，我們要找其合理的估計也要從 Xb 中去找；而如果要找的估計式是 Y 的線型函數 PY（雖然應找 c + PY 形的估計，但顯然 c 應取 0 向量。） 所以
PY = Xb; 且 PXβ = Xβ 對所有 β 成立。
上式等於要求 PY = Xb 且 PX = X，也就是說：矩陣 P 視為 R^n 上的一個線性變換，它把 R^n 的所有元素都映至 C(X)；並且所有 C(X) 上的元素，經 P 的變換將維持原樣。這樣的線性變換稱為 C(X) 上的投影算子 (projection operator) 或只稱投影，矩陣 P 就是投影矩陣。如果不關心投影的目的子空間，單論一個線性變換是投影算子的充要條件是 P^2 = P，或說 P 是自乘不變 (idempotent) 矩陣。所以，對於 Xβ 我們要的合理線性不偏估計是：
找一個 C(X) 上的投影矩陣 P，以 PY 估計 Xβ.
但最小平方法估計，由幾何上來看，是：
Xβ 的最小平方估計是 Y 在 C(X) 上的垂直投影。
資料向量在 C(X) 上的垂直投影是 X β^，則殘差向量 Y - X β^ 其實是 Y 在 C(X) 的正交補空間 (orthogonal complement) 上的垂直投影。以 M 代表 C(X) 上的垂直投影，則在 C(X)^⊥ 上的垂直投影算子是 I-M，故
對任意 x, y in R^n，(I-M)x ⊥ My
相當於 (I-M)'M = 0，或說 M = M'M，因此，
M 是垂直投影矩陣若且唯若 M 是對稱的而且 M^2 = M。
因為 M 是 C(X) 上的垂直投影矩陣，所以顯然 M 必須是 XAX' 形式，並且
XAX'XAX' = XAX', XAX'X = X
如果選擇 A 使 X'X 等於 A 的廣義反矩陣 (generalized inverse) A^-，則由廣義反矩陣之定義得 AX'XA = A，第一個式子成立。再者，如果第二式成立，則第一式也成立。前面已知：當 X'X 可逆時，Xβ 的 LSE 是 X(X'X)^{-1}X'Y，即 M = X(X'X)^{-1}X'；在 X'X 不可逆時，是否可以把 (X'X)^{-1} 換成廣義反矩陣，也就是取 A = (X'X)^-？
［定理］M = X (X'X)^- X' 是 C(X) 上的垂直投影。
［證］ 顯然 M 的值域在 C(X) 內，因此需要證明 M 對 C(X) 上的元素不會改變，即 MX = X。
設 y = u + v 是 R^n 中的任一元素，分解成 u = Xb 在 C(X) 中，v 在 C(X)^⊥ 中。則
y'MX = (b'X' + v') X (X'X)^- X'X
= b'X'X (X'X)^- X'X
= b'X'X
= (b'X' + v') X = y'X
因 y 任意，所以 MX = X，所以 M^2 = M。
設 A 是對稱矩陣，G 是 A 的任一廣義反矩陣，則
AGA = A = A' = (AGA)' = AG'A
這表示：(1) 存在對稱的 A^-，(2) A A^- A 是對稱的。
所以 M' = M。
當 X'X 可逆時，其反矩陣是唯一的廣義反矩陣；當 X'X 不可逆時，有無限多廣義反矩陣。但 C(X) 上的垂直投影只有唯一一個，這表示雖然 (X'X)^- 不唯一，但所有 M = X (X'X)^- X' 都是同一個投影算子。結論是：不論 X 是否行滿秩，Xβ 的 LSE 都是
(Xβ)^ = MY = [X (X'X)^- X'] Y
當 X 行滿秩時 (Xβ)^ = X β^ = X (X'X)^{-1} X'Y。
最小平方估計是「資料配適 (data fitting)」的估計方法，可以不顧 Y 或 ε 的機率分布特性；但在估計問題中，我們更關心估計方式或估計結果的評估，以決策的觀點，就是風險 (risk) 的評估。在線性模型，最常用的兩種評估：一是 Xβ^ 的共變異矩陣，二是可估實值線性函數 λ'β 之線性不偏估計 (linear unbiased estimator) 的變異數。
［定理］可估函數 λ'β 之 LSE 是 LUE 中變異數最小的；
可估函數 Λ'β 之 LSE 是 LUE 中共變異矩陣最小的。
［證］ 首先看實數值 λ'β 的估計問題。
設 λ = X'ρ, 則 λ'β 的 LSE 是 ρ'MY；而其 LUE 為 (Mρ+d)'Y，則 d'X = 0。
因假設 Cov[Y] = Cov[ε] = σ^2 I，由於 d'X = 0, 故 Md = 0, 故
Var[(Mρ+d)'Y] = σ^2 (Mρ+d)'(Mρ+d) = σ^2 (ρ'Mρ+d'd)
其最小值發生若且唯若 d = 0.
設 Λ = X'P, 則 Λ'β 的 LSE 是 P'MY；而其 LUE 為 (MP+D)'Y，則 D'X = 0。
類似賁數值可測函數情形，
Cov[(MP+D)'Y] = σ^2 (P'MP+D'D)
因為 D'D 是正確定 (positive definite) 或正半定 (positive semi-definite) 矩陣，除非 D = 0，再一次證明了 LSE 是 LUE 中共變異矩陣最小的。
在所有線性且不偏估計式中變異數或共變異矩陣最小的估計式，稱為「最佳線性不偏估計 (Best Linear Unbiased Estimator, BLUE)」。上列定理稱 Gauss-Markov Theorem，指出最小平方估計同時也是 BLUE，雖然只限制在線性且不偏的估計量裡面的最佳，也算是一種最佳性。如果誤差項的分布是 N(0, σ^2 I)，則最佳的範圍可擴充至所有不偏估計，而且最小平方估計就是 MLE；但如果比較的範圍不限於不偏估計，又對向量值參數函數估計量評價標準不是共變異或均方誤矩陣而是總和均方誤，甚至會發生 Stein 現象。
統計推論無非估計和檢定，可估函數 Λ'β 之 LSE 是 BLUE，而其估計誤差，以共變異矩陣來看是
Cov(Λ'β^) = Cov(P'MY) = σ^2 PMP = σ^2 Λ' (X'X)^- Λ
其中 σ^2 可用‵ MSE = Y'(I-M)Y/rank(I-M) 來估計，因為
E[Y'(I-M)Y] = σ^2 tr(I-M) = σ^2 rank(I-M) = σ^2(n-rank(M))
如果要做 H0: Λ'β = 0 的檢定，可以用
W = β'^Λ [Λ' (X'X)^- Λ]^- Λ'β^/MSE
= [Y'MP (PMP)^- P'MY]/MSE
另一方面，由 Λ'β = 0 是模型的縮減：
μ = E[X] = Xβ and Λ'β = P'Xβ = 0 <==> μ = Zγ
或者說：
μ in C(X) and μ ⊥ C(P) <==> u in C(Z)
首先，μ = Xβ ⊥ C(P) 若且唯若 μ ⊥ C(MP)，MP 是把 P 垂直投影到 C(X) = C(M)。也就是說，Λ'β = 0 對應到 C(X) 中與 C(MP) 垂直的子空間。所以，一個想法是取
Ｚ = (I-M_{MP})X；
另一個想法是取
Z = M - M_{MP}
前一個 Z 很容易證得
C(Z) = C(X)∩C(P)^⊥ = C(X)∩C(MP)^⊥
後一個想法可以直接由
M = M_{MP} ⊕ (M - M_{MP})
而 M, M_{MP} 及兩者之差都是垂直投影矩陣，所以
M_Z = M - M_{MP}
而假說 H0: Y = Zγ + ε 對假說 Ha: Y = Xβ + ε 的檢定可以用模型殘差平方和增量對大模型殘差平方和之比來達成：
F = {[Y'(I-M_Z)Y -Y'(I-M)Y]/(rank(I-M_Z)-rank(I-M))}/{Y'(I-M)Y/rank(I-M)}
= {Y'M_{MP}Y/rank(M_{MP})}/{Y'(I-M)Y/rank(I-M)}
此式分母即是 MSE，而 F 和前面的 W 僅差在 W 的分子沒有用 C(MP) 的維度，統計上所謂的自由度 (degrees of freedom)，調整。前面的 W 是用在大樣本，中央極限定理適用的情況，W 漸近服從卡方。而 F 適用於模型誤差項為常態，ε~N(0, σ^2 I)，小樣本情形。
第二類或第二級線性模型是 Cov(ε) = σ^2 V, 其中 V 已知，可逆，可參考「線性模型：當誤差非等幅變異時」一文。
設 V = QQ'，因 V 可逆，此處 Q 也是 n×n 可逆方陣。對原模型做下列變數變換：
Q^{-1}Y = Q^{-1}X β + Q^{-1}ε
則把問題化成第一類問題，所以普通最小平方法用於轉換後模型，得 Q^{-1}X β 的 LSE 為
Q^{-1}X β^ = (Q^{-1}X) [(Q^{-1}X)'(Q^{-1}X)]^- (Q^{-1}X)'Q^{-1}Y
或者，
X β^ = X (X'V^{-1}Ｘ)^- X'V^{-1}Y
從 LSE 來看，轉換模型的最小平方法是極小化
Q(β) = (Q^{-1}Y-Q^{-1}Xβ)'(Q^{-1}Y-Q^{-1}Xβ)
= (Y-Xβ)'V^{-1}(Y-Xβ)
這就是對原模型以 V^{-1} 加權的加權最小平方準則。而以 n 維空間幾何來看，Q^{-1}Y 在 C(Q^{-1}X) 上的正交投影，當逆轉回 Y 在 C(X) 之「投影」時是
Proj_{C(X)} Y ≡ Q (Proj_{C(Q^{-1}X)} Q^{-1}Y)
意思是將 Proj_{C(X)} Y 在 C(X) 的投影定義為 Q^{-1}Y 在 C(Q^{-1}X) 上的正交投影做逆轉換 Q。於是，此投影的矩陣是
A = Q M_{Q^{-1}X} Q^{-1}
= Q (Q^{-1}X) [(Q^{-1}X)'(Q^{-1}X)]^- (Q^{-1}X)' Q^{-1}
= X [X'V^{-1}X]^- X'V^{-1}
易證 A^2 = A, 所以 A 是一倨投影矩陣，但它卻不像 M 那樣是一個垂直投影——垂直投影矩陣是一個對稱矩陣，因為投影 A 和「殘差」垂直的充要條件是
A'(I-A) = 0 = A' - A'A
當 A 對稱時，A'(I-A) = 0 成立；反之，當 A'(I-A) = 0 時，因為 0 和 A'A 都是對稱方陣，所以 A 對稱。一般而言前面那投影矩陣 A 是不對稱的，可稱之為「斜交投影 (oblique projection)」, 或譯斜投影、斜軸投影。
與 V 相關的投影 A 矩陣很容易得 AV = VA' 及 A'V^{-1} = V^{-1}A；再者，
X' A'V^{-1}A X = X' A'V^{-1} X = X' V^{-1}A X = X' V^{-1} X
由此可得 X'(I-A)'V^{-1}(I-A)X = 0，因此證得 AX = X。當然也可由 V = QQ' 及經變數變換的模型證得 AX = X 的結果。
一個線性函數是否可估只涉及期望值，即：存在 P'Y 使 E[P'Y] = Λ'β 對任意 β in R^k 都成立，則線性函數 Λ'β 可估。但 E[P'Y] = P'Xβ，因此 Λ'β 可估的充要條件是存在 P 使 Λ' = P'X。這條件在各級線性模型都成立。當 Λ'β 可估，（加權）最小平方估計 Λ'β^ = P'AY 是 BLUE:
Cov(B'Y) = E[(B'Y-P'AY+P'AY-Λ'β)(B'Y-P'AY+P'AY-Λ'β)']
= E[(B'Y-P'AY)(B'Y-P'AY)'] + E[(P'AY-Λ'β)(P'AY-Λ'β)']
+ E[(B'Y-P'AY)(P'AY-Λ'β)'] + E[(P'AY-Λ'β)(B'Y-P'AY)']
因 B'Y 的不偏性，可證得 Cov((B'-P'A)Y, P'AY) = 0，故 Cov(B'Y) 等於 Cov(P'AY) 加一非負確定項。後者等於
Cov(P'AY) = σ^2 P'AVA' P = σ^2 Λ' (X'V^{-1}X)^- Λ
而 σ^2 可用
MSE = Y'(I-A)'V^(-1)(I-A)Y/tr(I-A)
估計。如欲做 Λ'β = 0 之檢定，則取
F = {Y'A'P (P'AVA' P)^- P'AY/rank(P'A)}/MSE
= {[β'^Λ [Λ' (X'V^{-1}X)^- Λ]^- Λ'β^]/rank(Λ)}/MSE
當 V = I 時，所有結果與先前一致。
留言列表
統計 (9)
{{ article.title }}