Risk-Evaluation: Breast Cancer Royston-Altman

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

op <- par(no.readonly = TRUE)
pander::panderOptions('digits', 3)
pander::panderOptions('table.split.table', 400)
pander::panderOptions('keep.trailing.zeros',TRUE)
source("C:/Users/jtame/Documents/GitHub/RISKPLOTS/CODE/auxplots.R")

Breast Cancer Royston-Altman data

data(gbsg, package=“survival”) and data(rotterdam, package=“survival”)

gbsgdata <- gbsg
rownames(gbsgdata) <- gbsgdata$pid
gbsgdata$pid <- NULL

odata <-rotterdam
rownames(odata) <- odata$pid
odata$pid <- NULL
odata$rfstime <- odata$rtime
odata$status <- odata$recur
odata$rtime <- NULL
odata$recur <- NULL

odata <- odata[,colnames(odata) %in% colnames(gbsgdata)]

odata$size <- 10*(odata$size=="<=20") + 
  35*(odata$size=="20-50") + 
  60*(odata$size==">50")

data <- as.data.frame(model.matrix(Surv(rfstime,status)~.*age,odata))

data$`(Intercept)` <- NULL

dataBrestCancerTrain <- cbind(time=odata[rownames(data),"rfstime"],status=odata[rownames(data),"status"],data)

colnames(dataBrestCancerTrain) <-str_replace_all(colnames(dataBrestCancerTrain),":","_")
colnames(dataBrestCancerTrain) <-str_replace_all(colnames(dataBrestCancerTrain)," ","")
colnames(dataBrestCancerTrain) <-str_replace_all(colnames(dataBrestCancerTrain),"\\.","_")
colnames(dataBrestCancerTrain) <-str_replace_all(colnames(dataBrestCancerTrain),"-","_")
colnames(dataBrestCancerTrain) <-str_replace_all(colnames(dataBrestCancerTrain),">","_")
dataBrestCancerTrain$time <- dataBrestCancerTrain$time/365 ## To years


pander::pander(table(odata[rownames(data),"status"]),caption="rotterdam")

rotterdam
0	1
1464	1518

data(gbsg, package=“survival”) data conditioning

gbsgdata <- gbsgdata[,colnames(odata)]
data <- as.data.frame(model.matrix(Surv(rfstime,status)~.*.,gbsgdata))

data$`(Intercept)` <- NULL

dataBrestCancerTest <- cbind(time=gbsgdata[rownames(data),"rfstime"],status=gbsgdata[rownames(data),"status"],data)

colnames(dataBrestCancerTest) <-str_replace_all(colnames(dataBrestCancerTest),":","_")
colnames(dataBrestCancerTest) <-str_replace_all(colnames(dataBrestCancerTest)," ","")
colnames(dataBrestCancerTest) <-str_replace_all(colnames(dataBrestCancerTest),"\\.","_")
colnames(dataBrestCancerTest) <-str_replace_all(colnames(dataBrestCancerTest),"-","_")
colnames(dataBrestCancerTest) <-str_replace_all(colnames(dataBrestCancerTest),">","_")
dataBrestCancerTest$time <- dataBrestCancerTest$time/365

pander::pander(table(odata[rownames(data),"status"]), caption="gbsg")

gbsg
0	1
499	183

Cox Modeling

ml <- BSWiMS.model(Surv(time,status)~.,data=dataBrestCancerTrain,loops=1,NumberofRepeats = 5)

—–.

sm <- summary(ml)
pander::pander(sm$coefficients)

	Estimate	lower	HR	upper	u.Accuracy	r.Accuracy	full.Accuracy	u.AUC	r.AUC	full.AUC	IDI	NRI	z.IDI	z.NRI	Delta.AUC	Frequency
age_nodes	5.33e-04	1.000	1.001	1.001	0.626	0.588	0.625	0.630	0.588	0.627	0.03069	0.4398	13.21	13.30	0.03894	1
nodes	4.40e-02	1.041	1.045	1.049	0.637	0.599	0.640	0.640	0.600	0.641	0.03197	0.4403	13.10	13.40	0.04147	1
size	7.33e-03	1.005	1.007	1.010	0.595	0.637	0.640	0.595	0.640	0.641	0.01103	0.2963	6.56	8.18	0.00124	1
age_size	7.57e-05	1.000	1.000	1.000	0.567	0.629	0.625	0.568	0.633	0.627	0.00762	0.3044	5.49	8.41	-0.00568	1
grade	2.17e-01	1.148	1.242	1.343	0.565	0.637	0.640	0.561	0.639	0.641	0.00777	0.2420	5.41	7.53	0.00214	1
age	-3.91e-03	0.995	0.996	0.998	0.513	0.626	0.640	0.513	0.626	0.641	0.00376	0.0931	4.97	2.54	0.01539	1
age_grade	-6.64e-04	0.999	0.999	1.000	0.508	0.616	0.625	0.509	0.619	0.627	0.00186	0.0574	3.44	1.57	0.00843	1
age_pgr	-2.60e-06	1.000	1.000	1.000	0.548	0.626	0.625	0.544	0.629	0.627	0.00267	0.1705	2.97	5.21	-0.00183	1

Cox Model Performance

Here we evaluate the model using the RRPlot() function.

The evaluation of the raw Cox model with RRPlot()

Here we will use the predicted event probability assuming a baseline hazard for events withing 5 years

timeinterval <- 5 # Five years

h0 <- sum(dataBrestCancerTrain$status & dataBrestCancerTrain$time <= timeinterval)
h0 <- h0/sum((dataBrestCancerTrain$time > timeinterval) | (dataBrestCancerTrain$status==1))

pander::pander(t(c(h0=h0,timeinterval=timeinterval)),caption="Initial Parameters")

Initial Parameters
h0	timeinterval
0.429	5

index <- predict(ml,dataBrestCancerTrain)
rdata <- cbind(dataBrestCancerTrain$status,ppoisGzero(index,h0))

rrAnalysisTrain <- RRPlot(rdata,atRate=c(0.90,0.80),
                     timetoEvent=dataBrestCancerTrain$time,
                     title="Train: Breast Cancer",
                     ysurvlim=c(0.00,1.0),
                     riskTimeInterval=timeinterval)

By Risk Categories

obsexp <- rrAnalysisTrain$OERatio$atThrEstimates

expObs(obsexp,"Train: Expected vs. Observed")

pander::pander(obsexp)

	Observed	L.CI	H.CI	Expected	pvalue
Total	1518	1443	1596	1451	8.08e-02
low	819	764	877	879	4.47e-02
90%	225	197	256	188	7.63e-03
80%	474	432	519	382	5.10e-06

Time to Event Analysis

rrAnalysisdata <- rrAnalysisTrain

