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