predict.coxph()
oblicza współczynnik ryzyka w stosunku do średniej próbki dla wszystkich zmiennych predykcyjnych . Czynniki są jak zwykle konwertowane na fałszywe predyktory, których średnią można obliczyć. Przypomnijmy, że model PH Coxa jest modelem liniowym dla log-hazardu ln h ( t ) :plnh(t)
lnh(t)=lnh0(t)+β1X1+⋯+βpXp=lnh0(t)+Xβ
Gdzie jest nieokreślonym ryzykiem odniesienia. Równoważnie, zagrożenie h ( t ) modeluje się jako h ( t ) = h 0 ( t ) ⋅ e β 1 X 1 + ⋯ + β p X p = h 0 ( t ) ⋅ e X β . Współczynnik ryzyka między dwiema osobami i i i ′ z wartościami predyktorah0(t)h(t)h ( t ) = h0( t ) ⋅ eβ1X1+ ⋯ + βpXp= godz0( t ) ⋅ eX βjaja′ i X i ' są zatem niezależne od podstawowego zagrożenia i niezależne od czasut:XjaXja′t
hi(t)hi′(t)=h0(t)⋅eXiβh0(t)⋅eXi′β=eXiβeXi′β
Dla stosunku szacunkowa zagrożenia między osobami i í ' , po prostu podłącz oszacowań współczynników b 1 , ... , b p dla p 1 , ... , β p , dając e X i b i e X i " b .ii′b1,…,bpβ1,…,βpeXibeXi′b
Jako przykład w R używam danych z dodatku Johna Foxa na modelu Cox-PH, który zapewnia bardzo ładny tekst wprowadzający. Najpierw pobieramy dane i budujemy prosty model Cox-PH dla czasu aresztowania uwolnionych więźniów ( fin
: czynnik - otrzymano pomoc finansową z kodowaniem pozorowanym "no"
-> 0, "yes"
-> 1 age
,: wiek w momencie uwolnienia, prio
: liczba wcześniejszych wyroków skazujących):
> 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
Teraz podłączamy średnie przykładowe dla naszych predyktorów do wzoru :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)
Teraz podłączamy wartości predykcyjne pierwszych 4 osób do wzoru .eXb
r1234 <- exp(coefCPH["finyes"]*(as.numeric(Rossi[1:4, "fin"])-1)
+ coefCPH["age"] *Rossi[1:4, "age"]
+ coefCPH["prio"] *Rossi[1:4, "prio"])
Teraz obliczyć względne ryzyko dla pierwszych 4 osób w stosunku do średniej próbki i porównać z wynikami 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
Jeśli masz model warstwowy, porównanie predict.coxph()
jest w stosunku do średnich warstw, można to kontrolować za pomocą reference
opcji opisanej na stronie pomocy.
meanFin <- mean(as.numeric(Rossi$fin) - 1)
nie ma większego sensu, ponieważfin
jest kategoryczny. Nie musiszmodeFin <- get_Mode(Rossi$fin)
w tym przypadku?fin
jest binarny, więc numeryczna reprezentacja współczynnika ma po prostu wartości 1 i 2. Odjęcie 1 daje nam zmienną kodowaną obojętnie o wartościach 0 i 1, która również pojawia się w macierzy obliczeniowej. Pamiętaj, że to nie zadziała w przypadku czynników z więcej niż 2 poziomami. Z pewnością można dyskutować, czy uśrednianie zmiennych zastępczych ma sens, ale właśnie topredict.coxph()
robi.