pander::pander(wilcox.test(rrAnalysisdata$timetoEventData$eTime,rrAnalysisdata$timetoEventData$expectedTime,paired = TRUE))

Wilcoxon signed rank test with continuity correction: `rrAnalysisdata$timetoEventData$eTime` and `rrAnalysisdata$timetoEventData$expectedTime`
Test statistic	P value	Alternative hypothesis
1738186	5.24e-25 * * *	two.sided

highrisk <- rrAnalysisdata$timetoEventData$class == 2
pander::pander(wilcox.test(rrAnalysisdata$timetoEventData$eTime[highrisk],rrAnalysisdata$timetoEventData$expectedTime[highrisk],paired = TRUE))

Wilcoxon signed rank test with continuity correction: `rrAnalysisdata$timetoEventData$eTime[highrisk]` and `rrAnalysisdata$timetoEventData$expectedTime[highrisk]`
Test statistic	P value	Alternative hypothesis
90068	0.129	two.sided

timesdata <- expObsTime(rrAnalysisdata,title="Train: Expected vs Observed")

pander::pander(timesdata)

	O:Low Risk < 0.364	O:0.364 <= Risk < 0.435	O:High Risk >= 0.435	E:Low Risk < 0.364	E:0.364 <= Risk < 0.435	E:High Risk >= 0.435
1Q	3.42	1.84	1.05	6.34	4.67	2.31
Median	6.83	4.27	2.23	7.55	4.99	3.15
3Q	9.42	7.97	4.78	8.58	5.29	3.86

The Time vs. Events are not calibrated. Lets do the calibration

Uncalibrated Performance Report

pander::pander(t(rrAnalysisTrain$keyPoints),caption="Threshold values")

Threshold values
	@:0.9	@:0.8	@MAX_BACC	@MAX_RR	@SPE100	p(0.5)
Thr	0.435	0.364	0.339	0.239	1.98e-01	0.5000
RR	1.723	1.712	1.777	2.230	9.22e+01	1.7119
RR_LCI	1.617	1.602	1.653	1.839	1.91e-01	1.6049
RR_UCI	1.835	1.829	1.910	2.704	4.46e+04	1.8261
SEN	0.312	0.460	0.586	0.947	1.00e+00	0.2134
SPE	0.899	0.800	0.704	0.172	1.23e-02	0.9426
BACC	0.606	0.630	0.645	0.559	5.06e-01	0.5780
NetBenefit	0.121	0.178	0.224	0.355	3.90e-01	0.0805

pander::pander(t(rrAnalysisTrain$OERatio$estimate),caption="O/E Ratio")

O/E Ratio
O/E	Low	Upper	p.value
1.05	0.994	1.1	0.0808

pander::pander(t(rrAnalysisTrain$OE95ci),caption="O/E Mean")

O/E Mean
mean	50%	2.5%	97.5%
1.16	1.16	1.15	1.17

pander::pander(t(rrAnalysisTrain$OAcum95ci),caption="O/Acum Mean")

O/Acum Mean
mean	50%	2.5%	97.5%
1.35	1.35	1.35	1.35

pander::pander(rrAnalysisTrain$c.index$cstatCI,caption="C. Index")

mean.C Index	median	lower	upper
0.678	0.678	0.664	0.691

pander::pander(t(rrAnalysisTrain$ROCAnalysis$aucs),caption="ROC AUC")

ROC AUC
est	lower	upper
0.693	0.674	0.712

pander::pander((rrAnalysisTrain$ROCAnalysis$sensitivity),caption="Sensitivity")

Sensitivity
est	lower	upper
0.312	0.289	0.336

pander::pander((rrAnalysisTrain$ROCAnalysis$specificity),caption="Specificity")

Specificity
est	lower	upper
0.899	0.882	0.914

pander::pander(t(rrAnalysisTrain$thr_atP),caption="Probability Thresholds")

Probability Thresholds
90%	80%
0.435	0.364

pander::pander(rrAnalysisTrain$surdif,caption="Logrank test")

Logrank test Chisq = 504.157535 on 2 degrees of freedom, p = 0.000000
	N	Observed	Expected	(O-E)^2/E	(O-E)^2/V
class=0	1990	819	1150	95.4	399.4
class=1	370	225	169	18.5	20.9
class=2	622	474	199	382.0	445.9

Cox Calibration

op <- par(no.readonly = TRUE)


calprob <- CoxRiskCalibration(ml,dataBrestCancerTrain,"status","time")

( 7.156166 , 29202.06 , 1.068545 , 1518 , 1795.957 )

pander::pander(c(h0=calprob$h0,
                 Gain=calprob$hazardGain,
                 DeltaTime=calprob$timeInterval),
               caption="Cox Calibration Parameters")

h0	Gain	DeltaTime
0.714	1.38	7.16

The RRplot() of the calibrated model

rcaldata <- cbind(dataBrestCancerTrain$status,calprob$prob)


rrAnalysisCalTrain <- RRPlot(rcaldata,atRate=c(0.90,0.80),
                     timetoEvent=dataBrestCancerTrain$time,
                     title="Cal. Train: Breast Cancer",
                     ysurvlim=c(0.00,1.0),
                     riskTimeInterval=calprob$timeInterval)

By Risk Categories

obsexp <- rrAnalysisCalTrain$OERatio$atThrEstimates

expObs(obsexp,"Cal. Expected vs. Observed")

pander::pander(obsexp)

	Observed	L.CI	H.CI	Expected	pvalue
Total	1518	1443	1596	1687	2.97e-05
low	819	764	877	1022	6.03e-11
90%	225	197	256	218	6.35e-01
80%	474	432	519	444	1.54e-01

Time to Event Analysis

rrAnalysisdata <- rrAnalysisCalTrain

pander::pander(wilcox.test(rrAnalysisdata$timetoEventData$eTime,rrAnalysisdata$timetoEventData$expectedTime,paired = TRUE))

Wilcoxon signed rank test with continuity correction: `rrAnalysisdata$timetoEventData$eTime` and `rrAnalysisdata$timetoEventData$expectedTime`
Test statistic	P value	Alternative hypothesis
2322576	0.0357 *	two.sided

highrisk <- rrAnalysisdata$timetoEventData$class == 2
pander::pander(wilcox.test(rrAnalysisdata$timetoEventData$eTime[highrisk],rrAnalysisdata$timetoEventData$expectedTime[highrisk],paired = TRUE))

Wilcoxon signed rank test with continuity correction: `rrAnalysisdata$timetoEventData$eTime[highrisk]` and `rrAnalysisdata$timetoEventData$expectedTime[highrisk]`
Test statistic	P value	Alternative hypothesis
104552	0.0869	two.sided

timesdata <- expObsTime(rrAnalysisdata,title="Train: Expected vs Observed")

