N∶M条件logistic回归分析在统计软件上的实现*

广州南方医科大学公共卫生与热带医学学院生物统计学系(510515) 郑一男 曹佩华 欧春泉

本文在STATA、R、SAS常用统计软件上实现N∶M条件logistic回归分析，并指出Cox风险比例模型进行条件logistic回归分析时结点调整方法的适用情况。

问题的提出

在有关疾病危险因素的流行病学调查研究中，常将病例与对照按年龄、性别等因素进行匹配，以消除某些混杂因素对研究结果的干扰作用〔1－3〕。大家熟悉的是1∶1或1∶M匹配，即给每个病例选择1个或M个对照与之匹配。匹配后往往再按匹配因素进行分层分析，这样可使每一个匹配层中都有一定数目的病例与对照。或者直接采用N∶M匹配，每一匹配层中有N个病例和M个对照，而且各层中病例数和对照数的比例可以不固定，研究实施更为灵活。配对的病例-对照研究数据均采用条件logistic回归分析。关于1∶1和1∶M条件logistic回归分析的软件应用有不少讨论，但对N∶M条件logistic回归分析却鲜有介绍，笔者在实际工作中遇到不少非统计专业的研究人员错误地运用SPSS软件的Cox模块进行N∶M条件logistic回归分析，导致结果出现较大偏差。本文将全面介绍STATA、R和SAS统计软件如何对N∶M配对病例-对照研究资料进行条件logistic回归分析，并详细阐述在使用Cox模块进行条件logistic回归分析时的注意事项。读者可根据自身的使用偏好选择任一软件进行分析。

原理及方法

*:广东省科技计划项目(项目号:0903109)△通讯作者:欧春泉，E-mail:ocq@fimmu．com

2．Cox比例风险模型

运用Cox比例风险模型可实现条件logistic回归分析〔5〕。由于同一层内的生存时间均相同，故生存时间的结点(ties)问题非常突出，如何处理结点是分析的一个关键。目前主流统计软件中有四种处理方法:Breslow、Efron、Exact marginal likelihood(精确边际似然法)与 Exact partial likelihood(精确偏似然法)。Breslow是最简单快捷的，Efron比Breslow更加精确，但计算速度略慢，两种方法在处理大量结点数据时均会使参数估计结果出现严重偏差。精确边际似然法实际上是利用15-point Gauss-Laguerre quadrature近似计算得到，对于结点数较少的数据，该法比Efron慢;随着结点数增加，其效率优势能迅速显露出来，但当平均结点数大于30时〔6〕，误差将逐渐加大。对于生存时间为离散变量的资料，由于产生的结点较多，应使用精确偏似然法进行处理。

在1∶1和1∶M资料中，每层中仅有一个结点，因此四种结点处理方法均等价可行。对于N∶M资料(N为病例数)，每层中有N个结点，使用精确偏似然法最为合适。

下面给出Breslow与精确偏似然法的对数似然函数 ln L〔6〕:

