Survival and hazard rate function estimation using wavelets was covered in Section 4.8 and Appendix C of my book. There were a few small typos, mainly in the code, in the book and these are explained here.
Alternatively, you can download all the necessary (and hopefully correct) code here
You can read this code into your R session via:
This code also contains the
fan data set as an example.
Here is a some example code. (see below if you don't know how to get packages)
library("muhaz") # Useful survivor function estimator library("wavethresh") # For the wavelet's we're going to do. kphaz.fan <- kphaz.fit(time=fan.z, status=fan.delta) # Do usual fit # # Now plot the kphaz estimation # plot(kphaz.fan$time, kphaz.fan$haz, type="l", ylim=c(0, 0.1)) # # Now compute and plot a wavelet estimate (on the same plot, in red) # hazest.fan <- hazest(z=fan.z, delta=fan.delta, nbins=32, levN=1, levD=1) lines(hazest.fan$bins[1:30], hazest.fan$haz[1:30], col=2) # # Note, that the previous example had minimal smoothing (controlled # by levN and levD. Also, we don't plot the last two values as they are # stupidly large (and the method needs improvement here). # # Let's compute another one with more smoothing (by increasing the # levN and levD values. The # # levN: value controls the smoothing for the sub-estimate (numerator) # levD: value controls the smoothing for the denominator. # hazest.fan3 <- hazest(z=fan.z, delta=fan.delta, nbins=32, levN=3, levD=3) lines(hazest.fan3$bins[1:30], hazest.fan3$haz[1:30], col=3)
Note, in general, this example code does not choose the smoothing parameters automatically (but someone might like to write a simple method to do this, e.g. cross-validation), it should be rewritten to more carefully modify/monitor values of the hazard function that get uncontrollably large (ie commensurately reducing the numerator). A better method would also investigate the use of wavelets smoother than Haar.
I look forward very much to seeing the result of someone who makes these improvements!
I would like to thank various citizens of the countries of Iraq and Iran who have helped me improve this code and explanation.
Getting packages. To get the
muhaz package, you might type:
install.packages("muhaz")you might have to answer questions about which repository/mirror you want to use, select one and then you should see the package install (then return to the instructions above. You also need to ensure you've got the