郭 威 于菲菲 秦嬰逸 何 倩 吳 騁

第二軍醫(yī)大學衛(wèi)生統(tǒng)計學教研室(200433)

·綜述·

觀察性研究中校正連續(xù)型混雜因素的非線性方法*

郭 威 于菲菲 秦嬰逸 何 倩 吳 騁△

第二軍醫(yī)大學衛(wèi)生統(tǒng)計學教研室(200433)

混雜偏倚(confounding bias)是觀察性研究中的一類重要偏倚，它是指由于混雜因素既與暴露因素又與結局存在相關關系，導致暴露與結局之間的真實關系受到了干擾而產(chǎn)生的偏倚[1]。因此，觀察性研究中如何控制混雜一直是研究人員所關注的重要問題。在統(tǒng)計分析階段一種常用的處理辦法是將混雜因素納入回歸模型中進行校正。實際問題中常常遇到這樣的情況，即混雜變量為連續(xù)型指標，該變量與結局變量間的真實的內(nèi)在關系有可能是非線性的，比如BMI與死亡率的關系可能為U型，即低于或超出BMI的參考值范圍，個人的死亡風險都會增加，BMI處于參考值范圍內(nèi)時個人的死亡風險最低。在過去幾十年中，主流的控制連續(xù)型混雜因素的做法有兩種，一是直接假定該變量與結局變量的關系為線性的，將其直接以線性形式納入回歸模型中加以校正，二是將連續(xù)型混雜變量進行分類(categorization)，以虛擬變量的形式納入回歸模型中。這兩種處理方法雖然操作簡單，但存在不少問題。近十年來，統(tǒng)計軟件的進一步開發(fā)和豐富使得許多復雜的非線性建模策略實現(xiàn)起來不再困難，越來越多的醫(yī)學科研人員開始將這些方法應用于自己的研究，就目前而言，劑量-反應分析和各種預后模型的建立就是非線性建模策略的重要應用方向[2-3]。在回歸模型中，若采用非線性方法對連續(xù)型混雜變量進行擬合，就能夠更充分地校正混雜從而獲得漸進無偏的暴露-結局關系。本文將首先對控制連續(xù)型混雜的兩種傳統(tǒng)回歸方法及其局限性進行探討，然后就流行的三種校正連續(xù)型混雜的非線性方法和相應的軟件實現(xiàn)作簡要介紹。

傳統(tǒng)線性回歸方法及其局限性

傳統(tǒng)的多元回歸模型中校正混雜因素的方法主要有兩種：一種是以變量的線性形式直接引入模型，另一種是將連續(xù)型變量轉(zhuǎn)化為多分類變量。在回歸模型中直接引入線性項的前提假定是該變量與結局變量之間為線性關系，如果這個線性假定與真實的關系偏差不大，則該方法能較好地控制混雜因素。然而，該簡化手段反映出研究人員對于事物內(nèi)部特定的真實關系存在一定程度的認知上的缺乏，因此，當線性假定嚴重背離了真實的關系時，這種回歸建模策略就損失了較多的信息，從而大大降低了統(tǒng)計效能。

將連續(xù)型變量轉(zhuǎn)化為二分類或多分類資料，然后納入回歸模型中進行校正，被廣泛地應用于臨床實踐和觀察性流行病學研究中。常見的做法是通過特定的閾值將人群分為具有某種屬性和不具有該屬性的兩個類別，比如吸煙和不吸煙人群，或者基于某個連續(xù)型變量的分布將數(shù)據(jù)劃分為多個類別，比如以分位數(shù)為臨界值將收入水平劃分為高、較高、中、較低、低五個類別等。分類化處理本身有許多優(yōu)勢：當采集資料時將某些連續(xù)型指標分段處理具有特定的臨床意義；本身的擬合優(yōu)度較高；當采用logistic或Cox回歸時結果可表示為OR或HR，容易解釋和相互比較。然而，分類處理也存在一定的局限性：當校正混雜因素時，某個連續(xù)型變量被分成若干類別，同一個類別內(nèi)部所有的數(shù)據(jù)點被視為具有同質(zhì)性，不可避免地損失了層內(nèi)變異信息，會產(chǎn)生殘余混雜(residual confounding)，有研究證實兩分類時殘余混雜最嚴重，此時的效應估計與未校正前相比(粗效應值)沒有太大變化[4]，并且增大了I類錯誤[5]，因此流行病學研究中較少采用兩分類法來校正混雜因素，而多采用四分類或五分類，如果繼續(xù)增加混雜變量的分類數(shù)，統(tǒng)計效能的增加不再明顯，反而大大增加了模型的復雜度。

