Transcript 模型误差独立性诊断

第三章
回归诊断
在前两章讨论一元与多元线性回归问题中，我们作了如下
一些假定：

回归函数的线性假设： E (Y ) 是 x1 , x2 ,  xt 的线性函数；

误差的独立性假设：随机误差  1 ,  2 ,  n 相互独立。

误差方差的齐性假设： E i

误差的正态性假设：  1 ,  2 ,  n 服从正态分布。
 0,Var i   2 ；
• 在实际中这些假定是否合理？如果实际数据与这些假设偏离比较大，那
么前面讨论的有关参数的区间估计，假设检验就不一定成立。如果经过
分析，已经确认对所研究的具体数据，上面的假设不成立，那么我们又
希望探讨对数据作怎样的修正后，能使它们满足或近似满足这些假设。
这些就是回归诊断中所要解决的第一个问题。
• 回归诊断的另一个研究的问题是对数据的诊断，探查对统计推断有较大
影响的样本点，这样的点称为强影响点;以及探查与其它数据不是来自
同一模型的样本点，这样的点称为异常点。
§3.1 残差及残差图
由第二章的讨论可知：

Ee  0,Var (e)   2 ( I  H )

Cov(Yˆ , e)  0

当  ~ N (0,  I n ) 时， e ~ N (0,  ( I  H ))
2
2
我们可以看出普通残差与 H 之间有着密切的关系，为
此进一步讨论 H  (hij ) 的性质如下：
定理 3.1 (1) 0  hii  1 ,且 hii  1 时， hij  0, j  i 。
n
(2)
h
i 1
ii
 t  1 , 其中 t 为自变量个数。
证明：（1）因为 H   H , H  H ，所以有
2
n
hii   hij2  hii2   hij2   hij2  0
j 1
j i
j i
故有： hii  hii2 ，由此可得。
n
（2）
1
1



h

tr
(
H
)

tr
(
X
(
X
X
)
X
)

tr
((
X
X
)
X X )  t  1
 ii
i 1
中心化模型下：
hii 
1
 ( xi  x)' L1 ( xi  x) ,
n
1 ( xi  x )2
对一元线性回归方程来说： hii  
.
n
lxx
通常称 hii 很大的点为高杠杆点。所谓 hii 很大，即
hii 很靠近于 1。但 hii 大到多少所对应的点才是高杠杆
点呢？一般说来很难给出一个处处适用的标准，
Sheather(2009)给出高杠杆点的标准是大于 2(t+1)/n,
其中 t 为自变量个数，n 为样本容量。高杠杆点对回归
拟合的影响很大(可参阅§8.1)，在回归诊断中应予以
特别注意。
一般来讲，e1 , e2 ,
, en 是相关的，且它们的方差不等。所以不宜直接用 ei
作比较来判别异常点(第八章会详细介绍)。
引进标准化残差
称：
ei
ri 
1  hii s
(i  1,2, , n)
为标准化残差，其中： hii 为 H 矩阵的第 i 个对角元， s

SSE
。
n  t 1
一般说来， ri 的确切分布由于 ei 与 s 是不独立的所以很
难求得，（ ri 的确切分布可见[4]）。但在模型假设成立时，
r1 , r2 ,, rn 近似独立，且近似服从 N (0,1) ，可以近似认为
r1 , r2 ,, rn 是来自 N (0,1) 的随机子样。依据标准化残差
r1 , r2 ,, rn 近似服从 N (0,1) 且近似相互独立这一结论，常
用残差图对模型假设的合理性进行检验。
残差图是一种直观的工具，它是以标准化残差 r 或
普通残差 e 为纵坐标，以其它的量为横坐标的散点图。
常用的横坐标有如下三种：

以拟合值 yˆ 为横坐标；

以 xi 为横坐标；

