6.5 一般线性模型的统计推断

本节将脱离实际背景, 从理论上来讨论一般线性模型的参数统计推断. 考虑模型依然为 y=Xβ+ε, 这里 X∈Rn×p, rankX≤p, n≥p. β∈Rp 是未知参数向量. 误差向量 ε:Eε=0.

1 可估参数函数及其估计

我们从 β 的最小二乘估计开始. 在 这里 我们得到它是正规方程 (1.1)XTXβ^=XTy 的解. 现在解释它的直观几何意义. 要使 ||y−Xβ||
最小, 则 Xβ 应该是 y 在平面 μ(X) 上的投影: Xβ^=PXy, 如图.
Pasted image 20260104152039.png|200
但是如果 rankX<p, 正规方程 (1.1) 有无穷多的解, 此时 β 是不可估的. 为此, 引入如下定义:

可估函数

设 aTβ 是参数 β 的线性函数. 如果有 y 的线性函数 cTy: E(cTy)=aTβ, 则 aTβ 是 β 的可估函数; 否则 aTβ 是不可估的.

定理 1.1

在 (1.1) 中, aTβ 可估等价于 a∈μ(XT).

推论

  1. Xβ 的每个分量 eiTXβ 都是可估的.
  2. 如果 rankX=p, 则 ∀a∈Rp, aTβ 都是可估的; 如果 rankX<p, 则 β 至少有一个分量不可估.

不可估还可以用正规方程的多解性说明. 此时 XTXβ=0 有非零解 β0. 如果 β∗=β+β0, 则必有 Xβ∗=Xβ, 而 β 在模型中的作用是通过 Xβ 体现的, 因此无法从模型推断 β, β∗.

如果 aTβ 可估, 记它线性无偏估计的全体为 EU(a)={cTy|E(cTy)=aTβ,∀β∈Rp}.

BLUE/MVLUE/Gauss-Markov 估计

如果 bTy∈EU(a) 满足 Var(bTy)=mincTy∈EU(a)Var(cTy), 则称 bTy 是 aTβ 的最优线性无偏估计(Best Linear Unbiased Estimation, BLUE), 或者极小方差线性无偏估计(MVLUE), 或者 Gauss-Markov 估计 (GME).

Gauss-Markov 定理

如果假设 Covε=σ2I (不相关、同方差), 且 aTβ 可估, β^ 是正规方程的任一解, 则 aTβ^ 是 aTβ 的唯一的 BLUE.

这表明, 如果 aTβ 可估(a∈μ(X)), 则 aTβ^ 不依赖 β^ 是正规方程的哪个解. 因此, 不管 X 是否满列秩, Xβ^ 总是唯一的.

而如果取消 Covε=σ2I 的假定, 虽然可估性和最小二乘估计没有改变, 但是 aTβ^ 就不一定是 BLUE 了. 此时, 假设 Covε=σ2G, G 是正定阵, 则模型为 y∼(Xβ,σ2G), 不难求出 BLUE: 令 z=G−12y, 原模型变为 z∼(G−12Xβ,σ2I). 此时正规方程 (XTG−1X)β^=XTG−1y 的解 β^^, 则 aTβ 的 BLUE 就是 aTβ^^.
称 β^^ 是新模型的加权最小二乘估计.

回到模型 y∼(Xβ,σ2I), 记 ε^=y−Xβ^ 为剩余/残差, ||ε^||2=||PX⊥y||2 为剩余平方和. 我们在 6.2中 看出 ||ε^||2n−r (r=rankX) 是 σ2 的无偏估计, 记为 σ^2.

从而

定理 1.2

对正态模型 y∼Nn(Xβ,σ2I), aTβ 是可估的, 则 aTβ^, σ^ 分别是 aTβ, σ2 的唯一的 UMVUE.

2 受约束线性模型的参数估计

实际问题中我们也要给 β 一个约束: Hβ=ξ. 这里 H∈Rk×p. 设 rankH=k. 则 H 有右逆 Hr, 即 HHr=Ik. 做变换 z=y−XHrξ, θ=β−Hrξ, 则原模型变为 z∼(Xθ,σ2I;Hθ=0). 由此, 可以不失一般性地讨论 Hβ=0. 下面就来讨论 (2.1)y∼(Xβ,σ2I;Hβ=0).
回顾 定理1.1. 注意到证明中用了 β 的任意性, 因此添加 Hβ=0 约束后, 可估的充要条件变为 ∃c: cTXβ=aTβ, ∀Hβ=0. 它等价于 ∃c: cTXβ=aTβ, ∀β∈μ⊥(HT). 又等价于 XTc−a∈μ(HT), 即 a∈μ((XTHT))=μ(XT)+μ(HT).
因为 Hβ=0, 则也理应有 Hβ^=0. 于是最小二乘推广为受约束最小二乘估计 β^H: ||y−Xβ^H||2=minHb=0||y−Xb||2. 注意到 ||y−Xb||2=||y−Xβ^H||2+||X(β^H−b)||2+2(β^H−b)TXT(y−Xβ^H). 可知 β^H 是方程 (β^H−b)TXT(y−Xβ^H)=0,∀Hb=0 的解. 这等价于 XT(y−Xβ^H)∈μ(HT). 即 ∃λ: XTy−XTXβ^H=HTλ. 从而 β^H 是方程 (2.2)(XTXHTH0)(XHλ)=(XTy0) 的解. 把它作为 (2.1) 的正规方程.

