大家好！我又出现了(*^__^*) 嘻嘻。刚结束PDE考试（不可避免的凉凉）我就赶紧过来完成了这一篇文章。

这一节我们会结束多元线性回归的内容，并且会努力结束下一个部分——违背基本假设的情况的相关内容。

提供之前的笔记：

  • 回归分析|笔记整理（1）——引入，一元线性回归（上）
  • 回归分析|笔记整理（2）——一元线性回归（下）
  • 回归分析|笔记整理（3）——多元正态分布理论（上）
  • 回归分析|笔记整理（4）——多元正态分布理论（中）
  • 回归分析|笔记整理（5）——多元正态分布理论（下）
  • 回归分析|笔记整理（6）——多元线性回归（上）

我们开始本节的内容。

目录

  • 多元线性回归（下）
    • 偏回归平方和
    • 部分系数显著性检验
  • 违背基本假设的情况
    • 异方差性
      • 异方差性检验
      • 异方差性问题处理方法
        • 一元加权最小二乘估计
        • 多元加权最小二乘估计
    • 自相关性
      • 自相关系数法
      • DW检验
      • 自相关性问题的处理方法
        • 迭代法
        • 差分法
    • 异常值与强影响点
      • 关于自变数 y 异常
        • 删除残差
      • 关于自变数 x 异常，强影响点
        • 强影响点的判定
        • 库克距离
    • Box-Cox变换

多元线性回归（下）

偏回归平方和

这一个概念其实是拟合优度里的内容，单独拉出来说的原因是它在理论上还有很多其余的重要的知识内容。

我们在上一节中说过，偏相关系数是一个不错的衡量回归系数显著性的指标。但是有没有发现我们有说过对方程做显著性检验，也有说过对每个系数做显著性检验，却好像没有说过对部分系数做显著性检验。这就是偏回归平方和可以用到的地方。事实上，它和偏相关系数差别也不大。

首先我们给出它的定义。

Definition 1: Sum of squares of partial regression

Delta SSR_{(j)}=SSR-SSR_{(j)}=SSE_{(j)}-SSE 称为自变数 x_j 的偏回归平方和，其中 SSE_{(j)} 表示在原本p个变数中剔除 x_j 后的剩下p-1个变数的回归出的残差平方和。 SSR_{(j)} 类似定义。

现在我们用它的思想来解决这个问题。

部分系数显著性检验

我们设 y=Xeta+epsilon=egin{bmatrix}X_1 & X_2end{bmatrix}egin{bmatrix}eta_1 \ eta_2end{bmatrix}+epsilon 。那么我们的目的是构造两个平方和的差，这样才有可能分离出部分变数的 SSR,SSE 等，进而进行假设检验。

因为 SSE=Y^T(I_n-H)Y （ H 为帽子矩阵），而这个 H 是由自变数的集合决定的。所以我们如果要构造只包含自变数集合 X_1 的残差平方和，只需要设 H=X_1(X_1^TX_1)^{-1}X_1^T 即可。

为了尊重Prof，我们遵守原来的标记，设 P_X=X(X^TX)^{-1}X^T,P_{X_1}=X_1(X_1^TX_1)^{-1}X_1^T （这个命名在《数值线性代数》里用的比较多，因为它们都是投影(projector)矩阵），那么容易得到

SSE_1-SSE=Y^T(P_X-P_{X_1})Y

回想一下上一节我们怎么进行方程的显著性检验的？对，关键的问题就是要研究统计量矩阵表示中，最中间的那个矩阵。如果我们说明了它是二次型，就可以使用 chi^2 分布相关的知识。所以我们需要看看 P_X-P_{X_1} 到底是个啥。

首先我们要计算 (X^TX)^{-1} ，在第四节中，我们计算过，这里我们直接用那里的公式。

(X^TX)^{-1}=[egin{bmatrix}X_1^T \ X_2^Tend{bmatrix}egin{bmatrix}X_1 & X_2end{bmatrix}]^{-1}=egin{bmatrix} X_1^TX_1 & X_1^TX_2 \ X_2^TX_1 & X_2^TX_2end{bmatrix}^{-1}