以观测时间或序号为横坐标。
一般残差图均要求 n 个点的散布是无规则的(若是标准
化残差，则还要求 n 个点大致落在 | r | 2 的水平带内)。
当残差图中的点呈现某种规律或趋向时，就可以对模型
的假设提出怀疑。利用残差图上点的散布规律作诊断的
方法是回归分析中对模型的诊断的最有效的方法之一，
下面在各节中将详述。
综合以上所述回归诊断有如下主要
内容：
• 识别、判定和检验异常点。
• 区分出对统计推断影响特别大的点（影响分析）。
• 残差分析和残差图能用于研究既定模型与实际数据是否能
很好拟合。其中包括：模型线性诊断、模型误差方差齐性
诊断、模型误差独立性诊断、模型误差正态性诊断等。
§3.2 回归诊断一（数据的诊
断）
(一)、数据诊断的两个基本概念
（1）异常点
在回归模型中，异常点是指对既定模型
偏离很大的数据点。但究竟偏离达到何种
程度才算是异常，这就必须对模型误差项
的分布有一定的假设（通常假定为正态分
布）。
目前对异常点有以下两种较为流行的看法：
把异常点看成是那些与数据集的主体明显不协
调，使得研究者大感惊讶的数据点。这时，异常
点可解释为所假定的分布中的极端点，即落在分
布的单侧或双侧分位点以外的点，而通常取很小
的值（如：0.05 ），致使观察者对数据中出现如
此极端的点感到意外。
（2）强影响点
数据集中的强影响点是指那些对统计
量的取值有非常大的影响力的点。在考虑
强影响点时，有几个基本问题需要考虑：
首先必须明确“是对哪个统计量的影响？”
例如，线性回归模型所考虑的是对回归系
数的估计量的影响；不是对误差方差的估
计影响；或是对拟合优度统计量的影响等
等。分析目标不同，所考虑的影响亦有所
不同。
• 其次，必须确定“度量影响的尺度是什么？”为了定量
地刻划影响的大小，迄今为止已提出多种尺度，基于置
信域的尺度，基于似然函数的尺度等等。在每一种类型
中又可能有不同的统计量，例如基于影响函数就已提出
多种“距离”来度量影响，有Cook距离、Welsch Kuh距离、Welsch距离等等。每一种度量都是着眼于某
一方面的影响，并在某种具体场合下较为有效。这一方
面反映了度量影响问题的复杂性，另一方面也说明了影
响分析的研究在统计诊断中是一个甚为活跃的方向，还
有大量有待解决的问题。
• 强影响点通常是数据集中更为重要的数据点，它往往能
提供比一般数据点更多的信息，因此需引起特别注意。
• 强影响点和异常点是两个不同的概念，它们之间既有联
系也有区别。强影响点可能同时又是异常点也可能不是；
反之，异常点可能同时又是强影响点也可能不是。
(二)、影响分析
设有模型：
 y i   0   1 xi1     t xit   t i  1,2, n

2
 i ~ N (0,  ), 独立同分布
记： IFi  ˆ (i)  ˆ ，其中： ˆ 为模型的最小二乘估计， ˆ (i) 为剔除
第 i 组数据后由剩下的（n-1)组数据求得的最小二乘估计。 IFi 为度量
数据 xi1 , xi 2 ,, xit , yi 
为经验影响函数。
对回归估计影响大小的一种重要统计量，称
( ˆ (i)  ˆ )M ( ˆ (i)  ˆ )
记： Di ( M , c) 
c
Di (M , c) 越大表明了第 i 组数据被剔除后， ˆ 移动的距离越
大。 Di (M , c) 度量了第 i 组数据对 ˆ 影响的大小。一般称
Di (M , c) 为 Cook 距离。
为了给出 Di (M , c) 进一步的表达式，记：
 1 x11

1 x21

X


 1 xn1
x1t   x1 
  
x2t   x2 
ˆ
,
  
  
xnt   xn 
n
n
j 1
j 1
X X   x j x j ， X Y   x j y j ；
剔除第 i 组数据后由剩下的（n-1)组数据得到的相应的系数矩阵及常数项矩阵
为，
n
X (i) X (i)   x j x j  xi xi  X X  xi xi ，
j 1
n
X (i)Y (i)   x j y j  xi yi  X Y  xi yi ；
j 1

 (i )   X (i) X (i)  X (i)Y (i)  X X  xi xi
1

1
 X Y  xi yi 
1
1





X
X
x
x
X
X



1
i i
  X Y  xi yi 
  X X  
1

1  xi  X X  xi 
X X  i

ˆ

1  hii
1
x xˆ
i
  X X  xi yi
1
X X 


1
X X  xi  yˆi  yi  ˆ  X X  xi ei

ˆ


1  hii
1  hii
1
1
xi xi  X X  xi yi
1  hii
1

