胡良平

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

1 回歸建模前的變量變換[1]

1.1 對多值名義自變量進行“優(yōu)化計分變換(Opscore)”或“單調(diào)變換(Monotone)”

1.1.1優(yōu)化計分變換

優(yōu)化計分變換就是采用Fisher[2]提出的給多值名義變量或多值有序變量各水平賦值的方法。Opscore(x)代表對變量x進行優(yōu)化計分變換，它可以被用于“多值名義變量”或“多值有序變量”。

1.1.2 單調(diào)變換

單調(diào)變換就是采用Kruskal提出的給多值有序變量各水平賦值的方法(被稱為次要最小平方單調(diào)變換)。Monotone(x)代表對變量x進行單調(diào)變換，它只能被用于“多值有序變量”。

1.1.3 “優(yōu)化計分變換”與“單調(diào)變換”的具體做法

兩者具體的變換方法是相同的，每一個不同的非缺失值形成一個不同的類，例如：“1,1,1,2,2,3”形成三類，即3個“1”算作一類、兩個“2”算作一類、一個“3”算作一類。具體賦值方法參見下面的實例：

設(shè)變量x的具體取值為( . . .A .A .B 1 1 1 2 2 3 3 3 4)' ；設(shè)變量y的具體取值為( 5 6 2 4 2 1 2 3 4 6 4 5 6 7)' 。

【說明】以上“x”與“y”都有14個取值，用“向量”形式表示，即這兩個向量都有14個分量。其中，變量x的前兩個分量都是缺失值“.”；第3和第4個分量都是“.A”；第5個分量是“.B”；第6～14個分量都是數(shù)字。而變量y的全部14個分量都是數(shù)字。

于是，以“y”為“因變量”，對“變量x”進行“優(yōu)化計分變換”和“單調(diào)變換”。用關(guān)鍵詞表示為：Opscore(x)和Monotone(x)。它們只有8個分量，即：

Opscore(x)=( 5 6 3 2 2 5 5 7)'

Monotone(x)=( 5 6 3 2 2 5 5 7)'

【說明】在“x”的向量中，最前面的兩個缺失值各形成一類，即有兩類；“.A”出現(xiàn)了兩次，形成一類；“.B”形成一類 ；三個“1”形成一類；兩個“2”形成一類；三個“3”形成一類；一個“4”形成一類，故總共形成8類。

在同一類中，求出“對應(yīng)于因變量y全部數(shù)值的算術(shù)平均值”，賦值給變換后的“變量x”，作為其“優(yōu)化計分變換”或“單調(diào)變換”的結(jié)果，例如：兩個“.A”對應(yīng)著變量y的兩個數(shù)值為2和4，其算術(shù)平均值為3；同理，三個“1”對應(yīng)著變量y的三個數(shù)值為“1、2、3”，其算術(shù)平均值為2；兩個“2”對應(yīng)著變量y的兩個數(shù)值為“4與6”，其算術(shù)平均值為5；三個“3”對應(yīng)著變量y的三個數(shù)值為“4、5、6”，其算術(shù)平均值為5；一個“4”對應(yīng)著變量y的一個數(shù)值為“7”，其算術(shù)平均值為7。

1.2 可對定量變量進行的變量變換的種類

1.2.1 變量擴展

1.2.2 非優(yōu)化變換

非優(yōu)化變換就是用一個變換后的“新變量”取代“原變量”。這個“新變量”在后續(xù)的迭代算法中不再被重新變換了(對具有缺失值估計所做的可能的線性變換除外)。常用的“非優(yōu)化變換”有如下幾種。

ARSIN(x)：此處的“x”為“百分率”數(shù)據(jù)，通常其取值區(qū)間為“0

EXP(x)：指數(shù)變換，通常為“ex”變換；若希望底數(shù)為任意的正實數(shù)a，即做“ax”變換，在使用SAS的“TRANSREG過程”中，需要另外再使用“Parameter=”來指定，即在“=”后寫出“a”的具體數(shù)值。

