predict.coxph()
вычисляет соотношение риска относительно среднего образца для всех p предикторов. Факторы обычно преобразуются в фиктивные предикторы, среднее значение которых можно рассчитать. Напомним, что модель Кокса РН является линейной моделью для логарифмической опасности lnh(t) :
lnh(t)=lnh0(t)+β1X1+⋯+βpXp=lnh0(t)+Xβ
Где - неопределенная базовая опасность. Эквивалентно, опасность ч ( т ) моделируется как ч ( т ) = ч 0 ( т ) ⋅ е β 1 X 1 + ⋯ + β р Х р = ч 0 ( т ) ⋅ е Х β . Соотношение рисков между двумя людьми i и i ′ с предикторамиh0(t)h(t)h(t)=h0(t)⋅eβ1X1+⋯+βpXp=h0(t)⋅eXβii′ иXiТаким образом, X i ' не зависят от базовой опасности и не зависят от времениt:Xi′t
hi(t)hi′(t)=h0(t)⋅eXiβh0(t)⋅eXi′β=eXiβeXi′β
Для оценки соотношения рисков между людьми и i ′ мы просто включаем оценки коэффициентов b 1 , … , b p для β 1 , … , β p , давая e X i b и e X i ′ b .ii′b1,…,bpβ1,…,βpeXibeXi′b
В качестве примера на R я использую данные из приложения Джона Фокса по модели Кокса-PH, которые дают очень хороший вводный текст. Сначала мы извлекаем данные и строим простую модель Кокса-PH для времени ареста освобожденных заключенных ( fin
: фактор - полученная финансовая помощь с фиктивным кодированием "no"
-> 0, "yes"
-> 1 age
,: возраст на момент освобождения, prio
: количество предыдущих приговоров):
> 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)
Теперь мы включаем значения предикторов первых 4 человек в формулу .eXb
r1234 <- exp(coefCPH["finyes"]*(as.numeric(Rossi[1:4, "fin"])-1)
+ coefCPH["age"] *Rossi[1:4, "age"]
+ coefCPH["prio"] *Rossi[1:4, "prio"])
Теперь вычислите относительный риск для первых 4 человек относительно среднего значения по выборке и сравните с выходом из 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
Если у вас стратифицированная модель, сравнение по сравнению predict.coxph()
со средними по стратам, это можно контролировать с помощью reference
опции, которая описана на странице справки.
meanFin <- mean(as.numeric(Rossi$fin) - 1)
не имеет большого смысла, так какfin
является категоричным. Вам не нужноmodeFin <- get_Mode(Rossi$fin)
в этом случае?fin
является двоичным, поэтому числовое представление фактора имеет только значения 1 и 2. Вычитание 1 дает нам фиктивную переменную со значениями 0 и 1, которая также появляется в матрице проекта. Обратите внимание, что это не будет работать для факторов с более чем 2 уровнями. Это, конечно , спорно , делает ли усреднение фиктивного смысла, но это то, чтоpredict.coxph()
делает.