pander::pander(timesdata)

	O:Low Risk < 0.529	O:0.529 <= Risk < 0.613	O:High Risk >= 0.613	E:Low Risk < 0.529	E:0.529 <= Risk < 0.613	E:High Risk >= 0.613
1Q	3.42	1.84	1.05	5.46	4.01	1.99
Median	6.83	4.27	2.23	6.49	4.29	2.71
3Q	9.42	7.97	4.78	7.38	4.55	3.32

Calibrated Train Performance

pander::pander(t(rrAnalysisTrain$keyPoints),caption="Threshold values")

Threshold values
	@:0.9	@:0.8	@MAX_BACC	@MAX_RR	@SPE100	p(0.5)
Thr	0.435	0.364	0.339	0.239	1.98e-01	0.5000
RR	1.723	1.712	1.777	2.230	9.22e+01	1.7119
RR_LCI	1.617	1.602	1.653	1.839	1.91e-01	1.6049
RR_UCI	1.835	1.829	1.910	2.704	4.46e+04	1.8261
SEN	0.312	0.460	0.586	0.947	1.00e+00	0.2134
SPE	0.899	0.800	0.704	0.172	1.23e-02	0.9426
BACC	0.606	0.630	0.645	0.559	5.06e-01	0.5780
NetBenefit	0.121	0.178	0.224	0.355	3.90e-01	0.0805

pander::pander(t(rrAnalysisTrain$OERatio$estimate),caption="O/E Ratio")

O/E Ratio
O/E	Low	Upper	p.value
1.05	0.994	1.1	0.0808

pander::pander(t(rrAnalysisTrain$OE95ci),caption="O/E Mean")

O/E Mean
mean	50%	2.5%	97.5%
1.16	1.16	1.15	1.17

pander::pander(t(rrAnalysisTrain$OAcum95ci),caption="O/Acum Mean")

O/Acum Mean
mean	50%	2.5%	97.5%
1.35	1.35	1.35	1.35

pander::pander(rrAnalysisTrain$c.index$cstatCI,caption="C. Index")

mean.C Index	median	lower	upper
0.678	0.678	0.664	0.691

pander::pander(t(rrAnalysisTrain$ROCAnalysis$aucs),caption="ROC AUC")

ROC AUC
est	lower	upper
0.693	0.674	0.712

pander::pander((rrAnalysisTrain$ROCAnalysis$sensitivity),caption="Sensitivity")

Sensitivity
est	lower	upper
0.312	0.289	0.336

pander::pander((rrAnalysisTrain$ROCAnalysis$specificity),caption="Specificity")

Specificity
est	lower	upper
0.899	0.882	0.914

pander::pander(t(rrAnalysisTrain$thr_atP),caption="Probability Thresholds")

Probability Thresholds
90%	80%
0.435	0.364

pander::pander(rrAnalysisTrain$surdif,caption="Logrank test")

Logrank test Chisq = 504.157535 on 2 degrees of freedom, p = 0.000000
	N	Observed	Expected	(O-E)^2/E	(O-E)^2/V
class=0	1990	819	1150	95.4	399.4
class=1	370	225	169	18.5	20.9
class=2	622	474	199	382.0	445.9

Performance on the external data set

index <- predict(ml,dataBrestCancerTest)
pp <- predictionStats_binary(cbind(dataBrestCancerTest$status,index),plotname="Breast Cancer")

par(op)


prob <- ppoisGzero(index,h0)
rdata <- cbind(dataBrestCancerTest$status,prob)
rrCoxTestAnalysis <- RRPlot(rdata,atThr=rrAnalysisTrain$thr_atP,
                     timetoEvent=dataBrestCancerTest$time,
                     title="Test: Breast Cancer",
                     ysurvlim=c(0.00,1.0),
                     riskTimeInterval=timeinterval)

par(op)

External Data Report

pander::pander(t(rrCoxTestAnalysis$keyPoints),caption="Threshold values")

Threshold values
	@:0.435	@:0.364	@MAX_BACC	@MAX_RR	@SPE100	p(0.5)
Thr	0.4349	0.364	0.358	0.250	2.12e-01	0.4999
RR	1.7921	1.746	1.773	2.737	2.64e+01	1.7885
RR_LCI	1.5300	1.472	1.490	1.446	5.65e-02	1.5235
RR_UCI	2.0991	2.070	2.110	5.181	1.23e+04	2.0997
SEN	0.3712	0.555	0.585	0.973	1.00e+00	0.2709
SPE	0.8475	0.690	0.667	0.103	1.55e-02	0.9044
BACC	0.6094	0.623	0.626	0.538	5.08e-01	0.5876
NetBenefit	0.0956	0.142	0.150	0.255	2.86e-01	0.0642

pander::pander(t(rrCoxTestAnalysis$OERatio$estimate),caption="O/E Ratio")

O/E Ratio
O/E	Low	Upper	p.value
1.38	1.23	1.54	1.36e-07

pander::pander(rrCoxTestAnalysis$c.index,caption="C. Index")

C Index: 0.668
Dxy: 0.336
S.D.: 0.0309
n: 686
missing: 0
uncensored: 299
Relevant Pairs: 266144
Concordant: 177775
Uncertain: 203702
cstatCI:

mean.C Index median lower upper

0.668 0.669 0.639 0.699

mean.C Index	median	lower	upper
0.668	0.669	0.639	0.699

pander::pander(t(rrCoxTestAnalysis$ROCAnalysis$aucs),caption="ROC AUC")

ROC AUC
est	lower	upper
0.669	0.629	0.71

pander::pander((rrCoxTestAnalysis$ROCAnalysis$sensitivity),caption="Sensitivity")

Sensitivity
est	lower	upper
0.371	0.316	0.429

pander::pander((rrCoxTestAnalysis$ROCAnalysis$specificity),caption="Specificity")

Specificity
est	lower	upper
0.85	0.811	0.884

pander::pander(t(rrCoxTestAnalysis$thr_atP),caption="Probability Thresholds")

Probability Thresholds
90%	80%
0.435	0.364

pander::pander(rrCoxTestAnalysis$surdif,caption="Logrank test")

Logrank test Chisq = 92.126440 on 2 degrees of freedom, p = 0.000000
	N	Observed	Expected	(O-E)^2/E	(O-E)^2/V
class=0	400	133	198.1	21.376	64.3
class=1	117	55	49.3	0.668	0.8
class=2	169	111	51.7	68.135	84.0

Calibrating the index on the test data

calprob <- CoxRiskCalibration(ml,dataBrestCancerTest,"status","time")

( 5.664115 , 3012.376 , 1.295745 , 299 , 329.4124 )

pander::pander(c(h0=calprob$h0,
                 Gain=calprob$hazardGain,
                 DeltaTime=calprob$timeInterval),
               caption="Cox Calibration Parameters")

h0	Gain	DeltaTime
0.579	1	5.66

