2017-04-21 41 views
0

我自教過我幾天後,我陷入了一個cox迴歸分析。使用類別作爲新變量?

我管理使用cut()函數到每個連續的變量我有分成2個分類組。 現在我想知道是否可以在cox迴歸分析中單獨合併這些類別。例如:我有3個變量(A,B,C),每個變量由2個類別組成:A-,A + & B-,B + & C-,C +。現在,如果我運行命名變量的coxph(),我只能得到「+」類別的結果(據我瞭解是因爲「 - 」類別被用作參考組)。然而,由於C-對生存似乎有負面影響,因此我會更有興趣接受該類別的結果。

此外,我不知道是否有一種方法,我可以定義每個類別作爲一個新的組/變量,分別結合他們看到每個在生存的影響?或者這是不必要的?

編輯:我會盡力給予更具體的例子,我希望它的工作原理:)

#example data: 
test<-structure(list(A = c(8, 6, 42, 97, 55, 1, 5, 7, 55, 4), B = c(93, 9, 65, 2, 51, 89, 1, 1, 5, 62), C = c(58, 99, 100, 98, 92, 100, 99, 95, 81, 67), time = c(1.6, 34.6, 1.5, 35.8, 7.7, 38.6, 40.2, 4.7, 37.6, 8.6), event= c(1, 0, 0, 0, 1, 0, 0, 1, 0, 1))) 
mydata<-as.data.frame(test) 

它應該看起來像這樣:

mydata 
    A B C time status 
1 8 93 58 1.6  1 
2 6 9 99 34.6  0 
3 42 65 100 1.5  0 
4 97 2 98 35.8  0 
5 55 51 92 7.7  1 
6 1 89 100 38.6  0 
7 5 1 99 40.2  0 
8 7 1 95 4.7  1 
9 55 5 81 37.6  0 
10 4 62 67 8.6  1 

而且關於上述問題,這是我迄今所做的:

#load survival package 
library("survival") 

#define variables 
A <- c(mydata[1:10,1]) 
B <- c(mydata[1:10,2]) 
C <- c(mydata[1:10,3]) 
OS <- c(mydata[1:10,4]) 
Event <- c(mydata[1:10,5]) 

# dependent and independent variables 
y <- Surv(OS, Event) 
x <- cbind(A, B, C) 
mydata <- data.frame(cbind(x,y)) 

#Cox proportional hazard model, with the "raw data" 
coxph1 <- coxph(y~x,data=mydata, method="breslow") 
summary(coxph1) 

#categorising the variables 

CA=cut(mydata$A, br=c(-1,20,101), labels = c("[A-]", "[A+]")) 
CB=cut(mydata$B, br=c(-1,20,101), labels = c("[B-]", "[B+]")) 
CC=cut(mydata$C, br=c(-1,96,101), labels = c("[C-]", "[C+]")) 

#Cox regression, combined with cut intervals 
coxph2=coxph(y~CA+CB+CC, data=mydata, method="breslow") 
summary(coxph2) 

和預期的輸出結果是:

coxph(formula = y ~ x, data = mydata, method = "breslow") 

    n= 10, number of events= 4 

     coef exp(coef) se(coef)  z Pr(>|z|) 
xA 0.0001443 1.0001443 0.0238329 0.006 0.995 
xB 0.0104826 1.0105378 0.0211830 0.495 0.621 
xC -0.0497490 0.9514682 0.0383305 -1.298 0.194 

    exp(coef) exp(-coef) lower .95 upper .95 
xA 1.0001  0.9999 0.9545  1.048 
xB 1.0105  0.9896 0.9694  1.053 
xC 0.9515  1.0510 0.8826  1.026 

Concordance= 0.769 (se = 0.167) 
Rsquare= 0.29 (max possible= 0.799) 
Likelihood ratio test= 3.43 on 3 df, p=0.33 
Wald test   = 3.3 on 3 df, p=0.3476 
Score (logrank) test = 4.24 on 3 df, p=0.2364 


coxph(formula = y ~ CA + CB + CC, data = mydata, method = "breslow") 

    n= 10, number of events= 4 

      coef exp(coef) se(coef)  z Pr(>|z|) 