=egin{bmatrix}(X_1^TX_1)^{-1} & 0 \ 0 & 0end{bmatrix}+egin{bmatrix}-(X_1^TX_1)X_1^TX_2 \ Iend{bmatrix}Sigma_{22 cdot 1}^{-1} egin{bmatrix}-X_2^TX_1(X_1^TX_1)^{-1} & Iend{bmatrix}

其中 Sigma_{22.1} = X_2^T(I-P_{X_1})X_2

因为计算 P_X 还需要之前加一个 X ，之后加一个 X^T 。所以把 X 按照之前的规则分块，写开乘上 (X^TX)^{-1} 的第一个部分，就是 P_{X_1} ，而第二项乘完相当于 (I-P_{X_1})X_2Sigma_{22cdot1}^{-1}X_2^T(I-P_{X_1}) （注意 P_{X_1} 为对称阵），所以加在一起我们容易得到最后的结果

P_X-P_{X_1}=(I-P_{X_1})X_2[X_2^T(I-P_{X_1})X_2]^{-1}X_2^T(I-P_{X_1})

下面我们试试看，能不能证明这是一个对称幂等阵，如果证明出来了，就说明这是一个二次型，就可以说明它服从一个 chi^2 分布了。

先试著把它平方一下看看

(P_X-P_{X_1})(P_X-P_{X_1})=P_X-P_XP_{X_1}-P_{X_1}P_X+P_{X_1}

所以我们关键要看 P_XP_{X_1},P_{X_1}P_X 分别是什么。事实上，你只需要注意到 P_X-P_{X_1} 它左右都有一个 I-P_{X_1} ，而这个矩阵无论是左乘 P_{X_1} ，还是右乘，都能让矩阵变为0。所以我们事实上就证明了 P_XP_{X_1}=P_{X_1}P_X=P_{X_1} ，也就自然说明了对称幂等。

现在，我们根据第五节的内容就可以得到

frac{SSE_1-SSE}{sigma^2} sim chi^2(r,lambda) ， r=tr(P_X-P_{X_1})=r(X)-r(X_1) ， lambda = frac1{sigma^2}eta^TX^T(P_X-P_{X_1})Xeta= frac1{sigma^2}eta_2^TX_2^T(I-P_{X_1})X_2eta_2 。

将 X 分块就可以得到上面的结果。

使用和上一节一样的分析思路，代入原假设的条件，看看非中心化参数是不是为0。有没有发现在 eta_2=0 的时候，确确实实 SSE_1-SSE 服从一个中心 chi^2 分布，这也一定程度上说明了我们确实通过分离出 X_2 的残差平方和实现了构造分布所需要的统计量。还有一个统计量，只需要根据 (I-P_X)(P_X-P_{X_1})=0 ，就可以得到 SSE_1-SSE 与 SSE 独立。所以这样的话，就可以得到下面的统计量

F=frac{(SSE_1-SSE)/(r(X)-r(X_1))}{SSE/(n-r(X))} sim F(r(X_2),n-r(X))

这是针对 H_0:eta_2=0 的一个检验统计量。有兴趣的同学可以看看，当 eta_2 是一维的时候，它与我们上一节说的t检验是一致的。

我们在这一节的推导中省略了较多的细节，如果发现推导出现困难，那就应该回去多看看之前的理论和思路了。

违背基本假设的情况

这是一个全新的部分，在这里我们会说明，在回归的三大基本假设不满足的情况下，会有什么可能的解决方案。为了防止大家有所遗忘，我这里再一次写一下所谓的Gauss-Markov条件。

Definition 2: Gauss-Markov

egin{cases}E(epsilon_i)=0 \ cov(epsilon_i,epsilon_j) = egin{cases}sigma^2 & i=j \ 0 & i e jend{cases}end{cases}

我们主要关注的就是违背G-M条件情况下，我们应该如何处理。

异方差性

