概率论拾遗：自由度

摘 要 在线性回归模型的方差分析中, 自由度一般被形象地表述为自由取值的变量个数, 即变量总数减去限制条件的数量. 这种表述虽然易于理解但是并不本质. 本文讲述自由度的一个严格定义: 方差分析中的每一个平方和都可以写成二次型的形式, 而自由度是这些二次型矩阵的秩.

1 多元线性回归模型

1.1 多元线性回归模型及其矩阵表示

  给定 $p-1$ 个解释变量 $x_1,\dots,x_{p-1}$ 和对应的响应变量 $y$. 假定有 $n$ 条观测数据

$${\rm OBS}=\{(x_1^{(i)},\dots,x_{p-1}^{(i)},y^{(i)})\}_{i=1}^n$$

建立多元线性回归模型:

$$y^{(i)}=b_0+b_1x_1^{(i)}+\dots+b_{p-1}x_{p-1}^{(i)}+e^{(i)}$$

其中 $e^{(i)}$ 是残差项. 我们希望 $\{e^{(i)}\}_{i=1}^n$ 是来自同一个正态随机变量 $\varepsilon\sim N(0,\sigma^2)$ 的简单随机样本.

  现在用矩阵语言重述上面的内容. 定义资料矩阵

$$X=\begin{pmatrix} 1 & x_1^{(1)} & \cdots & x_{p-1}^{(1)}\\ \vdots & \vdots & & \vdots\\ 1 & x_1^{(n)} & \cdots & x_{p-1}^{(n)} \end{pmatrix}$$

即每行是一个观测对象, 每列是一个解释变量. 定义响应变量向量 $\boldsymbol y=(y^{(1)},\dots,y^{(n)})^T$、系数向量 $\boldsymbol b=(b_0,b_1,\dots,b_{p-1})^T$、误差项向量 $\boldsymbol e=(e^{(1)},\dots,e^{(n)})^T$, 则原模型可以写成矩阵形式

$$\boldsymbol y=X\boldsymbol b+\boldsymbol e$$

1.2 最小二乘法

  最小二乘法的目的是找到一组 $\boldsymbol b$ 使得残差平方和最小. 定义残差平方和

$$\begin{aligned}{\rm SSE}:=\boldsymbol e^T\boldsymbol e&=(\boldsymbol y-X\boldsymbol b)^T(\boldsymbol y-X\boldsymbol b)\\&=\boldsymbol b^TX^TX\boldsymbol b-2(X^T\boldsymbol y)^T\boldsymbol b+\boldsymbol y^T\boldsymbol y\end{aligned}$$

这是一个关于 $\boldsymbol b$ 的二次函数.

一般的二次函数形式及其梯度为

$$f=\frac 12\boldsymbol x^TQ\boldsymbol x+\boldsymbol q^T\boldsymbol x+q\quad \implies \quad \nabla f=Q\boldsymbol x+\boldsymbol q$$

现在将 $\rm SSE$ 对 $\boldsymbol b$ 求梯度, 得到

$$\nabla {\rm SSE}=2X^TX\boldsymbol b-2X^T\boldsymbol y$$

令其为零, 得到最优的一组 $\boldsymbol b$

$$\nabla {\rm SSE}=0\quad\implies\quad\boldsymbol b=(X^TX)^{-1}X^T\boldsymbol y$$

1.3 帽子矩阵

  首先用矩阵语言表示拟合值 $\hat{\boldsymbol y}$. 由

$$\hat{\boldsymbol y}:=X\boldsymbol b=\underbrace{X(X^TX)^{-1}X^T}_H\boldsymbol y=:H\boldsymbol y$$

其中定义 $H:=X(X^TX)^{-1}X^T$ 为帽子矩阵. 显然帽子矩阵是对称的. 另外注意到 $H^2=H^T=H$, 即帽子矩阵是对称幂等的.

对称幂等阵 满足 $A^2=A^T=A$ 的矩阵 $A$ 称为对称幂等阵. 对称幂等阵有一些很好的性质. 令 $A$ 是一个对称幂等阵, 则:

