我们在 6.1 线性模型的概念和分类 开始讨论了一个因变量的线性模型. 现在推广到多个因变量.

1 多元线性模型

设计矩阵和之前相同: ${\displaystyle X=\left(\begin{array}{ccc}{x}_{11}& \cdots & x_{1q}\\ ⋮& & ⋮\\ x_{n1}& \cdots & x_{nq}\end{array}\right)}$. 而这里 ${\displaystyle n}$ 个观察值 ${\displaystyle y_1,\cdots ,y_n\in {\mathbb{R}}^p}$, 记为 ${\displaystyle Y=\left(\begin{array}{ccc}{y}_{11}& \cdots & y_{1p}\\ ⋮& & ⋮\\ y_{n1}& \cdots & y_{np}\end{array}\right)}$. 误差矩阵 ${\displaystyle \epsilon =\left(\begin{array}{ccc}{\epsilon }_{11}& \cdots & {\epsilon }_{1p}\\ ⋮& & ⋮\\ {\epsilon }_{n1}& \cdots & {\epsilon }_{np}\end{array}\right)}$. 未知参数矩阵 ${\displaystyle B=\left(\begin{array}{ccc}{\beta }_{11}& \cdots & {\beta }_{1p}\\ ⋮& & ⋮\\ {\beta }_{q1}& \cdots & {\beta }_{qp}\end{array}\right)}$. ${\displaystyle {\beta }_j}$ (第 ${\displaystyle j}$ 列) 影响实验结果的第 ${\displaystyle j}$ 个指标.
记 ${\displaystyle Y}$ 的第 ${\displaystyle j}$ 列为 ${\displaystyle Y_j}$, 则 $$Y_j=X{\beta }_j+{\epsilon }_j,j=1,\cdots ,p.$$
看起来多元线性模型就是 ${\displaystyle p}$ 个一元线性模型的混合, 但是不同之处是我们要考虑 ${\displaystyle p}$ 个指标之间存在相关关系.
设观察值矩阵 ${\displaystyle Y}$ 的各行互不相关, 有相同的协方差阵 ${\displaystyle \mathrm{\Sigma }}$: $$\mathrm{Cov}(y_{\alpha },y_{\beta })={\delta }_{\alpha \beta }\mathrm{\Sigma }.$$
因此我们修改记法为

它就是一元线性模型的推广, 把 ${\displaystyle \beta ,{\sigma }^2>0}$ 推广为 ${\displaystyle B,\mathrm{\Sigma }>\mathbf{0}}$.

2 参数估计及其分布

大多数情形下考虑列满秩情形: ${\displaystyle \mathrm{rank}X=q}$.
下面讨论 ${\displaystyle B}$ 的估计. 设估计为 ${\displaystyle \hat{B}}$, 则残差为 ${\displaystyle \hat{\epsilon }=Y-X\hat{B}}$. 可以定义残差阵 ${\displaystyle {\hat{\epsilon }}^{\mathrm{T}}\hat{\epsilon }}$. 我们希望好的 ${\displaystyle \hat{B}}$ 满足: $$(Y-XB{)}^{\mathrm{T}}(Y-XB)-(Y-X\hat{B}{)}^{\mathrm{T}}(Y-X\hat{B})\ge \mathbf{0}$$ (也即非负定), 或记为 $$(Y-XB{)}^{\mathrm{T}}(Y-XB)\ge (Y-X\hat{B}{)}^{\mathrm{T}}(Y-X\hat{B}),\mathrm{\forall }B.$$

这是一个很强的结果. 当 ${\displaystyle C\ge D\ge 0}$, 它们由大到小的特征值 ${\displaystyle {\lambda }_i(C),{\lambda }_i(D)}$ 满足 ${\displaystyle {\lambda }_i(C)\ge {\lambda }_i(D)}$, ${\displaystyle i=1,\cdots ,\mathrm{rank}D}$. 因此 ${\displaystyle \prod _{i=1}^r{\lambda }_i(C)\ge \prod _{i=1}^r{\lambda }_i(D)}$, ${\displaystyle \mathrm{tr}C\ge \mathrm{tr}D}$ 等性质成立.

类似 Gauss-Markov定理:

2.1 正态性假设下的讨论

