# Piece-wise constant estimator of 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 width argument is the bin or interval size.
 
pehaz<-function(times, status = NA, width = NA, min.time = 0, max.time = NA)
{
if(is.na(status[1]))
status <- rep(1, length(times))
if(is.na(max.time))
max.time <- max(times)
if(is.na(width)) {
nu <- sum(status)
width <- (max.time - min.time)/(8 * nu^0.2)
}
cat("\nmax.time=", max.time)
cat("\nwidth=", width)
nbins <- ceiling((max.time - min.time)/width)
cat("\nnbins=", nbins)
cuts <- rep(NA, nbins + 1)
hazfun <- rep(NA, nbins)
fusum <- rep(NA, nbins)
evnum <- rep(NA, nbins)
atrisk <- rep(NA, nbins)
cuts[1] <- min.time
for(i in 1:nbins) {
left <- min.time + (i - 1) * width

# bins are considered left-continuous

right <- min.time + i * width
cuts[i + 1] <- right
ind <- (times < right) & !(times < left)

# ind indicates pts terminating within interval

fu <- times[!(times < left)]
fu[(fu > right)] <- right
fu <- fu - left
sumfu <- sum(fu)
numevnt <- sum(status[ind])
haz <- numevnt/sumfu
fusum[i] <- sumfu
evnum[i] <- numevnt
atrisk[i] <- sum(!(times < left))
hazfun[i] <- haz
}

# save call

obj.call <- match.call() # organize results

# evnum = number of events in each bin
# atrisk = number at risk in each bin
# fusum = sum of f/u time in each bin

result <- list(Call = obj.call, Width = width, Cuts = cuts, Hazard = hazfun,
Events = evnum, At.Risk = atrisk, F.U.Time = fusum)
class(result) <- "pehaz"
result

#return(cuts, hazfun, evnum, atrisk, fusum)
 
}