基于正交化方法的回歸分析

胡良平

(1.軍事科學院研究生院，北京　100850；2.世界中醫(yī)藥學會聯(lián)合會臨床科研統(tǒng)計學專業(yè)委員會，北京　100029

1　概　　述

1.1　何為基于正交化方法的回歸分析

基于正交化方法的回歸分析是在構(gòu)建多重線性回歸模型時，先對數(shù)據(jù)進行“Gentleman-Givens 變換”[1-2]，但仍基于“線性最小二乘原理”推導出的公式估計回歸系數(shù)。

1.2　基于正交化方法的回歸分析應用的場合

當擬做多重線性回歸分析的原始數(shù)據(jù)存在嚴重的“病態(tài)”時，采用此方法構(gòu)建多重線性回歸模型，可以最大限度地消除“病態(tài)數(shù)據(jù)”對建模結(jié)果造成的影響。

1.3　基于正交化方法的回歸分析的原理

此法使用“Gentleman-Givens變換”[3]校正數(shù)據(jù)，并且在計算數(shù)據(jù)矩陣的QR分解[4]的上三角矩陣R時十分謹慎。相對于其他正交化方法(例如Householder變換[3])，此法的優(yōu)點是不需要將數(shù)據(jù)矩陣存儲在計算機的內(nèi)存中。

2　基于正交化方法的回歸分析解決實際問題

2.1　此法可否解決多重共線性問題

2.1.1問題與數(shù)據(jù)結(jié)構(gòu)

沿用文獻[5]中的“問題與數(shù)據(jù)”，并基于派生變量得到的“最優(yōu)回歸模型”所決定的“數(shù)據(jù)集”，來提出下面的“新問題”：即“weight”的回歸系數(shù)為“-88.00801”，這個“負值”表明：體重越重的人收縮壓(SBP)越低，這似乎不符合臨床專業(yè)知識。盡管計算出來的因變量的預測值在專業(yè)上都成立，而且模型的殘差方差=122.32418、R2= 0.9931，這些結(jié)果都提示所構(gòu)建的多重線性回歸模型很好。但畢竟存在回歸系數(shù)的正負號不符合專業(yè)知識的“嚴重瑕疵”，這是一個需要徹底解決的“疑難問題”！

2.1.2所需要的SAS程序

嘗試采用“基于正交化方法”解決上述提及的“疑難問題”。所需要的SAS程序如下：

data a1;

input id age height weight bmi sbp;

cards;

(此處輸入文獻[5]表1中50行6列數(shù)據(jù))

;

run;

/* 以上程序為了創(chuàng)建數(shù)據(jù)集a1 */

data a2;

set a1;

x1=age*age; x2=age*height; x3=age*weight;

x4=age*bmi; x5=height*height; x6=height*weight;

x7=height*bmi; x8=weight*weight; x9=weight*bmi;

x10=bmi*bmi;

run;

/* 以上程序是在數(shù)據(jù)集a1基礎上創(chuàng)建數(shù)據(jù)集a2，它增添了10個派生變量 */

proc reg data=a2;

model sbp=age weight x3 x6-x10/noint;

quit;

/* 以上程序是調(diào)用REG過程并基于數(shù)據(jù)集a2擬合文獻[5]中那個“最佳”回歸模型 */

proc orthoreg data=a2;

model sbp=age weight x3 x6-x10/noint;

quit;

/* 以上程序是調(diào)用ORTHOREG過程并基于數(shù)據(jù)集a2擬合文獻[5]中那個“最佳”回歸模型 */

【SAS程序說明】在以上的SAS程序中，都用“/* ……… */”注釋語句作了說明。

2.1.3SAS輸出結(jié)果及其解釋

以上第1個SAS過程步(REG過程)程序的主要輸出結(jié)果如下：

參數(shù)估計值

以上第2個SAS過程步(ORTHOREG過程)程序的主要輸出結(jié)果如下：

變量自由度參數(shù)估計值標準誤差t值Pr > |t|age11.821817155899260.49294033783.700.0006weight1-88.008006763805126.676358039-3.300.0020x31-0.00970853699460.0034186421-2.840.0069x610.645688011242420.19305159293.340.0017x714.324558948331461.30917299623.300.0020x81-0.058353118670680.01883696-3.100.0035x910.785296033542440.25836051573.040.0041x101-2.62458498623560.8771470203-2.990.0046

比較上述由“REG過程”與“ORTHOREG過程”對同一個具有嚴重共線性問題的資料進行多重回歸分析的結(jié)果可知：后者比前者所估計的“參數(shù)值”更精細，但仍沒有消除多重共線性對回歸系數(shù)的嚴重影響(尤其是weight前的回歸系數(shù)為負值且其絕對值還比較大，不符合臨床專業(yè)知識)。也就是說：SAS/STAT中的“ORTHOREG過程”并不能解決“多重共線性問題”。

2.2　此法可否解決高階多項式曲線擬合問題

2.2.1問題與數(shù)據(jù)結(jié)構(gòu)

假定有一個高階多項式資料，見表1。

表1　某資料中首尾各10對數(shù)據(jù)(n=101)

注：表1中省略編號為11到91的數(shù)據(jù)

繪出表1資料中(x，y)各點的散布圖，見圖1。

圖1　表1資料中(x，y)散布圖

【問題】試擬合圖1中y依賴x變化的回歸模型。

2.2.2所需要的SAS程序

data Polynomial;

do i=1 to 101;

x=(i-1)/(101-1);

y=10**(9/2);

do j=0 to 8;

y=y * (x - j/8);

end;

output;

end;

run;

/* 以上程序是為了產(chǎn)生表1中的全部101對數(shù)據(jù) */

proc gplot data=Polynomial;

plot y*x;

run;

/* 以上程序是為了繪制圖1中的散布圖 */

proc glm data=Polynomial;

model y=x|x|x|x|x|x|x|x|x;

run;

/* 以上程序是為了調(diào)用GLM過程擬合九次多項式曲線回歸模型 */

ods graphics on;

proc orthoreg data=Polynomial;

effect xMod=polynomial(x / degree=9);

model y=xMod;

effectplot fit / obs;

store OStore;

run;

ods graphics off;

/* 以上程序是為了調(diào)用ORTHOREG過程擬合九次多項式曲線回歸模型并繪圖 */

2.2.3SAS輸出結(jié)果及其解釋

以下是“GLM過程步”輸出的計算結(jié)果：

參數(shù)估計值標準誤差t值Pr > |t|截距0.448960.1254793.580.0006x24.610626.0548524.060.0001x*x-443.7403493.388508-4.75<0.0001x*x*x2626.90806B642.9727184.09<0.0001x*x*x*x-7371.23677B2327.022377-3.170.0021x*x*x*x*x10697.73145B4741.5988972.260.0264x*x*x*x*x*x-7749.24904B5468.101939-1.420.1598x*x*x*x*x*x*x2214.08419B3328.3556610.670.5076x*x*x*x*x*x*x*x-0.00610B830.439899-0.001.0000x*x*x*x*x*x*x*x*x0.00000B---

以上結(jié)果表明：用GLM擬合多項式曲線回歸模型并不太合適，其中，標志“B”的那些行上的參數(shù)估計值不準確。

以下是基于GLM過程擬合結(jié)果繪制出的曲線圖，見圖2。

圖2中的實線是基于GLM過程計算的結(jié)果，而圓圈是實際觀測到的結(jié)果。不難看出，擬合效果很不理想。

以下是“ORTHOREG過程步”輸出的計算結(jié)果：

圖2　基于GLM過程擬合結(jié)果繪制出的曲線圖

參數(shù)自由度參數(shù)估計值標準誤差t值Pr > |t|Intercept1-2.17224978239E-115.841067E-12-3.720.0003x175.99773124392843.526134E-102.16E11<0.0001x^21-1652.407813618026.8407273E-9-242E9<0.0001x^3114249.45397693735.9828349E-82.38E11<0.0001x^41-64932.46157501272.8072627E-7-231E9<0.0001x^51173315.3593603037.6781594E-72.26E11<0.0001x^61-280158.0364593531.2614251E-6-222E9<0.0001x^71269781.8128871421.2252919E-62.2E11<0.0001x^81-142302.4947098696.4807927E-7-22E10<0.0001x^9131622.77660222611.4379253E-72.2E11<0.0001

以上結(jié)果表明：擬合的九次多項式曲線回歸模型中各參數(shù)均具有統(tǒng)計學意義。

以下是基于ORTHOREG過程擬合結(jié)果繪制出的曲線圖，見圖3。

圖3　基于ORTHOREG過程擬合結(jié)果繪制出的曲線圖

圖3中的實線是基于ORTHOREG過程計算的結(jié)果，而圓圈是實際觀測到的結(jié)果。不難看出，擬合效果非常好。

【說明】本例表1中的資料屬于一種相當嚴重的“病態(tài)”資料，由圖3可知，當x在(0.00，0.18)之間變化時，y的變化形成了一個“尖峰”；當x在(0.18，0.82)之間變化時，y的變化形成了近似平行于x軸的一條“略有起伏的波浪線”；而當x在(0.82，1.00)之間變化時，y的變化形成了一個“低谷”。