kmrs.S
#
# kmrs.S
#
# S code for Kaplan-Meier and reduced sample
# estimates of a distribution function
# from _histograms_ of censored data.
#
# kaplan.meier()
# reduced.sample()
# km.rs()
#
# $Revision: 3.10 $ $Date: 2009/03/26 09:47:44 $
#
# The functions in this file produce vectors `km' and `rs'
# where km[k] and rs[k] are estimates of F(breaks[k+1]),
# i.e. an estimate of the c.d.f. at the RIGHT endpoint of the interval.
#
"kaplan.meier" <-
function(obs, nco, breaks, upperobs=0) {
# obs: histogram of all observations : min(T_i,C_i)
# nco: histogram of noncensored observations : T_i such that T_i <= C_i
# breaks: breakpoints (vector or 'breakpts' object, see breaks.S)
# upperobs: number of observations beyond rightmost breakpoint
#
breaks <- as.breakpts(breaks)
n <- length(obs)
if(n != length(nco))
stop("lengths of histograms do not match")
check.hist.lengths(nco, breaks)
#
#
# reverse cumulative histogram of observations
d <- revcumsum(obs) + upperobs
#
# product integrand
s <- ifelse(d > 0, 1 - nco/d, 1)
#
km <- 1 - cumprod(s)
# km has length n; km[i] is an estimate of F(r) for r=breaks[i+1]
#
widths <- diff(breaks$val)
lambda <- - log(ifelse(s > 0, s, 1))/widths
# lambda has length n; lambda[i] is an estimate of
# the average of \lambda(r) over the interval (breaks[i],breaks[i+1]).
#
return(list(km=km, lambda=lambda))
}
"reduced.sample" <-
function(nco, cen, ncc, show=FALSE, uppercen=0)
# nco: histogram of noncensored observations: T_i such that T_i <= C_i
# cen: histogram of all censoring times: C_i
# ncc: histogram of censoring times for noncensored obs:
# C_i such that T_i <= C_i
#
# Then nco[k] = #{i: T_i <= C_i, T_i \in I_k}
# cen[k] = #{i: C_i \in I_k}
# ncc[k] = #{i: T_i <= C_i, C_i \in I_k}.
#
# The intervals I_k must span an interval [0,R] beginning at 0.
# If this interval did not include all censoring times,
# then `uppercen' must be the number of censoring times
# that were not counted in 'cen'.
{
n <- length(nco)
if(n != length(cen) || n != length(ncc))
stop("histogram lengths do not match")
#
# denominator: reverse cumulative histogram of censoring times
# denom(r) = #{i : C_i >= r}
# We compute
# cc[k] = #{i: C_i > breaks[k]}
# except that > becomes >= for k=0.
#
cc <- revcumsum(cen) + uppercen
#
#
# numerator
# #{i: T_i <= r <= C_i }
# = #{i: T_i <= r, T_i <= C_i} - #{i: C_i < r, T_i <= C_i}
# We compute
# u[k] = #{i: T_i <= C_i, T_i <= breaks[k+1]}
# - #{i: T_i <= C_i, C_i <= breaks[k]}
# = #{i: T_i <= C_i, C_i > breaks[k], T_i <= breaks[k+1]}
# this ensures that numerator and denominator are
# comparable, u[k] <= cc[k] always.
#
u <- cumsum(nco) - c(0,cumsum(ncc)[1:(n-1)])
rs <- u/cc
#
# Hence rs[k] = u[k]/cc[k] is an estimator of F(r)
# for r = breaks[k+1], i.e. for the right hand end of the interval.
#
if(!show)
return(rs)
else
return(list(rs=rs, numerator=u, denominator=cc))
}
"km.rs" <-
function(o, cc, d, breaks) {
# o: censored lifetimes min(T_i,C_i)
# cc: censoring times C_i
# d: censoring indicators 1(T_i <= C_i)
# breaks: histogram breakpoints (vector or 'breakpts' object)
#
breaks <- as.breakpts(breaks)
# count any data outside range of breakpoints
maxbreak <- breaks$max
uppercen <- sum(cc > maxbreak)
upperobs <- sum(o > maxbreak)
# compile histograms (breakpoints may not span data)
obs <- whist( o, breaks=breaks$val, right=maxbreak)
nco <- whist( o[d], breaks=breaks$val, right=maxbreak)
cen <- whist( cc, breaks=breaks$val, right=maxbreak)
ncc <- whist( cc[d], breaks=breaks$val, right=maxbreak)
# go
km <- kaplan.meier(obs, nco, breaks, upperobs=upperobs)
rs <- reduced.sample(nco, cen, ncc, uppercen=uppercen)
#
return(list(rs=rs, km=km$km, hazard=km$lambda,
r=breaks$val[-1], breaks=breaks$val))
}