# Kernel-smoothed estimate of the hazard function from possibly right-censored data
# Only required argument is vector of survival times.
# The status argument is a vector of censoring indicators (1 = event, 0 = censored)
# The bandwidth argument controls the amount of smoothing.
# The estgrid argument is a vector of times for estimation.
 
kernhaz<-function(times,status=NA,estgrid=NA,bandwidth=NA)
{
if(is.na(status[1])) status <- rep(1, length(times))
if(is.na(estgrid[1])) estgrid <-seq(min(times),max(times),length=50)
n<-length(times) # num of times
d<-sum(status) # num of events
e<-length(estgrid) # num of points in estimation grid
ord<-order(times)
times<-times[ord]
status<-status[ord]
if(is.na(bandwidth)) bandwidth <- 3*(estgrid[e] - estgrid[1])/(8 * d^0.2)
hazest<-rep(NA,e)
for (i in 1:e)
{sum<-0
x0<-estgrid[i]
for (j in 1:n)
{arg<-(times[j]-x0)/bandwidth
arg<-arg^2
if (arg<1) sum<-sum+(1-arg)*status[j]/(n-j+1)}
hazest[i]<-sum*.75/bandwidth}
return(bandwidth,estgrid,hazest)
}
 
########################################################################
 
x1<-rexp(100)
c1<-c(rep(1,50),rep(0,50))
h1<-kernhaz(x1,c1,seq(0,3,.1))
pe1<-pehaz(x1,c1,max.time=3)
plot(stepfun(pe1$Cuts, pe1$Hazard), type = "s", xlab = "Time",
ylab = "Hazard Rate",xlim=c(0,3))
lines(h1$estgrid,h1$hazest)
 
########################################################################