• $A,I-A,\mathop{\rm adj}A,A^k(k\in\mathbb N)$ 都是对称幂等的;
• $A\sim\begin{pmatrix}I\\&O\end{pmatrix}$;
• $\mathop{\rm rank}A=\mathop{\rm tr}A$.

  帽子矩阵是一个投影矩阵.

投影矩阵 对于任意列满秩的矩阵 $A_{m\times n}(m>n)$, 满足形式 $P=A(A^TA)^{-1}A^T$ 形式的矩阵称为投影矩阵, 它在几何上将全空间 $\mathbb R^m$ 内的所有向量垂直地投影到 $A$ 的列空间, 即 $\mathop{\rm Im}P=\mathop{\rm Im}A$

等式两边取维数, 可以推出 $\mathop{\rm rank}P=\mathop{\rm rank}A=n$.

应用于帽子矩阵, 得到 $\mathop{\rm rank}H=\mathop{\rm rank}X=p$ (假设存在足够多的观测数据使得 $X$ 列满秩).

  有了帽子矩阵后, 残差向量就可以表示为

$$\boldsymbol e=\boldsymbol y-\hat{\boldsymbol y}=(I-H)\boldsymbol y$$

2 方差分析

  在多元线性回归模型中, 我们有 ANOVA 方差分析表:

	平方和	自由度	均方
Regression	${\rm SSR}=\sum(\hat y^{(i)}-\bar y)^2$	$p-1$	${\rm MSR}=\frac{\rm SSR}{p-1}$
Error	${\rm SSE}=\sum(y^{(i)}-\hat y^{(i)})^2$	$n-p$	${\rm MSE}=\frac{\rm SSE}{n-p}$
Total	${\rm SST}=\sum(y^{(i)}-\bar y)^2$	$n-1$

2.1 平方和的二次型表示

  首先回顾前面定义的残差平方和 $\rm SSE$

$$\begin{aligned}{\rm SSE}&:=\sum(y^{(i)}-\hat y^{(i)})^2\\&=\boldsymbol e^T\boldsymbol e=\boldsymbol y^T(I-H)^T(I-H)\boldsymbol y=\boldsymbol y^T(I-H)\boldsymbol y\end{aligned}$$

然后定义总平方和, 注意到

$$\bar y=\frac{\boldsymbol 1^T\boldsymbol y}{n}$$

则有

$$\begin{aligned} {\rm SST}&:=\sum(y^{(i)}-\bar y)^2=\sum y^{(i)}(y^{(i)}-\bar y)\\ &=\boldsymbol y^T\left(\boldsymbol y-\boldsymbol 1\frac{\boldsymbol 1^T\boldsymbol y}{n}\right)=\boldsymbol y^T\left(I-\underbrace{\boldsymbol 1\frac{\boldsymbol 1^T}{n}}_J\right)\boldsymbol y=\boldsymbol y^T(I-J)\boldsymbol y\end{aligned}$$

其中定义

$$J:=\frac{\boldsymbol 1\boldsymbol 1^T}{n}=\frac{\boldsymbol 1\boldsymbol 1^T}{\boldsymbol 1^T\boldsymbol 1}$$

最后定义回归平方和

$${\rm SSR}:=\sum(\hat y^{(i)}-\bar y)^2={\rm SST}-{\rm SSE}=\boldsymbol y^T(H-J)\boldsymbol y$$

如此一来, 三个平方和可以表示为 $\boldsymbol y$ 的二次型的形式:

$${\rm SSR}=\boldsymbol y^T(H-J)\boldsymbol y,\qquad {\rm SSE}=\boldsymbol y^T(I-H)\boldsymbol y,\qquad {\rm SST}=\boldsymbol y^T(I-J)\boldsymbol y$$

2.2 二次型矩阵的秩

  残差平方和的二次型矩阵 $I-H$ 是对称幂等的, 因为 $H$ 是对称幂等的. $I-H$ 的秩为

$$\mathop{\rm rank}(I-H)=\mathop{\rm tr}(I-H)=\mathop{\rm tr}I-\mathop{\rm tr}H=n-p$$

  然后是总平方和. 注意到 $J$ 是对称幂等的, 因为

