Powiedzmy, że mam zestaw danych „cewnika nerkowego”. Próbuję modelować krzywą przeżycia za pomocą modelu Coxa. Jeśli wezmę pod uwagę model Coxa: potrzebuję oszacowania podstawowego zagrożenia. Korzystając z wbudowanej funkcji pakietu R , mogę łatwo to zrobić w następujący sposób:
survival
basehaz()
library(survival)
data(kidney)
fit <- coxph(Surv(time, status) ~ age , kidney)
basehaz(fit)
Ale jeśli chcę napisać krok po kroku funkcję podstawowego zagrożenia dla danego oszacowania parametru, b
jak mogę kontynuować? Próbowałem:
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)
Ale to nie daje takiego samego rezultatu jak basehaz(fit)
. Jaki jest problem?
Odpowiedzi:
Najwyraźniej
basehaz()
faktycznie oblicza skumulowany wskaźnik zagrożenia, a nie sam wskaźnik zagrożenia. Wzór jest w godzinach : 0 ( y ( l ) ) = d ( L )Spróbujmy tego. (Poniższy kod służy wyłącznie do ilustracji i nie jest przeznaczony do bardzo dobrego pisania).
częściowa wydajność:
Podejrzewam, że niewielka różnica może wynikać z przybliżenia częściowego prawdopodobieństwa z
coxph()
powodu powiązań w danych ...źródło
kidney$time >= y[l]
może napotykać problemy numeryczne, gdy czas jest liczbowy z powodu tworzenia tabel w tabelachstatus=0
status=1
status=0
coxph
wywołaniafit<-coxph(Surv(time, status)~age, data=kidney, method="breslow")
naprawi różnicę metod.