R'de survfit'ten nasıl tehlikeler çıkartacağım?

Question

survfit nesnesine sahibim. Benim t=0:50 yıl faiz için bir özet survfit yeterince kolaydır.

summary(survfit, t=0:50) 

Her t'de sağkalımı verir.

her t tehlike elde etmek için bir yol (bu durumda, her bir t = 0, t için t-1 tehlike: 50) var mı? Kaplan Meier eğrisi ile ilgili tehlikeler için ortalama ve güven aralığı (veya standart hata) almak istiyorum.

Bu, bir dağıtım uygun olduğunda yapmak kolay gibi görünüyor (örneğin flexsurvreg içinde type="hazard"), ancak bunu düzenli bir survfit nesnesi için nasıl yapacağımı anlayamıyorum. Öneriler?

Answer 1

Tehlike, anlık bir olasılığa ilişkin bir tahmin olduğundan (ve bu, ayrık bir veri olduğundan) biraz zor olabilir, ancak basehaz işlevi bazı yardımcı olabilir, ancak yalnızca kümülatif tehlikeyi döndürür. Yani hala fazladan bir adım atmak zorunda kalacaktınız.

Ben de muhaz fonksiyonu ile şans oldu. onun belgelerden;

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

ben% 95 güven aralığında almanın en iyi yolu emin değilim, ama önyükleme bir yaklaşım olabilir.

#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 

Ahmet'in cevabı çalışmasına ek olarak

Answer 2

, tek bir Poisson dağılımının yerine normal dağılım ile etkinlik sayısını modellemek olabilir. Tehlike oranı daha sonra bir gama dağılımı ile hesaplanabilir. Kod olacaktı:

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

gibi şimdi gitti alt tehlike oranının bağlı olumsuz tahminler görülebilir .