在 r 中使用 predict.coxph 预测基线累积危险答案

【问题标题】：Predicting baseline cumulative hazard using predict.coxph in r在 r 中使用 predict.coxph 预测基线累积危险
【发布时间】：2020-09-02 08:48:41
【问题描述】：

我的目标是预测（从下面的拟合模型预测新观察的累积风险）从时间尺度 0 到拟合模型的开始时间的累积风险值。

我使用 2 次拟合 cox 模型（开始时间不等于 0 和结束时间）。那么我可以找到结束时间的累积风险（即从 0 到结束时间的累积风险，我已经使用相同的拟合模型计算）和开始时间的累积风险（即从 0 到结束时间，我想在这里计算）这将最终给出每次观察的开始和结束时间之间的 cum haz。

为了获得预期的事件数量，我使用了predict(coxph(), newdata, type= "expected")。

我使用的数据如下：

N <- 10^4 # population
H <- within(data.frame(start_time=runif(N, 0, 50), x1=rnorm(N, 2, 1), x2=rnorm(N, -2, 1)), {
  lp <-   0.05*x1 + 0.2*x2 
  Tm <- qweibull(runif(N,pweibull(start_time,shape = 7.5, scale = 84*exp(-lp/7.5)),1), shape=7.5, scale=84*exp(-lp/7.5))
  Cens1 <- 100
  event_time <- pmin(Tm,Cens1)
  status <- as.numeric(event_time == Tm)})

预测的代码是：

H$X <- rep(1,nrow(H))
D = coxph(Surv(start_time, event_time, status) ~ X, data =  H, x = TRUE )
pred2 <- predict(D, newdata = data.frame(start_time = rep(0,nrow(H)),event_time = H$start_time, status = rep(0,nrow(H)), X = rep(1, nrow(H))), type = "expected")

但pred2 只会产生“NA”值。有人可以指出我的想法或代码中是否有任何错误

如果需要进一步说明，请告诉我。

【问题讨论】：

您的代码不可重现。请提供d。
对不起，我已经更正了代码，我认为它现在可以重现了。你现在可以试试吗。让我知道是否需要更正
Cox 模型的重点是估计相对危害。它始终与具有协变量平均值的假设案例相关。由于没有协变量，您要求它估算相对危险，而没有可比较的风险。如果您想要的只是 KM 曲线，那么 cph 是错误的函数。

标签： r predict survival-analysis cox-regression hazard

【解决方案1】：

有两个问题。首先，您遇到了一个问题，因为当您指定 ~1 时，这意味着拟合没有系数的仅截距模型。所以你所有的预测都会给你一个价值？

library(survival)
D <- coxph(Surv(H$start_time, H$event_time, H$status) ~ 1, x = TRUE )
names(D)
 [1] "loglik"            "linear.predictors" "method"           
 [4] "residuals"         "n"                 "nevent"           
 [7] "terms"             "assign"            "concordance"      
[10] "x"                 "y"                 "timefix"          
[13] "formula"           "call"  

table(predict(D))

    0 
10000

我认为这没有多大意义，因此您会遇到所有错误。因此，您需要使用回归中使用的自变量进行预测，例如：

D <- coxph(Surv(start_time,event_time,status) ~ x1+x2, data=H )
pred2 <- predict(D, newdata = data.frame(t_0 = rep(0,nrow(H)),time = H$start_time, event_M = rep(0,nrow(H)), X = rep(1, nrow(H))), type = "expected")

predict(D,newdata=data.frame(x1=runif(10,0,1),x2=runif(10,-1,1)))
        1         2         3         4         5         6         7         8 
0.3033206 0.4213120 0.3952827 0.3879701 0.4798670 0.2170032 0.3385253 0.4141698 
        9        10 
0.3690579 0.4128084

当您拟合一个所有 X=1 的模型时，这会为您提供所有 NA，因为已经存在一个截距，这使得该变量变得多余。您可以检查：

H$X = 1
D <- coxph(Surv(start_time, event_time, status) ~ X,data=H)

Call:
coxph(formula = Surv(start_time, event_time, status) ~ X, data = H)

  coef exp(coef) se(coef)  z  p
X   NA        NA        0 NA NA

仅当 X 是拟合数据中的实际变量时才有效，因此我使用具有 2 个协变量的示例：

H$X = runif(nrow(H))
D <- coxph(Surv(start_time, event_time, status) ~ X + x1,data=H)

您可以通过例如将 X 固定为 1 并改变 x1 来进行预测：

predict(D,newdata=data.frame(X=1,x1=c(0.1,0.2,0.3)))
         1          2          3 
-0.1132548 -0.1084592 -0.1036637

或 X 在 2：

predict(D,newdata=data.frame(X=2,x1=c(0.1,0.2,0.3)))
                 1          2          3 
-0.1579480 -0.1531524 -0.1483568

【讨论】：

在问题中，我想在没有协变量影响的情况下估计基线累积风险，当我添加一个变量H$X = rep(1, nrow(H)) 并使用D = coxph(Surv(start_time, event_time, status) ~ X, data = H, x = TRUE ) 拟合模型时，代码运行没有任何错误但我得到“NA”值作为结果。我不确定我的概念是否正确。 PS：我刚刚修改了代码，我认为您不会在两者之间遇到错误（我已经在问题中更详细地解释了我想要的内容）
在你放入模型的数据中，X全为1。这与拟合模型 coxph(Surv(start_time, event_time, status) ~ 1 有何不同？你现在遇到了一个错误，因为你也有拦截，使 X 冗余
你真的能用合理的数据拟合模型吗？我已经更新了我的答案。它似乎不再是一个编程问题。也许您可以在发布之前先检查一下？ socialsciences.mcmaster.ca/jfox/Books/Companion/appendices/…

【解决方案2】：

我自己找到了答案，这只是一个快速技巧，我不确定它是否会一直有效。如果我在predict() 函数之前使用以下行：

D$coefficients["X"] <- 0

但是，我得到了使用 nelsonaalen() 函数检查的正确值，该函数不接受开始时间（或一次两个变量）

如果有其他合适的方法可以解决，请告诉我。

【讨论】：