校正連續(xù)型混雜因素的非線性方法

許多模擬及實證研究表明，回歸模型中的非線性方法在校正連續(xù)型混雜時能夠很好地降低殘余混雜[4,6-7]。尤其在某些情況下暴露對結局的影響較弱而混雜對結局的影響較強，例如，為研究空氣污染與肺癌發(fā)病率的關系，需要控制吸煙的影響，而吸煙對空氣污染而言是強混雜，必須加以很好地控制才能獲得漸近無偏地估計，此時采用非線性建模策略顯然具有明顯的優(yōu)勢。近年來，文獻中被經(jīng)常應用的方法主要有回歸樣條、廣義可加模型和分數(shù)多項式等，現(xiàn)依次介紹如下。

1.回歸樣條(regression spline)

(1)模型簡介 樣條函數(shù)是一種用于曲線擬合的特殊的分段多項式。最簡單的樣條函數(shù)是線性樣條，它將自變量分為多個區(qū)間，每個區(qū)間內(nèi)的數(shù)據(jù)用線性函數(shù)擬合，不同區(qū)間線段的斜率不同，在相鄰兩個區(qū)間的結合部，即節(jié)點處將區(qū)間內(nèi)的兩條直線相連接。線性樣條的優(yōu)點是形式簡潔，模型中的參數(shù)容易解釋，缺點是曲線在各個節(jié)點處的形狀往往不光滑，所以線性樣條模型在實際應用中并不常見。立方樣條是指在各個區(qū)間內(nèi)部的樣條函數(shù)表達式均為三次多項式Si(X)，并且在相鄰兩個區(qū)間的節(jié)點處二階連續(xù)可導，使曲線在節(jié)點處的形狀變得光滑。不過立方樣條也有其局限性，表現(xiàn)在曲線的兩個尾部，即第一個節(jié)點之前和最后一個節(jié)點之后，擬合效果不佳。限制性立方樣條(restricted cubic spline，RCS)[8]，也稱自然樣條，是在立方樣條的基礎上附加限制條件，令兩尾部區(qū)間內(nèi)的函數(shù)為線性形式。若令tk表示X在第k個節(jié)點處的值，則RCS如下式：

RCS(X)=β0+β1S1(X)+β2S2(X)+…+βk-1Sk-1(X)

上式中，S1(X)=X,j=1,2,…,k-2,

與立方樣條相比，RCS不僅優(yōu)化了數(shù)據(jù)兩尾部的估計，還使模型中的待估參數(shù)減少，對于用k個節(jié)點的RCS表示的變量，模型只需要估計k-1個參數(shù)。應用回歸樣條之前，需要根據(jù)變量間的內(nèi)在關系預先設定節(jié)點的數(shù)量k和位置t。RCS中的節(jié)點的數(shù)量對模型擬合影響較大，實際應用中一般取3～7個節(jié)點，在標準的統(tǒng)計軟件中，節(jié)點位置的設定默認為變量的等距百分位數(shù)。表1列出了選取特定節(jié)點數(shù)量時節(jié)點位置對應的百分位數(shù)[8]。另外，還可以采用AIC(Akaike information criteria)來選擇最佳的節(jié)點數(shù)量，對不同節(jié)點數(shù)建立不同的樣條回歸模型，得到不同的AIC值，選擇AIC值最小時的模型作為最優(yōu)的模型。

將樣條函數(shù)引入廣義線性模型如多元線性回歸、logistic回歸、Cox回歸等，即得到回歸樣條模型。Brenner認為當進行探索性分析時，需要精確地擬合劑量-反應關系，因此最好采用5個節(jié)點的樣條函數(shù)；當研究目的是校正連續(xù)型混雜時，為了模型的簡潔性，采用3個節(jié)點的樣條函數(shù)就足夠了[6]。國內(nèi)一些學者也對RCS進行了應用或評價，羅劍鋒等利用實例數(shù)據(jù)對比了logistic回歸中RCS和多分類處理的分析效果，認為RCS對數(shù)據(jù)的擬合效果更好[9]；余紅梅等將Cox回歸和RCS結合探索了生存數(shù)據(jù)中急性白血病患者持續(xù)緩解時間和預后之間的劑量-反應關系[10]。

