8

我有一个 Cox 比例风险模型，使用 R 中的以下代码来预测死亡率。添加协变量 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 个百分位数。我该怎么做？

我可以列出我尝试过的所有其他内容，但我不想混淆任何人！

非常感谢。

4

2 回答 2

6

这是取自这篇论文的一个例子。

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. 如本文档中所述，该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% 的频段。

于 2015-08-18T16:27:58.493 回答
3

X 的第 25 和第 75 个百分位数的估计生存曲线请求首先需要确定这些百分位数并指定数据框中所有其他协变量的值，以用作 survfit 的新数据参数：

可以使用 Fox 网站上其他响应者建议的数据，尽管在我的机器上它需要构建一个url-object：

 url <- url("http://socserv.mcmaster.ca/jfox/Books/Companion/data/Rossi.txt")
 Rossi <- read.table(url, header=TRUE)

这可能不是这个 wquestion 的最佳示例，但它确实有一个数字变量，我们可以计算四分位数：

> summary(Rossi$prio)
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  0.000   1.000   2.000   2.984   4.000  18.000

所以这将是模型拟合和 survfit 调用：

 mod.allison <- coxph(Surv(week, arrest) ~ 
                         fin + age + race + prio ,
                         data=Rossi)
 prio.fit <- survfit(mod.allison, 
                     newdata= data.frame(fin="yes", age=30, race="black", prio=c(1,4) ))
 plot(prio.fit, col=c("red","blue"))

在此处输入图像描述

于 2015-08-19T20:02:18.100 回答