# 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)
}