Jak wyodrębnić zagrożenia z przetrwania w R?

Question

Mam obiekt survfit. Podsumowanie dla mojego t=0:50 lat jest dość łatwe.

summary(survfit, t=0:50) 

Daje to przeżycie przy każdym t.

Czy istnieje sposób na uzyskanie zagrożenia dla każdego t (w tym przypadku zagrożenie od t-1 do t w każdym t = 0: 50)? Chcę uzyskać średnią i przedział ufności (lub błąd standardowy) dla zagrożeń związanych z krzywą Kaplana Meiera.

Wydaje się to łatwe do zrobienia, gdy rozkład jest odpowiedni (np. type="hazard" w flexsurvreg), ale nie mogę wymyślić, jak to zrobić dla zwykłego obiektu survfit. Propozycje?

Answer 1

Jest to nieco skomplikowane, ponieważ zagrożenie jest oszacowaniem chwilowego prawdopodobieństwa (i jest to dyskretne dane), ale funkcja basehaz może być pomocna, ale zwraca tylko skumulowane zagrożenie. Więc nadal musiałbyś wykonać dodatkowy krok.

Mam też szczęście z funkcją muhaz. Z dokumentacji:

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

nie jestem pewien, że najlepszym sposobem, aby dostać się w 95% przedziale ufności, ale ładowania może być jednym podejściem.

#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

Jako uzupełnienie odpowiedź Mike'a, można modelować liczbę zdarzeń przez rozkładu Poissona zamiast rozkładu normalnego. Stopień hazardu można następnie obliczyć za pomocą rozkładu gamma. Kod staną:

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

Jak widać negatywne prognozy dla dolnej granicy szybkości zagrożenia są już nie ma.