#' Plotting predicted risk
#'
#' Show predicted risk obtained by a risk prediction model as a function of
#' time.
#'
#' @aliases plot.riskRegression
#' @usage
#' \method{plot}{riskRegression}(x,
##' cause,
##' newdata,
##' xlab,
##' ylab,
##' xlim,
##' ylim,
##' lwd,
##' col,
##' lty,
##' axes=TRUE,
##' percent=TRUE,
##' legend=TRUE,
##' add=FALSE,
##' ...)
#' @param x Fitted object obtained with one of \code{ARR}, \code{LRR},
#' \code{riskRegression}.
#' @param cause For CauseSpecificCox models the cause of interest.
#' @param newdata A data frame containing predictor variable combinations for
#' which to compute predicted risk.
#' @param xlim See \code{plot}
#' @param ylim See \code{plot}
#' @param xlab See \code{plot}
#' @param ylab See \code{plot}
#' @param lwd A vector of line thicknesses for the regression coefficients.
#' @param col A vector of colors for the regression coefficients.
#' @param lty A vector of line types for the regression coefficients.
#' @param axes Logical. If \code{FALSE} then do not draw axes.
#' @param percent If true the y-axis is labeled in percent.
#' @param legend If true draw a legend.
#' @param add Logical. If \code{TRUE} then add lines to an existing plot.
#' @param \dots Used for transclusion of smart arguments for \code{plot},
#' \code{lines}, \code{axis} and \code{background}. See function
#' \code{\link{SmartControl}} from prodlim.
#' @author Thomas Alexander Gerds <tag@@biostat.ku.dk>
#' @keywords survival
##' @examples
##'
##' library(survival)
##' library(prodlim)
##' data(Melanoma)
##' fit.arr <- ARR(Hist(time,status)~invasion+age+strata(sex),data=Melanoma,cause=1)
##' plot(fit.arr,xlim=c(500,3000))
##'
##'
#' @export
plot.riskRegression <- function(x,
cause,
newdata,
xlab,
ylab,
xlim,
ylim,
lwd,
col,
lty,
axes=TRUE,
percent=TRUE,
legend=TRUE,
add=FALSE,
...){
# {{{ getting predicted risk
if ("CauseSpecificCox"%in%class(x))
plot.times <- x$eventTimes
else
plot.times <- x$timeVaryingEffects$coef[,"time"]
if (inherits(x=x,what="predictedRisk"))
Y <- split(x$risk,1:NROW(x$risk))
else{
if (missing(newdata)){
ff <- eval(x$call$formula)
xdat <- unique(eval(x$call$data)[all.vars(update(ff,NULL~.))])
if (NROW(xdat)<5){
if ("CauseSpecificCox"%in%class(x)){
p1 <- predictRisk(x,newdata=xdat,times=plot.times,cause=cause)
}
else{
p1 <- stats::predict(x,newdata=xdat,times=plot.times)$risk}
rownames(p1) <- paste("id",1:NROW(xdat))
}
else{
if ("CauseSpecificCox"%in%class(x)){
P1 <- predictRisk(x,newdata=eval(x$call$data),times=plot.times,cause=cause)
}
else{
P1 <- stats::predict(x,
newdata=eval(x$call$data),
times=plot.times)$risk
}
medianP1 <- P1[,prodlim::sindex(plot.times,median(plot.times))]
P1 <- P1[order(medianP1),]
p1 <- P1[round(quantile(1:NROW(P1))),]
rownames(p1) <- paste("Predicted risk",c("Min","q25","Median","q75","Max"),sep="=")
warning("Argument newdata is missing.\n",
"Shown are the cumulative incidence curves from the original data set.\nSelected are curves based on individual risk (min,q25,median,q75,max) at the median time:",
median(plot.times))
}
}
else{
p1 <- stats::predict(x,newdata=newdata,time=plot.times)$risk
}
Y <- lapply(1:NROW(p1),function(i){p1[i,]})
if (!is.null(rownames(p1)))
names(Y) <- rownames(p1)
}
nlines <- NROW(Y)
# }}}
# {{{ labels, limits etc
if (missing(ylab)) ylab <- "Cumulative incidence"
if (missing(xlab)) xlab <- "Time"
if (missing(xlim)) xlim <- c(0, max(plot.times))
if (missing(ylim)) ylim <- c(0, 1)
if (missing(lwd)) lwd <- rep(3,nlines)
if (missing(col)) col <- 1:nlines
if (missing(lty)) lty <- rep(1, nlines)
if (length(lwd) < nlines) lwd <- rep(lwd, nlines)
if (length(lty) < nlines) lty <- rep(lty, nlines)
if (length(col) < nlines) col <- rep(col, nlines)
# }}}
# {{{ defaults
plot.DefaultArgs <- list(x=0,y=0,type = "n",ylim = ylim,xlim = xlim,xlab = xlab,ylab = ylab)
lines.DefaultArgs <- list(type="s")
axis1.DefaultArgs <- list()
axis2.DefaultArgs <- list(at=seq(0,1,.25),side=2)
legend.DefaultArgs <- list(legend=names(Y),
lwd=lwd,
col=col,
lty=lty,
cex=1.5,
bty="n",
y.intersp=1.3,
x="topleft",
trimnames=TRUE)
# }}}
# {{{ smart control
smartA <- prodlim::SmartControl(call= list(...),
keys=c("plot","lines","legend","conf.int","marktime","axis1","axis2"),
ignore=c("x","type","cause","newdata","add","col","lty","lwd","ylim","xlim","xlab","ylab","legend","marktime","conf.int","automar","atrisk","timeOrigin","percent","axes","atrisk.args","conf.int.args","legend.args"),
ignore.case=TRUE,
defaults=list("plot"=plot.DefaultArgs,
"axis1"=axis1.DefaultArgs,
"axis2"=axis2.DefaultArgs,
"legend"=legend.DefaultArgs,
"lines"=lines.DefaultArgs),
forced=list("plot"=list(axes=FALSE),
"axis1"=list(side=1)),
verbose=TRUE)
# }}}
# {{{ empty plot
if (!add) {
do.call("plot",smartA$plot)
if (axes){
do.call("axis",smartA$axis1)
if (percent & is.null(smartA$axis1$labels))
smartA$axis2$labels <- paste(100*smartA$axis2$at,"%")
do.call("axis",smartA$axis2)
}
}
# }}}
# {{{ adding the lines
lines.type <- smartA$lines$type
nix <- lapply(1:nlines, function(s) {
lines(x = plot.times,
y = Y[[s]],
type = lines.type,
col = col[s],
lty = lty[s],
lwd = lwd[s])
})
# }}}
# {{{ legend
if(legend[[1]]==TRUE && !add[[1]] && !is.null(names(Y))){
if (all(smartA$legend$trimnames==TRUE) && all(sapply((nlist <- strsplit(names(Y),"=")),function(x)length(x))==2)){
smartA$legend$legend <- sapply(nlist,function(x)x[[2]])
smartA$legend$title <- unique(sapply(nlist,function(x)x[[1]]))
}
smartA$legend <- smartA$legend[-match("trimnames",names(smartA$legend))]
save.xpd <- par()$xpd
par(xpd=TRUE)
do.call("legend",smartA$legend)
par(xpd=save.xpd)
}
# }}}
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.