王純杰，　鄭　茜，　李　群，　林珊屹

(長春工業(yè)大學 基礎科學學院， 吉林 長春　130012)

?

基于SAS的回歸分析異方差的檢驗與消除

王純杰，鄭茜*，李群，林珊屹

(長春工業(yè)大學 基礎科學學院， 吉林 長春130012)

摘要：利用SAS軟件的REG過程步實現(xiàn)了回歸分析中異方差的White檢驗，利用加權最小二乘法消除了異方差的影響。通過與殘差圖法和等級相關系數(shù)法的結果比較，該檢驗效果更精確。利用2013年全國31個省市自治區(qū)保險賠付和保費數(shù)據(jù)進行實證分析。

關鍵詞：回歸分析； 異方差檢驗； 加權最小二乘估計

0引言

根據(jù)回歸分析中最小二乘估計的結果，可以得到自變量與因變量之間的關系，但是模型中可能存在整體模型或者個別變量不顯著的情況，由多種原因導致，如異常值與強影響點，多重共線性等。異方差是常見原因之一，消除異方差可以使模型顯著。異方差問題分為兩大類，即異方差的診斷和消除。關于異方差診斷，常見的方法有殘差圖檢驗法和斯皮爾曼等級相關系數(shù)檢驗法。這兩種方法可以判定是否存在異方差，但是它們有一定的局限性，可見探索其他的方法很有意義。

1回歸方程和普通最小二乘估計[1]

用y表示因變量，X=(X1,X2,…,Xp)T是對y有影響的p個自變量，回歸模型一般表述為:

(1)

(2)

即

(3)

回歸模型的基本假設如下：

1)E(ε)=0，var(ε)=σ2I;

2)E(X′ε)=0;

3)rank(X)=p+1。

于是有E(y)=Xβ,cov(εi,εj)=σ2，線性回歸模型簡記為(y,Xβ,σ2In)[1]。

定義離差平方和

(4)

回歸系數(shù)β的估計為

回歸方程的估計式為:

(5)

根據(jù)模型假設，回歸模型要檢驗3個方面：

1)回歸模型的顯著性；

2)每個自變量的顯著性；

3)殘差是否為“純隨機”。

2 異方差的檢驗

2.1異方差產(chǎn)生的原因和影響

經(jīng)典回歸假設中一個重要的假設是var(εi)=σ2，即擾動項的同方差性。當模型中隨機擾動項的方差不是常數(shù)時，稱為異方差或方差非齊性。記為：

(6)

隨機誤差項包含眾多因素對解釋變量的影響，如果它們隨著解釋變量觀測值的變化而變化，則會對解釋變量產(chǎn)生影響。異方差現(xiàn)象存在時，參數(shù)估計值雖然是無偏的，但不是有效的，即不再具有最小方差性。以一個簡單線性回歸模型為例，模型如下:

(7)

對回歸系數(shù)β1進行普通最小二乘(OLS)估計，可得:

(8)

(9)

(10)

(11)

則

(12)

則

2.2異方差的檢驗

關于異方差性的檢驗[2]，統(tǒng)計前輩們進行了大量的研究，提出的診斷方法很多，但是沒有一個公認的最好的方法。下面先介紹一下殘差圖分析法和等級相關系數(shù)法，這是兩種比較常見的方法[2]。

2.2.1殘差圖分析法

殘差圖分析法非常簡單、直觀，它是一種非正式的檢驗方法。如果在檢驗前不能判斷是否存在異方差，可以采用該法。做法是先不考慮存在異方差的假定下構造回歸模型,然后以殘差ei為縱坐標，以其他合適的變量為橫坐標，畫散點圖。常用的橫坐標有3種：

1)以解釋變量xi(i=1,2,…,p)為橫坐標；

3)以觀測值或序號為橫坐標。

如果回歸模型適合樣本數(shù)據(jù)，那么殘差ei應該反映εi所假定的性質，即殘差圖上的n個點應該是隨機分布和沒有任何規(guī)律的。如果回歸模型存在異方差，殘差圖上點的分布會呈現(xiàn)一定的趨勢，如殘差ei的值隨y值的增大而增大(或減小)，具有明顯的規(guī)律，由此可以判斷模型的隨機誤差項εi的方差是非齊性的，存在異方差。

