使用类别作为新变量？
Using categories as new variables?
几天以来我一直在自学 R，但我一直坚持使用 cox 回归分析。
我设法使用 cut() 函数将每个连续变量分成 2 个分类组。
现在我想知道我是否可以在 cox 回归分析中单独组合这些类别。
例如：我有 3 个变量（A、B、C），每个变量由 2 个类别组成：A-、A+ & B-、B+ & C-、C+。现在，如果我 运行 coxph() 命名变量，我只会得到“+”类别的结果（据我了解，这是因为“-”类别用作参考组）。
但是，由于 C- 似乎对生存有负面影响，我更感兴趣的是收到该类别的结果。
另外我想知道是否有一种方法可以将每个类别定义为一个新的 group/variable 并将它们单独组合以查看每个类别对生存的影响？或者这是不必要的？
编辑：我会试着给出一个更具体的例子，我希望它有用:)
#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
首先，有几次出现了不必要的复杂代码。测试数据可以这样写：
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 <- mydata$time 而不是 OS <- c(mydata[1:10,4])，但没有必要将它们作为我们的 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
希望这对您有所帮助:)
几天以来我一直在自学 R，但我一直坚持使用 cox 回归分析。
我设法使用 cut() 函数将每个连续变量分成 2 个分类组。
现在我想知道我是否可以在 cox 回归分析中单独组合这些类别。
例如：我有 3 个变量（A、B、C），每个变量由 2 个类别组成：A-、A+ & B-、B+ & C-、C+。现在，如果我 运行 coxph() 命名变量，我只会得到“+”类别的结果（据我了解，这是因为“-”类别用作参考组）。
但是，由于 C- 似乎对生存有负面影响，我更感兴趣的是收到该类别的结果。
另外我想知道是否有一种方法可以将每个类别定义为一个新的 group/variable 并将它们单独组合以查看每个类别对生存的影响？或者这是不必要的？
编辑：我会试着给出一个更具体的例子，我希望它有用:)
#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
首先，有几次出现了不必要的复杂代码。测试数据可以这样写：
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 <- mydata$time 而不是 OS <- c(mydata[1:10,4])，但没有必要将它们作为我们的 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
希望这对您有所帮助:)