sm.survival - S+ code by Bowman & Azzalini
 
DESCRIPTION:
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.
 
USAGE:
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, ...)
 
REQUIRED ARGUMENTS:
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.
 
OPTIONAL ARGUMENTS:
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.
 
VALUE:
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.
 
SIDE EFFECTS:
none.
 
DETAILS:
see Section 3.5 of the reference below.
 
REFERENCES:
Bowman, A.W. and Azzalini, A. (1997). Applied Smoothing Techniques for
Data Analysis: the Kernel Approach with S-Plus Illustrations. Oxford University Press, Oxford.
 
SEE ALSO:
sm.regression
 
EXAMPLES:
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(is.na(x))
if(absent(eval.points))
eval.points <- seq(min(x), max(x), length = ngrid)
if(is.na(xlab))
xlab <- deparse(substitute(x))
if(is.na(ylab))
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))
}