Riesgo basal de Cox

Question

Supongamos que tengo un conjunto de datos de "catéter renal". Intento modelar una curva de supervivencia utilizando un modelo de Cox. Si considero un modelo de Cox: $$h(t,Z) = h_0 \exp(b'Z),$$ Necesito la estimación del riesgo base. Utilizando el survival función del paquete R basehaz() puedo hacerlo fácilmente así:

library(survival)

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

Pero si quiero escribir una función paso a paso del peligro base para una estimación dada del parámetro b ¿cómo puedo proceder? Lo he intentado:

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)

Pero esto no da el mismo resultado que basehaz(fit) . ¿Cuál es el problema?

Answer 1

Aparentemente, basehaz() en realidad calcula una tasa de peligrosidad acumulada, en lugar de la tasa de peligrosidad propiamente dicha. La fórmula es la siguiente: $$ \hat{H}_0(t) = \sum_{y_{(l)} \leq t} \hat{h}_0(y_{(l)}), $$ con $$ \hat{h}_0(y_{(l)}) = \frac{d_{(l)}}{\sum_{j \in R(y_{(l)})} \exp(\mathbf{x}^{\prime}_j \mathbf{\beta})} $$ donde $y_{(1)} < y_{(2)} < \cdots$ denotan los distintos tiempos de los sucesos, $d_{(l)}$ es el número de eventos en $y_{(l)}$ y $R(y_{(l)})$ es el riesgo fijado en $y_{(l)}$ que contiene a todos los individuos aún susceptibles de sufrir el suceso en $y_{(l)}$ .

Intentemos esto. (El siguiente código es sólo ilustrativo y no pretende estar muy bien 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))

salida 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

Sospecho que la ligera diferencia puede deberse a la aproximación de la probabilidad parcial en coxph() debido a los empates en los datos...