Почему этот код дает ошибочные P-значения?

Question

Я пытаюсь вычислить P -значения, связанные с точечными оценками, полученными из модели Кокса PH с изменяющимися во времени коэффициентами. Функция, которую я написал, не обеспечивает правильные значения P . Я проиллюстрирую это, используя данные NCCTG по раку легкого из пакета на выживание.

# Setup
require(survival)

# Effect of Karnofsky score, linear
fit <- coxph(Surv(time/365.25, status == 2) ~ ph.karno + tt(ph.karno), 
             lung, tt=function(x, t, ...) {x*t})

Функция:

# Same function but now with a P-value in the output
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 <- pnorm(-abs(z.at.time))

  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

}

# Point estimates after 1.05*365.25 = 383.5 days of follow-up
calculate.timeDependentHazard.P(fit,1.05)

Выход:

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

Очевидно, значение P должно быть> 0,05, но почему-то это не так. P -значения, рассчитанные с помощью этого подхода, кажутся слишком низкими. Кто-нибудь, кто может обнаружить недостаток?

Answer 1

Кажется, вы хотите двустороннюю альтернативу, поэтому умножьте pnorm(-abs(z.at.time)) на два.То есть сделать 2*pnorm(-abs(z.at.time)).