# R/plotCalibration.R In riskRegression: Risk Regression Models and Prediction Scores for Survival Analysis with Competing Risks

#### Documented in plotCalibration

### plotCalibration.R ---
#----------------------------------------------------------------------
## author: Thomas Alexander Gerds
## created: Feb 23 2017 (11:15)
## Version:
## last-updated: Mar  9 2022 (15:48)
##           By: Thomas Alexander Gerds
##     Update #: 384
#----------------------------------------------------------------------
##
### Commentary:
##
### Change Log:
#----------------------------------------------------------------------
##
### Code:
##' Plot Calibration curves for risk prediction models
##'
##' In uncensored data, the observed frequency of the outcome event is calculated locally at the predicted risk.
##' In right censored data with and without competing risks, the actual risk
##' is calculated using the Kaplan-Meier and the Aalen-Johansen method, respectively,
##' locally at the predicted risk.
##' @title Plot Calibration curve
##' @export
##' @param x Object obtained with function \code{Score}
##' @param models Choice of models to plot
##' @param times Time point specifying the prediction horizon.
##' @param method The method for estimating the calibration curve(s):
##' \itemize{
##' \item{\code{"quantile"}}{The observed proportion at predicted risk value 'p'
##' is obtained in groups
##' defined by quantiles of the predicted event probabilities of all subjects.
##' The number of groups is controlled by argument \code{q}.}
##' \item{\code{"nne"}}: {The observed proportion at predicted risk value 'p' is obtained based
##' on the subjects whose predicted risk is inside a nearest
##' neighborhood around the value 'p'. The larger the
##' bandwidth the more subjects are included in the current neighborhood.
##' }}
##' @param cens.method For right censored data only. How observed proportions are calculated. Either \code{"jackknife"} or \code{"local"}:
##' \itemize{
##' \item{\code{"jackknife"}}{Compute a running mean of the jackknife pseudovalues across neighborhoods/groups of the predicted risks.
##' Here we rely on the
##' assumption that censoring is independent of the event time and the covariates, see References. }
##' \item{\code{"local"}}{Compute the Kaplan-Meier estimator in absence of competing risks and the Aalen-Johansen estimator in presence of competing risks
##' locally like a running mean in neighborhoods of the predicted risks. The widths of the neighborhoods
##' are defined according to method.
##' }}
##' @param round If \code{TRUE} predicted probabilities are rounded to
##'     two digits before smoothing. This may have a considerable
##'     effect on computing efficiency in large data sets.
##' @param bandwidth The bandwidth for \code{method="nne"}
##' @param q The number of quantiles for \code{method="quantile"} and
##'     \code{bars=TRUE}.
##' @param bars If \code{TRUE}, use barplots to show calibration.
##' @param hanging Barplots only. If \code{TRUE}, hang bars
##'     corresponding to observed frequencies (estimated actual risk)  at the value of the
##'     corresponding prediction.
##' @param names Barplots only. Names argument passed to
##'     \code{names.arg} of \code{barplot}.
##' @param pseudo If \code{TRUE} show pseudo values (only for right
##'     censored data).
##' @param rug If \code{TRUE} show rug plot at the predictions
##' @param show.frequencies Barplots only. If \code{TRUE}, show
##'     frequencies above the bars.
##' @param plot If \code{FALSE}, do not plot the results, just return
##'     a plottable object.
##' @param add If \code{TRUE} the line(s) are added to an existing
##'     plot.
##' @param diag If \code{FALSE} no diagonal line is drawn.
##' @param legend Logical. If \code{TRUE} draw legend.
##' @param auc.in.legend Logical. If \code{TRUE} add AUC to legend.
##' @param brier.in.legend Logical. If \code{TRUE} add Brier score to
##'     legend.
##' @param axes If \code{FALSE} no axes are drawn.
##' @param xlim Limits of x-axis.
##' @param ylim Limits of y-axis.
##' @param xlab Label for y-axis.
##' @param ylab Label for x-axis.
##' @param col Vector with colors, one for each element of
##' @param lwd Vector with line widths, one for each element of
##' @param lty lwd Vector with line style, one for each element of
##' @param pch Passed to \code{\link{lines}}.
##' @param type Passed to \code{\link{lines}}.
##' @param percent If TRUE axes labels are multiplied by 100 and thus
##'     interpretable on a percent scale.
##' @param na.action what to do with NA values. Passed to
##' @param cex Default cex used for legend and labels.
##' @param ... Used to control the subroutines: plot, axis, lines,
##'     barplot, legend, addtable2plot, points (pseudo values), rug. See
##' @examples
##' library(prodlim)
##' # binary
##' set.seed(10)
##' db=sampleData(100,outcome="binary")
##' fb1=glm(Y~X1+X5+X7,data=db,family="binomial")
##' fb2=glm(Y~X1+X3+X6+X7,data=db,family="binomial")
##' xb=Score(list(model1=fb1,model2=fb2),Y~1,data=db,
##'           plots="cal")
##' plotCalibration(xb,brier.in.legend=TRUE)
##' plotCalibration(xb,bars=TRUE,model="model1")
##' plotCalibration(xb,models=1,bars=TRUE,names.cex=1.3)
##'
##' # survival
##' library(survival)
##' library(prodlim)
##' dslearn=sampleData(56,outcome="survival")
##' dstest=sampleData(100,outcome="survival")
##' fs1=coxph(Surv(time,event)~X1+X5+X7,data=dslearn,x=1)
##' fs2=coxph(Surv(time,event)~strata(X1)+X3+X6+X7,data=dslearn,x=1)
##' xs=Score(list(Cox1=fs1,Cox2=fs2),Surv(time,event)~1,data=dstest,
##'           plots="cal",metrics=NULL)
##' plotCalibration(xs)
##' plotCalibration(xs,cens.method="local",pseudo=1)
##' plotCalibration(xs,method="quantile")
##'
##'
##' # competing risks
##'
##' \dontrun{
##' data(Melanoma)
##' f1 <- CSC(Hist(time,status)~age+sex+epicel+ulcer,data=Melanoma)
##' f2 <- CSC(Hist(time,status)~age+sex+logthick+epicel+ulcer,data=Melanoma)
##' x <- Score(list(model1=f1,model2=f2),Hist(time,status)~1,data=Melanoma,
##'            cause= 2,times=5*365.25,plots="cal")
##' plotCalibration(x)
##' }
##'
plotCalibration <- function(x,
models,
times,
method="nne",
cens.method,
round=TRUE,
bandwidth=NULL,
q=10,
bars=FALSE,
hanging=FALSE,
names="quantiles",
pseudo=FALSE,
rug,
## boxplot=FALSE,
show.frequencies=FALSE,
plot=TRUE,
auc.in.legend,
brier.in.legend,
xlim=c(0,1),
ylim=c(0,1),
xlab=ifelse(bars,"Risk groups","Predicted risk"),
ylab,
col,
lwd,
lty,
pch,
type,
percent=TRUE,
na.action=na.fail,
cex=1,
...){
if (x$response.type!="binary" && missing(cens.method)){ cens.method <- "local" message("The default method for estimating calibration curves based on censored data has changed for riskRegression version 2019-9-8 or higher\nSet cens.method=\"jackknife\" to get the estimate using pseudo-values.\nHowever, note that the option \"jackknife\" is sensititve to violations of the assumption that the censoring is independent of both the event times and the covariates.\nSet cens.method=\"local\" to suppress this message.") } # {{{ plot frame model=risk=event=status=NULL if (missing(ylab)){ if (bars==TRUE){ ylab="" } else{ if (x$cens.type=="rightCensored"){
ylab="Estimated actual risk"
} else{
ylab="Observed frequency"
}
}
}
if (missing(auc.in.legend))
auc.in.legend <- ("auc" %in% tolower(x$metrics)) if (missing(brier.in.legend)) brier.in.legend <- ("auc" %in% tolower(x$metrics))
if (missing(pseudo) & missing(rug))
if (x$cens.type=="rightCensored"){ showPseudo <- FALSE showRug <- FALSE } else{ showPseudo <- FALSE showRug <- TRUE } else{ if (missing(pseudo) || pseudo[[1]]==FALSE) showPseudo <- FALSE else showPseudo <- TRUE if (missing(rug) || rug[[1]]==FALSE) showRug <- FALSE else showRug <- TRUE } pframe <- x$Calibration$plotframe if (is.null(pframe)) stop("Object has no information for calibration plot.\nYou should call the function \"riskRegression::Score\" with plots=\"calibration\".") Rvar <- grep("^(ReSpOnSe|pseudovalue)$",names(pframe),value=TRUE)
if (!missing(models)){
fitted.models <- pframe[,unique(model)]
if (all(models%in%fitted.models)){
pframe <- pframe[model%in%models]
} else{
if (all(is.numeric(models)) && (max(models)<=length(fitted.models))){
models <- fitted.models[models]
pframe <- pframe[model%in%models]
}else{
stop(paste0("The requested models: ",
models,
"\ndo not all match the fitted models: ",
paste0(fitted.models,collapse=", ")))
}
}
}
data.table::setkey(pframe,model)
if (x$response.type!="binary"){ if (missing(times)){ tp <- max(pframe[["times"]]) if (length(unique(pframe$times))>1)
warning("Time point not specified, use max of the available times: ",tp)
} else{ ## can only do one time point
tp <- times[[1]]
if (!(tp%in%unique(pframe$times))) stop(paste0("Requested time ",times[[1]]," is not in object")) } pframe <- pframe[times==tp] }else tp <- NULL ## plot(0,0,type="n",ylim = 0:1,xlim = 0:1,axes=FALSE,xlab = xlab,ylab = ylab) NF <- pframe[,length(unique(model))] if (bars){ method="quantile" if (!(NF==1)) stop(paste0("Barplots work only for one risk prediction model at a time. Provided are ",NF, "models.")) } # }}} # {{{ lines if (missing(lwd)) lwd <- rep(3,NF) if (missing(col)) { if (bars) col <- c("grey90","grey30") else{ cbbPalette <- c("#000000", "#E69F00", "#56B4E9", "#009E73", "#D55E00", "#0072B2", "#CC79A7", "#F0E442") if (NF>length(cbbPalette)) col <- 1:NF else col <- cbbPalette[1:NF] } } if (missing(type)){ if (method=="quantile"){ type <- rep("b",length.out=NF) } else{ type <- rep("l",length.out=NF) } } if (missing(lty)) lty <- rep(1, length.out=NF) if (missing(pch)) pch <- rep(1, length.out=NF) if (length(lwd) < NF) lwd <- rep(lwd, length.out=NF) if (length(type) < NF) type <- rep(type,length.out=NF) if (length(lty) < NF) lty <- rep(lty,length.out=NF) if (length(col) < NF) col <- rep(col,length.out=NF) if (length(pch) < NF) pch <- rep(pch,length.out=NF) # }}} # {{{ SmartControl modelnames <- pframe[,unique(model)] axis1.DefaultArgs <- list(side=1,las=1,at=seq(0,xlim[2],xlim[2]/4)) axis2.DefaultArgs <- list(side=2,las=2,at=seq(0,ylim[2],ylim[2]/4),mgp=c(4,1,0)) if (is.character(legend[[1]])|| legend[[1]]==TRUE||auc.in.legend==TRUE|| brier.in.legend==TRUE){ legend.data <- getLegendData(object=x, models=models, times=tp, auc.in.legend=auc.in.legend, brier.in.legend=brier.in.legend, drop.null.model=TRUE) if (is.character(legend)) legend.text <- legend else{ legend.text <- unlist(legend.data[,1]) } nrows.legend <- NROW(legend.data) if (nrows.legend==1){ legend.lwd <- NA } else{ legend.lwd <- lwd } legend.DefaultArgs <- list(legend=legend.text,lwd=legend.lwd,col=col,ncol=1,lty=lty,cex=cex,bty="n",y.intersp=1,x="topleft",title="") if (NCOL(legend.data)>1){ addtable2plot.DefaultArgs <- list(yjust=1.18, cex=cex, table=legend.data[,-1,drop=FALSE]) }else{ addtable2plot.DefaultArgs <- NULL } }else{ legend.DefaultArgs <- NULL addtable2plot.DefaultArgs <- NULL } if (bars){ legend.DefaultArgs <- list(legend=modelnames,col=col,cex=cex,bty="n",x="topleft") names.DefaultArgs <- list(cex=.7*par()$cex,y=c(-abs(diff(ylim))/15,-abs(diff(ylim))/25))
frequencies.DefaultArgs <- list(cex=.7*par()$cex,percent=FALSE,offset=0) } else{ if (length(modelnames)<=1 && sum(auc.in.legend+brier.in.legend)==0){ legend=FALSE } } if(bars){ if (x$cens.type=="rightCensored")
legend.DefaultArgs$legend <- c("Predicted risk", "Estimated actual risk") else legend.DefaultArgs$legend <- c("Predicted risk",
"Observed frequency")
}
lines.DefaultArgs <- list(pch=pch,type=type,cex=cex,lwd=lwd,col=col,lty=lty)
abline.DefaultArgs <- list(lwd=1,col="red")
if (missing(ylim)){
if (showPseudo[1] && !bars[1]){
ylim <- range(pframe[[Rvar]])
}
else
ylim <- c(0,1)
}
if (missing(xlim)){
xlim <- c(0,1)
}
plot.DefaultArgs <- list(x=0,
y=0,
type = "n",
ylim = ylim,
xlim = xlim,
ylab=ylab,
xlab=xlab)
rug.DefaultArgs <- list(col=col,lwd=lwd,side=1,ticksize=0.03)
pseudo.DefaultArgs <- list(col=col,cex=cex,density=55)
barplot.DefaultArgs <- list(ylim = ylim,
col=col,
axes=FALSE,
ylab=ylab,
xlab=xlab,
beside=TRUE,
legend.text=NULL,
cex.axis=cex,
cex.lab=par()$cex.lab, cex.names=cex) if (bars){ control <- prodlim::SmartControl(call= list(...), keys=c("barplot","legend","addtable2plot","axis2","abline","names","frequencies"), ignore=NULL, ignore.case=TRUE, defaults=list("barplot"=barplot.DefaultArgs, "addtable2plot"=addtable2plot.DefaultArgs, "abline"=abline.DefaultArgs, "legend"=legend.DefaultArgs, "names"=names.DefaultArgs, "frequencies"=frequencies.DefaultArgs, "axis2"=axis2.DefaultArgs), forced=list("abline"=list(h=0)), verbose=TRUE) }else{ control <- prodlim::SmartControl(call= list(...), keys=c("plot","rug","pseudo","lines","legend","addtable2plot","axis1","axis2"), ignore=NULL, ignore.case=TRUE, defaults=list("plot"=plot.DefaultArgs,"pseudo"=pseudo.DefaultArgs,"rug"=rug.DefaultArgs, "addtable2plot"=addtable2plot.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) } # }}} # {{{ smoothing method <- match.arg(method,c("quantile","nne")) getXY <- function(f){ risk=NULL p <- pframe[model==f,risk] jackF <- pframe[model==f][[Rvar]] switch(method, "quantile"={ if (length(q)==1) groups <- quantile(p,seq(0,1,1/q)) else{ groups <- q } xgroups <- (groups[-(length(groups))]+groups[-1])/2 pcut <- cut(p,groups,include.lowest=TRUE) ## if (x$cens.type=="rightCensored"){
if (x$response.type=="binary"){ plotFrame=data.frame(Pred=tapply(p,pcut,mean), Obs=pmin(1,pmax(0,tapply(jackF,pcut,mean)))) }else{ if(x$response.type=="survival"){
censcode <- pframe[status==0,status[1]]
qfit <- prodlim::prodlim(prodlim::Hist(time,status,cens.code=censcode)~pcut,data=pframe)
plotFrame=data.frame(Pred=tapply(p,pcut,mean),
Obs=predict(qfit,
newdata=data.frame(pcut=levels(pcut)),
cause=1,
mode="matrix",
times=tp,type="cuminc"))
}else{
censcode <- pframe[status==0,event[1]]
qfit <- prodlim::prodlim(prodlim::Hist(time,event,cens.code=censcode)~pcut,data=pframe)
n.cause <- match(x$cause,x$states)
plotFrame=data.frame(Pred=tapply(p,pcut,mean),
Obs=predict(qfit,
newdata=data.frame(pcut=levels(pcut)),
cause=n.cause,
mode="matrix",
times=tp,type="cuminc"))
}
}
attr(plotFrame,"quantiles") <- groups
plotFrame
},
"nne"={
## Round probabilities to 2 digits
## to avoid memory explosion ...
## a difference in the 3 digit should
## not play a role for the single subject.
if (round==TRUE){
if (!is.null(bandwidth) && bandwidth[1]>=1){
## message("No need to round predicted probabilities to calculate calibration in the large")
} else{
p <- round(p,2)
}
}
p <- na.omit(p)
if (no <- length(attr(p,"na.action")))
warning("calPlot: removed ",no," missing values in risk prediction.",call.=FALSE,immediate.=TRUE)
if (is.null(bandwidth)){
## if (length(p)>length(apppred[,f+1])){
## bw <- prodlim::neighborhood(apppred[,f+1])$bandwidth ## }else{ if (length(unique(p))==1) bw <- 1 else bw <- prodlim::neighborhood(p)$bandwidth
## }
} else{
bw <- bandwidth
}
if (bw>=1){
## calibration in the large
plotFrame <- data.frame(Pred=mean(p),Obs=mean(jackF))
} else{
if (x$response.type=="binary"){ nbh <- prodlim::meanNeighbors(x=p,y=jackF,bandwidth=bw) plotFrame <- data.frame(Pred=nbh$uniqueX,Obs=nbh$averageY) ## plotFrame=data.frame(Pred=tapply(p,pcut,mean), ## Obs=pmin(1,pmax(0,tapply(jackF,pcut,mean)))) }else{ ## local Kaplan-Meier/Aalen-Johansen if (cens.method=="local"){ if(x$response.type=="survival"){
censcode <- pframe[status==0,status[1]]
pfit <- prodlim::prodlim(prodlim::Hist(time,status,cens.code=censcode)~p,data=pframe,bandwidth=bandwidth)
plotFrame=data.frame(Pred=sort(unique(p)),
Obs=predict(pfit,
newdata=data.frame(p=sort(unique(p))),
cause=1,
mode="matrix",
times=tp,type="cuminc"))
}else{
censcode <- pframe[status==0,event[1]]
## pframe[,Event:=factor(event,levels=1:(length(x$states)+1),labels=c(x$states,censcode))]
pfit <- prodlim::prodlim(prodlim::Hist(time,event,cens.code=censcode)~p,data=pframe,bandwidth=bandwidth)
n.cause <- match(x$cause,x$states)
plotFrame=data.frame(Pred=sort(unique(p)),
Obs=predict(pfit,
newdata=data.frame(p=sort(unique(p))),
cause=n.cause,
mode="matrix",
times=tp,type="cuminc"))
}
} else{
## jackknife pseudo values
nbh <- prodlim::meanNeighbors(x=p,y=jackF,bandwidth=bw)
plotFrame <- data.frame(Pred=nbh$uniqueX,Obs=nbh$averageY)
}
}
}
attr(plotFrame,"bandwidth") <- bw
plotFrame
})
}
plotFrames <- lapply(modelnames,function(f){getXY(f)})
names(plotFrames) <- modelnames
# }}}
# {{{ plot and/or invisibly output the results
if (bars){
if ((is.logical(names[[1]]) && names[[1]]==TRUE)|| names[[1]] %in% c("quantiles.labels","quantiles")){
qq <- attr(plotFrames[[1]],"quantiles")
if (names[1]=="quantiles.labels"){
pp <- seq(0,1,1/q)
names <- paste0("(",
sprintf("%1.0f",100*pp[-length(pp)]),",",
sprintf("%1.0f",100*pp[-1]),
")\n",
sprintf("%1.1f",100*qq[-length(qq)])," - ",
sprintf("%1.1f",100*qq[-1]))
}
else
names <- paste0(sprintf("%1.1f",100*qq[-length(qq)])," - ",
sprintf("%1.1f",100*qq[-1]))
}
}
out <- list(plotFrames=plotFrames,
times=tp,
control=control,
legend=legend,
bars=bars,
diag=diag,
names=names,
method=method,
axes=axes,
percent=percent,
hanging=hanging,
show.frequencies=show.frequencies,
col=col,
ylim=ylim,
xlim=xlim,
ylab=ylab,
xlab=xlab,
lwd=lwd,
lty=lty,
pch=pch,
lty=lty,
type=type,
NF=NF)
if (method=="nne")
out <- c(out,list(bandwidth=sapply(plotFrames,
function(x)attr(x,"bandwidth"))))
# }}}
# {{{ do the actual plot
if (plot){
## if (boxplot){
## nbox <- length(control$lines$col)
## layout(matrix(c(1:nbox, nbox, 1), widths = 100, heights = c(100-nbox*5,rep(5,nbox))))
## for (m in 1:nbox){
## pframe[model==m,boxplot(risk,col=control$lines$col[m],ylim=c(0,1),horizontal=1L)]
## }
## }
if (out$add[1]==FALSE && !out$bars[1]){
do.call("plot",control$plot) } if (out$diag[1] && !out$bars[1]){ segments(x0=0,y0=0,x1=1,y1=1,col="gray77",lwd=2,xpd=FALSE) } if (out$diag[1] && !out$bars[1]){ segments(x0=0,y0=0,x1=1,y1=1,col="gray77",lwd=2,xpd=FALSE) } if (out$bars[1]) {
stopifnot(NF==1)
pf <- na.omit(plotFrames[[1]])
Pred <- pf$Pred Obs <- pf$Obs
if (is.character(legend[[1]]) || legend[[1]]==FALSE){
control$barplot$legend.text <- NULL
}else{
if (is.null(control$barplot$legend.text)){
control$barplot$legend.text <- control$legend$legend
}
## }else{
control$barplot$args.legend <- control$legend ## } } if (is.null(control$barplot$space)) control$barplot$space <- rep(c(1,0),length(Pred)) PredObs <- c(rbind(Pred,Obs)) control$barplot$height <- PredObs if (out$hanging){
control$barplot$offset <- c(rbind(0,Pred-Obs))
minval <- min(Pred-Obs)
if (minval<0)
negY.offset <- 0.05+seq(0,1,0.05)[prodlim::sindex(jump.times=seq(0,1,0.05),eval.times=abs(minval))]
else
negY.offset <- 0
control$barplot$ylim[1] <- min(control$barplot$ylim[1],-negY.offset)
control$names$y <- control$names$y-negY.offset
}
coord <- do.call("barplot",control$barplot) if (length(out$names)>0 && (out$names[[1]]!=FALSE) && is.character(out$names)){
if (out$names[[1]]!=FALSE && length(out$names)==(length(coord)/2)){
mids <- rowMeans(matrix(coord,ncol=2,byrow=TRUE))
text(x=mids,
## x=coord,
y=control$names$y,
## c(rbind(out$names,rbind(rep("",length(coord)/2)))), out$names,
xpd=NA,
cex=control$names$cex)
}
}
## if (out$legend) print(control$barplot$args.legend)n ## message(paste0("Bars are located at ",paste(coord,collapse=","))) if (out$hanging){
do.call("abline",control$abline) } if (out$show.frequencies){
if(out$hanging){ text(x=coord, cex=control$frequencies$cex, pos=3, y=(as.vector(rbind(Pred,Pred)) +rep(control$frequencies$offset,times=length(as.vector(coord))/2)), paste(round(100*c(rbind(Pred,Obs)),0),ifelse(control$frequencies$percent,"%",""),sep=""),xpd=NA) }else{ text(coord, pos=3, c(rbind(Pred,Obs))+control$frequencies$offset, cex=control$frequencies$cex, paste(round(100*c(rbind(Pred,Obs)),0),ifelse(control$frequencies$percent,"%",""),sep=""),xpd=NA) } } coords <- list(xcoord=coord[,1],ycoord=PredObs,offset=control$barplot$offset) out <- c(out,coords) }else{ nix <- lapply(1:NF,function(f){ pf <- out$plotFrames[[f]]
pf <- na.omit(pf)
## rug plot predictions
if (showRug){
do.call(graphics::rug,c(list(x=pf$Pred,col=control$rug$col[f],lwd=control$rug$lwd[f]*0.5))) } ## show pseudo values if (showPseudo) { if (!is.null(control$pseudo$density) & control$pseudo$density>0){ control$pseudo$col <- prodlim::dimColor(control$pseudo$col[f], control$pseudo$density) } if ((gotcha <- match("density",names(control$pseudo),nomatch=0))>0){
control$pseudo <- control$pseudo[-gotcha]
}
do.call(points, c(list(x=pframe[model==modelnames[f]][,risk],y=pframe[model==modelnames[f]][[Rvar]]),control$pseudo)) } ## add the lines lineF <- lapply(control$lines,function(x)x[[min(f,length(x))]])
lineF$x <- pf$Pred
lineF$y <- pf$Obs
do.call("lines",lineF)
})
if (is.character(out$legend[[1]]) || out$legend[[1]]==TRUE){
legend.coords <- do.call("legend",control$legend) if (!is.null(addtable2plot.DefaultArgs)){ if (is.null(control$addtable2plot[["x"]]))
control$addtable2plot[["x"]] <- legend.coords$rect$left+legend.coords$rect$w ## strange error:$y does not work as yjust is matched, thus use [["y"]]
if (is.null(control$addtable2plot[["y"]])) control$addtable2plot[["y"]] <- legend.coords$rect$top-legend.coords$rect$h
do.call(plotrix::addtable2plot,control$addtable2plot) } } } if (out$axes){
if (out$percent){ if (is.null(control$axis1$labels)){ control$axis1$labels <- paste(100*control$axis1$at,"%") } if (is.null(control$axis2$labels)){ control$axis2$labels <- paste(100*control$axis2$at,"%") } } if (!out$bars)
do.call("axis",control$axis1) do.call("axis",control$axis2)
}
}
# }}}
class(out) <- "calibrationPlot"
invisible(out)
}

#----------------------------------------------------------------------
### plotCalibration.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 March 23, 2022, 5:07 p.m.