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