rdata <- cbind(dataBrestCancerTest$status,calprob$prob)

rrAnalysisTest <- RRPlot(rdata,atRate=c(0.90,0.80),
                     timetoEvent=dataBrestCancerTest$time,
                     title="Cal. Test: Breast Cancer",
                     ysurvlim=c(0.00,1.0),
                     riskTimeInterval=calprob$timeInterval)

By Risk Categories

obsexp <- rrAnalysisTest$OERatio$atThrEstimates

expObs(obsexp,"Cal. Expected vs. Observed")

pander::pander(obsexp)

	Observed	L.CI	H.CI	Expected	pvalue
Total	299	266.1	334.9	256.6	0.00954
low	171	146.3	198.6	143.7	0.02692
90%	43	31.1	57.9	29.9	0.02174
80%	85	67.9	105.1	80.8	0.61637

Time to Event Analysis

rrAnalysisdata <- rrAnalysisTest

pander::pander(wilcox.test(rrAnalysisdata$timetoEventData$eTime,rrAnalysisdata$timetoEventData$expectedTime,paired = TRUE))

Wilcoxon signed rank test with continuity correction: `rrAnalysisdata$timetoEventData$eTime` and `rrAnalysisdata$timetoEventData$expectedTime`
Test statistic	P value	Alternative hypothesis
26492	3.01e-69 * * *	two.sided

highrisk <- rrAnalysisdata$timetoEventData$class == 2
pander::pander(wilcox.test(rrAnalysisdata$timetoEventData$eTime[highrisk],rrAnalysisdata$timetoEventData$expectedTime[highrisk],paired = TRUE))

Wilcoxon signed rank test with continuity correction: `rrAnalysisdata$timetoEventData$eTime[highrisk]` and `rrAnalysisdata$timetoEventData$expectedTime[highrisk]`
Test statistic	P value	Alternative hypothesis
3703	0.669	two.sided

timesdata <- expObsTime(rrAnalysisdata,title="Train: Expected vs Observed")

pander::pander(timesdata)

	O:Low Risk < 0.508	O:0.508 <= Risk < 0.599	O:High Risk >= 0.599	E:Low Risk < 0.508	E:0.508 <= Risk < 0.599	E:High Risk >= 0.599
1Q	1.95	1.37	1.08	4.95	3.34	1.58
Median	3.35	2.35	1.72	5.72	3.65	2.30
3Q	4.83	3.69	3.01	6.50	3.86	2.69

After Calibration Report

pander::pander(t(rrAnalysis$keyPoints),caption="Threshold values")

Quitting from lines 373-383 (BreastCancerRoyAltman.Rmd) Error in t(rrAnalysis$keyPoints) : object ‘rrAnalysis’ not found Calls: … eval_with_user_handlers -> eval -> eval -> -> t In addition: There were 45 warnings (use warnings() to see them)

pander::pander(t(rrAnalysis$OERatio$estimate),caption="O/E Ratio")

Quitting from lines 373-383 (BreastCancerRoyAltman.Rmd) Error in t(rrAnalysis$OERatio$estimate) : object ‘rrAnalysis’ not found Calls: … eval_with_user_handlers -> eval -> eval -> -> t

pander::pander(rrAnalysis$c.index,caption="C. Index")

Quitting from lines 373-383 (BreastCancerRoyAltman.Rmd) Error in pander::pander(rrAnalysis$c.index, caption = “C. Index”) : object ‘rrAnalysis’ not found Calls: … eval_with_user_handlers -> eval -> eval ->

pander::pander(t(rrAnalysis$ROCAnalysis$aucs),caption="ROC AUC")

Quitting from lines 373-383 (BreastCancerRoyAltman.Rmd) Error in t(rrAnalysis$ROCAnalysis$aucs) : object ‘rrAnalysis’ not found Calls: … eval_with_user_handlers -> eval -> eval -> -> t

pander::pander((rrAnalysis$ROCAnalysis$sensitivity),caption="Sensitivity")

Quitting from lines 373-383 (BreastCancerRoyAltman.Rmd) Error in pander::pander((rrAnalysis$ROCAnalysis$sensitivity), caption = “Sensitivity”) : object ‘rrAnalysis’ not found Calls: … eval_with_user_handlers -> eval -> eval ->

pander::pander((rrAnalysis$ROCAnalysis$specificity),caption="Specificity")

Quitting from lines 373-383 (BreastCancerRoyAltman.Rmd) Error in pander::pander((rrAnalysis$ROCAnalysis$specificity), caption = “Specificity”) : object ‘rrAnalysis’ not found Calls: … eval_with_user_handlers -> eval -> eval ->

pander::pander(t(rrAnalysis$thr_atP),caption="Probability Thresholds")

Quitting from lines 373-383 (BreastCancerRoyAltman.Rmd) Error in t(rrAnalysis$thr_atP) : object ‘rrAnalysis’ not found Calls: … eval_with_user_handlers -> eval -> eval -> -> t

pander::pander(rrAnalysis$surdif,caption="Logrank test")

Quitting from lines 373-383 (BreastCancerRoyAltman.Rmd) Error in pander::pander(rrAnalysis$surdif, caption = “Logrank test”) : object ‘rrAnalysis’ not found Calls: … eval_with_user_handlers -> eval -> eval ->

Logistic Model

Here we train a logistic model on the same data set

## Only label subjects that present event withing five years

dataBrestCancerR <- subset(dataBrestCancerTrain, time>=5 | status==1)
dataBrestCancerR$status <- dataBrestCancerR$status * (dataBrestCancerR$time < 5)
dataBrestCancerR$time <- NULL

#ml <- BSWiMS.model(status~1,data=dataBrestCancerR,loops=20,NumberofRepeats = 5)
mlog <- BSWiMS.model(status~1,data=dataBrestCancerR,loops=1,NumberofRepeats = 5)

—–.

sm <- summary(mlog)
pander::pander(sm$coefficients)

	Estimate	lower	OR	upper	u.Accuracy	r.Accuracy	full.Accuracy	u.AUC	r.AUC	full.AUC	IDI	NRI	z.IDI	z.NRI	Delta.AUC	Frequency
