Perigo da linha de base de Cox

19

Digamos que eu tenha um conjunto de dados de "cateter renal". Estou tentando modelar uma curva de sobrevivência usando um modelo de Cox. Se eu considerar um modelo de Cox: preciso da estimativa do risco da linha de base. Usando a função interna do pacote R , posso fazer isso da seguinte maneira:

h(t,Z)=h0exp(bZ),
survivalbasehaz()
library(survival)

data(kidney)
fit <- coxph(Surv(time, status) ~ age , kidney)
basehaz(fit)

Mas se eu quiser escrever uma função passo a passo do risco da linha de base para uma determinada estimativa de parâmetro, bcomo posso proceder? Eu tentei:

bhaz <- function(beta, time, status, x) {

    data <- data.frame(time,status,x)
    data <- data[order(data$time), ]
    dt   <- data$time
    k    <- length(dt)
    risk <- exp(data.matrix(data[,-c(1:2)]) %*% beta)
    h    <- rep(0,k)

    for(i in 1:k) {
        h[i] <- data$status[data$time==dt[i]] / sum(risk[data$time>=dt[i]])          
    }

    return(data.frame(h, dt))
}

h0 <- bhaz(fit$coef, kidney$time, kidney$status, kidney$age)

Mas isso não dá o mesmo resultado que basehaz(fit). Qual é o problema?

Dihan
fonte
@gung você poderia ajudar com esta pergunta ? Lutei durante dois dias ...
Haitao Du

Respostas:

21

Aparentemente, basehaz()na verdade calcula uma taxa de risco cumulativa, em vez da própria taxa de risco. A fórmula é como se com h 0 ( y ( l ) ) = d ( l )

H^0(t)=y(l)th^0(y(l)),
ondey(1)<y(2)<otedenota os tempos distintos dos eventos,d(l)é o número de eventos emy(l), eR(y(l))é o risco definido emy(l)
h^0(y(l))=d(l)jR(y(l))exp(xjβ)
y(1)<y(2)<d(l)y(l)R(y(l))y(l)contendo todos os indivíduos ainda suscetíveis ao evento em .y(l)

Vamos tentar isso. (O código a seguir existe apenas para ilustração e não pretende ser muito bem escrito.)

#------package------
library(survival)

#------some data------
data(kidney)

#------preparation------
tab <- data.frame(table(kidney[kidney$status == 1, "time"])) 
y <- as.numeric(levels(tab[, 1]))[tab[, 1]] #ordered distinct event times
d <- tab[, 2]                               #number of events

#------Cox model------
fit<-coxph(Surv(time, status)~age, data=kidney)

#------cumulative hazard obtained from basehaz()------
H0 <- basehaz(fit, centered=FALSE)
H0 <- H0[H0[, 2] %in% y, ] #only keep rows where events occurred

#------my quick implementation------
betaHat <- fit$coef

h0 <- rep(NA, length(y))
for(l in 1:length(y))
{
  h0[l] <- d[l] / sum(exp(kidney[kidney$time >= y[l], "age"] * betaHat))
}

#------comparison------
cbind(H0, cumsum(h0))

saída parcial:

       hazard time cumsum(h0)
1  0.01074980    2 0.01074980
5  0.03399089    7 0.03382306
6  0.05790570    8 0.05757756
7  0.07048941    9 0.07016127
8  0.09625105   12 0.09573508
9  0.10941921   13 0.10890324
10 0.13691424   15 0.13616338

Suspeito que a pequena diferença possa ser devida à aproximação da probabilidade parcial coxph()devido a laços nos dados ...

ocram
fonte
Muito obrigado. Sim, há uma pequena diferença no método de aproximação. Mas há 76 pontos no tempo, se eu quiser encontrar o risco da linha de base para cada ponto no tempo. O que eu posso fazer? Que tipo de modificação no código R é necessário?
Dihan
1
O risco discretizado é zero, exceto nos horários dos eventos. Isso realmente dá a maior contribuição para a probabilidade se uma função discreta de risco é suposta. Você pode querer interpolar entre duas estimativas, supondo, por exemplo, que o risco permaneça constante.
Ocram 26/12/12
Método de Breslow (1974)
tomka
kidney$time >= y[l]ystatus=0status=1d=2d=1status=0
Como @tomka mencionou. Substituir a coxphchamada por fit<-coxph(Surv(time, status)~age, data=kidney, method="breslow")corrigirá a diferença de métodos.
Mr.bjerre 16/05/19