R - cox.zph() 不返回 rho 值，p 值与示例不同

Question

当我在Cox模型上运行cox.zph时，我得到的回报类型与我在其他地方看到的相同。

我尝试运行以下代码：

library(survival) #version 3.1-7
library(survminer) #version 0.4.6
res.cox <- coxph(Surv(time, status) ~ age + sex + wt.loss, data = lung)
#lung data is in the survival package and loads from there.
( test.ph <- cox.zph(res.cox) )

这给了我以下回报：

chisq df    p
age     0.5077  1 0.48
sex     2.5489  1 0.11
wt.loss 0.0144  1 0.90
GLOBAL  3.0051  3 0.39

但是，其他地方的示例(包括我在这里尝试遵循的示例(返回一个带有"rho"列的表，如下所示：

rho chisq     p
age     -0.0483 0.378 0.538
sex      0.1265 2.349 0.125
wt.loss  0.0126 0.024 0.877
GLOBAL       NA 2.846 0.416

此外，无论它抛出什么，似乎也会改变我的 chi-sq 和 p 值。

此外，当我随后尝试使用ggcoxzph(test.ph)绘制勋菲尔德残差时，我得到以下图：

我的勋菲尔德残差图

与示例相比：

相同情节的STHDA版本

这些问题对我的团队在当前项目上的努力造成了巨大的阻碍，任何提供的帮助将不胜感激。

提前感谢！

Answer 1

尝试更新survival包的版本。它应该对我有用(打印方法显示"rho"(。

我使用的是 2.43-3 版

编辑：42是对的，我是从过时的镜像安装的。我刚刚更新到最新版本。

当然，一种解决方案是将安装回滚到旧版本。但假设你不想这样做...

看来你需要自己计算你想要什么。我不知道这是否是最简洁的方法，但我只是将旧版本包的源代码改编成一个函数，您可以在其中传递适合的对象和 cox.zph 测试对象。

library(survival) #version 3.1-7
library(survminer) #version 0.4.6
res.cox <- coxph(Surv(time, status) ~ age + sex + wt.loss, data = lung)
#lung data is in the survival package and loads from there.
test.ph <- cox.zph(res.cox) 
your_func <- function(fit, transform = "km", new_cox.zph = NULL) {
sresid <- resid(fit, "schoenfeld")
varnames <- names(fit$coefficients)
nvar <- length(varnames)
ndead <- length(sresid)/nvar
if (nvar == 1) {
times <- as.numeric(names(sresid)) 
} else {
times <- as.numeric(dimnames(sresid)[[1]])
}
if (is.character(transform)) {
tname <- transform
ttimes <- switch(transform, identity = times, rank = rank(times), 
log = log(times), km = {
temp <- survfitKM(factor(rep(1, nrow(fit$y))), 
fit$y, se.fit = FALSE)
t1 <- temp$surv[temp$n.event > 0]
t2 <- temp$n.event[temp$n.event > 0]
km <- rep(c(1, t1), c(t2, 0))
if (is.null(attr(sresid, "strata"))) 1 - km else (1 - 
     km[sort.list(sort.list(times))])
}, stop("Unrecognized transform"))
}
else {
tname <- deparse(substitute(transform))
if (length(tname) > 1) 
tname <- "user"
ttimes <- transform(times)
}
xx <- ttimes - mean(ttimes)
r2 <- sresid %*% fit$var * ndead
test <- xx %*% r2
corel <- c(cor(xx, r2))
cbind(rho = c(corel,NA), new_cox.zph$table)
}

调用函数

your_func(fit = res.cox, new_cox.zph = test.ph)
rho      chisq df         p
age     -0.04834523 0.50774082  1 0.4761185
sex      0.12651872 2.54891522  1 0.1103700
wt.loss  0.01257876 0.01444092  1 0.9043482
GLOBAL           NA 3.00505434  3 0.3908466

更新

我不知道版本之间的 chisq 和 p 计算有何不同。您可能需要查看发行说明以了解版本之间的更改。但就您的目的而言，我不确定"rho"是否会有差异，这似乎只是 1 之间的梨逊相关性。与观测时间和平均时间[ttimes - 平均值(ttimes(]和2的区别。拟合值乘以勋菲尔德残差(按死数缩放(。从代码..

sresid <- resid(fit, "schoenfeld")
xx <- ttimes - mean(ttimes)
r2 <- sresid %*% fit$var * ndead
test <- xx %*% r2
corel <- c(cor(xx, r2))
##corel is rho