Disons que j'ai un ensemble de données «cathéter rénal». J'essaie de modéliser une courbe de survie en utilisant un modèle de Cox. Si je considère un modèle de Cox: j'ai besoin de l'estimation du danger de base. En utilisant la fonction intégrée du package R , je peux facilement le faire comme ceci:
survival
basehaz()
library(survival)
data(kidney)
fit <- coxph(Surv(time, status) ~ age , kidney)
basehaz(fit)
Mais si je veux écrire une fonction étape par étape du danger de base pour une estimation donnée du paramètre, b
comment puis-je procéder? J'ai essayé:
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)
Mais cela ne donne pas le même résultat que basehaz(fit)
. Quel est le problème?
Réponses:
Apparemment,
basehaz()
calcule en fait un taux de risque cumulé, plutôt que le taux de risque lui-même. La formule est la suivante: avec où désignent les temps d'événement distincts, est le nombre d'événements à et est le risque fixé à contenant tous les individus encore sensibles à l'événement à .Essayons ça. (Le code suivant est uniquement à titre d'illustration et n'est pas destiné à être très bien écrit.)
sortie partielle:
Je soupçonne que la légère différence pourrait être due à l'approximation de la probabilité partielle en
coxph()
raison de liens dans les données ...la source
kidney$time >= y[l]
status=0
status=1
status=0
coxph
appel parfit<-coxph(Surv(time, status)~age, data=kidney, method="breslow")
fixera la différence de méthodes.