LOG(x)：對數(shù)變換，通常取以“e”為底數(shù)的對數(shù)(e=2.718281828)，即取自然對數(shù)變換；若希望底數(shù)為任意的正實數(shù)a，即做“l(fā)ogax”變換，在使用SAS的“TRANSREG過程”中，需要另外再使用“Parameter=”來指定，即在“=”后寫出“a”的具體數(shù)值。

LOGIT(x)：logit變換，即做“l(fā)oga[x/(1-x)]”變換，這里的“a”通常為“e=2.718281828”。注意：x為數(shù)值型變量，其定義域為(0.0，1.0)。

POWER(x)：冪變換，即做“xa”變換。例如：希望做“x1.5”變換，關(guān)鍵詞寫為POWER(x/PARAMETER=1.5)。顯然，當(dāng)“PARAMETER=”后數(shù)值為“2”時，就是平方變換；為“-1”時，就是倒數(shù)變換；為“0.5”時，就是開平方根變換。

RANK(x)：秩變換。將變量x由小到大排序，再分別賦值為“1、2、3、4、5、……”。

1.2.3 非線性擬合變換

BOXCOX(x)：這是Box等[5]提出的變量變換方法，實際應(yīng)用參見文獻[6]。此變換只適用于回歸分析中的“定量因變量”。

PBSPLINE(x)：這是非迭代懲罰B樣條變換，使用時的特殊要求參見文獻[1]。此變換只適用于回歸分析中的“定量自變量”。

SMOOTH(x)：這是一種非迭代光滑樣條變換，由Reinsch[7]提出的變量變換方法，使用時的特殊要求參見文獻[1]。此變換只適用于回歸分析中的“定量自變量”。

1.2.4 優(yōu)化變換

優(yōu)化變換就是采用迭代計算導(dǎo)出的變量變換方法。在SAS的“TRANSREG過程”中，涉及到以下六種優(yōu)化變換方法，即：LINEAR(x)、MONOTONE(x)、MSPLINE(x)、OPSCORE(x)、SPLINE(x)和UNTIE(x)。以上六種優(yōu)化變換在使用時都有具體要求，詳見文獻[1]，此處從略。

1.3 對定量變量進行其他變量變換

1.3.1 恒等變換[IDENTITY(x)]

恒等變換，就是沒有改變變量的取值。通常情況下，就是“不做變換”。應(yīng)注意：在使用SAS的“TRANSREG過程”時，寫在“model語句”中的所有“因變量”和“自變量”之前必需要有“變量變換”的關(guān)鍵詞。若不想對“x4-x10”進行任何變量變換，必須按如下的方式寫：IDENTITY(x4-x10)。

使用此種變換，可以產(chǎn)生“交互項”，具體使用方法參見文獻[1] ，此處從略。

1.3.2 迭代光滑樣條變換[SSPLINE(x)]

這是對變量x進行迭代光滑樣條變換，此變量變換方法通常不會使誤差平方達到最小值。具體使用方法參見文獻[1] ，此處從略。

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

沿用本期科研方法專題第一篇文章《提高回歸模型擬合優(yōu)度的策略(Ⅰ)——啞變量變換與其他變量變換》中的“實際問題與數(shù)據(jù)結(jié)構(gòu)”[1]，此處從略。

3 解決問題的思路和做法

3.1 對多值名義自變量和多值有序自變量進行變量變換

3.1.1所需要的SAS程序

對“燃油種類(fuel)” 這個“6值名義自變量”進行“優(yōu)化計分變換”；對“壓縮比(cpratio)” 這個“5值有序自變量”進行“單調(diào)變換”。所需要的SAS程序如下：

ods graphics on;

title' Gasoline Example';

title2' Iteratively Estimate NOx, CpRatio, EqRatio, and Fuel';

* Fit the Nonparametric Model;

proc transreg data=sashelp.Gas solve test nomiss plots=all;

ods exclude where=(_path_ ? ' MV' );

model mspline(NOx / nknots=9) =

spline(EqRatio / nknots=9)

monotone(CpRatio) opscore(Fuel);

