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

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

1 可估参数函数及其估计

我们从 β 的最小二乘估计开始. 在 这里 我们得到它是正规方程 (1.1)XTXβ^=XTy 的解. 现在解释它的直观几何意义. 要使 ||yXβ||
最小, 则 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, 则 aRp, aTβ 都是可估的; 如果 rankX<p, 则 β 至少有一个分量不可估.

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

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

BLUE/MVLUE/Gauss-Markov 估计

如果 bTyEU(a) 满足 Var(bTy)=mincTyEU(a)Var(cTy), 则称 bTyaTβ最优线性无偏估计(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=G12y, 原模型变为 z(G12Xβ,σ2I). 此时正规方程 (XTG1X)β^=XTG1y 的解 β^^, 则 aTβ 的 BLUE 就是 aTβ^^.
β^^ 是新模型的加权最小二乘估计.

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

从而

定理 1.2

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

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

实际问题中我们也要给 β 一个约束: Hβ=ξ. 这里 HRk×p. 设 rankH=k. 则 H 有右逆 Hr, 即 HHr=Ik. 做变换 z=yXHrξ, θ=β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). 又等价于 XTcaμ(HT), 即 aμ((XTHT))=μ(XT)+μ(HT).
因为 Hβ=0, 则也理应有 Hβ^=0. 于是最小二乘推广为受约束最小二乘估计 β^H: ||yXβ^H||2=minHb=0||yXb||2. 注意到 ||yXb||2=||yXβ^H||2+||X(β^Hb)||2+2(β^Hb)TXT(yXβ^H). 可知 β^H 是方程 (β^Hb)TXT(yXβ^H)=0,Hb=0 的解. 这等价于 XT(yXβ^H)μ(HT). 即 λ: XTyXTXβ^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λ=XTyLT(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β^HaTβ 的唯一 BLUE.
σ^H2=||yXβ^H||2ns, 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 秩至多为 pr.
在 (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 区间估计

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

定理 3.1

cTRm×p, 满足 cp×m=XTKn×m+HTSk×m, K,S, β^H(2.2) 的解, SSHε=||yXβ^H||2, 则

  1. cTβ^HNm(KTPXLXβ,σ2KTPXLK).
  2. SSHεσ2χns2(δ), s=rank((XTHT))rankHT, δ2=βTXTP(XL)Xβσ2.
  3. cTβ^HSSHε.

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

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

根据定理 (3.1) 立即得 F=(β^β)Tc(KTPXK)1cT(β^β)SSεnrmFm,nr, 这是因为 cT(β^β)=KTX(β^β)=KTPX(yXβ)(β^β)Tc(KTPXK)1cT(β^β)σ2χm2.

P{FFm,nr(α)}=1α.
Rm 中椭球 G(cTβ^)={z|(zcTβ^)T(KTPXK)1(zcTβ^)mSSεnrFm,nr(α)},P{cTβG(cTβ^)}=1α, 从而 G(cTβ^) 就是 cTβ1α 置信系数的置信椭球.

特别地对 m=1, aTβ 可估, 故 b:a=XTb, 于是 aTβ^N1(aTβ,σ2bTPXb), 得 T=aT(β^β)nrSSεbTPXbtnr. 不难算出置信区间 aTβ^±SSεbTPXbnrtnr(α2).

4 一般线性假设检验

此时依然有 yNn(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||yXβ||2}. 记SSε=min{||yXβ||2|βRp},SSHε=min{||yXβ||2|Hβ=0}. 前面得到 SSε=||yXβ^||2, SSHε=||yXβ^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εnrrsFrs,nr,δ, 其中 r=rankX, s=rank((XTHT))rankHT, δ2=||XβXEβ^H||2σ2.
当假设成立, FFrs,nr, 从而拒绝域为 {FFrs,nr(α)}.

5 其他讨论

5.1 方差分析可行的条件

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

引理 5.1

yNn(μ,I), 记 ξ=yTPLy, PL 是到子空间 L 的正交补空间的投影阵, 则 ξχd2(δ), d=ndimL, δ2=μTPLμ, 并且 Eξ=d+μTPLμ.μLξ 有偏大趋向.

接下来讨论对设计矩阵 X 要求. 设 yNn(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=(X1rX2c)=(1pc1p1pc1p1pc1p1pc1p).
a=(a(1)T,,a(r)T)Tμ(X1), 则 a(i)T1pc=0, i=1,,r. 而 μ(X) 中的一般元可记为 t=X(uv)=(1pcu11pcur)+(1pv1 1pvc 1pvc)T,tμ(X)μ1(X)uipc+pj=1cvj=0t=1nu0+(1pv1 1pvc 1pvc)T, 这里 cu0=j=1cvj, n=rcp. 类似地 sμ(X)μ2(X)s=(1pcw11pcwr)+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+ε~.