ˆ  X X 
IFi  ˆ (i)  
1
xi ei
1  hii
进一步，可以把 Cook 距离表示为：
( ˆ (i)  ˆ )M ( ˆ (i)  ˆ )
Di ( M , c) 
c
  X X 1 xi ei  M
 
 
1  hii

 c
  X X 1 xi ei 
 

1  hii


1
1
2  




x
X
X
M
X
X
xi 



ei
s
i 


  
1  hii

 1  hii  s  c 
2
2
s
 ri 2   Pi ( M )
c
进一步，可以把 Cook 距离表示为：
s2
Di ( M , c)  ri
Pi ( M )
c
2

ri2 是第 i 组数据的学生化残差，表示在 xi1 , xi 2 ,, xit , yi 
处拟合的好坏；

s2
是与 i 无关的常数；
c

Pi (M ) 刻划了点 xi1 , xi 2 ,, xit  在自变量空间 R t 中的位
置。
M，c 常用的选择： M  X X , c  (t  1)s 2 ，此时，有：
hii
1
Di ( M , c)  ri
t  1 1  hii 
2
其中： hii 为矩阵 H  X ( X X ) X  的对角线第 i 个元素。
1

Di (M , c) 相对较大的对应点称为强影响点。

hii 较大的点称为高杠杆点。

| ri |较大(大于 2)的点称为异常点。
§3.3 回归诊断二（模型的诊断）
一、
回归函数线性的诊断
（1）模型诊断
诊断回归函数是否是 x1 , x2 ,  xt 的线性函数的主要工具
是以拟合值 yˆ 为横坐标的残差图。t=1 时可用散点图。
例 3.1
I
x
y
yˆ
eˆ
1
80
0.6
1.6625
-1.0625
2
220
6.7
7.75
-1.05
3
140
5.3
4.27143
1.02857
4
120
4
3.40179
0.59821
5
180
6.55
6.01071
0.53929
6
100
2.15
2.53214
-0.38214
7
200
6.6
6.88036
-0.28036
8
160
5.75
5.14107
0.60893
由上述数据，可得 y 关于 x 的一元线性回归方程
yˆ  1.82  0.0435x
回归方程所对应的  的估计值及 x,y 的样本相关系数分别为：
s  0.869
r  0.935
其拟合值及残差见上表。其散点图(a)及残差图(b)如下

从散点图可以看出，x 与 y 之间的关系用线性函数去描述不太
适合，最好看成是 x 的二次函数；

从残差图发现，点的散布是有规律的， yˆ 小的及 yˆ 大的残差
为负，而 yˆ 介于中间的点的残差是正的。
从上述的分析，我们有理由怀疑回归函数线性这一假定是不成
立的。
（2）模型修正
为了修改模型，我们再作以xt为横坐标的残差图
从上图可见， x 较小或较大时，残差小于 0， x 介于中间值时，残
2
差大于 0，从而设想改变回归模型，建立 y 关于 x 、 x 的回归方
程。其基本模型为
yi   0  1 xi   2 xi2   i ,
i  1,2,, n
经计算可得回归方程如下：
yˆ  10.028 0.16424x  0.0040205x 2
其拟合值和残差见表 3.3.2，残差图见图 3.3.3。从残差图可见点的
分布已无规律，说明用二次回归方程是合适的。
模型修改后的预测值及残差
I
yˆ
eˆ
1
0.53542
0.06458
2
6.62292
0.07708
3
5.07649
0.22351
4
3.88482
0.11518
5
6.49375
0.05625
6
2.37113
-0.22113
7
6.71935
-0.11935
8
5.94613
-0.19613
模型修改后的残差图
二、 误差方差齐性诊断
对多元线性回归模型：

 y i   0   1 xi1     t xit   t i  1,2, n

2

 i ~ N (0,  i ),并且独立
要诊断  i2 是否相等。即要检验：
H0 : 12   22 
  n2 。
（1）模型诊断
常用的对误差方差的诊断的方法有：
1． 残差图
误差方差非齐性时，以拟合值 yˆ 为横坐标的残差图一般是
呈现喇叭型或倒喇叭型，
1． 若有重复试验， 可采用如下三种检验方法：

