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

Question

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

我试过运行以下代码：

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) )

这给了我以下 return:

         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

但是，其他地方的示例（包括我正在尝试遵循的示例 here）return 带有 "rho" 列的 table，如下所示：

            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

此外，无论它抛出什么，似乎也在改变我的卡方值和 p 值。

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

My Schoenfeld residual plots

对比例子：

sthda version of the same plots

这些问题严重阻碍了我的团队在当前项目上的工作，我们将不胜感激提供的任何帮助。

提前致谢！

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 - mean(ttimes)] 的差异和 2. 之间的皮尔逊相关性。拟合值乘以 schoenfeld 残差（按死亡人数缩放）。从代码..

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

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

R - cox.zph() not returning rho values, p-values differ from examples

r

cox-regression

survival