不难证明上述方程是相容(有解)的.

假设 (2.2) 解为 β^H=Ly. 由 Hβ^H=0, 得 HLy=0, 从而 Cov(HLy)=σ2HLLTHT=0, 得 HL=0. 根据正规方程 (2.3)(XTX)Ly+HTλ=XTy⇒LT(XTX)Ly=LTXTy.

又有 LTXTXL=LTXT, 得 XL 是幂等矩阵, 即 PXL. 从而 Xβ^H=PXLy.

对 XL 给出如下引理:

引理 2.1

设 Ly 是 (2.2) 的任一解, R 满足 μ(R)=Ker(H). 则 μ(XL)=μ(XR), 且 rank(XR)=rank(XH)−rankH.

定理 2.1

设 y∼(Xβ,σ2I;Hβ=0), β^H 是受约束最小二乘估计 (不必唯一), aTβ 可估, 则 aTβ^H 是 aTβ 的唯一 BLUE.
记 σ^H2=||y−Xβ^H||2n−s, s 定义见 (2.4), 则 σ^H2 是 σ2 的无偏估计.

2.1 如何附加约束

如何附加约束让原模型不缩小? 设模型为 y=Xβ+ε, β∈Rp. 现在希望 {Xβ|β∈Rp}={Xβ|Hβ=0}={Xβ|β∈Ker(H)}. 根据 引理2.1, {Xβ|β∈Ker(H)}={Xβ|β∈μ(R)}={XRβ|β∈Rp}. 因此模型不缩小的充要条件是 rank(XR)=rank((XTHT))−rankHT=rankXT, 等价于 (2.5)μ(XT)∩μ(HT)={0}.
如果 rankX=r<p, 则 H 秩至多为 p−r.
在 (2.5) 下的约束, 对模型的统计推断没有实质性的影响. 从估计的角度, 我们有下列结果:

定理 2.2

设模型是 (2.1) 且满足 (2.5). 设可估 aTβ 的 BLUE 是 aTβ^H, 则 aTβ^H=aTβ^, 这里 β^ 是 y∼(Xβ,σ2I) 的正规方程 XTXβ^=XTy 的满足 Hβ^=0 的解. 进而, 如果 rank((XTHT))=p, 上面的 β^H,β^ 唯一.

3 区间估计

接下来我们只讨论 y∼Nn(Xβ,σ2I). 先给出一个结果:

定理 3.1

设 cT∈Rm×p, 满足 cp×m=XTKn×m+HTSk×m, ∀K,S, β^H 是 (2.2) 的解, SSHε=||y−Xβ^H||2, 则

  1. cTβ^H∼Nm(KTPXLXβ,σ2KTPXLK).
  2. SSHεσ2∼χn−s2(δ), s=rank((XTHT))−rankHT, δ2=βTXTP(XL)⊥Xβσ2.
  3. cTβ^H⊥⊥SSHε.

特别地, 当 H=0: c=XTK, μ(XL)=μ(X), β^H=β^, 则 cTβ^∼Nm(Xβ,σ2KTPXK), SSεσ2∼χn−r2, cTβ^⊥⊥SSε. 其中 SSε=||y−Xβ^||2, r=rankX.
又如果 Hβ=0, 则 E(cTβ^H)=cTβ, Xβ^∈μ(XL), δ2=0.

现在对于无约束模型 y∼Nn(Xβ,σ2I), 设 cTβ 各分量均可估, 且 rankc=m. 则有 K:c=XTK, 则 KTPXK 是满秩的.

根据定理 (3.1) 立即得 F=(β^−β)Tc(KTPXK)−1cT(β^−β)SSε⋅n−rm∼Fm,n−r, 这是因为 cT(β^−β)=KTX(β^−β)=KTPX(y−Xβ)⇒(β^−β)Tc(KTPXK)−1cT(β^−β)σ2∼χm2.

则 P{F≤Fm,n−r(α)}=1−α.
记 Rm 中椭球 G(cTβ^)={z|(z−cTβ^)T(KTPXK)−1(z−cTβ^)≤mSSεn−rFm,n−r(α)}, 则 P{cTβ∈G(cTβ^)}=1−α, 从而 G(cTβ^) 就是 cTβ 的 1−α 置信系数的置信椭球.

特别地对 m=1, aTβ 可估, 故 ∃b:a=XTb, 于是 aTβ^∼N1(aTβ,σ2bTPXb), 得 T=aT(β^−β)n−rSSε⋅bTPXb∼tn−r. 不难算出置信区间 aTβ^±SSεbTPXbn−r⋅tn−r(α2).

4 一般线性假设检验

