如何从 R 中的 survfit 中提取危害？

Question

我有一个survfit 对象。对我的t=0:50 年感兴趣的摘要生存很容易。

summary(survfit, t=0:50)

它给出了每个 t 的生存率。

有没有办法获得每个 t 的危险（在这种情况下，每个 t=0:50 中从 t-1 到 t 的危险）？我想获得与 Kaplan Meier 曲线相关的风险的平均值和置信区间（或标准误差）。

当分布适合时，这似乎很容易做到（例如flexsurvreg 中的type="hazard"），但我不知道如何为常规的 survfit 对象执行此操作。建议？

Answer 1

这有点棘手，因为风险是对瞬时概率的估计（这是离散数据），但basehaz 函数可能会有所帮助，但它只返回累积风险。所以你仍然需要执行一个额外的步骤。

muhaz 函数我也很幸运。从其文档中：

library(muhaz)
?muhaz
data(ovarian, package="survival")
attach(ovarian)
fit1 <- muhaz(futime, fustat)
plot(fit1)

我不确定获得 95% 置信区间的最佳方法，但自举可能是一种方法。

#Function to bootstrap hazard estimates
haz.bootstrap <- function(data,trial,min.time,max.time){
  library(data.table)
  data <- as.data.table(data)
  data <- data[sample(1:nrow(data),nrow(data),replace=T)]
  fit1 <- muhaz(data$futime, data$fustat,min.time=min.time,max.time=max.time)
  result <- data.table(est.grid=fit1$est.grid,trial,haz.est=fit1$haz.est)
  return(result)
}

#Re-run function to get 1000 estimates
haz.list <- lapply(1:1000,function(x) haz.bootstrap(data=ovarian,trial=x,min.time=0,max.time=744))
haz.table <- rbindlist(haz.list,fill=T)

#Calculate Mean,SD,upper and lower 95% confidence bands
plot.table <- haz.table[, .(Mean=mean(haz.est),SD=sd(haz.est)), by=est.grid]
plot.table[, u95 := Mean+1.96*SD]
plot.table[, l95 := Mean-1.96*SD]

#Plot graph
library(ggplot2)
p <- ggplot(data=plot.table)+geom_smooth(aes(x=est.grid,y=Mean))
p <- p+geom_smooth(aes(x=est.grid,y=u95),linetype="dashed")
p <- p+geom_smooth(aes(x=est.grid,y=l95),linetype="dashed")
p

Answer 2

作为对 Mike 的回答的补充，可以通过泊松分布而不是正态分布来模拟事件的数量。然后可以通过伽马分布计算危险率。代码将变为：

library(muhaz)
library(data.table)
library(rGammaGamma)
data(ovarian, package="survival")
attach(ovarian)
fit1 <- muhaz(futime, fustat)
plot(fit1)

#Function to bootstrap hazard estimates
haz.bootstrap <- function(data,trial,min.time,max.time){
  library(data.table)
  data <- as.data.table(data)
  data <- data[sample(1:nrow(data),nrow(data),replace=T)]
  fit1 <- muhaz(data$futime, data$fustat,min.time=min.time,max.time=max.time)
  result <- data.table(est.grid=fit1$est.grid,trial,haz.est=fit1$haz.est)
  return(result)
}

#Re-run function to get 1000 estimates
haz.list <- lapply(1:1000,function(x) haz.bootstrap(data=ovarian,trial=x,min.time=0,max.time=744))
haz.table <- rbindlist(haz.list,fill=T)

#Calculate Mean, gamma parameters, upper and lower 95% confidence bands
plot.table <- haz.table[, .(Mean=mean(haz.est),
                            Shape = gammaMME(haz.est)["shape"],
                            Scale = gammaMME(haz.est)["scale"]), by=est.grid]
plot.table[, u95 := qgamma(0.95,shape = Shape + 1, scale = Scale)]
# The + 1 is due to the discrete character of the poisson distribution.
plot.table[, l95 := qgamma(0.05,shape = Shape, scale = Scale)]

#Plot graph
ggplot(data=plot.table) + 
  geom_line(aes(x=est.grid, y=Mean),col="blue") + 
  geom_ribbon(aes(x=est.grid, y=Mean, ymin=l95, ymax=u95),alpha=0.5, fill= "lightblue")

可以看出，对危险率下限的负估计现在已经消失了。

Answer 3

为了性能，我们可以使用更精简的引导函数。

## define custom times
t0 <- 0
t1 <- 744

## bootstrap fun
boot_fun <- function(x) {
  n <- dim(x)[1]
  x <- x[sample.int(n, n, replace=TRUE), ]
  muhaz::muhaz(x$futime, x$fustat, min.time=t0, max.time=t1)
}

# bootstrap
set.seed(42)
R <- 999
B <- replicate(R, boot_fun(ovarian))

结果可以手动计算。

## extract matrix from bootstrap
r <- `colnames<-`(t(array(unlist(B[3, ]), dim=c(101, R))), B[2, ][[1]])

## calculate result
library(matrixStats)  ## for fast matrix calculations
r <- cbind(x=as.numeric(colnames(r)), 
           y=colMeans2(r),
           shape=(colMeans2(r)/colSds(r))^2, 
           scale=colVars(r)/colMeans2(r))
r <- cbind(r[, 1:2], 
           lower=qgamma(0.025, shape=r[, 'shape'] + 1, scale=r[, 'scale']),
           upper=qgamma(0.975, shape=r[, 'shape'], scale=r[, 'scale']))

给予

head(r)
#          x            y        lower       upper
# [1,]  0.00 0.0003836816 9.359267e-05 0.001400539
# [2,]  7.44 0.0003913992 9.746868e-05 0.001387551
# [3,] 14.88 0.0003997275 1.018897e-04 0.001374005
# [4,] 22.32 0.0004087439 1.069353e-04 0.001360212
# [5,] 29.76 0.0004178464 1.123697e-04 0.001346187
# [6,] 37.20 0.0004275531 1.184685e-04 0.001332237

range(r[, 'y'])
# [1] 0.0003836816 0.0011122910

情节

matplot(r[, 1], r[, -1], type='l', lty=c(1, 2, 2), col=4, 
        xlab='Time', ylab='Hazard Rate', main='Hazard Estimates')
legend('topleft', legend=c('estimate', '95% CI'), col=4, lty=1:2, cex=.8)

---

数据

data(cancer, package="survival")  ## loads `ovarian` data set