R Code “Development of nomograms to predict seizure outcome after antiepileptic drug withdrawal following epilepsy surgery”.
knitr::opts_chunk$set(warning=FALSE, size = 8)
#fig.width=8 Setup
Load necessary R libraries:
library(survival)
library(survminer)
library(rms)
library(readxl)
library(MASS)
library(ggplot2)
library(DynNom)Load data from spreadsheet:
data <- read_excel("Only_sz_free_all_cohort_(all variables).xlsx")## readxl works best with a newer version of the tibble package.
## You currently have tibble v1.4.2.
## Falling back to column name repair from tibble <= v1.4.2.
## Message displays once per session.ddist <- datadist(data)
options(datadist = 'ddist')Model 1 (Cox regression)
Cohort: Seizure free before reduction Outcome: Time to seizure relapse
km <- survfit(Surv(data$time_sz, data$sz) ~ begin, data = data, type =
"kaplan-meier", conf.type = "plain")
km## Call: survfit(formula = Surv(data$time_sz, data$sz) ~ begin, data = data,
## type = "kaplan-meier", conf.type = "plain")
##
## n events median 0.95LCL 0.95UCL
## 350.0 106.0 21.5 14.0 NAsummary(km, times = c(0,1,2,5,10, 15))## Call: survfit(formula = Surv(data$time_sz, data$sz) ~ begin, data = data,
## type = "kaplan-meier", conf.type = "plain")
##
## time n.risk n.event survival std.err lower 95% CI upper 95% CI
## 0 350 0 1.000 0.0000 1.000 1.000
## 1 303 29 0.915 0.0151 0.885 0.945
## 2 268 17 0.862 0.0190 0.824 0.899
## 5 179 38 0.730 0.0255 0.680 0.780
## 10 84 15 0.653 0.0298 0.595 0.712
## 15 25 6 0.575 0.0408 0.495 0.655ggsurvplot(km, xlim = c(0,15), data = data, xlab = "Time(years)", ylab = "Proportion seizure free patients", color = 'black')
model1 <- coxph(Surv(time = data$time_sz, event = data$sz) ~ data$auras + data$time_begin + data$duration + data$drugs + data$gtcs + data$nosz + data$temp + data$MRI + data$complete + data$pathology_HS + data$age_at_surgery, x = TRUE, y = TRUE)
summary(model1)## Call:
## coxph(formula = Surv(time = data$time_sz, event = data$sz) ~
## data$auras + data$time_begin + data$duration + data$drugs +
## data$gtcs + data$nosz + data$temp + data$MRI + data$complete +
## data$pathology_HS + data$age_at_surgery, x = TRUE, y = TRUE)
##
## n= 350, number of events= 106
##
## coef exp(coef) se(coef) z Pr(>|z|)
## data$auras 1.399e+00 4.052e+00 3.244e-01 4.314 1.6e-05 ***
## data$time_begin -3.757e-02 9.631e-01 2.733e-02 -1.375 0.169
## data$duration 1.314e-02 1.013e+00 1.117e-02 1.177 0.239
## data$drugs 1.860e-01 1.204e+00 1.239e-01 1.502 0.133
## data$gtcs 3.576e-01 1.430e+00 2.267e-01 1.577 0.115
## data$nosz -2.112e-04 9.998e-01 2.044e-03 -0.103 0.918
## data$temp -3.252e-01 7.224e-01 4.204e-01 -0.774 0.439
## data$MRI -5.579e-02 9.457e-01 6.175e-01 -0.090 0.928
## data$complete -1.391e+01 9.068e-07 1.895e+03 -0.007 0.994
## data$pathology_HS -1.049e-01 9.004e-01 2.779e-01 -0.377 0.706
## data$age_at_surgery -9.888e-03 9.902e-01 1.291e-02 -0.766 0.444
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## data$auras 4.052e+00 2.468e-01 2.1460 7.653
## data$time_begin 9.631e-01 1.038e+00 0.9129 1.016
## data$duration 1.013e+00 9.869e-01 0.9913 1.036
## data$drugs 1.204e+00 8.303e-01 0.9448 1.535
## data$gtcs 1.430e+00 6.994e-01 0.9169 2.230
## data$nosz 9.998e-01 1.000e+00 0.9958 1.004
## data$temp 7.224e-01 1.384e+00 0.3169 1.647
## data$MRI 9.457e-01 1.057e+00 0.2819 3.173
## data$complete 9.068e-07 1.103e+06 0.0000 Inf
## data$pathology_HS 9.004e-01 1.111e+00 0.5223 1.552
## data$age_at_surgery 9.902e-01 1.010e+00 0.9654 1.016
##
## Concordance= 0.639 (se = 0.03 )
## Rsquare= 0.062 (max possible= 0.961 )
## Likelihood ratio test= 22.29 on 11 df, p=0.02
## Wald test = 25.47 on 11 df, p=0.008
## Score (logrank) test = 27.76 on 11 df, p=0.004Tests for assumptions
Proportional hazards (Schoenfeld residuals)
ph<-cox.zph(model1)
print(ph)## rho chisq p
## data$auras -0.05978 3.93e-01 0.531
## data$time_begin 0.14118 2.07e+00 0.151
## data$duration -0.14262 1.86e+00 0.173
## data$drugs 0.02498 6.91e-02 0.793
## data$gtcs -0.00500 2.52e-03 0.960
## data$nosz -0.03425 1.08e-01 0.743
## data$temp -0.03673 1.48e-01 0.701
## data$MRI -0.05249 2.84e-01 0.594
## data$complete 0.00195 1.86e-11 1.000
## data$pathology_HS 0.03258 1.12e-01 0.738
## data$age_at_surgery 0.09603 1.11e+00 0.292
## GLOBAL NA 5.11e+00 0.926ggcoxzph(ph)
####Multicollinearity (Variance inflation factor)
vif(model1)## data$auras data$time_begin data$duration
## 1.350470 1.242102 1.844854
## data$drugs data$gtcs data$nosz
## 1.084586 1.030591 1.058922
## data$temp data$MRI data$complete
## 1.452994 1.108288 1.000000
## data$pathology_HS data$age_at_surgery
## 1.657564 1.803111Linearity (Martingale residuals)
res <- residuals(model1, type = 'martingale')
plot(data$time_begin, res)
lines(loess.smooth(data$time_begin, res), lwd=2)
plot(data$duration, res)
lines(loess.smooth(data$duration, res), lwd=2)
plot(data$age_at_surgery, res)
lines(loess.smooth(data$age_at_surgery, res), lwd=2)
plot(data$nosz, res)
lines(loess.smooth(data$nosz, res), lwd=2)
plot(data$duration, res)
lines(loess.smooth(data$duration, res), lwd=2)
Stepwise reduction using AIC
step(model1, data = data, direction = "both", criterion = "AIC")## Start: AIC=1139.16
## Surv(time = data$time_sz, event = data$sz) ~ data$auras + data$time_begin +
## data$duration + data$drugs + data$gtcs + data$nosz + data$temp +
## data$MRI + data$complete + data$pathology_HS + data$age_at_surgery
##
## Df AIC
## - data$MRI 1 1137.2
## - data$nosz 1 1137.2
## - data$pathology_HS 1 1137.3
## - data$age_at_surgery 1 1137.8
## - data$complete 1 1137.8
## - data$temp 1 1137.8
## - data$duration 1 1138.6
## <none> 1139.2
## - data$time_begin 1 1139.2
## - data$drugs 1 1139.4
## - data$gtcs 1 1139.8
## - data$auras 1 1152.0
##
## Step: AIC=1137.17
## Surv(time = data$time_sz, event = data$sz) ~ data$auras + data$time_begin +
## data$duration + data$drugs + data$gtcs + data$nosz + data$temp +
## data$complete + data$pathology_HS + data$age_at_surgery
##
## Df AIC
## - data$nosz 1 1135.2
## - data$pathology_HS 1 1135.3
## - data$age_at_surgery 1 1135.8
## - data$complete 1 1135.8
## - data$temp 1 1135.8
## - data$duration 1 1136.6
## <none> 1137.2
## - data$time_begin 1 1137.2
## - data$drugs 1 1137.4
## - data$gtcs 1 1137.8
## + data$MRI 1 1139.2
## - data$auras 1 1150.0
##
## Step: AIC=1135.18
## Surv(time = data$time_sz, event = data$sz) ~ data$auras + data$time_begin +
## data$duration + data$drugs + data$gtcs + data$temp + data$complete +
## data$pathology_HS + data$age_at_surgery
##
## Df AIC
## - data$pathology_HS 1 1133.3
## - data$age_at_surgery 1 1133.8
## - data$complete 1 1133.8
## - data$temp 1 1133.9
## - data$duration 1 1134.6
## <none> 1135.2
## - data$time_begin 1 1135.2
## - data$drugs 1 1135.4
## - data$gtcs 1 1135.9
## + data$nosz 1 1137.2
## + data$MRI 1 1137.2
## - data$auras 1 1148.0
##
## Step: AIC=1133.35
## Surv(time = data$time_sz, event = data$sz) ~ data$auras + data$time_begin +
## data$duration + data$drugs + data$gtcs + data$temp + data$complete +
## data$age_at_surgery
##
## Df AIC
## - data$age_at_surgery 1 1131.9
## - data$temp 1 1131.9
## - data$complete 1 1132.0
## - data$duration 1 1132.6
## <none> 1133.3
## - data$time_begin 1 1133.5
## - data$drugs 1 1133.5
## - data$gtcs 1 1134.0
## + data$pathology_HS 1 1135.2
## + data$MRI 1 1135.3
## + data$nosz 1 1135.3
## - data$auras 1 1146.0
##
## Step: AIC=1131.86
## Surv(time = data$time_sz, event = data$sz) ~ data$auras + data$time_begin +
## data$duration + data$drugs + data$gtcs + data$temp + data$complete
##
## Df AIC
## - data$temp 1 1130.3
## - data$complete 1 1130.5
## - data$duration 1 1130.6
## - data$time_begin 1 1131.8
## <none> 1131.9
## - data$drugs 1 1131.9
## - data$gtcs 1 1132.5
## + data$age_at_surgery 1 1133.3
## + data$pathology_HS 1 1133.8
## + data$MRI 1 1133.8
## + data$nosz 1 1133.9
## - data$auras 1 1144.0
##
## Step: AIC=1130.27
## Surv(time = data$time_sz, event = data$sz) ~ data$auras + data$time_begin +
## data$duration + data$drugs + data$gtcs + data$complete
##
## Df AIC
## - data$complete 1 1129.0
## - data$duration 1 1129.2
## - data$drugs 1 1130.1
## - data$time_begin 1 1130.1
## <none> 1130.3
## - data$gtcs 1 1131.0
## + data$temp 1 1131.9
## + data$age_at_surgery 1 1131.9
## + data$MRI 1 1132.2
## + data$nosz 1 1132.3
## + data$pathology_HS 1 1132.3
## - data$auras 1 1142.4
##
## Step: AIC=1129.01
## Surv(time = data$time_sz, event = data$sz) ~ data$auras + data$time_begin +
## data$duration + data$drugs + data$gtcs
##
## Df AIC
## - data$duration 1 1128.0
## - data$time_begin 1 1128.8
## - data$drugs 1 1128.9
## <none> 1129.0
## - data$gtcs 1 1129.6
## + data$complete 1 1130.3
## + data$temp 1 1130.5
## + data$age_at_surgery 1 1130.6
## + data$MRI 1 1131.0
## + data$pathology_HS 1 1131.0
## + data$nosz 1 1131.0
## - data$auras 1 1141.2
##
## Step: AIC=1127.98
## Surv(time = data$time_sz, event = data$sz) ~ data$auras + data$time_begin +
## data$drugs + data$gtcs
##
## Df AIC
## <none> 1128.0
## - data$drugs 1 1128.1
## - data$time_begin 1 1128.2
## - data$gtcs 1 1129.0
## + data$duration 1 1129.0
## + data$complete 1 1129.2
## + data$temp 1 1129.2
## + data$pathology_HS 1 1129.8
## + data$nosz 1 1130.0
## + data$MRI 1 1130.0
## + data$age_at_surgery 1 1130.0
## - data$auras 1 1141.2## Call:
## coxph(formula = Surv(time = data$time_sz, event = data$sz) ~
## data$auras + data$time_begin + data$drugs + data$gtcs, x = TRUE,
## y = TRUE)
##
## coef exp(coef) se(coef) z p
## data$auras 1.3748 3.9542 0.3091 4.45 8.7e-06
## data$time_begin -0.0376 0.9631 0.0264 -1.43 0.153
## data$drugs 0.1727 1.1885 0.1177 1.47 0.142
## data$gtcs 0.3761 1.4566 0.2242 1.68 0.093
##
## Likelihood ratio test=19.47 on 4 df, p=6e-04
## n= 350, number of events= 106Final model
model1f <- coxph(Surv(time = data$time_sz, event = data$sz) ~ data$auras + data$time_begin + data$drugs + data$gtcs, x = TRUE, y = TRUE)
summary(model1f)## Call:
## coxph(formula = Surv(time = data$time_sz, event = data$sz) ~
## data$auras + data$time_begin + data$drugs + data$gtcs, x = TRUE,
## y = TRUE)
##
## n= 350, number of events= 106
##
## coef exp(coef) se(coef) z Pr(>|z|)
## data$auras 1.37479 3.95424 0.30908 4.448 8.67e-06 ***
## data$time_begin -0.03762 0.96308 0.02636 -1.427 0.1535
## data$drugs 0.17270 1.18851 0.11771 1.467 0.1423
## data$gtcs 0.37611 1.45661 0.22419 1.678 0.0934 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## data$auras 3.9542 0.2529 2.1576 7.247
## data$time_begin 0.9631 1.0383 0.9146 1.014
## data$drugs 1.1885 0.8414 0.9436 1.497
## data$gtcs 1.4566 0.6865 0.9387 2.260
##
## Concordance= 0.626 (se = 0.03 )
## Rsquare= 0.054 (max possible= 0.961 )
## Likelihood ratio test= 19.47 on 4 df, p=6e-04
## Wald test = 23.82 on 4 df, p=9e-05
## Score (logrank) test = 25.36 on 4 df, p=4e-05AIC(model1f)## [1] 1127.981Internal validation
Discrimination:
fitmodel1f <- cph(Surv(time = time_sz, event = sz) ~ auras + time_begin + drugs + gtcs, data = data, x = TRUE, y = TRUE, surv = TRUE)
summary(fitmodel1f)## Effects Response : Surv(time = time_sz, event = sz)
##
## Factor Low High Diff. Effect S.E. Lower 0.95
## auras 0.00000 1.0000 1.0000 1.37480 0.309080 0.769020
## Hazard Ratio 0.00000 1.0000 1.0000 3.95430 NA 2.157700
## time_begin 0.87569 2.9951 2.1194 -0.07974 0.055863 -0.189230
## Hazard Ratio 0.87569 2.9951 2.1194 0.92336 NA 0.827600
## drugs 2.00000 3.0000 1.0000 0.17270 0.117710 -0.058017
## Hazard Ratio 2.00000 3.0000 1.0000 1.18850 NA 0.943630
## gtcs 0.00000 1.0000 1.0000 0.37611 0.224190 -0.063282
## Hazard Ratio 0.00000 1.0000 1.0000 1.45660 NA 0.938680
## Upper 0.95
## 1.98060
## 7.24700
## 0.02975
## 1.03020
## 0.40341
## 1.49690
## 0.81551
## 2.26030rms::validate(fitmodel1f, dxy = TRUE, B = 1000)## index.orig training test optimism index.corrected n
## Dxy 0.2527 0.2509 0.2253 0.0256 0.2271 1000
## R2 0.0563 0.0652 0.0481 0.0172 0.0391 1000
## Slope 1.0000 1.0000 0.9031 0.0969 0.9031 1000
## D 0.0162 0.0192 0.0137 0.0055 0.0107 1000
## U -0.0018 -0.0018 0.0017 -0.0035 0.0017 1000
## Q 0.0180 0.0209 0.0119 0.0090 0.0090 1000
## g 0.3985 0.4357 0.3750 0.0607 0.3379 1000Calibration for year 2:
fitmodel1f_2y <- cph(Surv(time = time_sz, event = sz) ~ auras + time_begin + drugs + gtcs , data = data, x = TRUE, y = TRUE, surv = TRUE, time.inc = 2)
calibrate1_2y <- calibrate(fitmodel1f_2y, u = 2, cmethod = 'KM', B = 1000, m = 80)## Using Cox survival estimates at 2 Dayssummary(calibrate1_2y)## index.orig training test
## Min. :-0.024432 Min. :-0.090672 Min. :-0.067770
## 1st Qu.:-0.017645 1st Qu.:-0.028307 1st Qu.:-0.037796
## Median :-0.001278 Median : 0.007582 Median :-0.016245
## Mean : 0.001495 Mean :-0.012935 Mean :-0.021076
## 3rd Qu.: 0.017862 3rd Qu.: 0.022954 3rd Qu.: 0.000476
## Max. : 0.032968 Max. : 0.023767 Max. : 0.015956
## mean.optimism mean.corrected n
## Min. :-0.0058728 Min. :-0.027432 Min. : 964.0
## 1st Qu.: 0.0007814 1st Qu.:-0.013990 1st Qu.: 978.2
## Median : 0.0046925 Median :-0.001534 Median : 986.5
## Mean : 0.0033117 Mean :-0.001817 Mean : 984.2
## 3rd Qu.: 0.0072229 3rd Qu.: 0.010639 3rd Qu.: 992.5
## Max. : 0.0097346 Max. : 0.023233 Max. :1000.0
## mean.predicted KM KM.corrected std.err
## Min. :0.7876 Min. :0.7722 Min. :0.7781 Min. :0.03216
## 1st Qu.:0.8431 1st Qu.:0.8210 1st Qu.:0.8202 1st Qu.:0.03300
## Median :0.8723 Median :0.8765 Median :0.8702 Median :0.04045
## Mean :0.8609 Mean :0.8624 Mean :0.8591 Mean :0.04319
## 3rd Qu.:0.8901 3rd Qu.:0.9180 3rd Qu.:0.9091 3rd Qu.:0.05064
## Max. :0.9116 Max. :0.9245 Max. :0.9181 Max. :0.05969plot(calibrate1_2y)
plot(calibrate1_2y, xlim = c(0,1), ylim = c(0,1))
Calibration for year 5:
fitmodel1f_5y <- cph(Surv(time = time_sz, event = sz) ~ auras + time_begin +drugs + gtcs , data = data, x = TRUE, y = TRUE, surv = TRUE, time.inc = 5)
calibrate1_5y <- calibrate(fitmodel1f_5y, u = 5, cmethod='KM', B = 1000, m = 80)## Using Cox survival estimates at 5 Daysplot(calibrate1_5y)
plot(calibrate1_5y, xlim = c(0,1), ylim = c(0,1), df = 5)
Nomogram
ddist <- datadist(data)
options(datadist = 'ddist')
fitmodel1f$Design$label <- c("Auras before withdrawal", "Time to begin withdrawal (years from surgery)", " Number of AEDs at time of surgery","GTCS before surgery")
surv.fitmodel1f <- Survival(fitmodel1f)
nom.cox1 <- nomogram(fitmodel1f, fun = list(function(x) surv.fitmodel1f(1, x), function(x) surv.fitmodel1f(5, x)), funlabel = c("2-year seizure freedom", "5-year seizure freedom"), lp = T, fun.at = c(0.99,0.95,0.9,0.8,0.7,0.6,0.5,0.4,0.3,0.2))
plot(nom.cox1)
data$auras <- as.factor(data$auras)
levels(data$auras) <- c("No", "Yes")
data$gtcs <- as.factor(data$gtcs)
levels(data$gtcs) <- c("No", "Yes")//{r} web <- coxph(Surv(time_sz, sz) ~ auras + time_begin + drugs + gtcs, data = data) DNbuilder.coxph(web, data = Only_sz_free_all_cohort, clevel = 0.95, m.summary = c("formatted"), covariate = c("numeric")) setwd("..") //
Model 2 (Cox regression)
Cohort: Seizure free before reduction Outcome: Time to withdrawal all
model2 <- coxph(Surv(time = begin_all_follow, event = all) ~ auras +age_at_surgery + drugs + time_begin + pathology_HS + gtcs + MRI + psychiatric_pre_any + duration, data = data, x = TRUE, y = TRUE)
print(model2)## Call:
## coxph(formula = Surv(time = begin_all_follow, event = all) ~
## auras + age_at_surgery + drugs + time_begin + pathology_HS +
## gtcs + MRI + psychiatric_pre_any + duration, data = data,
## x = TRUE, y = TRUE)
##
## coef exp(coef) se(coef) z p
## aurasYes -0.49809 0.60769 0.43974 -1.13 0.2573
## age_at_surgery -0.02830 0.97209 0.01128 -2.51 0.0121
## drugs -0.35980 0.69782 0.11731 -3.07 0.0022
## time_begin -0.06649 0.93567 0.03253 -2.04 0.0410
## pathology_HS -0.08387 0.91955 0.19476 -0.43 0.6667
## gtcsYes -0.04517 0.95583 0.18056 -0.25 0.8024
## MRI 0.69351 2.00072 0.59904 1.16 0.2470
## psychiatric_pre_any -0.31818 0.72747 0.17762 -1.79 0.0732
## duration 0.00933 1.00938 0.01004 0.93 0.3525
##
## Likelihood ratio test=31.19 on 9 df, p=3e-04
## n= 350, number of events= 141summary(model2)## Call:
## coxph(formula = Surv(time = begin_all_follow, event = all) ~
## auras + age_at_surgery + drugs + time_begin + pathology_HS +
## gtcs + MRI + psychiatric_pre_any + duration, data = data,
## x = TRUE, y = TRUE)
##
## n= 350, number of events= 141
##
## coef exp(coef) se(coef) z Pr(>|z|)
## aurasYes -0.498091 0.607690 0.439735 -1.133 0.25734
## age_at_surgery -0.028302 0.972095 0.011278 -2.510 0.01209 *
## drugs -0.359798 0.697817 0.117311 -3.067 0.00216 **
## time_begin -0.066489 0.935674 0.032535 -2.044 0.04099 *
## pathology_HS -0.083867 0.919553 0.194757 -0.431 0.66674
## gtcsYes -0.045175 0.955830 0.180556 -0.250 0.80243
## MRI 0.693506 2.000718 0.599037 1.158 0.24699
## psychiatric_pre_any -0.318181 0.727471 0.177617 -1.791 0.07323 .
## duration 0.009333 1.009376 0.010038 0.930 0.35253
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## aurasYes 0.6077 1.6456 0.2567 1.4387
## age_at_surgery 0.9721 1.0287 0.9508 0.9938
## drugs 0.6978 1.4330 0.5545 0.8782
## time_begin 0.9357 1.0687 0.8779 0.9973
## pathology_HS 0.9196 1.0875 0.6278 1.3470
## gtcsYes 0.9558 1.0462 0.6710 1.3617
## MRI 2.0007 0.4998 0.6184 6.4727
## psychiatric_pre_any 0.7275 1.3746 0.5136 1.0304
## duration 1.0094 0.9907 0.9897 1.0294
##
## Concordance= 0.648 (se = 0.026 )
## Rsquare= 0.085 (max possible= 0.987 )
## Likelihood ratio test= 31.19 on 9 df, p=3e-04
## Wald test = 28.35 on 9 df, p=8e-04
## Score (logrank) test = 29.01 on 9 df, p=6e-04Tests for assumptions
Proportional hazards (Schoenfeld residuals)
ph2<-cox.zph(model2)
print(ph2)## rho chisq p
## aurasYes 0.19506 6.5562 0.0105
## age_at_surgery -0.08675 1.0676 0.3015
## drugs 0.16194 4.1516 0.0416
## time_begin -0.02325 0.0727 0.7875
## pathology_HS -0.10230 1.4964 0.2212
## gtcsYes 0.06150 0.5633 0.4529
## MRI 0.00333 0.0016 0.9681
## psychiatric_pre_any 0.12537 2.2753 0.1314
## duration 0.10211 1.3995 0.2368
## GLOBAL NA 16.5797 0.0557ggcoxzph(ph2)
##### Log(-log(S(t))) test
check <- npsurv(Surv(begin_all_follow, all) ~ auras , data = data)
survplot(check, loglog = T)
Multicollinearity (Variance inflation factor)
vif(model2)## aurasYes age_at_surgery drugs
## 1.109297 1.676398 1.077035
## time_begin pathology_HS gtcsYes
## 1.109815 1.170926 1.030007
## MRI psychiatric_pre_any duration
## 1.053272 1.025763 1.807382Linearity (Martingale residuals)
res2 <- residuals(model2, type = 'martingale')
plot(data$time_begin, res2)
lines(loess.smooth(data$time_begin, res2), lwd=2)
plot(data$duration, res2)
lines(loess.smooth(data$duration, res2), lwd=2)
plot(data$age_at_surgery, res2)
lines(loess.smooth(data$age_at_surgery, res2), lwd=2)
Stepwise reduction using AIC
stepAIC(model2, direction = "both")## Start: AIC=1505.7
## Surv(time = begin_all_follow, event = all) ~ auras + age_at_surgery +
## drugs + time_begin + pathology_HS + gtcs + MRI + psychiatric_pre_any +
## duration
##
## Df AIC
## - gtcs 1 1503.8
## - pathology_HS 1 1503.9
## - duration 1 1504.6
## - auras 1 1505.2
## - MRI 1 1505.3
## <none> 1505.7
## - psychiatric_pre_any 1 1507.0
## - time_begin 1 1508.5
## - age_at_surgery 1 1510.5
## - drugs 1 1513.6
##
## Step: AIC=1503.76
## Surv(time = begin_all_follow, event = all) ~ auras + age_at_surgery +
## drugs + time_begin + pathology_HS + MRI + psychiatric_pre_any +
## duration
##
## Df AIC
## - pathology_HS 1 1502.0
## - duration 1 1502.6
## - auras 1 1503.2
## - MRI 1 1503.4
## <none> 1503.8
## - psychiatric_pre_any 1 1505.0
## + gtcs 1 1505.7
## - time_begin 1 1506.6
## - age_at_surgery 1 1508.5
## - drugs 1 1511.8
##
## Step: AIC=1501.96
## Surv(time = begin_all_follow, event = all) ~ auras + age_at_surgery +
## drugs + time_begin + MRI + psychiatric_pre_any + duration
##
## Df AIC
## - duration 1 1500.6
## - auras 1 1501.3
## - MRI 1 1501.5
## <none> 1502.0
## - psychiatric_pre_any 1 1503.4
## + pathology_HS 1 1503.8
## + gtcs 1 1503.9
## - time_begin 1 1504.8
## - age_at_surgery 1 1506.5
## - drugs 1 1509.9
##
## Step: AIC=1500.64
## Surv(time = begin_all_follow, event = all) ~ auras + age_at_surgery +
## drugs + time_begin + MRI + psychiatric_pre_any
##
## Df AIC
## - auras 1 1500.0
## - MRI 1 1500.4
## <none> 1500.6
## - psychiatric_pre_any 1 1501.9
## + duration 1 1502.0
## + gtcs 1 1502.6
## + pathology_HS 1 1502.6
## - time_begin 1 1503.2
## - age_at_surgery 1 1505.1
## - drugs 1 1508.0
##
## Step: AIC=1499.98
## Surv(time = begin_all_follow, event = all) ~ age_at_surgery +
## drugs + time_begin + MRI + psychiatric_pre_any
##
## Df AIC
## - MRI 1 1499.8
## <none> 1500.0
## + auras 1 1500.6
## + duration 1 1501.3
## - psychiatric_pre_any 1 1501.7
## + pathology_HS 1 1501.9
## + gtcs 1 1502.0
## - time_begin 1 1505.2
## - age_at_surgery 1 1505.4
## - drugs 1 1506.7
##
## Step: AIC=1499.76
## Surv(time = begin_all_follow, event = all) ~ age_at_surgery +
## drugs + time_begin + psychiatric_pre_any
##
## Df AIC
## <none> 1499.8
## + MRI 1 1500.0
## + auras 1 1500.4
## + duration 1 1500.9
## - psychiatric_pre_any 1 1501.2
## + pathology_HS 1 1501.7
## + gtcs 1 1501.7
## - time_begin 1 1505.0
## - age_at_surgery 1 1505.3
## - drugs 1 1506.7## Call:
## coxph(formula = Surv(time = begin_all_follow, event = all) ~
## age_at_surgery + drugs + time_begin + psychiatric_pre_any,
## data = data, x = TRUE, y = TRUE)
##
## coef exp(coef) se(coef) z p
## age_at_surgery -0.02394 0.97635 0.00891 -2.69 0.0072
## drugs -0.33250 0.71713 0.11335 -2.93 0.0034
## time_begin -0.07435 0.92834 0.03059 -2.43 0.0151
## psychiatric_pre_any -0.32048 0.72580 0.17585 -1.82 0.0684
##
## Likelihood ratio test=27.13 on 4 df, p=2e-05
## n= 350, number of events= 141Final model
model2f <- coxph(Surv(time = begin_all_follow, event = all) ~ age_at_surgery + drugs + time_begin + psychiatric_pre_any, data = data)
summary(model2f)## Call:
## coxph(formula = Surv(time = begin_all_follow, event = all) ~
## age_at_surgery + drugs + time_begin + psychiatric_pre_any,
## data = data)
##
## n= 350, number of events= 141
##
## coef exp(coef) se(coef) z Pr(>|z|)
## age_at_surgery -0.023937 0.976348 0.008912 -2.686 0.00723 **
## drugs -0.332498 0.717130 0.113345 -2.933 0.00335 **
## time_begin -0.074352 0.928345 0.030590 -2.431 0.01507 *
## psychiatric_pre_any -0.320483 0.725799 0.175845 -1.823 0.06837 .
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## exp(coef) exp(-coef) lower .95 upper .95
## age_at_surgery 0.9763 1.024 0.9594 0.9936
## drugs 0.7171 1.394 0.5743 0.8955
## time_begin 0.9283 1.077 0.8743 0.9857
## psychiatric_pre_any 0.7258 1.378 0.5142 1.0245
##
## Concordance= 0.639 (se = 0.026 )
## Rsquare= 0.075 (max possible= 0.987 )
## Likelihood ratio test= 27.13 on 4 df, p=2e-05
## Wald test = 25.2 on 4 df, p=5e-05
## Score (logrank) test = 25.64 on 4 df, p=4e-05Internal validation
Discrimination:
fitmodel2f <- cph( Surv(time = begin_all_follow, event = all) ~ age_at_surgery + strat(drugs) + time_begin + psychiatric_pre_any, data = data, surv=TRUE, x = TRUE, y = TRUE)
rms::validate(fitmodel2f, dxy = TRUE, B = 1000, u = 5)##
## Divergence or singularity in 384 samples## index.orig training test optimism index.corrected n
## Dxy 0.2800 0.2983 0.2606 0.0377 0.2423 616
## R2 0.0425 0.0510 0.0367 0.0143 0.0282 616
## Slope 1.0000 1.0000 0.8957 0.1043 0.8957 616
## D 0.0114 0.0140 0.0097 0.0043 0.0072 616
## U -0.0017 -0.0017 0.0011 -0.0028 0.0011 616
## Q 0.0131 0.0157 0.0087 0.0070 0.0061 616
## g 0.4162 0.4522 0.3816 0.0706 0.3457 616Calibration plot for year 2:
fitmodel2f <- cph( Surv(time = begin_all_follow, event = all) ~ age_at_surgery + strat(drugs) + time_begin + psychiatric_pre_any, data = data, surv=TRUE, x = TRUE, y = TRUE, time.inc = 2)
calibrate2_2y <- calibrate(fitmodel2f, u = 2, cmethod = 'KM', B = 1000, m = 80)## Using Cox survival estimates at 2 Days
##
## Divergence or singularity in 2 samples
##
## Divergence or singularity in 2 samples
##
## Divergence or singularity in 4 samples
##
## Divergence or singularity in 4 samples
##
## Divergence or singularity in 1 samples
##
## Divergence or singularity in 3 samples
##
## Divergence or singularity in 3 samples
##
## Divergence or singularity in 6 samples
##
## Divergence or singularity in 5 samples
##
## Divergence or singularity in 5 samples
##
## Divergence or singularity in 3 samples
##
## Divergence or singularity in 3 samples
##
## Divergence or singularity in 6 samples
##
## Divergence or singularity in 3 samples
##
## Divergence or singularity in 5 samples
##
## Divergence or singularity in 5 samples
##
## Divergence or singularity in 3 samples
##
## Divergence or singularity in 5 samples
##
## Divergence or singularity in 5 samples
##
## Divergence or singularity in 6 samples
##
## Divergence or singularity in 3 samples
##
## Divergence or singularity in 4 samples
##
## Divergence or singularity in 6 samples
##
## Divergence or singularity in 3 samples
##
## Divergence or singularity in 7 samples
##
## Divergence or singularity in 4 samples
##
## Divergence or singularity in 4 samples
##
## Divergence or singularity in 4 samples
##
## Divergence or singularity in 5 samples
##
## Divergence or singularity in 3 samples
##
## Divergence or singularity in 5 samples
##
## Divergence or singularity in 1 samples
##
## Divergence or singularity in 6 samples
##
## Divergence or singularity in 3 samples
##
## Divergence or singularity in 4 samples
##
## Divergence or singularity in 2 samples
##
## Divergence or singularity in 5 samples
##
## Divergence or singularity in 7 samples
##
## Divergence or singularity in 7 samples
##
## Divergence or singularity in 4 samples
##
## Divergence or singularity in 3 samples
##
## Divergence or singularity in 5 samples
##
## Divergence or singularity in 5 samples
##
## Divergence or singularity in 3 samples
##
## Divergence or singularity in 7 samples
##
## Divergence or singularity in 4 samples
##
## Divergence or singularity in 5 samples
##
## Divergence or singularity in 3 samples
##
## Divergence or singularity in 1 samples
##
## Divergence or singularity in 2 samples
##
## Divergence or singularity in 4 samples
##
## Divergence or singularity in 2 samples
##
## Divergence or singularity in 1 samples
##
## Divergence or singularity in 4 samples
##
## Divergence or singularity in 2 samples
##
## Divergence or singularity in 5 samples
##
## Divergence or singularity in 3 samples
##
## Divergence or singularity in 5 samples
##
## Divergence or singularity in 6 samples
##
## Divergence or singularity in 4 samples
##
## Divergence or singularity in 2 samples
##
## Divergence or singularity in 6 samples
##
## Divergence or singularity in 4 samples
##
## Divergence or singularity in 3 samples
##
## Divergence or singularity in 2 samples
##
## Divergence or singularity in 4 samples
##
## Divergence or singularity in 6 samples
##
## Divergence or singularity in 2 samples
##
## Divergence or singularity in 5 samples
##
## Divergence or singularity in 4 samples
##
## Divergence or singularity in 5 samples
##
## Divergence or singularity in 5 samples
##
## Divergence or singularity in 6 samples
##
## Divergence or singularity in 2 samples
##
## Divergence or singularity in 4 samples
##
## Divergence or singularity in 3 samples
##
## Divergence or singularity in 6 samples
##
## Divergence or singularity in 3 samples
##
## Divergence or singularity in 4 samples
##
## Divergence or singularity in 6 samples
##
## Divergence or singularity in 3 samples
##
## Divergence or singularity in 5 samples
##
## Divergence or singularity in 6 samples
##
## Divergence or singularity in 5 samples
##
## Divergence or singularity in 6 samples
##
## Divergence or singularity in 4 samples
##
## Divergence or singularity in 2 samples
##
## Divergence or singularity in 2 samples
##
## Divergence or singularity in 5 samples
##
## Divergence or singularity in 1 samples
##
## Divergence or singularity in 2 samples
##
## Divergence or singularity in 6 samples
##
## Divergence or singularity in 8 samples
##
## Divergence or singularity in 6 samples
##
## Divergence or singularity in 4 samples
##
## Divergence or singularity in 4 samples
##
## Divergence or singularity in 1 samples
##
## Divergence or singularity in 2 samples
##
## Divergence or singularity in 5 samplessummary(calibrate2_2y)## index.orig training test
## Min. :-0.047180 Min. :-0.036416 Min. :-0.05235
## 1st Qu.:-0.006462 1st Qu.:-0.018844 1st Qu.:-0.04327
## Median : 0.008006 Median : 0.006278 Median :-0.03590
## Mean :-0.002834 Mean : 0.001535 Mean :-0.01763
## 3rd Qu.: 0.011634 3rd Qu.: 0.026657 3rd Qu.:-0.01027
## Max. : 0.019833 Max. : 0.030000 Max. : 0.05362
## mean.optimism mean.corrected n mean.predicted
## Min. :-0.025102 Min. :-0.042383 Min. :601 Min. :0.5149
## 1st Qu.:-0.009873 1st Qu.:-0.022725 1st Qu.:601 1st Qu.:0.6163
## Median : 0.004303 Median :-0.011232 Median :601 Median :0.7034
## Mean : 0.002144 Mean :-0.004978 Mean :601 Mean :0.6942
## 3rd Qu.: 0.016320 3rd Qu.: 0.006515 3rd Qu.:601 3rd Qu.:0.7813
## Max. : 0.025073 Max. : 0.044934 Max. :601 Max. :0.8551
## KM KM.corrected std.err
## Min. :0.5347 Min. :0.5598 Min. :0.04643
## 1st Qu.:0.5859 1st Qu.:0.5957 1st Qu.:0.05834
## Median :0.6834 Median :0.6791 Median :0.07724
## Mean :0.6914 Mean :0.6892 Mean :0.07575
## 3rd Qu.:0.7889 3rd Qu.:0.7726 3rd Qu.:0.09464
## Max. :0.8640 Max. :0.8389 Max. :0.10209plot(calibrate2_2y)
plot(calibrate2_2y, xlim = c(0,1), ylim = c(0,1))
Calibration plot for year 5:
fitmodel2f <- cph( Surv(time = begin_all_follow, event = all) ~ age_at_surgery + strat(drugs) + time_begin + psychiatric_pre_any, data = data, surv=TRUE, x = TRUE, y = TRUE, time.inc = 5)
calibrate2_5y <- calibrate(fitmodel2f, u = 5, cmethod = 'KM', B = 1000, m = 80)## Using Cox survival estimates at 5 Days
##
## Divergence or singularity in 2 samples
##
## Divergence or singularity in 4 samples
##
## Divergence or singularity in 3 samples
##
## Divergence or singularity in 3 samples
##
## Divergence or singularity in 5 samples
##
## Divergence or singularity in 6 samples
##
## Divergence or singularity in 6 samples
##
## Divergence or singularity in 6 samples
##
## Divergence or singularity in 5 samples
##
## Divergence or singularity in 3 samples
##
## Divergence or singularity in 4 samples
##
## Divergence or singularity in 3 samples
##
## Divergence or singularity in 5 samples
##
## Divergence or singularity in 4 samples
##
## Divergence or singularity in 3 samples
##
## Divergence or singularity in 2 samples
##
## Divergence or singularity in 4 samples
##
## Divergence or singularity in 5 samples
##
## Divergence or singularity in 5 samples
##
## Divergence or singularity in 5 samples
##
## Divergence or singularity in 3 samples
##
## Divergence or singularity in 2 samples
##
## Divergence or singularity in 2 samples
##
## Divergence or singularity in 3 samples
##
## Divergence or singularity in 3 samples
##
## Divergence or singularity in 4 samples
##
## Divergence or singularity in 2 samples
##
## Divergence or singularity in 1 samples
##
## Divergence or singularity in 7 samples
##
## Divergence or singularity in 6 samples
##
## Divergence or singularity in 4 samples
##
## Divergence or singularity in 5 samples
##
## Divergence or singularity in 3 samples
##
## Divergence or singularity in 4 samples
##
## Divergence or singularity in 6 samples
##
## Divergence or singularity in 2 samples
##
## Divergence or singularity in 5 samples
##
## Divergence or singularity in 2 samples
##
## Divergence or singularity in 5 samples
##
## Divergence or singularity in 3 samples
##
## Divergence or singularity in 4 samples
##
## Divergence or singularity in 6 samples
##
## Divergence or singularity in 1 samples
##
## Divergence or singularity in 6 samples
##
## Divergence or singularity in 4 samples
##
## Divergence or singularity in 4 samples
##
## Divergence or singularity in 3 samples
##
## Divergence or singularity in 3 samples
##
## Divergence or singularity in 6 samples
##
## Divergence or singularity in 5 samples
##
## Divergence or singularity in 3 samples
##
## Divergence or singularity in 5 samples
##
## Divergence or singularity in 5 samples
##
## Divergence or singularity in 5 samples
##
## Divergence or singularity in 2 samples
##
## Divergence or singularity in 3 samples
##
## Divergence or singularity in 7 samples
##
## Divergence or singularity in 4 samples
##
## Divergence or singularity in 4 samples
##
## Divergence or singularity in 7 samples
##
## Divergence or singularity in 2 samples
##
## Divergence or singularity in 4 samples
##
## Divergence or singularity in 2 samples
##
## Divergence or singularity in 5 samples
##
## Divergence or singularity in 4 samples
##
## Divergence or singularity in 5 samples
##
## Divergence or singularity in 2 samples
##
## Divergence or singularity in 4 samples
##
## Divergence or singularity in 2 samples
##
## Divergence or singularity in 3 samples
##
## Divergence or singularity in 1 samples
##
## Divergence or singularity in 3 samples
##
## Divergence or singularity in 5 samples
##
## Divergence or singularity in 3 samples
##
## Divergence or singularity in 3 samples
##
## Divergence or singularity in 4 samples
##
## Divergence or singularity in 5 samples
##
## Divergence or singularity in 3 samples
##
## Divergence or singularity in 4 samples
##
## Divergence or singularity in 2 samples
##
## Divergence or singularity in 7 samples
##
## Divergence or singularity in 4 samples
##
## Divergence or singularity in 1 samples
##
## Divergence or singularity in 4 samples
##
## Divergence or singularity in 4 samples
##
## Divergence or singularity in 4 samples
##
## Divergence or singularity in 4 samples
##
## Divergence or singularity in 4 samples
##
## Divergence or singularity in 4 samples
##
## Divergence or singularity in 5 samples
##
## Divergence or singularity in 5 samples
##
## Divergence or singularity in 2 samples
##
## Divergence or singularity in 5 samples
##
## Divergence or singularity in 4 samples
##
## Divergence or singularity in 4 samples
##
## Divergence or singularity in 1 samples
##
## Divergence or singularity in 4 samples
##
## Divergence or singularity in 4 samplessummary(calibrate2_5y)## index.orig training test
## Min. :-0.0305527 Min. :-0.039522 Min. :-0.0020676
## 1st Qu.:-0.0149086 1st Qu.:-0.018107 1st Qu.:-0.0003758
## Median :-0.0091930 Median :-0.002436 Median : 0.0017203
## Mean :-0.0054368 Mean :-0.004124 Mean : 0.0118068
## 3rd Qu.: 0.0002788 3rd Qu.: 0.011547 3rd Qu.: 0.0139030
## Max. : 0.0271914 Max. : 0.027898 Max. : 0.0458543
## mean.optimism mean.corrected n mean.predicted
## Min. :-0.0229590 Min. :-0.034499 Min. :621.0 Min. :0.3606
## 1st Qu.:-0.0090543 1st Qu.:-0.028225 1st Qu.:621.0 1st Qu.:0.4727
## Median :-0.0006058 Median :-0.019017 Median :621.5 Median :0.5796
## Mean : 0.0001587 Mean :-0.005595 Mean :621.5 Mean :0.5737
## 3rd Qu.: 0.0086072 3rd Qu.: 0.003613 3rd Qu.:622.0 3rd Qu.:0.6806
## Max. : 0.0248053 Max. : 0.050150 Max. :622.0 Max. :0.7752
## KM KM.corrected std.err
## Min. :0.3878 Min. :0.4107 Min. :0.06812
## 1st Qu.:0.4565 1st Qu.:0.4656 1st Qu.:0.08289
## Median :0.5599 Median :0.5605 Median :0.10501
## Mean :0.5683 Mean :0.5681 Mean :0.10566
## 3rd Qu.:0.6717 3rd Qu.:0.6630 3rd Qu.:0.12778
## Max. :0.7655 Max. :0.7407 Max. :0.14452plot(calibrate2_2y)
plot(calibrate2_5y, xlim = c(0,1), ylim = c(0,1))
Nomogram
ddist <- datadist(data)
options(datadist = 'ddist')
fitmodel2f$Design$label <- c("Age at surgery", " Number of AEDs at time of surgery", "Time to begin withdrawal (years from surgery)", "Psychiatric comorbidity")
surv.fitmodel2f <- Survival(fitmodel2f)
nom.fitmodel2f<- nomogram(fitmodel2f, fun = list(function(x) surv.fitmodel2f(2,x), function(x) surv.fitmodel2f(5,x)), funlabel = c("2-years AED freedom (= 1-value)", "5-years AED freedom (= 1-value)"), lp = T, fun.at = c(0.99,0.95,0.9,0.8,0.7,0.6,0.5,0.4,0.3,0.2))
plot(nom.fitmodel2f)
//{r} web2 <- coxph(Surv(begin_all_follow, all) ~ age_at_surgery + drugs + time_begin + psychiatric_pre_any, data = data) getwd() DNbuilder.coxph(web2, data = Only_sz_free_all_cohort, clevel = 0.95, m.summary = c("formatted"), covariate = c("numeric"), ptype = c("st", "1-st")) setwd("..") //