Zagrożenie podstawowe Coxa

20

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:

h(t,Z)=h0exp(bZ),
survivalbasehaz()
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, bjak 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?

Dihan
źródło
@ gung, czy możesz pomóc z tym pytaniem ? Walczyłem przez kilka dni ...
Haitao Du

Odpowiedzi:

22

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 )

H.^0(t)=y(l)th^0(y(l)),
gdziey(1)<y(2)<oznaczają różne czasy zdarzeń,d(l)jest liczbą zdarzeń wy(l), aR(y(l))oznacza ryzyko ustalone nay(l)
h^0(y(l))=re(l)jotR(y(l))exp(xjotβ)
y(1)<y(2))<re(l)y(l)R(y(l))y(l)zawierający wszystkie osoby wciąż podatne na zdarzenie w punkcie .y(l)

Spróbujmy tego. (Poniższy kod służy wyłącznie do ilustracji i nie jest przeznaczony do bardzo dobrego pisania).

#------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))

częściowa wydajność:

       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

Podejrzewam, że niewielka różnica może wynikać z przybliżenia częściowego prawdopodobieństwa z coxph()powodu powiązań w danych ...

ocram
źródło
Wielkie dzięki. Tak, istnieje niewielka różnica w metodzie aproksymacji. Ale jest 76 punktów czasowych z remisami, jeśli chcę znaleźć podstawowe zagrożenie dla każdego punktu czasowego. Co mogę zrobić? Jakiego rodzaju modyfikacja kodu R jest potrzebna?
Dihan
1
Zagrożenie dyskrecjonalne wynosi zero, z wyjątkiem czasu zdarzenia. To rzeczywiście ma największy udział w prawdopodobieństwie, jeśli zakłada się dyskretną funkcję hazardu. Możesz interpolować dowolne dwa oszacowania, zakładając na przykład, że zagrożenie pozostaje stałe.
ocram
Method of Breslow (1974)
tomka
Muszę zauważyć pewne problemy z tym wdrożeniem. Używanie kidney$time >= y[l]może napotykać problemy numeryczne, gdy czas jest liczbowy z powodu tworzenia tabel w tabelachystatus=0status=1re=2)re=1status=0
Jak wspomniano @tomka. Zastąpienie coxphwywołania fit<-coxph(Surv(time, status)~age, data=kidney, method="breslow")naprawi różnicę metod.
mr.bjerre