output out=aaa tdprefix=td tiprefix=ti;

run;

proc freq data=aaa;

tables cpratio ticpratio fuel tifuel;

run;

3.1.2 與所需變量變換有關(guān)的SAS輸出結(jié)果

Compression RatioCpRatio頻數(shù)百分比累積頻數(shù)累積百分比7.59354.399354.399179.9411064.33122414.0413478.36152011.7015490.0618179.94171100.00

以上是“壓縮比(cpratio)”的原始取值及其各水平的頻數(shù)分布。

Compression Ratio TransformationtiCpRatio頻數(shù)百分比累積頻數(shù)累積百分比7.15982275949354.399354.3910.953020684179.9411064.3312.9445525662414.0413478.3614.201537262011.7015490.0617.433541388179.94171100.00

以上是“壓縮比(cpratio)”經(jīng)過“單調(diào)變換”后的取值及其各水平的頻數(shù)分布。

由以上兩部分輸出結(jié)果可知：“單調(diào)變換”是將“多值有序變量”的水平做一一對應(yīng)的變換，各水平出現(xiàn)的“頻數(shù)”不會發(fā)生變化。

Fuel頻數(shù)百分比累積頻數(shù)累積百分比095.2695.2612514.623419.8829052.6312472.513137.6013780.1242212.8715992.985127.02171100.00

以上是“燃油種類(fuel)”的“6種水平代碼(分別為0、1、2、3、4、5)”及其各水平的頻數(shù)分布。

Fuel頻數(shù)百分比累積頻數(shù)累積百分比82rongas95.2695.2694%Eth2514.623419.88Ethanol9052.6312472.51Gasohol137.6013780.12Indolene2212.8715992.98Methanol127.02171100.00

以上是“燃油種類(fuel)”的“6種水平的真實含義”及其各水平的頻數(shù)分布。

Fuel TransformationtiFuel頻數(shù)百分比累積頻數(shù)累積百分比1.3742357556127.02127.021.55720545762514.623721.641.59754088319052.6312774.274.209335677695.2613679.534.36641535932212.8715892.404.4654059829137.60171100.00

以上是“燃油種類(fuel)” 經(jīng)過“優(yōu)化計分變換”后的取值及其各水平的頻數(shù)分布。

3.2 對定量自變量進行變量變換

對定量自變量可以進行“樣條變換(spline)”“單調(diào)樣條變換(mspline)”或“變量擴展變換(pspline)(即產(chǎn)生‘派生變量’，例如變量的平方項、立方項等)”。例如：對“eqratio”進行“樣條變換”，可以寫成spline(eqratio)或spline(eqratio/degree=3 nknots=0)。后一種寫法涉及到兩個參數(shù)，一是“degree=”，它是指“樣條函數(shù)(通常是多項式)”中多項式的“次數(shù)”，“degree=3”代表“三次多項式”；二是“nknots=”，它是指“結(jié)點(或?qū)懗伞?jié)點’)個數(shù)”，實際上就是“斷點個數(shù)”。其具體含義是：在該定量變量的取值區(qū)間內(nèi)，確定“幾個斷點(如nknots=9)”，于是，就將該區(qū)間劃分成“10段”。在每一段上，擬合一個“多項式”。就整體而言，被稱為“分段多項式(a piecewise polynomial)”。

3.3 對定量因變量進行變量變換

若需要變換的定量變量為“因變量”，除了可以進行“定量自變量”的某些變量變換(注意：“非迭代懲罰B樣條變換”只能用于“定量自變量”)之外，還可以進行其他一些變量變換，如“BOX-COX變換(注意：該變換僅適用于‘定量因變量’)”等。

3.4 針對實際問題，尋找合適的變量變換方法

3.4.1尋找最優(yōu)回歸模型的策略