2.2.2等級相關系數(shù)法

等級相關系數(shù)法又稱斯皮爾曼(Spearman)檢驗，是一種應用比較廣泛的方法。進行等級相關系數(shù)檢驗有3個步驟：

1)求出y關于x的普通最小二乘估計，求出εi的估計值，即ei的值。

2)取ei的絕對值，即|ei|，把xi和|ei|按遞增(或遞減)的次序排列，分成等級，按下式計算等級相關系數(shù):

(13)

式中:n----樣本容量;

di----對應xi和|ei|的等級的差數(shù)。

3)做等級相關系數(shù)的顯著性檢驗。在n>8的情況下，用下式對樣本等級相關系數(shù)rs做t檢驗，檢驗統(tǒng)計量為:

(14)

2.2.3White檢驗

(15)

(16)式中：n----觀測值的個數(shù);

p----回歸參數(shù)的個數(shù)，包括常數(shù)項。

(17)

(18)

研究發(fā)現(xiàn)，當樣本量小于250時，HC3的性質最優(yōu)。

在SAS中,可以通過在model語句后面加入HCCMETHOD=0(,1,2,3)選項來選擇異方差性相容協(xié)方差矩陣的估計方法，默認選項是HC0。

model語句后面的ACOV選項的作用是在參數(shù)估計的結果中顯示異方差性相容協(xié)方差矩陣，而且增加了異方差性標準差，就是White標準差。如果在model語句后面選擇了HCC選項或者WHITE選項，但是沒有選擇ACOV選項，結果中會顯示異方差性標準差，但是不會顯示異方差性相容協(xié)方差矩陣。SPEC選項是對模型做一個檢驗，原假設是回歸估計的誤差項,是獨立同方差的。備擇假設是誤差項具有異方差。可以指定model語句后面的SPEC，ACOV，HCC，White選項，來選擇檢驗的類型，檢驗給出的一致協(xié)方差矩陣都是漸近的。ACOV選項和SPEC選項可以用在model語句或者print語句中。

3異方差的消除

如果模型存在異方差，就違背了回歸模型的基本假定，不能使用普通最小二乘法估計參數(shù)了。如果對原有的模型進行變換，使其滿足同方差性的假設，再進行參數(shù)估計，可以得到理想的模型。最常用的方法是加權最小二乘法[4]。

3.1加權最小二乘估計

假設觀測數(shù)據(jù)y與參數(shù)β1,β2，…，βp之間服從如下線性關系:

(19)

其中,ε是擾動項，進行n次觀測可以得到n個線性方程:

(20)

式(20)可以用矩陣表達，即:

(21)

(22)

(23)

令式(23)為0，解得:

(24)

式(24)就是最小二乘估計的表達式，它與觀測數(shù)據(jù)y和系數(shù)陣X有關[5]。

(25)

(26)

解得:

(27)

3.2加權最小二乘估計權數(shù)的確定

加權最小二乘估計誤差的方差推導過程如下：

(28)

我們的目標是求出使式(28)達到最小的w, 即加權系數(shù)，由此引入矩陣型許瓦茲不等式。

定理1設A,B分別為m×n和n×l矩陣，且AAT是可逆矩陣，則有：

(29)

當且僅當存在一個m×l維矩陣P，使得B=ATP時式(29)取等。

證明：

BTB-2BTAT(AAT)-1AB+BTAT(AAT)-1AAT(AAT)-1AB=

(30)

所以BTB≥(AB)T(AAT)-1AB得證。

當且僅當B=AT(AAT)-1AB時,式(30)右端為零矩陣。令P=(AAT)-1(AB)，有B=ATP；反之,如果B=ATP，則有AT(AAT)-1AB=AT(AAT)-1AATP=ATP=B，由此證明了式(28)中等號成立的充要條件是B=ATP，證畢。