采用 (2.1) 的记号, 设 $$\epsilon \sim {\mathcal{N}}_{np}(0,I,\mathrm{\Sigma }),$$ 接下来讨论这样的正态模型.

和 这个定理 类似, 可以证明 ${\displaystyle (\hat{B},\frac{R_0}{n-q})}$ 是 ${\displaystyle (B,\mathrm{\Sigma })}$ 的 UMVUE, 方法完全类似. 在没有正态性的假设下, 依然有 ${\displaystyle \frac{R_0}{n-q}}$ 是 ${\displaystyle \mathrm{\Sigma }}$ 的无偏估计.
事实上, 因为 ${\displaystyle P_{X^{\perp }}}$ 是 ${\displaystyle n-q}$ 秩的正投影阵, 必然存在正交 ${\displaystyle U}$: $$U^{\mathrm{T}}{P}_{X^{\perp }}U=\left(\begin{array}{cc}{I}_{n-q}& 0\\ 0& 0\end{array}\right),$$ 于是 $$R_0=Y^{\mathrm{T}}U{U}^{\mathrm{T}}{P}_{X^{\perp }}U{U}^{\mathrm{T}}Y=Y^{\mathrm{T}}U\left(\begin{array}{cc}{I}_{n-q}& 0\\ 0& 0\end{array}\right)U^{\mathrm{T}}Y.$$
注意到 ${\displaystyle Z=U^{\mathrm{T}}Y\sim {\mathcal{N}}_{np}(0,I,\mathrm{\Sigma })}$, 记 ${\displaystyle Z^{\mathrm{T}}=(Z_1,\cdots ,Z_n)}$, 有 ${\displaystyle R_0=\sum _{\alpha =1}^{n-q}{Z}_{\alpha }{Z}_{\alpha }^{\mathrm{T}}}$, 因此 $$\mathbb{E}{R}_0=\sum _{\alpha =1}^{n-q}\mathbb{E}[Z_{\alpha }{Z}_{\alpha }^{\mathrm{T}}]=(n-q)\mathrm{\Sigma }.$$

2.2 多元情形的推广

设因变量 ${\displaystyle Y_1,\cdots ,Y_p}$ 对 ${\displaystyle X_1,\cdots ,X_q}$ 的统计依赖表现为: 给定 ${\displaystyle X_1,\cdots ,X_q}$, 有 $$\begin{array}{}\text{(2.2)}& (Y_1,\cdots ,Y_p)=(1,X_1,\cdots ,X_q)B+({\epsilon }_1,\cdots ,{\epsilon }_p),\end{array}$$ 其中 ${\displaystyle B=\left(\begin{array}{ccc}{\beta }_{01}& \cdots & {\beta }_{0p}\\ {\beta }_{11}& \cdots & {\beta }_{1p}\\ ⋮& & ⋮\\ {\beta }_{q1}& \cdots & {\beta }_{qp}\end{array}\right)}$. 称 ${\displaystyle {\beta }_{ij}}$ 是自变量 ${\displaystyle X_i}$ 对因变量 ${\displaystyle Y_j}$ 的效应. (2.2) 称为总体回归模型.
根据 性质5, 当 ${\displaystyle (Y_1,\cdots ,Y_p,X_1,\cdots ,X_q)=Z\sim {\mathcal{N}}_{p+q}(v,V)}$ 时, 记 $$v=\left(\begin{array}{c}\mathbb{E}Y\\ \mathbb{E}X\end{array}\right),V=\left(\begin{array}{cc}{\mathrm{\Sigma }}_{YY}& {\mathrm{\Sigma }}_{YX}\\ {\mathrm{\Sigma }}_{XY}& {\mathrm{\Sigma }}_{XX}\end{array}\right),V>0,$$ 则 ${\displaystyle X_1,\cdots ,X_q}$ 给定时, $$Y|X\sim {\mathcal{N}}_p(\mathbb{E}Y+{\mathrm{\Sigma }}_{YX}{\mathrm{\Sigma }}_{XX}^{-1}(X-\mathbb{E}X),{\mathrm{\Sigma }}_{YY}-{\mathrm{\Sigma }}_{YX}{\mathrm{\Sigma }}_{XX}^{-1}{\mathrm{\Sigma }}_{XY}),$$ 这里 ${\displaystyle Y=(Y_1,\cdots ,Y_p{)}^{\mathrm{T}}}$. 可见总体回归函数是 $$\begin{array}{}\text{(2.3)}& \mathbb{E}Y+{\mathrm{\Sigma }}_{YX}{\mathrm{\Sigma }}_{XX}^{-1}(X-\mathbb{E}X).\end{array}$$ 不难看出, 以回归方程 $$\hat{Y}=\mathbb{E}Y+{\mathrm{\Sigma }}_{YX}{\mathrm{\Sigma }}_{XX}^{-1}(X-\mathbb{E}X)$$ 做 ${\displaystyle Y}$ 的预测, 预测误差的协方差阵是 $${\mathrm{\Sigma }}_{YY}-{\mathrm{\Sigma }}_{YX}{\mathrm{\Sigma }}_{XX}^{-1}{\mathrm{\Sigma }}_{XY}.$$ 这些都是一元情形的自然推广.

