Wie ist die Ausgabe von predict.coxph zu interpretieren?

17

Nach dem Anpassen eines Steuermodells ist es möglich, Vorhersagen zu treffen und das relative Risiko neuer Daten abzurufen. Was ich nicht verstehe, ist, wie das relative Risiko für eine Person berechnet wird und in Bezug auf was es sich handelt (dh den Durchschnitt der Bevölkerung)? Irgendwelche Empfehlungen für Ressourcen zum besseren Verständnis (ich bin in der Überlebensanalyse nicht sehr weit fortgeschritten, je einfacher desto besser)?

user4673
quelle

Antworten:

32

predict.coxph()berechnet die Hazard Ratio relativ zum Stichprobenmittel für alle Prädiktorvariablen. Faktoren werden wie gewohnt in Dummy-Prädiktoren konvertiert, deren Durchschnitt berechnet werden kann. Recall , dass das PH Cox - Modell ein lineares Modell für das log-hazard ln h ( t ) :plnh(t)

lnh(t)=lnh0(t)+β1X1++βpXp=lnh0(t)+Xβ

Dabei ist die nicht spezifizierte Grundgefahr. In äquivalenter Weise die Gefahr h ( t ) wird als modellierten h ( t ) = h 0 ( t ) e β 1 X 1 + + β p X p = h 0 ( t ) e X β . Das Gefährdungsverhältnis zwischen zwei Personen i und i ' mit Prädiktorwertenh0(t)h(t)h(t)=h0(t)eβ1X1++βpXp=h0(t)eXβii und XXi ist somit unabhängig von der Grundgefahr und unabhängig von der Zeitt:Xit

hi(t)hi(t)=h0(t)eXiβh0(t)eXiβ=eXiβeXiβ

Für das geschätzte Gefährdungsverhältnis zwischen Personen und i ' fügen wir einfach die Koeffizientenschätzungen b 1 , , b p für die β 1 , , β p ein , wobei sich e X i b und e X i ' b ergeben .iib1,,bpβ1,,βpeXibeXib

Als Beispiel in R verwende ich die Daten aus dem Anhang von John Fox zum Cox-PH-Modell , das einen sehr schönen Einführungstext enthält. Zuerst holen wir die Daten und bauen ein einfaches Cox-PH-Modell für die Zeit bis zur Festnahme von freigelassenen Gefangenen auf ( fin: Faktor - erhaltene finanzielle Unterstützung mit Dummy-Codierung "no"-> 0, "yes"-> 1 age,: Alter zum Zeitpunkt der Freilassung, prioAnzahl früherer Verurteilungen):

> URL   <- "http://socserv.mcmaster.ca/jfox/Books/Companion/data/Rossi.txt"
> Rossi <- read.table(URL, header=TRUE)                  # our data
> Rossi[1:3, c("week", "arrest", "fin", "age", "prio")]  # looks like this
  week arrest fin age prio
1   20      1  no  27    3
2   17      1  no  18    8
3   25      1  no  19   13

> library(survival)                                      # for coxph()    
> fitCPH <- coxph(Surv(week, arrest) ~ fin + age + prio, data=Rossi)    # Cox-PH model
> (coefCPH <- coef(fitCPH))                              # estimated coefficients
     finyes         age        prio 
-0.34695446 -0.06710533  0.09689320 

eXb

meanFin  <- mean(as.numeric(Rossi$fin) - 1)   # average of financial aid dummy
    meanAge  <- mean(Rossi$age)                   # average age
meanPrio <- mean(Rossi$prio)                  # average number of prior convictions
rMean <- exp(coefCPH["finyes"]*meanFin        # e^Xb
           + coefCPH["age"]   *meanAge
           + coefCPH["prio"]  *meanPrio)

eXb

r1234 <- exp(coefCPH["finyes"]*(as.numeric(Rossi[1:4, "fin"])-1)
           + coefCPH["age"]   *Rossi[1:4, "age"]
           + coefCPH["prio"]  *Rossi[1:4, "prio"])

Berechnen Sie nun das relative Risiko für die ersten 4 Personen anhand des Stichprobenmittelwerts und vergleichen Sie es mit der Ausgabe von predict.coxph().

> r1234 / rMean
[1] 1.0139038 3.0108488 4.5703176 0.7722002

> relRisk <- predict(fitCPH, Rossi, type="risk")   # relative risk
> relRisk[1:4]
        1         2         3         4 
1.0139038 3.0108488 4.5703176 0.7722002

Wenn Sie ein geschichtetes Modell haben, erfolgt der Vergleich mit predict.coxph()den Schichtendurchschnitten. Dies kann über die referenceOption gesteuert werden, die auf der Hilfeseite erläutert wird.

caracal
quelle
2
+1 weil es nicht offensichtlich ist, was predict.coxph genau auf der Hilfeseite macht!
2.
das war großartig! Sehr einfach zu verstehen!
user4673
meanFin <- mean(as.numeric(Rossi$fin) - 1)macht nicht viel sinn, da finkategorisch ist. Müssen Sie modeFin <- get_Mode(Rossi$fin)in diesem Fall nicht?
Zhubarb
1
@Zhubarb finist binär, daher hat die numerische Darstellung des Faktors nur die Werte 1 und 2. Durch Subtrahieren von 1 erhalten wir die Dummy-codierte Variable mit den Werten 0 und 1, die auch in der Entwurfsmatrix erscheint. Beachten Sie, dass dies für Faktoren mit mehr als 2 Ebenen nicht funktioniert. Es ist sicherlich fraglich, ob die Mittelung von Dummy-Variablen sinnvoll ist, aber genau das predict.coxph()tut es.
Caracal
Wie würden Sie in Worten eine Hazard Ratio von 3,01 interpretieren (zB relRisk [2])?
RNB