Nothing
### plotBrier.R ---
#----------------------------------------------------------------------
## author: Thomas Alexander Gerds
## created: Feb 23 2017 (11:07)
## Version:
## last-updated: Mar 9 2022 (15:49)
## By: Thomas Alexander Gerds
## Update #: 74
#----------------------------------------------------------------------
##
### Commentary:
##
### Change Log:
#----------------------------------------------------------------------
##
### Code:
##' Plot Brier score curves
##' @title Plot Brier curve
#' @param x Object obtained with \code{Score}
#' @param models Choice of models to plot
#' @param which Character. Either \code{"score"} to show AUC or
#' \code{"contrasts"} to show differences between AUC.
#' @param xlim Limits for x-axis
#' @param ylim Limits for y-axis
#' @param xlab Label for x-axis
#' @param ylab Label for y-axis
#' @param col line color
#' @param lwd line width
#' @param lty line style
#' @param cex point size
#' @param pch point style
#' @param type line type
#' @param axes Logical. If \code{TRUE} draw axes.
#' @param percent Logical. If \code{TRUE} scale y-axis in percent.
#' @param conf.int Logical. If \code{TRUE} draw confidence shadows.
#' @param legend Logical. If \code{TRUE} draw legend.
#' @param ... Used for additional control of the subroutines: plot,
#' axis, lines, legend. See \code{\link{SmartControl}}.
##' @examples
##' # survival
##' library(survival)
##' library(prodlim)
##' set.seed(10)
##' ds1=sampleData(40,outcome="survival")
##' ds2=sampleData(40,outcome="survival")
##' f1 <- coxph(Surv(time,event)~X1+X3+X5+X7+X9,data=ds1,x=TRUE)
##' f2 <- coxph(Surv(time,event)~X2+X4+X6+X8+X10,data=ds1,x=TRUE)
##' xscore <- Score(list(f1,f2),formula=Hist(time,event)~1,data=ds2,times=0:12,metrics="brier")
##' plotBrier(xscore)
#' @export
#'
#'
plotBrier <- function(x,
models,
which="score",
xlim,
ylim,
xlab,
ylab,
col,
lwd,
lty=1,
cex=1,
pch=1,
type="l",
axes=1L,
percent=1L,
conf.int=0L,
legend=1L,
...){
times=contrast=model=se=Brier=lower=upper=delta.Brier=reference=NULL
## cbPalette <- c("#999999", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")
cbbPalette <- c("#000000", "#E69F00", "#56B4E9", "#009E73", "#F0E442", "#0072B2", "#D55E00", "#CC79A7")
which <- tolower(which[1])
pframe <- switch(which,
"score"={copy(x$Brier$score)},
"ipa"={copy(x$Brier$score)},
"contrasts"={copy(x$Brier$contrasts)},
{stop("argument 'which' has to be either 'score' for Brier, 'ipa' for IPA, or 'contrasts' for differences in Brier.")})
if (length(pframe$times)<2) stop(paste("Need at least two time points for plotting time-dependent Brier. Object has only ",length(pframe$times),"times"))
if (!missing(models)) pframe <- pframe[model %in% models]
if (which=="score"){
mm <- unique(pframe$model)
if ("se"%in%names(pframe)){
pframe[is.na(se)×==0,lower:=0]
pframe[is.na(se)×==0,upper:=0]
}
}else{
if (which =="ipa"){
mm <- unique(pframe[model!="Null model",model[1]])
pframe <- pframe[model%in%mm]
if ("se"%in%names(pframe)){
pframe[is.na(se)×==0,lower:=0]
pframe[is.na(se)×==0,upper:=0]
}
}else{
pframe[,contrast:=factor(paste(model,reference,sep=" - "))]
mm <- unique(pframe$contrast)
if ("se"%in%names(pframe)){
pframe[is.na(se)×==0,lower:=0]
pframe[is.na(se)×==0,upper:=0]
}
}
}
lenmm <- length(mm)
if(missing(xlab)) xlab <- "Time"
if(missing(ylab)) ylab <- switch(which,
"score"="Brier score",
"ipa"="Index of prediction accuracy",
expression(paste(Delta, " Brier score")))
if(missing(col)) col <- rep(cbbPalette,length.out=lenmm)
names(col) <- mm
if(missing(lwd)) lwd <- 2
lwd <- rep(lwd,length.out=lenmm)
names(lwd) <- mm
pch <- rep(pch,length.out=lenmm)
names(pch) <- mm
type <- rep(type,length.out=lenmm)
names(type) <- mm
if(missing(lwd)) lty <- 1
lty <- rep(lty,length.out=lenmm)
names(lty) <- mm
if (missing(xlim)) xlim <- pframe[,range(times)]
if (missing(ylim)){
if (which%in%c("score","ipa")) {
ylim <- c(0,.3)
if (which=="ipa"){
ylim <- c(0,max(pframe$IPA,na.rm=0))
}
axis2.DefaultArgs <- list(side=2,
las=2,
at=seq(0,ylim[2],ylim[2]/4),
mgp=c(4,1,0))
} else{
ylim <- c(floor(10*min(pframe$lower))/10,ceiling(10*max(pframe$upper))/10)
yat <- seq(ylim[1],ylim[2],0.05)
## this is a strange behaviour of R: seq(-0.6,.1,0.05)
## [1] -6.000000e-01 -5.500000e-01 -5.000000e-01 -4.500000e-01 -4.000000e-01 -3.500000e-01 -3.000000e-01 -2.500000e-01
## [9] -2.000000e-01 -1.500000e-01 -1.000000e-01 -5.000000e-02 1.110223e-16 5.000000e-02 1.000000e-01
yat <- round(100*yat)/100
## axis2.DefaultArgs <- list(side=2,las=2,at=seq(ylim[1],ylim[2],abs(ylim[2]-ylim[1])/4),mgp=c(4,1,0))
axis2.DefaultArgs <- list(side=2,las=2,at=yat,mgp=c(4,1,0))
}
}else{
axis2.DefaultArgs <- list(side=2,las=2,at=seq(ylim[1],ylim[2],abs(ylim[2]-ylim[1])/4),mgp=c(4,1,0))
}
lines.DefaultArgs <- list(pch=pch,type=type,cex=cex,lwd=lwd,col=col,lty=lty)
axis1.DefaultArgs <- list(side=1,las=1,at=seq(0,xlim[2],xlim[2]/4))
if (which%in%c("score","ipa")){
legend.DefaultArgs <- list(legend=mm,lwd=lwd,col=col,lty=lty,cex=cex,bty="n",y.intersp=1.3,x="topleft")
if (which=="ipa") legend.DefaultArgs$x="bottomleft"
} else{
legend.DefaultArgs <- list(legend=as.character(unique(pframe$contrast)),lwd=lwd,col=col,lty=lty,cex=cex,bty="n",y.intersp=1.3,x="topleft")
}
plot.DefaultArgs <- list(x=0,y=0,type = "n",ylim = ylim,xlim = xlim,ylab=ylab,xlab=xlab)
control <- prodlim::SmartControl(call= list(...),
keys=c("plot","lines","legend","axis1","axis2"),
ignore=NULL,
ignore.case=TRUE,
defaults=list("plot"=plot.DefaultArgs,
"lines"=lines.DefaultArgs,
"legend"=legend.DefaultArgs,
"axis1"=axis1.DefaultArgs,
"axis2"=axis2.DefaultArgs),
forced=list("plot"=list(axes=FALSE),
"axis1"=list(side=1)),
verbose=TRUE)
if (which%in%c("score","ipa")){
## Brier
do.call("plot",control$plot)
pframe[,{thisline <- control$line
thisline$col=thisline$col[[as.character(model[1])]]
thisline$lwd=thisline$lwd[[as.character(model[1])]]
thisline$lty=thisline$lty[[as.character(model[1])]]
thisline$pch=thisline$pch[[as.character(model[1])]]
thisline$type=thisline$type[[as.character(model[1])]]
thisline$x=times
if (which =="ipa"){
thisline$y=IPA
}else{
thisline$y=Brier
}
do.call("lines",thisline)},by=model]
}else{
## delta Brier
do.call("plot",control$plot)
pframe[,{thisline <- control$line;
thisline$col=thisline$col[[as.character(contrast[1])]];
thisline$lwd=thisline$lwd[[as.character(contrast[1])]];
thisline$lty=thisline$lty[[as.character(contrast[1])]];
thisline$pch=thisline$pch[[as.character(contrast[1])]];
thisline$type=thisline$type[[as.character(contrast[1])]];
thisline$x=times;
thisline$y=delta.Brier;
do.call("lines",thisline)},by=contrast]
}
## legend
if (!(is.logical(legend[[1]]) && legend[[1]]==FALSE)){
do.call("legend",control$legend)
}
## x-axis
if (conf.int==TRUE){
dimcol <- sapply(col,function(cc){prodlim::dimColor(cc)})
names(dimcol) <- names(col)
if (which=="score"){
pframe[,polygon(x=c(times,rev(times)),y=c(lower,rev(upper)),col=dimcol[[as.character(model)]],border=NA),by=model]
}else{
pframe[,polygon(x=c(times,rev(times)),y=c(lower,rev(upper)),col=dimcol[[as.character(contrast)]],border=NA),by=contrast]
}
}
if (axes){
control$axis2$labels <- paste(100*control$axis2$at,"%")
do.call("axis",control$axis1)
do.call("axis",control$axis2)
}
invisible(pframe)
}
#----------------------------------------------------------------------
### plotBrier.R ends here
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.