sm.survival - S+ code by Bowman & Azzalini
This function creates a smooth, nonparametric estimate of the quantile of the
distribution of survival data as a function of a single covariate. A weighted
Kaplan-Meier survivor function is obtained by smoothing across the covariate scale.
A small amount of smoothing is then also applied across the survival time scale in order
to achieve a smooth estimate of the quantile.
sm.survival(x, y, status, h, hv=0.05, p=0.5, status.code=1,
eval.points=NA, ngrid=50, display="lines", xlab=NA, ylab=NA,
lty=1, add=F, ...)
x: a vector of covariate values.

y: a vector of survival times.
status: an indicator of a complete survival time or a censored value.
             The value of status.code defines a complete survival time.

h: the smoothing parameter applied to the covariate scale. A normal kernel
     function is used and h is its standard deviation.
hv: a smoothing parameter applied to the weighted Kaplan-Meier
functions derived from the smoothing procedure in the covariate scale.
This ensures that a smooth estimate is obtained.

p: the quantile to be estimated at each covariate value.

status.code: the value of status which defines a complete survival time.

eval.points: the points at which the estimate will be evaluated.

ngrid: the number of points in a regular grid over the covariate scale at which
            the estimate will be evaluated, if eval.points is set to NA.

display: The setting "none" will prevent any graphical output from being produced.
                The default setting "lines" (or indeed any other value) will produce a plot of
                the data and estimate.

xlab: the label attached to the x-axis.

ylab: the label attached to the y-axis.

lty: the line type used to plot the estimate.

add: a logical value which controls whether the estimate is added to the current plot.
         Its default value is F, which creates a new plot.

...: additional graphical parameters.
a list containing the values of the estimate at the evaluation points and the values of the
smoothing parameters for the covariate and survival time scales.
see Section 3.5 of the reference below.
Bowman, A.W. and Azzalini, A. (1997). Applied Smoothing Techniques for
Data Analysis: the Kernel Approach with S-Plus Illustrations. Oxford University Press, Oxford.
x <- runif(50, 0, 10)
y <- rexp(50, 2)
z <- rexp(50, 1)
status <- rep(1, 50)
status[z<y] <- 0
y <- pmin(z, y)
sm.survival(x, y, status, h=2)
sm.survival <- function(x, y, status, h, hv = 0.05, p = 0.5, status.code = 1,
eval.points = NA,  ngrid = 50, display = "lines", xlab = NA, ylab = NA, lty = 1,
add = F, ...)
absent <- function(x)
missing(x) | any(
eval.points <- seq(min(x), max(x), length = ngrid)
xlab <- deparse(substitute(x))
ylab <- deparse(substitute(y))
if(!(display == "none" | add == T)) {
plot(x, y, type = "n", xlab = xlab, ylab = ylab, ...)
text(x[status == status.code], y[status == status.code], "x")
text(x[status != status.code], y[status != status.code], "o")
n <- length(x)
ne <- length(eval.points)
xr <- x[order(y)]
statusr <- status[order(y)]
yr <- sort(y)
w <- matrix(rep(eval.points, rep(n, ne)), ncol = n, byrow = T)
w <- w - matrix(rep(xr, ne), ncol = n, byrow = T)
w <- exp(-0.5 * (w/h)^2)
wf <- t(apply(w, 1, rev))
wf <- t(apply(wf, 1, cumsum))
wf <- t(apply(wf, 1, rev))
w <- w/wf
st <- rep(0, n)
st[statusr == status.code] <- 1
w <- 1 - w * matrix(rep(st, ne), ncol = n, byrow = T)
w <- w[, st == 1]
if(ne == 1)
w <- matrix(w, ncol = length(w))
yw <- yr[st == 1]
w <- t(apply(w, 1, cumprod))
w <- cbind(rep(1, ne), w)
j <- - t(apply(w, 1, diff))
J <- t(apply(j, 1, cumsum))
wd <- J - p
w <- exp(-0.5 * (wd/hv)^2) # * j
ns <- length(yw)
s0 <- w %*% rep(1, ns)
s1 <- (w * wd) %*% rep(1, ns)
s2 <- (w * wd^2) %*% rep(1, ns)
w <- w * (matrix(rep(s2, ns), ncol = ns) - wd * matrix(rep(s1, ns), ncol = ns))
w <- w/(matrix(rep(s2, ns), ncol = ns) * matrix(rep(s0, ns), ncol = ns) - matrix(rep(s1, ns), ncol = ns)^2)
estimate <- w %*% yw
if(!(display == "none"))
lines(eval.points, estimate, lty = lty)
invisible(list(estimate = estimate, eval.points = eval.points, h = h, hv = hv))