TABLE OF CONTENTS
S3Methods/plot.splinesurv [ Functions ]
NAME
plot.splinesurv --- plot method for splinesurv objects
FUNCTION
Plots either a hazard, frailty, survival curve estimate, or coefficient density. There are a large number of options and settings, consult the package documentation for details on the options.
This function mainly works by calling the predict function to construct the curve fits. Most of the rest is just handling inputs and keeping track of the various parameters needed by different plotting routines.
SYNOPSIS
3772 plot.splinesurv <- function(x, which = c("hazard", "survival", "frailty", "coef", "all"), 3773 newdata = NULL, iter = NULL, fn = mean, plotknots = TRUE, npoints = 100, npost = 100, 3774 alpha=.05, legend = NULL, lty = 1, col = 2, lwd = 2, lty.knots = 1, col.knots = 8, 3775 lwd.knots = 1, xlab = NULL, ylab = NULL, main = NULL, xlim = NULL, ylim = NULL, 3776 tk = FALSE, ...)
SOURCE
3779 { 3780 3781 if(tk) { 3782 # call the routine that uses tcltk via the tkrplot to plot an interactive 3783 # curve viewer 3784 splinesurvtkplot(x, newdata = newdata, plotknots = plotknots, npoints = npoints, 3785 legend = legend, lty = lty, col = col, lwd = lwd, lty.knots = lty.knots, 3786 col.knots = col.knots, lwd.knots = lwd.knots, xlab = xlab, ylab = ylab, 3787 main = main, xlim = xlim, ylim = ylim, ...) 3788 return(invisible()); 3789 } 3790 # save old parameters 3791 oldask <- par("ask") 3792 which <- match.arg(which) 3793 if(which == "all") par(ask = TRUE) 3794 if(!is.null(iter) && iter <= 0) iter <- NULL 3795 # Hazard and survival curve plotting 3796 if(which == "hazard" | which == "survival" | which == "all"){ 3797 type <- if(which == "all") "hazard" else which 3798 # get the set of knots to use. If an iteration is given, use knots for 3799 # that iteration, otherwise use only boundary knots 3800 if(is.null(x$hazard$spline.knots)){ 3801 knots <- NULL 3802 }else{ 3803 if(is.null(iter)) { 3804 knots <- x$history$hazard.spline.knots[1, ];knots <- knots[knots>-Inf] 3805 knots <- c(min(knots), max(knots)) 3806 }else{ 3807 knots <- x$history$hazard.spline.knots[iter, ];knots <- knots[knots>-Inf] 3808 } 3809 } 3810 if(is.null(knots)) knots <- range(x$data$time) 3811 if(is.null(xlim)) xlim1 = range(x$data$time) else xlim1 <- xlim 3812 # set the points at which the curve should be predicted 3813 times = seq(from = max(xlim1[1], min(knots)), to = min(xlim1[2], max(knots)), 3814 length = npoints) 3815 if(is.null(xlab)) xlab1 <- "Time" else xlab1 <- xlab 3816 # set axis labels 3817 if(type == "hazard"){ 3818 if(is.null(ylab)) ylab1 <- "Hazard" else ylab1 <- ylab 3819 if(is.null(main)){ 3820 if(is.null(newdata)) main1 <- "Baseline hazard" 3821 else main1 <- "Hazard" 3822 }else{ main1 <- main } 3823 } 3824 if(type == "survival"){ 3825 if(is.null(ylab)) ylab1 <- "Survival" else ylab1 <- ylab 3826 if(is.null(main)){ 3827 if(is.null(newdata)) main1 <- "Baseline survival" 3828 else main1 <- "Survival" 3829 }else{ main1 <- main } 3830 } 3831 # a single curve is plotted if either no newdata is given, or newdata has only one row 3832 if(is.null(newdata) | (!is.null(newdata) && dim(newdata)[1] < 2)){ 3833 # estimate the hazard curve 3834 haz <- predict(x, type = type, x = times, newdata = newdata, iter = iter, fn = fn, 3835 npost = npost, alpha = alpha) 3836 plot(haz[, 1:2], type = "l", lty = lty, col = col, lwd = lwd, main = main1, 3837 xlab = xlab1, ylab = ylab1, xlim = xlim1, ylim = ylim, ...) 3838 # plot knots as vertical lines 3839 if(plotknots) { 3840 abline(v = knots, col = col.knots, lty = lty.knots, lwd = lwd.knots, ...) 3841 lines(haz, lty = lty, col = col, lwd = lwd) 3842 } 3843 # plot confidence intervals 3844 if(!is.null(alpha) & is.null(iter)){ 3845 lines(haz[, c(1, 3)], lty = 2, col = col) 3846 lines(haz[, c(1, 4)], lty = 2, col = col) 3847 } 3848 # if newdata has more than one row, multiple lines have to be predicted and plotted 3849 }else{ 3850 haz <- times 3851 # predict a curve for each row of newdata 3852 for(i in 1:dim(newdata)[1]) haz <- cbind(haz, predict(x, type = type, x = times, 3853 newdata = newdata[i, ,drop = FALSE], iter = iter, fn = fn, npost = npost)[, 2]) 3854 if(length(col) == 1 & length(lty) == 1 & length(lwd) == 1) col = 1:dim(newdata)[1] 3855 # plot all the new lines 3856 matplot(haz[, 1], haz[, -1], type = "l", col = col, lwd = lwd, lty = lty, 3857 main = main1, xlab = xlab1, ylab = ylab1, xlim = xlim1, ylim = ylim, ...) 3858 if(is.null(legend)) legend <- rownames(newdata) 3859 # plot knots as vertical lines 3860 if(plotknots) abline(v = knots, col = col.knots, lty = lty.knots, 3861 lwd = lwd.knots, ...) 3862 # add a legend 3863 if(!is.na(legend)) legend("topright", legend = legend, col = col, lty = lty, 3864 lwd = lwd) 3865 } 3866 if(which == "all") plot(x, which = "survival", newdata = newdata, iter = iter, 3867 plotknots = plotknots, npoints = npoints, npost = npost, alpha = alpha, 3868 legend = legend, lty = lty, col = col, lwd = lwd, lty.knots = lty.knots, 3869 col.knots = col.knots, xlab = xlab, ylab = ylab, main = main, xlim = xlim, 3870 ylim = ylim, tk = tk, ...) 3871 } 3872 # Frailty density plotting 3873 if(which == "frailty" | which == "all"){ 3874 knots <- x$frailty$spline.knots 3875 # construct knots to use 3876 if(is.null(iter)){ 3877 knots <- x$history$frailty.spline.knots[1, ];knots <- knots[knots>-Inf] 3878 knots <- c(min(knots), max(knots)) 3879 } else{ 3880 knots <- x$history$frailty.spline.knots[iter, ];knots <- knots[knots>-Inf] 3881 } 3882 if(is.null(knots)) knots <- range(x$posterior.mean$frailty) 3883 # limits on the graph 3884 if(is.null(xlim)) { 3885 if(hasspline(x$frailty)) xlim1 = attr(knots, "b") 3886 else xlim1 <- range(x$posterior.mean$frailty) 3887 }else{xlim1 <- xlim} 3888 # construct the set of values at which the curve shoult be estimated 3889 Ui = seq(from = max(xlim1[1], min(knots)), to = min(xlim1[2], max(knots)), 3890 length = npoints) 3891 if(is.null(xlab)) xlab1 <- "x" else xlab1 <- xlab 3892 if(is.null(ylab)) ylab1 <- "Density" else ylab1 <- ylab 3893 if(is.null(main)) main1 <- "Frailty density" else main1 <- main 3894 # comput the curve 3895 dens <- predict(x, type = "frailty", x = Ui, iter = iter, fn = fn, npost = npost, 3896 alpha = alpha) 3897 # plot the density curve 3898 plot(dens[, 1:2], type = "l", lty = lty, col = col, lwd = lwd, main = main1, 3899 xlab = xlab1, ylab = ylab1, xlim = xlim1, ylim = ylim, ...) 3900 # plot knots as vertical lines 3901 if(plotknots){ 3902 abline(v = knots, col = col.knots, lty = lty.knots, lwd = lwd.knots, ...) 3903 lines(dens, type = "l", lty = lty, col = col, lwd = lwd) 3904 } 3905 # plot CIs 3906 if(!is.null(alpha) & is.null(iter)){ 3907 lines(dens[, c(1, 3)], lty = 2, col = col) 3908 lines(dens[, c(1, 4)], lty = 2, col = col) 3909 } 3910 } 3911 # plot coefficient density 3912 if(which == "coef" | which == "all"){ 3913 burnin <- x$control$burnin 3914 if(is.null(burnin)) burnin <- x$call$control$burnin 3915 if(is.null(burnin)) burnin <- dim(x$history$coefficients)[2] 3916 if(is.null(xlab)) xlab1 <- "x" else xlab1 <- xlab 3917 if(is.null(ylab)) ylab1 <- "Posterior density" else ylab1 <- ylab 3918 coefs <- x$history$coefficients 3919 coefnames <- colnames(coefs) 3920 if(length(coefnames) > 1) par(ask = TRUE) 3921 for(i in 1:dim(coefs)[2]){ 3922 if(is.null(main)) main1 <- paste("Coefficient of", coefnames[i]) else main1 <- main 3923 betai <- coefs[burnin:(dim(coefs)[1]), i] 3924 plot(density(betai), lty = lty, col = col, lwd = lwd, main = main1, xlab = xlab1, ylab = ylab1, xlim = xlim, ylim = ylim, ...) 3925 } 3926 } 3927 par(ask = oldask) 3928 }