类似地, 由最小二乘估计得回归系数的估计$$\begin{array}{rl}\hat{B}& =[(1X{)}^{\mathrm{T}}(1X){]}^{-1}(1X{)}^{\mathrm{T}}Y\\ & ={\left(\begin{array}{cc}n& 1^{\mathrm{T}}X\\ X^{\mathrm{T}}1& X^{\mathrm{T}}X\end{array}\right)}^{-1}(1X{)}^{\mathrm{T}}Y\\ & =\left(\begin{array}{cc}\frac{1}{n}+\frac{{\mathbf{1}}^{\mathrm{T}}X}{n}(X^{\mathrm{T}}{P}_{1^{\perp }}X{)}^{-1}\frac{X^{\mathrm{T}}\mathbf{1}}{n}& -\frac{{\mathbf{1}}^{\mathrm{T}}X}{n}(X^{\mathrm{T}}{P}_{{\mathbf{1}}^{\perp }}X{)}^{-1}\\ -(X^{\mathrm{T}}{P}_{{\mathbf{1}}^{\perp }}X{)}^{-1}\frac{X^{\mathrm{T}}\mathbf{1}}{n}& (X^{\mathrm{T}}{P}_{1^{\perp }}X{)}^{-1}\end{array}\right)\left(\begin{array}{c}{\mathbf{1}}^{\mathrm{T}}Y\\ X^{\mathrm{T}}Y\end{array}\right)\\ & =\left(\begin{array}{c}\frac{{\mathbf{1}}^{\mathrm{T}}Y}{n}-\frac{{\mathbf{1}}^{\mathrm{T}}X}{n}(X^{\mathrm{T}}{P}_{1^{\perp }}X{)}^{-1}(X^{\mathrm{T}}{P}_{1^{\perp }}Y)\\ (X^{\mathrm{T}}{P}_{{\mathbf{1}}^{\perp }}X{)}^{-1}(X^{\mathrm{T}}{P}_{1^{\perp }}Y)\end{array}\right).\end{array}$$
这和总体回归情形用样本矩代替总体矩得到的结果完全一致, 和一元情形也相仿.

3 线性假设检验

3.1 检验参数矩阵