前面提及的“實際問題”涉及“定量因變量(nox)”“定量自變量(eqratio)”“多值有序自變量(cpratio)”和“多值名義自變量(fuel)”。在擬合多重回歸模型的過程中，涉及到如何對上述4個變量進行“變量變換”。當(dāng)然，對每個變量選取不同的“變量變換”方法，將會得到不同的擬合效果。下面將固定“多值名義自變量(fuel)”的變量變換方法為：優(yōu)化計分變換，即“OPSCORE(fuel)”，對其他變量依次采取各種可能的變量變換，呈現(xiàn)出各種情況下回歸模型的擬合優(yōu)度，最終選擇“擬合優(yōu)度最好”的回歸模型。

3.4.2不同變量變換方法組合下的回歸模型擬合優(yōu)度

在對“多值名義自變量(fuel)”做優(yōu)化計分變換“OPSCORE(fuel)”的前提條件下，對其他變量依次采取各種可能的變量變換，對應(yīng)的回歸模型擬合優(yōu)度見表1。

表1 在OPSCORE(fuel)的前提下對“nox”“eqratio”和“cpratio”進行不同變量變換后對應(yīng)的回歸模型的擬合優(yōu)度

注：“Root MSE”代表“均方誤差的平方根”，此值越小越好；“Nox_B”代表對“Nox”進行變量變換的方法，表中第3列和第4列符號代表的含義與“Nox_B”相同，此處從略

由表1可知：模型1和模型7的結(jié)果相同且擬合優(yōu)度最好，對應(yīng)的SAS過程步程序如下：

proc transreg data=sashelp.Gas solve test nomiss

plots=all;

ods exclude where=(_path_ ?' MV' );

model spline(NOx / nknots=9)

=spline(EqRatio / nknots=9)

monotone(CpRatio) opscore(Fuel);

run;

3.4.3 模型1的主要輸出結(jié)果

由于“模型1”是一個“非參數(shù)模型”，只能給出總模型的有關(guān)信息和與“擬合優(yōu)度”有關(guān)的結(jié)果，不能給出“回歸系數(shù)”等信息。但可以以“圖形”方式給出各變量變換的結(jié)果，以提示對不同變量采取什么變量變換方法最有效。下面先給出與“擬合優(yōu)度”有關(guān)的結(jié)果：

Root MSE0.30935R-Square0.9586Dependent Mean2.34593Adj R-Sq0.9527Coeff Var13.18661

表1中的第1行就是摘錄了以上的結(jié)果。

在前面給出的“SAS過程步”的“MODEL語句”中，對定量因變量(nox)、定量自變量(EqRatio)、多值有序自變量(CpRatio)和多值名義自變量(Fuel)分別做了“樣條變換”“樣條變換”“單調(diào)變換”和“優(yōu)化計分變換”，變換的結(jié)果用“圖形”呈現(xiàn)出來，見圖1。

圖1 基于MODEL語句中的要求對4個變量做的變換結(jié)果

在圖1中，左上方一圖表明：定量因變量(nox)呈現(xiàn)出“對數(shù)”變化趨勢，提示可對其進行“對數(shù)變換”；右上角一圖表明：定量自變量(EqRatio)呈現(xiàn)出“二次拋物線”變化趨勢，提示可對其進行“二次多項式變換”；左下角一圖表明：多值有序自變量(CpRatio)呈現(xiàn)出近似“直線”變化趨勢，提示可對其進行“恒等變換(即不做變換)”；右下角一圖表明：多值名義自變量(Fuel)呈現(xiàn)出“兩種級別”，其中，水平代碼為“1、2、5”的燃料種類對結(jié)果的影響數(shù)量大約在“1與2”之間；而水平代碼為“0、3、4”的燃料種類對結(jié)果的影響數(shù)量大約在“4與5”之間，縱軸上的“數(shù)量”就是“優(yōu)化計分”的結(jié)果。

基于以上的“非參數(shù)多重非線性回歸分析結(jié)果(以圖1反映其變量變換的效果)”，可以構(gòu)建參數(shù)多重非參數(shù)回歸分析模型如下：

log(nox) = b0+b1*EqRatio+b2*EqRatio**2+b3*CpRatio+Sum b(j)*Fuel(j)+Error