在时间区间(t0i，ti］中的第 i个观测(i=1，…，N)

有:

j为死亡时间 t(j)的顺序编号(j=1，…，D);dj，rj分别为t(j)时刻的死亡数与存活数;Dj，Rj分别为在t(j)时刻死亡病人的集合与存活病人的集合(危险集);k为观测编号使得t0k＜t(j)≤tk;δij为0-1变量，表示病人死亡(=1)或截尾(=0)。在病例-对照研究中，上述“死亡”与“存活”的概念分别对应“病例”与“对照”。比照似然函数(2)和(4)可知，使用精确偏似然法调整结点的Cox比例风险模型与条件logistic回归模型完全等价。

实例分析

1．资料概述

为研究婴儿出现低出生体重(体重≤2500g)的危险因素，Hosmer和Lemeshow按照母亲年龄采取配对的病例-对照研究〔7〕，调查共涉及189名女性，其中59名产有低出生体重婴儿，130名所产的婴儿体重正常。数据文件可在南方医科大学生物统计学系网站(http://www．echobelt．org)下载。

表1 变量说明

2．条件logistic模型的软件实现

STATA 9．0和R语言均提供了条件logistic分析模块，即clogit模块。

(1)STATA 9．0对实例资料的分析过程为:

clogit low lwt smoke ht ut，group(age)/* 输出参数估计值及置信域*/clogit，or/*输出OR值及置信域*/

(2)R2．11．1对实例资料的分析过程为:

library(splines)

library(survival)#载入自带软件包镜像

clogit(Low～LWT+Smoke+HT+UT+strata(Age))

STATA和R软件的logit模块分析结果完全相同，整理结果于表2。表中最后两列展示的是比数比(Odds Ratio，OR)的95%置信区间。

表2 clogit模块结果整理

3．Cox比例风险模型的软件实现

各软件处理结点的方法大体相同，但需要注意区分各软件命名与含义的差异:STATA中有breslow(默认)、efron、exactm(精确边际似然法)与exactp(精确偏似然法);R中“method”选项里包括 efron(默认)、breslow与exact(精确偏似然法);SAS中“ties”选项里包括breslow(默认)、efron、exact(精确边际似然法)与discrete(精确偏似然法);SPSS没有结点处理的可选项，运算过程使用内置的breslow法。

使用Cox比例风险模型前，需在原始数据基础上增加生存时间变量Time=2－Low(注:对照组生存时间大于病例组即可)。

(1)STATA 9．0的Cox比例风险模型为stcox模块，实例资料的分析过程为:

stset time，failure(low)/*定义时间变量与反应变量*/

stcox lwt smoke ht ut，strata(age)exactp/* 输出HR值及置信域*/

stcox，nohr/*输出参数估计值及置信域*/

(2)R2．11．1的 Cox比例风险模型为 coxph模块，实例资料的分析过程为:

coxph(Surv(Time，Low)～ LWT+Smoke+HT+UT+strata(Age)，method="exact")

(3)SAS9．2对于1∶1资料的条件logistic回归可以使用logistic过程(proc logistic)来完成，但对于1∶M和N∶M资料，就必须使用Cox比例风险模型(proc phreg)，实例资料的分析过程为:

data LBW;

input Age Low LWT Smoke HT UT;

*［在此输入数据］

Time=2－Low;

run;

proc phreg data=LBW;

model Time*Low(0)=LWT Smoke HT UT/ties=discrete;

strata Age;

run;

上述三种方法的结果汇总于表3，可见与表2中条件logistic回归分析结果完全一致，但输出结果通常称为风险比(Hazard Ratio，HR)，而非比数比。表中展示的是风险比的95%置信区间。

表3 各软件Cox模块结果整理

4．SPSS17．0 软件的误用

在SPSS17．0软件的菜单操作中没有提供处理条件logistic回归的模块，如果使用SPSS软件的Cox模块对实例资料进行分析，结果如表4:

表4 SPSS软件Cox Regression模块结果整理

可以看出上述估计值低估了真实结果，原因就在于结点的处理使用的是Breslow法。因此SPSS软件目前在菜单操作中尚无法实现N∶M资料的条件logistic回归分析。但对于1∶1与1∶M的资料是可行的。

讨 论

STATA软件与R软件中均有条件logistic回归的专用模块(clogit)，非常简便灵活;SAS软件中没有类似的专用模块，而是使用与之等价的Cox风险比例回归模块(phreg)进行分析。使用前需要在原数据基础上加入生存时间变量(time)，并且使用精确偏似然法处理结点，这在STATA中的stcox模块与R软件中coxph模块也得到了印证。1∶1和1∶M配对设计是N∶M配对的特例，故本文介绍的N∶M条件logistic回归分析的方法同样适合于1∶1或1∶M条件logistic回归。

1．方亚，胡海兰．女性乳腺癌危险因素及其变化．中国卫生统计，2009，26(3):241-246．

2．刘静，王洁贞，薛付忠，等．肾综合征出血热发病率与气象因素关系的研究．中国卫生统计，2006，23(4):326-329．

3．陆云霞，施侣元，余红平，等．武汉市居民饮食模式与食管癌发病关系的条件logistic回归分析．中国卫生统计，2005，22(3):146-148．

4．STATA Base Reference Manual Volume 1 Release 10．Texas，USA:STATA，2007，286．

5．张业武．Cox比例风险模型对条件logistic回归参数估计原理和方法．中国卫生统计，2002，19(1):23-25．

6．STATA Survival Analysis and Epidemiological Table Reference Manual Release 10．Texas，USA:STATA，2007，153．

7．Hosmer D，Lemeshow S．Applied Logistic Regression．2nd edit．New York，USA:John Wiley ＆ Sons，1989，319．