Peligro de referencia de Cox

Digamos que tengo un conjunto de datos de "catéter renal". Estoy tratando de modelar una curva de supervivencia usando un modelo de Cox. Si considero un modelo de Cox: necesito la estimación del riesgo de referencia. Al usar la función integrada del paquete R , puedo hacerlo fácilmente de esta manera:

$$h(t,Z) = h_0 \exp(b'Z),$$ survivalbasehaz()

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 de la línea de base para una estimación dada del parámetro, b¿cómo puedo proceder? Lo intenté:

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?

Aparentemente, en basehaz()realidad calcula una tasa de riesgo acumulativo, en lugar de la tasa de riesgo en sí. La fórmula es la con h 0 ( y ( l ) ) = d ( l

$$\hat{H}_0(t) = \sum_{y_{(l)} \leq t} \hat{h}_0(y_{(l)}),$$ dondey(1)<y(2)<⋯denotan los tiempos de eventos distintos,d(l)es el número de eventos eny(l), yR(y(l))es el riesgo establecido eny(l $$\hat{h}_0(y_{(l)}) = \frac{d_{(l)}}{\sum_{j \in R(y_{(l)})} \exp(\mathbf{x}^{\prime}_j \mathbf{\beta})}$$ $y_{(1)} < y_{(2)} < \cdots$$d_{(l)}$$y_{(l)}$$R(y_{(l)})$$y_{(l)}$que contiene todos los individuos aún susceptibles al evento en .$y_{(l)}$

Intentemos esto. (El siguiente código está ahí solo a modo ilustrativo y no está destinado a 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 podría deberse a la aproximación de la probabilidad parcial coxph()debido a lazos en los datos ...