Angenommen, ich habe einen "Nierenkatheter" -Datensatz. Ich versuche, eine Überlebenskurve mit einem Cox-Modell zu modellieren. Wenn ich ein Cox-Modell betrachte: brauche ich die Schätzung der . Mit der eingebauten Paket-R-Funktion kann ich das ganz einfach so machen:
survival
basehaz()
library(survival)
data(kidney)
fit <- coxph(Surv(time, status) ~ age , kidney)
basehaz(fit)
Aber wenn ich für eine gegebene Parameterschätzung eine Schritt-für-Schritt-Funktion der Grundgefahr schreiben möchte, b
wie kann ich dann vorgehen? Ich habe es versucht:
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)
Dies ergibt jedoch nicht das gleiche Ergebnis wie basehaz(fit)
. Was ist das Problem?
Antworten:
Anscheinend wird
basehaz()
tatsächlich eine kumulative Gefährdungsrate berechnet und nicht die Gefährdungsrate selbst. Die Formel ist wie mit wobei bezeichnen die verschiedenen Ereigniszeiten, ist die Anzahl der Ereignisse bei undLass uns das versuchen. (Der folgende Code dient nur zur Veranschaulichung und ist nicht als sehr gut geschrieben gedacht.)
Teilleistung:
Ich vermute, dass der geringfügige Unterschied auf die Annäherung der Teilwahrscheinlichkeit
coxph()
aufgrund von Unstimmigkeiten in den Daten zurückzuführen ist ...quelle
kidney$time >= y[l]
status=0
status=1
status=0
coxph
Anrufs durchfit<-coxph(Surv(time, status)~age, data=kidney, method="breslow")
wird der Methodenunterschied behoben.