NIK and RRPlots
The libraries
library(survival)
library(FRESA.CAD)
## Loading required package: Rcpp
## Loading required package: stringr
## Loading required package: miscTools
## Loading required package: Hmisc
##
## Attaching package: 'Hmisc'
## The following objects are masked from 'package:base':
##
## format.pval, units
## Loading required package: pROC
## Type 'citation("pROC")' for a citation.
##
## Attaching package: 'pROC'
## The following objects are masked from 'package:stats':
##
## cov, smooth, var
library(Biobase)
## Loading required package: BiocGenerics
##
## Attaching package: 'BiocGenerics'
## The following object is masked from 'package:pROC':
##
## var
## The following objects are masked from 'package:stats':
##
## IQR, mad, sd, var, xtabs
## The following objects are masked from 'package:base':
##
## anyDuplicated, append, as.data.frame, basename, cbind, colnames,
## dirname, do.call, duplicated, eval, evalq, Filter, Find, get, grep,
## grepl, intersect, is.unsorted, lapply, Map, mapply, match, mget,
## order, paste, pmax, pmax.int, pmin, pmin.int, Position, rank,
## rbind, Reduce, rownames, sapply, setdiff, sort, table, tapply,
## union, unique, unsplit, which.max, which.min
## Welcome to Bioconductor
##
## Vignettes contain introductory material; view with
## 'browseVignettes()'. To cite Bioconductor, see
## 'citation("Biobase")', and for packages 'citation("pkgname")'.
##
## Attaching package: 'Biobase'
## The following object is masked from 'package:Hmisc':
##
## contents
## The following object is masked from 'package:miscTools':
##
## rowMedians
op <- par(no.readonly = TRUE)
pander::panderOptions('digits', 3)
pander::panderOptions('table.split.table', 400)
pander::panderOptions('keep.trailing.zeros',TRUE)
if (!require("BiocManager", quietly = TRUE))
{
install.packages("BiocManager")
BiocManager::install("seventyGeneData")
}
## Bioconductor version '3.15' is out-of-date; the current release version '3.17'
## is available with R version '4.3'; see https://bioconductor.org/install
library(seventyGeneData)
data(vanDeVijver)
class(vanDeVijver)
## [1] "ExpressionSet"
## attr(,"package")
## [1] "Biobase"
Getting the clinical data
pdata <- pData(vanDeVijver)
ROC Plots
table(pdata$Posnodes)
n y 151 144
pander::pander(table(pdata$TTMevent))
pmroc <- plotModels.ROC(cbind(pdata$TTMevent,-pdata$C1used),name="NIK",thr= -0.4) ## Using paper threshold
par(op)
Node positive data
Node Negative data
RRPlot Cox Model
timeinterval <- 5 # Five years
h0 <- sum(pdata$TTMevent & pdata$RFS <= timeinterval)
h0 <- h0/sum((pdata$RFS > timeinterval) | (pdata$TTMevent==1))
mcox <- coxph(Surv(RFS,TTMevent)~C1used,pdata)
pander::pander(summary(mcox)$coefficients)
| C1used |
-1.5 |
0.224 |
0.263 |
-5.69 |
1.3e-08 |
index <- predict(mcox,pdata)
rdata <- cbind(pdata$TTMevent,ppoisGzero(index,h0))
RRAnalysisCox <- RRPlot(rdata,atRate=c(0.90,0.95),
timetoEvent=pdata$RFS,
title="NIK: Breast Cancer",
ysurvlim=c(0.00,1.0),
riskTimeInterval=timeinterval)
par(op)
Expected time to event
toinclude <- rdata[,1]==1 | pdata$RFS > 2.0*timeinterval
obstiemToEvent <- pdata[toinclude,"RFS"]
tmin<-min(obstiemToEvent)
sum(toinclude)
[1] 167
timetoEvent <- meanTimeToEvent(rdata[toinclude,2],timeinterval)
tmax<-max(c(obstiemToEvent,timetoEvent))
lmfit <- lm(obstiemToEvent~0+timetoEvent)
sm <- summary(lmfit)
pander::pander(sm)
| timetoEvent |
0.736 |
0.0364 |
20.2 |
1.24e-46 |
Fitting linear model: obstiemToEvent ~ 0 +
timetoEvent
| 167 |
4.86 |
0.711 |
0.709 |
plot(timetoEvent,obstiemToEvent,
col=1+rdata[toinclude,1],
xlab="Expected time",
ylab="Observed time",
main="Expected vs. Observed",
xlim=c(tmin,tmax),
ylim=c(tmin,tmax),
log="xy")
lines(x=c(tmin,tmax),y=lmfit$coefficients*c(tmin,tmax),lty=1,col="blue")
text(tmin+0.01*(tmax-tmin),tmax,sprintf("R.sq=%3.2f",sm$r.squared),cex=0.7)
MADerror2 <- mean(abs(timetoEvent-obstiemToEvent))
pander::pander(MADerror2)
4.67
RRPlot Cox Adjusted Model
This time we will include Lymph node status from pathology report and
Estrogen receptor alpha expression measurement from microarray
mcox <- coxph(Surv(RFS,TTMevent)~C1used*(ESR1 + Posnodes),pdata)
pander::pander(summary(mcox)$coefficients)
| C1used |
-0.403 |
0.668 |
0.629 |
-0.640 |
0.52186 |
| ESR1 |
0.123 |
1.131 |
0.255 |
0.481 |
0.63079 |
| Posnodesy |
-0.305 |
0.737 |
0.217 |
-1.401 |
0.16112 |
| C1used:ESR1 |
-1.913 |
0.148 |
0.739 |
-2.588 |
0.00966 |
| C1used:Posnodesy |
0.378 |
1.460 |
0.583 |
0.649 |
0.51661 |
index <- predict(mcox,pdata)
rdata <- cbind(pdata$TTMevent,ppoisGzero(index,h0))
RRAnalysisAdCox <- RRPlot(rdata,atRate=c(0.90,0.95),
timetoEvent=pdata$RFS,
title="Adjusted: Breast Cancer",
ysurvlim=c(0.00,1.0),
riskTimeInterval=timeinterval)
par(op)
Expected time to event
timetoEvent <- meanTimeToEvent(rdata[toinclude,2],timeinterval)
tmax<-max(c(obstiemToEvent,timetoEvent))
lmfit <- lm(obstiemToEvent~0+timetoEvent)
sm <- summary(lmfit)
pander::pander(sm)
| timetoEvent |
0.639 |
0.0339 |
18.9 |
3.63e-43 |
Fitting linear model: obstiemToEvent ~ 0 +
timetoEvent
| 167 |
5.09 |
0.682 |
0.68 |
plot(timetoEvent,obstiemToEvent,
col=1+rdata[toinclude,1],
xlab="Expected time",
ylab="Observed time",
main="Expected vs. Observed",
xlim=c(tmin,tmax),
ylim=c(tmin,tmax),
log="xy")
lines(x=c(tmin,tmax),y=lmfit$coefficients*c(tmin,tmax),lty=1,col="blue")
text(tmin+0.01*(tmax-tmin),tmax,sprintf("R.sq=%3.2f",sm$r.squared),cex=0.7)
MADerror2 <-c(MADerror2,mean(abs(timetoEvent-obstiemToEvent)))
pander::pander(MADerror2)
4.67 and 5.51
Calibrating the index
calprob <- CoxRiskCalibration(mcox,pdata,"TTMevent","RFS")
pander::pander(c(h0=calprob$h0,
Gain=calprob$hazardGain,
DeltaTime=calprob$timeInterval),
caption="Cox Calibration Parameters")
h0 <- calprob$h0
timeinterval <- calprob$timeInterval;
rdata <- cbind(pdata$TTMevent,calprob$prob)
RRAnalysisCalAdCox <- RRPlot(rdata,atRate=c(0.90,0.95),
timetoEvent=pdata$RFS,
title="Cal. NIK: Breast Cancer",
ysurvlim=c(0.00,1.0),
riskTimeInterval=timeinterval)
par(op)
Expected time to event
timetoEvent <- meanTimeToEvent(rdata[toinclude,2],timeinterval)
tmax<-max(c(obstiemToEvent,timetoEvent))
lmfit <- lm(obstiemToEvent~0+timetoEvent)
sm <- summary(lmfit)
pander::pander(sm)
| timetoEvent |
0.623 |
0.033 |
18.9 |
3.63e-43 |
Fitting linear model: obstiemToEvent ~ 0 +
timetoEvent
| 167 |
5.09 |
0.682 |
0.68 |
plot(timetoEvent,obstiemToEvent,
col=1+rdata[toinclude,1],
xlab="Expected time",
ylab="Observed time",
main="Expected vs. Observed",
xlim=c(tmin,tmax),
ylim=c(tmin,tmax),
log="xy")
lines(x=c(tmin,tmax),y=lmfit$coefficients*c(tmin,tmax),lty=1,col="blue")
text(tmin+0.01*(tmax-tmin),tmax,sprintf("R.sq=%3.2f",sm$r.squared),cex=0.7)
MADerror2 <-c(MADerror2,mean(abs(timetoEvent-obstiemToEvent)))
pander::pander(MADerror2)
4.67, 5.51 and 5.69
Comparing Risks
Comparing concordance Index
## Correlation Index
cindex <- RRAnalysisCI$c.index$cstatCI
## Cox Index
cindex <- rbind(cindex,RRAnalysisCox$c.index$cstatCI)
## Adjusted Cox Index
cindex <- rbind(cindex,RRAnalysisAdCox$c.index$cstatCI)
## Adjusted and Calibrated Cox Index
cindex <- rbind(cindex,RRAnalysisCalAdCox$c.index$cstatCI)
rownames(cindex) <- c("CI","Cox","Adj. Cox","Cal. Adj. Cox")
pander::pander(cindex)
| CI |
0.698 |
0.698 |
0.647 |
0.745 |
| Cox |
0.698 |
0.698 |
0.649 |
0.745 |
| Adj. Cox |
0.707 |
0.707 |
0.659 |
0.751 |
| Cal. Adj. Cox |
0.707 |
0.706 |
0.660 |
0.751 |
Comparing Risk Ratios Index
## Correlation Index
RRratio <- c(RR=RRAnalysisCI$keyPoints$RR[1],
LCI=RRAnalysisCI$keyPoints$RR_LCI[1],
UCI=RRAnalysisCI$keyPoints$RR_UCI[1])
## Cox Index
RRratio <- rbind(RRratio,c(RR=RRAnalysisCox$keyPoints$RR[1],
LCI=RRAnalysisCox$keyPoints$RR_LCI[1],
UCI=RRAnalysisCox$keyPoints$RR_UCI[1]))
## Adjusted Cox Index
RRratio <- rbind(RRratio,c(RR=RRAnalysisAdCox$keyPoints$RR[1],
LCI=RRAnalysisAdCox$keyPoints$RR_LCI[1],
UCI=RRAnalysisAdCox$keyPoints$RR_UCI[1]))
## Adjusted and Calibrated Cox Index
RRratio <- rbind(RRratio,c(RR=RRAnalysisCalAdCox$keyPoints$RR[1],
LCI=RRAnalysisCalAdCox$keyPoints$RR_LCI[1],
UCI=RRAnalysisCalAdCox$keyPoints$RR_UCI[1]))
rownames(RRratio) <- c("CI","Cox","Adj. Cox","Cal. Adj. Cox")
pander::pander(RRratio)
| CI |
3.81 |
2.08 |
6.96 |
| Cox |
3.81 |
2.08 |
6.96 |
| Adj. Cox |
3.42 |
1.93 |
6.07 |
| Cal. Adj. Cox |
3.42 |
1.93 |
6.07 |
Comparing logRank values
## Correlation Index
SurvDif <- c(chisq=RRAnalysisCI$surdif$chisq,pvalue=RRAnalysisCI$surdif$pvalue)
## Cox Index
SurvDif <- rbind(SurvDif,c(chisq=RRAnalysisCox$surdif$chisq,pvalue=RRAnalysisCox$surdif$pvalue))
## Adjusted Cox Index
SurvDif <- rbind(SurvDif,c(chisq=RRAnalysisAdCox$surdif$chisq,pvalue=RRAnalysisAdCox$surdif$pvalue))
## Adjusted and Calibrated Cox Index
SurvDif <- rbind(SurvDif,c(chisq=RRAnalysisCalAdCox$surdif$chisq,pvalue=RRAnalysisCalAdCox$surdif$pvalue))
rownames(SurvDif) <- c("CI","Cox","Adj. Cox","Cal. Adj. Cox")
pander::pander(SurvDif)
| CI |
28.1 |
7.97e-07 |
| Cox |
28.1 |
7.97e-07 |
| Adj. Cox |
28.5 |
6.46e-07 |
| Cal. Adj. Cox |
28.5 |
6.46e-07 |
Comparing Sensitivity
## Correlation Index
sensi <- RRAnalysisCI$ROCAnalysis$sensitivity
## Cox Index
sensi <- rbind(sensi,RRAnalysisCox$ROCAnalysis$sensitivity)
## Adjusted Cox Index
sensi <- rbind(sensi,RRAnalysisAdCox$ROCAnalysis$sensitivity)
## Adjusted and Calibrated Cox Index
sensi <- rbind(sensi,RRAnalysisCalAdCox$ROCAnalysis$sensitivity)
rownames(sensi) <- c("CI","Cox","Adj. Cox","Cal. Adj. Cox")
pander::pander(sensi)
| CI |
0.891 |
0.813 |
0.944 |
| Cox |
0.891 |
0.813 |
0.944 |
| Adj. Cox |
0.891 |
0.813 |
0.944 |
| Cal. Adj. Cox |
0.891 |
0.813 |
0.944 |
Comparing Specificity
## Correlation Index
speci <- RRAnalysisCI$ROCAnalysis$specificity
## Cox Index
speci <- rbind(speci,RRAnalysisCox$ROCAnalysis$specificity)
## Adjusted Cox Index
speci <- rbind(speci,RRAnalysisAdCox$ROCAnalysis$specificity)
## Adjusted and Calibrated Cox Index
speci <- rbind(speci,RRAnalysisCalAdCox$ROCAnalysis$specificity)
rownames(speci) <- c("CI","Cox","Adj. Cox","Cal. Adj. Cox")
pander::pander(speci)
| CI |
0.397 |
0.328 |
0.469 |
| Cox |
0.397 |
0.328 |
0.469 |
| Adj. Cox |
0.397 |
0.328 |
0.469 |
| Cal. Adj. Cox |
0.397 |
0.328 |
0.469 |
Comparing O/E
OERatio <- NULL
## Cox Index
OERatio <- rbind(OERatio,RRAnalysisCox$OERatio$estimate)
## Adjusted Cox Index
OERatio <- rbind(OERatio,RRAnalysisAdCox$OERatio$estimate)
## Adjusted and Calibrated Cox Index
OERatio <- rbind(OERatio,RRAnalysisCalAdCox$OERatio$estimate)
rownames(OERatio) <- c("Cox","Adj. Cox","Cal. Adj. Cox")
pander::pander(OERatio)
| Cox |
0.871 |
0.709 |
1.06 |
0.178 |
| Adj. Cox |
0.904 |
0.737 |
1.10 |
0.344 |
| Cal. Adj. Cox |
0.923 |
0.752 |
1.12 |
0.445 |
Comparing O/Acum
OARatio <- NULL
## Cox Index
OARatio <- rbind(OARatio,RRAnalysisCox$OARatio$estimate)
## Adjusted Cox Index
OARatio <- rbind(OARatio,RRAnalysisAdCox$OARatio$estimate)
## Adjusted and Calibrated Cox Index
OARatio <- rbind(OARatio,RRAnalysisCalAdCox$OARatio$estimate)
rownames(OARatio) <- c("Cox","Adj. Cox","Cal. Adj. Cox")
pander::pander(OARatio)
| Cox |
1.263 |
1.029 |
1.53 |
0.0217 |
| Adj. Cox |
1.290 |
1.051 |
1.57 |
0.0128 |
| Cal. Adj. Cox |
0.997 |
0.812 |
1.21 |
1.0000 |
Comparing NetBenefit
NetBen <- NULL
## Cox Index
NetBen <- rbind(NetBen,RRAnalysisCox$keyPoints$NetBenefit)
## Adjusted Cox Index
NetBen <- rbind(NetBen,RRAnalysisAdCox$keyPoints$NetBenefit)
## Adjusted and Calibrated Cox Index
NetBen <- rbind(NetBen,RRAnalysisCalAdCox$keyPoints$NetBenefit)
colnames(NetBen) <- rownames(RRAnalysisCox$keyPoints)
rownames(NetBen) <- c("Cox","Adj. Cox","Cal. Adj. Cox")
pander::pander(NetBen)
| Cox |
0.223 |
0.232 |
0.222 |
0.239 |
0.252 |
-0.02795 |
| Adj. Cox |
0.230 |
0.249 |
0.202 |
0.252 |
0.274 |
0.00617 |
| Cal. Adj. Cox |
0.196 |
0.217 |
0.168 |
0.220 |
0.244 |
0.01345 |
Compare the ROC AUC
pander::pander(pROC::roc.test(RRAnalysisCI$ROCAnalysis$ROC.analysis$roc.predictor,
RRAnalysisAdCox$ROCAnalysis$ROC.analysis$roc.predictor))
DeLong’s test for two correlated ROC curves:
RRAnalysisCI$ROCAnalysis$ROC.analysis$roc.predictor and
RRAnalysisAdCox$ROCAnalysis$ROC.analysis$roc.predictor
| -0.988 |
0.323 |
two.sided |
0.692 |
0.702 |
