使用R中的flexsurvreg预测样本外数据

Question

我在R语言中有以下模型

library(flexsurv)

data(ovarian)

model = flexsurvreg(Surv(futime, fustat) ~ ecog.ps + rx, data = ovarian, dist='weibull')

model

predict(model,data = ovarian, type = 'response')

模型概述如下：flexsurvreg模型输出 我正在尝试使用R中的预测函数来预测生存时间，但是遇到了以下错误：在尝试预测时出现错误 如何使用这个flexsurvreg模型来预测预期寿命？

我知道文档中提到了一个totlos.fs函数，但是这个数据似乎没有totlos.fs需要提供输出的trans变量。

如果没有其他替代方法，如何在这个数据中创建一个trans变量，并将其与现有协变量一起处理？

请给予建议。

- Nik

Answer 1

补充示例文档的第3节介绍了一个示例，其中使用模型方程直接计算预测值。由于您正在使用Weibull分布（具有n=2个参数），我相信这应该是可行的：

pred.model <- model.matrix(model) %*% model$res[-(1:n),"est"]

干杯

- Kivevi

---

谢谢，我会尝试一下！ - Nik

Answer 2

Nik，

我知道你的问题是一个老问题，但是请看下面我是如何通过一种方法来解决它的。它涉及到从测试数据的拟合中检索形状和速率参数，然后不使用predict，而是使用flexsurv中的qgompertz()函数。请原谅我使用了自己封装的示例代码，但你应该能够跟上。

# generate the training data "lung1" from data(lung) in survival package
# hacked way for truncating the lung data to 2 years of follow up
require(survival)
lung$yrs <- lung$time/365
lung1 <- lung[c("status", "yrs")]
lung1$status[ lung1$yrs >2] <- 1
lung1$yrs[ lung1$yrs >2]  <- 2

# from the training data build KM to obtain survival %s
s <- Surv(time=lung1$yrs, event=lung1$status)
km.lung <- survfit(s ~ 1, data=lung1)
plot(km.lung)

# generate dataframe to use later for plotting 
cut.length <- sum((km.lung$time <= 2)) # so I can create example test data
test.data <- data.frame(yrs = km.lung$time[1:cut.length] , surv=round(km.lung$surv[1:cut.length], 3))


##
##  doing the same as above with gompertz
##
require(flexsurv) #needed to run gompertz model
s <- Surv(time=lung1$yrs, event=lung1$status)
gomp <- flexsurvreg(s ~ 1, data=lung1, dist="gompertz") # run this to get shape and rate estimates for gompertz
gomp # notice the shape and rate values 

# create variables for these values
g.shape <- 0.5866
g.rate <- 0.5816


##
##  plot data and vizualize the gomperts
##
# vars for plotting
df1 <- test.data 
xvar <- "yrs"
yvar <- "surv"

extendedtime <- 3 # 
ylim1 <- c(0,1)
xlim1 <- c(0, extendedtime)

# plot the survival % for training data
plot(df1[,yvar]~df1[,xvar], type="S", ylab="", xlab="", lwd=3, xlim=xlim1, ylim=ylim1)
# Nik--here is where the magic happens... pay special attention to: qgompertz(seq(.01,.99,by=.01), shape=0.58656, rate = .5816) 
lines (qgompertz(seq(.01,.99,by=.01), shape=0.58656, rate = .5816) ,  seq(.99,.01,by=-.01) , col="red", lwd=2, lty=2  )

# generate a km curve from the testing data
s <- Surv(time=lung$yrs, event=lung$status)
km.lung <- survfit(s ~ 1, data=lung)
par(new=T)
# now draw remaining survival curve from the testing section
plot(km.lung$surv[(cut.length+1):length(km.lung$time)]~km.lung$time[(cut.length+1):length(km.lung$time)], type="S", col="blue", ylab="", xlab="", lwd=3, xlim=xlim1, ylim=ylim1)

- Chris