Hartley 检验（最大 F 比检验）要求有相等的重复试验三次或三次以上，统
计量为：
Fmax
max(s12 , s22 , sn2 )
2

,
其中
:
s
i 为组内方差
2
2
2
min(s1 , s2 , sn )
拒绝域为：
W  Fmax  Fmax, (n, m 1)，
其中：n 表示不同条件数，m 表示重复次数。

Cochran 检验（最大方差检验）要求有相等的重复试验二次或二次以
上，统计量为
Gmax 
max(s12 , s 22 , s n2 )
n

i 1
, 其中 : si2为组内方差
si2
拒绝域为：
W  Gmax  Gmax, (n, m 1)，
其中：n 表示不同条件数，m 表示重复次数。

Bartlett 检验（  2 检验）只要有重复就可以，这是三个检验中要求条件最弱的，
其统计量为：
1

 2   f e ln s 2   (mi  1) ln si2 
c


i
1
1

i m  1 f
i
e
其中： f e   (mi  1) ， s 2   (mi  1)si2 / f e ， c 
1
3(n  1)
i
i
当 H0 :  12   22     n2 为真时，上述统计量近似服从自由度为 n  1 的  分
2
布，对给定的显著水平  ，其拒绝域为


W   2  1 (n 1)
（1）
模型修正
当模型经检验表明方差非齐性时，常有两种方法可修改模型。
方法一：加权最小二乘估计
若我们可以得到 Var (Y )   G 中 G 的表达式（此时， G 为对
2
角矩阵，对角线上的元素不相等）
，那么，可采用加权最小二乘估计
去求模型的参数估计。
方法二、方差稳定性变换
设 Y 是 一 个 随 机 变 量 ， E (Y )   ， 这 里  可 能 是 自 变 量
x1 , x2 ,  xt 的函数，又假定 Var (Y )  g ( ) ，其中 g 是一个已知函
数。设
Z  f (Y )
使 Var (Z )  c （其中 c 为常数）
。下面我们要找出 f 表达式。
记 z  f ( y ) ，并在 y   处作泰勒展开，取近似式为：
f ( y)  f ( )  f ( )( y   )
将 y 改为随机变量 Y 则有：
Z  f (Y )  f ( )  f ( )(Y   )
其方差
D(Z )  ( f ( ))2 g ( )
由 Var (Z )  c 知， ( f ( )) g ( )  c ，有
2
f ()  c / g() (这里取算术平方根)
我们得到方差稳定性变换
f ( y )   c / g ( y )dy
常见的方差稳定性变换(忽略常数)
N
g ( )
变换
1
g (  )  a
Z Y
2
g (  )  a (1   )
Z  arcsin Y
3
g ( )  a 2
Z  ln Y
4
g ( )  a 3
Z  1/ Y
5
g ( )  a
Z  Y 1
4
三、误差的独立性诊断
• 在不少有关时间序列问题中，观测值往往
呈相关的趋势。如河流的水位总有一个变
化过程，当一场暴雨使河流水位上涨后往
往需要几天才能使水位降低，因而当我们
逐日测定河流最高水位时，相邻两天的观
测间就不一定独立。
我们考虑相邻观测间存在的一种最简单的相关情况——一阶
自相关。其基本模型为：
 yi   0  1 xi1 

 i 1   i  ui 1
u , iid ~ N (0,1)
 i
  i xit   i
i  1, 2,
i  1, 2,
, n 1
当   0 时，  i ,  2 ,,  n 之间存在一阶自相关。
这时检验误差的独立性就变成了检验假设：
H 0 :   0  H1 :   0
,n
（1）
模型诊断
1． 残差图
对模型独立性的诊断的残差图一般用以序号为横坐标的残差
图。当残差符号的改变非常频繁，大致有正负相间的趋势时，这时
序列是负相关的   0 ;相反，当残差符号具有“集团”性，即有一
段是下号，另一段是负号，然后又有一段全是下号，这时序列是正
相关的   0 。
 0
 0
判断结果？
判断结果？
2.
DW检验
对假设 H 0 :   0  H1 :   0
Durbin与Watson提出了DW检验，其检验用
的统计量为：
n
DW 
2
(
e

e
)
 i i 1