age_nodes	0.001339	1.001	1.001	1.001	0.678	0.610	0.683	0.642	0.583	0.652	0.07157	0.5528	14.27	15.66	0.069089	1.0
nodes	0.087232	1.082	1.091	1.100	0.676	0.632	0.691	0.639	0.622	0.663	0.07210	0.5673	14.21	16.08	0.040845	1.0
size	0.010407	1.007	1.010	1.014	0.611	0.686	0.691	0.618	0.653	0.663	0.01523	0.3186	6.26	8.38	0.009924	1.0
age_size	0.000149	1.000	1.000	1.000	0.608	0.679	0.683	0.577	0.644	0.652	0.01440	0.3009	6.11	7.90	0.008129	1.0
grade	0.379818	1.284	1.462	1.665	0.571	0.683	0.691	0.500	0.653	0.663	0.01179	0.1809	5.72	4.76	0.010365	1.0
age_meno	-0.003201	0.995	0.997	0.998	0.571	0.684	0.683	0.500	0.652	0.652	0.00775	0.0649	4.60	1.69	0.000242	1.0
pgr	-0.000253	1.000	1.000	1.000	0.571	0.682	0.683	0.500	0.650	0.652	0.00403	0.1949	3.33	5.65	0.001397	1.0
age_grade	-0.001887	0.997	0.998	1.000	0.574	0.690	0.691	0.507	0.662	0.663	0.00273	0.0713	2.62	1.85	0.001270	1.0
age_hormon	-0.000729	0.999	0.999	1.000	0.578	0.684	0.682	0.522	0.653	0.650	0.00131	0.1436	2.17	3.93	-0.002961	0.4

Logistic Model Performance

op <- par(no.readonly = TRUE)


cprob <- predict(mlog,dataBrestCancerTrain)

rdata <- cbind(dataBrestCancerTrain$status,cprob)
rrAnalysisLogTrain <- RRPlot(rdata,atRate=c(0.90,0.80),
                     timetoEvent=dataBrestCancerTrain$time,
                     title="Logistic Train: Breast Cancer",
                     ysurvlim=c(0.00,1.0),
                     riskTimeInterval=5.0)

par(op)

By Risk Categories

obsexp <- rrAnalysisLogTrain$OERatio$atThrEstimates

expObs(obsexp,"Logistic: Expected vs. Observed")

pander::pander(obsexp)

	Observed	L.CI	H.CI	Expected	pvalue
Total	1518	1443	1596	1857	4.90e-16
low	802	747	859	1000	1.19e-10
90%	226	197	257	258	4.96e-02
80%	490	448	535	596	8.71e-06

Time to Event Analysis

rrAnalysisdata <- rrAnalysisLogTrain

pander::pander(wilcox.test(rrAnalysisdata$timetoEventData$eTime,rrAnalysisdata$timetoEventData$expectedTime,paired = TRUE))

Wilcoxon signed rank test with continuity correction: `rrAnalysisdata$timetoEventData$eTime` and `rrAnalysisdata$timetoEventData$expectedTime`
Test statistic	P value	Alternative hypothesis
2351656	0.00656 * *	two.sided

highrisk <- rrAnalysisdata$timetoEventData$class == 2
pander::pander(wilcox.test(rrAnalysisdata$timetoEventData$eTime[highrisk],rrAnalysisdata$timetoEventData$expectedTime[highrisk],paired = TRUE))

Wilcoxon signed rank test with continuity correction: `rrAnalysisdata$timetoEventData$eTime[highrisk]` and `rrAnalysisdata$timetoEventData$expectedTime[highrisk]`
Test statistic	P value	Alternative hypothesis
136782	3.7e-14 * * *	two.sided

timesdata <- expObsTime(rrAnalysisdata,title="Logistic: Expected vs Observed")

pander::pander(timesdata)

	O:Low Risk < 0.445	O:0.445 <= Risk < 0.558	O:High Risk >= 0.558	E:Low Risk < 0.445	E:0.445 <= Risk < 0.558	E:High Risk >= 0.558
1Q	3.44	2.00	1.06	5.37	3.39	1.41
Median	6.89	4.38	2.25	6.50	3.68	2.02
3Q	9.43	7.98	4.70	8.21	4.06	2.53

Training Report

pander::pander(t(rrAnalysisLogTrain$keyPoints),caption="Threshold values")

Threshold values
	@:0.9	@:0.8	@MAX_BACC	@MAX_RR	@SPE100	p(0.5)
Thr	0.558	0.445	0.380	0.244	1.89e-01	0.500
RR	1.755	1.746	1.818	2.197	9.73e+01	1.762
RR_LCI	1.648	1.634	1.686	1.767	2.01e-01	1.654
RR_UCI	1.868	1.865	1.959	2.731	4.71e+04	1.877
SEN	0.323	0.472	0.626	0.958	1.00e+00	0.381
SPE	0.900	0.800	0.673	0.134	1.30e-02	0.867
BACC	0.611	0.636	0.649	0.546	5.06e-01	0.624
NetBenefit	0.102	0.161	0.220	0.351	3.96e-01	0.129

pander::pander(t(rrAnalysisLogTrain$OERatio$estimate),caption="O/E Ratio")

O/E Ratio
O/E	Low	Upper	p.value
0.817	0.777	0.859	4.9e-16

pander::pander(rrAnalysisLogTrain$c.index,caption="C. Index")

C Index: 0.678
Dxy: 0.357
S.D.: 0.014
n: 2982
missing: 0
uncensored: 1518
Relevant Pairs: 6184528
Concordant: 4196033
Uncertain: 2703838
cstatCI:

mean.C Index median lower upper

0.678 0.678 0.665 0.692

pander::pander(t(rrAnalysisLogTrain$ROCAnalysis$aucs),caption="ROC AUC")

ROC AUC
est	lower	upper
0.695	0.676	0.714

pander::pander((rrAnalysisLogTrain$ROCAnalysis$sensitivity),caption="Sensitivity")

Sensitivity
est	lower	upper
0.323	0.299	0.347

pander::pander((rrAnalysisLogTrain$ROCAnalysis$specificity),caption="Specificity")

Specificity
est	lower	upper
0.9	0.883	0.915

pander::pander(t(rrAnalysisLogTrain$thr_atP),caption="Probability Thresholds")

Probability Thresholds
90%	80%
0.558	0.445

pander::pander(rrAnalysisLogTrain$surdif,caption="Logrank test")

Logrank test Chisq = 531.385386 on 2 degrees of freedom, p = 0.000000
	N	Observed	Expected	(O-E)^2/E	(O-E)^2/V
class=0	1973	802	1142	101.1	413.6
class=1	372	226	173	16.1	18.2
class=2	637	490	203	405.1	474.8

Results on the validation set using Logistic model

pre <- predict(mlog,dataBrestCancerTest)
rdata <- cbind(dataBrestCancerTest$status,pre)

rrAnalysisLogTest <- RRPlot(rdata,atThr=rrAnalysisTrain$thr_atP,
                     timetoEvent=dataBrestCancerTest$time,
                     title="Logistic Test: Breast Cancer",
                     ysurvlim=c(0.00,1.0),
                     riskTimeInterval=5)

par(op)

By Risk Categories

obsexp <- rrAnalysisLogTest$OERatio$atThrEstimates

expObs(obsexp,"Logistic Test: Expected vs. Observed")

pander::pander(obsexp)

	Observed	L.CI	H.CI	Expected	pvalue