下面利用這個定理求加權系數(shù)。因為var(ε)=σε是對稱的正定矩陣,令σε=PTP,其中P為n×n的可逆矩陣。取l=n,令A=XTP-1,B=PwX(XTwX)-1,可以看出AB=1,所以對任意的加權矩陣可以得到如下不等式：

BTB≥

(31)

比較式(31)的兩邊，有

(32)

此時加權最小二乘估計的誤差方差達到最小,即：

(33)

此時的估計為：

(34)

3.3加權最小二乘估計的步驟

1)將回歸方程的各項除以比例因子Zi，得到新的回歸方程，同時滿足誤差項的同方差性；

2)對新回歸方程進行最小二乘估計，求得回歸參數(shù)。

4實證分析

SAS系統(tǒng)(Statistics Analysis System)的重要組成部分和核心功能是統(tǒng)計分析，最新版本是9.4版。 SAS可以提供多個統(tǒng)計分析過程，每個過程均含有豐富的選項，滿足客戶的不同需求[6]。

若想對某一實際現(xiàn)象作回歸分析，首先要檢驗數(shù)據(jù)是否存在異方差現(xiàn)象，然后選擇合適的方法進行消除，最后觀察模型擬合數(shù)據(jù)的程度。2013年全國31個省市自治區(qū)保險保費收入和賠付支出情況見表1。

表1　2013年全國31個省市自治區(qū)保險保費收入和賠付支出情況

文中利用的數(shù)據(jù)是2013年全國31個省市自治區(qū)保險賠付支出和保費收入的數(shù)據(jù)，該數(shù)據(jù)是橫截面數(shù)據(jù)(來自中國統(tǒng)計年鑒，2014)。設保費收入為自變量x，賠付支出為因變量y，單位是億元，由此建立回歸模型，分析自變量對因變量產(chǎn)生的影響。

在進行回歸分析之前，需要先進行相關分析，利用SAS軟件的CORR過程步可以實現(xiàn)。因變量支出與自變量收入的相關系數(shù)是0.989 75，在顯著性水平為0.05下，通過檢驗，說明自變量與因變量之間有很大的相關性, 可以進行回歸分析[7]。 SAS程序如下：

data expend2013;

input expend income;

cards;

318.17994.44

102.00276.80

…

24.0472.70

106.59273.49

;;

proc corr data=expend2013;

run;

proc reg data=expend2013;

model expend=income;

output out=expend2013_r r=residual;

run;

quit;

普通回歸分析結果見表2。

表2　方差分析表

因為數(shù)據(jù)的單位小，每個變量的值很大，所以導致總平方和、離差平方和、回歸平方和都很大，F(xiàn)值為1 393.24，P值<0.000 1，表明回歸方程高度顯著，說明自變量對因變量有高度顯著影響。通過復決定系數(shù)R2=0.979 6，也可以判定回歸方程高度顯著。在顯著性水平為0.05以下，自變量的t檢驗通過了，說明它對因變量有顯著影響，常數(shù)項也通過了檢驗，下面開始檢驗異方差[8]，SAS程序如下：

proc gplot data=expend2013_r;

plot residual*income=1;

symbol1 c=blue v=circle i=none;

label income="保費" residual="殘差";

run;

殘差關于自變量收入的散點圖如圖1所示。

圖1殘差圖

可以看出殘差不是隨機分布的，有向外擴張的趨勢，初步判定存在異方差現(xiàn)象。下面進行等級相關系數(shù)檢驗，程序如下：

data r2013_abs;

set expend2013_r;

abs_r=abs(residual);

run;

proc corr data=r2013_abs spearman;

var abs_r income;

label income="保費" abs_r="殘差絕對值";

run;

等級相關系數(shù)檢驗結果見表3。

表3通過SAS軟件得到自變量收入與殘差絕對值|ei|的相關系數(shù)是0.572 18，P值為0.000 8，自變量收入與殘差絕對值|ei|顯著相關，即存在異方差現(xiàn)象，與殘差圖的結果吻合。接下來進行white檢驗，程序如下：

proc reg data=expend2013;

model expend=income / white;

run;

quit;

