北京大學生育健康研究所/衛(wèi)生部生育健康重點實驗室(100191) 李宏田 袁 悅 李智文 王琳琳 劉建蒙

GAM識別非線性相關及其在醫(yī)學統(tǒng)計建模中的應用*

北京大學生育健康研究所/衛(wèi)生部生育健康重點實驗室(100191) 李宏田 袁 悅 李智文 王琳琳 劉建蒙△

目的 介紹廣義相加模型(GAM)識別非線性相關及其在醫(yī)學統(tǒng)計建模中的應用。方法 應用SAS軟件PROC GAM模塊識別實例數據結局變量與自變量之間的非線性相關，通過比較考慮該非線性相關和不考慮該非線性相關時多元線性回歸和logistic回歸模型的擬合和預測效果，闡明GAM識別非線性相關在統(tǒng)計建模中的重要性。結果 與不考慮非線性相關的模型相比，考慮非線性相關的模型擬合和預測效果更優(yōu)。結論 合理使用GAM，在模型中納入非線性成分，可改善回歸模型的建模效果和預測精度。

廣義相加模型 非線性相關 統(tǒng)計建模

*:國家自然科學基金面上項目(編號:81072372)和科技部973項目(編號:2007CB5119001)資助

△通訊作者:劉建蒙，E-mail:liujm@pku．edu．cn

廣義相加模型GAM，于1986年由Hastie和Tibshirani提出〔1，2〕。GAM 是對傳統(tǒng)廣義線性模型(包括多元線性回歸和logistic回歸模型)的擴展。廣義線性模型一般形式為E(Y|X1，X2，…，Xp)= β0+ β1X1+β2X2+… +βpXp，而 GAM 一般形式是E(Y|X1，X2，…，Xp)=β0+f1(X1)+f2(X2)+… +fp(Xp)。fp(Xp)是關于Xp的非指定類別的非參數函數，其估計方法有平滑樣條法(smoothing splines)、局部加權回歸散點平滑法(LOESS)和薄盤平滑樣條法(thin-plate smoothing spline);平滑參數選擇的方法有交叉驗證(cross validation)或廣義交叉驗證(generalized cross validation)。SAS軟件設有專門的GAM模塊，是GAM建模常用軟件之一〔3〕。本文將采用SAS軟件PROC GAM模塊識別實例數據結局變量與自變量之間的非線性相關，通過比較考慮該非線性相關和不考慮該非線性相關時多元線性回歸和logistic回歸模型的擬合和預測效果，闡明GAM識別非線性相關在統(tǒng)計建模中的重要性。

實例數據

實例數據是關于兒童智商(IQ)影響因素的研究資料，于2000年收集，樣本量為7340;變量及其分布見表1。建模時文化程度(EDU)按啞變量進入模型;以初中及以下組為參照，大專及以上、高中或中專和不詳組對應的啞變量依次為EDU1、EDU2和EDU3。

兒童智商(CIQ) 99．0±16．3 兒童高智商(CIQTOP) 9．4%)兒童月齡(CAGE) 67．7±7．4 母親文化程度(EDU)母親智商(MIQ) 94．5±17．0 大專及以上 4．3%母親年齡(MAGE)25．5±3．1 高中或中專 17．6%初中及以下 76．7%不詳 1．4%

模型擬合和評價方法

7340名兒童隨機分成數據集IQSAMPLE(n=3687)和IQTEST(n=3653)。IQSAMPLE用于建模，IQTEST用于預測評價。先以CIQTOP為因變量，以CAGE，MAGE，MIQ和EDU為自變量建立logistic回歸模型1;進而通過GAM識別CIQTOP與CAGE，MAGE和MIQ之間是否有非線性相關以及非線性相關的具體類型，并將其引入logistic回歸模型，建立模型2。利用赤池信息準則(AIC)比較模型1和2的建模效果，AIC值越小，建模效果越好。比較模型1和2用于IQSAMPLE預測時KAPPA統(tǒng)計量的最大取值，KAPPA值越大，建模效果越好;將IQTEST分別回代至模型1和2，以前述KAPPA值最大時的判別概率為標準，比較兩個模型用于IQTEST預測的KAPPA值。再以 CIQ為因變量，以 CAGE，MAGE，MIQ和 EDU為自變量建立多元線性回歸模型1和2，建模過程與logistic回歸建模類似，也利用AIC比較模型1和2的建模效果。將IQTEST分別回代至模型1和2，比較兩個模型殘差平方和的大小，殘差平方和越小，建模效果越好。