数学上说就是 D(epsilon_i) e D(epsilon_j)(i e j) 。现实中这样的例子也有很多，比方说收入模型，贫穷如我的人整天就会想怎么才能吃饱，而富有的人就会在想双十一要买多少东西，这消费就不可能是一个量级。在异方差出现的时候，会有很多问题。比方说参数不再是最佳线性无偏估计（但依然无偏），显著性检验也失效了。所以回归的效果也很不理想。所以统计学家要想办法去侦测到它，并且努力去消除它。

异方差性检验

因为正常情况下， var(e_i) = sigma^2 ，所以异方差性是可以通过残差看出来的。这就是残差图检验的由来。比方说我的残差长下面这个样子

有没有看到随著观测的进行，它的方差越来越大了？所以这个时候就可以大概的判断出，它是存在异方差性的。但是这样的方法未免过于主观，而统计学家显然更喜欢用检验和p值的方法，所以就有了下面的秩相关检验。请注意，与前六节不同，我们这里提供的检验大部分都只是应用的用途，并不会提供相关的证明。

Definition 3: Spearman test

定义 r_s = 1 - frac{6}{n(n^2-1)}sum_{i=1}^n d_i^2 ，那么斯皮尔曼检验量为 t= frac{sqrt{n-2}}{sqrt{1-r_s^2}}r_s 。

我们要强调的是这个是等级相关系数检验。在做检验之前需要先对模型做一次回归（虽然我们这里已经知道异方差性存在的情况下，回归没啥用了。但是如果你不做回归测试异方差性，你又怎么确定呢？）。得到 {e_i} 。然后取残差的绝对值 |e_i| ，把 x_i 与 |e_i| 都按照从高到低或者从低到高排序，最后标记二者的排位（就是第几大或者第几小），算出二者等级的对应差值计算出来就是 d_i 。比方说一个数据的自变数值 x_i 是第8大的，但是它的对应的残差绝对值 |e_i| 是第3大的，那么对应的 d_i^2=25 。

这个检验量在 n > 8 的时候是近似服从t分布的，因此如果检验量的值 le t_{alpha/2}(n-2) ，就可以认为没有异方差。

异方差问题处理方法

一元加权最小二乘估计(WLS)

这是一种解决异方差性的方法。一般来说，如果我们什么都不管，那么实际上就是要最小化 sum(y_i-eta_0-eta_1x_i)^2 。要注意到的是这个和式的每一项的期望都是 sigma_i^2 （因为异方差性假设存在，所以我们不再使用 sigma^2 ）。所以如果某一项方差越大，实际上这一项所占的比重就很大，那么为了最小化我们的离差平方和，就必须要让回归直线「尽量偏向」这个方差很大的数据点。实际生活中的情况也符合这样的预知。

为了解决这个问题，我们把我们的平方和改一下，写成下面这个样子

Q_omega(eta_0,eta_1)=sum omega_i(y_i-eta_0-eta_1x_i)^2

按照相同的方法进行回归，可得

hat eta_{0omega}=ar y_omega- hat eta_{1omega}ar x_omega

hat eta_{1omega} = frac{sum_{i=1}^{n}(x_i-ar x_omega)(y_i-ar y_w)}{sum_{i=1}^{n}(x_i- ar x_omega)^2}

其中 ar x_omega=frac{1}{sumomega_i}sum_{i=1}^{n}omega_ix_i, ar y_omega = frac1 {sum omega_i}sum_{i=1}^{n}omega_iy_i 。

所以这个回归的关键就是如何选择我们的 omega_i 。直观上来看，因为每一项的期望是 sigma_i^2 ，所以只需要让 omega_i=frac1{sigma_i^2} 即可。是不是就这么简单的结束了？

不好意思，没那么简单。理论可行，可是 sigma_i 是啥你不知道啊。所以如果没有电脑，我们一般是通过残差图去「猜测」应该用什么权。比方说如果 x_i 与 sigma_i 成正比，那么这个时候可以考虑拿 omega_i=frac1{kx_i^2} 去作为权函数。实际上我们也是一般用诸如 frac1{x_i^m} 这样的自变数的幂函数来构造权函数。