(2)軟件實現(xiàn)Desquilbet等在2009年編寫了回歸樣條RCS的SAS Macro，該程序的特點有：①可以對連續(xù)型變量建立RCS函數(shù)；②可以對一個連續(xù)型暴露和一個結局的曲線及95%置信區(qū)間帶圖形展示；③適用范圍包括線性模型、logistic模型、Cox模型和廣義估計方程；④可以給出全模型及其中的非線性關系的統(tǒng)計學檢驗值。Ruifeng Li等在2010年開發(fā)了SAS LGTPHCURV9 Macro[11]，該程序可以在非條件logistic回歸、條件logistic回歸、pooled logistic回歸和比例風險模型中擬合RCS，在控制混雜變量的同時研究暴露與結局(OR或IRR)的關系；根據(jù)用戶指定的節(jié)點數(shù)來自動選擇節(jié)點位置等。

2.廣義可加模型(generalized additive model，GAM)

(1)模型簡介 廣義可加模型是廣義線性模型(GLM)的擴展，由Hastie和Tibshirani于1990年提出。GAM保留了GLM中反應變量的分布和連接函數(shù)的多樣性的特性，不同的是，它的預測變量采用非參數(shù)形式。它不需要預先對模型進行線性假定，唯一需要的假定是各函數(shù)項是可加且光滑的，克服了維度的影響，通過“加性”假設，GAM能將一些與因變量存在復雜非線性關系的自變量以不同函數(shù)加和的形式進入模型，從而可以探索到變量間非單調(diào)和非線性的關系，具有較高的靈活性。GAM的模型表達式為：

其中fj(xj),(j=1,…,p)為自變量xj的光滑函數(shù)。從GAM的形式上可以看出，它對自變量的形式?jīng)]有規(guī)定，具有較好的靈活性；連接函數(shù)可根據(jù)資料的分布類型的不同而不同，比如資料為正態(tài)分布時，連接函數(shù)為probit；資料為二項分布時，連接函數(shù)為logit等。光滑函數(shù)Sj的擬合方法有多種，常用的有核光滑函數(shù)法、局部加權散點圖平滑法(LOESS或LOWESS)和光滑樣條(smoothing spline)等。光滑樣條GAM的估計方法通常為懲罰最小二乘法。所謂懲罰最小二乘法就是在最小二乘法的基礎上增加了一個懲罰項來保證樣條函數(shù)擬合的預測變量在節(jié)點處的光滑性，如下式：

其中，前一項為最小二乘項，后一項為懲罰項，λ是光滑參數(shù)，使上式最小就可得到fj。光滑參數(shù)的設定不僅要使觀測值和估計值之間的距離之和達到最小，即達到較好的曲線擬合優(yōu)度，還要控制回歸曲線的光滑度，所以理想的λ是曲線擬合優(yōu)度和光滑度的一種折中。實際操作中模型的估計采用局部記分(local-scoring)算法，該算法是由迭代再加權最小二乘法與backfitting過程合并而成。光滑函數(shù)fj的選擇通常根據(jù)廣義交叉驗證的偏差(generalized cross-validated deviance)和AIC等。

國內(nèi)很多學者對廣義可加模型進行了研究和應用。陳長生等較早地對光滑樣條非參數(shù)回歸進行了部分理論探索和實例應用[12]。馮國雙等用實例說明了使用SAS 8.2中GAM模塊擬合廣義可加模型的過程[13]。類似于多元線性回歸中的共線性問題，GAM也可能存在共曲線性(concurvity)問題，當有共曲線性存在時，它會低估模型參數(shù)項的標準誤，增大I類錯誤和導致模型的解不唯一。近年來發(fā)展的非參數(shù)條件自助抽樣法是克服共曲線性影響的方法之一[14]。另外，在實際擬合數(shù)據(jù)時，模型的估計結果可能受到離群點的影響而產(chǎn)生偏差，王彤等將穩(wěn)健估計的思想方法引入到GAM中，通過模擬對Y方向的存在離群點的情況進行了討論，導出穩(wěn)健估計較一般估計的結果更加可靠[15]。由于GAM靈活性強，并且可以有效控制與時間相關的混雜因素的影響，目前國內(nèi)外多將其應用于探索環(huán)境污染物和人體疾病之間關系的環(huán)境流行病學領域[16-17]，此類數(shù)據(jù)多為時間序列數(shù)據(jù)，殘差的自相關性可能增大I類錯誤，針對這個問題余松林等提出在GAM基礎上通過增加反應變量函數(shù)的勻滑函數(shù)的方法，有效地校正了時間序列中殘差的自相關性對參數(shù)假設檢驗的影響[18]。