SAS程序與建模評價

1．logistic回歸建模

PROC LOGISTIC DATA=IQSAMPLE DESC;MODEL CIQTOP=EDU1 EDU2 EDU3 MIQ CAGE MAGE;RUN;

logistic回歸模型 1的 AIC值為 1984.02，對IQSAMPLE預測的最大KAPPA值為0．305，相應的判別概率為0．23;據此概率值，模型1用于IQTEST預測的KAPPA值為0．307。參數估計結果見表2。

(2)GAM識別非線性相關

①程序及主要結果

PROC GAM DATA=IQSAMPLE;MODEL CIQTOP=PARAM(EDU1 EDU2 EDU3)SPLINE(MIQ)SPLINE(MAGE)SPLINE(CAGE)/DIST=BINOMINAL;QUIT;

表2 logistic回歸模型1參數估計

調用 GAM程序，擬合 IQSAMPLE數據集，以CIQTOP為因變量，EDU以啞變量形式按參數函數〔PARAM(變量名)〕進行擬合，MIQ、MAGE和CAGE按非參數函數〔SPLINE(變量名)〕進行擬合。DIST指定CIQTOP呈二項分布(BINOMINAL)。GAM程序擬合非參數函數默認自由度為4，線性部分為1，非線性部分為3。SAS主要輸出結果見表3－5，第1部分與模型1參數估計基本一致，僅CAGE檢驗的P值由0．11變?yōu)?．07。第3部分MAGE非線性部分檢驗有統(tǒng)計學意義，即MAGE與IQTOP呈非線性相關。

表3 GAM參數函數及非參函數線性部分估計結果

表4 GAM非參數函數非線性部分平滑擬合結果

表5 GAM非參數函數非線性部分假設檢驗結果

②程序及主要結果

2017年互聯(lián)網期刊出版行業(yè)的主要出版商仍然是以同方知網（北京）技術有限公司（以下簡稱同方知網）、萬方數據科技有限公司（以下簡稱萬方數據）、重慶維普資訊有限公司（以下簡稱維普資訊）、龍源數字傳媒集團（以下簡稱龍源數媒）四家出版企業(yè)占市場最大份額，還有其他出版企業(yè)也開始接觸互聯(lián)網期刊業(yè)務。

PROC GAM DATA=IQSAMPLE;MODEL CIQTOP=PARAM(EDU1 EDU2 EDU3)SPLINE(MIQ)SPLINE(MAGE)SPLINE(CAGE)/METHOD=GCV DIST=BINOMINAL;QUIT;

程序①擬合非參數函數默認自由度為4;程序②增加了METHOD=GCV語句，指定參數估計方法為廣義交叉驗證法，不限定自由度。參數函數及非參數函數線性部分的擬合結果與程序1基本一致，非參數函數非線性部分的假設檢驗仍顯示MAGE與CIQTOP呈非線性相關，CAGE非線性部檢驗的P值減小至0.053，提示CAGE與CIQTOP呈非線性相關(表6)。

表6 廣義交叉驗證法GAM假設檢驗結果

③程序及主要結果

ODS HTML;ODS GRAPHICS ON;PROC GAM DATA=IQSAMPLE PLOT(CLM);MODEL CIQTOP=PARAM(EDU1 EDU2 EDU3)SPLINE(MIQ)SPLINE(MAGE)SPLINE(CAGE)/METHOD=GCV LINK=LOGIT DIST=BINOMINAL;QUIT;ODS GRAPHICS OFF;ODS HTML CLOSE;

