Nothing
R/plotBrier.R
In riskRegression: Risk Regression Models and Prediction Scores for Survival Analysis with Competing Risks
Defines functions
plotBrier
Documented in
plotBrier
### plotBrier.R ---
#----------------------------------------------------------------------
## author: Thomas Alexander Gerds
## created: Feb 23 2017 (11:07)
## Version:
## last-updated: Sep 6 2023 (09:53)
## By: Thomas Alexander Gerds
## Update #: 77
#----------------------------------------------------------------------
##
### 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(7)
##' 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"={data.table::copy(x$Brier$score)},
"ipa"={data.table::copy(x$Brier$score)},
"contrasts"={data.table::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)) {
if (any(!(models %in% unique(pframe$models))))
stop(paste0("Fitted object does not contain models named: ",paste0(models,collapse = ", "),
"\nAvailable models are named: ",paste0(unique(pframe$models),collapse = ", ")))
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)
do.call("plot",control$plot)
if (which%in%c("score","ipa")){
## Brier
# not a very nice solution but this fixes the problem with the plot
model.list <- unique(pframe[["model"]])
times <- unique(pframe[["times"]])
for (mod in model.list){
thisline <- control$line
thisline$col=thisline$col[[as.character(mod)]]
thisline$lwd=thisline$lwd[[as.character(mod)]]
thisline$lty=thisline$lty[[as.character(mod)]]
thisline$pch=thisline$pch[[as.character(mod)]]
thisline$type=thisline$type[[as.character(mod)]]
thisline$x=times
if (which =="ipa"){
thisline$y=pframe[model==mod][["IPA"]] ## does this work?
}else{
thisline$y=pframe[model==mod][["Brier"]]
}
do.call("lines",thisline)}
}
else{
## delta Brier
# not a very nice solution but this fixes the problem with the plot
contrast.list <- unique(pframe[["contrast"]])
times <- unique(pframe[["times"]])
for (con in contrast.list){
thisline <- control$line
thisline$col=thisline$col[[as.character(con)]]
thisline$lwd=thisline$lwd[[as.character(con)]]
thisline$lty=thisline$lty[[as.character(con)]]
thisline$pch=thisline$pch[[as.character(con)]]
thisline$type=thisline$type[[as.character(con)]]
thisline$x=times
thisline$y=pframe[contrast==con][["delta.Brier"]]
do.call("lines",thisline)
}
}
## 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
Try the riskRegression package in your browser
Any scripts or data that you put into this service are public.
riskRegression documentation built on Sept. 8, 2023, 6:12 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions plotBrier
Documented in plotBrier
### plotBrier.R ---
#----------------------------------------------------------------------
## author: Thomas Alexander Gerds
## created: Feb 23 2017 (11:07)
## Version:
## last-updated: Sep 6 2023 (09:53)
## By: Thomas Alexander Gerds
## Update #: 77
#----------------------------------------------------------------------
##
### 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(7)
##' 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"={data.table::copy(x$Brier$score)},
"ipa"={data.table::copy(x$Brier$score)},
"contrasts"={data.table::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)) {
if (any(!(models %in% unique(pframe$models))))
stop(paste0("Fitted object does not contain models named: ",paste0(models,collapse = ", "),
"\nAvailable models are named: ",paste0(unique(pframe$models),collapse = ", ")))
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)
do.call("plot",control$plot)
if (which%in%c("score","ipa")){
## Brier
# not a very nice solution but this fixes the problem with the plot
model.list <- unique(pframe[["model"]])
times <- unique(pframe[["times"]])
for (mod in model.list){
thisline <- control$line
thisline$col=thisline$col[[as.character(mod)]]
thisline$lwd=thisline$lwd[[as.character(mod)]]
thisline$lty=thisline$lty[[as.character(mod)]]
thisline$pch=thisline$pch[[as.character(mod)]]
thisline$type=thisline$type[[as.character(mod)]]
thisline$x=times
if (which =="ipa"){
thisline$y=pframe[model==mod][["IPA"]] ## does this work?
}else{
thisline$y=pframe[model==mod][["Brier"]]
}
do.call("lines",thisline)}
}
else{
## delta Brier
# not a very nice solution but this fixes the problem with the plot
contrast.list <- unique(pframe[["contrast"]])
times <- unique(pframe[["times"]])
for (con in contrast.list){
thisline <- control$line
thisline$col=thisline$col[[as.character(con)]]
thisline$lwd=thisline$lwd[[as.character(con)]]
thisline$lty=thisline$lty[[as.character(con)]]
thisline$pch=thisline$pch[[as.character(con)]]
thisline$type=thisline$type[[as.character(con)]]
thisline$x=times
thisline$y=pframe[contrast==con][["delta.Brier"]]
do.call("lines",thisline)
}
}
## 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
Try the riskRegression package in your browser
Any scripts or data that you put into this service are public.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.