此模型的含義是：對定量因變量(nox)進行自然對數(shù)變換，對定量自變量(EqRatio)進行二次多項式變換，對多值有序自變量(CpRatio)進行恒等變換，而對多值名義自變量(Fuel)進行優(yōu)化計分變換并將其6個計分值[Fuel(j)]與其回歸系數(shù)[b(j)]分別相乘后求和(j=0、1、2、3、4、5)(注意：相當(dāng)于把“6值名義自變量”視為6個新變量，就有6個回歸系數(shù)；但它們又屬于同一個多值名義變量，因此，其自由度為6-1=5，在本質(zhì)上相當(dāng)于是5個啞變量)。

3.4.4 構(gòu)建參數(shù)多重非線性回歸模型

采用SAS擬合上述參數(shù)多重非線性回歸模型所需要的SAS過程步程序如下：

title2' Now fit log(nox) = b0 + b1*EqRatio + b2*EqRatio**2 +';

title3' b3*CpRatio + Sum b(j)*Fuel(j) + Error';

*-Fit the Parametric Model Suggested by the Nonparametric Analysis-;

proc transreg data=sashelp.Gas solve ss2 short

nomiss plots=all cldetail;

model log(nox) = pspline(EqRatio / deg=2)

identity(CpRatio) opscore(Fuel);

run;

3.4.5 參數(shù)多重非線性回歸模型輸出結(jié)果

Univariate ANOVA Table Based on the Usual Degrees of FreedomSourceDFSum of SquaresMean SquareFPr > FModel879.338389.917298213.09<0.0001Error1607.446590.046541Corrected Total16886.78498

以上是總模型的假設(shè)檢驗結(jié)果，F(xiàn)=213.09，P<0.0001，表明總模型具有統(tǒng)計學(xué)意義。

Root MSE0.21573R-Square0.9142Dependent Mean0.63130Adj R-Sq0.9099Coeff Var34.17294

以上是擬合優(yōu)度評價結(jié)果：均方誤差平方根為0.21573，R2為0.9142，調(diào)整R2為0.9099，相對于表1中“模型1”的三個對應(yīng)的數(shù)值(0.30935、0.9586、0.9527)，均方誤差平方根略微變小了一點，然而，R2和調(diào)整R2的數(shù)值卻下降得比較多(因為模型中自變量的項數(shù)減少得很多，現(xiàn)在模型的自由度df=8；而模型1的df=21)。

Univariate Regression Table Based on the Usual Degrees of FreedomVariableDFCoefficientType II Sum of SquaresMean SquareFPr > FLabelIntercept1-15.27464957.133857.13381227.60<0.0001InterceptPspline.EqRatio_1135.10291462.747862.74781348.22<0.0001Equivalence Ratio 1Pspline.EqRatio_21-19.38646864.643064.64301388.94<0.0001Equivalence Ratio 2Identity(CpRatio)10.0320581.44451.444531.04<0.0001Compression RatioOpscore(Fuel)50.1583885.56191.112423.90<0.0001Fuel

以上為回歸模型的回歸系數(shù)與假設(shè)檢驗結(jié)果(因輸出結(jié)果過寬，將其中參數(shù)的置信區(qū)間等信息省略掉了)。其中，第2、3兩行上分別為“定量自變量(EqRatio)”的“一次項”與“二次項”的計算結(jié)果；第4行上為“多值有序自變量(CpRatio)”的計算結(jié)果；而第5行上為“多值名義自變量(Fuel)”的計算結(jié)果，具有5個自由度，但只有一個回歸系數(shù)。

基于上述回歸模型計算出定量因變量的預(yù)測值作為橫坐標(biāo)的數(shù)值，其原始觀測值作為縱坐標(biāo)的數(shù)值，繪制出散布圖，以直觀的方式呈現(xiàn)模型對資料的擬合效果。見圖2。

圖2 反映因變量(nox)的觀測值與預(yù)測值的散布圖

由圖2可知：此回歸模型對資料的擬合效果比較令人滿意。

4 總 結(jié)

4.1 本文方法的小結(jié)

