为何Cox时变系数模型自定义函数计算的P值存在偏差？

我一眼就看到问题出在P值的计算逻辑上——你现在用的是单侧检验的P值，而我们在Cox比例风险模型中，针对系数是否为0的假设检验通常用的是双侧检验，这就是为什么你的P值普遍偏低的原因。

具体错误点 #

在你的函数里，这一行代码计算的是单侧P值（对应备择假设H₁: β < 0）：

p.value <- pnorm(-abs(z.at.time))

但我们需要的是双侧P值（对应备择假设H₁: β ≠ 0），所以应该把结果乘以2，修正为：

p.value <- 2 * pnorm(-abs(z.at.time))

验证修正结果 #

用你提供的例子，随访时间1.05年时，z值的绝对值约为1.67，单侧P值是0.047，乘以2之后得到≈0.094，这就大于0.05了，和你预期的一致。

其他步骤的正确性确认 #

顺便检查了函数里的其他逻辑：

• 协方差矩阵的提取是正确的，通过索引定位ph.karno和tt(ph.karno)的子矩阵没问题；
• 时变系数的log(HR)和方差的二次型计算逻辑正确；
• 置信区间的计算符合常规方法，指数化后得到HR的区间也没问题。

修正后的完整函数 #

calculate.timeDependentHazard.P <- function(model,time) {
  index.1 <- which(names(model$coef)=="ph.karno")
  index.2 <- which(names(model$coef)=="tt(ph.karno)")
  coef <- model$coef[c(index.1,index.2)]
  var <- rbind(c(model$var[index.1,index.1],model$var[index.1,index.2]),
               c(model$var[index.2,index.1],model$var[index.2,index.2]))
  var.at.time <- t(c(1,time)) %*% var %*% c(1,time)
  hazard.at.time <- t(c(1,time)) %*% coef
  lower.95 <- hazard.at.time - 1.96*sqrt(var.at.time)
  upper.95 <- hazard.at.time + 1.96*sqrt(var.at.time)
  z.at.time <- hazard.at.time/(sqrt(var.at.time))
  p.value <- 2 * pnorm(-abs(z.at.time))  # 修正为双侧P值
  results <- c(exp(c(hazard.at.time,lower.95,upper.95)),p.value)
  names(results) <- c("hazard ratio","95% lower","95% upper","P.value")
  options(scipen = 999)
  results
}

运行修正后的函数，会得到符合预期的结果：

calculate.timeDependentHazard.P(fit,1.05)
#  hazard ratio   95% lower   95% upper      P.value 
#   0.98913256    0.97654719    1.00188013    0.09442684 

内容的提问来源于stack exchange，提问作者Dion Groothof