主成分回歸克服多重共線性的R語言實(shí)現(xiàn)

汪朋

摘要：多重共線性是回歸分析中容易出現(xiàn)的一類重要問題，現(xiàn)有的克服多重共線性的方法有很多，這其中主成分回歸是非常有效的一種，但該方法計(jì)算復(fù)雜，必須借助于計(jì)算軟件才能完成。為此，本文在已有R函數(shù)的基礎(chǔ)上，通過自編一定的R函數(shù)和代碼，探討了應(yīng)用R語言實(shí)現(xiàn)主成分回歸的過程。最后的案例表明，通過R語言實(shí)現(xiàn)主成分回歸來克服模型的多重共線性，過程簡(jiǎn)單，效果明顯，且容易被學(xué)習(xí)者和應(yīng)用者掌握。

關(guān)鍵詞：主成分回歸 多重共線性 R語言

中圖分類號(hào)： 文獻(xiàn)標(biāo)識(shí)碼：A 文章編號(hào)：1672-3791（2015）03（b）-0000-00

1引言

回歸模型是應(yīng)用最為廣泛的統(tǒng)計(jì)模型之一，常被用來處理多變量間相依關(guān)系問題，而多重共線性問題是回歸分析中一個(gè)非常棘手的問題，最早是R.Frisch提出。大量實(shí)踐表明，當(dāng)模型存在嚴(yán)重的多重共線性時(shí)，模型的參數(shù)估計(jì)變得不夠穩(wěn)定，甚至出現(xiàn)估計(jì)的回歸系數(shù)與其實(shí)際意義不相符，模型的顯著性檢驗(yàn)和預(yù)測(cè)功能失效等嚴(yán)重后果。正因如此，從20世紀(jì)六七十年代至今，人們提出了多種方法來克服模型的多重共線性，比較常用的有逐步回歸法、嶺回歸法和主成分回歸法等。逐步回歸法是剔除變量，容易造成信息損失和結(jié)構(gòu)分析的預(yù)期目的難以實(shí)現(xiàn)的問題；嶺回歸法要求選擇比較準(zhǔn)確的嶺參數(shù)，而嶺參數(shù)的選擇主觀性很強(qiáng)；相比之下，主成分回歸法通過將自變量分解成幾個(gè)不相關(guān)的主成分，能完全克服模型的多重共線性，且選擇的主成分能基本上概括原變量的信息，使得主成分回歸法成為克服多重共線性中的一種較為有效的方法，但該方法計(jì)算復(fù)雜，需要借助于計(jì)算軟件才能完成。目前主流的統(tǒng)計(jì)軟件雖一般提供了主成分分析的功能，都基本上沒有直接實(shí)現(xiàn)主成分回歸的模塊，因此，本文在R語言提供的主成分分析的函數(shù)的基礎(chǔ)上，通過自編相關(guān)的函數(shù)，來探討如何快速簡(jiǎn)便地運(yùn)用主成分回歸法克服回歸模型的多重共線性。

二、主成分回歸法的R語言實(shí)現(xiàn)過程

主成分回歸法是利用主成分分析找到解釋變量的若干個(gè)主成份，由于主成份之間是互不相關(guān)的，因此考慮將被解釋變量關(guān)于這些主成份進(jìn)行回歸，可以消除模型的多重共線性。

應(yīng)用R語言克服模型的多重共線性，具體步驟如下：

（1）運(yùn)用princomp（ ）函數(shù)進(jìn)行主成分分析，使用方法如下：

pr = princomp（～X1+X2+…4+Xk， cor=TRUE ）

這里cor =TRUE表示用樣本的相關(guān)矩陣R作主成分分析，cor=FALSE（默認(rèn)值），表示用樣本的協(xié)方差作主成分分析。

在這一步中還可以使用summmary（ ）函數(shù)來提取因子載荷矩陣等模型的各項(xiàng)信息。

（2）根據(jù)（1）中得到的結(jié)果選定主成分z1，z2，…，zm，然后將被解釋變量關(guān)于這些主成分進(jìn)行回歸，代碼如下：

pre=predict（pr）

z1=pre[，1]

z2=pre[，2]

…………

zm=pre[，m]

model.pr=lm（Y～z1+z2+…+zm）

summary（model.pr）

（3）運(yùn)用自編函數(shù)將主成分回歸模型轉(zhuǎn)換為關(guān)于原變量的回歸模型，自編函數(shù)的代碼如下：

prc.trans=function（formula， prc.object）{

object=lm（formula）； X=as.matrix（object$model）[，-1]

pr=princomp（～X， cor=TRUE）； A=loadings（pr）

apha=coef（prc.object）； n=length（apha）

if （n<=2）{

beta=A[，1：n-1]*apha[2：n] }

else {

beta=A[，1：n-1]%*%apha[2：n] }

xbar=apply（X，2，mean）； s=sd（X）

b=beta/s； b0=apha[1]-sum（b*xbar）

coefficients=c（b0，b）

coefficients=as.matrix（coefficients）

rownames（coefficients）=rownames（summary（object）$coefficients）

colnames（coefficients）="Estimate"

coefficients

}