(2)軟件實現(xiàn) SAS軟件設有專門的GAM模塊，作為SAS軟件中非參數(shù)回歸建模的重要過程，PROC GAM具有優(yōu)良的多維適應性和結果的可解釋性。在SAS 9.3中，它的主要特點有：①支持一元光滑樣條、二元薄板光滑樣條和局部回歸平滑法；②能夠擬合非參數(shù)可加及半?yún)?shù)可加模型；③支持多個SCORE語句；④允許用戶自定義光滑參數(shù)或根據(jù)GCV自動選擇光滑參數(shù)；⑤可以通過ODS圖形系統(tǒng)進行圖形展示等[19]。

3.分數(shù)多項式(fractional polynomials，F(xiàn)P)

(1)模型簡介 FP由Royston和Altman兩位統(tǒng)計學家提出[20]，它是二次和立方多項式的擴展，與傳統(tǒng)多項式不同之處在于，F(xiàn)P的冪可以是整數(shù)，也可以是分數(shù)，故稱為分數(shù)多項式，又譯作分式多項式。FP模型的形式如下：

針對實際中常需處理多個暴露或混雜變量的問題，Sauerbrei和Royston提出了多元分數(shù)多項式(multivariable fractional polynomials，MFP)[2]。為確定需要進入模型的變量，研究人員往往依靠專業(yè)知識和文獻報道，當研究新事物時，可以獲得的背景知識很少，這就為篩選進入模型的變量造成了很大的困難。MFP的一大優(yōu)勢是它能同時篩選重要的暴露或混雜變量和確定FP模型的函數(shù)形式，它將向后剔除法和FP自適應算法結合起來構建多元模型，是一種基于數(shù)據(jù)的建模方法[22]。Royston等提出當候選變量個數(shù)很多時，對于連續(xù)型暴露變量和連續(xù)型混雜變量應分別設置剔除標準，暴露變量的標準宜嚴，推薦0.01或0.05，混雜變量的標準宜寬，推薦0.10或0.20。

FP模型的思想是采用單個函數(shù)形式來擬合某個變量的樣本數(shù)據(jù)，這一點使得FP對數(shù)據(jù)的局部特征不敏感，與樣條函數(shù)相比，如果FP擬合暴露變量，其靈活性稍差。不過，如果待擬合的變量為混雜變量，那么它產(chǎn)生的殘余混雜很小，校正混雜的效果不會受到太大影響。

(2)軟件實現(xiàn) STATA、SAS、R軟件中逐步加入了FP和MFP的程序或命令，見表2。關于程序的更多細節(jié)描述可參考相關文獻[23]。

結 語

觀察性研究中混雜因素可在回歸模型中加以校正。對于特定情況下的某些連續(xù)型的混雜變量，傳統(tǒng)回歸方法——線性擬合和分類處理并不能完全控制其混雜效應，導致殘余混雜的產(chǎn)生。本文介紹的三種方法——回歸樣條、廣義可加模型、分數(shù)多項式是目前較為流行的校正連續(xù)型混雜的非線性建模方法。盡管它們的回歸系數(shù)不便于解釋，不過混雜變量與結局變量的具體的數(shù)量關系不是主要的關注點，實際上人們對暴露結局的關聯(lián)關系更感興趣，此時非線性方法扮演的角色避免了這個局限性。鑒于目前研究人員對連續(xù)型混雜因素的重視程度還不高，本文對三種控制連續(xù)型混雜因素的非線性方法的基本理論和軟件實現(xiàn)作了簡要的介紹，供廣大醫(yī)學科研工作者參考。

[1]胡永華,耿直.關于混雜概念的討論.中華流行病學雜志,2001,22(6):459-461.

