[R] plot(cox.zph()): customize xlab & ylab

[Dan Bebber]

Dear all,

I've modified the plot.cox.zph function to allow
customized xlab and ylab (see below). Someone might
like to confirm that it works.

Thanks for all the assistance.

Dan
___________________________________

plot.cox.zph <- function (x, resid = TRUE, se = TRUE,
df = 4, nsmo = 40, var, 
    xlab="Time",ylab = paste("Beta(t) for",
dimnames(yy)[[2]]),...) 
{
    xx <- x$x
    yy <- x$y
    d <- nrow(yy)
    df <- max(df)
    nvar <- ncol(yy)
    pred.x <- seq(from = min(xx), to = max(xx), length
= nsmo)
    temp <- c(pred.x, xx)
    lmat <- ns(temp, df = df, intercept = TRUE)
    pmat <- lmat[1:nsmo, ]
    xmat <- lmat[-(1:nsmo), ]
    qmat <- qr(xmat)
    if (se) {
        bk <- backsolve(qmat$qr[1:df, 1:df], diag(df))
        xtx <- bk %*% t(bk)
        seval <- d * ((pmat %*% xtx) * pmat) %*%
rep(1, df)
    }
    if (missing(var)) 
        var <- 1:nvar
    else {
        if (is.character(var)) 
            var <- match(var, dimnames(yy)[[2]])
        if (any(is.na(var)) || max(var) > nvar ||
min(var) < 
            1) 
            stop("Invalid variable requested")
    }
    if (x$transform == "log") {
        xx <- exp(xx)
        pred.x <- exp(pred.x)
    }
    else if (x$transform != "identity") {
        xtime <- as.numeric(dimnames(yy)[[1]])
        apr1 <- approx(xx, xtime, seq(min(xx),
max(xx), length = 17)[2 * 
            (1:8)])
        temp <- signif(apr1$y, 2)
        apr2 <- approx(xtime, xx, temp)
        xaxisval <- apr2$y
        xaxislab <- rep("", 8)
        for (i in 1:8) xaxislab[i] <- format(temp[i])
    }
    for (i in var) {
        y <- yy[, i]
        yhat <- pmat %*% qr.coef(qmat, y)
        if (resid) 
            yr <- range(yhat, y)
        else yr <- range(yhat)
        if (se) {
            temp <- 2 * sqrt(x$var[i, i] * seval)
            yup <- yhat + temp
            ylow <- yhat - temp
            yr <- range(yr, yup, ylow)
        }
        if (x$transform == "identity") 
            plot(range(xx), yr, type = "n", xlab =
xlab, ylab = ylab[i], 
                ...)
        else if (x$transform == "log") 
            plot(range(xx), yr, type = "n", xlab =
xlab, ylab = ylab[i], 
                log = "x", ...)
        else {
            plot(range(xx), yr, type = "n", xlab =
xlab, ylab = ylab[i], 
                axes = FALSE, ...)
            axis(1, xaxisval, xaxislab)
            axis(2)
            box()
        }
        if (resid) 
            points(xx, y)
        lines(pred.x, yhat)
        if (se) {
            lines(pred.x, yup, lty = 2)
            lines(pred.x, ylow, lty = 2)
        }
    }
}


--- Thomas Lumley <tlumley at u.washington.edu> wrote:

> On Mon, 11 Jul 2005, Adaikalavan Ramasamy wrote:
> 
> > I am not sure if there is an easy way around this.
> An ugly hack is to
> > make a copy the function "survival:::plot.cox.zph"
> and make your
> > modified function. But there are others in the
> list who might know
> > neater solutions.
> 
> If you then send a patch to the package maintainer
> it stops being an ugly 
> hack and turns into an example of collaborative
> open-source development :)
> 
>  	-thomas
> 
>