如何在 R 中的 ggplot 中绘制风险比 + CI 随着时间的推移生存数据?
How to plot the Hazard Ratio + CI over time of survival data in ggplot in R?
背景
我想绘制生存数据集随时间变化的风险比,包括其置信区间。作为示例,我将从 survival
包中获取一个简化的数据集:冒号数据集。
library(survival)
library(tidyverse)
# Colon survival dataset
data <- colon %>%
filter(etype == 2) %>%
select(c(id, rx, status, time)) %>%
filter(rx == "Obs" | rx == "Lev+5FU") %>%
mutate(rx = factor(rx))
数据集包含接受治疗的患者(即“Lev+5FU”)和未接受治疗的患者(即“Obs”)。生存曲线如下:
fit <- survfit(Surv(time, status) ~ rx, data = data )
plot(fit)
尝试
使用 cox.zph
函数,您可以绘制 cox 模型的风险比。
cox <- coxph(Surv(time, status) ~ rx, data = data)
plot(cox.zph(cox))
但是,我想使用 ggplot
.
绘制此生存数据集的风险比,包括 95% CI
问题
- 如何从这个 cox.zph 对象中提取风险比数据和 95% CIs 以将它们绘制在
ggplot
中?
- 是否有其他
R
软件包能够以更方便的方式实现同样的功能?
survminer
包 will do this for you:
library(survminer)
ggcoxzph(cox.zph(cox))
注意:认清 Dion Groothof 的修正很重要。线条和 CIs 并不是真正的风险比。它们是随时间变化 log-hazard-ratios 的估计值和界限。您需要取幂才能获得 HR。
这些值在 cox.zph
:
返回的结果中
str(cox.zph(cox))
#----------------------
List of 7
$ table : num [1:2, 1:3] 1.188 1.188 1 1 0.276 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : chr [1:2] "rx" "GLOBAL"
.. ..$ : chr [1:3] "chisq" "df" "p"
$ x : num [1:291] 0 0.00162 0.00323 0.00485 0.00646 ...
$ time : num [1:291] 23 34 45 52 79 113 125 127 138 141 ...
$ y : num [1:291, 1] 2.09 2.1 2.1 2.1 2.11 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : chr [1:291] "23" "34" "45" "52" ...
.. ..$ : chr "rx"
$ var : num [1, 1] 4.11
$ transform: chr "km"
$ call : language cox.zph(fit = cox)
- attr(*, "class")= chr "cox.zph"
要使用任何范例(基础、点阵或 ggplot2)绘制图表,请使用 time
作为 x 轴,使用 x
作为实线,使用“点”处的 y “
z <- cox.zph(cox)
ggdf <- data.frame( unclass(z)[c("time", "x","y")])
ggplot(data=ggdf, aes(x=time, y=-x))+
geom_line()+ ylim(range(z$y))+
geom_point(aes(x=time,y=z$y) )
要获得 CI 查看 getAnywhere(plot.cox.zph)
xx <- x$x
yy <- x$y
df <- max(df)
nvar <- ncol(yy)
pred.x <- seq(from = min(xx), to = max(xx), length = nsmo)
#------------
if (se) {
bk <- backsolve(qmat$qr[1:df, 1:df], diag(df))
xtx <- bk %*% t(bk)
seval <- ((pmat %*% xtx) * pmat) %*% rep(1, df)
temp <- 2 * sqrt(x$var[i, i] * seval)
yup <- yhat + temp
ylow <- yhat - temp
yr <- range(yr, yup, ylow)
#---------------
if (se) {
lines(pred.x, exp(yup), col = col[2], lty = lty[2],
lwd = lwd[2])
lines(pred.x, exp(ylow), col = col[2], lty = lty[2],
lwd = lwd[2])
}
背景
我想绘制生存数据集随时间变化的风险比,包括其置信区间。作为示例,我将从 survival
包中获取一个简化的数据集:冒号数据集。
library(survival)
library(tidyverse)
# Colon survival dataset
data <- colon %>%
filter(etype == 2) %>%
select(c(id, rx, status, time)) %>%
filter(rx == "Obs" | rx == "Lev+5FU") %>%
mutate(rx = factor(rx))
数据集包含接受治疗的患者(即“Lev+5FU”)和未接受治疗的患者(即“Obs”)。生存曲线如下:
fit <- survfit(Surv(time, status) ~ rx, data = data )
plot(fit)
尝试
使用 cox.zph
函数,您可以绘制 cox 模型的风险比。
cox <- coxph(Surv(time, status) ~ rx, data = data)
plot(cox.zph(cox))
但是,我想使用 ggplot
.
问题
- 如何从这个 cox.zph 对象中提取风险比数据和 95% CIs 以将它们绘制在
ggplot
中? - 是否有其他
R
软件包能够以更方便的方式实现同样的功能?
survminer
包 will do this for you:
library(survminer)
ggcoxzph(cox.zph(cox))
注意:认清 Dion Groothof 的修正很重要。线条和 CIs 并不是真正的风险比。它们是随时间变化 log-hazard-ratios 的估计值和界限。您需要取幂才能获得 HR。
这些值在 cox.zph
:
str(cox.zph(cox))
#----------------------
List of 7
$ table : num [1:2, 1:3] 1.188 1.188 1 1 0.276 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : chr [1:2] "rx" "GLOBAL"
.. ..$ : chr [1:3] "chisq" "df" "p"
$ x : num [1:291] 0 0.00162 0.00323 0.00485 0.00646 ...
$ time : num [1:291] 23 34 45 52 79 113 125 127 138 141 ...
$ y : num [1:291, 1] 2.09 2.1 2.1 2.1 2.11 ...
..- attr(*, "dimnames")=List of 2
.. ..$ : chr [1:291] "23" "34" "45" "52" ...
.. ..$ : chr "rx"
$ var : num [1, 1] 4.11
$ transform: chr "km"
$ call : language cox.zph(fit = cox)
- attr(*, "class")= chr "cox.zph"
要使用任何范例(基础、点阵或 ggplot2)绘制图表,请使用 time
作为 x 轴,使用 x
作为实线,使用“点”处的 y “
z <- cox.zph(cox)
ggdf <- data.frame( unclass(z)[c("time", "x","y")])
ggplot(data=ggdf, aes(x=time, y=-x))+
geom_line()+ ylim(range(z$y))+
geom_point(aes(x=time,y=z$y) )
要获得 CI 查看 getAnywhere(plot.cox.zph)
xx <- x$x
yy <- x$y
df <- max(df)
nvar <- ncol(yy)
pred.x <- seq(from = min(xx), to = max(xx), length = nsmo)
#------------
if (se) {
bk <- backsolve(qmat$qr[1:df, 1:df], diag(df))
xtx <- bk %*% t(bk)
seval <- ((pmat %*% xtx) * pmat) %*% rep(1, df)
temp <- 2 * sqrt(x$var[i, i] * seval)
yup <- yhat + temp
ylow <- yhat - temp
yr <- range(yr, yup, ylow)
#---------------
if (se) {
lines(pred.x, exp(yup), col = col[2], lty = lty[2],
lwd = lwd[2])
lines(pred.x, exp(ylow), col = col[2], lty = lty[2],
lwd = lwd[2])
}