当然了，如果你有电脑那就方便多了，直接考虑使用极大似然估计，找到让这个值最大的那个 m 就好。

多元加权最小二乘估计

思路是完全类似的。

我们设 cov(epsilon) = diag(sigma_1^2,sigma_2^2,ldots,sigma_n^2) ，因为在一元的情况下，我们相当于每一个元素都加了一个权 omega_i ，那么在这里，我们也想这么干，那矩阵的意义下就应该两边乘一个对角阵。因此我们设 omega=diag(sqrt{omega_1},sqrt{omega_2},ldots,sqrt{omega_n}) ，那么注意到这个时候，如果原来的方程是 Y=Xeta+epsilon ，那么两边乘上矩阵 omega 就可以得到 ilde Y = ilde X eta+ ilde epsilon ，其中 cov(epsilon)=diag(omega_1sigma_1^2,omega_2sigma_2^2,ldots,omega_nsigma_n^2) 。那么这个时候，对应的估计值

hat eta = ( ilde X^T ilde X)^{-1} ilde X^T ilde Y=(X^Tomega^2X)^{-1}X^Tomega^2Y

注意 omega^2=diag(omega_1,omega_2,ldots,omega_n) ，所以这自然就有办法可以让残差 epsilon 满足G-M条件了。

当然了，多元本身也有一些其余的需要处理的地方，比如说我们的权函数。在一元中我们可以用自变数的幂函数构造。但是多元的情况，如果我们用每一个自变数的幂函数构造，那么对应的计算量可能就是 O(n^2) 级别的，所以一般都只使用其中一个自变数。所以我们用哪一个自变数呢？

这也是有一个法则的，一般来说需要计算每一个自变数与变通残差（ ilde epsilon ）的等级相关系数，取最大的那个构造即可。

自相关性

自相关性的意思就是 cov(epsilon_i,epsilon_j) e 0( i e j) 。这在现实生活中也是很常见的。比如说金融危机一般都是要延后两三年才会有很显著的负面影响。又比如说回归中遗漏了关键变数。时间序列模型本质上也就是一种自相关的模型。

在回归中，因为自相关相当于违背了G-M条件，所以会出现很多问题。比如参数估计值不再是BLUE了，MSE也会严重低估误差项的方差，这就会导致 t 过高，所以F,t检验也就失效了。另外，最小二乘估计量也会对抽样的波动很敏感，意思是说 hat eta 虽然无偏，但是估计出来的值却可能严重歪曲真实情况。

一样，如何检测自相关性？有一种方法是根据时间的顺序画出 e_t,t 的残差关系图，这里我拿PPT上的图片来举例子。

注意这里的残差是根据时间顺序画出来的。可以看出这相当于是说随著时间的推移，残差并不是散乱，而是有序，或者说以一个函数形式出现的。这就说明存在自相关性了。

但是，使用图片总是不能够说服别人。所以我们需要一些更好的办法。

自相关系数法

我们先给出它的定义。

Definition 4: Autocorrelation coefficient

定义误差序列 epsilon_1,epsilon_2,cdots,epsilon_n 的自相关系数为 ho = frac{sum_{t=2}^{n}epsilon_tepsilon_{t-1}}{sqrt{sum_{t=2}^{n}epsilon_t^2}sqrt{sum_{t=2}^{n}{epsilon_{t-1}^2}}}

这也是时间序列中一个很重要的统计量。和简单相关系数对比容易得到它的范围是 [-1,1] 。

因为误差是未知的，所以一般用残差 {e_i} 代替。这样的到的 hat ho 因为和样本量有关，所以需要定义其余的量来处理它。

DW检验

这是一种检验自相关性的方法，但是它只能检验随机扰动项具有一阶自回归形式的序列相关问题。也就是说，它的误差满足下面的式子

epsilon_t= ho epsilon_{t-1}+u_t

其中 u_t 服从G-M条件。

并且定义DW统计量为