[2]Sauerbrei W,Royston P.Building multivariable prognostic and diagnostic models:transformation of the predictors by using fractional polynomials.Journal of the Royal Statistical Society:Series A(Statistics in Society),1999,162(1):71-94.

[3]Desquilbet L,Mariotti F.Dose-response analyses using restricted cubic spline functions in public health research.Statistics in medicine,2010,29(9):1037-1057.

[4]Groenwold RH,Klungel OH,Altman DG,et al.Adjustment for continuous confounders:an example of how to prevent residual confounding.Canadian Medical Association Journal,2013,185(5):401.

[5]Austin PC,Brunner LJ.Inflation of the type I error rate when a continuous confounding variable is categorized in logistic regression analyses.Statistics in medicine,2004,23(7):1159-1178.

[6]Brenner H,Blettner M.Controlling for continuous confounders in epidemiologic research.Epidemiology,1997,8(4):429-434.

[7]Benedetti A,Abrahamowicz M.Using generalized additive models to reduce residual confounding.Statistics in medicine,2004,23(24):3781-3801.

[8]Harrell FE.Regression modeling strategies:with applications to linear models,logistic regression,and survival analysis.Springer,2001,20-23.

[9]羅劍鋒,金歡,李寶月.限制性立方樣條在非線性回歸中的應用研究.中國衛(wèi)生統(tǒng)計,2010,27(3):229-232.

[10]余紅梅,羅艷虹,吳燕萍.基于三次樣條Cox回歸的劑量-反應關系分析.中國衛(wèi)生統(tǒng)計,2012,29(5):721-722.

[11]Li R,Hertzmark E,Louie M,et al.The SAS LGTPHCURV9 Macro.Boston,MA Channing Laboratory,2011.

[12]陳生長,徐勇勇.光滑樣條非參數(shù)回歸方法及醫(yī)學應用.中國衛(wèi)生統(tǒng)計,1999,16(6):342-345.

[13]馮國雙,陳景武,張國英.用GAM程序擬合光滑樣條非參數(shù)回歸.數(shù)理醫(yī)藥學雜志,2005,18(5):403-405.

[14]賈彬,王彤,王琳娜.廣義可加模型共曲線性及其在空氣污染問題研究中的應用.第四軍醫(yī)大學學報,2005,26(3):280-283.

[15]王彤,賈彬,王琳娜.廣義可加模型穩(wěn)健估計在空氣污染對健康影響評價中的應用.中國衛(wèi)生統(tǒng)計,2007,24(3):245-247.

[16]Dominici F,Mcdermott A,Zeger SL,et al.On the use of generalized additive models in time-series studies of air pollution and health.American journal of epidemiology,2002,156(3):193-203.

[17]莫運政,鄭亞安,陶輝.日均氣溫與呼吸系統(tǒng)疾病急診人次相關性的時間序列分析.北京大學學報(醫(yī)學版),2012,44(3):416-420.

[18]余松林,彭曉武.廣義加性模型配合時間序列資料時消除殘差自相關性的一種方法.中國衛(wèi)生統(tǒng)計,2010,27(5):450-454.

[19]Cai W.Fitting generalized additive models with the GAM procedure in SAS 9.2.SAS Institute Inc,Cary NC(USA),2008.

[20]Royston P,Altman DG.Regression using fractional polynomials of continuous covariates:parsimonious parametric modelling.Applied Statistics,1994,429-467.

[21]Royston P,Sauerbrei W.Multivariable model-building:a pragmatic approach to regression anaylsis based on fractional polynomials for modelling continuous variables.John Wiley & Sons,2008.

[22]Sauerbrei W,Royston P,Binder H.Selection of important variables and determination of functional form for continuous predictors in multivariable model building.Stat Med,2007,26(30):5512-5528.

[23]Sauerbrei W,Meier-Hirmer C,Benner A,et al.Multivariable regression model building by using fractional polynomials:Description of SAS,STATA and R programs.Computational Statistics & Data Analysis,2006,50(12):3464-3485.

(責任編輯：郭海強)

*：上海市公共衛(wèi)生重點學科建設項目(12GWZX0602)；上海市軟科學研究重點項目(14692101700)；總后優(yōu)秀青年科技人才扶持對象項目；第二軍醫(yī)大學衛(wèi)勤系基金項目(2014WK02)資助

△通信作者：吳騁，Email：wucheng_wu@126.com