UA MATH571A 一元线性回归II 统计推断

在上一篇的例子中，我们讨论到仅使用系数的估计值无法进行稳健的推断。因为系数的估计量服从某个随机分布，给定样本下系数的估计值只是这个随机分布的一个实现，据此无法窥测到系数估计量的整体分布情况。所以本篇从系数估计量的分布出发，试图得出一些稳健的统计推断。

在一元线性回归模型 $Y_i = \beta_0 + \beta_1 X_i + \epsilon_i$Yi​=β0​+β1​Xi​+ϵi​ 中，回归系数$\beta_1$β1​与$\beta_0$β0​的估计量为
$\hat{\beta}_1= \frac{\sum_{i=1}^{N} (X_i -\bar{X}) (Y_i - \bar{Y})}{\sum_{i=1}^{N} (X_i - \bar{X})^2} \\ \hat{\beta}_0 = \bar{Y} - \hat{\beta}_1\bar{X}$β^​1​=∑i=1N​(Xi​−Xˉ)2∑i=1N​(Xi​−Xˉ)(Yi​−Yˉ)​β^​0​=Yˉ−β^​1​Xˉ
注意到
$\sum_{i=1}^{N} (X_i -\bar{X}) (Y_i - \bar{Y}) = \sum_{i=1}^{N} (X_i -\bar{X}) Y_i -\bar{Y} \sum_{i=1}^{N} (X_i -\bar{X}) = \sum_{i=1}^{N} (X_i -\bar{X}) Y_i$i=1∑N​(Xi​−Xˉ)(Yi​−Yˉ)=i=1∑N​(Xi​−Xˉ)Yi​−Yˉi=1∑N​(Xi​−Xˉ)=i=1∑N​(Xi​−Xˉ)Yi​
因此定义
$k_i = \frac{(X_i -\bar{X}) }{\sum_{i=1}^{N} (X_i - \bar{X})^2}$ki​=∑i=1N​(Xi​−Xˉ)2(Xi​−Xˉ)​
$\beta_1$β1​可以写成$Y_i$Yi​的线性组合
$\hat{\beta}_1=\sum_{i=1}^{N} \frac{(X_i -\bar{X}) Y_i }{\sum_{i=1}^{N} (X_i - \bar{X})^2} = \sum_{i=1}^{N} k_i Y_i$β^​1​=i=1∑N​∑i=1N​(Xi​−Xˉ)2(Xi​−Xˉ)Yi​​=i=1∑N​ki​Yi​
同样地，$\beta_0$β0​ 也可以写成$Y_i$Yi​的线性组合
$\hat{\beta}_0 = \bar{Y} - \hat{\beta}_1\bar{X} = \frac{1}{N} \sum_{i=1}^{N} (1- k_i X_i) Y_i$β^​0​=Yˉ−β^​1​Xˉ=N1​i=1∑N​(1−ki​Xi​)Yi​

$\beta_1$β1​的假设检验与置信区间

从$\beta_1$β1​关于$Y_i$Yi​的线性组合出发，
$\hat{\beta}_1= \sum_{i=1}^{N} k_i Y_i = \beta_0 \sum_{i=1}^{N} k_i + \beta_1 \sum_{i=1}^{N} k_i X_i + \sum_{i=1}^{N} k_i \epsilon_i$β^​1​=i=1∑N​ki​Yi​=β0​i=1∑N​ki​+β1​i=1∑N​ki​Xi​+i=1∑N​ki​ϵi​
其中
$\sum_{i=1}^{N} k_i = \sum_{i=1}^{N} \frac{(X_i -\bar{X}) }{\sum_{i=1}^{N} (X_i - \bar{X})^2} = \frac{(N\bar{X}-N\bar{X}) }{\sum_{i=1}^{N} (X_i - \bar{X})^2} =0$i=1∑N​ki​=i=1∑N​∑i=1N​(Xi​−Xˉ)2(Xi​−Xˉ)​=∑i=1N​(Xi​−Xˉ)2(NXˉ−NXˉ)​=0
$\sum_{i=1}^{N} k_i X_i= \sum_{i=1}^{N} \frac{(X_i -\bar{X})X_i }{\sum_{i=1}^{N} (X_i - \bar{X})^2} = \frac{ \sum_{i=1}^{N} (X_i -\bar{X})(X_i -\bar{X}) }{\sum_{i=1}^{N} (X_i - \bar{X})^2} =1$i=1∑N​ki​Xi​=i=1∑N​∑i=1N​(Xi​−Xˉ)2(Xi​−Xˉ)Xi​​=∑i=1N​(Xi​−Xˉ)2∑i=1N​(Xi​−Xˉ)(Xi​−Xˉ)​=1
因此
$\hat{\beta}_1= \beta_1 + \sum_{i=1}^{N} k_i \epsilon_i \sim N(\beta_1, \sum_{i=1}^{N} k_i^2 \sigma^2 )$β^​1​=β1​+i=1∑N​ki​ϵi​∼N(β1​,i=1∑N​ki2​σ2)
从$\hat{\beta}_1$β^​1​的分布中可以发现：1）用最小二乘法与最大似然估计得到的$\beta_1$β1​的估计量都是无偏的；2）估计量服从正态分布。

