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:Xi′t
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 .ii′b1,…,bpβ1,…,βpeXibeXi′b
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, prio
Anzahl 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 reference
Option gesteuert werden, die auf der Hilfeseite erläutert wird.
meanFin <- mean(as.numeric(Rossi$fin) - 1)
macht nicht viel sinn, dafin
kategorisch ist. Müssen SiemodeFin <- get_Mode(Rossi$fin)
in diesem Fall nicht?fin
ist 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 daspredict.coxph()
tut es.