|
混雜是生物醫(yī)學(xué)研究中最為棘手的問(wèn)題之一，混雜可以在設(shè)計(jì)階段，采用配對(duì)方法將可能的混雜因素加以控制，以提高研究效率和可靠性。然而這樣的資料確不能采用經(jīng)典的logistic回歸，需要使用條件logistic回歸。
理論部分為何采用條件logistic回歸？這里以1:1配對(duì)設(shè)計(jì)為例，1:1配對(duì)設(shè)計(jì)的特點(diǎn)是對(duì)子內(nèi)部控制的混雜因素一致，不同對(duì)子的這些因素不同。從統(tǒng)計(jì)學(xué)角度來(lái)看，把每個(gè)配對(duì)試做一個(gè)成組病例對(duì)照研究，將配對(duì)層設(shè)置為啞變量。按照傳統(tǒng)的logistic回歸方法建立模型，估計(jì)每個(gè)自變量的比值比，這樣做有的主要困難有： 其一、每個(gè)配對(duì)層僅有2個(gè)觀察樣品； 其二、啞變量樹(shù)是對(duì)子數(shù)-1，大大增加了估計(jì)參數(shù)的數(shù)量。
為克服原有方法解決配對(duì)資料參數(shù)估計(jì)的缺陷，往往通過(guò)構(gòu)造特殊的條件似然函數(shù)，仍然采用極大似然估計(jì)法進(jìn)行參數(shù)估計(jì)。 構(gòu)造條件似然函數(shù)1:1病例對(duì)照的資料常常整理成如下格式。設(shè)有n對(duì)獨(dú)立的觀察，每個(gè)對(duì)子包含兩個(gè)人，第1個(gè)已經(jīng)患病，第2個(gè)沒(méi)有患病。自變量為X，第i層第一個(gè)人自變量記為Xi1Xi1，第二個(gè)人自變量記為Xi2Xi2。 | 配對(duì)號(hào) | 病例 |
| 對(duì)照 |
|
|---|
| X | Y | X | Y | | 1 | X11X11 | 1 | X10X10 | 0 | | 2 | X21X21 | 1 | X20X20 | 0 | | …… | …… | …… | …… | …… | | n | Xn1Xn1 | 1 | Xn0Xn0 | 0 |
任何一層中，第一個(gè)人患病的概率和未患病的概率分別為： π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β) 第二個(gè)人患病的概率與不患病的概率分別為： π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β) 現(xiàn)在假定同一層中的2個(gè)人中，只有1個(gè)人患病。在只有1個(gè)人患病的情況下，恰好第1個(gè)人患病而第2個(gè)人不患病的條件概率為 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個(gè)配對(duì)層的聯(lián)合概率，即似然函數(shù)為 L=P1×P2×…Pn=∏11+exp(?(XT1?XT0)β)L=P1×P2×…Pn=∏11+exp(?(X1T?X0T)β) 對(duì)上面的式子求對(duì)數(shù)，得到對(duì)數(shù)似然方程，通過(guò)迭代法求得偏回歸系數(shù)的值。假設(shè)檢驗(yàn)同樣采用似然比、score檢驗(yàn)以及wald檢驗(yàn)。 需要注意，上面的模型與經(jīng)典的logistic回歸有兩點(diǎn)不同： 與偏回歸系數(shù)相乘的是病例與對(duì)照相應(yīng)變量之差； 模型不含常數(shù)項(xiàng)。
實(shí)戰(zhàn)部分在多種統(tǒng)計(jì)軟件如SPSS、SAS等中，配對(duì)logistic回歸通常采用分層COX回歸模型來(lái)實(shí)現(xiàn)。這是因?yàn)椋謱覥ox回歸假設(shè)不同層之間的基線風(fēng)險(xiǎn)函數(shù)完全無(wú)關(guān)，只需估計(jì)協(xié)變量的偏回歸系數(shù)值。條件logistic回歸同樣如此，模型中不存在截距項(xiàng)，不同層之間有相同的偏回歸系數(shù)。可見(jiàn)用分層Cox回歸來(lái)擬合條件logistic回歸完全是一種取巧。 在調(diào)用分層Cox回歸時(shí)，會(huì)將病例算作發(fā)生終點(diǎn)事件，對(duì)照算作刪失。下面以一個(gè)例子來(lái)說(shuō)明： 為探討女性乳腺癌危險(xiǎn)因素，研究者在某市1996-1997年間確診的女性乳腺癌患者中隨機(jī)抽取350名病例，對(duì)每一病例以一名性別相同、年兩差別不超過(guò)±2.5歲的對(duì)照。收集的信息包括： | 變量 | 變量含義 | 值標(biāo)簽 |
|---|
| X1X1 | 文化程度 | 0：大專以下，1：大專及以上 | | X2X2 | 體質(zhì)指數(shù) | 0：小于等于27，1：大于27 | | X3X3 | 近年精神壓抑 | 0：無(wú)，1：有 | | X4X4 | 乳腺良性疾病史 | 0：無(wú)，1：有 | | X5X5 | 惡性腫瘤家族史 | 0：無(wú)，1：有 | | X6X6 | 初潮年齡 | 0：大于等于14歲，1：小于14歲 | | X7X7 | 哺乳史 | 0：有，1：無(wú) |
部分?jǐn)?shù)據(jù)如下，數(shù)據(jù)中以id號(hào)來(lái)標(biāo)識(shí)層變量，同一個(gè)id中第一個(gè)為病例，后一個(gè)為對(duì)照。： | 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 |
現(xiàn)在調(diào)用R程序進(jìn)行條件logistic回歸，使用的函數(shù)是survival包中的clogit函數(shù)，該函數(shù)實(shí)際上默認(rèn)調(diào)用了coxph函數(shù)。下面的程序首先讀取了數(shù)據(jù)，然后創(chuàng)建結(jié)局變量，進(jìn)而將分類變量轉(zhuǎn)化為因子。 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('無(wú)','有'))
df$x4<-factor(df$x4,levels = 0:1,labels = c('無(wú)','有'))
df$x5<-factor(df$x5,levels = 0:1,labels = c('無(wú)','有'))
df$x6<-factor(df$x6,levels = 0:1,labels = c('≤14歲','>14歲'))
df$x7<-factor(df$x7,levels = 0:1,labels = c('有','無(wú)'))
head(df)
|
clogit具有和其他建模函數(shù)同樣的用法，只需要使用strata參數(shù)指定分層變量即可
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無(wú) -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無(wú) 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不是顯著的，現(xiàn)在將其剔除， 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
|
剩下的變量都是顯著的，輸出的結(jié)果最后一部分是該模型與空模型的似然比檢驗(yàn)，wald檢驗(yàn)以及score檢驗(yàn)，結(jié)果同樣有統(tǒng)計(jì)學(xué)意義。 上述模型中，X1、X2、X3、X4、X5均有統(tǒng)計(jì)學(xué)意義，OR值均大于1。因此可以認(rèn)為高文化程度、肥胖、精神壓抑、乳腺良性疾病史和惡性腫瘤家族史是女性乳腺癌的危險(xiǎn)因素。
|
