Case study on survival data
Chapter 23 illustrates survival modeling with the SMART data (Principal investigator: Prof Yolanda van der Graaf). The R script file is shown below, which was constructed with support from Martijn Gondrie (Utrecht).
# EXAMPLE_SYNTAX paper prof. dr. Y. van der Graaf + Prof. dr. E. Steyerberg
# DATE: 29May2009
library(Design)
library(mice)
library(foreign)
# Importation of SMART-data-file with correct variables and endpoints
# Variables are already truncated to 1th en 99th percentiles
SMART <- read.spss('F:/Book/Applications/Graaf/SMARTst.sav',use.value.labels=F, to.data.frame=T)
# Original SMARTfile not truncated
SMARTo <- read.spss('F:/Book/Applications/Graaf/SMARTs.sav',use.value.labels=F, to.data.frame=T)
#-------------------------------------------------------------------------------
## PAPER STEP 1: Preliminary#
# See the first row of the data-file; also to get an idea of which variables are included
SMART[1,]
# Retrieve the dimensions of your file: (How many rows(cases) and how many colums (variables))
dim(SMART)
# Retrieve basic information per variable (e.g: mean, highest value, lowest value, missings etc.)
describe(SMART)
# Generate an extra window for 1 plot
par(mfrow=c(1,1))
# Selecting 18 variables out of 29 in the loaded SMART-file and immediatly retrieve info again.
SMARTm <- SMART[,Cs(SYSTBP, DIASTBP, SYSTH, DIASTH,
CHOL,HDL,LDL,TRIG, HOMOC, GLUT, CREAT, IMT, ALBUMIN,
STENOSIS, DIABETES, SMOKING, PACKYRS, ALCOHOL)]
describe(SMARTm)
# With the following function you determine the fraction of NA's (missing values)
na.patterns <- naclus(SMARTm)
# Plot the 18 variables you just selected in a graph: Fig 23.2.
plot(na.patterns, ylab="Fraction of NAs in common")
# A different approach to show the NA's is the following:
# First you generate a screen in which two graphs can be shown: Fig 23.1.
par(mfrow=c(2,2))
naplot(na.patterns, which=c('na per var'))
naplot(na.patterns, which=c('na per obs'))
#-------------------------------------------------------------------------------
## PAPER STEP 2: Coding of predictors. ##
# Start of IMT-variable-example just as in the paper
# Showing the maximum value of IMT in the truncated file.
max(SMART$IMT,na.rm=T)
# Show how many times this maximum value is exceeded in the original file
sum(SMARTo$IMTO>1.83,na.rm=T)
# Same, but now for the minimum value
min(SMART$IMT,na.rm=T)
sum(SMARTo$IMTO<.47,na.rm=T)
# show relation between IMT in both the not-truncated file and the truncated file with the cox proportional hazard function
cph(Surv(TEVENT,EVENT) ~ IMT, data=SMART)
cph(Surv(TEVENT,EVENT) ~ IMTO, data=SMARTo)
# again an extra screen
par(mfrow=c(1,2),mar=c(5,5,3,1))
# make an new matrix with the the original and truncated IMT variable. Then name the two colums
IMTdata <- as.matrix(cbind(SMARTo$IMTO, SMART$IMT))
IMTlist <- split(IMTdata, col(IMTdata))
names(IMTlist) <- Cs(Original, Truncated)
# make a boxplot of this matirx to show the difference of the data when it's truncated
boxplot(IMTlist, ylab= "IMT (mm)", cex.lab=1.2) # Fig 23.4
# Now the effect on outcome is shown with the cph-function in three ways
#to show how the variable which fit the model best:
#1: Not truncated with restricted cubic spline with 5 knots
dd <- datadist(SMARTo)
options(datadist="dd")
IMTOfit <- cph(Surv(TEVENT,EVENT) ~ rcs(IMTO,5), data=SMARTo)
plot(IMTOfit, xlab="IMT (mm)", lty=2, lwd=2, conf.int=F)
#2: Truncated with restricted cubic spline with 5 knots
dd <- datadist(SMART)
options(datadist="dd")
IMTfit <- cph(Surv(TEVENT,EVENT) ~ rcs(IMT,5), data=SMART)
plot(IMTfit, xlab="", add=T, lty=3, lwd=2,conf.int=F)
#3: Truncated linear
IMTfitlin <- cph(Surv(TEVENT,EVENT) ~ IMT, data=SMART)
plot(IMTfitlin, xlab="", add=T, lty=1, lwd=2, conf.int=F)
# Putting text in the graphs
text(x=2,y=.8,"Original,\nspline",adj=0)
text(x=1.8,y=1.1,"Truncated,\nspline",adj=0)
text(x=1.4,y=.8,"Truncated,\nlinear",adj=1)
## End IMT illustration, Fig 23.4
# Start of Age-example just as in the paper
# Example transformation of age to get best fit.
# With cph-function the X2 and dof are measured when the variable in transformed.
dd <- datadist(SMART)
options(datadist="dd")
par(mfrow=c(2,2), mar=c(5,5,3,1))
# age related to outcome: age linear --> X2=97
AGEfit <- cph(Surv(TEVENT,EVENT) ~ AGE, data=SMART)
plot(AGEfit, xlab="Age (years)", lty=1, lwd=2,conf.int=T, ylim=c(-1.5, 1.5), xlim=c(20,80))
anova(AGEfit)
# age related to outcome: age + age squared --> X2= 125
AGEfit2 <- cph(Surv(TEVENT,EVENT) ~ pol(AGE,2), data=SMART) # age + age squared
anova(AGEfit2)
plot(AGEfit2, xlab="Age (years)", add=T, lty=2, lwd=2,conf.int=T)
title(substitute(paste("Age, 1 df, ", chi^2," ", 97,
"; Age+",Age^2,", 2 df, ", chi^2," ", 125,sep="")), cex.main=1)
# Age related to outcome: 125truncated under 50 / 55 and then linear --> X2= 119
AGEfit3 <- cph(Surv(TEVENT,EVENT) ~ ifelse(AGE>55, (AGE-55),0), data=SMART) # age linear
anova(AGEfit3)
plot(AGEfit3, xlab="Age (years)", lty=1, lwd=2,conf.int=T, ylim=c(-1.5, 1.5), xlim=c(20,80))
# Age related to outcome: 125truncated under 50 / 55 and then squared --> X2= 130
AGEfit4 <- cph(Surv(TEVENT,EVENT) ~ ifelse(AGE>50, (AGE-50)^2,0), data=SMART) # age + age squared
anova(AGEfit4)
plot(AGEfit4, xlab="Age (years)", add=T, lty=2, lwd=2,conf.int=T)
title(substitute(paste("Age>55, 1 df, ", chi^2," ", 119,
"; ",(Age>50)^2,", 2 df, ", chi^2," ", 130,sep="")), cex.main=1)
# age related to outcome with Spline and 3 df (4 knots) --> X2=125
AGEfit5 <- cph(Surv(TEVENT,EVENT) ~ rcs(AGE,4), data=SMART)
anovaAge <- anova(AGEfit5)
plot(AGEfit5, xlab="Age (years)", lty=3, lwd=2,conf.int=T,
ylim=c(-1.5, 1.5), xlim=c(20,80))
title(substitute(paste("Spline, 3 df, ", chi^2," ", 125,sep="")), cex.main=1)
# Age related to outcome in quartiles --> X2=93
tapply(SMART$AGE, ifelse(SMART$AGE<50,1,ifelse(SMART$AGE<60,2,ifelse(SMART$AGE<70,3,4))), mean)
SMART$AGEcat <- as.factor(ifelse(SMART$AGE<50,43,ifelse(SMART$AGE<60,55,ifelse(SMART$AGE<70,65,74))))
AGEfit6 <- cph(Surv(TEVENT,EVENT) ~ as.factor(AGEcat), data=SMART)
anova (AGEfit6)
plot(AGEfit6, xlab="Age (years)", lty=3, lwd=2,conf.int=T,
ylim=c(-1.5, 1.5))
title(substitute(paste("Categorized, 3 df, ", chi^2," ", 93,sep="")), cex.main=1)
# Age related to outcome dichotomization --> X2=72
tapply(SMART$AGE, cut2(SMART$AGE,g=2), mean)
AGEfit7 <- cph(Surv(TEVENT,EVENT) ~ ifelse(SMART$AGE<61,51,69), data=SMART)
anova(AGEfit7)
##End AGE-illustration Fig 23.5 and Tab 23.4
# Transformations of other continuous vars: RCS with 3 or 4 knots and normal linear is checked
# BMI
fit <- cph(Surv(TEVENT,EVENT) ~ rcs(BMI,4), data=SMART)
anova(fit)
fit <- cph(Surv(TEVENT,EVENT) ~ rcs(BMI,3), data=SMART)
anova(fit)
# Blood pressure: diastolic
fit <- cph(Surv(TEVENT,EVENT) ~ rcs(DIASTH,4), data=SMART)
anova(fit)
fit <- cph(Surv(TEVENT,EVENT) ~ rcs(DIASTH,3), data=SMART)
anova(fit)
fit <- cph(Surv(TEVENT,EVENT) ~ DIASTH, data=SMART)
anova(fit)
# Blood pressure: Systolic
fit <- cph(Surv(TEVENT,EVENT) ~ rcs(SYSTH,4), data=SMART)
anova(fit)
fit <- cph(Surv(TEVENT,EVENT) ~ rcs(SYSTH,3), data=SMART)
anova(fit)
fit <- cph(Surv(TEVENT,EVENT) ~ SYSTH, data=SMART)
anova(fit)
# Blood pressure: other type of measurement: diastolic
fit <- cph(Surv(TEVENT,EVENT) ~ rcs(DIASTBP,4), data=SMART)
anova(fit)
fit <- cph(Surv(TEVENT,EVENT) ~ rcs(DIASTBP,3), data=SMART)
anova(fit)
plot(fit)
fit <- cph(Surv(TEVENT,EVENT) ~ DIASTBP, data=SMART)
anova(fit)
# Blood pressure: other typeof measurment: systolic
fit <- cph(Surv(TEVENT,EVENT) ~ rcs(SYSTBP,4), data=SMART)
anova(fit)
fit <- cph(Surv(TEVENT,EVENT) ~ rcs(SYSTBP,3), data=SMART)
anova(fit)
fit <- cph(Surv(TEVENT,EVENT) ~ SYSTBP, data=SMART)
anova(fit)
# Lipids: Cholesterol
fit <- cph(Surv(TEVENT,EVENT) ~ rcs(CHOL,4), data=SMART)
anova(fit)
fit <- cph(Surv(TEVENT,EVENT) ~ rcs(CHOL,3), data=SMART)
anova(fit)
fit <- cph(Surv(TEVENT,EVENT) ~ CHOL, data=SMART)
anova(fit)
# Lipids: HDL
fit <- cph(Surv(TEVENT,EVENT) ~ rcs(HDL,4), data=SMART)
anova(fit)
fit <- cph(Surv(TEVENT,EVENT) ~ rcs(HDL,3), data=SMART)
anova(fit)
fit <- cph(Surv(TEVENT,EVENT) ~ HDL, data=SMART)
anova(fit)
# Lipids: LDL
fit <- cph(Surv(TEVENT,EVENT) ~ rcs(LDL,4), data=SMART)
anova(fit)
fit <- cph(Surv(TEVENT,EVENT) ~ rcs(LDL,3), data=SMART)
anova(fit)
fit <- cph(Surv(TEVENT,EVENT) ~ LDL, data=SMART)
anova(fit)
# Lipids: Triglyceriden
fit <- cph(Surv(TEVENT,EVENT) ~ rcs(TRIG,4), data=SMART)
anova(fit)
fit <- cph(Surv(TEVENT,EVENT) ~ rcs(TRIG,3), data=SMART)
anova(fit)
fit <- cph(Surv(TEVENT,EVENT) ~ TRIG, data=SMART)
anova(fit)
# Only HDL seems to have an effect
# Other: Homocysteine
fit <- cph(Surv(TEVENT,EVENT) ~ rcs(HOMOC,4), data=SMART)
anova(fit)
fit <- cph(Surv(TEVENT,EVENT) ~ rcs(HOMOC,3), data=SMART)
anova(fit)
fit <- cph(Surv(TEVENT,EVENT) ~ HOMOC, data=SMART)
anova(fit)
# Other: glucose
fit <- cph(Surv(TEVENT,EVENT) ~ rcs(GLUT,4), data=SMART)
anova(fit)
fit <- cph(Surv(TEVENT,EVENT) ~ rcs(GLUT,3), data=SMART)
anova(fit)
fit <- cph(Surv(TEVENT,EVENT) ~ GLUT, data=SMART)
anova(fit)
# Other: Diabetes
fit <- cph(Surv(TEVENT,EVENT) ~ DIABETES, data=SMART)
anova(fit)
# Other: Diabetes + Glucose
fit <- cph(Surv(TEVENT,EVENT) ~ DIABETES + GLUT, data=SMART)
anova(fit)
#Best is to keep diabetes, no glucose
#Other: Creatinine: Also visualiised in plots, because different transformation is done (logistic)
#Plot 4 graphs for CREAT: Fig 23.6
par(mfrow=c(2,2), mar=c(5,5,3,1))
fit <- cph(Surv(TEVENT,EVENT) ~ CREAT, data=SMART)
anova(fit)
plot(fit, xlab="Creatinin (mmol/l)", lty=1, lwd=2,conf.int=T, ylim=c(-.5, 2.5))
title(substitute(paste("Linear, 1 df, ", chi^2," ", 93,sep="")), cex.main=1.2)
fit <- cph(Surv(TEVENT,EVENT) ~ rcs(CREAT,4), data=SMART)
anova(fit)
plot(fit, xlab="Creatinin (mmol/l)", lty=3, lwd=2,conf.int=T, ylim=c(-.5, 2.5))
title(substitute(paste("Spline, 3 df, ", chi^2," ", 116,sep="")), cex.main=1.2)
fit <- cph(Surv(TEVENT,EVENT) ~ rcs(CREAT,3), data=SMART)
anova(fit)
plot(fit, xlab="Creatinin (mmol/l)", lty=2, lwd=2,conf.int=T, ylim=c(-.5, 2.5))
title(substitute(paste("Spline, 2 df, ", chi^2," ", 99,sep="")), cex.main=1.2)
fit <- cph(Surv(TEVENT,EVENT) ~ log(CREAT), data=SMART)
anova(fit)
plot(fit, xlab="Creatinin (mmol/l)", lty=1, lwd=2,conf.int=T, ylim=c(-.5, 2.5))
title(substitute(paste("Log transformation, 1 df, ", chi^2," ", 131,sep="")), cex.main=1.2)
# Hence: log(CREAT) good choice, Fig 23.6 and Tab 23.4 ##
# Start of History-example just as in the paper
# Score with individual terms, 4 dof
fitsum <- cph(Surv(TEVENT,EVENT) ~ CEREBRAL+CARDIAC+AAA+PERIPH, data=SMART)
anova (fitsum)
# Score with each affected organ seperately in a score, assuming equal weight
SMART$HISTCARD <- SMART$CEREBRAL+SMART$CARDIAC+SMART$AAA+SMART$PERIPH
fitsum2 <- cph(Surv(TEVENT,EVENT) ~ HISTCARD, data=SMART)
anova(fitsum2)
# Score with each affected organ seperately in a score, AAA doubled
SMART$HISTCAR2 <- SMART$CEREBRAL+SMART$CARDIAC+2*SMART$AAA+SMART$PERIPH
describe(SMART$HISTCAR2)
fitsum3 <- cph(Surv(TEVENT,EVENT) ~ HISTCAR2, data=SMART)
anova(fitsum3)
## End history of cardiovascular disease ##
#-------------------------------------------------------------------------------
## PAPER STEP 3 + 4: Model Specification AND model estimation##
# imputation model with outcome, all potential predictors, and auxilliary variables
SMARTM <-aregImpute(~I(TEVENT)+EVENT+SEX+I(AGE)+
SYSTBP+DIASTBP+SYSTH+DIASTH+
DIABETES+CEREBRAL+CARDIAC+AAA+PERIPH+I(HISTCARD)+STENOSIS+
I(LENGTH)+I(WEIGHT)+I(BMI)+
I(CHOL)+I(HDL)+I(LDL)+I(TRIG)+I(HOMOC)+I(GLUT)+I(CREAT)+I(IMT)+
as.factor(ALBUMIN) +as.factor(SMOKING) + I(PACKYRS) + as.factor(ALCOHOL),
n.impute=5, data=SMART)
# Check results of imputation, Fig 23.3
SMARTM
# Use transcan to make 1 singly imputed file.
SMARTc <- SMART
imputed <-impute.transcan(SMARTM, imputation=1, data=SMART, list.out=TRUE, pr=FALSE, check=FALSE)
SMARTc[names(imputed)] <- imputed
SMARTc$HISTCAR2 <- SMARTc$CEREBRAL+SMARTc$CARDIAC+2*SMARTc$AAA+SMARTc$PERIPH
SMARTc$HISTCAR3 <- ifelse(SMARTc$HISTCARD==4, 3, SMARTc$HISTCARD)
SMARTc[1:4,]
#Check and investigate your new file, no missings anymore!!!
describe(SMARTc)
# 3 models: full model, LASSO-model, Backward-step model
#1: Full model combined
dd <- datadist(SMARTc)
options(datadist="dd")
cox01 <-cph(Surv(TEVENT,EVENT) ~ ifelse(AGE>50, (AGE-50)^2,0) + SEX +
as.factor(SMOKING) + as.factor(ALCOHOL)+BMI + SYSTH +
HDL+DIABETES +HISTCAR2 +HOMOC + log(CREAT)+ as.factor(ALBUMIN)+STENOSIS+IMT,
data = SMARTc, x=T, y=T, surv=T)
anova(cox01)
summary(cox01)
cox01
# See Table 23.6: slightly different numbers because of new impuation
# Interactions with sex, 10 df, chi2 14.8
anova(cph(Surv(TEVENT,EVENT) ~ SEX * (ifelse(AGE>55, (AGE-55),0) +
BMI + HDL+DIABETES +HISTCAR2 + log(CREAT)+ as.factor(ALBUMIN)+STENOSIS+IMT),data = SMARTc))
# specific sex*HISTCAR2, 1 df
anova(cph(Surv(TEVENT,EVENT) ~ ifelse(AGE>55, (AGE-55),0) +
BMI + HDL+DIABETES +SEX * HISTCAR2 + log(CREAT)+ as.factor(ALBUMIN)+STENOSIS+IMT,data = SMARTc))
# HISTCAR2 * other predictors, 9 df
anova(cph(Surv(TEVENT,EVENT) ~ SEX + HISTCAR2 * (ifelse(AGE>55, (AGE-55),0) +
BMI + HDL+DIABETES + log(CREAT)+ as.factor(ALBUMIN)+STENOSIS+IMT),data = SMARTc))
# specific sex* +CEREBRAL+CARDIAC+AAA+PERIPH
anova(cph(Surv(TEVENT,EVENT) ~ ifelse(AGE>55, (AGE-55),0) +
BMI + HDL+DIABETES +SEX*(CEREBRAL+CARDIAC+AAA+PERIPH) + log(CREAT)+ as.factor(ALBUMIN)+STENOSIS+IMT,data = SMARTc))
# Plot interaction history * sex
par(mfrow=c(1,2))
# levels(SMARTc$SEX) <- Cs(Male, Female)
dd <- datadist(SMARTc)
options(datadist="dd")
plot(cph(Surv(TEVENT,EVENT) ~ ifelse(AGE>55, (AGE-55),0) +
BMI + HDL+DIABETES +SEX*HISTCAR2 +log(CREAT)+ as.factor(ALBUMIN)+STENOSIS+IMT,data = SMARTc),
HISTCAR2=NA, SEX=NA, xlim=c(1,5.8), xlab="Sumscore previous symptomatic atherosclerosis")
title('sex*sumscore')
text(x=c(5.2,5.2), y=c(0.93,2.56), Cs(Male, Female), adj=0)
# C-statistic of the full model:
rcorr.cens(-cox01$linear.pred,cox01$y)
#2: Stepwise selection
backwardstepwise <- fastbw(cox01, rule="aic", type="individual")
backwardstepwise
cox01.step <- update(cox01, ~ cox01$x[,backwardstepwise$parms.kept])
cox01.step
#3: LASSO-method
## glmpath analysis; 'lasso'
library(glmpath)
memory.size(size = 2048)
## SMART data analysis
SMARTdata <- list(x=cox01$x, time=SMARTc$TEVENT, status=SMARTc$EVENT)
summary(SMARTdata)
## VERY LONG computing time ...
# predictie corrector methode voor berekenen van CPH-model met L1 penalty
SMARTpath <- coxpath(data=SMARTdata)
par(mar=c(5,4,4,6))
plot.coxpath(SMARTpath, cex.lab=1.2)
plot.coxpath(SMARTpath, cex.lab=1.2,type = "aic")
SMARTpath
SMARTpath$aic
round(SMARTpath$b.predictor[SMARTpath$aic==min(SMARTpath$aic),],4)
####################CHECK SMARTHPATH AND CHOOSE ROW IN WHICH HOMOCYST IS DELETED FROM THE MODEL############################
SMARTpath$b.predictor
# model ?? is without homocyst
Lasso.coefs <- SMARTpath$b.predictor[#CHOOSE ROW,]
Lasso.coefs <- Lasso.coefs[Lasso.coefs > 0 | Lasso.coefs < 0]
Lasso.coefs
# Alternsative: Goeman's library
library(penalized)
memory.size(size = 2048)
## new SMART data analysis
opt <- optL1(Surv(SMARTc$TEVENT, SMARTc$EVENT), penalized = cox01$x, startbeta=cox01$coefficients, fold = 10, model='cox', maxlambda1=20, standardize=T)
coefficients(opt$fullfit)
cox01.lasso <- penalized(Surv(SMARTc$TEVENT, SMARTc$EVENT) ~ cox01$x, model='cox', standardize=T, lambda1=18)
coef(cox01.lasso)
cox01.lasso.path <- penalized(Surv(SMARTc$TEVENT, SMARTc$EVENT) ~ cox01$x, model='cox', standardize=T, lambda1=1, steps=10)
plotpath(cox01.lasso.path, standardize=T)
title('Standardized coefficients')
#-------------------------------------------------------------------------------
## PAPER STEP 5: Model performance ##
#COX01DEF:
cox01def <-cph(Surv(TEVENT,EVENT) ~ ifelse(AGE>50, (AGE-50)^2,0) + BMI +
HDL+DIABETES +HISTCAR2 + log(CREAT)+ as.factor(ALBUMIN)+STENOSIS+IMT,data = SMARTc, x=T, y=T, surv=T)
dd <- datadist(SMARTc)
options(datadist="dd")
anova(cox01def)
summary(cox01def)
#C-statistic and KM-curves of cox01def
#C-statistic:
rcorr.cens(-cox01def$linear.pred,cox01def$y)
#KM-curves:
validate(cox01def, B=20, dxy=T)
cox01def$linear.pred
summary(cox01def$linear.predictor)
SMARTc$lp4 <- cut2(as.numeric(cox01def$linear.predictor), g=4)
dd <- datadist(SMARTc)
options(datadist="dd")
table(SMARTc$lp4)
levels(SMARTc$lp4) <- as.character(1:4)
KMrisk4 <- survfit(Surv(TEVENT,EVENT) ~ lp4, data=SMARTc)
KMrisk4 # Fig 23.8
par(las=1)
survplot(KMrisk4, conf="none", n.risk=T, xlab="Follow-up (days)", cex.n.risk=1, adj=.5,
ylab="Fraction free of cardiovascular event", time.inc=365)
#-------------------------------------------------------------------------------
## PAPER STEP 6: Model performance ##
#Bootstrap:
val.step <- validate(cox01, bw=T, rule="aic", type="individual")
vcox01def <- validate.cph(cox01def, method='boot', B=200, bw=T,type="individual",dxy=T)
coef.sh <- cox01def$coefficients*vcox01def["Slope","index.corrected"]
lp.s <-cox01def$x%*%coef.sh
rcorr.cens (-lp.s,Surv(SMARTc$TEVENT,SMARTc$EVENT))
# SHRINKAGE FACTORS
Lasso.coefs <- SMARTpath$b.predictor[SMARTpath$aic==min(SMARTpath$aic),]
Lasso.coefs <- Lasso.coefs[SMARTpath$b.predictor[SMARTpath$aic==min(SMARTpath$aic),] > 0 |
SMARTpath$b.predictor[SMARTpath$aic==min(SMARTpath$aic),] < 0]
Lasso.coefs
cox01def$coef
round(Lasso.coefs/cox01def$coef,3)
mean(Lasso.coefs/cox01def$coef) # 0.80; stronger shrinkage for weaker effects
# v Houwelingen shrinkage factor, with 17 df --> 0.94
(cox01def$stats[3] - cox01$stats[4]) / cox01def$stats[3]
#-------------------------------------------------------------------------------
## PAPER STEP 7: Model performance ##
# Nomogram with lasso coefs as in the original paper ( EERST cox01def draaien, daarvoor moet worden geimputeerd).
lassocoefs <- c(AGE=0.0012,BMI=-0.001,HDL=-0.16,DIABETES=0.11,HISTCAR2=0.33,CREAT=0.71,'ALBUMIN=2'=0.13,'ALBUMIN=3'=0.20,STENOSIS=0.16,IMT=0.50)
cox01def$coefficients <- lassocoefs
surv <- Survival(cox01def)
surv3 <- function(x) surv(3*365.25,lp=x)
surv5 <- function(x) surv(5*365.25,lp=x)
nomogram(cox01def, fun=list(surv3, surv5), lp=F, AGE=c(15,seq(50,85,5)), CREAT=c(60,80,100,150,200,400),
funlabel=c('3-year survival', '5-year survival'), maxscale=10, lp.at=c(0,1,2,2.4), fun.at=c(.93,.9,.85,.8,.7,.6,.5))
# nomogram ( if you run this before you run the above, you get the nomogram based on the coefficients of the cph (cox01def) fit)
surv <- Survival(cox01def)
surv1 <- function(x) surv(1*365.25,lp=x)
surv3 <- function(x) surv(3*365.25,lp=x)
surv5 <- function(x) surv(5*365.25,lp=x)
par(mfrow=c(1,1))
nomogram(cox01def, fun=list(surv3,surv5), lp=F, AGE=c(15,seq(50,85,5)), CREAT=c(60,80,100,150,200,400), funlabel=c('3-year survival', '5-year survival'), maxscale=10,lp.at=c(0,1,2,2.4), fun.at=c(.93,.9,.85,.8,.7,.6,.5))
#-------------------------------------------------------------------------------
## EXTRA
# how I computed the absolute risks per patient
#make baseline database
baseline <- data.frame(AGE=0,BMI=0,HDL=0,DIABETES=0,HISTCAR2=0,CREAT=1,ALBUMIN=1,STENOSIS=0,IMT=0)
summary(survfit(cox01def, newdata=baseline))
plot(survfit(cox01def, newdata=baseline))
# make new variables in which tranformed variables are transformed.
SMARTc_new<-SMARTc
#New columns in SMARTc, columns 35-38:
SMARTc_new$ALBUMIN2<-as.numeric(SMARTc$ALBUMIN==2)
SMARTc_new$ALBUMIN3<-as.numeric(SMARTc$ALBUMIN==3)
SMARTc_new$logCREAT<-log(SMARTc$CREAT)
SMARTc_new$newAGE<-ifelse(SMARTc$AGE>50, (SMARTc$AGE-50)^2,0)
##IMPORTANT: NUMBERS ARE THE NUMBERS OF THE COLOMNS IN ORDER OF THE COEFFICIENTS### THIS CAN CHANGE PER DATASET
SMARTc_new[1,]
cox01def$coefficients
# you can give the coefficients every value you like it to have : see above ...$coefficients<- .....
linpred<-as.matrix(SMARTc_NEW[,c(38,17,19,5,32,37,35,36,10,25)])%*%cox01def$coefficients # depends wheterh you have overwritten the coeffcients if lasso-beta's are used or not
#compute risks
base1yr<-min(survfit(cox01def, newdata=baseline)$surv[survfit(cox01def, newdata=baseline)$time<=365.25])
base3yr<-min(survfit(cox01def, newdata=baseline)$surv[survfit(cox01def, newdata=baseline)$time<=1095.75])
base5yr<-min(survfit(cox01def, newdata=baseline)$surv[survfit(cox01def, newdata=baseline)$time<=1826.25])
SMARTc$oneyrsurv <- base1yr^(exp(linpred))
SMARTc$threeyrsurv <- base3yr^(exp(linpred))
SMARTc$fiveyrsurv <- base5yr^(exp(linpred))
SMARTc$linpred <- linpred
## SAVE IMAGE ##
save(list = ls(all=TRUE), file = "#choose name#.RData")
