时间戳

首页 > 解决方案 > R - cox.zph() 不返回 rho 值，p 值与示例不同

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

问题描述

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

我尝试运行以下代码：

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

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

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

我的 Schoenfeld 残差图

与示例相比：

sthda 版本相同的情节

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

提前致谢！

标签： rcox-regressionsurvival

解决方案

尝试更新您的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

推荐阅读