混杂是生物医学研究中最为棘手的问题之一,混杂可以在设计阶段,采用配对方法将可能的混杂因素加以控制,以提高研究效率和可靠性。然而这样的资料确不能采用经典的logistic回归,需要使用条件logistic回归。
理论部分
为何采用条件logistic回归?
这里以1:1配对设计为例,1:1配对设计的特点是对子内部控制的混杂因素一致,不同对子的这些因素不同。从统计学角度来看,把每个配对试做一个成组病例对照研究,将配对层设置为哑变量。按照传统的logistic回归方法建立模型,估计每个自变量的比值比,这样做有的主要困难有:
其一、每个配对层仅有2个观察样品;
其二、哑变量树是对子数-1,大大增加了估计参数的数量。
为克服原有方法解决配对资料参数估计的缺陷,往往通过构造特殊的条件似然函数,仍然采用极大似然估计法进行参数估计。
构造条件似然函数
1:1病例对照的资料常常整理成如下格式。设有n对独立的观察,每个对子包含两个人,第1个已经患病,第2个没有患病。自变量为X,第i层第一个人自变量记为Xi1Xi1,第二个人自变量记为Xi2Xi2。
配对号病例
对照
XYXY
1X11X111X10X100
2X21X211X20X200
…………………………
nXn1Xn11Xn0Xn00
任何一层中,第一个人患病的概率和未患病的概率分别为:
π1=exp(β0+XT1β)1+exp(β0+XT1β)1−π1=11+exp(β0+XT1β)π1=exp(β0+X1Tβ)1+exp(β0+X1Tβ)1−π1=11+exp(β0+X1Tβ)
第二个人患病的概率与不患病的概率分别为:
π0=exp(β0+XT0β)1+exp(β0+XT0β)1−π0=11+exp(β0+XT0β)π0=exp(β0+X0Tβ)1+exp(β0+X0Tβ)1−π0=11+exp(β0+X0Tβ)
现在假定同一层中的2个人中,只有1个人患病。在只有1个人患病的情况下,恰好第1个人患病而第2个人不患病的条件概率为
P=π1(1−π0)π1(1−π0)+π0(1−π1)=11+exp(−(XT1−XT0)β)P=π1(1−π0)π1(1−π0)+π0(1−π1)=11+exp(−(X1T−X0T)β)
n个配对层的联合概率,即似然函数为
L=P1×P2×…Pn=∏11+exp(−(XT1−XT0)β)L=P1×P2×…Pn=∏11+exp(−(X1T−X0T)β)
对上面的式子求对数,得到对数似然方程,通过迭代法求得偏回归系数的值。假设检验同样采用似然比、score检验以及wald检验。
需要注意,上面的模型与经典的logistic回归有两点不同:
与偏回归系数相乘的是病例与对照相应变量之差;
模型不含常数项。
实战部分
在多种统计软件如SPSS、SAS等中,配对logistic回归通常采用分层COX回归模型来实现。这是因为,分层Cox回归假设不同层之间的基线风险函数完全无关,只需估计协变量的偏回归系数值。条件logistic回归同样如此,模型中不存在截距项,不同层之间有相同的偏回归系数。可见用分层Cox回归来拟合条件logistic回归完全是一种取巧。
在调用分层Cox回归时,会将病例算作发生终点事件,对照算作删失。下面以一个例子来说明:
为探讨女性乳腺癌危险因素,研究者在某市1996-1997年间确诊的女性乳腺癌患者中随机抽取350名病例,对每一病例以一名性别相同、年两差别不超过±2.5岁的对照。收集的信息包括:
变量变量含义值标签
X1X1文化程度0:大专以下,1:大专及以上
X2X2体质指数0:小于等于27,1:大于27
X3X3近年精神压抑0:无,1:有
X4X4乳腺良性疾病史0:无,1:有
X5X5恶性肿瘤家族史0:无,1:有
X6X6初潮年龄0:大于等于14岁,1:小于14岁
X7X7哺乳史0:有,1:无
部分数据如下,数据中以id号来标识层变量,同一个id中第一个为病例,后一个为对照。:
idx1x2x3x4x5x6x7
10000010
10000001
20100101
20000001
30000001
30000001
40010001
40000011
51111010
现在调用R程序进行条件logistic回归,使用的函数是survival包中的clogit函数,该函数实际上默认调用了coxph函数。下面的程序首先读取了数据,然后创建结局变量,进而将分类变量转化为因子。
1
2
3
4
5
6
7
8
9
10
11
df<-read.delim('clipboard',header = T)
df$result<-rep(c(1,0),times=350)
head(df)
df$x1<-factor(df$x1,labels = c('大专以下','大专及以上'),levels=0:1)
df$x2<-factor(df$x2,levels = 0:1,labels = c('≤27','>27'))
df$x3<-factor(df$x3,levels = 0:1,labels = c('无','有'))
df$x4<-factor(df$x4,levels = 0:1,labels = c('无','有'))
df$x5<-factor(df$x5,levels = 0:1,labels = c('无','有'))
df$x6<-factor(df$x6,levels = 0:1,labels = c('≤14岁','>14岁'))
df$x7<-factor(df$x7,levels = 0:1,labels = c('有','无'))
head(df)
clogit具有和其他建模函数同样的用法,只需要使用strata参数指定分层变量即可
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
> library(survival)
> clogit.full<-clogit(result~.-id+strata(id),data = df)
> summary(clogit.full)
Call:
coxph(formula = Surv(rep(1, 700L), result) ~ . - id + strata(id),
data = df, method = "exact")
n= 700, number of events= 350
coef exp(coef) se(coef) z Pr(>|z|)
x1大专及以上 0.61655 1.85252 0.24368 2.530 0.011403 *
x2>27 1.07581 2.93237 0.47854 2.248 0.024570 *
x3有 1.38708 4.00313 0.26438 5.247 1.55e-07 ***
x4有 1.90482 6.71821 0.53637 3.551 0.000383 ***
x5有 0.60911 1.83879 0.26039 2.339 0.019326 *
x6>14岁 0.13795 1.14792 0.17979 0.767 0.442893
x7无 -0.08851 0.91530 0.23665 -0.374 0.708407
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
exp(coef) exp(-coef) lower .95 upper .95
x1大专及以上 1.8525 0.5398 1.1490 2.987
x2>27 2.9324 0.3410 1.1478 7.491
x3有 4.0031 0.2498 2.3843 6.721
x4有 6.7182 0.1488 2.3480 19.223
x5有 1.8388 0.5438 1.1038 3.063
x6>14岁 1.1479 0.8711 0.8070 1.633
x7无 0.9153 1.0925 0.5756 1.455
Rsquare= 0.102 (max possible= 0.5 )
Likelihood ratio test= 75.58 on 7 df, p=1.091e-13
Wald test = 50.46 on 7 df, p=1.173e-08
Score (logrank) test = 65.08 on 7 df, p=1.447e-11
变量x6和x7不是显著的,现在将其剔除,
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
> clogit.mod<-clogit(result~.-id+strata(id)-x6-x7,data = df)
> summary(clogit.mod)
Call:
coxph(formula = Surv(rep(1, 700L), result) ~ . - id + strata(id) -
x6 - x7, data = df, method = "exact")
n= 700, number of events= 350
coef exp(coef) se(coef) z Pr(>|z|)
x1大专及以上 0.6105 1.8414 0.2433 2.509 0.012092 *
x2>27 1.0975 2.9965 0.4775 2.299 0.021529 *
x3有 1.3690 3.9313 0.2631 5.204 1.95e-07 ***
x4有 1.9001 6.6864 0.5327 3.567 0.000361 ***
x5有 0.6069 1.8348 0.2588 2.345 0.019039 *
---
Signif. codes: 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
exp(coef) exp(-coef) lower .95 upper .95
x1大专及以上 1.841 0.5431 1.143 2.966
x2>27 2.997 0.3337 1.175 7.639
x3有 3.931 0.2544 2.348 6.583
x4有 6.686 0.1496 2.354 18.993
x5有 1.835 0.5450 1.105 3.047
Rsquare= 0.101 (max possible= 0.5 )
Likelihood ratio test= 74.83 on 5 df, p=1.01e-14
Wald test = 49.86 on 5 df, p=1.483e-09
Score (logrank) test = 64.61 on 5 df, p=1.35e-12
剩下的变量都是显著的,输出的结果最后一部分是该模型与空模型的似然比检验,wald检验以及score检验,结果同样有统计学意义。
上述模型中,X1、X2、X3、X4、X5均有统计学意义,OR值均大于1。因此可以认为高文化程度、肥胖、精神压抑、乳腺良性疾病史和恶性肿瘤家族史是女性乳腺癌的危险因素。