i 2
n
e
i 1
2
i
若记 r 为序列 e1 , e2 ,, en1 与序列 e2 , e3 ,, en 的样本相关系数(把 r 看成  的
一个估计)，则可以证明(证明细节可见教材)：
DW  2  2r
当 DW  2 过大时，拒绝假设 H 0 :   0  H1 :   0 ，据DW的值
可按下面规则下结论：

DW  d L ，认为  i ,  2 ,,  n 之间存在正相关；

dU  DW  4  dU ，认为  i ,  2 ,,  n 之间不相关；

DW  4  d L ，认为  i ,  2 ,,  n 之间存在负相关；

d L  DW  dU 或 4  dU  DW  4  d L ，对  i ,  2 ,,  n 是
否相关不下结论。
d L ， dU 的大小可由附表7查得。
（2） 模型修正
经检验发现模型的随机误差序列是相关时，常可采用下面的方法对
模型进行修正。
方法一：加权最小二乘估计
若我们可以得到 D(Y )   G 中 G 的表达式（此时， G 是一般的
2
正定矩阵，而非对角矩阵）
，那么，可采用加权最小二乘估计去求模型
的参数估计。
方法二：可用差分法
令： y  y
i
i 1
y ,
i
x  x
x ,
i
i 1
i
i  1,2,, n  1
建立 y 关于 x 的回归方程。由此可得：
y   0  1x
对此模型再进行DW检验，若认为 y i 之间是不相关，则可采用下列方程为最
终形式：
yˆ
i 1
 y   0  1 ( x
x )
i
i 1
i
方法三：可用迭代法
其步骤如下：

先建立 y 关于 x 的回归方程，并求出相应各点的残差 e1 , e2 ,
, en ，
经DW检验，若认为存在序列相关，则进行下述各步骤;

利用已求出的残差，求序列 e1 , e2 ,, en1 与序列 e2 , e3 ,, en 的相关
系数，记为 r

(1)
令：  y  y
(1)
i
；
 r (1) y , (1) x  x
 r (1) x , i  1,2,, n  1 ；
i 1
i
i
i 1
i

(1)
先建立  y 关于  x 的回归方程
(1)
(1) y   0  1(1) x ，
(1)
(1)
(1)
并求出相应各点的残差 e1 , e2 ,, en1 ，经DW检验若认为已不存在序列相关，
则把上述方程认为是最终所求方程；若经检验仍存在序列相关，则求出序列
e1(1) , e2(1) ,, en(1)2 与序列 e2(1) , e3(1) ,, en(1)1 的相关系数，记为 r ( 2 ) ；

令

现重新建立回归方程，并进行DW检验，一直对DW检验表明序列不相关为
(2)
止。
(2) (1)
(2)
(1)
(2) (1)
y  yi(1)

r
y
,

x

x

r
xi , i  1, 2,
1
i
i 1
i
i
,n  2 ；
四、误差的正态性诊断
（一）模型诊断
检验总体分布是否为正态的方法有很多。常用的有如下几种：
方法一：在第本章的第二节已讨论过，若 Y ~ N ( X ,  I n ) 时，标准化
2
残差 r1 , r2 ,, rn 可近似看成来自 N (0,1) 的随机子样，从而可粗略地统计
一下 r1 , r2 ,, rn 中正负个数是否则各占一半左右，介于(-1,1)之间的比
例是否约为68%,介于（-2，2）之间的比例是否约为95%,介于(-3,3)之间
的比例是否约为99%.
方法二：正态检验方法。如数理统计中介绍过的柯莫
哥洛夫检验，小样本的夏皮诺和威尔克的W检验，大样
本的D检验等。
(二) 模型修正
当模型的随机误差序列为非正态时，我们希望能找到一个变
换，使变换后有：
 y1   g ( y1 ) 
  

 y 2   g ( y 2 ) 
2

~
N
(
X

,

In )
n
     
  

 y   g ( y ) 
n 
 n 
Box与Cox提出了一类应用较广的Box-Cox变换：
y
( )
 y 1

 
ln y

0
 0
要求变换后
Y ( )
 y1(  ) 
 ( ) 
 y2 
2

~
N
(
X

,

In )
n

  
 y ( ) 
 n 
从上式可以看出，Box-Cox变换不仅在误差非正态场合可用，
而且在回归函数非线性，误差方差非齐性，观测值间不独立
时均可选用，因为通过这个变换，线性回归中四条假定都满
足。
下面分两步来找出参数  ,  ,  的极大似然估计。
2
Y (  ) 的似然函数为：
LL(  ,  2 ,  )  (2 2 )