DW=frac{sum_{t=2}^{n}(e_t-e_{t-1})^2}{sum_{t=2}^{n}e_t^2}

显然，这里序列相关性的原假设为 H_0: ho=0 。现在问题来了，DW值可以取到哪些呢？

为了解决这个问题，我们只需要把它的分子展开一下。

frac{sum_{t=2}^{n}(e_t-e_{t-1})^2}{sum_{t=2}^{n}e_t^2}=frac{sum_{t=2}^{n}e_t^2+sum_{t=2}^{n}e_{t-1}^{2}-2sum_{t=2}^ne_te_{t-1}}{sum_{t=2}^{n}e_t^2}

分子的第一项和第二项在 n 比较大的时候是几乎相同的（所以一般来说DW检验要求 n>15 ），而第三项与分母的比就是我们的 hat ho 。所以这个值其实大约为 2(1-hat ho)，所以这样的话，我们大概可以推出它的范围为 [0,4] 。

一般来说，我们会根据样本数 n 和解释变数的数目 k ，查DW分布表来确定两个值 d_L,d_U 再将 [0,4] 分块，大家可以在网上查阅更多相关细节，这里我们略去。

自相关问题的处理方法

其实我们可以通过修改回归模型，增加自变数等方法解决相关问题。但是如果这二者都不行的话，就只能考虑下面的方法了。

迭代法

还是要以误差项存在一阶自相关来举例。并且我们假设模型为 y_t=eta_0+eta_1x_t+epsilon_t 。

根据这个模型，让时间倒退一秒，就可以得到 y_{t-1}=eta_0+eta_1x_{t-1}+epsilon_{t-1} 。所以为了消除自相关性，归根到底就是要让误差项回归到 u_t ，这就需要我们得到 epsilon_t- ho epsilon_{t-1} 。所以我们来计算 y_t- ho y_{t-1}

y_t- ho y_{t-1}=(eta_0- hoeta_0)+eta_1(x_t- ho x_{t-1})+(epsilon_t- ho epsilon_{t-1})

对应的变数做换元就能得到 y_t=eta_0+eta_1x_t+u_t 。这个时候可以看出误差项就满足G-M条件了。

那么这样的方法可以看出如果真的误差项存在一阶自相关的话，那么很明显是有效的。但是实际情况并不总是如此，所以我们的方法是不停的迭代，直到我们的DW检验能够说明它没有自相关了为止，可以说是简单粗暴啊。

差分法

差分法的适用范围就更窄了，它是适用于原模型存在较高程度一阶自相关的情况才可使用。在迭代法的模型中我们设 ho=1 ，就可以得到一个差分法的模型

Delta{y_t}=eta_1Delta x_t+u_t

对它做一个回归可以得到 hat eta_1=frac{sum_{t=2}^{n}Delta y_tDelta x_t}{sum_{t=2}^n Delta x_t^2} （注意从 t=2 开始的原因是差分肯定只能从第二项开始才会有数据）。

因为差分法更简单，且迭代法需要用样本估计自相关系数，而它的估计误差又会影响迭代法的效率，所以如果自相关系数很高的时候，一般优先考虑差分法。

异常值与强影响点

有的时候数据会有一些很诡异的点，这个时候也会导致对应的数据点产生较大的残差，同样是会让回归直线「尽量偏向」这些异常点，自然也会带来很大的麻烦。

一般来说我们会分 x,y 两个维度讨论异常值。

关于因变数 y 异常

在数据分析中，刚开始总是要看有没有特别特别「高」的点。一般来说会认为超过 pm 3hat sigma 的残差的话它就是异常值。但是问题在于，多元回归中 D(e_i)=(1-h_{ii})sigma^2 ，其中 h_{ii} 为帽子矩阵 H=X(X^TX)^{-1}X^T 的主对角线元素，这也就说明每一个数据点的误差是不相同的。那么单纯的因为它「特别高」就认为数据异常就不合适了。因为这很有可能是残差导致的，换句话说这个数据「特别高」不是因为它异常，而是因为它「就完全有可能这么高」。换句话说，因为误差是每一个数据点的固有性质，所以如果是因为残差特别大，导致某一个数据点像异常值，那么即使你剔除掉这个异常值，也不会对回归有任何帮助。