$$J^2=\frac{\boldsymbol 1\boldsymbol 1^T\boldsymbol 1\boldsymbol 1^T}{n^2}=\frac{\boldsymbol 1\boldsymbol 1^T}{n}=J$$

所以总平方和的二次型矩阵 $I-J$ 是对称幂等的. 它的秩为

$$\mathop{\rm rank}(I-J)=\mathop{\rm tr}(I-J)=\mathop{\rm tr}I-\mathop{\rm tr}J=n-1$$

  最后是回归平方和. 注意到 $HX=X(X^TX)^{-1}X^TX=X$, 对 $X$ 作列分块得

$$H\begin{pmatrix}\boldsymbol 1 & * & \cdots & *\end{pmatrix}=\begin{pmatrix}\boldsymbol 1 & * & \cdots & *\end{pmatrix}$$

关注第一列, 得到 $H\boldsymbol 1=\boldsymbol 1$. 所以有

$$HJ=\frac 1n H\boldsymbol 1\boldsymbol 1^T=\frac 1n \boldsymbol 1\boldsymbol 1^T=J,\quad JH=\frac 1n \boldsymbol 1\boldsymbol 1^TH=\frac 1n \boldsymbol 1(H^T\boldsymbol 1)^T=\frac 1n \boldsymbol 1\boldsymbol 1^T=J$$

可以得到

$$(H-J)^2=H^2-HJ-JH+J^2=H-J$$

这样就证明了 $H-J$ 的幂等性. 它的秩为

$$\mathop{\rm rank}(H-J)=\mathop{\rm tr}(H-J)=\mathop{\rm tr}H-\mathop{\rm tr}J=p-1$$

  以上内容可以整理为以下表格:

平方和	二次型表示	二次型的秩
$\rm SSR$	$\boldsymbol y^T(H-J)\boldsymbol y$	$p-1$
$\rm SSE$	$\boldsymbol y^T(I-H)\boldsymbol y$	$n-p$
$\rm SST$	$\boldsymbol y^T(I-J)\boldsymbol y$	$n-1$

2.3 平方和服从的分布

  有如下定理: 设 $\boldsymbol x=(x_1,\dots,x_n)$ 是标准正态分布的一组简单随机样本, $A$ 是一个秩为 $\nu$ 的对称幂等阵, 则有

$$\boldsymbol x^TA\boldsymbol x\sim \chi^2(\nu)$$

这是因为存在正交变换 $\boldsymbol x=Q\boldsymbol y$ 使得 $\boldsymbol y$ 也是相互独立的标准正态样本, 且在该正交变换下原二次型有标准型

$$\boldsymbol x^TA\boldsymbol x=y_1^2+\cdots+y_\nu^2\sim \chi^2(\nu)$$

  据此也可以得到数理统计中的经典结论: 设 $x_1,\dots,x_n$ 是正态变量 $N(\mu,\sigma^2)$ 的简单随机样本, 则有

$$\frac{\sum (x_i-\bar x)^2}{\sigma^2}\sim \chi^2(n-1)$$

因为如果归一化 $y_i:=(x_i-\mu)/\sigma$, 则上式变为

$$\frac{\sum (x_i-\bar x)^2}{\sigma^2}=\sum(y_i-\bar y)^2={\rm SST}\sim \chi^2(n-1)$$

这是我们上面已经导出过的结论.

  平方和除以方差后服从各自自由度的 $\chi^2$ 分布, 即

$$\frac{\rm SSR}{\sigma^2}\sim\chi^2(p-1),\quad \frac{\rm SSE}{\sigma^2}\sim\chi^2(n-p),\quad \frac{\rm SST}{\sigma^2}\sim\chi^2(n-1)$$

均方之比亦服从对应的 $F$ 分布, 例如

$$\frac{\rm MSR}{\rm MSE}=\frac{\rm SSR}{(p-1)\sigma^2}\Big/\frac{\rm SSE}{(n-p)\sigma^2}\sim F(p-1,n-p)$$