Come si estrapolano i rischi di sopravvivenza in R?

Question

Ho un oggetto survfit. Un riassunto sopravvissuto per il mio anno di interesse t=0:50 è abbastanza facile.

summary(survfit, t=0:50) 

Dà la sopravvivenza a ogni t.

C'è un modo per ottenere il rischio per ogni t (in questo caso, il pericolo da t-1 a t in ogni t = 0: 50)? Voglio ottenere la media e l'intervallo di confidenza (o errore standard) per i rischi relativi alla curva di Kaplan Meier.

Questo sembra facile da fare quando una distribuzione è in forma (ad esempio type="hazard" in flexsurvreg) ma non riesco a capire come farlo per un normale oggetto di sopravvivenza. Suggerimenti?

Answer 1

È un po 'complicato poiché il pericolo è una stima di una probabilità istantanea (e si tratta di dati discreti), ma la funzione basehaz potrebbe essere di aiuto, ma restituisce solo il rischio cumulativo. Quindi avresti ancora dovuto eseguire un ulteriore passaggio.

Ho anche avuto fortuna con la funzione muhaz. Dalla sua documentazione:

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

Non sono sicuro che il modo migliore per ottenere con l'intervallo di confidenza del 95%, ma il bootstrap potrebbe essere un approccio.

#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

Come supplemento alla risposta di Mike, si potrebbe modellare il numero di eventi da una distribuzione di Poisson, invece di una distribuzione normale. Il tasso di pericolo può quindi essere calcolato tramite una distribuzione gamma. Il codice sarebbe diventato:

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") 

Come si può vedere le stime negative per il limite inferiore del tasso di rischio ora sono andati.