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

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

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

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

廣義相加模型 非線性相關(guān) 統(tǒng)計(jì)建模

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

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

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

實(shí)例數(shù)據(jù)

實(shí)例數(shù)據(jù)是關(guān)于兒童智商(IQ)影響因素的研究資料，于2000年收集，樣本量為7340;變量及其分布見表1。建模時文化程度(EDU)按啞變量進(jìn)入模型;以初中及以下組為參照，大專及以上、高中或中專和不詳組對應(yīng)的啞變量依次為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名兒童隨機(jī)分成數(shù)據(jù)集IQSAMPLE(n=3687)和IQTEST(n=3653)。IQSAMPLE用于建模，IQTEST用于預(yù)測評價。先以CIQTOP為因變量，以CAGE，MAGE，MIQ和EDU為自變量建立logistic回歸模型1;進(jìn)而通過GAM識別CIQTOP與CAGE，MAGE和MIQ之間是否有非線性相關(guān)以及非線性相關(guān)的具體類型，并將其引入logistic回歸模型，建立模型2。利用赤池信息準(zhǔn)則(AIC)比較模型1和2的建模效果，AIC值越小，建模效果越好。比較模型1和2用于IQSAMPLE預(yù)測時KAPPA統(tǒng)計(jì)量的最大取值，KAPPA值越大，建模效果越好;將IQTEST分別回代至模型1和2，以前述KAPPA值最大時的判別概率為標(biāo)準(zhǔn)，比較兩個模型用于IQTEST預(yù)測的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預(yù)測的最大KAPPA值為0．305，相應(yīng)的判別概率為0．23;據(jù)此概率值，模型1用于IQTEST預(yù)測的KAPPA值為0．307。參數(shù)估計(jì)結(jié)果見表2。

(2)GAM識別非線性相關(guān)

①程序及主要結(jié)果

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

表2 logistic回歸模型1參數(shù)估計(jì)

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

表3 GAM參數(shù)函數(shù)及非參函數(shù)線性部分估計(jì)結(jié)果

表4 GAM非參數(shù)函數(shù)非線性部分平滑擬合結(jié)果

表5 GAM非參數(shù)函數(shù)非線性部分假設(shè)檢驗(yàn)結(jié)果

②程序及主要結(jié)果

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

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

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

表6 廣義交叉驗(yàn)證法GAM假設(shè)檢驗(yàn)結(jié)果

③程序及主要結(jié)果

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均呈非線性相關(guān)，程序③增加了ODS GRAPHICS和PLOT(CLM)語句，該語句會輸出非參數(shù)函數(shù)非線性部分對CIQTOP影響的效應(yīng)圖，見圖1。MAGE和CAGE非參數(shù)函數(shù)非線性部分對CIQTOP影響近似于二次方曲線，MAGE曲線開口向下，CAGE曲線開口向上。基于此曲線，預(yù)期在Logistic回歸模型中增加MAGE和CAGE的二次方項(xiàng)會改善建模和預(yù)測效果。

圖1 GAM SAS程序③的部分輸出結(jié)果

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

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

表7 logistic回歸模型2參數(shù)估計(jì)

多元線性回歸建模

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

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

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

表8 多元線性回歸模型1參數(shù)估計(jì)

2．GAM 非線性相關(guān)識別

(1)程序及主要結(jié)果

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

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

表9 GAM參化函數(shù)及非參數(shù)函數(shù)線性部分估計(jì)結(jié)果

表10 GAM非參數(shù)函數(shù)非線性部分平滑擬合結(jié)果

表11 GAM非參數(shù)函數(shù)非線性部分假設(shè)檢驗(yàn)結(jié)果

(2)程序及主要結(jié)果

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

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

表12 廣義交叉驗(yàn)證法假設(shè)檢驗(yàn)結(jié)果

(3)程序及主要結(jié)果

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非參數(shù)函數(shù)非線性部分對IQ影響效應(yīng)曲線(圖2)。盡管CAGE和MIQ圖像近似二次方曲線，但其效應(yīng)值(縱坐標(biāo))遠(yuǎn)小于MAGE，自由度遠(yuǎn)小于0，檢驗(yàn)P值無法估計(jì)，提示此類曲線無實(shí)際意義。MAGE曲線較為復(fù)雜，但母親分娩年齡在22～35歲之間時，近似呈二次方曲線，而這部分人群占總?cè)巳旱谋壤_(dá)91．5%，提示兩側(cè)曲線的穩(wěn)定性弱。

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

(4)程序及主要結(jié)果

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;

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

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

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

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

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

表13 多元線性回歸模型2參數(shù)估計(jì)

討 論

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

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