回归分析中异常值的诊断与处理

回归分析中异常值的诊断与处理1 引言对调查得到的大量原始数据进行加工处理，提取其中有用的信息，即统计整理，是进一步进行统计分析的前提．但是通常的统计整理方法往往都对样本数据有一个前提假设，即样本数据是来自同一个总体，而这个假设有时却不能成立．原因一是由客观因素造成的，如总体条件的突然变化或人们未知的某个因素的突然出现等等；二是由主观方面的因素造成的，即人为的差错如调查人员读错或抄错数据，不小心把另一些不同条件下的样本数据混杂进来．当样本中的个别数据明显的偏离样本中其余数据时，这些数据可能是来自不同的总体，我们称这样的数据为异常数据．若对混有异常数据的样本按常规进行统计整理、分析、推断，往往会得出不符合实际的结论．本文就样本中的异常数据提出了诊断方法和处理方法，并结合实例说明了简单回归中上述方法的具体应用．2 异常值的概念所谓异常值就是在所获统计数据中相对误差较大的观察数据，也称奇异值(Outliers)．狭义地定义异常值就是一批数据中有部分数据与其余数据相比明显不一致的数据，也称离群值．社会经济统计中一切失实数据统称为异常值．由于人为或随机因素的影响，失实的数据随时都有可能出现，因而统计数据中的任何一个都有可能成为异常值，而狭义界定的异常值是指离群值，如果把统计数据按由小到大排列，若有异常值，它必位于其数据的两端，左端称为异常小值，右端的称为异常大值．残差：考虑线性回归模型y X e β=+，()0E e =，()2n Cov e I σ=．()1 其中 ⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=n y y y y 21，⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡=---1,11,2211,111111p n n p p x x x x x x X ，⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡=-110p ββββ ，⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=n e e e e 21． 分量形式为011,11i i i p p i y x x e βββ--=++++，n i ，， 1=．()2 定义ˆˆˆey y y X β=-=-为残差向量，其中ˆˆy X β=称为拟合值向量，βˆ为β在模型()1下的最小二乘估计．如果用n x x '',1表示X 的n 个行向量，则称 ˆˆi i i ey x β'=-， n i ，， 1=， 为第i 次试验或观测的残差．对简单回归，12111n x x X x ⎡⎤⎢⎥⎢⎥=⎢⎥⎢⎥⎣⎦，⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=n y y y y 21，⎥⎦⎤⎢⎣⎡=10βββ．则ˆˆiiie y x β'=-，n i ，， 1=，其中()1,iix x '=．于是，异常值就是在回归分析中，一组数据()i i y x ,'如果它的残差i eˆ较其它组数据的残差大的多，则称此数据为异常值．异常值的出现有主客观的原因．主观上抽样调查技术有问题，疏忽大意记错，或人为的虚报，谎报数据等都可能导致异常数据，这直接影响了统计数据的质量；客观上某些样品由于特定原因在某些变量上的确表现突出，明显超出平均水平，这也可能产生异常值．异常值的存在必将导致相应统计分析误差增大，会对分析结果(平均值与标准差)产生重要影响，会降低测量的精度．如果不预先处理它们，用通常的统计整理方法所得出的结论可靠性差．而异常值的诊断与处理是保证原始数据可靠性，平均值与标准差计算准确性的前提．3 异常值的诊断方法在讨论异常值诊断问题时，通常要假设所得样本观测值在某中意义下遵从一定的分布规律．拿到一批数据，若能从其实际背景中明确看出它服从某中分布形式时，一般的做法是在这种分布假设下，导出能较好反映异常值与正常值差异的统计量，在没有异常值的原假设下作假设检验．以下给出两种检验方法．3.1 F 分布检验法学生化残差：考虑线性回归模型()1，记ˆˆyX β=，称y ˆ为拟合值向量，称其第i 个分量βˆˆi i x y'=为第i 个拟合值，则 ()1ˆy X X X X y Hy -''==，这里()1H X X X X -''=．文献中通常称H 为帽子矩阵．前面已经定义了ˆˆi i i ey x β'=-, n i ，， 1=， 为第i 次试验或观测的残差．将其标准化为iii h e-1ˆσ，再用σˆ代替σ，得到所谓学生化残差ii i i h e r -=1ˆˆσ， n i ， ,1=，这里ii h 为H 的第i 个对角元，pn eni i-=∑=122ˆˆσ．把正态线性回归模型()1改写成分量形式i i i e x y +'=β，()2,0~σN e i ，n i ，， 1=，这里i e ，)1(n i ，， =相互独立．如果第j 组数据()j j y x ,'是一个异常点，那么它的残差就很大．它的残差之所以很大是因为它的均值()j y E 发生了非随机性漂移η．从而()ηβ+'=j j x y E ．这产生了一个新的模型i i ijj j y x e y x e ββη'=+⎧⎪⎨'=++⎪⎩ ()2~0,,i i j e N σ≠ ()3 记[]0,,0,1,0,,0j d '=．将模型()3改写成矩阵的形式j y X d e βη=++， ()I N e 2,0~σ， ()4模型()3和()4称为均值漂移线性回归模型．要判定()j j y x ,'不是异常点，等价于检验假设0=η：H ．引理1 用()i y ，()i X 和()i e 分别表示从Y 、X 和e 剔除第i 行所得到的向量或矩．从线性回归模型()1剔除第i 组数据后，剩余的1-n 组数据的线性回归模型为()()()i i i y X e β=+，()()0i E e =，()()21n i Cov e I σ-=． ()5 将从这个模型求到的β的最小二乘估计记为()i βˆ，则 ()()i iii i x X X h e11ˆˆˆ-'--=ββ． 证明：因为()()()()()()1ˆi i i i i X X X y β-''=． ()6设A 为n n ⨯可逆阵，v u 和均为1⨯n 向量．用恒等式()vA u A v u A A v u A 111111-----'-'+='-有()()()()()()()111111i i i i iiiiX X x x X X X X X X x x X X h -----'''''''=-=+-， ()7这里i x '为X 的第i 行．将上式两边右乘X y '，并利用()()i i i i X y X y y x ''=+以及()6式，有()()()()()()iii i ii i i i h x x X X x X X y -''-'+=--1ˆˆˆ11βββ．()8将()7式右乘i x ，可以得到如下关系式()()()()i iii iix X X h x X X 1111--'-='． 将其代入()8式，得到()()i iii i x X X h e11ˆˆˆ-'--=ββ． 引理2 对均值漂移线性回归模型()4，β和η的最小二乘估计分别为()j ββˆ=*和j jje h ˆ11-=*η，其中()j βˆ为从非均值漂移线性回归模型()1剔除第j 组数据后得到的β的最小二乘估计．()()X X X X h H jj ''==-1，jj h 为H 的第j 个对角元．j eˆ为从模型()2导出的第j 个残差． 证明：显然，j j y y d ='，1='j j d d ．记()12,,,n X x x x '=．则j j x d X ='．于是，根据定义()111j j j j j j X X X Xx X y Xd y d d x y βη-*-*''''⎡⎤⎛⎫⎡⎤⎡⎤⎛⎫⎡⎤==⎢⎥ ⎪ ⎪⎢⎥⎢⎥⎢⎥'''⎢⎥⎣⎦⎝⎭⎣⎦⎣⎦⎝⎭⎣⎦． 根据分块矩阵的逆矩阵公式(见附录1)，以及()j j jj x X X x h 1-''=．()()()()()()()()111111111111111111ˆˆ1111ˆ111ˆˆ11ˆ1j j j jj jj j j jj jjj j j j jj jjj j jj jj j j jjj jj X X X X x x X X X X x h h X y y x X X h h X X x x X X x y h h x y h h X X x e h e h βηββββ----**----⎡⎤'''''+-⎢⎥--'⎡⎤⎡⎤⎢⎥=⎢⎥⎢⎥⎢⎥⎣⎦⎣⎦''-⎢⎥--⎣⎦⎛⎫'''+- ⎪--⎪= ⎪'-+ ⎪ ⎪--⎝⎭⎛'--=-⎝⎫ ⎪⎪ ⎪ ⎪ ⎪⎭． 再由引理1知命题得证．现在应用引理来求检验0=η：H 的检验统计量．注意到，对现在的情形，在约束条件0=η下，模型()4就化为模型()1，于是=H RSS 模型()1无约束情形下的残差平方和 y X y y ''-'=βˆ． 而模型()4的无约束残差平方和y d y X y y RSS j '-'-'=**'ηβ． ()9利用引理2得()2ˆˆˆ1ˆˆ111H jj j j j j jj jj jjRSS RSS X y d y e y e e x h h h ββηβ**'''-=-+'=-+=---．这里βˆˆj j j x y e'-=为第j 组数据的残差． 利用*β和*η的具体表达式将()9式作进一步化简：jj j j jjj j h y e h y e y X y y RSS ---+''-'=1ˆ1ˆˆˆβ=()jjj h e p n ---1ˆˆ22σ，其中pn eni i-=∑=122ˆˆσ．根据引理2，所求的检验统计量为()()()()22222111ˆ1ˆ1ˆ1jj jj j jj jH r p n r p n h p n e p n p n h e p n RSS RSSRSS F ----=--------=---=σ． 于是，我们证明了如下事实：定理[]2 对于均值漂移线性回归模型()4，如果假设0=η：H 成立，则()1,122~1------=p n j j j F r p n r p n F ．据此，我们得到如下检验：对给定的()10<<αα，若()()α1,1221-->----=p n jj j F rp n r p n F ，则判定第j 组数据()j j y x ,'为异常点．当然，这种检验会犯“判无为有”的错误，也就是()jjy x ,'可能不是异常点，而被误判为异常点．但我们犯这种错误的概率只有α，事先我们可以把它控制的很小．显然，根据t 分布与F 分布的关系，我们也可以用t 检验法完成上面的检验．若定义()()212211⎥⎥⎦⎤⎢⎢⎣⎡----==j j j j r p n p n r F t ． 对给定的α，当⎪⎭⎫⎝⎛>--21αp n j t t时，我们拒绝假设0=η：H ．即判定第j 组数据()j j y x ,'为异常点．3.2 残差及残差图检验异常值前面定义了βˆˆX y e-=，称为残差向量，其分量形式ˆˆi i i e y x β'=-，()1,,i n =，称为第i 次试验或观测的残差．特别地，对简单回归，()i i x x ，1='，n i ，， 1=．⎪⎪⎪⎪⎪⎪⎭⎫ ⎝⎛---=⎥⎦⎤⎢⎣⎡=∑∑==n i i n i i i x n x y x n y x x y 1221110ˆˆˆˆββββ．所以，()i i i x y e 10ˆˆˆββ+-=，n i ，， 1=． 残差是最重要的一种回归诊断量，它蕴涵了有关模型基本假设的许多重要信息．残差分析就是对残差进行统计处理，从中提炼出这些信息的方法．而残差图就是残差分析中使用的基本工具．所谓残差图就是残差i eˆ对因变量y 或自变量n X X X ，， 21,，或其它导出统计量（如拟合值i yˆ）的点子图，有时候也用残差对时间或对数据序数的点子图．最简单的图，尤其在简单回归中，为残差i eˆ对拟合值i y ˆ的图． 所谓异常数据就是相对于其它观测值来说，具有大的残差的数据点．利用残差及残差图检验异常值的方法是用所给数据计算出残差i eˆ，与其余观测值的残差进行比较，具有大的残差的数据点被怀疑为异常值．然后作出自变量与因变量的散点图，残差i eˆ对拟合值i y ˆ的残差图以及残差i eˆ对自变量i x 的残差图，从图中观察，那些远离大多数观测点的孤立的点有理由被认为是异常点．然后从数据中删除这些点，再次估计回归方程，作出X 与Y 的散点图以及i eˆ对i y ˆ的残差图，计算标准差，与删除前进行比较． 4 提出两种处理方法4.1 采用虚拟变量消除异常值的影响（1） 虚拟变量某类变量（如性别、种族、颜色、宗教、国家、战争、地震等）常常表示某属性是否存在，如男或女，黑种人或白种人，教徒或非教徒，对于这类变量可用“数量化”方法设计人为变量来表示，如用 “0”或“1”，“0”表示不存在某种属性，“1”表示存在该属性．例如：用“1”表示男性，“0”表示女性，或用“1”表示大学毕业，“0”表示没有大学毕业，等等．这样假设的变量，如“0”或“1”称为虚拟变量，也称为类型变量或属性变量，它的主要特点就是将一些可以划分为不同类型或属性的变量用“0”和“1”分别表示，即属于某一类型的变量用“1”表示，不属于这一类型的变量用“0”表示．这里用D 表示虚拟变量．（2） 处理异常值的方法在这里，我们用虚拟变量（0或1）把试验数据划分为两类，属于异常值一类用“1”表示，属于正常值一类用“0”表示．引入虚拟变量D ，建立回归方程：01122221ˆˆˆˆˆp p p Y X X X D βββββ---=+++++， 其中⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=n y y y Y 21，111,11212,121,1111p p n n p n x x d x x d X x x d ---⎡⎤⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎣⎦，011ˆˆˆˆp ββββ-⎡⎤⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎣⎦，12n d d D d ⎡⎤⎢⎥⎢⎥=⎢⎥⎢⎥⎣⎦．特别地，对简单回归，建立回归方程：012ˆˆˆY X D βββ=++， 其中，⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=n y y y Y 21，1122111n n x d x d Z x d ⎡⎤⎢⎥⎢⎥=⎢⎥⎢⎥⎣⎦，()0112ˆˆˆˆZ Z Z y ββββ-⎡⎤⎢⎥''==⎢⎥⎢⎥⎢⎥⎣⎦． 利用所给数据，作出回归估计．这样，异常值的影响被虚拟变量D 的系数吸收，从而估计更接近于现实．4.2 剔除异常值的方法对于测定中的异常值的剔除，我们必须持慎重态度，不能贸然从事．否则会出现误删有效数据或保留异常数据的错误．通常处理步骤大致如下：（1）初分析：首先利用所给数据作出回归估计：βˆˆX y=． 其中111,1212,11,1111p p n n p x x xx x x x ---⎡⎤⎢⎥⎢⎥=⎢⎥⎢⎥⎢⎥⎣⎦，()y X X X p ''=⎥⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎢⎣⎡=--1110ˆˆˆˆββββ ．特别地，对简单回归，x y 10ˆˆˆββ+=，其中x y 10ˆˆββ-=，∑∑==--=ni ini ii x n xy x n yx 12211ˆβ，并计算标准差：s =值i yˆ、残差i e ˆ，并列入表．然后作出i x 与i y 的散点图，拟合值i y ˆ与残差i e ˆ的残差图，分析判别出异常值．（2）判别出异常值后，从测量数据中删除异常值，然后利用余下的观测值再次作回归估计，同第一步，计算标准差、拟合值以及残差，作出删除异常数据后的散点图与残差图．（3）对比删除异常值前后的标准差以及残差图．对于单一可疑异常数据，直接按上述步骤进行剔除异常值．对于多个可疑异常值的情况，上述剔除步骤反复逐次进行至无可剔除为止．通常可采用两种方法：向前逐一剔除法（从极大值依次向次大值逐一剔除）和向后逐一剔除法（从可疑的j 个最大值依次从小到大逐一剔除）．5 实例分析5.1 实例一假定某调查数据X 和Y 由表1中给出．表1 数据X 和Y根据表中数据作回归估计，得回归方程：ˆ0.6610.837yx =+．并算得 ()()457.0ˆ112212=---=∑∑==n i ini i iy yy yR ．表2 各组数据的诊断统计量()120.05 3.65t =．从表2最后一列可以看出只有841.96=t 超过这个值，于是我们诊断出第六号数据为异常值．然后引入虚拟变量D ，建立回归方程：D x Y 210ˆˆˆβββ++=，作回归估计．令 βˆZ Y =，⎥⎥⎥⎥⎦⎤⎢⎢⎢⎢⎣⎡=15152211111d x d xd x Z ，()y Z Z Z ''=⎥⎥⎥⎦⎤⎢⎢⎢⎣⎡=-1210ˆˆˆˆββββ． 得0ˆ0.303β=-，1ˆ 1.007β=，2ˆ 4.985β=．所以，回归模型为：ˆ0.303 1.007 4.985y x D =-++．并且算得993.02=R ．从2R 的值来看，这个回归方程［见图2中()2］要比未采用虚拟变量所估计的回归方程［见图2中()1］其拟合程度好多了．这是因为异常值的影响被虚拟变量D 的系数吸收，所以比处理前求得的估计值更接近于现实．因此虚拟变量是消除异常值的影响，探求变量之间真正关系的一种有效方法．01234567812345675.2 实例二表3给出了Forbes 数据．表3 Forbes 数据根据表中数据，(气压用100log ⨯的数据)，估计回归方程，算得x y 895.0131.42ˆ+-=．计算出Forbes 数据的拟合值及残差，如表4，进一步作出散点图（图3）、残差i eˆ对拟合值i y ˆ的残差图（图4）．表4 Forbes 数据的拟合值及残差-0.6-0.4-0.200.20.40.60.811.21.41.6比较表中所列的残差，并观察残差图，我们发现序号12的残差比其它的残差大的多．其它残差的绝对值都小于0.35，而第12号数据的残差为1.36．并且从残差图（图4）中观察到第12号数据的跳跃度比较大，远离其它点．由残差及残差图法诊断出第12号数据为异常数据．然后剔除第12号数据．再估计回归方程，重新拟合数据，检查参数估计，拟合值，残差方差的变化．删除第12号数据后，得到回归方程：x y891.0302.41ˆ+-=． 表5 Forbes 数据中删除12号数据后的拟合值及残差图４ Forbes 数据的残差图ˆi yˆi e-2.5-2-1.5-1-0.50.511.522.5表6 删除异常值前后的比较量对删除异常值前后的结果进行对比，我们得到对于删除前后获得的参数估计基本是相同的，所以第12号数据是无关的．但是，第12号数据在标准误差上的影响要显著些．删除后标准误差减小约3.1倍，方差减小约10倍．并且，从删除异常值获得的残差图()图5和散点图（图6）显示出，对余下的16个数据无明显的拟合失败．图5 Forbes 数据中删除12号数据后的残差图ˆi yˆi e6 结束语异常值的诊断与处理问题是一个很有实际意义的问题，是一个逐步研究逐步实践的问题．但迄今为止，还没有一种广泛适用的方法，本文所介绍的方法也是一种尝试，对于异常值的诊断与处理问题的研究还有待于进一步完善．致谢本文在撰写过程中得到了李惠东老师的悉心指导和大力支持．在论文撰写过程中遇到的难点和疑点李老师都给我作了详细的讲解，并提供了许多有参考价值的资料和专业软件，使我学到了很多知识和方法．在此表示深深的敬意与感谢．参考文献[1] 魏立力.概率论与数理统计[M].银川:宁夏人民出版社,1999.[2] 王松桂等.线性模型引论[M].北京:科学出版社,2004.[3] 何平.剔除测量数据中异常值的若干方法[J].数理统计与管理,1995,(1):19-22.[4] 王静龙,梁小筠,等译.应用线性回归[M].北京:中国统计出版社,1998.[5] 刘宗鹤,赵明强译.计量经济学概论[M].北京:农业出版社,1988.附录1：分块矩阵的逆矩阵：设0>A ，将其分块为：⎥⎦⎤⎢⎣⎡=22211211A A A AA ，则它的逆矩阵 ⎥⎦⎤⎢⎣⎡--+=-⋅--⋅-⋅---⋅---1122111211122112212111111211122*********A A A A A A A A A A A A A A=⎥⎦⎤⎢⎣⎡+----⋅---⋅---⋅-⋅122121************211211221221212111211A A A A A A A A A A A A A ． 这里 121112122122A A A A A -⋅-=，211221211211A A A A A -⋅-=．附录2：SPSS处理实例1数据：1．未用虚拟变量的数据处理：RegressionModel Summary(b)b Dependent Variable: YCoefficients(a)Casewise Diagnostics(a)2．引进虚拟变量后的数据处理：Model Summary(b)b Dependent Variable: YCoefficients(a)Casewise Diagnostics(a)附录3：SPSS处理实例2的数据：1．删除前的数据处理：Coefficients(a)a Dependent Variable: YModel Summary(b)a Predictors: (Constant), Xb Dependent Variable: YCasewise Diagnostics(a)a Dependent Variable: Y2．删除后的数据处理：a Dependent Variable: YModel Summary(b)b Dependent Variable: YCasewise Diagnostics(b)a Missing Caseb Dependent Variable: Y。