表3　等級相關系數(shù)檢驗結果

white檢驗結果見表4。

表4　white檢驗結果

由表4可以判斷存在異方差現(xiàn)象(P值<0.000 1)。

下面開始進行加權最小二乘的處理，分別采用殘差絕對值和殘差平方項作為權重，經(jīng)比較殘差平方項作權重的效果更好。利用殘差平方項加權的SAS程序如下：

procregdata=resid;

modelsqresid=income;

outputout=v_weightsp=v_hat;

run;

quit;

datav_weights;

setv_weights;

v_weight=/v_hat;

labelv_weight="weightsusingsquaredresiduals";

run;

procregdata=v_weights;

weightv_weight;

modelexpend=income/noint;

outputout=e2013_r_vr=residual;

labelincome="保費"expend="賠付";

run;

quit;

datar2013_abs_v;

sete2000_r_v;

abs_r_v=abs(residual2);

run;

proccorrdata=r2013_abs_vspearman;

varabs_r_vincome;

run;

加權之后，需要去掉常數(shù)項，加權之后的等級相關系數(shù)的結果是自變量收入與殘差絕對值|ei|的相關系數(shù)降為0.243 16，P值為0.111 75，可以判定異方差現(xiàn)象消除了。

表5　加權后的的參數(shù)估計結果

圖2回歸擬合圖

5結語

通過SAS軟件可以實現(xiàn)殘差圖法、等級相關系數(shù)法、white檢驗，還可以實現(xiàn)加權最小二乘法消除異方差。SAS軟件功能強大，處理異方差的方法的還有很多，SAS軟件在這個方面的應用還有待開發(fā)。

參考文獻：

[1]何曉群.應用回歸分析[M].北京:中國人民大學出版社,2011:96-116.

[2]龔秀芳.幾種異方差檢驗方法的比較[J].菏澤師范?？茖W校學報,2003(25):19-22.

[3]WhiteH.Aheteroskedasticity-consistentcovariancematrixestimatorandadirecttestforheteroskedastivity[J].Econometrica,1980(48):817-838.

[4]楊波.加權最小二乘估計中加權系數(shù)的確定[J].現(xiàn)代電子技術,2002(12):45-46.

[5]王純杰.基于分位數(shù)回歸的長春市職工工資水平的分析[J].長春工業(yè)大學學報:自然科學版,2010,31(4):367-373.

[6]董大鈞.SAS統(tǒng)計分析應用[M].北京：電子工業(yè)出版社,2008.

[7]朱世武.SAS編程技術教程[M].北京：清華大學出版社,2013.

[8]InstituteS.SAS/STATuser’sguide[M]. [S.l.]:SASInstitute,1990.

Inspection and elimination of heteroscedasticity based on SAS

WANG Chunjie,ZHENG Xi*,LI Qun,LIN Shanyi

(School of Basic Sciences, Changchun University of Technology, Changchun 130012, China)

Abstract：Heteroscedasticity phenomenon is tested with white test of REG procedure in SAS, and the weighted least square estimation is used to eliminate heteroscedasticity. The test is compared with both the residual graph and Spearman correlation coefficient method to show more reasonable results. The data of insurance and premium in 2013 is applied to the example analysis.

Key words：regression analysis; heteroscedasticity; weighted least square estimate.

收稿日期：2016-01-21

基金項目：國家自然科學基金資助項目(11301031，11571051)； 吉林省教育廳“十三五”規(guī)劃項目(2016317)

作者簡介：王純杰(1978-)，女，漢族，遼寧遼陽人，長春工業(yè)大學副教授，博士，主要從事數(shù)理統(tǒng)計方向研究,E-mail:wangchunjie@ccut.edu.cn. *通訊作者：鄭茜(1989-)，女，漢族，吉林遼源人，長春工業(yè)大學碩士研究生，主要從事數(shù)理統(tǒng)計方向研究,E-mail:13069008450@163.com.

DOI:10.15923/j.cnki.cn22-1382/t.2016.3.01

中圖分類號：O 212.1

文獻標志碼：A

文章編號：1674-1374(2016)03-0209-08