Total	299	266.1	334.9	291.7	0.6605
low	73	57.2	91.8	59.9	0.0926
90%	48	35.4	63.6	44.0	0.5458
80%	178	152.8	206.2	185.2	0.6328

Time to Event Analysis

rrAnalysisdata <- rrAnalysisLogTest

pander::pander(wilcox.test(rrAnalysisdata$timetoEventData$eTime,rrAnalysisdata$timetoEventData$expectedTime,paired = TRUE))

Wilcoxon signed rank test with continuity correction: `rrAnalysisdata$timetoEventData$eTime` and `rrAnalysisdata$timetoEventData$expectedTime`
Test statistic	P value	Alternative hypothesis
43507	1.85e-46 * * *	two.sided

highrisk <- rrAnalysisdata$timetoEventData$class == 2
pander::pander(wilcox.test(rrAnalysisdata$timetoEventData$eTime[highrisk],rrAnalysisdata$timetoEventData$expectedTime[highrisk],paired = TRUE))

Wilcoxon signed rank test with continuity correction: `rrAnalysisdata$timetoEventData$eTime[highrisk]` and `rrAnalysisdata$timetoEventData$expectedTime[highrisk]`
Test statistic	P value	Alternative hypothesis
24472	0.36	two.sided

timesdata <- expObsTime(rrAnalysisdata,title="Logistic Test: Expected vs Observed")

pander::pander(timesdata)

	O:Low Risk < 0.364	O:0.364 <= Risk < 0.435	O:High Risk >= 0.435	E:Low Risk < 0.364	E:0.364 <= Risk < 0.435	E:High Risk >= 0.435
1Q	2.10	1.75	1.30	6.11	4.66	1.83
Median	3.67	3.24	2.27	6.59	4.93	2.73
3Q	4.98	4.72	3.93	7.77	5.29	3.68

Validation Report

pander::pander(t(rrAnalysisLogTest$keyPoints),caption="Threshold values")

Threshold values
	@:0.435	@:0.364	@MAX_BACC	@MAX_RR	@SPE100	p(0.5)
Thr	0.4352	0.364	0.443	0.290	2.17e-01	0.5007
RR	1.6630	1.601	1.781	2.741	2.64e+01	1.7135
RR_LCI	1.3951	1.294	1.495	1.659	5.65e-02	1.4549
RR_UCI	1.9822	1.980	2.121	4.529	1.23e+04	2.0180
SEN	0.5953	0.759	0.592	0.957	1.00e+00	0.4749
SPE	0.6279	0.411	0.661	0.163	1.55e-02	0.7545
BACC	0.6116	0.585	0.627	0.560	5.08e-01	0.6147
NetBenefit	0.0978	0.141	0.106	0.224	2.82e-01	0.0682

pander::pander(t(rrAnalysisLogTest$OERatio$estimate),caption="O/E Ratio")

O/E Ratio
O/E	Low	Upper	p.value
1.03	0.912	1.15	0.66

pander::pander(rrAnalysisLogTest$c.index,caption="C. Index")

C Index: 0.665
Dxy: 0.33
S.D.: 0.0309
n: 686
missing: 0
uncensored: 299
Relevant Pairs: 266144
Concordant: 177019
Uncertain: 203702
cstatCI:

mean.C Index median lower upper

0.665 0.665 0.634 0.698

pander::pander(t(rrAnalysisLogTest$ROCAnalysis$aucs),caption="ROC AUC")

ROC AUC
est	lower	upper
0.666	0.625	0.707

pander::pander((rrAnalysisLogTest$ROCAnalysis$sensitivity),caption="Sensitivity")

Sensitivity
est	lower	upper
0.595	0.537	0.651

pander::pander((rrAnalysisLogTest$ROCAnalysis$specificity),caption="Specificity")

Specificity
est	lower	upper
0.628	0.578	0.676

pander::pander(t(rrAnalysisLogTest$thr_atP),caption="Probability Thresholds")

Probability Thresholds
90%	80%
0.435	0.364

pander::pander(rrAnalysisLogTest$surdif,caption="Logrank test")

Logrank test Chisq = 50.463152 on 2 degrees of freedom, p = 0.000000
	N	Observed	Expected	(O-E)^2/E	(O-E)^2/V
class=0	232	73	118.3	17.34	28.90
class=1	132	48	61.9	3.12	3.94
class=2	322	178	118.8	29.46	49.41

Logistic Model Poisson Calibration

riskdata <- cbind(dataBrestCancerTrain$status,predict(mlog,dataBrestCancerTrain,type="prob"),dataBrestCancerTrain$time)
calprob <- CalibrationProbPoissonRisk(riskdata)

( 7.468681 , 29202.06 , 1.115209 , 1518 , 1836.231 )

pander::pander(c(h0=calprob$h0,
                 Gain=calprob$hazardGain,
                 DeltaTime=calprob$timeInterval),
               caption="Logistic Calibration Parameters")

h0	Gain	DeltaTime
0.689	1.33	7.47

timeinterval <- calprob$timeInterval;
gain <- calprob$hazardGain

rdata <- cbind(dataBrestCancerTrain$status,calprob$prob)


rrAnalysisLogCalTrain <- RRPlot(rdata,atRate=c(0.90,0.80),
                     timetoEvent=dataBrestCancerTrain$time,
                     title="Cal. Logistic Train: Breast Cancer",
                     ysurvlim=c(0.00,1.0),
                     riskTimeInterval=timeinterval)

par(op)

By Risk Categories

obsexp <- rrAnalysisLogCalTrain$OERatio$atThrEstimates

expObs(obsexp,"Logistic Cal: Expected vs. Observed")

pander::pander(obsexp)

	Observed	L.CI	H.CI	Expected	pvalue
Total	1518	1443	1596	1653	0.000821
low	802	747	859	890	0.003165
90%	226	197	257	229	0.868816
80%	490	448	535	530	0.082339

Time to Event Analysis

rrAnalysisdata <- rrAnalysisLogCalTrain

pander::pander(wilcox.test(rrAnalysisdata$timetoEventData$eTime,rrAnalysisdata$timetoEventData$expectedTime,paired = TRUE))

Wilcoxon signed rank test with continuity correction: `rrAnalysisdata$timetoEventData$eTime` and `rrAnalysisdata$timetoEventData$expectedTime`
Test statistic	P value	Alternative hypothesis
1927275	2.85e-10 * * *	two.sided

highrisk <- rrAnalysisdata$timetoEventData$class == 2
pander::pander(wilcox.test(rrAnalysisdata$timetoEventData$eTime[highrisk],rrAnalysisdata$timetoEventData$expectedTime[highrisk],paired = TRUE))