此时依然有 y∼Nn(Xβ,σ2I). 考虑一般线性假设 H0:Hβ=0,rankHk×p=k.
考虑 似然比 λ=MMH=max{L(y;β,σ2)|β∈Rp,σ2>0}max{L(y;β,σ2)|Hβ=0,σ2>0}, 其中 L(y;β,σ2)=σ−nexp⁡{−12σ2||y−Xβ||2}. 记SSε=min{||y−Xβ||2|β∈Rp},SSHε=min{||y−Xβ||2|Hβ=0}. 前面得到 SSε=||y−Xβ^||2, SSHε=||y−Xβ^H||2. 然后可以容易得到 M=(SSεn)−n2exp⁡(−n2),MH=(SSHεn)−n2exp⁡(−n2), 从而 λ=MMH=(SSεSSHε)−n2.
对 λ 变形: λ=(SSHε−SSεSSε+1)n2=(SSHSSε+1)n2.
注意到 λ 关于 SSHSSε 严格增加, 似然比检验的否定域 {λ≥c} 易于用这个量表示. 这样再由 定理3.1, F=SSHSSε⋅n−rr−s∼Fr−s,n−r,δ, 其中 r=rankX, s=rank((XTHT))−rankHT, δ2=||Xβ−XEβ^H||2σ2.
当假设成立, F∼Fr−s,n−r, 从而拒绝域为 {F≥Fr−s,n−r(α)}.

5 其他讨论

5.1 方差分析可行的条件

对于上面似然比引出的检验, 也可以用方差分析来解释: ||y||2=||PXLy||2+||(P(XL)⊥−PX⊥)y||2+||PX⊥y||2≡||PXLy||2+SSH+SSε.
这里 SSε 是原模型的剩余平方和, 反映了模型的精确程度; SSH 是附加约束后和原模型的剩余平方和之差, 反映了约束带来的误差情况. 如果客观上 Hβ=0, 约束存在, 则 SSH 应该偏小.

引理 5.1

设 y∼Nn(μ,I), 记 ξ=yTPLy, PL⊥ 是到子空间 L 的正交补空间的投影阵, 则 ξ∼χd2(δ), d=n−dim⁡L, δ2=μTPL⊥μ, 并且 Eξ=d+μTPL⊥μ. 故 μ∈―L 时 ξ 有偏大趋向.

接下来讨论对设计矩阵 X 要求. 设 y∼Nn(Xβ,σ2I). 剖分 X=(X1X2), 相应地 βT=(β(1)Tβ(2)T). 记 μi=μ(Xi), μ=μ(X). 给出方差分析可行的条件:

定理 5.1

||y||2 可以分解为独立的二次型的和 (5.1)||y||2=SSε+SS1+SS2+SSg, 且 SS1,SS2 可以分别解释为由第一/二个因子引起的平方和的充要条件是 (5.2)(μ∩μ1⊥)⊥(μ∩μ2⊥).

推论

设 μ(X1)∩μ(X2)={0}, 则 (5.2) 成立的充要条件是 μ1⊥μ2.

这表明 rankX=p (列满秩) 时, 如果 μ1,μ2 不正交, 将无法分解为 (5.1) 的样子. 由于在回归分析中我们都假定 X 列满秩, 则 XTX 需要有分快对角形: XTX=(X1TX100X2TX2).
而对于方差分析模型, 要求不必如此严格. 例如对 两向分类模型、每格试验次数相同的模型, 我们有 X=(X1⏟rX2⏟c)=(1pc1p⋱⋱1pc1p⋱⋮1pc1p⋱⋱1pc1p).
设 a=(a(1)T,⋯,a(r)T)T∈μ⊥(X1), 则 a(i)T1pc=0, i=1,⋯,r. 而 μ(X) 中的一般元可记为 t=X(uv)=(1pcu1⋮1pcur)+(1pv1 1pvc 1pvc)T, 故 t∈μ(X)∩μ1⊥(X) 时 uipc+p∑j=1cvj=0⇒t=1nu0+(1pv1 1pvc 1pvc)T, 这里 cu0=−∑j=1cvj, n=rcp. 类似地 s∈μ(X)∩μ2⊥(X) 时 s=(1pcw1⋮1pcwr)+1nv0,rv0=−∑i=1rwi. 从而 tTs=0.

5.2 非线性模型

非线性模型有些时候可以转换为线性模型. 考虑 y=F(x1,⋯,xp;β1,⋯,βp)+ε, 这里 β1,⋯,βp 是未知参数, ε 是随机误差. 如果 ∃ 可逆的连续函数 f: f(F(x1,⋯,xp;β1,⋯,βp))=∑i=1pgi(x1,⋯,xp)φi(β1,⋯,βp), 且满足 (x1,⋯,xp)↦(g1,⋯,gp),(β1,⋯,βp)↦(φ1,⋯,φp) 一一对应, 则记 z=f(y), x~i=gi, β~i=φi, i=1,⋯,p. 原模型可以近似看成 z=∑i=1px~iβ~i+ε~.