我在R中使用以下代码建立了Cox比例风险模型,该模型可预测死亡率。添加协变量A,B和C只是为了避免混淆(即年龄,性别,种族),但我们对预测变量X确实很感兴趣.X是连续变量。
cox.model <- coxph(Surv(time, dead) ~ A + B + C + X, data = df)
现在,我很难为此绘制Kaplan-Meier曲线。我一直在寻找如何创建这个数字的方法,但是运气并不好。我不确定是否可以为Cox模型绘制Kaplan-Meier? Kaplan-Meier是否针对我的协变量进行调整,或者不需要它们?
下面是我尝试过的内容,但有人告诉我这是不正确的。
plot(survfit(cox.model), xlab = 'Time (years)', ylab = 'Survival Probabilities')
我还试图绘制一个数字,显示死亡的累积危害。我不知道我是否做对了,因为我尝试了几种不同的方法并获得了不同的结果。理想情况下,我想绘制两条线,一条线显示X的第75个百分点的死亡风险,另一条线显示X的25个百分点的死亡风险。我该怎么做?
我可以列出我尝试过的所有其他内容,但是我不想让任何人感到困惑!
非常感谢。
最佳答案
这是取自this paper的示例。
url <- "http://socserv.mcmaster.ca/jfox/Books/Companion/data/Rossi.txt"
Rossi <- read.table(url, header=TRUE)
Rossi[1:5, 1:10]
# week arrest fin age race wexp mar paro prio educ
# 1 20 1 no 27 black no not married yes 3 3
# 2 17 1 no 18 black no not married yes 8 4
# 3 25 1 no 19 other yes not married yes 13 3
# 4 52 0 yes 23 black yes married yes 1 5
# 5 52 0 no 19 other yes not married yes 3 3
mod.allison <- coxph(Surv(week, arrest) ~
fin + age + race + wexp + mar + paro + prio,
data=Rossi)
mod.allison
# Call:
# coxph(formula = Surv(week, arrest) ~ fin + age + race + wexp +
# mar + paro + prio, data = Rossi)
#
#
# coef exp(coef) se(coef) z p
# finyes -0.3794 0.684 0.1914 -1.983 0.0470
# age -0.0574 0.944 0.0220 -2.611 0.0090
# raceother -0.3139 0.731 0.3080 -1.019 0.3100
# wexpyes -0.1498 0.861 0.2122 -0.706 0.4800
# marnot married 0.4337 1.543 0.3819 1.136 0.2600
# paroyes -0.0849 0.919 0.1958 -0.434 0.6600
# prio 0.0915 1.096 0.0286 3.194 0.0014
#
# Likelihood ratio test=33.3 on 7 df, p=2.36e-05 n= 432, number of events= 114
请注意,该模型使用
fin, age, race, wexp, mar, paro, prio
预测arrest
。如this document中所述,survfit()
函数将Kaplan-Meier估计用于生存率。plot(survfit(mod.allison), ylim=c(0.7, 1), xlab="Weeks",
ylab="Proportion Not Rearrested")
我们得到了一个生存率图(置信区间为95%)。对于累积危害率,您可以执行
# plot(survfit(mod.allison)$cumhaz)
但这没有给出置信区间。但是,不用担心!我们知道H(t)= -ln(S(t)),我们有S(t)的置信区间。我们要做的就是
sfit <- survfit(mod.allison)
cumhaz.upper <- -log(sfit$upper)
cumhaz.lower <- -log(sfit$lower)
cumhaz <- sfit$cumhaz # same as -log(sfit$surv)
然后画这些
plot(cumhaz, xlab="weeks ahead", ylab="cumulative hazard",
ylim=c(min(cumhaz.lower), max(cumhaz.upper)))
lines(cumhaz.lower)
lines(cumhaz.upper)
您需要使用
survfit(..., conf.int=0.50)
来获得75%和25%的频段,而不是97.5%和2.5%的频段。关于r - 绘制Kaplan-Meier进行Cox回归,我们在Stack Overflow上找到一个类似的问题:https://stackoverflow.com/questions/32076310/