How can I extract hazards from residues in R?

Question

I have a survfit object. The summary for my t=0:50 years of interest is quite simple.

 summary(survfit, t=0:50) 

This gives survival at every t.

Is there a way to get a danger for each t (in this case, a danger from t-1 to t in each t = 0: 50)? I want to get the average and confidence intervals (or standard error) for the dangers associated with the Kaplan Meyer curve.

This is easy to do when the distribution is suitable (for example, type="hazard" in flexsurvreg ), but I can’t figure out how to do this for a regular surfactant object. Suggestions?

Answer 1

This is a little complicated, since the danger is an estimate of the instantaneous probability (and this is discrete data), but the basehaz function can provide some help, but it returns only the cumulative danger. Thus, you still have to perform an additional step.

I was also lucky with the muhaz function. From the documentation:

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

I'm not sure the best way to get a 95% confidence interval, but bootstraping might be one approach.

 #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

As a complement to Mike’s answer, one could model the number of events by the Poisson distribution instead of the normal distribution. The degree of danger can then be calculated by distributing gamma radiation. The code will look like this:

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

As can be seen, negative assessments of the lower boundary of the danger rate now disappeared.