先讨论 (2.2) 的假设检验. 此时假设为 $$H_0:B_1=\left(\begin{array}{ccc}{\beta }_{11}& \cdots & {\beta }_{1p}\\ ⋮& & ⋮\\ {\beta }_{q1}& \cdots & {\beta }_{qp}\end{array}\right)=\mathbf{0}.$$
对正态模型, 令似然比 $$\lambda =\frac{max\{L(Y,B,\mathrm{\Sigma })|B,\mathrm{\Sigma }>\mathbf{0}\}}{max\{L(Y,{\beta }_0,\mathrm{\Sigma })|{\beta }_0,\mathrm{\Sigma }>\mathbf{0}\}}\equiv \frac{M}{M_H}.$$ 其中 ${\displaystyle {\beta }_0}$ 是 ${\displaystyle B}$ 的第一个行向量. 由 定理2.3, $$M=C{[det\left(\frac{R_0}{n}\right)]}^{-\frac{n}{2}}\exp \{-\frac{n}{2}\}.$$
而 ${\displaystyle H_0}$ 成立时 $$L(Y,{\beta }_0,\mathrm{\Sigma })=C(det\mathrm{\Sigma }{)}^{-\frac{n}{2}}\exp \{-\frac{1}{2}\mathrm{tr}({\mathrm{\Sigma }}^{-1}\sum _{\alpha =1}^n(y_{\alpha }-{\beta }_0)(y_{\alpha }-{\beta }_0{)}^{\mathrm{T}})\}.$$
注意到$$\begin{array}{rl}\sum _{\alpha =1}^n(y_{\alpha }-{\beta }_0)(y_{\alpha }-{\beta }_0{)}^{\mathrm{T}}& =\sum _{\alpha =1}^n(y_{\alpha }-\bar{y})(y_{\alpha }-\bar{y}{)}^{\mathrm{T}}+n(\bar{y}-{\beta }_0)(\bar{y}-{\beta }_0{)}^{\mathrm{T}}\\ & =Y^{\mathrm{T}}{P}_{{\mathbf{1}}^{\perp }}Y+n(\bar{y}-{\beta }_0)(\bar{y}-{\beta }_0{)}^{\mathrm{T}},\end{array}$$ 得 $$M_H=C{[det\left(\frac{Y^{\mathrm{T}}{P}_{{\mathbf{1}}^{\perp }}Y}{n}\right)]}^{-\frac{n}{2}}\exp \{-\frac{n}{2}\}.$$
记 ${\displaystyle Y^{\mathrm{T}}{P}_{{\mathbf{1}}^{\perp }}Y=R_1}$. 有 ${\displaystyle \lambda ={\left(\frac{det{R}_0}{det{R}_1}\right)}^{-\frac{n}{2}}}$. 它是 ${\displaystyle \mathrm{\Lambda }=\frac{det{R}_0}{det{R}_1}}$ 的严格下降函数, 因此拒绝域为 ${\displaystyle \{\mathrm{\Lambda }\le C_{\alpha }\}}$. 由于 ${\displaystyle R_0\sim W_p(n-q-1,\mathrm{\Sigma })}$,^[2] 当 ${\displaystyle H_0}$ 成立时, ${\displaystyle R_1\sim W_p(n-1,\mathrm{\Sigma })}$, 得 ${\displaystyle R_1-R_0\perp \perp R_0}$, 且 ${\displaystyle R_1-R_0\sim W_p(q,\mathrm{\Sigma })}$. 故 $$\mathrm{\Lambda }\sim \mathrm{\Lambda }(p,n-q-1,q).$$

3.2 检验单一效应

当 ${\displaystyle H_0}$ 被拒绝, 接受回归模型 (2.2). 类似一元情形, 检验 ${\displaystyle X_i}$ 效应的显著性: $$H_{0i}={\beta }_{i1}=\cdots ={\beta }_{ip}=0,i=1,\cdots ,q.$$ 类似上面, 用似然比检验, 容易验证, ${\displaystyle H_{0i}}$ 成立时, ${\displaystyle R_{1i}=Y^{\mathrm{T}}{P}_{{\mu }_i^{\perp }}Y}$, 这里 ${\displaystyle {\mu }_i}$ 是 ${\displaystyle (1X)}$ 划去 ${\displaystyle X}$ 的第 ${\displaystyle i}$ 列后长成的线性空间, 于是有 ${\displaystyle R_{1i}\sim W_p(n-q,\mathrm{\Sigma },\mathrm{\Delta })}$. 当假设成立时, 有 ${\displaystyle R_{1i}\sim W_p(n-q,\mathrm{\Sigma })}$, ${\displaystyle R_{1i}-R_0\sim W_p(1,\mathrm{\Sigma })}$. 于是有 $$\mathrm{\Lambda }=\frac{det{R}_0}{det{R}_{1i}}\sim \mathrm{\Lambda }(p,n-q-1,1).$$ 令 ${\displaystyle F=\frac{1-\mathrm{\Lambda }}{\mathrm{\Lambda }}\cdot \frac{n-q-p}{p}}$, 当 ${\displaystyle H_{0i}}$ 成立时 $$\frac{1-\mathrm{\Lambda }}{\mathrm{\Lambda }}\cdot \frac{n-q-p}{p}\sim F_{p,n-p-q},$$ 故得拒绝域 ${\displaystyle \{{\mathrm{\Lambda }}^{-1}\ge 1+\frac{p}{n-p-q}\cdot F_{p,n-p-q,\alpha }\}}$.^[2:1]

