胡良平
(1.軍事醫學科學院研究生院,北京 100850;2.世界中醫藥學會聯合會臨床科研統計學專業委員會,北京 100029
1.1.1 選擇合適變量變換方法的必要性
在進行回歸分析時,若自變量中有“多值名義變量”(如職業、血型、儀器品牌等),其具體的“表現或水平”不能用“文字”或“字母”表示,也不能簡單地賦值“1、2、3……”前者無法參與統計計算,而后者將會導致計算結果錯誤。那么,究竟應該對“多值名義變量”進行什么樣的變量變換呢?本文將介紹常規做法,即進行“啞變量變換”。
在回歸分析中應如何處置“多值有序變量”?在統計學上,人們認為:直接采用多值有序變量各水平的數值為其取值,例如:假定x代表“腫瘤分級”,依據臨床專業知識,已知它可以分為五級,于是,認為 x的取值就是“1、2、3、4、5”。依據基本常識可知,這樣的做法是不妥的。因為當腫瘤處于不同等級,其對結果的影響可能不是“線性關系”,很可能是較復雜的“非線性關系”。因此,應將“多值有序自變量”視為“多值名義自變量”,采用合適的變量變換方法。
1.1.2 對多值名義變量進行“啞變量變換”
所謂啞變量變換,就是將一個具有k個水平的多值名義變量轉換成(k-1)個新變量,每個新變量都是一個“二值變量(即僅有兩個不同取值的變量)”。這些新變量像“啞巴”一樣,其中的每一個都攜帶著原變量的一部分信息,在計算中發揮一定的作用,但又不能完全取代原變量,故它們都被形象地稱為“啞變量”。
實施啞變量變換的方法是:選擇一個頻率高的水平作為“基準水平”,其他水平都與該基準水平作比較而產生一個“比較變量”(即啞變量)。例如:在ABO血型系統中,假定在樣本資料中屬于O型血的人數最多,就可以以“O型血的人”為“基準水平”,其他三種血型的人相對于O型血的人分別產生一個“啞變量”。簡化形式呈現如下:

個體編號血型XA|OXB|OXAB|O1A1002B0103AB0014O000
在上面的簡化形式中,“XA|O、XB|O、XAB|O”這三個變量都是與“血型”這個4值名義變量對應的“啞變量”,它們分別代表“是否為A型血”“是否為B型血”和“是否為AB型血”。
1.1.3 對多值名義變量進行“其他變量變換”
在進行回歸分析中,上面的“啞變量變換”已經成為統計學界處置“多值名義自變量”的“金標準”。是否還有更合理的“變量變換”方法可以取代“啞變量變換”呢?此問題將在本期“科研方法專題”的另三篇文章中深入討論。
1.2.1 選擇合適變量變換方法的必要性
通常情況下,人們在進行回歸分析時,對于定量的自變量和/或因變量不作任何變換。然而,由基本常識可知,前述做法是不切實際的,通常情況下,效果是不夠好的。因為變量之間的關系往往是錯綜復雜的,它們之間永遠以“一次方”形式存在聯系的可能性是非常罕見的。因變量Y可能與某個自變量之間是拋物線關系、指數曲線關系或對數曲線關系;因變量Y本身可能偏離正態分布很遠,而很多統計模型要求因變量必須服從正態分布。因此,需要對定量因變量作合適的變量變換,以使其符合特定統計模型的基本要求;需要對某些定量自變量作合適的變量變換,以更真實地呈現其與定量因變量之間的變化趨勢。
1.2.2 對定量自變量進行兩方面的變量變換
第一方面的變量變換就是對某定量自變量作了某種變量變換后,丟棄原先的那個自變量,而僅采用變換后的變量。例如:建模時,只用“log(x1)”,而丟棄“x1”。第二方面的變量變換就是不僅用變換后的變量,還保留未變換的原變量。這樣做的結果會使自變量的數目大大增加,常稱為產生“派生變量”。例如:假定有10個定量變量,可以給它們都取對數變換,就會增加10個新變量;也可以對10個變量進行平方變換或平方根變換;還可以基于10個定量變量產生交叉乘積項等。
1.2.3 對定量因變量進行變量變換
在通常情況下,人們進行的是“一元多重回歸分析”,因此,若對定量因變量進行變量變換,在回歸建模時,只使用變換后的因變量,而不會同時使用原先的“因變量”與變換后的因變量(因為這樣做已經把“一元”問題轉變成“二元”問題了)。
何時需要對定量因變量進行變換呢?通常在以下兩種情況之一:其一,已知因變量與自變量之間呈某種函數關系,就選擇相應的變量變換方法。例如:當因變量與自變量之間呈“指數函數”變化關系時,就可以對因變量取對數變換;其二,當定量因變量(嚴格地說,應該是模型的誤差項)偏離正態分布很遠時,需要選擇一種合適的變量變換方法,目的是使變換后的因變量服從模型所要求的某種概率分布,如正態分布、指數分布或威布爾分布等。
研究者關心的定量結果變量為“氧化氮釋放量(nox)”,該定量指標的數值測自單缸發動機。已知影響因素有:燃油種類(fuel)、壓縮比(cpratio)和等值比(eqratio)。其中,燃油種類(fuel)是多值名義變量,而氧化氮釋放量(nox)、壓縮比(cpratio)和等值比(eqratio)都是計量變量。該資料來自SAS軟件中的“幫助”數據庫,數據集名為:sashelp.gas。
試以“氧化氮釋放量(nox)”為因變量,以“燃油種類(fuel)、壓縮比(cpratio)和等值比(eqratio)”為自變量,創建一元多重回歸模型。
【說明】該實際問題和對應的數據來源于“SAS/STAT的TRANSREG過程中的樣例及SASHELP數據庫,其數據集名為sashelp.gas”[1]。
利用以下SAS程序可以顯示該例的數據結構:
proc print data=sashelp.gas;
run;
【燃油資料的數據結構】

ObsFuelCpRatioEqRatioNOx1Ethanol120.9073.7412Ethanol120.7612.2953Ethanol121.1081.4984Ethanol121.0162.8815Ethanol121.1890.760
以上顯示出數據集的前5個觀測,全部資料共171個觀測。其中,在結果變量nox上有兩個缺失值。
利用如下SAS程序可以顯示三個自變量(一個為多值名義自變量、一個為多值有序自變量、一個為定量自變量)及定量結果變量(nox)的頻數分布情況:
proc freq data=sashelp.gas;
tables fuel eqratio cpratio nox;
run;
【燃油種類的頻數分布】

Fuel頻數百分比累積頻數累積百分比82rongas95.2695.2694%Eth2514.623419.88Ethanol9052.6312472.51Gasohol137.6013780.12Indolene2212.8715992.98Methanol127.02171100.00
以上結果表明:共有6種燃油,其中,頻數最多的是“Ethanol”,涉及此種燃油的觀測共有90個。
【壓縮比的頻數分布】

Compression RatioCpRatio頻數百分比累積頻數累積百分比7.59354.399354.399179.9411064.33122414.0413478.36152011.7015490.0618179.94171100.00
以上結果表明:壓縮比只有5種,屬于“多值有序”變量(注意:以下簡稱為“定量變量”)。其中,頻數最多的是“7.5”,涉及此種壓縮比的觀測共有93個。
等值比(eqratio)與氧化氮釋放量(nox)的取值都很多,其頻數分布表此處從略;但利用下面的SAS程序可以顯示這兩個變量的頻數分布直方圖,同時,還可以對它們進行正態性檢驗:
proc univariate data=sashelp.gas normal;
var eqratio nox;
histogram eqratio nox/normal;
run;
【等值比的正態性檢驗結果】

正態性檢驗檢驗統計量PShapiro-WilkW0.969774Pr
以上結果表明:等值比不服從正態分布。
等值比的頻數分布直方圖見圖1。由圖1可知,等值比呈“負偏態分布”
【氧化氮釋放量的正態性檢驗結果】

正態性檢驗檢驗統計量PShapiro-WilkW0.945485Pr
以上結果表明:氧化氮釋放量不服從正態分布。
氧化氮釋放量的頻數分布直方圖見圖2。由圖2可知:氧化氮釋放量呈“正偏態分布”。

圖1 等值比的頻數分布直方圖 圖2 氧化氮釋放量的頻數分布直方圖
選擇出現頻數最多的水平“Ethanol”為“基準”,產生5個啞變量:g1到g5。實現此任務的SAS程序如下:
data a1;
set sashelp.gas;
g1=0;g2=0;g3=0;g4=0;g5=0;
if fuel=' 82rongas' then g1=1;
else if fuel=' 94%Eth' then g2=1;
else if fuel=' Gasohol' then g3=1;
else if fuel=' Indolene' then g4=1;
else if fuel=' Methanol' then g5=1;
run;
g1到g5分別代表:“是否為82rongas燃油”“是否為94%Eth燃油”“是否為Gasohol燃油”“是否為Indolene燃油”和“是否為Methanol燃油”。
在數據集a1基礎上增加由定量自變量派生出來的13個新變量,產生數據集a2。SAS程序如下:
data a2;
set a1;
x1=eqratio**2;x2=eqratio*cpratio;
x3=cpratio**2;x4=x1*eqratio;
x5=x3*cpratio;x6=x1*cpratio;
x7=x3*eqratio;x8=sqrt(eqratio);
x9=sqrt(cpratio);x10=log(eqratio);
x11=log(cpratio);x12=exp(eqratio);
x13=exp(cpratio);
run;
【說明】cpratio和eqratio是資料中兩個原始的定量自變量;x1、x4、x8、x10、x12分別是“eqratio”的平方變換、立方變換、平方根變換、自然對數變換和指數變換的結果;x3、x5、x9、x11、x13分別是“cpratio”的平方變換、立方變換、平方根變換、自然對數變換和指數變換的結果;x2是“eqratio”與“cpratio”的交叉乘積項;x6是“eqratio”的平方項與“cpratio”的交叉乘積項;而x7是“cpratio”的平方項與“eqratio”的交叉乘積項。
在數據集a2基礎上同時增加定量因變量的對數變換y1、平方根變換y2、指數變換y3、倒數變換y4和Logistic變換y5,產生數據集a3。SAS程序如下:
data a3;
set a2;
y1=log(nox);y2=sqrt(nox);y3=exp(nox);
y4=1/nox;y5=exp(nox)/(1+exp(nox));
run;
對一個“多值名義自變量”采取“啞變量變換”,以其為基礎,再分別選取定量因變量(nox)的6種不同“表現”為每次建模的“因變量”,并對定量自變量在“不做變量變換”和“引入13個派生變量”且分別在回歸模型中假定“包含截距項”與“不含截距項”的條件下,采取“前進法”“后退法”和“逐步法”篩選自變量。
定量因變量(nox)的6種不同“表現”分別是:①定量因變量(nox),即對“定量因變量(nox)”不做變量變換;②定量因變量[y1=log(nox)],即對“定量因變量(nox)”做自然對數變換;③定量因變量[y2=SQRT(nox)],即對“定量因變量(nox)”做平方根變換;④定量因變量[y3=exp(nox)],即對“定量因變量(nox)”做指數變換;⑤定量因變量(y4=1/nox),即對“定量因變量(nox)”做倒數變換;⑥定量因變量{y5=exp(nox)/[1+exp(nox)]},即對“定量因變量(nox)”做Logistic變換。
在定量因變量(nox)每種“表現”且分別在定量自變量“不做變換”與“引入派生變量”的條件下,再在回歸模型中假定“包含截距項”與“不含截距項”時,分別采取“前進法”“后退法”和“逐步法”篩選自變量。這實際上就有“2×2×3=12”個回歸模型,它們分屬于4種情形:①“定量自變量不做變換”且假定“包含截距項”;②“定量自變量不做變換”且假定“不含截距項”;③“定量自變量做變換”且假定“包含截距項”;④“定量自變量做變換”且假定“不含截距項”。每種情形都涉及3種篩選自變量的方法,最多有3種不同的回歸模型,從中選取一個擬合最好的回歸模型。
所以,在每種特定的因變量條件下,就對應著4個“最優回歸模型”;故在因變量的6種條件下,一共有24個“最優回歸模型”。見表1。

表1 反映24個多重回歸模型擬合優度的計算結果
續表1:

第5組模型:對定量因變量做倒數變換170.0891 0.0781 0.37187 2.5112 2 有180.58300.5780 0.37923 2.2523 2 無190.8285 0.8199 0.0726513.5416 8 有200.9243 0.9185 0.0732018.7606 12 無第6組模型:對定量因變量做Logistic變換210.0856 0.0746 0.01436 7.06592 2 有220.9545 0.9525 0.03543 7.00000 7 無230.9539 0.9504 0.0007715.4067 12 有240.9991 0.9990 0.0007616.1852 16 無
注:第1組模型對應的因變量為“氧化氮釋放量(nox)”;第2組模型對應的因變量為“氧化氮釋放量的自然對數變換結果(y1)”;第3組模型對應的因變量為“氧化氮釋放量的平方根變換結果(y2)”;第4組模型對應的因變量為“氧化氮釋放量的指數變換結果(y3)”;第5組模型對應的因變量為“氧化氮釋放量的倒數變換結果(y4)”;第6組模型對應的因變量為“氧化氮釋放量的Logistic變換結果(y5)”
一般來說,當模型中包含的自變量數目相等且都包含截距項或都不含截距項時,R2值越大越好;此時,Cp值越接近自變量個數越好;當保留在模型中的自變量個數相差較多時,在前述判斷方法基礎上,再加上“均方誤差”(越小越好)和“調整R2”(越大越好),則更好。
5.2.1 第1組模型的擬合效果評價
第1組模型對應的因變量為“氧化氮釋放量”,模型1與模型2都是基于“5個啞變量加上2個定量自變量”進行變量篩選,其區別在于模型1假定包含截距項,而模型2假定不含截距項;模型3與模型4都是基于“5個啞變量加上2個定量自變量及其13個派生變量”進行變量篩選,其區別在于模型3假定包含截距項,而模型4假定不含截距項。由表1中前4行結果可知:模型2優于模型1、模型4優于模型3,即在相同情況下,假定不含截距項的擬合結果優于假定包含截距項的擬合結果;進一步比較可知:模型4優于模型2,即引入派生變量的擬合結果優于不引入派生變量的擬合結果。
5.2.2 第2組模型的擬合效果評價
第2組模型對應的因變量為“氧化氮釋放量的自然對數變換結果(y1)”,模型5與模型6都是基于“5個啞變量加上2個定量自變量”進行變量篩選,其區別在于模型5假定包含截距項,而模型6假定不包含截距項;模型7與模型8都是基于“5個啞變量加上2個定量自變量及其13個派生變量”進行變量篩選,其區別在于模型7假定包含截距項,而模型8假定不包含截距項。由表1中第5~8行結果可知:模型6優于模型5、模型8優于模型7,即在相同情況下,假定不含截距項的擬合結果優于假定包含截距項的擬合結果;進一步比較可知:模型8優于模型6,即引入派生變量的擬合結果優于不引入派生變量的擬合結果。
5.2.3 第3組模型的擬合效果評價
第3組模型對應的因變量為“氧化氮釋放量的平方根變換結果(y2)”,模型9與模型10都是基于“5個啞變量加上2個定量自變量”進行變量篩選,其區別在于模型9假定包含截距項,而模型10假定不包含截距項;模型11與模型12都是基于“5個啞變量加上2個定量自變量及其13個派生變量”進行變量篩選,其區別在于模型11假定包含截距項,而模型12假定不包含截距項。由表1中第9~12行結果可知:模型10優于模型9、模型12優于模型11,即在相同情況下,假定不含截距項的擬合結果優于假定包含截距項的擬合結果;進一步比較可知:模型12優于模型10,即引入派生變量的擬合結果優于不引入派生變量的擬合結果。
5.2.4 第4組模型的擬合效果評價
第4組模型對應的因變量為“氧化氮釋放量的指數變換結果(y3)”,模型13與模型14都是基于“5個啞變量加上2個定量自變量”進行變量篩選,其區別在于模型13假定包含截距項,而模型14假定不包含截距項;模型15與模型16都是基于“5個啞變量加上2個定量自變量及其13個派生變量”進行變量篩選,其區別在于模型15假定包含截距項,而模型16假定不包含截距項。由表1中第13~16行結果可知:模型14優于模型13、模型16優于模型15,即在相同情況下,假定不含截距項的擬合結果優于假定包含截距項的擬合結果;進一步比較可知:模型16優于模型14,即引入派生變量的擬合結果優于不引入派生變量的擬合結果。
5.2.5 第5組模型的擬合效果評價
第5組模型對應的因變量為“氧化氮釋放量的倒數變換結果(y4)”,模型17與模型18都是僅基于“3個定量自變量”進行變量篩選,其區別在于模型17假定包含截距項,而模型18假定不包含截距項;模型19與模型20都是基于“3個定量自變量及其18個派生變量”進行變量篩選,其區別在于模型19假定包含截距項,而模型20假定不包含截距項。由表1中第17~20行結果可知:模型18優于模型17、模型20優于模型19,即在相同情況下,假定不含截距項的擬合結果優于假定包含截距項的擬合結果;進一步比較可知:模型20優于模型18,即引入派生變量的擬合結果優于不引入派生變量的擬合結果。
5.2.6 第6組模型的擬合效果評價
第6組模型對應的因變量為“氧化氮釋放量的Logistic變換結果(y5)”,模型21與模型22都是僅基于3個定量自變量進行變量篩選,其區別在于模型21假定包含截距項,而模型22假定不包含截距項;模型23與模型24都是基于3個定量自變量及其18個派生變量進行變量篩選,其區別在于模型23假定包含截距項,而模型24假定不包含截距項。由表1中第21~24行結果可知:模型22優于模型21、模型24優于模型23,即在相同情況下,假定不含截距項的擬合結果優于假定包含截距項的擬合結果;進一步比較可知:模型24優于模型22,即引入派生變量的擬合結果優于不引入派生變量的擬合結果。
5.2.7各組模型中最優模型擬合優度總評價
從以上的“評價結果”可知:模型4、模型8、模型12、模型16、模型20和模型24分別是從6組模型中挑選出來的“最優模型”,現將它們從表1中摘錄出來,以便直觀比較和判斷。見表2。

表2 各組挑選出來的6個“最優”多重回歸模型擬合優度的計算結果
由表2可知:模型24是6個“最優”模型中“最佳”的。該模型的因變量為“氧化氮釋放量(nox)的Logistic變換結果(y5)”,從全部(5+2+13=20個)自變量中篩選出了16個具有統計學意義的自變量,模型中不含截距項。具體計算結果如下:

方差分析源自由度平方和均方FPr > F模型16126.039047.8774410431.6<0.0001誤差1530.115540.00075515未校正合計169126.15458

變量參數估計值標準誤差II 型 SSFPr > Fg10.053670.009990.0217928.86<0.0001g30.060210.008660.0365048.33<0.0001g40.059570.007130.0527569.85<0.0001EqRatio2915.10929665.612980.0144819.18<0.0001CpRatio-932.92081221.155370.0134417.79<0.0001x1-591.67642128.766190.0159421.11<0.0001x2-0.096580.043720.003694.880.0287x329.765427.057010.0134317.79<0.0001x493.2702919.619880.0170722.60<0.0001x5-0.558000.132320.0134317.78<0.0001x60.072000.019400.0104013.770.0003x7-0.002270.001120.003074.070.0454x8-5597.266531310.349450.0137818.25<0.0001x93191.32640756.404990.0134417.80<0.0001x10785.94985188.216330.0131717.44<0.0001x136.991368E-71.657578E-70.0134317.79<0.0001
輸出以上結果的“SAS過程步程序”如下:
/*模型24:R2=0.9991,調整R2=0.9990,MSE=0.00075515,Cp=16.1852,niv=16,無截距項*/
proc reg data=a3;
model y5=g1-g5 eqratio cpratio x1-x13/noint
selection=backward sls=0.05 r;
/*模型24*/
run;
應注意:全部啞變量共有5個(它們之間不是互相對立的),采用篩選自變量的方法,保留下來其中的3個。嚴格地說,由一個多值名義自變量產生的全部啞變量應當同時被保留在回歸模型中或同時被排除出回歸模型,但這兩種結局都存在局限性;而將有關聯性的5個啞變量視為“獨立”的,根據假設檢驗結果保留其中的3個,這個結果也存在弊端。如何更妥善地處置“多值名義自變量”,將在本期科研方法專題后續文章中繼續討論。