該函數(shù)的prc.object參數(shù)為做主成分回歸所返回的對(duì)象。

有了上面的prc.trans（ ）函數(shù)，就可、可以很方便地將進(jìn)行系數(shù)轉(zhuǎn)換了，轉(zhuǎn)換代碼如下：

prc.trans（Y～X1+X2+…+Xk， model.pr）

三、案例分析

為了驗(yàn)證上節(jié)內(nèi)容的可行性，下面選取1991—2013年西藏第三產(chǎn)業(yè)及一些影響因素的數(shù)據(jù)來驗(yàn)證上節(jié)提出的方法和結(jié)論。

根據(jù)經(jīng)濟(jì)學(xué)的有關(guān)理論和經(jīng)驗(yàn)知，第三產(chǎn)業(yè)的發(fā)展主要與第三產(chǎn)業(yè)的規(guī)模、最終消費(fèi)、資本形成總額、本國(guó)交通運(yùn)輸規(guī)模等因素有關(guān)。為此考慮

——第三產(chǎn)業(yè)從業(yè)人員數(shù)（反映第三產(chǎn)業(yè)的規(guī)模）；

——財(cái)政支出（最終消費(fèi)的構(gòu)成因素）；

——全區(qū)居民消費(fèi)水平（最終消費(fèi)的影響因素）；

——全社會(huì)的固定資產(chǎn)投資（GDP的影響因素）；

——旅客周轉(zhuǎn)量（反映交通運(yùn)輸規(guī)模）。

對(duì)西藏第三產(chǎn)業(yè)總產(chǎn)值Y的影響，從而尋找影響西藏第三產(chǎn)業(yè)發(fā)展的關(guān)鍵因素。有關(guān)數(shù)據(jù)見表1。

表1 西藏1991—2012年第三產(chǎn)業(yè)產(chǎn)值及其影響因素的數(shù)據(jù)

年份

第三產(chǎn)業(yè)總值（億元）

第三產(chǎn)業(yè)從業(yè)人員數(shù)（萬人）

財(cái)政支出（萬元）

全區(qū)居民消費(fèi)水平（元/人）

全社會(huì)固定資產(chǎn)投資（萬元）

旅客周轉(zhuǎn)量（億人公里）

1991

10.03

16.67

129242

735

76105

3.90

1992

10.86

18.28

150018

839

105665

3.65

1993

12.24

19.47

166120

903

133297

4.08

1994

13.63

18.80

216012

931

181458

4.11

1995

16.97

21.92

302998

1110

211718

4.65

1996

19.39

19.96

348749

1202

369492

4.80

1997

26.46

22.14

381195

1312

303605

5.24

1998

31.13

23.08

390961

1471

345495

5.29

1999

39.99

24.08

461966

1551

427457

5.53

2000

47.86

25.26

544223

1669

566030

5.78

2001

54.37

25.85

616108

1823

665044

6.20

2002

69.65

28.52

1062067

1939

857725

6.42

2003

89.56

32.46

1398904

2725

1089868

6.71

2004

96.76

35.31

1481966

2825

1386165

6.49

2005

123.30

38.15

1360690

2950

1684361

7.09

2006

137.24

43.61

1891612

3019

1961916

18.40

2007

159.76

46.60

2023024

2990

2323503

22.00

2008

188.06

52.45

2793631

3215

2711811

27.00

2009

218.67

57.00

3840173

3504

3099304

30.39

2010

240.85

58.72

4711288

4027

3794158

30.00

2011

274.82

61.58

5625834

4326

4632585

42.94

2012

322.57

69.57

7756827

4730

5492690

45.17

2013

377.8

81.36

9339713

5340

7099822

47.29

數(shù)據(jù)來源：《西藏統(tǒng)計(jì)年鑒（2014）》

作出Y關(guān)于各解釋變量之間的散點(diǎn)圖如下：

圖1 被解釋變量Y與各解釋變量之間的散點(diǎn)圖

根據(jù)圖1，可以發(fā)現(xiàn)，Y與各變量之間都大致成線性的關(guān)系，因此總體回歸模型可設(shè)為

（8）

為檢驗(yàn)?zāi)Ｐ褪欠翊嬖诙嘀毓簿€性，可計(jì)算得到各解釋變量之間相關(guān)系數(shù)矩陣如下：

表2 各解釋變量之間的相關(guān)系數(shù)矩陣

X1

X2

X3

X4

X5

X1

1.0000000

0.9574153

0.9806276

0.9826848

0.9649556

