文章/答案/技术大牛

发布

社区首页 >问答首页 >R- cox.zph()不返回rho值，p-值与示例不同。

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

Stack Overflow用户

提问于 2019-11-23 20:12:40

回答 1查看 1.2K关注 0票数 1

当我在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)绘制Schoenfeld残差时，我得到了以下图：

我的Schoenfeld残差图

相对于示例：

相同情节的sthda版本

这些问题正在造成一个巨大的障碍，我的小组在当前的项目，并提供任何帮助将是非常感谢的。

提前感谢！

survival

cox-regression

回答 1

Stack Overflow用户

回答已采纳

发布于 2019-11-23 20:36:31

尝试更新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“中是否会有差异，这似乎只是一个pearson相关性。1.观测时间和平均时间的差异--平均( time )和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

票数 0

页面原文内容由Stack Overflow提供。腾讯云小微IT领域专用引擎提供翻译支持

原文链接：

https://stackoverflow.com/questions/59011818

复制

相似问题

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

回答 1

Stack Overflow用户

社区

活动

圈层

关于

腾讯云开发者

热门产品

热门推荐

更多推荐

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

回答 1

Stack Overflow用户

社区

活动

圈层

关于

腾讯云开发者

热门产品

热门推荐

更多推荐

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