本文對“多值名義自變量”進行“優(yōu)化計分變換(Opscore)”，取代了經(jīng)典統(tǒng)計學(xué)中常規(guī)的“啞變量變換”。在此基礎(chǔ)上，對“定量因變量(nox)”分別進行了“樣條變換(Spline)”“單調(diào)樣條變換(Mspline)”和“Box-Cox變換”三種變量變換；對“定量自變量(EqRatio)”分別進行了“樣條變換(Spline)”和“單調(diào)樣條變換(Mspline)”兩種變量變換；對“多值有序自變量(CpRatio)”分別進行了“單調(diào)變換(Monotone)”和“優(yōu)化計分變換(Opscore)”。從而，獲得“最優(yōu)模型”的擬合優(yōu)度評價指標(biāo)為“均方誤差平方根=0.30935、R2=0.9586、調(diào)整R2=0.9527”。由于此模型相對比較復(fù)雜(自變量有21項且寫不出具體的回歸系數(shù))；簡化后的參數(shù)多重非線性回歸模型僅含8個自變量(注：包括派生變量)，R2為0.9142，調(diào)整R2為0.9099。

4.2 其他方法的回顧

在本期科研方法專題的前三篇文章中，分別對“多值名義自變量”進行“啞變量變換”“算術(shù)均值變換”和“校正均值變換”，對“定量因變量(nox)”進行了“不變換”“自然對數(shù)變換”“平方根變換”“倒數(shù)變換”“指數(shù)變換”和“Logistic變換”；對“定量自變量(EqRatio)”和“多值有序自變量(CpRatio)(被視為‘定量的’)”引入了“派生變量”。還在模型中假定“包含截距項”和“不含截距項”兩種情形下，分別采取“前進法”“后退法”和“逐步法”篩選自變量。在前三篇文章中的每一篇中，都從大約72個模型中選出了24個“最優(yōu)模型”且以“模型24”為“最佳模型”。

4.3 四組方法的總結(jié)

“四組方法”指本期科研方法專題中四篇文章所介紹的方法，即分別對“多值名義自變量”進行“啞變量變換”“算術(shù)均值變換”“校正均值變換”和“優(yōu)化計分變換”。在此基礎(chǔ)上，再對其他變量采用多種變量變換方法，最后以“擬合優(yōu)度評價指標(biāo)”的取值為判斷標(biāo)準(zhǔn)，總結(jié)于下表2。

表2 四篇文章中“最佳”回歸模型擬合優(yōu)度評價指標(biāo)匯總

注：前三篇文章原先給出的是“均方誤差(MSE)”，分別為“0.00076”“0.00060”和“0.00061”，將它們分別開平方根，得到“0.02757”“0.02449”和“0.02470”。由表2可知，前三篇文章回歸模型的擬合效果非常接近，均優(yōu)于第4篇文章回歸模型的擬合效果；然而，第一篇文章回歸模型中保留了“部分啞變量”，嚴(yán)格來說，它不是非常理想的結(jié)果

由此可以得出如下結(jié)論：在對“定量因變量”構(gòu)建多重非線性回歸模型中，對“多值名義自變量”是采取“啞變量變換”“優(yōu)化計分變換”還是“均值或校正均值變換”，對結(jié)果的影響都不是特別大，關(guān)鍵在于以下四點：第一，應(yīng)對“定量因變量”“定量自變量”和“多值有序自變量”采取合適的變量變換；第二，應(yīng)盡可能引入較多的“派生變量”；第三，應(yīng)假定“模型中不含截距項”；第四，應(yīng)盡可能多采取一些篩選自變量的策略，如前進法、后退法和逐步法。

最后還需要注意：引入派生變量后獲得的“最佳”回歸模型通常具有“嚴(yán)重多重共線性”，若最終的回歸模型中的某些回歸系數(shù)的“正負號”與專業(yè)知識不符合，此時，再基于“最佳回歸模型”所確定的“自變量組合”進行“嶺回歸分析”。從而獲得滿意的結(jié)果[4]。