CA[A+] -1.036e+00 3.549e-01 1.262e+00 -0.821 0.412 
CB[B+] 4.294e-01 1.536e+00 1.274e+00 0.337 0.736 
CC[C+] -2.162e+01 4.095e-10 2.094e+04 -0.001 0.999 

     exp(coef) exp(-coef) lower .95 upper .95 
CA[A+] 3.549e-01 2.818e+00 0.0299  4.213 
CB[B+] 1.536e+00 6.509e-01 0.1266 18.653 
CC[C+] 4.095e-10 2.442e+09 0.0000  Inf 

Concordance= 0.904 (se = 0.165) 
Rsquare= 0.542 (max possible= 0.799) 
Likelihood ratio test= 7.8 on 3 df, p=0.05031 
Wald test   = 1.15 on 3 df, p=0.7653 
Score (logrank) test = 6.42 on 3 df, p=0.09288 
+3

歡迎SO。請閱讀[如何問一個好問題](http://stackoverflow.com/help/how-to-ask),[如何給出一個可重現的例子](http://stackoverflow.com/questions/ 5963269/how-to-make-a-great-r-reproducible-example)併發布您的預期輸出。 – Sotos

+0

我編輯它,以便現在有一個例子。你可能有一個想法,我可以如何使用單獨的組?或者這只是一個愚蠢的問題? – Alex

+0

看起來「錯了」。應該。是一個因素變量。 –

回答

0

首先,幾次有不必要的複雜代碼。測試數據可以寫爲:

mydata <- data.frame(A = c(8, 6, 42, 97, 55, 1, 5, 7, 55, 4), 
        B = c(93, 9, 65, 2, 51, 89, 1, 1, 5, 62), 
        C = c(58, 99, 100, 98, 92, 100, 99, 95, 81, 67), 
        time = c(1.6, 34.6, 1.5, 35.8, 7.7, 38.6, 40.2, 4.7, 37.6, 8.6), 
        event= c(1, 0, 0, 0, 1, 0, 0, 1, 0, 1)) 

而不是OS <- c(mydata[1:10,4])你可以寫OS <- mydata$time,但沒有必要把他們與我們的mydata。正如你所看到的,該模型是一樣的你的問題的第一個:

> (coxph1 <- coxph(Surv(time, event) ~ ., data=mydata, method="breslow")) 
Call: 
coxph(formula = Surv(time, event) ~ ., data = mydata, method = "breslow") 

     coef exp(coef) se(coef)  z p 
A 0.000144 1.000144 0.023833 0.01 1.00 
B 0.010483 1.010538 0.021183 0.49 0.62 
C -0.049749 0.951468 0.038331 -1.30 0.19 

Likelihood ratio test=3.43 on 3 df, p=0.33 
n= 10, number of events= 4 

關於你提到的個別相結合的協變量的問題 - 使用所有其他協變量~.手段。您可以使用~ A + B + C或其他任意組合專門指定它們。

關於改變參考類別 - 這僅需要具有> 2級的類別。係數的含義是類別和參考類別之間的差異。如果只有2個類別存在,比改變參考將給出相同的係數,位帶有' - '符號。 要更改因子參考類別使用relevel功能:

mydata$CA <- cut(mydata$A, br=c(-1,20,101), labels = c("[A-]", "[A+]")) 
mydata$CB <- cut(mydata$B, br=c(-1,20,101), labels = c("[B-]", "[B+]")) 
mydata$CC <- cut(mydata$C, br=c(-1,96,101), labels = c("[C-]", "[C+]")) 

mydata$CA <- relevel(mydata$CA, 2) 
> (coxph1 <- coxph(Surv(time, event) ~ CA, data=mydata, method="breslow")) 
Call: 
coxph(formula = Surv(time, event) ~ CA, data = mydata, method = "breslow") 

     coef exp(coef) se(coef) z p 
CA[A-] 0.559  1.749 1.158 0.48 0.63 

希望這有助於:)

+0

非常感謝! :) – Alex