3.3 变量的选择

在多元情形, 就算自变量 ${\displaystyle X_i}$ 对整个因变量集的效应是显著的, 也不见得它对单个因变量都是显著的. 从而我们可以检验 ${\displaystyle B}$ 的列向量是否为 ${\displaystyle \mathbf{0}}$. 另外就算自变量对整个 ${\displaystyle Y}$ 不显著, 对单个因变量也可能显著. 所以问题变得很复杂.

在这里我们讨论一般线性假设 ${\displaystyle H_0:HBG=\mathbf{0}}$.
这里 ${\displaystyle H\in {\mathbb{R}}^{k\times q}}$, ${\displaystyle G\in {\mathbb{R}}^{p\times r}}$, 且 ${\displaystyle \mathrm{rank}H=k}$, ${\displaystyle \mathrm{rank}G=r}$.

先讨论一般线性假设 ${\displaystyle H_0:HB=\mathbf{0}}$. 此约束下的 ${\displaystyle {\hat{B}}_H}$, 相当于在约束 ${\displaystyle H{\beta }_j=\mathbf{0}}$ 下的最小二乘估计, 故有 ${\displaystyle X{\hat{B}}_H=P_{XQ}Y}$, 这里 ${\displaystyle Q}$ 是 ${\displaystyle q}$ 阶方针: ${\displaystyle \mathrm{Im}(Q)=\mathrm{Ker}(H)}$. 实际上可取 $$Q=I-H^{\mathrm{T}}(H{H}^{\mathrm{T}}{)}^{-1}H\equiv P_{{H^{\mathrm{T}}}^{\perp }},$$ 它是到 ${\displaystyle \mathrm{Im}(H^{\mathrm{T}})}$ 的正交补空间的投影阵. 记 ${\displaystyle {\hat{B}}_H}$ 的残积阵为 ${\displaystyle R_H}$, 则 $$R_H=Y^{\mathrm{T}}{P}_{(XQ{)}^{\perp }}Y,$$ 从而仿照前面, 似然比 ${\displaystyle \lambda =\frac{M}{M_H}}$ 是 ${\displaystyle \mathrm{\Lambda }=\frac{det{R}_0}{det{R}_H}}$ 的严格下降函数. 因此 ${\displaystyle HB=\mathbf{0}}$ 的拒绝域为 ${\displaystyle \{\mathrm{\Lambda }\le C_{\alpha }\}}$. 根据 (2.4), 知 $$\mathrm{rank}{P}_{(XQ{)}^{\perp }}=n-[\mathrm{rank}((X^{\mathrm{T}}{H}^{\mathrm{T}}))-\mathrm{rank}{H}^{\mathrm{T}}]\equiv n-s,$$ 因此 ${\displaystyle H_0}$ 成立时 ${\displaystyle \mathrm{\Lambda }\sim \mathrm{\Lambda }(p,n-q,q-s)}$. 只有在特殊情况下才等价于 ${\displaystyle F}$ 检验.

回到 ${\displaystyle H_0:HBG=\mathbf{0}}$. 令 ${\displaystyle Z=YG,\mathrm{\Theta }=BG}$, 则模型变为 $$Z=X\mathrm{\Theta }+\epsilon G,\epsilon G\sim {\mathcal{N}}_{nr}(0,I,G^{\mathrm{T}}\mathrm{\Sigma }G).$$ 假设变为 ${\displaystyle {\tilde{H}}_0:H\mathrm{\Theta }=\mathbf{0}}$. 从而类似上面, 检验统计量为 $$\mathrm{\Lambda }=\frac{det(Z^{\mathrm{T}}{P}_{X^{\perp }}Z)}{det(Z^{\mathrm{T}}{P}_{(XQ{)}^{\perp }}Z)}=\frac{det(G^{\mathrm{T}}{Y}^{\mathrm{T}}{P}_{X^{\perp }}YG)}{det(G^{\mathrm{T}}{Y}^{\mathrm{T}}{P}_{(XQ{)}^{\perp }}YG)}.$$${\displaystyle H_0}$ 成立时 ${\displaystyle \mathrm{\Lambda }\sim \mathrm{\Lambda }(r,n-q,q-s)}$.