上一篇我们给出了$\sigma^2$σ2的无偏估计为$MSE$MSE，此处给出简单论证。
$\frac{(N-2)MSE}{\sigma^2} = \sum_{i=1}^N (\frac{Y_i-\hat{Y}_i}{\sigma})^2 \sim \chi^2 (N-2)$σ2(N−2)MSE​=i=1∑N​(σYi​−Y^i​​)2∼χ2(N−2)
简单解释一下，由于$\frac{Y_i-\hat{Y}_i}{\sigma}$σYi​−Y^i​​ 服从标准正态分布，但$\hat{Y}_i$Y^i​使用了两个系数的估计量，所以有两个*度损失，因此总*度为$N-2$N−2，而标准正态分布平方为卡方分布。
$E(MSE) = \sigma^2 \Longleftrightarrow E(\frac{(N-2)MSE}{\sigma^2}) = N-2$E(MSE)=σ2⟺E(σ2(N−2)MSE​)=N−2
因此$\beta_1$β1​的标准差的无偏估计为
$se(\hat{\beta_1})=\hat{\sigma} \{\hat{\beta}_1\} = \sqrt{MSE\sum_{i=1}^{N} k_i^2 }$se(β1​^​)=σ^{β^​1​}=MSEi=1∑N​ki2​​
其中
$\sum_{i=1}^{N} k_i^2 = \frac{\sum_{i=1}^{N} (X_i-\bar{X})^2}{[\sum_{i=1}^{N} (X_i-\bar{X})^2]^2} = \frac{1}{\sum_{i=1}^{N} (X_i-\bar{X})^2}$i=1∑N​ki2​=[∑i=1N​(Xi​−Xˉ)2]2∑i=1N​(Xi​−Xˉ)2​=∑i=1N​(Xi​−Xˉ)21​
从而可以构造t分布
$t = \frac{\hat{\beta}_1 - \beta_1}{se(\hat{\beta_1})} = \frac{\hat{\beta}_1 - \beta_1}{\sqrt{\frac{MSE}{\sum_{i=1}^{N} (X_i-\bar{X})^2}}} \sim t(N-2)$t=se(β1​^​)β^​1​−β1​​=∑i=1N​(Xi​−Xˉ)2MSE​​β^​1​−β1​​∼t(N−2)

Gauss-Markov定理

回归系数的最小二乘估计是最优线性无偏估计（Best Linear Unbiased Estimate， BLUE）。这表示最小二乘估计在所有线性无偏估计中方差最小。假设某估计量
$\tilde{\beta}_1 = \sum_{i=1}^{N} c_i Y_i$β~​1​=i=1∑N​ci​Yi​
是无偏估计，则
$E(\tilde{\beta}_1) = \sum_{i=1}^{N} c_i E(Y_i) = \beta_0 \sum_{i=1}^{N} c_i + \beta_1 \sum_{i=1}^{N} c_i X_i = \beta_1$E(β~​1​)=i=1∑N​ci​E(Yi​)=β0​i=1∑N​ci​+β1​i=1∑N​ci​Xi​=β1​
从而
$\sum_{i=1}^{N} c_i = 0 \\ \sum_{i=1}^{N} c_i X_i = 1$i=1∑N​ci​=0i=1∑N​ci​Xi​=1
不妨假设$c_i = k_i + d_i$ci​=ki​+di​，
$Var(\tilde{\beta}_1) = \sigma^2 \sum_{i=1}^{N} c_i^2 = \sigma^2 \sum_{i=1}^{N} (k_i+d_i)^2 = \sigma^2 \sum_{i=1}^{N} (k_i^2 + d_i^2 + 2k_i d_i)$Var(β~​1​)=σ2i=1∑N​ci2​=σ2i=1∑N​(ki​+di​)2=σ2i=1∑N​(ki2​+di2​+2ki​di​)
其中交叉项为零
$\sum_{i=1}^{N} k_i d_i = \sum_{i=1}^{N} k_i (c_i - k_i) = \frac{ \sum_{i=1}^{N}c_i (X_i - \bar{X}) - 1}{\sum_{i=1}^{N}(X_i-\bar{X})^2}= \frac{( \sum_{i=1}^{N}c_i X_i - \bar{X}\sum_{i=1}^{N}c_i ) - 1}{\sum_{i=1}^{N}(X_i-\bar{X})^2}=0$i=1∑N​ki​di​=i=1∑N​ki​(ci​−ki​)=∑i=1N​(Xi​−Xˉ)2∑i=1N​ci​(Xi​−Xˉ)−1​=∑i=1N​(Xi​−Xˉ)2(∑i=1N​ci​Xi​−Xˉ∑i=1N​ci​)−1​=0
所以方差可以简化为
$Var(\tilde{\beta}_1) = \sigma^2 (\sum_{i=1}^{N} k_i^2 +\sum_{i=1}^{N} d_i^2 ) \ge \sigma^2 \sum_{i=1}^{N} k_i^2$Var(β~​1​)=σ2(i=1∑N​ki2​+i=1∑N​di2​)≥σ2i=1∑N​ki2​
因此最小二乘估计是BLUE。

