观察性研究中校正连续型混杂因素的非线性方法*

郭 威 于菲菲 秦婴逸 何 倩 吴 骋

第二军医大学卫生统计学教研室(200433)

·综述·

观察性研究中校正连续型混杂因素的非线性方法*

郭 威 于菲菲 秦婴逸 何 倩 吴 骋△

第二军医大学卫生统计学教研室(200433)

混杂偏倚(confounding bias)是观察性研究中的一类重要偏倚，它是指由于混杂因素既与暴露因素又与结局存在相关关系，导致暴露与结局之间的真实关系受到了干扰而产生的偏倚[1]。因此，观察性研究中如何控制混杂一直是研究人员所关注的重要问题。在统计分析阶段一种常用的处理办法是将混杂因素纳入回归模型中进行校正。实际问题中常常遇到这样的情况，即混杂变量为连续型指标，该变量与结局变量间的真实的内在关系有可能是非线性的，比如BMI与死亡率的关系可能为U型，即低于或超出BMI的参考值范围，个人的死亡风险都会增加，BMI处于参考值范围内时个人的死亡风险最低。在过去几十年中，主流的控制连续型混杂因素的做法有两种，一是直接假定该变量与结局变量的关系为线性的，将其直接以线性形式纳入回归模型中加以校正，二是将连续型混杂变量进行分类(categorization)，以虚拟变量的形式纳入回归模型中。这两种处理方法虽然操作简单，但存在不少问题。近十年来，统计软件的进一步开发和丰富使得许多复杂的非线性建模策略实现起来不再困难，越来越多的医学科研人员开始将这些方法应用于自己的研究，就目前而言，剂量-反应分析和各种预后模型的建立就是非线性建模策略的重要应用方向[2-3]。在回归模型中，若采用非线性方法对连续型混杂变量进行拟合，就能够更充分地校正混杂从而获得渐进无偏的暴露-结局关系。本文将首先对控制连续型混杂的两种传统回归方法及其局限性进行探讨，然后就流行的三种校正连续型混杂的非线性方法和相应的软件实现作简要介绍。

传统线性回归方法及其局限性

传统的多元回归模型中校正混杂因素的方法主要有两种：一种是以变量的线性形式直接引入模型，另一种是将连续型变量转化为多分类变量。在回归模型中直接引入线性项的前提假定是该变量与结局变量之间为线性关系，如果这个线性假定与真实的关系偏差不大，则该方法能较好地控制混杂因素。然而，该简化手段反映出研究人员对于事物内部特定的真实关系存在一定程度的认知上的缺乏，因此，当线性假定严重背离了真实的关系时，这种回归建模策略就损失了较多的信息，从而大大降低了统计效能。

将连续型变量转化为二分类或多分类资料，然后纳入回归模型中进行校正，被广泛地应用于临床实践和观察性流行病学研究中。常见的做法是通过特定的阈值将人群分为具有某种属性和不具有该属性的两个类别，比如吸烟和不吸烟人群，或者基于某个连续型变量的分布将数据划分为多个类别，比如以分位数为临界值将收入水平划分为高、较高、中、较低、低五个类别等。分类化处理本身有许多优势：当采集资料时将某些连续型指标分段处理具有特定的临床意义；本身的拟合优度较高；当采用logistic或Cox回归时结果可表示为OR或HR，容易解释和相互比较。然而，分类处理也存在一定的局限性：当校正混杂因素时，某个连续型变量被分成若干类别，同一个类别内部所有的数据点被视为具有同质性，不可避免地损失了层内变异信息，会产生残余混杂(residual confounding)，有研究证实两分类时残余混杂最严重，此时的效应估计与未校正前相比(粗效应值)没有太大变化[4]，并且增大了I类错误[5]，因此流行病学研究中较少采用两分类法来校正混杂因素，而多采用四分类或五分类，如果继续增加混杂变量的分类数，统计效能的增加不再明显，反而大大增加了模型的复杂度。

校正连续型混杂因素的非线性方法

许多模拟及实证研究表明，回归模型中的非线性方法在校正连续型混杂时能够很好地降低残余混杂[4,6-7]。尤其在某些情况下暴露对结局的影响较弱而混杂对结局的影响较强，例如，为研究空气污染与肺癌发病率的关系，需要控制吸烟的影响，而吸烟对空气污染而言是强混杂，必须加以很好地控制才能获得渐近无偏地估计，此时采用非线性建模策略显然具有明显的优势。近年来，文献中被经常应用的方法主要有回归样条、广义可加模型和分数多项式等，现依次介绍如下。

1.回归样条(regression spline)

(1)模型简介 样条函数是一种用于曲线拟合的特殊的分段多项式。最简单的样条函数是线性样条，它将自变量分为多个区间，每个区间内的数据用线性函数拟合，不同区间线段的斜率不同，在相邻两个区间的结合部，即节点处将区间内的两条直线相连接。线性样条的优点是形式简洁，模型中的参数容易解释，缺点是曲线在各个节点处的形状往往不光滑，所以线性样条模型在实际应用中并不常见。立方样条是指在各个区间内部的样条函数表达式均为三次多项式Si(X)，并且在相邻两个区间的节点处二阶连续可导，使曲线在节点处的形状变得光滑。不过立方样条也有其局限性，表现在曲线的两个尾部，即第一个节点之前和最后一个节点之后，拟合效果不佳。限制性立方样条(restricted cubic spline，RCS)[8]，也称自然样条，是在立方样条的基础上附加限制条件，令两尾部区间内的函数为线性形式。若令tk表示X在第k个节点处的值，则RCS如下式：

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

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