程序②顯示MAGE和CAGE與CIQTOP均呈非線性相關，程序③增加了ODS GRAPHICS和PLOT(CLM)語句，該語句會輸出非參數函數非線性部分對CIQTOP影響的效應圖，見圖1。MAGE和CAGE非參數函數非線性部分對CIQTOP影響近似于二次方曲線，MAGE曲線開口向下，CAGE曲線開口向上?；诖饲€，預期在Logistic回歸模型中增加MAGE和CAGE的二次方項會改善建模和預測效果。

圖1 GAM SAS程序③的部分輸出結果

(3)logistic回歸模型2的SAS程序及主要結果PROC LOGISTIC DATA=IQSAMPLE DESC;MODEL CIQTOP=EDU1 EDU2 EDU3 MIQ CAGE CAGE*CAGE MAGE MAGE*MAGE;RUN;

logistic回歸模型 2的AIC值為 1970.66，對IQSAMPLE預測的最大KAPPA值為0.324，相應的判別概率為0.22;據此概率值，模型2用于IQTEST預測的KAPPA值為0.349。參數估計結果見表7。模型2的AIC值小于模型1，對IQSAMPLE和IQTEST預測的最大KAPPA值均大于模型1，表明模型2優(yōu)于模型1。MAGE二次方項檢驗有統(tǒng)計學意義，CAGE二次方項檢驗的P值為0.06，接近有統(tǒng)計學意義;回歸系數符號所反映的開口方向與GAM輸出的MAGE和CAGE非線性部分效應圖相吻合。

表7 logistic回歸模型2參數估計

多元線性回歸建模

1．多元線性回歸模型1的SAS程序及主要結果

PROC REG DATA=IQSAMPLE;MODEL CIQ=EDU1 EDU2 EDU3 MIQ CAGE MAGE;QUIT;

多元線性回歸模型1的AIC值為19614.33，對IQSAMPLE預測的殘差平方和為 750603.39，對IQTEST預測的殘差平方和為728649．26;參數估計結果見表8。

表8 多元線性回歸模型1參數估計

2．GAM 非線性相關識別

(1)程序及主要結果

PROC GAM DATA=IQSAMPLE;MODEL CIQ=PARAM(EDU1 EDU2 EDU3)SPLINE(MIQ)SPLINE(MAGE)SPLINE(CAGE)/DIST=GAUSSIAN;QUIT;

調用GAM程序，擬合IQSAMPLE數據集，以CIQ為因變量，EDU以啞變量形式按參數函數進行擬合，按默認自由度(df=4)對MIQ、MAGE和CAGE進行非參數函數擬合。DIST指定CIQ的分布為高斯分布(GAUSSIAN)，默認的連接函數為IDENTITY。結果見表9－11。表9與多元線性回歸模型1參數估計基本一致。表11顯示MAGE非線性部分檢驗有統(tǒng)計學意義，CAGE檢驗P值為0.06。

表9 GAM參化函數及非參數函數線性部分估計結果

表10 GAM非參數函數非線性部分平滑擬合結果

表11 GAM非參數函數非線性部分假設檢驗結果

(2)程序及主要結果

PROC GAM DATA=IQSAMPLE;MODEL CIQ=PARAM(EDU1 EDU2 EDU3)SPLINE(MIQ)SPLINE(MAGE)SPLINE(CAGE)/METHOD=GCV DIST=GAUSSIAN;QUIT;

程序增加GCV語句，參數函數及非參數函數線性部分的擬合結果與程序(1)基本一致，非參數函數非線性部分的假設檢驗仍顯示MAGE與IQ呈非線性相關，MIQ和CAGE的自由度遠小于0，檢驗P值無法估計，提示MIQ和CAGE與IQ基本無非線性相關(表12)。

表12 廣義交叉驗證法假設檢驗結果

(3)程序及主要結果

ODS HTML;ODS GRAPHICS ON;PROC GAM DATA=IQSAMPLE;MODEL CIQ=PARAM(EDU1 EDU2 EDU3)SPLINE(MIQ)SPLINE(MAGE)SPLINE(CAGE)/METHOD=GCV IST=GAUSSIAN;QUIT;ODSGRAPHICSOFF;ODSHTML CLOSE;