检验的势

如果我们想要检验$\beta_1$β1​是否等于猜测值$\beta_{10}$β10​，那么可以使用检验方差未知时正态分布均值的t检验方法。用假设检验的语言描述如下：
$H_0: \beta_1 = \beta_{10} \\ H_a: \beta_1 \ne \beta_{10}$H0​:β1​=β10​Ha​:β1​​=β10​
其中原假设$H_0$H0​为$\beta_1$β1​等于猜测值$\beta_{10}$β10​，备择假设$H_a$Ha​为$\beta_1$β1​不等于猜测值$\beta_{10}$β10​。显然原假设和备择假设包含了$\beta_1$β1​所有可能的值，因此原假设与备择假设总是有且仅有一个成立。假设检验有两种可能的错误，弃真、取伪。弃真（第一类错误）指的是拒绝了应该是正确的原假设，取伪（第二类错误）指的是接受了应该是错误的原假设。弃真的概率用$\alpha$α表示，取伪的概率用$\beta$β表示
$\alpha = P(reject\ H_0|H_0\ is\ true) \\ \beta = P(accept\ H_0|H_0\ is\ false)$α=P(reject H0​∣H0​ is true)β=P(accept H0​∣H0​ is false)
检验的势（Power）的含义是当原假设错误时，能准确拒绝原假设的概率，被定义为
$1-\beta =1- P(accept\ H_0|H_0\ is\ false) = P(reject\ H_0|H_0\ is\ false)$1−β=1−P(accept H0​∣H0​ is false)=P(reject H0​∣H0​ is false)
我们希望犯这两类错误的概率都尽可能低，而在给定样本与统计模型时这两个概率总是此消彼长的。因为取伪的后果更加严重，因此在做假设检验时总是在控制$\alpha$α的基础上让$\beta$β尽可能小，基于这种思想，假设检验其实是一个优化问题。

双边检验，单边检验与置信区间

置信区间

基于真实的系数构造的t分布为
$t = \frac{\hat{\beta}_1 - \beta_1}{se(\hat{\beta_1})} = \frac{\hat{\beta}_1 - \beta_1}{\sqrt{\frac{MSE}{\sum_{i=1}^{N} (X_i-\bar{X})^2}}} \sim t(N-2)$t=se(β1​^​)β^​1​−β1​​=∑i=1N​(Xi​−Xˉ)2MSE​​β^​1​−β1​​∼t(N−2)
根据该分布可以给出面的关系式，其中$1-\alpha$1−α是置信水平
$1-\alpha = P(t(\frac{\alpha}{2},N-2)< t<t(1-\frac{\alpha}{2},N-2))$1−α=P(t(2α​,N−2)<t<t(1−2α​,N−2))
从而
$t(\frac{\alpha}{2},N-2)< t<t(1-\frac{\alpha}{2},N-2) \\ t(\frac{\alpha}{2},N-2)< \frac{\hat{\beta}_1 - \beta_1}{se(\hat{\beta_1})} <t(1-\frac{\alpha}{2},N-2) \\ \beta_1+se(\hat{\beta_1})t(\frac{\alpha}{2},N-2)< \hat{\beta}_1 < \beta_1+se(\hat{\beta_1})t(1-\frac{\alpha}{2},N-2) \\ \beta_1-se(\hat{\beta_1})t(1-\frac{\alpha}{2},N-2)< \hat{\beta}_1 < \beta_1+se(\hat{\beta_1})t(1-\frac{\alpha}{2},N-2)$t(2α​,N−2)<t<t(1−2α​,N−2)t(2α​,N−2)<se(β1​^​)β^​1​−β1​​<t(1−2α​,N−2)β1​+se(β1​^​)t(2α​,N−2)<β^​1​<β1​+se(β1​^​)t(1−2α​,N−2)β1​−se(β1​^​)t(1−2α​,N−2)<β^​1​<β1​+se(β1​^​)t(1−2α​,N−2)
上式给出了回归系数估计量的置信水平为$1-\alpha$1−α的置信区间，如果根据根据样本计算得到的回归系数的估计值在置信区间之内，那么我们可以相信这个估计值是合理的，否则我们可以不认可系数的估计值。