那么应该如何去做呢？我们在之前介绍过一个学生化残差 SRE_i = frac{e_i}{hat sigma sqrt{1-h_{ii}}} 。看似通过把杠杆值的影响去除掉可以解决方差不等的问题，但是如果观测数据中真的存在异常值，学生化残差也没有什么卵用。这是因为这个时候，异常值的存在会使得回归线「偏向」它，进而使得回归的标准差 hat sigma 实际上是偏大的。那么这样子的话，实际的学生化残差是偏小的，这就不能使用 |SRE_i|>3 的原则来判断残差了。

为了解决异常值的问题，我们需要别的办法。

删除残差

我们这么构造删除残差：针对第 i 个观测值，我们计算它的残差时，用其余 n-1 个观测值拟合回归方程，计算出第 i 个观测值的删除拟合值 hat y_{(i)} ，那么这个值显然不会受到第 i 个值是否是异常值的影响。所以我们定义删除残差为

e_{(i)}=y_i-hat y_{(i)}

下面我们回到理论的模型下观察下这样的残差应该如何使用和检验。

如果我们设原模型为 y=Xeta+epsilon ，那么对应的估计为 hat eta=(X^TX)^{-1}X^Ty, hat sigma^2=frac{SSE}{n-p-1} ，并且残差为 e=y-hat y ，标准化残差为 frac{e_i}{hat sigmasqrt{1-h_{ii}}} 。现在我们删除第 i 个观察值，让模型变为

Y_{(i)}=X_{(i)}eta+epsilon_{(i)}

那么对应的估计即为 hat eta_{(i)}=(X^T_{(i)}X_{(i)})^{-1}X_{(i)}^TY_{(i)},hat sigma^2=frac{SSE_{(i)}}{n-p-2} 。下面注意到

(X^T_{(i)}X_{(i)})^{-1}=(X^TX-x_ix_i^T)^{-1}

=(X^TX)^{-1}+frac{(X^TX)^{-1}x_ix_i^T(X^TX)^{-1}}{1-x_i^T(X^TX)^{-1}x_i} （这一步我也不知道怎么来的，Prof和我们说不要求证明……） =(X^TX)^{-1}+frac{(X^TX)^{-1}x_ix_i^T(X^TX)^{-1}}{1-h_{ii}}

且有

X^T_{(i)}Y_{(i)}=X^TY-x_iy_i （这可以通过矩阵分块看出来，注意设计矩阵按行分块就是每一个样本的数据，按列分块就是每一个自变数对应的数据）

所以我们可以算出

hat eta_{(i)}=[(X^TX)^{-1}+frac{(X^TX)^{-1}x_ix_i^T(X^TX)^{-1}}{1-h_{ii}}](X^TY-x_iy_i)

=hat eta+frac{(X^TX)^{-1}x_ix_i^That eta}{1-h_{ii}}-(X^TX)^{-1}x_iy_i-frac{(X^TX)^{-1}x_ix_i^T(X^TX)^{-1}x_iy_i}{1-h_{ii}}

注意到最后一项中间有一个 x_i^T(X^TX)^{-1}x_i=h_{ii} ，而剩下的部分是 (X^TX)^{-1}x_iy_i 又正好与第三项重合。所以化简一下会有

hat eta_{(i)}= hat eta+frac{(X^TX)^{-1}x_ix_i^That eta}{1-h_{ii}}-frac{(X^TX)^{-1}x_iy_i}{1-h_{ii}} =hat eta-frac{(X^TX)^{-1}x_ie_i}{1-h_{ii}}

你应该还记得我们的残差是与 SSE 有关系的，所以我们需要弄到 SSE_{(i)} 与 SSE 的关系，这样才能构造出全新的残差表达式。那么注意到

SSE_{(i)}=|Y_{(i)}-X_{(i)}hat eta_{(i)}|^2