与立方样条相比，RCS不仅优化了数据两尾部的估计，还使模型中的待估参数减少，对于用k个节点的RCS表示的变量，模型只需要估计k-1个参数。应用回归样条之前，需要根据变量间的内在关系预先设定节点的数量k和位置t。RCS中的节点的数量对模型拟合影响较大，实际应用中一般取3～7个节点，在标准的统计软件中，节点位置的设定默认为变量的等距百分位数。表1列出了选取特定节点数量时节点位置对应的百分位数[8]。另外，还可以采用AIC(Akaike information criteria)来选择最佳的节点数量，对不同节点数建立不同的样条回归模型，得到不同的AIC值，选择AIC值最小时的模型作为最优的模型。

将样条函数引入广义线性模型如多元线性回归、logistic回归、Cox回归等，即得到回归样条模型。Brenner认为当进行探索性分析时，需要精确地拟合剂量-反应关系，因此最好采用5个节点的样条函数；当研究目的是校正连续型混杂时，为了模型的简洁性，采用3个节点的样条函数就足够了[6]。国内一些学者也对RCS进行了应用或评价，罗剑锋等利用实例数据对比了logistic回归中RCS和多分类处理的分析效果，认为RCS对数据的拟合效果更好[9]；余红梅等将Cox回归和RCS结合探索了生存数据中急性白血病患者持续缓解时间和预后之间的剂量-反应关系[10]。

(2)软件实现Desquilbet等在2009年编写了回归样条RCS的SAS Macro，该程序的特点有：①可以对连续型变量建立RCS函数；②可以对一个连续型暴露和一个结局的曲线及95%置信区间带图形展示；③适用范围包括线性模型、logistic模型、Cox模型和广义估计方程；④可以给出全模型及其中的非线性关系的统计学检验值。Ruifeng Li等在2010年开发了SAS LGTPHCURV9 Macro[11]，该程序可以在非条件logistic回归、条件logistic回归、pooled logistic回归和比例风险模型中拟合RCS，在控制混杂变量的同时研究暴露与结局(OR或IRR)的关系；根据用户指定的节点数来自动选择节点位置等。

2.广义可加模型(generalized additive model，GAM)

(1)模型简介 广义可加模型是广义线性模型(GLM)的扩展，由Hastie和Tibshirani于1990年提出。GAM保留了GLM中反应变量的分布和连接函数的多样性的特性，不同的是，它的预测变量采用非参数形式。它不需要预先对模型进行线性假定，唯一需要的假定是各函数项是可加且光滑的，克服了维度的影响，通过“加性”假设，GAM能将一些与因变量存在复杂非线性关系的自变量以不同函数加和的形式进入模型，从而可以探索到变量间非单调和非线性的关系，具有较高的灵活性。GAM的模型表达式为：

其中fj(xj),(j=1,…,p)为自变量xj的光滑函数。从GAM的形式上可以看出，它对自变量的形式没有规定，具有较好的灵活性；连接函数可根据资料的分布类型的不同而不同，比如资料为正态分布时，连接函数为probit；资料为二项分布时，连接函数为logit等。光滑函数Sj的拟合方法有多种，常用的有核光滑函数法、局部加权散点图平滑法(LOESS或LOWESS)和光滑样条(smoothing spline)等。光滑样条GAM的估计方法通常为惩罚最小二乘法。所谓惩罚最小二乘法就是在最小二乘法的基础上增加了一个惩罚项来保证样条函数拟合的预测变量在节点处的光滑性，如下式：

其中，前一项为最小二乘项，后一项为惩罚项，λ是光滑参数，使上式最小就可得到fj。光滑参数的设定不仅要使观测值和估计值之间的距离之和达到最小，即达到较好的曲线拟合优度，还要控制回归曲线的光滑度，所以理想的λ是曲线拟合优度和光滑度的一种折中。实际操作中模型的估计采用局部记分(local-scoring)算法，该算法是由迭代再加权最小二乘法与backfitting过程合并而成。光滑函数fj的选择通常根据广义交叉验证的偏差(generalized cross-validated deviance)和AIC等。

