在 r 中使用 predict.coxph 预测基线累积危害
Predicting baseline cumulative hazard using predict.coxph in r
我的目标是预测(从下面的拟合模型预测新观察的累积风险)从拟合模型的时间尺度 0 到开始时间的累积风险值。
我已经使用2次(不等于0的开始时间和结束时间)拟合了cox模型。这样我就可以找到结束时间的累积危险(即从 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”值。谁能指出我的想法或代码是否有错误
如果需要进一步说明,请告诉我。
有两个问题。首先,您 运行 遇到了一个问题,因为当您指定 ~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
我自己找到了答案,这只是一个小技巧,但我不确定是否总能奏效。
如果我在 predict()
函数之前使用以下行:
D$coefficients["X"] <- 0
但是,我得到了使用不接受开始时间(或一次两个变量)的 nelsonaalen()
函数检查的正确值
如果有任何其他正确的解决方法,请告诉我。
我的目标是预测(从下面的拟合模型预测新观察的累积风险)从拟合模型的时间尺度 0 到开始时间的累积风险值。
我已经使用2次(不等于0的开始时间和结束时间)拟合了cox模型。这样我就可以找到结束时间的累积危险(即从 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”值。谁能指出我的想法或代码是否有错误
如果需要进一步说明,请告诉我。
有两个问题。首先,您 运行 遇到了一个问题,因为当您指定 ~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
我自己找到了答案,这只是一个小技巧,但我不确定是否总能奏效。
如果我在 predict()
函数之前使用以下行:
D$coefficients["X"] <- 0
但是,我得到了使用不接受开始时间(或一次两个变量)的 nelsonaalen()
函数检查的正确值
如果有任何其他正确的解决方法,请告诉我。