程序輸出了MAGE、CAGE和MIQ非參數函數非線性部分對IQ影響效應曲線(圖2)。盡管CAGE和MIQ圖像近似二次方曲線，但其效應值(縱坐標)遠小于MAGE，自由度遠小于0，檢驗P值無法估計，提示此類曲線無實際意義。MAGE曲線較為復雜，但母親分娩年齡在22～35歲之間時，近似呈二次方曲線，而這部分人群占總人群的比例達91．5%，提示兩側曲線的穩(wěn)定性弱。

圖2 GAM程序(3)部分輸出結果

(4)程序及主要結果

ODS HTML;ODS GRAPHICS ON;PROC GAM DATA=IQSAMPLE;MODEL CIQ=PARAM(EDU1 EDU2 EDU3 MIQ CAGE)SPLINE(MAGE，DF=3)/DIST=GAUSSIAN;QUIT;ODS GRAPHICS OFF;ODS HTML CLOSE;

基于以上輸出結果，對CAGE和MIQ按參數函數擬合，并限定 MAGE的總自由度為2和3，以簡化MAGE非線性部分的效應曲線。簡化后的圖像均呈二次方曲線(圖3)。

3．多元回歸模型2的SAS程序及主要結果

PROC REG DATA=IQSAMPLE;MODEL CIQ=EDU1 EDU2 EDU3 MIQ CAGE MAGE MAGE_SQUARE;QUIT;

MAGE_SQUARE是新生成的變量，是MAGE的平方項。多元線性回歸模型2的AIC值為19597.41，對IQSAMPLE預測的殘差平方和為746762.65，對IQTEST預測的殘差平方和為725704.04;參數估計結果見表13。模型2的AIC值以及對IQSAMPLE和IQTEST預測的殘差平方和均小于模型1，提示模型2優(yōu)于模型1。MAGE二次方項檢驗有統(tǒng)計學意義，回歸系數符號所反映的開口方向與GAM輸出的MAGE非線性部分效應曲線相吻合。

圖3 GAM程序(4)的部分輸出結果

表13 多元線性回歸模型2參數估計

討 論

簡要介紹了GAM有關知識及其SAS程序，通過實例數據說明了GAM識別變量間非線性相關對統(tǒng)計建模的重要性。強調了如何使用GAM識別變量間非線性相關，并將識別出的非線性相關引入經典的多元線性回歸和logistic回歸模型，進而對比評價了引入非線性成分和未引入線性成分的模型。GAM用于識別變量間非線性相關的特點是直觀性好，以統(tǒng)計學檢驗為基礎，可同時考察因變量與諸多自變量間的關系;合理使用GAM可改善多元線性回歸和logistic回歸模型的建模效果和預測精度。

1．Hastie T，Tibshirani R．Generalized additive models．Stat Sci，1986，1(3):297-318．

2．Hastie TJ，Tibshirani RJ．Generalized additive models．New York，NY:Chapman and Hall，Inc，1990．

3．SAS Institute Inc．SAS/STAT User's Guide，Version 9．2．Cary，NC:SAS Institute Inc，2008．

An Introduction of GAM in Identifying Non-linear Correlations and its Application in Statistical Modeling

Li Hongtian，Yuan Yue，Li Zhiwen，et al．Institute of Reproductive and Child Health/Ministry ofHealth Key Laboratory ofReproductive Health，Peking University Health Science Center(100191)，Beijing

ObjectiveTo introduce Generalized Additive Models(GAM)in identifying non-linear correlations and its application in statistical modeling for medical research data.MethodsA dataset was used for modeling with SAS PROC GAM．Goodness of fit and prediction precision were compared between models with and without non-linear components.ResultsA non-linear correlation could be identified by GAM．Compared with models without non-linear components，goodness of fit and prediction precision were improved by involving non-linear components.ConclusionModels with non-linear components reflect a true relationship between dependent and independent variables and hence improve the predictive ability．

Generalized additive models;Non-linear correlations;Statistical modeling