X2

0.9574153

1.0000000

0.9282874

0.9906462

0.9600453

X3

0.9806276

0.9282874

1.0000000

0.9582267

0.9164908

X4

0.9826848

0.9906462

0.9582267

1.0000000

0.9707032

X5

0.9649556

0.9600453

0.9164908

0.9707032

1.0000000

由相關(guān)系數(shù)矩陣可以看出，各解釋變量相互之間的相關(guān)系數(shù)較高，證實(shí)原模型確實(shí)存在嚴(yán)重多重共線性。

下面在R語言中應(yīng)用主成分回歸法來克服模型的多重共線性。

首先，進(jìn)行主成分分析，運(yùn)用的函數(shù)為princomp（ ），對(duì)于本案例數(shù)據(jù)，代碼如下：

pr=princomp（～X1+X2+X3+X4+X5， cor=TRUE）

summary（pr，loadings=TRUE）

運(yùn)行結(jié)果如下：

Importance of components：

Comp.1 Comp.2 Comp.3 Comp.4 Comp.5

Proportion of Variance 0.9689 0.0198 0.0092 0.0015 0.0006

Cumulative Proportion 0.9689 0.9887 0.9979 0.9994 1.0000

從這個(gè)結(jié)果中可以看到，前兩個(gè)主成分的累計(jì)貢獻(xiàn)率已達(dá)到98.6%，包含了原5個(gè)變量的絕大部分信息，因此可以舍去另外三個(gè)主成分?？紤]Y與前兩個(gè)主成分之間的回歸，并將其轉(zhuǎn)換為Y關(guān)于原解釋變量的回歸，所需代碼如下：

pre=predict（pr）

z1=pre[，1]

z2=pre[，2]

model.pr=lm（Y～z1+z2）

prc.trans（Y～X1+X2+X3+X4+X5，model.pr）

運(yùn)行結(jié)果為：

prc.trans（Y～X1+X2+X3+X4+X5，model.pr）

運(yùn)行結(jié)果為：

Estimate

（Intercept） -4.049499e+01

X1 1.373448e+00

X2 6.272358e-06

X3 2.369798e-02

X4 1.038333e-05

X5 1.010673e+00

于是得到用主成分回歸法得到的最終模型為：

（11）

從模型的結(jié)果上看，得到的回歸系數(shù)基本上在合理的經(jīng)濟(jì)意義范圍內(nèi)。

四、結(jié)束語

主成分回歸法是近些年應(yīng)用得較多的一種新方法，相比于其他方法，能完全消除模型多重共線性，可以比較精準(zhǔn)的作出關(guān)于所有解釋變量的經(jīng)濟(jì)結(jié)構(gòu)分析，但該方法相對(duì)較復(fù)雜，尤其是在將關(guān)于主成分的回歸轉(zhuǎn)換為關(guān)于原變量的回歸上。因此，主成分回歸的實(shí)現(xiàn)不應(yīng)用計(jì)算軟件是很難實(shí)現(xiàn)的，而大多數(shù)軟件并沒有提供直接實(shí)現(xiàn)主成分回歸的功能，應(yīng)用時(shí)較為困難。本文在已有R語言函數(shù)的基礎(chǔ)上，通過自編系數(shù)轉(zhuǎn)換函數(shù)prc.trans（ ），較為完整地實(shí)現(xiàn)了主成分回歸的過程。最后的案例表明，應(yīng)用主成分回歸來克服模型的多重共線性，結(jié)果十分理想，而應(yīng)用R語言來實(shí)現(xiàn)主成分回歸，過程簡(jiǎn)單，結(jié)果明確，易學(xué)易懂，且通過對(duì)R語言實(shí)現(xiàn)主成分回歸的過程的把握和對(duì)相應(yīng)結(jié)果的理解，有助于計(jì)量經(jīng)濟(jì)學(xué)習(xí)者和應(yīng)用者把握主成分回歸背后的思想和原理，這是其他計(jì)算軟件，尤其是菜單式界面操作軟件所不能比擬的。

參考文獻(xiàn)

[1] 李子奈.計(jì)量經(jīng)濟(jì)學(xué)（第三版）[M].北京：高等教育出版社，2010：134-143.

[2] 何曉群.《多元統(tǒng)計(jì)分析》（第三版）[M].中國(guó)人民大學(xué)出版社，2004：114-125.

[3] 劉國(guó)旗.多重共線性的產(chǎn)生原因及其診斷處理[J].合肥工業(yè)大學(xué)學(xué)報(bào)：自然科學(xué)版，2001，24 （4）：607-610.

[4] 馬雄威.線性回歸方程中多重共線性診斷方法及其實(shí)證分析[J].華中農(nóng)業(yè)大學(xué)學(xué)報(bào)（社會(huì)科學(xué)版），2008，2：49-54.