双边检验

在一元线性回归中进行的如下检验是双边检验：
$H_0: \beta_1 = 0 \\ H_a: \beta_1 \ne 0$H0​:β1​=0Ha​:β1​​=0
原假设即我们认为不存在X对Y的效应，备择假设的含义是X对Y存在非零的效应。因为回归分析总是想要去验证某种效应是否存在，以及是正向还是负向的效应，而错误拒绝原假设的后果更小并且犯错的概率（$\alpha$α）是被控制在某个之下的。所以假设检验是想看能否拒绝原假设，进而证明某种效应是存在的。在上面的叙述中，我们已经知道了估计量$\hat{\beta}_1$β^​1​服从均值为$\beta_1$β1​的正态分布，因此这个检验其实就是方差未知时对正态分布均值的检验。构造t统计量
$t^* = \frac{\hat{\beta}_1 }{se(\hat{\beta_1})} \sim t(N-2)$t∗=se(β1​^​)β^​1​​∼t(N−2)
$t^*$t∗相当于在原假设下对t的一个猜测值，如果希望将弃真的概率控制为$\alpha$α （检验水平），若
$|t^*| \le t(1-\frac{\alpha}{2},N-2)$∣t∗∣≤t(1−2α​,N−2)
接收原假设，若
$|t^*| > t(1-\frac{\alpha}{2},N-2)$∣t∗∣>t(1−2α​,N−2)
拒绝原假设，接受备择假设。考虑
$|t^*| \le t(1-\frac{\alpha}{2},N-2) \Longleftrightarrow \\ -se(\hat{\beta_1})t(1-\frac{\alpha}{2},N-2)< \hat{\beta}_1 < se(\hat{\beta_1})t(1-\frac{\alpha}{2},N-2)$∣t∗∣≤t(1−2α​,N−2)⟺−se(β1​^​)t(1−2α​,N−2)<β^​1​<se(β1​^​)t(1−2α​,N−2)
上式为原假设的接受域，显然检验水平与置信水平互补时，如果真实系数为0，接受域与置信区间完全一致。检验的p值为如下概率：
$p=P(|t|>t^*)=2P(t>t^*) = 2(1-P(t \le t^*))$p=P(∣t∣>t∗)=2P(t>t∗)=2(1−P(t≤t∗))
所以
$p \le \alpha \Longleftrightarrow |t^*| > t(1-\frac{\alpha}{2},N-2)$p≤α⟺∣t∗∣>t(1−2α​,N−2)

单边检验

在之前的数值例子中，我们想要检验的是年龄对女性肌肉量的效应是否为负，因此我们需要单边检验：
$H_0: \beta_1 \ge 0 \\ H_a: \beta_1 < 0$H0​:β1​≥0Ha​:β1​<0
单边检验如双边检验仅在一些细节上有所不同。若
$t^* \ge -t(1-\frac{\alpha}{2},N-2)$t∗≥−t(1−2α​,N−2)
接收原假设，若
$t^*< -t(1-\frac{\alpha}{2},N-2)$t∗<−t(1−2α​,N−2)
拒绝原假设，接受备择假设。此检验的p值为
$p=P(t<-t^*)$p=P(t<−t∗)
如果检验的是系数是否为正，则.
$H_0: \beta_1 \le 0 \\ H_a: \beta_1 > 0$H0​:β1​≤0Ha​:β1​>0
若
$t^* \le t(1-\frac{\alpha}{2},N-2)$t∗≤t(1−2α​,N−2)
接收原假设，若
$t^*> t(1-\frac{\alpha}{2},N-2)$t∗>t(1−2α​,N−2)
拒绝原假设，接受备择假设。此检验的p值为
$p=P(t>t^*)$p=P(t>t∗)

$\beta_0$β0​的分布

