Участок Каплан-Мейер для регрессии Кокса

Question

У меня есть модель пропорциональных рисков Кокса, созданная с использованием следующего кода в R, который прогнозирует смертность. Ковариат A, B и C добавляются просто, чтобы избежать смешения (например, возраста, пола, расы), но нас действительно интересует предсказатель X. X - непрерывная переменная.

cox.model <- coxph(Surv(time, dead) ~ A + B + C + X, data = df) 

Теперь у меня возникают проблемы с графикой кривой Каплана-Мейера. Я искал, как создать эту фигуру, но мне не повезло. Я не уверен, возможно ли создание Каплана-Мейера для модели Кокса? Каплан-Мейер корректирует мои ковариаты или не нуждается в них?

То, что я пробовал, приведен ниже, но мне сказали, что это неправильно.

plot(survfit(cox.model), xlab = 'Time (years)', ylab = 'Survival Probabilities') 

Я также попытался построить рисунок, показывающий совокупную опасность смертности. Я не знаю, делаю ли я это правильно, так как я пробовал несколько разных способов и получал разные результаты. В идеале я хотел бы построить две строки, которые показывают риск смертности для 75-го процентиля X и один, который показывает 25-й процентиль X. Как я могу это сделать?

Я мог бы перечислить все, что я пробовал, но я никого не хочу смущать!

Большое спасибо.

Answer 1

Вот пример, взятый из this paper.

url <- "http://socserv.mcmaster.ca/jfox/Books/Companion/data/Rossi.txt" 
Rossi <- read.table(url, header=TRUE) 
Rossi[1:5, 1:10] 

# week arrest fin age race wexp   mar paro prio educ 
# 1 20  1 no 27 black no not married yes 3 3 
# 2 17  1 no 18 black no not married yes 8 4 
# 3 25  1 no 19 other yes not married yes 13 3 
# 4 52  0 yes 23 black yes  married yes 1 5 
# 5 52  0 no 19 other yes not married yes 3 3 

mod.allison <- coxph(Surv(week, arrest) ~ 
         fin + age + race + wexp + mar + paro + prio, 
         data=Rossi) 
mod.allison 

# Call: 
# coxph(formula = Surv(week, arrest) ~ fin + age + race + wexp + 
# mar + paro + prio, data = Rossi) 
# 
# 
#     coef exp(coef) se(coef)  z  p 
# finyes   -0.3794  0.684 0.1914 -1.983 0.0470 
# age   -0.0574  0.944 0.0220 -2.611 0.0090 
# raceother  -0.3139  0.731 0.3080 -1.019 0.3100 
# wexpyes  -0.1498  0.861 0.2122 -0.706 0.4800 
# marnot married 0.4337  1.543 0.3819 1.136 0.2600 
# paroyes  -0.0849  0.919 0.1958 -0.434 0.6600 
# prio   0.0915  1.096 0.0286 3.194 0.0014 
# 
# Likelihood ratio test=33.3 on 7 df, p=2.36e-05 n= 432, number of events= 114  

Обратите внимание, что модель использует fin, age, race, wexp, mar, paro, prio предсказать arrest. Как упоминалось в this document, функция survfit() использует оценку Каплана-Мейера для коэффициента выживаемости.

plot(survfit(mod.allison), ylim=c(0.7, 1), xlab="Weeks", 
    ylab="Proportion Not Rearrested") 

Мы получаем участок (с 95% доверительным интервалом) для выживаемости. Для кумулятивной нормы опасности вы можете сделать

# plot(survfit(mod.allison)$cumhaz) 

но это не дает доверительных интервалов. Однако, не беспокойтесь! Мы знаем, что H (t) = -ln (S (t)), и мы имеем доверительные интервалы для S (t). Все, что нам нужно сделать, это

sfit <- survfit(mod.allison) 
cumhaz.upper <- -log(sfit$upper) 
cumhaz.lower <- -log(sfit$lower) 
cumhaz <- sfit$cumhaz # same as -log(sfit$surv) 

Тогда просто построить эти

plot(cumhaz, xlab="weeks ahead", ylab="cumulative hazard", 
    ylim=c(min(cumhaz.lower), max(cumhaz.upper))) 
lines(cumhaz.lower) 
lines(cumhaz.upper) 

Вы хотите использовать survfit(..., conf.int=0.50), чтобы получить полосы на 75% и 25% вместо 97,5% и 2,5%.

Answer 2

Запрос на оценочную кривую выживаемости на 25-м и 75-м процентилях для X сначала требует определения этих процентилей и определения значений для всех остальных ковариаций в фрейме данных, который будет использоваться в качестве аргумента newdata для сохранения.:

Можно использовать данные, предложенные другим resondent с сайта Фокса, хотя на моей машине это требовало построения url -объекта:

url <- url("http://socserv.mcmaster.ca/jfox/Books/Companion/data/Rossi.txt") 
Rossi <- read.table(url, header=TRUE) 

Это, вероятно, не лучший пример для этого wquestion, но у него есть числовая переменная, которую мы можем вычислить квартили:

> summary(Rossi$prio) 
    Min. 1st Qu. Median Mean 3rd Qu. Max. 
    0.000 1.000 2.000 2.984 4.000 18.000 

Так это будет модель покрой и survfit вызовы:

mod.allison <- coxph(Surv(week, arrest) ~ 
         fin + age + race + prio , 
         data=Rossi) 
prio.fit <- survfit(mod.allison, 
        newdata= data.frame(fin="yes", age=30, race="black", prio=c(1,4))) 
plot(prio.fit, col=c("red","blue"))