Cómo estimar la función de riesgo basal en el modelo de Cox con R

Question

Necesito estimar la función de riesgo de base $\lambda_0(t)$ en un modelo de Cox dependiente del tiempo

$\lambda(t) = \lambda_0(t) \exp(Z(t)'\beta)$

Cuando hice el curso de Supervivencia, recuerdo que la derivada directa de la función de riesgo acumulativa ( $\lambda_0(t) dt = d\Lambda_0(t)$ ) no sería un buen estimador porque el estimador de Breslow da una función escalonada.

Entonces, ¿hay alguna función en R que pueda utilizar directamente? ¿O alguna referencia sobre este tema?

No estoy seguro de si vale la pena abrir otra pregunta, así que sólo añado algunos antecedentes de por qué la función de riesgo de referencia es importante para mí. La fórmula siguiente estima la probabilidad de que el tiempo de supervivencia de un sujeto sea mayor que el de otro. En el marco de un modelo de Cox, la función de riesgo inicial $\lambda_0(t)$ es necesario.

$P(T_1 > T_2 ) = - \int_0^\infty S_1(t) dS_2(t) = - \int_0^\infty S_1(t)S_2(t)\lambda_2(t)dt $

Answer 1

Se diseñó explícitamente un modelo de Cox para poder estimar las razones de riesgo sin tener que estimar la función de riesgo de base. Esto es un punto fuerte y un punto débil. El punto fuerte es que no se pueden cometer errores en las funciones que no se estiman. Esta es una verdadera fortaleza y es la razón por la que la gente se refiere a ella como "semiparamétrica" y es en gran medida responsable de su popularidad. Sin embargo, también es una debilidad real, en el sentido de que una vez que se quiere saber algo más que la razón de riesgo, a menudo se necesitará la función de riesgo de referencia y eso anula el propio propósito de un modelo de Cox.

Así que suelo utilizar los modelos de Cox sólo cuando me interesan los cocientes de riesgo y nada más. Si quiero saber otras cosas, suelo pasar a otros modelos como los que se comentan aquí: http://www.stata.com/bookstore/flexible-parametric-survival-analysis-stata/

Answer 2

La función de peligro de la línea de base se puede estimar en R utilizando la función "basehaz". El archivo de "ayuda" afirma que se trata de la función de "supervivencia predicha", lo que claramente no es así. Si uno inspecciona el código, es claramente la función de peligro acumulativo de un survfit objeto. Para más tonterías, la configuración por defecto es centered=TRUE que a) no es una función de riesgo de base (como su nombre indica), y b) emplea la predicción en la media, que está ampliamente desacreditada como válida en cualquier sentido práctico.

Y en cuanto a tu punto anterior: sí, esta función hace uso de la función escalonada. Puedes transformar esa salida en una función de riesgo utilizando el suavizado. Lo peor de todo, ¿cuál es el intervalo de incertidumbre para esa predicción? Puedes conseguir una medalla Fields si puedes derivarlo. Creo que ni siquiera sabemos si el bootstrapping funciona o no.

Como ejemplo:

set.seed(1234)
x <- rweibull(1000, 2, 3)
coxfit <- coxph(Surv(x) ~ 1)
bhest <- basehaz(coxfit)
haz <- exp(diff(bhest[, 1])*diff(bhest[, 2]))
time <- (bhest[-1,2] + bhest[-1000, 2])/2
b <- 2^-3

curve(3*b*x, from=0, to=max(x), xlab='Survival time', ylab='Weibull hazard')
points(t <- bhest[-1,2], h <- diff(bhest[, 1])/diff(bhest[, 2]), col='grey')
smooth <- loess.smooth(t, h)
lines(smooth$x, smooth$y, col='red')
legend('topright', lty=c(1,1,0), col=c('black', 'red', 'grey'), pch=c(NA,NA,1), c('Actual hazard fun', 'Smoothed hazard fun', 'Stepped discrete-time hazards'), bg='white')