将$\beta_0$β0​关于$Y_i$Yi​的线性组合进一步展开
$\hat{\beta}_0 = \frac{1}{N}\sum_{i=1}^{N} ( 1- k_i X_i) Y_i = \frac{1}{N}\sum_{i=1}^{N} ( 1- k_i X_i) (\beta_0+\beta_1 X_i + \epsilon_i)$β^​0​=N1​i=1∑N​(1−ki​Xi​)Yi​=N1​i=1∑N​(1−ki​Xi​)(β0​+β1​Xi​+ϵi​)
其中
$\beta_0 (1-\sum_{i=1}^{N} k_i X_i) = \beta_0 \\ \beta_1 \sum_{i=1}^{N} ( 1- k_i X_i) X_i=0 \\$β0​(1−i=1∑N​ki​Xi​)=β0​β1​i=1∑N​(1−ki​Xi​)Xi​=0
说明
$\sum_{i=1}^{N} k_i X_i^2 = \sum_{i=1}^{N} \frac{(X_i-\bar{X})X_i^2}{\sum_{i=1}^{N} (X_i-\bar{X})^2}=\frac{\sum_{i=1}^{N} (X_i-\bar{X})X_i \sum_{i=1}^{N} X_i}{\sum_{i=1}^{N} (X_i-\bar{X})^2} = \sum_{i=1}^{N} X_i$i=1∑N​ki​Xi2​=i=1∑N​∑i=1N​(Xi​−Xˉ)2(Xi​−Xˉ)Xi2​​=∑i=1N​(Xi​−Xˉ)2∑i=1N​(Xi​−Xˉ)Xi​∑i=1N​Xi​​=i=1∑N​Xi​
因此
$\hat{\beta}_0 = \beta_0 + \frac{1}{N}\sum_{i=1}^{N} ( 1- k_i X_i) \epsilon_i \\ E(\hat{\beta}_0 ) = \beta_0 \\ Var(\hat{\beta}_0 ) = \sigma^2 \sum_{i=1}^{N} \frac{( 1- k_i X_i) ^2}{N^2} = \sigma^2 (\frac{1}{N}+\sum_{i=1}^{N} k_i^2 \bar{X}^2) \\ \hat{\beta}_0 \sim N(\beta_0, \sigma^2 (\frac{1}{N}+\sum_{i=1}^{N} k_i^2 \bar{X}^2))$β^​0​=β0​+N1​i=1∑N​(1−ki​Xi​)ϵi​E(β^​0​)=β0​Var(β^​0​)=σ2i=1∑N​N2(1−ki​Xi​)2​=σ2(N1​+i=1∑N​ki2​Xˉ2)β^​0​∼N(β0​,σ2(N1​+i=1∑N​ki2​Xˉ2))
知道$\hat{\beta}_0$β^​0​的分布后，可以像对$\hat{\beta}_1$β^​1​做统计推断那样，对$\hat{\beta}_0$β^​0​进行推断。

数值例子：女性肌肉量与年龄的关系

上一篇我们已经建立了女性肌肉量与年龄的一元线性回归模型
$Y_i = \beta_0 + \beta_1 X_i + \epsilon_i$Yi​=β0​+β1​Xi​+ϵi​
其中$Y_i$Yi​表示女性个体的肌肉量，$X_i$Xi​表示女性个体的年龄。现在我们按假设检验的思路对女性个体肌肉量会随着年龄增长而减少的猜想进行验证。
$H_0: \beta_1 \ge 0 \\ H_a: \beta_1 < 0$H0​:β1​≥0Ha​:β1​<0
原假设的含义是女性个体的肌肉量会随着年龄增长而变多或是保持不变，备择假设的含义是女性个体的肌肉量会随着年龄增长变少。从summary()的结果中读取统计量$t^* = \frac{\hat{\beta}_1 }{se(\hat{\beta_1})}$t∗=se(β1​^​)β^​1​​的值：

红框内的结果是$se(\hat{\beta_1})$se(β1​^​)，黄框中的结果是$t^*$t∗，单边检验中$t^*$t∗需要和$t(1-\frac{\alpha}{2},N-2)$t(1−2α​,N−2)比较，假设检验水平为1%

> -qt(1-(.01/2),58)
[1] -2.663287

显然$t^*<-2.663287$t∗<−2.663287，拒绝原假设，接受备择假设：女性个体的肌肉量会随着年龄增长变少。蓝框中的值并非是这个检验的p值，而是双边检验的p值。可以根据上面叙述的结论计算该检验的p值

> pt(-13.19,58)
[1] 2.084381e-19

灰框中是$\beta_0$β0​相关的量，可以用来对$\beta_0$β0​的推断。