n
2
 ( )
 1

( )
exp
Y

X

Y

X


2
 2




利用变换的形式可知 Y 的似然函数为：
L(  ,  ,  )  (2 )
2
其中： J 
2
n

i 1

n
2
 ( )
 1
( )
exp
Y  X Y  X
2
 2
n
dyi(  )
  yi 1
dyi
i 1




 J

对固定的  ，可得  ,  2 的极大似然估计(MLE)为：
1
( )
ˆ


  ( X X ) X Y
 1

1
ˆ 2  Y (  ) [1  X ( X X ) 1 X ]Y (  )  SSE ( , y )
n
n
对固定的  ，有：

n
2
Lmax ( )  (2ˆ 2 ) e

n
2
J
若略去常数不计，仍记其值为 Lmax ( ) ，并取对数有：

( )
( ) 
n
n Y
Y

ln Lmax ( )   ln[SSE( , y )]  ln J   ln  1 n [1  X ( X X ) 1 X ] 1 n 
2
2 J
J 


若令： Z
( )

zi
Y ( )
J
( )
1n

，则有：
yi
J
( )
1n
 1
n
 ( )
( yi ) n ,
 yi

i 1

 1
n
(ln y )( y ) n ,
i  i

i 1
0
i  1, 2,
 0
,n
并记： SSE ( , Z )  Z
则：
( )

[1  X ( X X ) 1 X ]Z (  )
n
Lmax ( )   ln[ SSE ( , Z )]
2
为找出  使 Lmax ( ) 达到最大，只要使 SSE( , Z ) 达到最小即可。其
解析解是比较难找的，通常的做法是给出一系列的  值，画出  与
SSE( , Z ) 的曲线，从图上找出使 SSE( , Z ) 达到最小的近似 ˆ 。并
用此 ˆ 作Box-Cox变换并求出回归方程。
回归诊断在SAS上的实现
• 用语句plot r.*p. (r是residual的缩写，p是
predicted的缩写)可以作残差r相对于拟合值
p之间的散点图。如果此散点图在0水平线
上下均匀散布，且对p没有趋向性，则可认
为误差满足方差齐性假设、独立性假设，
和回归函数线性假设合理。
• model y=x/dw r;
选项里加上dw表示计算DW检验的值。
r表示计算学生化残差，并计算Cook距离，
若Cook距离相对较大,则认为是强影响点。
若学生化残差的绝对值大于2，则可认为是
异常点。从学生化残差也可判断误差的正
态性假设是否满足。
• 例子.
给10只大白鼠注射内霉素(30mg/kg)后，测
得每只大鼠红细胞x与血红蛋白含量Y数据
(见下页SAS文件)，试对X和Y进行回归分
析。
data mouse;
input x y;
cards;
654 130
786 168
667 143
605 130
761 158
642 129
652 151
706 153
602 151
539 109
;
proc reg;
model y=x;
run;
proc reg;
model y=x/noint dw r cli clm;
plot r.*p.;
run;
注意：开始的时候并没有noint，
分析后才加入noint
残差图
y = 0.2148 x
25
N
10
Rsq
0.9958
AdjRsq
0.9953
RMSE
9.7628
20
15
10
5
0
-5
-10
-15
115
120
125
130
135
140
145
Predicted Value
150
155
160
165
170
误差的独立性诊断
第九个为异常点、强影响点
Output Statistics
Std Error Student
Obs Residual Residual -2-1 0 1 2
1
2
3
4
5
6
7
8
9
10
9.279
9.056
9.259
9.350
9.101
9.297
9.282
9.196
9.354
9.437
-1.129
-0.0913
-0.0288
0.00543
-0.600
-0.957
1.180
0.148
2.319
-0.718
|
|
|
|
|
|
|
|
|
|
**|
|
|
|
|
|
|
|
*|
|
*|
|
|** |
|
|
|**** |
*|
|
Cook's
D
0.136
0.001
0.000
0.000
0.054
0.094
0.148
0.003
0.480
0.036
Sum of Residuals
1.36513
Sum of Squared Residuals
857.80435
Predicted Residual SS (PRESS) 1031.32594