Построение простого сурбрега

Question

Это более простая вариация вопроса, на который был дан ответ на Как построить кривую выживания, генерируемую с помощью суррег (выживание пакета R)?

# Create simple Weibull survival fit using library(survival)
  surmo<-survreg( Surv(validtimes, status)~1, dist="weibull")

# Getting Kaplan-Meier
  fKM<-survfit( Surv(validtimes, status)~1)

# Plot Kaplan-Meier
  plot(fKM,xlab="Time,Days",conf.int=TRUE,mark.time=TRUE,ylab="Fraction",main="Kaplan-Meier Plot")

Все доэто работало нормально без каких-либо проблем.

Проблемы возникли, когда я захотел наложить прогнозируемое соответствие Вейбулла на данные.Основываясь на примере, который я использовал.

 pct <- seq(.01,.99,by=.01)
 maxvalidtimes<-max(validtimes)
# Getting the Weibull lines to overlay
  lines(predict(surmo,newdata=list(1:maxvalidtimes),type="quantile",p=pct),1-pct,col="red")

Я получаю ошибку

Ошибка в xy.coords (x, y): длины 'x' и 'y' различаются

Я предположил, что проблема связана с термином: newdata = list (1: maxvalidtimes)

Я попытался удалить термин newdata, а также установил newdata = list (1:99) вбезрезультатно.

Я попробовал то же самое в пакете flexsurv, и я получил точные графики, которые хотел, без особых усилий.

 # Using flexsurv package here
  surmof  <- flexsurvreg( Surv(validtimes, status)~1,dist='weibull')
  plot(surmof,mark.time=TRUE,xlab="Time,Days",ylab="Fraction",main="FlexSurv Plot")

Answer 1

Поскольку вы не предлагали никаких данных, я собираюсь изменить последний пример на странице ?predict.survreg, который использует набор данных lung.Вам не нужны никакие новые данные, так как вам нужен только график квантильного типа, для которого требуется векторный аргумент, заданный p.

lfit <- survreg(Surv(time, status) ~ 1, data=lung)
pct <- 1:98/100   # The 100th percentile of predicted survival is at +infinity
ptime <- predict(lfit,  type='quantile',
                  p=pct, se=TRUE)
 str(ptime)
#------------
List of 2
 $ fit   : num [1:228, 1:98] 12.7 12.7 12.7 12.7 12.7 ...
 $ se.fit: num [1:228, 1:98] 2.89 2.89 2.89 2.89 2.89 ...

Таким образом, у вас на самом деле слишком много точек данных, и если вы посмотритев 228 строках данных в ptime вы обнаружите, что все строки одинаковы, поэтому просто используйте первый ряд.

identical( ptime$fit[1,], ptime$fit[2,])
#[1] TRUE


 str(ptime$fit[1,])
# num [1:98] 12.7 21.6 29.5 36.8 43.8 ...

Таким образом, у вас есть прогнозируемое время для каждого квантиля и помните, что выживаниефункция равна всего 1 минус функция квантиля и что значения y являются заданными квентилями, а ее времена, которые формируют значения x:

 plot(x=ptime$fit[1,], y=1-pct, type="l")