Cox-Basisrisiko

19

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:

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

Dihan
quelle
@gung könntest du bei dieser Frage helfen ? Ich habe ein paar Tage gekämpft ...
Haitao Du

Antworten:

21

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 und

H^0(t)=y(l)th^0(y(l)),
h^0(y(l))=d(l)jR(y(l))exp(xjβ)
y(1)<y(2)<d(l)y(l)R(y(l))ist das Risiko bei das alle Personen enthält, die für das Ereignis bei noch anfällig sind .y(l)y(l)

Lass uns das versuchen. (Der folgende Code dient nur zur Veranschaulichung und ist nicht als sehr gut geschrieben gedacht.)

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

Teilleistung:

       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

Ich vermute, dass der geringfügige Unterschied auf die Annäherung der Teilwahrscheinlichkeit coxph()aufgrund von Unstimmigkeiten in den Daten zurückzuführen ist ...

Ocram
quelle
Vielen Dank. Ja, es gibt geringfügige Unterschiede bei der Näherungsmethode. Es gibt jedoch 76 Zeitpunkte mit Bindungen, wenn ich die Grundgefahr für jeden Zeitpunkt ermitteln möchte. Was kann ich tun? Welche Art der Änderung im R-Code ist erforderlich?
Dihan
1
Die diskretisierte Gefahr ist null, außer zu Ereigniszeiten. Dies liefert in der Tat den größten Beitrag zur Wahrscheinlichkeit, wenn eine diskrete Gefahrenfunktion angenommen wird. Möglicherweise möchten Sie zwischen zwei Schätzungen interpolieren, wenn beispielsweise die Gefahr konstant bleibt.
26.
Methode von Breslow (1974)
Tomka
kidney$time >= y[l]ystatus=0status=1d=2d=1status=0
Wie @tomka erwähnt. Durch Ersetzen des coxphAnrufs durch fit<-coxph(Surv(time, status)~age, data=kidney, method="breslow")wird der Methodenunterschied behoben.
mr.bjerre