直观意义上，这相当于计算了「去掉了第 i 个观测值之后的向量与去掉了第 i 个预测值之差的2-范数」。那么如果我把 Y_{(i)} 改为 Y ， X_{(i)} 改为 X ，就相当于增加了「第 i 个观测值与第 i 个通过未使用第 i 个观测值的信息得到的预测值的差」。这可能有些绕，但是也许你通过下面的结果就可以帮助你理解。

|Y_{(i)}-X_{(i)}hat eta_{(i)}|^2=|Y-X hat eta_{(i)}|^2-(y_i-x_i^That eta_{(i)})^2

=|Y-Xhat eta+X(hat eta- hat eta_{(i)})|^2-(y_i-x_i^That eta_{(i)})^2 （这是为了弄出 SSE ）

注意到

(Y-Xhat eta)^{T}X(hat eta-hat eta_{(i)})=(Y^T-hat eta^TX^T)frac{X(X^TX)^{-1}x_ie_i}{1-h_{ii}}

=(Y^TX(X^TX)^{-1}-hat eta^T)x_ie_i/(1-h_{ii})=0

所以两个向量正交，平方和就是和的平方，可以得到

SSE_{(i)}=SSE+|X(hat eta-hat eta_{(i)}|^2-(y_i-x_i^That eta_{(i)})^2

简单计算容易得到

|X(hat eta-hat eta_{(i)})|^2=frac{h_{ii}e_i^2}{(1-h_{ii})^2}

y_i-x_i^T hat eta_{(i)}=y_i-x_i^That eta+x_i^T(hat eta-hat eta_{(i)})=frac{e_i}{1-h_{ii}}

所以最终可以得到残差平方和

SSE_{(i)}=SSE-frac{e_i^2}{1-h_{ii}}

有了残差平方和，我们就可以考虑计算我们的删除学生化残差了。学生化残差是使用残差和对应的 hat sigma^2 来构造的统计量（相关内容可以查看第2节一元或者第6节多元的情况）。所以删除学生化残差就自然是使用预测残差和对应的 hat sigma_{(i)}^2 了。而因为 hat sigma_{(i)}^2 和 SSE_{(i)} 有关，上面已经给出了表达式，所以不难得到

hat sigma_{(i)}^2=frac{n-p-1-SRE_{i}^2}{n-p-2}hat sigma^2

下面自然是考虑预测残差的方差，进而构造分布。而预测残差我们上面也已经给出表达式，所以有

D(y_i-x_i^Thateta_{(i)})=frac{sigma^2}{1-h_{ii}}

因为我们是使用的部分数据（删除了第 i 个观测值），所以自然是考虑使用 hat sigma_{(i)} 作为 sigma 的估计（因为虽然去掉了一个数据，但是注意到G-M条件是不变的，所以对应的误差的方差依然是 sigma^2 ）。这样就可以得到最终的结果

SRE_{(i)}=frac{y_i-x_i^That eta_{(i)}}{hat sigma_{(i)}/sqrt{1-h_{ii}}}=frac{e_i}{sqrt{1-h_{ii}}hat sigma_{(i)}}

总结一下这些漫长的证明，就是下面的结果

Proposition 1:

e_{(i)}=frac{e_i}{1-h_{ii}} SRE_{(i)}=SRE_i(frac{n-p-2}{n-p-1-SRE_i^2})^{1/2}

一般来说，认为 |SRE_{(i)}|>3 的时候就是异常值点。

关于自变数 x 的异常值，强影响点

这里主要涉及强影响点的概念。还是关于残差的方差式 D(e_i)=sigma^2(1-h_{ii}) ，可以看出 h_{ii} 大的点残差小，因此如果观测值的杠杆值（ h_{ii} ）大，就会使得回归方程偏移产生影响。所以一般来说杠杆值大的点我们叫做强影响点，注意它不一定是 y 的异常值。

实际上，强影响点还是很需要被关注的，这是因为有可能它会暗示我们回归方程不再是线性的。同时也会使得回归方程「偏向」自身。所以不可避免的会产生一些影响。这就使得检测强影响点变得很有必要。

不过另一方面，如果它不是 y 的异常值，那么就不一定会说对回归产生影响。所以归根到底我们还是需要判断它是不是 y 的异常值。所以总结一下就是两步：判断强影响点，判断是否是 y 的异常值。

强影响点的判定

直觉上很明显，超过均值太多肯定就是不太正常的。根据这个思路，我们设 ch_i=h_{ii}-frac1n ，那么这样的话注意到 tr(H)=r(H)=p+1=sum_{i=1}^{n}h_{ii} ，就可以得到 sum_{i=1}^{n}ch_i=sum_{i=1}^{n}h_{ii}-1=p 。所以均值就是 overline{ch}=frac p n 。实际情况下，如果某一个自变数对应的杠杆值 h_{ii} >2overline{ch} ，就认为它是杠杆值。

库克距离

这也是一个较为常用的判定的公式

D_i=frac{e_i^2}{(p+1)hat sigma^2}cdot frac{h_{ii}}{(1-h_{ii})^2}

实际情况是很复杂的，所以一般使用一个粗略的标准，认为 D_i > 1 就是异常值， D_i<0.5 就是非异常值。

Box-Cox变换

也叫方差稳定性变换。这是我们涉及的最后一个消除异方差的变换。

y_i^{(lambda)}=egin{cases}frac{(y_i+a)^lambda-1}{lambda} & lambda e 0 \ log(y_i+a) & lambda=0end{cases}

因为要求 y_i+a>0 ，所以这里我们要取 a>-min_i{y_i} 。变换根据 lambda 的不同而形式不同。

这个变换几乎就是万能的了。比方说在 lambda=0 的时候，其实数据就可以很好的消除异常值。比方说 y 很大，但是在经过对数变换后，只要你不是高的很厉害，基本上都会让所有的 y 保持差不多的高度。

当然了，最终的目的自然是要 Y^{(lambda)} sim N_n(Xeta,sigma^2I_n) 。而实际情况下，我们也会取 lambda 使得似然函数的值最大。

注意到似然函数为

f(Y|lambda,eta,sigma^2)=(2pisigma^2)^{-n/2}|J|exp{-frac1{2sigma^2}(Y^{(lambda)}-Xeta)^T(Y^{(lambda)}-Xeta)}

其中 J=prod_{i=1}^{n}y_i^{lambda-1} 为变换对应的Jacobian行列式。代入最小二乘回归的结果可以得到

f(Y|lambda,hat eta_lambda,hat sigma_lambda^2)=(2pifrac{SSE_lambda}{n})^{-n/2}|J|exp{-frac n 2}

取对数即可。剩下的事情交给计算机吧。

小结

本节的内容量很大，我们主要是结束了多元线性回归的一些显著性检验的内容，并且相对简要的介绍了一下现实中常出现的基本假设不满足的情况会有什么常用的解决方案。注意到这里的很多内容我都没有贴出理论的证明，是因为它们一般都是用在计算机上的，重要性也没有一元，多元回归那些大，所以我们就只是把结论贴在那里，主要的操作还是交给计算机去做的。

当然，只要涉及到理论和矩阵的东西，一般都是比较复杂的，这可能也是统计应该有的面目吧。

——————————————————广告——————————————————

本专栏为我的个人专栏，也是我学习笔记的主要生产地。任何笔记都具有著作权，不可随意转载和剽窃。

个人微信公众号：cha-diary，你可以通过它来有效的快速的获得最新文章更新的通知。

本人最近在寻找与数据科学，计算数学，统计有关的科研和实习机会。希望各路大佬神仙如果有看得上我的可以和我联系下~谢谢你们！

专栏目录：笔记专栏|目录

想要更多方面的知识分享吗？欢迎关注专栏：一个大学生的日常笔记。我鼓励和我相似的同志们投稿于此，增加专栏的多元性，让更多相似的求知者受益~


推荐阅读：
相关文章