Wilcoxon signed rank test with continuity correction: `rrAnalysisdata$timetoEventData$eTime[highrisk]` and `rrAnalysisdata$timetoEventData$expectedTime[highrisk]`
Test statistic	P value	Alternative hypothesis
125667	2.23e-07 * * *	two.sided

timesdata <- expObsTime(rrAnalysisdata,title="Logistic Cal: Expected vs Observed")

pander::pander(timesdata)

	O:Low Risk < 0.543	O:0.543 <= Risk < 0.662	O:High Risk >= 0.662	E:Low Risk < 0.543	E:0.543 <= Risk < 0.662	E:High Risk >= 0.662
1Q	3.44	2.00	1.06	6.03	3.81	1.59
Median	6.89	4.38	2.25	7.30	4.13	2.26
3Q	9.43	7.98	4.70	9.22	4.56	2.84

Report of the calibrated logistic: training

pander::pander(t(rrAnalysisLogCalTrain$keyPoints),caption="Threshold values")

Threshold values
	@:0.9	@:0.8	@MAX_BACC	@MAX_RR	@SPE100	p(0.5)
Thr	0.6626	0.543	0.470	0.310	2.43e-01	0.500
RR	1.7547	1.746	1.818	2.197	9.73e+01	1.770
RR_LCI	1.6484	1.634	1.686	1.767	2.01e-01	1.647
RR_UCI	1.8679	1.865	1.959	2.731	4.71e+04	1.902
SEN	0.3228	0.472	0.626	0.958	1.00e+00	0.580
SPE	0.8996	0.800	0.673	0.134	1.30e-02	0.709
BACC	0.6112	0.636	0.649	0.546	5.06e-01	0.644
NetBenefit	0.0675	0.123	0.176	0.297	3.54e-01	0.152

pander::pander(t(rrAnalysisLogCalTrain$OERatio$estimate),caption="O/E Ratio")

O/E Ratio
O/E	Low	Upper	p.value
0.918	0.873	0.966	0.000821

pander::pander(rrAnalysisLogCalTrain$c.index,caption="C. Index")

C Index: 0.678
Dxy: 0.357
S.D.: 0.014
n: 2982
missing: 0
uncensored: 1518
Relevant Pairs: 6184528
Concordant: 4196034
Uncertain: 2703838
cstatCI:

mean.C Index median lower upper

0.678 0.679 0.665 0.692

pander::pander(t(rrAnalysisLogCalTrain$ROCAnalysis$aucs),caption="ROC AUC")

ROC AUC
est	lower	upper
0.695	0.676	0.714

pander::pander((rrAnalysisLogCalTrain$ROCAnalysis$sensitivity),caption="Sensitivity")

Sensitivity
est	lower	upper
0.323	0.299	0.347

pander::pander((rrAnalysisLogCalTrain$ROCAnalysis$specificity),caption="Specificity")

Specificity
est	lower	upper
0.9	0.883	0.915

pander::pander(t(rrAnalysisLogCalTrain$thr_atP),caption="Probability Thresholds")

Probability Thresholds
90%	80%
0.662	0.543

pander::pander(rrAnalysisLogCalTrain$surdif,caption="Logrank test")

Logrank test Chisq = 531.385386 on 2 degrees of freedom, p = 0.000000
	N	Observed	Expected	(O-E)^2/E	(O-E)^2/V
class=0	1973	802	1142	101.1	413.6
class=1	372	226	173	16.1	18.2
class=2	637	490	203	405.1	474.8

probLog <- predict(mlog,dataBrestCancerTest)
aprob <- adjustProb(probLog,gain)

rdata <- cbind(dataBrestCancerTest$status,aprob)
rrAnalysisTestLogistic <- RRPlot(rdata,atThr=rrAnalysisTrain$thr_atP,
                     timetoEvent=dataBrestCancerTest$time,
                     title="Cal. Logistic Test: Breast Cancer",
                     ysurvlim=c(0.00,1.0),
                     riskTimeInterval=timeinterval)

par(op)

By Risk Categories

obsexp <- rrAnalysisTestLogistic$OERatio$atThrEstimates

expObs(obsexp,"Logistic Test: Expected vs. Observed")

pander::pander(obsexp)

	Observed	L.CI	H.CI	Expected	pvalue
Total	299	266.07	334.9	259.6	0.01679
low	13	6.92	22.2	16.3	0.53401
90%	44	31.97	59.1	29.0	0.00891
80%	242	212.47	274.5	212.0	0.04263

Time to Event Analysis

rrAnalysisdata <- rrAnalysisTestLogistic

pander::pander(wilcox.test(rrAnalysisdata$timetoEventData$eTime,rrAnalysisdata$timetoEventData$expectedTime,paired = TRUE))

Wilcoxon signed rank test with continuity correction: `rrAnalysisdata$timetoEventData$eTime` and `rrAnalysisdata$timetoEventData$expectedTime`
Test statistic	P value	Alternative hypothesis
28930	1.07e-65 * * *	two.sided

highrisk <- rrAnalysisdata$timetoEventData$class == 2
pander::pander(wilcox.test(rrAnalysisdata$timetoEventData$eTime[highrisk],rrAnalysisdata$timetoEventData$expectedTime[highrisk],paired = TRUE))

Wilcoxon signed rank test with continuity correction: `rrAnalysisdata$timetoEventData$eTime[highrisk]` and `rrAnalysisdata$timetoEventData$expectedTime[highrisk]`
Test statistic	P value	Alternative hypothesis
26784	3.21e-26 * * *	two.sided

timesdata <- expObsTime(rrAnalysisdata,title="Logistic Test: Expected vs Observed")

pander::pander(timesdata)

	O:Low Risk < 0.364	O:0.364 <= Risk < 0.435	O:High Risk >= 0.435	E:Low Risk < 0.364	E:0.364 <= Risk < 0.435	E:High Risk >= 0.435
1Q	2.65	2.08	1.45	8.78	6.96	2.59
Median	4.08	3.67	2.45	9.31	7.18	4.17
3Q	5.51	4.66	4.35	10.21	7.67	5.35

Report of the calibrated validation

pander::pander(t(rrAnalysisTestLogistic$keyPoints),caption="Threshold values")

Threshold values
	@:0.435	@:0.364	@MAX_BACC	@MAX_RR	@SPE100	p(0.5)
Thr	0.4360	0.364	0.5401	0.366	2.77e-01	0.4999
RR	1.7104	2.660	1.7806	2.741	2.64e+01	1.7914
RR_LCI	1.3502	1.612	1.4950	1.659	5.65e-02	1.4791
RR_UCI	2.1667	4.389	2.1207	4.529	1.23e+04	2.1696
SEN	0.8094	0.957	0.5920	0.957	1.00e+00	0.6823
SPE	0.3618	0.158	0.6615	0.163	1.55e-02	0.5607
BACC	0.5856	0.557	0.6267	0.560	5.08e-01	0.6215
NetBenefit	0.0745	0.145	0.0338	0.145	2.23e-01	0.0497