4 广义方差分析

回顾 6.3 方差分析, 我们把某个值进行平方和分解. 在这里自然推广为 ${\displaystyle Y^{\mathrm{T}}Y=\sum _{i=1}^k{W}_i}$. 这里 ${\displaystyle W_i}$ 是相互独立的遵从 Wishart分布 的随机矩阵. 根据 Cochran定理在矩阵的推广, 利用 Wilks 统计量 ${\displaystyle \mathrm{\Lambda }}$ 可以按照分解式进行一些检验.
例如两向分类模型 $$y_{ij}={\theta }_0+{\beta }_i+{\gamma }_j+{\epsilon }_{ij},i=1,\cdots ,r;j=1,\cdots ,c,$$ 其中 ${\displaystyle y_{ij},{\theta }_0,{\beta }_i,{\gamma }_j,{\epsilon }_{ij}\in {\mathbb{R}}^p}$. 容易得到$$\begin{array}{rl}{Y}^{\mathrm{T}}Y=& \bar{y}{\bar{y}}^{\mathrm{T}}+\sum _{i=1}^{r}C({\bar{y}}_{i\cdot }-\bar{y})({\bar{y}}_{i\cdot }-\bar{y}{)}^{\mathrm{T}}+\sum _{j=1}^{c}r({\bar{y}}_{\cdot j}-\bar{y})({\bar{y}}_{\cdot j}-\bar{y}{)}^{\mathrm{T}}\\ & +\sum _{i=1}^r\sum _{j=1}^c(y_{ij}-{\bar{y}}_{i\cdot }-{\bar{y}}_{\cdot j}+\bar{y})(y_{ij}-{\bar{y}}_{i\cdot }-{\bar{y}}_{\cdot j}+\bar{y}{)}^{\mathrm{T}}\\ =& W_0+W_r+W_c+W_{\epsilon },\end{array}$$ 其中 ${\displaystyle \bar{y}=\frac{1}{rc}\sum _{i=1}^r\sum _{j=1}^c{y}_{ij}}$, ${\displaystyle {\bar{y}}_{i\cdot }=\frac{1}{c}\sum _{j=1}^c{y}_{ij}}$, ${\displaystyle {\bar{y}}_{\cdot j}=\frac{1}{r}\sum _{i=1}^r{y}_{ij}}$. 想要检验$$\begin{array}{rl}& H_{01}:{\beta }_1=\cdots ={\beta }_r=\mathbf{0},\\ & H_{02}:{\gamma }_1=\cdots ={\gamma }_c=\mathbf{0}.\end{array}$$
${\displaystyle H_{01}}$ 的统计量可以取 ${\displaystyle {\mathrm{\Lambda }}_1=\frac{det{W}_{\epsilon }}{det(W_{\epsilon }+W_r)}}$; ${\displaystyle H_{02}}$ 的可以取 ${\displaystyle {\mathrm{\Lambda }}_2=\frac{det{W}_{\epsilon }}{det(W_{\epsilon }+W_c)}}$. 零假设成立时它们分别遵从 ${\displaystyle \mathrm{\Lambda }}$ 分布.

一般地, 广义方差分析总可以从方差分析导出, 这是因为 Wishart 矩阵和 ${\displaystyle {\chi }^2}$ 变量存在内在联系. 如果想分解 ${\displaystyle Y^{\mathrm{T}}Y=\sum _{i=1}^k{Y}^{\mathrm{T}}{P}_{{\mu }_i}Y}$, 可以考虑 ${\displaystyle Ya=XBa+\epsilon a}$. 记 ${\displaystyle Z=Ya}$, 有分解 $$Z^{\mathrm{T}}Z=\sum _{i=1}^k{Z}^{\mathrm{T}}{P}_{{\mu }_i}Z,$$ 满足 ${\displaystyle R^n={\mu }_1\oplus \cdots \oplus {\mu }_k}$ (正交直和).