Как построить кривую выживаемости, вызванную выживанием (пакетная выживаемость R)?

Im пытается подогнать и построить модель Вейбулла для данных о выживании. Данные имеют только одну ковариацию, когорту, которая работает с 2006 по 2010 год. Итак, какие-либо идеи о том, что добавить к двум строкам кода, которые следует, чтобы построить кривую выживаемости когорты 2010 года?

library(survival)
s <- Surv(subSetCdm$dur,subSetCdm$event)
sWei <- survreg(s ~ cohort,dist='weibull',data=subSetCdm)

Выполнение этого же с моделью Cox PH довольно просто, со следующими строками. Проблема в том, что survif() не принимает объекты типа revreg.

sCox <- coxph(s ~ cohort,data=subSetCdm)
cohort <- factor(c(2010),levels=2006:2010)
sfCox <- survfit(sCox,newdata=data.frame(cohort))
plot(sfCox,col='green')

Используя данные легкие (из пакета выживания), вот что я пытаюсь выполнить.

#create a Surv object
s <- with(lung,Surv(time,status))

#plot kaplan-meier estimate, per sex
fKM <- survfit(s ~ sex,data=lung)
plot(fKM)

#plot Cox PH survival curves, per sex
sCox <- coxph(s ~ as.factor(sex),data=lung)
lines(survfit(sCox,newdata=data.frame(sex=1)),col='green')
lines(survfit(sCox,newdata=data.frame(sex=2)),col='green')

#plot weibull survival curves, per sex, DOES NOT RUN
sWei <- survreg(s ~ as.factor(sex),dist='weibull',data=lung)
lines(survfit(sWei,newdata=data.frame(sex=1)),col='red')
lines(survfit(sWei,newdata=data.frame(sex=2)),col='red')

Ответ 1

Надеюсь, это поможет, и я не сделал какую-то вводящую в заблуждение ошибку:

скопировано сверху:

    #create a Surv object
    s <- with(lung,Surv(time,status))

    #plot kaplan-meier estimate, per sex
    fKM <- survfit(s ~ sex,data=lung)
    plot(fKM)

    #plot Cox PH survival curves, per sex
    sCox <- coxph(s ~ as.factor(sex),data=lung)
    lines(survfit(sCox,newdata=data.frame(sex=1)),col='green')
    lines(survfit(sCox,newdata=data.frame(sex=2)),col='green')

для Weibull, используйте предсказать, комментарий от Vincent:

    #plot weibull survival curves, per sex,
    sWei <- survreg(s ~ as.factor(sex),dist='weibull',data=lung)

    lines(predict(sWei, newdata=list(sex=1),type="quantile",p=seq(.01,.99,by=.01)),seq(.99,.01,by=-.01),col="red")
    lines(predict(sWei, newdata=list(sex=2),type="quantile",p=seq(.01,.99,by=.01)),seq(.99,.01,by=-.01),col="red")

Трюк здесь заключался в том, чтобы изменить порядок квантилей для построения прогноза. Вероятно, лучший способ сделать это, но он работает здесь. Удачи!

Ответ 2

Альтернативным вариантом является использование пакета flexsurv. Это предлагает некоторые дополнительные функции над пакетом survival - в том числе, что функция параметрической регрессии flexsurvreg() имеет хороший метод построения, который делает то, что вы просите.

Использование легких, как указано выше;

#create a Surv object
s <- with(lung,Surv(time,status))

require(flexsurv)
sWei  <- flexsurvreg(s ~ as.factor(sex),dist='weibull',data=lung)
sLno  <- flexsurvreg(s ~ as.factor(sex),dist='lnorm',data=lung)   

plot(sWei)
lines(sLno, col="blue")

Вы можете построить график совокупной опасности или опасности с помощью аргумента type и добавить доверительные интервалы с аргументом ci.

Ответ 3

Это всего лишь примечание, поясняющее ответ Тима Риффа, в котором используется следующий код:

lines(predict(sWei, newdata=list(sex=1),type="quantile",p=seq(.01,.99,by=.01)),seq(.99,.01,by=-.01),col="red")
lines(predict(sWei, newdata=list(sex=2),type="quantile",p=seq(.01,.99,by=.01)),seq(.99,.01,by=-.01),col="red")

Причина для двух последовательностей зеркального изображения seq(.01,.99,by=.01) и seq(.99,.01,by=-.01) заключается в том, что метод pred() дает квантили для распределения событий f (t) - то есть значения обратного CDF f (t) - в то время как кривая выживаемости рисует 1- (CDF по f) по сравнению с t. Другими словами, если вы построите p по сравнению с прогнозом (p), вы получите CDF, и если вы построите 1-p против прогноза (p), вы получите кривую выживаемости, которая равна 1-CDF. Следующий код более прозрачен и обобщает на произвольные векторы значений p:

pct <- seq(.01,.99,by=.01)
lines(predict(sWei, newdata=list(sex=1),type="quantile",p=pct),1-pct,col="red")
lines(predict(sWei, newdata=list(sex=2),type="quantile",p=pct),1-pct,col="red")