pander::pander(t(rrAnalysisTestLogistic$OERatio$estimate),caption="O/E Ratio")

O/E Ratio
O/E	Low	Upper	p.value
1.15	1.03	1.29	0.0168

pander::pander(rrAnalysisTestLogistic$c.index,caption="C. Index")

C Index: 0.665
Dxy: 0.33
S.D.: 0.0309
n: 686
missing: 0
uncensored: 299
Relevant Pairs: 266144
Concordant: 177019
Uncertain: 203702
cstatCI:

mean.C Index median lower upper

0.665 0.665 0.634 0.694

pander::pander(t(rrAnalysisTestLogistic$ROCAnalysis$aucs),caption="ROC AUC")

ROC AUC
est	lower	upper
0.666	0.625	0.707

pander::pander((rrAnalysisTestLogistic$ROCAnalysis$sensitivity),caption="Sensitivity")

Sensitivity
est	lower	upper
0.809	0.76	0.852

pander::pander((rrAnalysisTestLogistic$ROCAnalysis$specificity),caption="Specificity")

Specificity
est	lower	upper
0.362	0.314	0.412

pander::pander(t(rrAnalysisTestLogistic$thr_atP),caption="Probability Thresholds")

Probability Thresholds
90%	80%
0.435	0.364

pander::pander(rrAnalysisTestLogistic$surdif,caption="Logrank test")

Logrank test Chisq = 37.907344 on 2 degrees of freedom, p = 0.000000
	N	Observed	Expected	(O-E)^2/E	(O-E)^2/V
class=0	75	13	44.1	21.93	26.01
class=1	122	44	59.9	4.22	5.29
class=2	489	242	195.0	11.32	32.86

Comparing the COX and Logistic Models on the Independent Data

pander::pander(t(rrCoxTestAnalysis$OAcum95ci))

mean	50%	2.5%	97.5%
1.24	1.24	1.24	1.25

pander::pander(t(rrAnalysisTestLogistic$OAcum95ci))

mean	50%	2.5%	97.5%
0.931	0.931	0.928	0.934

pander::pander(t(rrCoxTestAnalysis$OE95ci))

mean	50%	2.5%	97.5%
1.17	1.17	1.14	1.2

pander::pander(t(rrAnalysisTestLogistic$OE95ci))

mean	50%	2.5%	97.5%
0.977	0.977	0.951	1

maxobs <- sum(dataBrestCancerTest$status)

par(mfrow=c(1,2),cex=0.75)

plot(rrCoxTestAnalysis$CumulativeOvs[,1:2],type="l",lty=1,
     main="Cumulative Probability",
     xlab="Observed",
     ylab="Cumulative Probability",
     ylim=c(0,maxobs),
     xlim=c(0,maxobs))
lines(rrAnalysisTestLogistic$CumulativeOvs[,1:2],lty=2,col="red")
lines(x=c(0,maxobs),y=c(0,maxobs),lty=3,col="gray")
legend("topleft",legend = c("Cox","Logistic","Ideal"),
       col=c("black","red","gray"),
       lty=c(1,2,3),
       cex=0.75
)

dxcox <- rrCoxTestAnalysis$CumulativeOvs$Cumulative-
       rrCoxTestAnalysis$CumulativeOvs$Observed

dxlogit <- rrAnalysisTestLogistic$CumulativeOvs$Cumulative-
       rrAnalysisTestLogistic$CumulativeOvs$Observed

miny <- min(c(dxcox,dxlogit))
maxy <- max(c(dxcox,dxlogit))
plot(rrCoxTestAnalysis$CumulativeOvs$Observed,
     dxcox,
     main="Cumulative Risk Difference",
     xlab="Observed",
     ylab="Cumulative Risk - Observed",
     type="l",
     ylim=c(miny,maxy),
     lty=1)
lines(rrAnalysisTestLogistic$CumulativeOvs$Observed,
     dxlogit,
     lty=2,
     col="red")
legend("topleft",legend = c("Cox","Logistic"),
       col=c("black","red"),
       lty=c(1,2),
       cex=0.75
)

plot(rrCoxTestAnalysis$OEData[,2:3],type="l",lty=1,
     main="Expected over Time",
     xlab="Observed",
     ylab="Expected",
     ylim=c(0,maxobs),
     xlim=c(0,maxobs))
lines(rrAnalysisTestLogistic$OEData[,2:3],lty=2,col="red")
lines(x=c(0,maxobs),y=c(0,maxobs),lty=3,col="gray")
legend("topleft",legend = c("Cox","Logistic","Ideal"),
       col=c("black","red","gray"),
       lty=c(1,2,3),
       cex=0.75
)

coxdif <- rrCoxTestAnalysis$OEData$Expected-
       rrCoxTestAnalysis$OEData$Observed

logdif <- rrAnalysisTestLogistic$OEData$Expected-
       rrAnalysisTestLogistic$OEData$Observed

miny <- min(c(coxdif,logdif))
maxy <- max(c(coxdif,logdif))

plot(rrCoxTestAnalysis$OEData$Observed,
     coxdif,
      ylim=c(miny,maxy),

     main="Expected vs Observed Difference",
     xlab="Observed",
     ylab="Cumulative - Observed",
     type="l",
     lty=1)
lines(rrAnalysisTestLogistic$OEData$Observed,
     logdif,
     lty=2,col="red")

legend("bottomleft",legend = c("Cox","Logistic"),
       col=c("black","red"),
       lty=c(1,2),
       cex=0.75
)

par(op)

Risk-Evaluation: Breast Cancer Royston-Altman

Jose Tamez

2023-05-22

Evaluation of RISK survival models

The libraries

Breast Cancer Royston-Altman data

data(gbsg, package=“survival”) and data(rotterdam, package=“survival”)

data(gbsg, package=“survival”) data conditioning

Cox Modeling

Cox Model Performance

The evaluation of the raw Cox model with RRPlot()

By Risk Categories

Time to Event Analysis

Uncalibrated Performance Report

Cox Calibration

The RRplot() of the calibrated model

By Risk Categories

Time to Event Analysis

Calibrated Train Performance

Performance on the external data set

External Data Report

Calibrating the index on the test data

By Risk Categories

Time to Event Analysis

After Calibration Report

Logistic Model

Logistic Model Performance

By Risk Categories

Time to Event Analysis

Training Report

Results on the validation set using Logistic model

By Risk Categories

Time to Event Analysis

Validation Report

Logistic Model Poisson Calibration

By Risk Categories

Time to Event Analysis

Report of the calibrated logistic: training

By Risk Categories

Time to Event Analysis

Report of the calibrated validation

Comparing the COX and Logistic Models on the Independent Data