Я использую книгу Мура «Прикладной анализ выживания с использованием R», чтобы попытаться смоделировать некоторые данные о времени до события. Проблема, с которой я сталкиваюсь, - это построение расчетных кривых выживаемости по модели Кокса. Из-за этого мне интересно, неверно ли мое понимание модели. Мои данные просты: столбец времени t , столбец индикатора события (1 для события 0 для цензора) i и столбец предиктора с 6 уровнями факторов p .
Я считаю, что могу построить расчетные кривые выживаемости для модели Кокса, как показано ниже. Но я не понимаю, как использовать Survfit и Baseplot или функции из Survminer для достижения той же цели. Вот некоторый общий код c для пояснения моего вопроса. Я буду использовать набор данных PharmcoSmoking, чтобы продемонстрировать свою проблему.
library(survival)
library(asaur)
t<-pharmacoSmoking$longestNoSmoke
i<-pharmacoSmoking$relapse
p<-pharmacoSmoking$levelSmoking
data<-as.data.frame(cbind(t,i,p))
model <- coxph(Surv(data$t, data$i) ~ p, data=data)
Насколько я понимаю, со следующими фрагментами кода, смоделированными на основе книжных примеров, базовая (совокупная) опасность на моем уровне эталонного фактора для p может быть получено из
base<-basehaz(model, centered=F)
Оценка кривой выживаемости дается как
s<-exp(-base$hazard)
t<-base$time
plot(s~t, typ = "l")
Кривая выживаемости, связанная с другой уровень фактора может быть затем задан как
beta_n<-model$coefficients #only one coef in this case
s_n <- s^(exp(beta_n))
lines(s_n~t)
, где beta_n - коэффициент для n-го уровня фактора из модели Кокса. Приведенный выше код дает то, что я считаю оценочными кривыми выживаемости для тяжелых и легких курильщиков в наборе данных PharmcoSmokers.
Так как это небольшой код, который я искал для пакетов для однострочного решения, мне было трудно с документацией по Survival (примеров в документации было немного), а также попробовал Survminer. Для последнего я пробовал:
library(survminer)
ggadjustedcurves(model, variable ="p" , data=data)
Это дает мне нечто иное, чем мой предыдущий код, хотя он похож. Метод, который я использовал ранее, неверен? Или существует другая методология, объясняющая разницу? Код выжившего не работает с моими данными (я получаю ошибку «не могу выделить вектор размера yada yada, а мои данные составляют ~ 1 млн строк»), что кажется странным, учитывая, что я могу без проблем создавать графики, используя то, что делал раньше. Это основная причина, по которой мне интересно, понимаю ли я, как построить кривые выживания для моей модели.