国内很多学者对广义可加模型进行了研究和应用。陈长生等较早地对光滑样条非参数回归进行了部分理论探索和实例应用[12]。冯国双等用实例说明了使用SAS 8.2中GAM模块拟合广义可加模型的过程[13]。类似于多元线性回归中的共线性问题，GAM也可能存在共曲线性(concurvity)问题，当有共曲线性存在时，它会低估模型参数项的标准误，增大I类错误和导致模型的解不唯一。近年来发展的非参数条件自助抽样法是克服共曲线性影响的方法之一[14]。另外，在实际拟合数据时，模型的估计结果可能受到离群点的影响而产生偏差，王彤等将稳健估计的思想方法引入到GAM中，通过模拟对Y方向的存在离群点的情况进行了讨论，导出稳健估计较一般估计的结果更加可靠[15]。由于GAM灵活性强，并且可以有效控制与时间相关的混杂因素的影响，目前国内外多将其应用于探索环境污染物和人体疾病之间关系的环境流行病学领域[16-17]，此类数据多为时间序列数据，残差的自相关性可能增大I类错误，针对这个问题余松林等提出在GAM基础上通过增加反应变量函数的匀滑函数的方法，有效地校正了时间序列中残差的自相关性对参数假设检验的影响[18]。

(2)软件实现 SAS软件设有专门的GAM模块，作为SAS软件中非参数回归建模的重要过程，PROC GAM具有优良的多维适应性和结果的可解释性。在SAS 9.3中，它的主要特点有：①支持一元光滑样条、二元薄板光滑样条和局部回归平滑法；②能够拟合非参数可加及半参数可加模型；③支持多个SCORE语句；④允许用户自定义光滑参数或根据GCV自动选择光滑参数；⑤可以通过ODS图形系统进行图形展示等[19]。

3.分数多项式(fractional polynomials，FP)

(1)模型简介 FP由Royston和Altman两位统计学家提出[20]，它是二次和立方多项式的扩展，与传统多项式不同之处在于，FP的幂可以是整数，也可以是分数，故称为分数多项式，又译作分式多项式。FP模型的形式如下：

针对实际中常需处理多个暴露或混杂变量的问题，Sauerbrei和Royston提出了多元分数多项式(multivariable fractional polynomials，MFP)[2]。为确定需要进入模型的变量，研究人员往往依靠专业知识和文献报道，当研究新事物时，可以获得的背景知识很少，这就为筛选进入模型的变量造成了很大的困难。MFP的一大优势是它能同时筛选重要的暴露或混杂变量和确定FP模型的函数形式，它将向后剔除法和FP自适应算法结合起来构建多元模型，是一种基于数据的建模方法[22]。Royston等提出当候选变量个数很多时，对于连续型暴露变量和连续型混杂变量应分别设置剔除标准，暴露变量的标准宜严，推荐0.01或0.05，混杂变量的标准宜宽，推荐0.10或0.20。

FP模型的思想是采用单个函数形式来拟合某个变量的样本数据，这一点使得FP对数据的局部特征不敏感，与样条函数相比，如果FP拟合暴露变量，其灵活性稍差。不过，如果待拟合的变量为混杂变量，那么它产生的残余混杂很小，校正混杂的效果不会受到太大影响。

(2)软件实现 STATA、SAS、R软件中逐步加入了FP和MFP的程序或命令，见表2。关于程序的更多细节描述可参考相关文献[23]。

结 语

观察性研究中混杂因素可在回归模型中加以校正。对于特定情况下的某些连续型的混杂变量，传统回归方法——线性拟合和分类处理并不能完全控制其混杂效应，导致残余混杂的产生。本文介绍的三种方法——回归样条、广义可加模型、分数多项式是目前较为流行的校正连续型混杂的非线性建模方法。尽管它们的回归系数不便于解释，不过混杂变量与结局变量的具体的数量关系不是主要的关注点，实际上人们对暴露结局的关联关系更感兴趣，此时非线性方法扮演的角色避免了这个局限性。鉴于目前研究人员对连续型混杂因素的重视程度还不高，本文对三种控制连续型混杂因素的非线性方法的基本理论和软件实现作了简要的介绍，供广大医学科研工作者参考。

[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]罗剑锋,金欢,李宝月.限制性立方样条在非线性回归中的应用研究.中国卫生统计,2010,27(3):229-232.

[10]余红梅,罗艳虹,吴燕萍.基于三次样条Cox回归的剂量-反应关系分析.中国卫生统计,2012,29(5):721-722.

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

[12]陈生长,徐勇勇.光滑样条非参数回归方法及医学应用.中国卫生统计,1999,16(6):342-345.

[13]冯国双,陈景武,张国英.用GAM程序拟合光滑样条非参数回归.数理医药学杂志,2005,18(5):403-405.

[14]贾彬,王彤,王琳娜.广义可加模型共曲线性及其在空气污染问题研究中的应用.第四军医大学学报,2005,26(3):280-283.

[15]王彤,贾彬,王琳娜.广义可加模型稳健估计在空气污染对健康影响评价中的应用.中国卫生统计,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]莫运政,郑亚安,陶辉.日均气温与呼吸系统疾病急诊人次相关性的时间序列分析.北京大学学报(医学版),2012,44(3):416-420.

[18]余松林,彭晓武.广义加性模型配合时间序列资料时消除残差自相关性的一种方法.中国卫生统计,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.

(责任编辑：郭海强)

*：上海市公共卫生重点学科建设项目(12GWZX0602)；上海市软科学研究重点项目(14692101700)；总后优秀青年科技人才扶持对象项目；第二军医大学卫勤系基金项目(2014WK02)资助

△通信作者：吴骋，Email：wucheng_wu@126.com