0

我的目标是预测（从下面的拟合模型预测新观察的累积风险）从时间尺度 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”值。有人可以指出我的想法或代码中是否有任何错误

如果需要进一步澄清，请告诉我。

4

2 回答 2

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
于 2020-09-02T14:45:51.530 回答
0

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

D$coefficients["X"] <- 0

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

让我知道是否有任何其他适当的方法来解决它。

于 2020-09-22T11:32:41.477 回答