examples/ch6.R

In tagteam/MedicalRiskPredictionModels: R-codes for book: Medical Risk Prediction Models by Gerds & kattan

library(MedicalRiskPredictionModels)
prepareExamples()

# Chunk1
# conventional model
fit1 <- lrm(ohss~ant.foll+cyclelen+smoking+age,data=ivftrain)
# experimental model
fit2 <- lrm(ohss~rcs(ant.foll,3)*smoking+cyclelen+age+fsh+bmi+ovolume,data=ivftrain,penalty=10)
# head to head comparison in test dataset
x <- Score(list("Conventional"=fit1,"Experimental"=fit2),
           data=ivftest, formula=ohss~1, summary=c("risks","ipa"),
           plots=c("roc","cal"))
# scatterplot showing predicted risks of the rival models  
plotRisk(x,
         col=c("gray22","black"),
         xlab="Conventional model: risk of OHSS",
         ylab="Experimental model: risk of OHSS")

# Chunk2
summary(x,what="score")

# Chunk3
summary(x,what="contrasts")

# Chunk4
plotCalibration(x,auc.in.legend=1,brier.in.legend=1)

# Chunk5
plotROC(x)

# Chunk6
dd <- datadist(octrain.cc)
options(datadist="dd")
fit1 <- cph(Surv(survtime,survstatus)~rcs(age,3)+rcs(tumorthickness,3)+gender+tobacco+deep.invasion+site+race+x.posnodes+tumormaxdimension+vascular.invasion,data=octrain.cc,x=TRUE,surv=TRUE)
surv <- Survival(fit1)
plot(nomogram(fit1,fun=list(function(x) 1-surv(60, x),
                           function(x) 1-surv(120, x)),
              funlabel=c("5 year risk","10 year risk")))

# Chunk7
set.seed(1972)
fit2 <- rfsrc(Surv(survtime,survstatus)~ age+tumorthickness+gender+tobacco+deep.invasion+site+race+x.posnodes+tumormaxdimension+vascular.invasion,data=octrain.cc)

# Chunk8
fit1 <- cph(Surv(survtime,survstatus)~rcs(age,3)+tumorthickness+gender+tobacco+deep.invasion+site+race+x.posnodes+tumormaxdimension+vascular.invasion,
            data=octrain.cc, x=TRUE, y=TRUE, surv=TRUE)
set.seed(1972)
fit2 <- rfsrc(Surv(survtime,survstatus)~ age+tumorthickness+gender+tobacco+deep.invasion+site+race+x.posnodes+tumormaxdimension+vascular.invasion,data=octrain.cc)
x <- Score(list("Conventional"=fit1,"Experimental"=fit2),
           data=octest.cc,
           formula=Surv(survtime,survstatus)~1,
           times=c(12,60,120,180),
           summary=c("risks","IPA"),
           plots=c("cali","roc"))
summary(x,what="score")[[1]]

# Chunk9
summary(x,what="contrasts")

# Chunk10
library(bootstrap)
set.seed(8)
# time of death or censored in the complete cases
# of the training set
boot.max <- bootstrap(octrain.cc$survtime, 
                      100000, # number of bootstrap samples
                      theta=function(x){max(x)}) # maximum time
# minimum of the bootstrap maxima, divided by 12
# to convert from months to years
min(boot.max$thetastar)/12

# Chunk11
x <- Score(list("Conventional"=fit1,"Experimental"=fit2),
           data=octest.cc,
           formula=Surv(survtime,survstatus)~1,
           times=seq(12,180,6))
plotBrier(x,
          legend.x="bottomleft",
          xlab="Years from surgery",
          axis1.at=seq(0,180,12),
          axis1.labels=seq(0,15,1))

# Chunk12
x <- Score(list("Conventional"=fit1,"Experimental"=fit2),
           data=octest.cc,
           formula=Surv(survtime,survstatus)~1,
           times=seq(12,180,6))
plotAUC(x, legend.x="bottomleft",
        xlab="Years from surgery",
        col=c("#E69F00", "#56B4E9"),
        axis1.at=seq(0,180,12),
        axis1.labels=seq(0,15,1))

# Chunk13
fit1 <- FGR(Hist(asprogtime,asprog)~age+psa+ct1+diaggs+ppb5,data=astrain,cause="progression")
fit1

# Chunk14
fit2 <- CSC(list(Hist(asprogtime,asprog)~age+psa+ct1+diaggs+ppb5,Hist(asprogtime,asprog)~age),data=astrain,cause="progression")
fit2

# Chunk15
x <- Score(list("Conventional"=fit1,"Experimental"=fit2),
           data=astest,
           formula=Hist(asprogtime,asprog)~1,
           times=3,summary="risk",
           cause="progression")
plotRisk(x,
         times=3,
         plot.main="Risk of progression within 3-years",
         xlab="Conventional model prediction",
         ylab="Experimental model prediction")

# Chunk16
x <- Score(list("Conventional"=fit1,"Experimental"=fit2),
           data=astest,
           formula=Hist(asprogtime,asprog)~1,
           times=3,
           cause="progression")
summary(x,what="score")

# Chunk17
summary(x,what="contrasts")

# Chunk18
x <- Score(list("Conventional"=fit1,"Experimental"=fit2),
           data=astest, formula=Hist(asprogtime,asprog)~1,
           cause="progression", times=seq(.5,5,.5))
plotBrier(x)

# Chunk19
x <- Score(list("Conventional"=fit1,"Experimental"=fit2),
           data=astest, formula=Hist(asprogtime,asprog)~1,
           cause="progression", times=seq(.5,5,.5))
plotAUC(x)

# Chunk20
set.seed(1972)
fit1 <- rfsrc(Surv(survtime,survstatus)~ age+tumorthickness+gender+tobacco+deep.invasion+site+race+x.posnodes+tumormaxdimension+vascular.invasion,data=octrain.cc)
set.seed(1972)
fit2 <- rfsrc(Surv(survtime,survstatus)~ age+tumorthickness+gender+tobacco+deep.invasion+site+race+x.posnodes+tumormaxdimension+vascular.invasion+Grade,data=octrain.cc)
x <- Score(list("Conventional"=fit1,"New marker"=fit2),
           data=octest.cc,
           formula=Surv(survtime,survstatus)~1,
           times=60,
           summary=c("risks","IPA"))
plotRisk(x,times=60,preclipse.shade=1,legend.ncol=3,legend.x=0,legend.y=1.1,legend.xpd=NA)
mtext("Predicted risk of 5-year mortality",line=2)

# Chunk21
fit1 <- CSC(list(Hist(asprogtime,asprog)~age+psa+ct1+diaggs+ppb5,Hist(asprogtime,asprog)~age),data=astrain,cause="progression")
fit2 <- CSC(list(Hist(asprogtime,asprog)~age+psa+ct1+diaggs+ppb5+erg.status,Hist(asprogtime,asprog)~age+erg.status),data=astrain,cause="progression")
x <- Score(list("Conventional"=fit1,"New marker"=fit2),formula=Hist(asprogtime,asprog)~1,cause="progression",data=astest,times=3,summary="risks")
plotRisk(x,times=3,xlim=c(0,1),ylim=c(0,1),preclipse.shade=1,legend.ncol=3,legend.x=0,legend.y=1,legend.xpd=NA)
mtext("Risk of progression within 3-years",side=3,line=2.3)