|
混杂是生物医学研究中最为棘手的问题之一，混杂可以在设计阶段，采用配对方法将可能的混杂因素加以控制，以提高研究效率和可靠性。然而这样的资料确不能采用经典的logistic回归，需要使用条件logistic回归。
理论部分为何采用条件logistic回归？这里以1:1配对设计为例，1:1配对设计的特点是对子内部控制的混杂因素一致，不同对子的这些因素不同。从统计学角度来看，把每个配对试做一个成组病例对照研究，将配对层设置为哑变量。按照传统的logistic回归方法建立模型，估计每个自变量的比值比，这样做有的主要困难有： 其一、每个配对层仅有2个观察样品； 其二、哑变量树是对子数-1，大大增加了估计参数的数量。
为克服原有方法解决配对资料参数估计的缺陷，往往通过构造特殊的条件似然函数，仍然采用极大似然估计法进行参数估计。 构造条件似然函数1:1病例对照的资料常常整理成如下格式。设有n对独立的观察，每个对子包含两个人，第1个已经患病，第2个没有患病。自变量为X，第i层第一个人自变量记为Xi1Xi1，第二个人自变量记为Xi2Xi2。 | 配对号 | 病例 |
| 对照 |
|
|---|
| X | Y | X | Y | | 1 | X11X11 | 1 | X10X10 | 0 | | 2 | X21X21 | 1 | X20X20 | 0 | | …… | …… | …… | …… | …… | | n | Xn1Xn1 | 1 | Xn0Xn0 | 0 |
任何一层中，第一个人患病的概率和未患病的概率分别为： π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中第一个为病例，后一个为对照。： | id | x1 | x2 | x3 | x4 | x5 | x6 | x7 |
|---|
| 1 | 0 | 0 | 0 | 0 | 0 | 1 | 0 | | 1 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | | 2 | 0 | 1 | 0 | 0 | 1 | 0 | 1 | | 2 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | | 3 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | | 3 | 0 | 0 | 0 | 0 | 0 | 0 | 1 | | 4 | 0 | 0 | 1 | 0 | 0 | 0 | 1 | | 4 | 0 | 0 | 0 | 0 | 0 | 1 | 1 | | 5 | 1 | 1 | 1 | 1 | 0 | 1 | 0 |
现在调用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。因此可以认为高文化程度、肥胖、精神压抑、乳腺良性疾病史和恶性肿瘤家族史是女性乳腺癌的危险因素。
|
