Last updated: 2018-04-16
Code version: 899ad0d
library(dplyr)
library(data.table)
library(magrittr)
library(purrr)
library(here) # for tracking working directory
library(ggplot2)
library(epistats)
library(broom)
This excercise was intented for SPSS
In this practical exercise we will study the prognosis of patients with traumatic brain injury. We will assess the individual prognostic strength in univariable and multivariable analyses. The aim of the multivariable analysis is to adjust for correlation between prognostic factors, either to adjust for confounding (1) or to predict the prognosis with multiple predictors (2). The data are in SPSS format: “TBI.sav”. See below for a description of the dataset.
Data set traumatic brain injury (TBI.sav, n=2159) Patients are from the International and US Tirilazad trials; distributed here for didactic purposes only. The primary outcome was 6 months Glasgow Outcome Scale (range 1 to 5).
Name Description (coding: no/yes is coded as 0/1) Development (n=2159) Trial Study identification: 74 = Tirilazad international (n=1118) 75 = US (n=1041)
52% 48% d.gos Glasgow Outcome Scale at 6 months: 1 = dead 2 = vegetative 3 = severe disability 4 = moderate disability 5 = good recovery
23% 4% 12% 16% 44% d.mort Mortality at 6 months (0/1) 23% d.unfav Unfavorable outcome at 6 months (0/1) 39% cause Cause of injury 3 = Road traffic accident 4 = Motorbike 5 = Assault 6 = Domestic/fall 9 = Other
39% 20% 6% 17% 18% age Age, in years (median [interquartile range]) 29 [21 - 42] d.motor Admission motor score, range: 1 - 6 (median) 4 d.pupil Pupillary reactivity 1=both reactive 2=one reactive 3= no reactive pupils
70% 14% 16% pupil.i Single imputed pupillary reactivity, 1;2;3 70%/14%/16% hypoxia Hypoxia before / at admission, 1=yes 22% hypotens Hypotension before / at admission, 1=yes 19% ctclass Marshall CT classification, 1 - 6 (median) 2 tsah tSAH at CT, 1=yes 46% edh EDH at CT, 1=yes 13% cisterns Compressed cisterns at CT, 0=no;1=slightly;2=fully 57%/26%/10% shift Midline shift > 5 mm at CT, 1=yes 18% glucose Glucose at admission, mmol/l (median [interquartile range]) 8.2 [6.7 - 10.4] glucoset ph Truncated glucose values (median [interquartile range]) pH (median [interquartile range]) 8.2 [6.7 - 10.4] 7.4 [7.3 - 7.5] sodium Sodium, mmol/l (median [interquartile range]) 140 [137 - 142] sodiumt Truncated sodium (median [interquartile range]) 140 [137 - 142] hb Hb, g/dl (median [interquartile range]) 12.8 [10.9 - 14.3] hbt Truncated hb (median [interquartile range]) 12.8 [10.9 - 14.3] * d. variables denoted ‘derived’.
Exercises
Load in data
require(haven)
tbi <- read_spss(here("data", "TBI.sav"))
str(tbi)
Classes 'tbl_df', 'tbl' and 'data.frame': 2159 obs. of 24 variables:
$ trial :Class 'labelled' atomic [1:2159] 74 74 74 74 ...
.. ..- attr(*, "label")= chr "Study identification"
.. ..- attr(*, "format.spss")= chr "A6"
.. ..- attr(*, "display_width")= int 4
.. ..- attr(*, "labels")= Named chr [1:2] "74" "75"
.. .. ..- attr(*, "names")= chr [1:2] "Tirilazad International" "Tirilazad US"
$ d.gos :Class 'labelled' atomic [1:2159] 5 5 5 4 5 5 5 5 5 5 ...
.. ..- attr(*, "label")= chr "GOS at 6 months"
.. ..- attr(*, "format.spss")= chr "F2.0"
.. ..- attr(*, "labels")= Named num [1:5] 1 2 3 4 5
.. .. ..- attr(*, "names")= chr [1:5] "dead" "vegetative" "severe disability" "moderate disability" ...
$ d.mort : atomic 0 0 0 0 0 0 0 0 0 0 ...
..- attr(*, "label")= chr "Mortality at 6 months"
..- attr(*, "format.spss")= chr "F2.0"
$ d.unfav : atomic 0 0 0 0 0 0 0 0 0 0 ...
..- attr(*, "label")= chr "Unfavorable outcome at 6 months"
..- attr(*, "format.spss")= chr "F2.0"
$ cause :Class 'labelled' atomic [1:2159] 4 4 6 4 4 4 3 4 4 6 ...
.. ..- attr(*, "label")= chr "Cause of injury recoded"
.. ..- attr(*, "format.spss")= chr "F2.0"
.. ..- attr(*, "display_width")= int 10
.. ..- attr(*, "labels")= Named num [1:5] 3 4 5 6 9
.. .. ..- attr(*, "names")= chr [1:5] "Road traffic accident" "Motorbike" "Assault" "domestic/fall" ...
$ age : atomic 14 14 14 14 14 14 14 14 14 14 ...
..- attr(*, "label")= chr "Age in years"
..- attr(*, "format.spss")= chr "F2.0"
$ d.motor : atomic 5 4 4 4 5 3 5 5 4 5 ...
..- attr(*, "label")= chr "Admission motor score"
..- attr(*, "format.spss")= chr "F2.0"
$ d.pupil :Class 'labelled' atomic [1:2159] 1 1 1 1 1 1 1 1 1 2 ...
.. ..- attr(*, "label")= chr "Pupillary reactivity"
.. ..- attr(*, "format.spss")= chr "F2.0"
.. ..- attr(*, "labels")= Named num [1:3] 1 2 3
.. .. ..- attr(*, "names")= chr [1:3] "both reactive" "one reactive" "no reactive pupils"
$ pupil.i :Class 'labelled' atomic [1:2159] 1 1 1 1 1 1 1 1 1 2 ...
.. ..- attr(*, "label")= chr "Single imputed pupillary reactivity"
.. ..- attr(*, "format.spss")= chr "F2.0"
.. ..- attr(*, "labels")= Named num [1:3] 1 2 3
.. .. ..- attr(*, "names")= chr [1:3] "both reactive" "one reactive" "no reactive pupils"
$ hypoxia : atomic 0 0 1 0 0 0 0 0 0 0 ...
..- attr(*, "label")= chr "Hypoxia before / at admission"
..- attr(*, "format.spss")= chr "F2.0"
$ hypotens: atomic 0 0 0 0 0 0 0 0 0 0 ...
..- attr(*, "label")= chr "Hypotension before / at admission"
..- attr(*, "format.spss")= chr "F2.0"
$ ctclass : atomic 2 2 4 2 2 2 2 2 2 1 ...
..- attr(*, "label")= chr "CT classification according to Marshall"
..- attr(*, "format.spss")= chr "F2.0"
$ tsah : atomic 0 0 1 0 0 NA 1 1 NA 0 ...
..- attr(*, "label")= chr "tSAH at CT"
..- attr(*, "format.spss")= chr "F2.0"
$ edh : atomic 0 0 0 0 0 0 0 0 0 0 ...
..- attr(*, "label")= chr "EDH at CT"
..- attr(*, "format.spss")= chr "F2.0"
$ cisterns: atomic 1 1 NA 1 1 1 1 2 1 1 ...
..- attr(*, "label")= chr "Compressed cisterns at CT"
..- attr(*, "format.spss")= chr "F2.0"
$ shift : atomic 0 0 NA 1 0 0 0 1 0 0 ...
..- attr(*, "label")= chr "Midline shift > 5 mm at CT"
..- attr(*, "format.spss")= chr "F2.0"
$ d.sysbpt: atomic 119 131 136 137 138 ...
..- attr(*, "label")= chr "Systolic blood pressure (truncated, mm Hg)"
..- attr(*, "format.spss")= chr "F4.0"
$ glucose : atomic 7.7 7.4 7.56 5.7 6 ...
..- attr(*, "label")= chr "Glucose at admission (mmol/l)"
..- attr(*, "format.spss")= chr "F4.2"
$ glucoset: atomic 7.7 7.4 7.56 5.7 6 ...
..- attr(*, "label")= chr "Truncated glucose values"
..- attr(*, "format.spss")= chr "F4.2"
$ ph : atomic 7.35 7.33 7.49 NA NA ...
..- attr(*, "label")= chr "pH"
..- attr(*, "format.spss")= chr "F4.2"
$ sodium : atomic 143 143 141 144 142 138 141 135 137 137 ...
..- attr(*, "label")= chr "Sodium (mmol/l)"
..- attr(*, "format.spss")= chr "F4.1"
$ sodiumt : atomic 143 143 141 144 142 138 141 135 137 137 ...
..- attr(*, "label")= chr "Truncated sodium"
..- attr(*, "format.spss")= chr "F4.1"
$ hb : atomic 15 15 12.8 7.1 8.5 12.3 12.8 9.2 12.8 11 ...
..- attr(*, "label")= chr "hb (g/dl)"
..- attr(*, "format.spss")= chr "F4.1"
$ hbt : atomic 15 15 12.8 7.1 8.5 12.3 12.8 9.2 12.8 11 ...
..- attr(*, "label")= chr "Truncated hb"
..- attr(*, "format.spss")= chr "F4.1"
Coerce labelled
variables into factors, as R works with factors and labelled
variables are foreign to R.
Use the package haven
with the function as_factor
to get this done while preserving factor labels.
tbi %<>%
mutate_if(is.labelled, as_factor)
str(tbi)
Classes 'tbl_df', 'tbl' and 'data.frame': 2159 obs. of 24 variables:
$ trial : Factor w/ 2 levels "Tirilazad International",..: 1 1 1 1 1 1 1 1 1 1 ...
..- attr(*, "label")= chr "Study identification"
$ d.gos : Factor w/ 5 levels "dead","vegetative",..: 5 5 5 4 5 5 5 5 5 5 ...
..- attr(*, "label")= chr "GOS at 6 months"
$ d.mort : atomic 0 0 0 0 0 0 0 0 0 0 ...
..- attr(*, "label")= chr "Mortality at 6 months"
..- attr(*, "format.spss")= chr "F2.0"
$ d.unfav : atomic 0 0 0 0 0 0 0 0 0 0 ...
..- attr(*, "label")= chr "Unfavorable outcome at 6 months"
..- attr(*, "format.spss")= chr "F2.0"
$ cause : Factor w/ 5 levels "Road traffic accident",..: 2 2 4 2 2 2 1 2 2 4 ...
..- attr(*, "label")= chr "Cause of injury recoded"
$ age : atomic 14 14 14 14 14 14 14 14 14 14 ...
..- attr(*, "label")= chr "Age in years"
..- attr(*, "format.spss")= chr "F2.0"
$ d.motor : atomic 5 4 4 4 5 3 5 5 4 5 ...
..- attr(*, "label")= chr "Admission motor score"
..- attr(*, "format.spss")= chr "F2.0"
$ d.pupil : Factor w/ 3 levels "both reactive",..: 1 1 1 1 1 1 1 1 1 2 ...
..- attr(*, "label")= chr "Pupillary reactivity"
$ pupil.i : Factor w/ 3 levels "both reactive",..: 1 1 1 1 1 1 1 1 1 2 ...
..- attr(*, "label")= chr "Single imputed pupillary reactivity"
$ hypoxia : atomic 0 0 1 0 0 0 0 0 0 0 ...
..- attr(*, "label")= chr "Hypoxia before / at admission"
..- attr(*, "format.spss")= chr "F2.0"
$ hypotens: atomic 0 0 0 0 0 0 0 0 0 0 ...
..- attr(*, "label")= chr "Hypotension before / at admission"
..- attr(*, "format.spss")= chr "F2.0"
$ ctclass : atomic 2 2 4 2 2 2 2 2 2 1 ...
..- attr(*, "label")= chr "CT classification according to Marshall"
..- attr(*, "format.spss")= chr "F2.0"
$ tsah : atomic 0 0 1 0 0 NA 1 1 NA 0 ...
..- attr(*, "label")= chr "tSAH at CT"
..- attr(*, "format.spss")= chr "F2.0"
$ edh : atomic 0 0 0 0 0 0 0 0 0 0 ...
..- attr(*, "label")= chr "EDH at CT"
..- attr(*, "format.spss")= chr "F2.0"
$ cisterns: atomic 1 1 NA 1 1 1 1 2 1 1 ...
..- attr(*, "label")= chr "Compressed cisterns at CT"
..- attr(*, "format.spss")= chr "F2.0"
$ shift : atomic 0 0 NA 1 0 0 0 1 0 0 ...
..- attr(*, "label")= chr "Midline shift > 5 mm at CT"
..- attr(*, "format.spss")= chr "F2.0"
$ d.sysbpt: atomic 119 131 136 137 138 ...
..- attr(*, "label")= chr "Systolic blood pressure (truncated, mm Hg)"
..- attr(*, "format.spss")= chr "F4.0"
$ glucose : atomic 7.7 7.4 7.56 5.7 6 ...
..- attr(*, "label")= chr "Glucose at admission (mmol/l)"
..- attr(*, "format.spss")= chr "F4.2"
$ glucoset: atomic 7.7 7.4 7.56 5.7 6 ...
..- attr(*, "label")= chr "Truncated glucose values"
..- attr(*, "format.spss")= chr "F4.2"
$ ph : atomic 7.35 7.33 7.49 NA NA ...
..- attr(*, "label")= chr "pH"
..- attr(*, "format.spss")= chr "F4.2"
$ sodium : atomic 143 143 141 144 142 138 141 135 137 137 ...
..- attr(*, "label")= chr "Sodium (mmol/l)"
..- attr(*, "format.spss")= chr "F4.1"
$ sodiumt : atomic 143 143 141 144 142 138 141 135 137 137 ...
..- attr(*, "label")= chr "Truncated sodium"
..- attr(*, "format.spss")= chr "F4.1"
$ hb : atomic 15 15 12.8 7.1 8.5 12.3 12.8 9.2 12.8 11 ...
..- attr(*, "label")= chr "hb (g/dl)"
..- attr(*, "format.spss")= chr "F4.1"
$ hbt : atomic 15 15 12.8 7.1 8.5 12.3 12.8 9.2 12.8 11 ...
..- attr(*, "label")= chr "Truncated hb"
..- attr(*, "format.spss")= chr "F4.1"
Give the frequencies of the outcome (d.gos). What is the most commonly observed outcome?
tabl(tbi$d.gos)
dead vegetative severe disability
503 86 262
moderate disability good recovery <NA>
351 957 0
Check the categorization in favorable vs unfavorable outcome (d.unfav variable). What is the overall risk of an unfavorable outcome? How was d.gos dichotomized?
tabl(tbi$d.gos, tbi$d.unfav)
0 1 <NA>
dead 0 503 0
vegetative 0 86 0
severe disability 0 262 0
moderate disability 351 0 0
good recovery 957 0 0
<NA> 0 0 0
Overall risk of unfavorable outcome:
mean(tbi$d.unfav)
[1] 0.394164
Give the frequencies of cause of injury. What is the most common cause of injury?
tabl(tbi$cause)
Road traffic accident Motorbike Assault
848 420 134
domestic/fall other <NA>
370 387 0
Study the univariable effect of the prognostic factor ‘cause of injury’ on ‘unfavorable outcome’ (with crosstabs). What is the risk of an unfavorable outcome for each cause of injury? (use option
)
gmodels::CrossTable(tbi$cause, tbi$d.unfav, prop.chisq = F, prop.t = F, prop.c = F)
Cell Contents
|-------------------------|
| N |
| N / Row Total |
|-------------------------|
Total Observations in Table: 2159
| tbi$d.unfav
tbi$cause | 0 | 1 | Row Total |
----------------------|-----------|-----------|-----------|
Road traffic accident | 530 | 318 | 848 |
| 0.625 | 0.375 | 0.393 |
----------------------|-----------|-----------|-----------|
Motorbike | 277 | 143 | 420 |
| 0.660 | 0.340 | 0.195 |
----------------------|-----------|-----------|-----------|
Assault | 87 | 47 | 134 |
| 0.649 | 0.351 | 0.062 |
----------------------|-----------|-----------|-----------|
domestic/fall | 200 | 170 | 370 |
| 0.541 | 0.459 | 0.171 |
----------------------|-----------|-----------|-----------|
other | 214 | 173 | 387 |
| 0.553 | 0.447 | 0.179 |
----------------------|-----------|-----------|-----------|
Column Total | 1308 | 851 | 2159 |
----------------------|-----------|-----------|-----------|
Quantify the effect with a logistic regression model.
Specify cause of injury as a categorical covariate.
Let’s make sure that ‘other’ is the reference category
tbi %<>% mutate(cause = relevel(cause, ref = "other"))
fit <- glm(d.unfav ~ cause, data = tbi, family = binomial("logit"))
summary(fit)
Call:
glm(formula = d.unfav ~ cause, family = binomial("logit"), data = tbi)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.1092 -0.9695 -0.9124 1.2690 1.4679
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.21268 0.10224 -2.080 0.0375 *
causeRoad traffic accident -0.29814 0.12444 -2.396 0.0166 *
causeMotorbike -0.44849 0.14511 -3.091 0.0020 **
causeAssault -0.40308 0.20790 -1.939 0.0525 .
causedomestic/fall 0.05017 0.14607 0.343 0.7313
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 2895.5 on 2158 degrees of freedom
Residual deviance: 2877.0 on 2154 degrees of freedom
AIC: 2887
Number of Fisher Scoring iterations: 4
By default SPSS will use the last category (‘Other’) as the reference category. Which causes give the highest risk of unfavorable outcome, and which the lowest risk, according to the regression result?
Motorbike is lowest, domestic/fall is highest. These match with the cross-table. Let’s look at the predicted probabilities for each category
cbind(levels(tbi$cause), predict(fit, newdata = data.frame(cause = levels(tbi$cause)),
type = "response"))
[,1] [,2]
1 "other" "0.447028423772611"
2 "Road traffic accident" "0.375000000000002"
3 "Motorbike" "0.340476190476219"
4 "Assault" "0.350746268656731"
5 "domestic/fall" "0.459459459459461"
Verify that the intercept estimate in e) corresponds to the risk for the “Other” cause category as noted in d). Is the exp(intercept) an “odds ratio” or simply an “odds”?
With our fit the reference category is other
o_int = exp(coef(fit)[1])
o_int
(Intercept)
0.8084112
To probability
o_int / (1 + o_int)
(Intercept)
0.4470284
This matches the observed probability for other.
This is an actual ‘odds’
Verify also that the risk for the category “Domestic/fall” is only slightly higher than that of the “Other” cause category, both according to the crosstable (d)) and the regression result in e).
Is the pattern of risk in d) and e) what you would expect? Can you think of confounders? Hint: What is the mean age for each cause of injury?
You would expect the risk for traffic accidents to be higher. Let’s include age in the summary
tbi %>%
group_by(cause) %>%
summarize(mean_age = mean(age), prop_unfavouroble = mean(d.unfav))
# A tibble: 5 x 3
cause mean_age prop_unfavouroble
<fct> <dbl> <dbl>
1 other 35.1 0.447
2 Road traffic accident 29.5 0.375
3 Motorbike 31.4 0.340
4 Assault 35.3 0.351
5 domestic/fall 41.0 0.459
Motorbike and traffic have low age and low unfavourable outcomes
Now fit a multivariable model to adjust the effect of cause of injury for age.
fit2 <- glm(d.unfav ~ cause + age, data = tbi, family = binomial("logit"))
summary(fit2)
Call:
glm(formula = d.unfav ~ cause + age, family = binomial("logit"),
data = tbi)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.6204 -0.9600 -0.8317 1.2374 1.7103
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -1.318676 0.162146 -8.133 4.2e-16 ***
causeRoad traffic accident -0.129030 0.128343 -1.005 0.3147
causeMotorbike -0.350214 0.148844 -2.353 0.0186 *
causeAssault -0.421374 0.211557 -1.992 0.0464 *
causedomestic/fall -0.137041 0.151044 -0.907 0.3643
age 0.031325 0.003498 8.955 < 2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 2895.5 on 2158 degrees of freedom
Residual deviance: 2794.1 on 2153 degrees of freedom
AIC: 2806.1
Number of Fisher Scoring iterations: 4
This model did not fit, the residual deviance is higher than the degrees of freedom
Let’s look if non-linear transformations of age are in place
Hmisc::rcspline.plot(x = tbi$age, y = tbi$d.unfav,
model = "logistic", adj = tbi %>% select(hbt, shift))
cbind(xe, lower, upper)
xe
[1,] 14.00000 -9.444175e-02 1.186191
[2,] 14.09349 -9.071695e-02 1.181569
[3,] 14.18698 -8.704456e-02 1.176999
[4,] 14.28047 -8.342560e-02 1.172483
[5,] 14.37396 -7.986113e-02 1.168021
[6,] 14.46745 -7.635220e-02 1.163615
[7,] 14.56093 -7.289988e-02 1.159265
[8,] 14.65442 -6.950526e-02 1.154974
[9,] 14.74791 -6.616942e-02 1.150741
[10,] 14.84140 -6.289347e-02 1.146567
[11,] 14.93489 -5.967853e-02 1.142455
[12,] 15.02838 -5.652571e-02 1.138405
[13,] 15.12187 -5.343615e-02 1.134418
[14,] 15.21536 -5.041099e-02 1.130496
[15,] 15.30885 -4.745136e-02 1.126639
[16,] 15.40234 -4.455841e-02 1.122849
[17,] 15.49583 -4.173330e-02 1.119127
[18,] 15.58932 -3.897719e-02 1.115473
[19,] 15.68280 -3.629121e-02 1.111890
[20,] 15.77629 -3.367654e-02 1.108378
[21,] 15.86978 -3.113433e-02 1.104939
[22,] 15.96327 -2.866571e-02 1.101573
[23,] 16.05676 -2.627185e-02 1.098282
[24,] 16.15025 -2.395388e-02 1.095067
[25,] 16.24374 -2.171293e-02 1.091928
[26,] 16.33723 -1.955013e-02 1.088868
[27,] 16.43072 -1.746657e-02 1.085888
[28,] 16.52421 -1.546337e-02 1.082987
[29,] 16.61770 -1.354159e-02 1.080168
[30,] 16.71119 -1.170230e-02 1.077432
[31,] 16.80467 -9.946554e-03 1.074779
[32,] 16.89816 -8.275365e-03 1.072210
[33,] 16.99165 -6.689733e-03 1.069727
[34,] 17.08514 -5.190610e-03 1.067331
[35,] 17.17863 -3.778756e-03 1.065024
[36,] 17.27212 -2.454739e-03 1.062807
[37,] 17.36561 -1.218951e-03 1.060683
[38,] 17.45910 -7.160723e-05 1.058653
[39,] 17.55259 9.872584e-04 1.056720
[40,] 17.64608 1.957799e-03 1.054884
[41,] 17.73957 2.840355e-03 1.053146
[42,] 17.83306 3.635458e-03 1.051508
[43,] 17.92654 4.343830e-03 1.049971
[44,] 18.02003 4.966390e-03 1.048534
[45,] 18.11352 5.504246e-03 1.047199
[46,] 18.20701 5.958700e-03 1.045966
[47,] 18.30050 6.331246e-03 1.044834
[48,] 18.39399 6.623568e-03 1.043804
[49,] 18.48748 6.837537e-03 1.042876
[50,] 18.58097 6.975207e-03 1.042048
[51,] 18.67446 7.038814e-03 1.041320
[52,] 18.76795 7.030769e-03 1.040691
[53,] 18.86144 6.953652e-03 1.040161
[54,] 18.95492 6.810212e-03 1.039727
[55,] 19.04841 6.603353e-03 1.039389
[56,] 19.14190 6.336134e-03 1.039145
[57,] 19.23539 6.011758e-03 1.038993
[58,] 19.32888 5.633568e-03 1.038932
[59,] 19.42237 5.205033e-03 1.038959
[60,] 19.51586 4.729750e-03 1.039073
[61,] 19.60935 4.211424e-03 1.039271
[62,] 19.70284 3.653870e-03 1.039552
[63,] 19.79633 3.060997e-03 1.039912
[64,] 19.88982 2.436802e-03 1.040350
[65,] 19.98331 1.785363e-03 1.040862
[66,] 20.07679 1.110826e-03 1.041446
[67,] 20.17028 4.173979e-04 1.042100
[68,] 20.26377 -2.906600e-04 1.042821
[69,] 20.35726 -1.009045e-03 1.043605
[70,] 20.45075 -1.733417e-03 1.044451
[71,] 20.54424 -2.459414e-03 1.045355
[72,] 20.63773 -3.182651e-03 1.046314
[73,] 20.73122 -3.898739e-03 1.047327
[74,] 20.82471 -4.603284e-03 1.048388
[75,] 20.91820 -5.291901e-03 1.049497
[76,] 21.01169 -5.960219e-03 1.050650
[77,] 21.10518 -6.603889e-03 1.051844
[78,] 21.19866 -7.218593e-03 1.053077
[79,] 21.29215 -7.800051e-03 1.054346
[80,] 21.38564 -8.344026e-03 1.055647
[81,] 21.47913 -8.846334e-03 1.056979
[82,] 21.57262 -9.302848e-03 1.058339
[83,] 21.66611 -9.709510e-03 1.059724
[84,] 21.75960 -1.006233e-02 1.061132
[85,] 21.85309 -1.035740e-02 1.062560
[86,] 21.94658 -1.059090e-02 1.064007
[87,] 22.04007 -1.075910e-02 1.065470
[88,] 22.13356 -1.085836e-02 1.066946
[89,] 22.22705 -1.088516e-02 1.068434
[90,] 22.32053 -1.083609e-02 1.069933
[91,] 22.41402 -1.070784e-02 1.071439
[92,] 22.50751 -1.049726e-02 1.072953
[93,] 22.60100 -1.020129e-02 1.074471
[94,] 22.69449 -9.817052e-03 1.075994
[95,] 22.78798 -9.341786e-03 1.077518
[96,] 22.88147 -8.772897e-03 1.079045
[97,] 22.97496 -8.107951e-03 1.080571
[98,] 23.06845 -7.344900e-03 1.082098
[99,] 23.16194 -6.483891e-03 1.083625
[100,] 23.25543 -5.526330e-03 1.085153
[101,] 23.34891 -4.473794e-03 1.086681
[102,] 23.44240 -3.328010e-03 1.088212
[103,] 23.53589 -2.090839e-03 1.089747
[104,] 23.62938 -7.642641e-04 1.091285
[105,] 23.72287 6.496148e-04 1.092829
[106,] 23.81636 2.148597e-03 1.094381
[107,] 23.90985 3.730390e-03 1.095940
[108,] 24.00334 5.392617e-03 1.097510
[109,] 24.09683 7.132829e-03 1.099091
[110,] 24.19032 8.948515e-03 1.100685
[111,] 24.28381 1.083711e-02 1.102295
[112,] 24.37730 1.279599e-02 1.103921
[113,] 24.47078 1.482252e-02 1.105566
[114,] 24.56427 1.691402e-02 1.107231
[115,] 24.65776 1.906778e-02 1.108918
[116,] 24.75125 2.128109e-02 1.110629
[117,] 24.84474 2.355124e-02 1.112367
[118,] 24.93823 2.587551e-02 1.114132
[119,] 25.03172 2.825120e-02 1.115927
[120,] 25.12521 3.067563e-02 1.117753
[121,] 25.21870 3.314612e-02 1.119613
[122,] 25.31219 3.566004e-02 1.121508
[123,] 25.40568 3.821481e-02 1.123440
[124,] 25.49917 4.080788e-02 1.125410
[125,] 25.59265 4.343674e-02 1.127421
[126,] 25.68614 4.609896e-02 1.129473
[127,] 25.77963 4.879216e-02 1.131569
[128,] 25.87312 5.151403e-02 1.133710
[129,] 25.96661 5.426234e-02 1.135897
[130,] 26.06010 5.703493e-02 1.138132
[131,] 26.15359 5.982974e-02 1.140415
[132,] 26.24708 6.264478e-02 1.142749
[133,] 26.34057 6.547817e-02 1.145133
[134,] 26.43406 6.832813e-02 1.147569
[135,] 26.52755 7.119297e-02 1.150058
[136,] 26.62104 7.407110e-02 1.152601
[137,] 26.71452 7.696108e-02 1.155197
[138,] 26.80801 7.986154e-02 1.157848
[139,] 26.90150 8.277124e-02 1.160554
[140,] 26.99499 8.568906e-02 1.163316
[141,] 27.08848 8.861400e-02 1.166133
[142,] 27.18197 9.154518e-02 1.169005
[143,] 27.27546 9.448184e-02 1.171933
[144,] 27.36895 9.742336e-02 1.174916
[145,] 27.46244 1.003692e-01 1.177954
[146,] 27.55593 1.033190e-01 1.181046
[147,] 27.64942 1.062725e-01 1.184192
[148,] 27.74290 1.092296e-01 1.187391
[149,] 27.83639 1.121902e-01 1.190643
[150,] 27.92988 1.151545e-01 1.193945
[151,] 28.02337 1.181227e-01 1.197298
[152,] 28.11686 1.210951e-01 1.200700
[153,] 28.21035 1.240722e-01 1.204150
[154,] 28.30384 1.270545e-01 1.207646
[155,] 28.39733 1.300428e-01 1.211187
[156,] 28.49082 1.330378e-01 1.214771
[157,] 28.58431 1.360404e-01 1.218397
[158,] 28.67780 1.390516e-01 1.222062
[159,] 28.77129 1.420725e-01 1.225765
[160,] 28.86477 1.451042e-01 1.229503
[161,] 28.95826 1.481481e-01 1.233276
[162,] 29.05175 1.512053e-01 1.237079
[163,] 29.14524 1.542774e-01 1.240912
[164,] 29.23873 1.573658e-01 1.244772
[165,] 29.33222 1.604720e-01 1.248656
[166,] 29.42571 1.635977e-01 1.252562
[167,] 29.51920 1.667444e-01 1.256488
[168,] 29.61269 1.699138e-01 1.260431
[169,] 29.70618 1.731077e-01 1.264388
[170,] 29.79967 1.763278e-01 1.268358
[171,] 29.89316 1.795759e-01 1.272336
[172,] 29.98664 1.828537e-01 1.276321
[173,] 30.08013 1.861629e-01 1.280311
[174,] 30.17362 1.895047e-01 1.284303
[175,] 30.26711 1.928797e-01 1.288296
[176,] 30.36060 1.962887e-01 1.292288
[177,] 30.45409 1.997322e-01 1.296280
[178,] 30.54758 2.032108e-01 1.300269
[179,] 30.64107 2.067247e-01 1.304254
[180,] 30.73456 2.102744e-01 1.308236
[181,] 30.82805 2.138601e-01 1.312212
[182,] 30.92154 2.174820e-01 1.316182
[183,] 31.01503 2.211401e-01 1.320146
[184,] 31.10851 2.248345e-01 1.324103
[185,] 31.20200 2.285651e-01 1.328052
[186,] 31.29549 2.323318e-01 1.331993
[187,] 31.38898 2.361345e-01 1.335925
[188,] 31.48247 2.399730e-01 1.339848
[189,] 31.57596 2.438469e-01 1.343762
[190,] 31.66945 2.477559e-01 1.347667
[191,] 31.76294 2.516996e-01 1.351562
[192,] 31.85643 2.556777e-01 1.355447
[193,] 31.94992 2.596895e-01 1.359322
[194,] 32.04341 2.637346e-01 1.363187
[195,] 32.13689 2.678124e-01 1.367042
[196,] 32.23038 2.719223e-01 1.370887
[197,] 32.32387 2.760636e-01 1.374722
[198,] 32.41736 2.802356e-01 1.378547
[199,] 32.51085 2.844377e-01 1.382363
[200,] 32.60434 2.886691e-01 1.386169
[201,] 32.69783 2.929289e-01 1.389966
[202,] 32.79132 2.972164e-01 1.393754
[203,] 32.88481 3.015307e-01 1.397533
[204,] 32.97830 3.058709e-01 1.401304
[205,] 33.07179 3.102362e-01 1.405067
[206,] 33.16528 3.146256e-01 1.408821
[207,] 33.25876 3.190382e-01 1.412568
[208,] 33.35225 3.234731e-01 1.416308
[209,] 33.44574 3.279292e-01 1.420042
[210,] 33.53923 3.324055e-01 1.423769
[211,] 33.63272 3.369012e-01 1.427490
[212,] 33.72621 3.414150e-01 1.431205
[213,] 33.81970 3.459461e-01 1.434916
[214,] 33.91319 3.504933e-01 1.438622
[215,] 34.00668 3.550556e-01 1.442324
[216,] 34.10017 3.596320e-01 1.446023
[217,] 34.19366 3.642214e-01 1.449718
[218,] 34.28715 3.688226e-01 1.453410
[219,] 34.38063 3.734347e-01 1.457101
[220,] 34.47412 3.780566e-01 1.460790
[221,] 34.56761 3.826872e-01 1.464478
[222,] 34.66110 3.873253e-01 1.468165
[223,] 34.75459 3.919700e-01 1.471851
[224,] 34.84808 3.966201e-01 1.475539
[225,] 34.94157 4.012746e-01 1.479227
[226,] 35.03506 4.059324e-01 1.482916
[227,] 35.12855 4.105924e-01 1.486607
[228,] 35.22204 4.152536e-01 1.490301
[229,] 35.31553 4.199150e-01 1.493997
[230,] 35.40902 4.245756e-01 1.497696
[231,] 35.50250 4.292342e-01 1.501399
[232,] 35.59599 4.338899e-01 1.505106
[233,] 35.68948 4.385416e-01 1.508817
[234,] 35.78297 4.431885e-01 1.512533
[235,] 35.87646 4.478295e-01 1.516254
[236,] 35.96995 4.524637e-01 1.519981
[237,] 36.06344 4.570901e-01 1.523714
[238,] 36.15693 4.617079e-01 1.527453
[239,] 36.25042 4.663160e-01 1.531199
[240,] 36.34391 4.709137e-01 1.534952
[241,] 36.43740 4.755001e-01 1.538712
[242,] 36.53088 4.800743e-01 1.542479
[243,] 36.62437 4.846356e-01 1.546254
[244,] 36.71786 4.891831e-01 1.550037
[245,] 36.81135 4.937160e-01 1.553828
[246,] 36.90484 4.982337e-01 1.557627
[247,] 36.99833 5.027353e-01 1.561435
[248,] 37.09182 5.072203e-01 1.565251
[249,] 37.18531 5.116880e-01 1.569075
[250,] 37.27880 5.161376e-01 1.572908
[251,] 37.37229 5.205686e-01 1.576750
[252,] 37.46578 5.249804e-01 1.580600
[253,] 37.55927 5.293725e-01 1.584459
[254,] 37.65275 5.337443e-01 1.588327
[255,] 37.74624 5.380954e-01 1.592203
[256,] 37.83973 5.424252e-01 1.596087
[257,] 37.93322 5.467333e-01 1.599980
[258,] 38.02671 5.510193e-01 1.603880
[259,] 38.12020 5.552829e-01 1.607788
[260,] 38.21369 5.595237e-01 1.611704
[261,] 38.30718 5.637414e-01 1.615628
[262,] 38.40067 5.679356e-01 1.619558
[263,] 38.49416 5.721062e-01 1.623495
[264,] 38.58765 5.762529e-01 1.627438
[265,] 38.68114 5.803755e-01 1.631387
[266,] 38.77462 5.844738e-01 1.635342
[267,] 38.86811 5.885478e-01 1.639302
[268,] 38.96160 5.925972e-01 1.643266
[269,] 39.05509 5.966221e-01 1.647235
[270,] 39.14858 6.006224e-01 1.651207
[271,] 39.24207 6.045980e-01 1.655182
[272,] 39.33556 6.085490e-01 1.659159
[273,] 39.42905 6.124753e-01 1.663138
[274,] 39.52254 6.163772e-01 1.667119
[275,] 39.61603 6.202546e-01 1.671100
[276,] 39.70952 6.241076e-01 1.675080
[277,] 39.80301 6.279365e-01 1.679060
[278,] 39.89649 6.317413e-01 1.683038
[279,] 39.98998 6.355222e-01 1.687013
[280,] 40.08347 6.392795e-01 1.690986
[281,] 40.17696 6.430134e-01 1.694954
[282,] 40.27045 6.467242e-01 1.698917
[283,] 40.36394 6.504120e-01 1.702875
[284,] 40.45743 6.540773e-01 1.706825
[285,] 40.55092 6.577202e-01 1.710769
[286,] 40.64441 6.613412e-01 1.714704
[287,] 40.73790 6.649406e-01 1.718629
[288,] 40.83139 6.685186e-01 1.722544
[289,] 40.92487 6.720758e-01 1.726448
[290,] 41.01836 6.756125e-01 1.730340
[291,] 41.11185 6.791291e-01 1.734219
[292,] 41.20534 6.826262e-01 1.738084
[293,] 41.29883 6.861044e-01 1.741936
[294,] 41.39232 6.895643e-01 1.745773
[295,] 41.48581 6.930064e-01 1.749595
[296,] 41.57930 6.964312e-01 1.753402
[297,] 41.67279 6.998393e-01 1.757194
[298,] 41.76628 7.032312e-01 1.760970
[299,] 41.85977 7.066072e-01 1.764730
[300,] 41.95326 7.099680e-01 1.768474
[301,] 42.04674 7.133140e-01 1.772201
[302,] 42.14023 7.166456e-01 1.775911
[303,] 42.23372 7.199632e-01 1.779605
[304,] 42.32721 7.232672e-01 1.783280
[305,] 42.42070 7.265581e-01 1.786939
[306,] 42.51419 7.298362e-01 1.790579
[307,] 42.60768 7.331020e-01 1.794201
[308,] 42.70117 7.363556e-01 1.797805
[309,] 42.79466 7.395976e-01 1.801391
[310,] 42.88815 7.428282e-01 1.804959
[311,] 42.98164 7.460478e-01 1.808507
[312,] 43.07513 7.492566e-01 1.812037
[313,] 43.16861 7.524549e-01 1.815548
[314,] 43.26210 7.556431e-01 1.819040
[315,] 43.35559 7.588213e-01 1.822512
[316,] 43.44908 7.619900e-01 1.825966
[317,] 43.54257 7.651492e-01 1.829400
[318,] 43.63606 7.682993e-01 1.832815
[319,] 43.72955 7.714404e-01 1.836210
[320,] 43.82304 7.745728e-01 1.839586
[321,] 43.91653 7.776967e-01 1.842943
[322,] 44.01002 7.808123e-01 1.846280
[323,] 44.10351 7.839198e-01 1.849597
[324,] 44.19699 7.870193e-01 1.852895
[325,] 44.29048 7.901109e-01 1.856174
[326,] 44.38397 7.931950e-01 1.859433
[327,] 44.47746 7.962715e-01 1.862672
[328,] 44.57095 7.993407e-01 1.865892
[329,] 44.66444 8.024027e-01 1.869093
[330,] 44.75793 8.054575e-01 1.872274
[331,] 44.85142 8.085054e-01 1.875436
[332,] 44.94491 8.115463e-01 1.878579
[333,] 45.03840 8.145805e-01 1.881703
[334,] 45.13189 8.176079e-01 1.884808
[335,] 45.22538 8.206287e-01 1.887894
[336,] 45.31886 8.236429e-01 1.890961
[337,] 45.41235 8.266507e-01 1.894009
[338,] 45.50584 8.296519e-01 1.897038
[339,] 45.59933 8.326468e-01 1.900049
[340,] 45.69282 8.356353e-01 1.903042
[341,] 45.78631 8.386175e-01 1.906017
[342,] 45.87980 8.415935e-01 1.908973
[343,] 45.97329 8.445631e-01 1.911911
[344,] 46.06678 8.475265e-01 1.914832
[345,] 46.16027 8.504837e-01 1.917735
[346,] 46.25376 8.534346e-01 1.920621
[347,] 46.34725 8.563793e-01 1.923489
[348,] 46.44073 8.593177e-01 1.926340
[349,] 46.53422 8.622499e-01 1.929175
[350,] 46.62771 8.651757e-01 1.931992
[351,] 46.72120 8.680953e-01 1.934793
[352,] 46.81469 8.710086e-01 1.937578
[353,] 46.90818 8.739154e-01 1.940347
[354,] 47.00167 8.768159e-01 1.943100
[355,] 47.09516 8.797099e-01 1.945837
[356,] 47.18865 8.825974e-01 1.948559
[357,] 47.28214 8.854783e-01 1.951265
[358,] 47.37563 8.883527e-01 1.953956
[359,] 47.46912 8.912203e-01 1.956633
[360,] 47.56260 8.940812e-01 1.959295
[361,] 47.65609 8.969353e-01 1.961943
[362,] 47.74958 8.997825e-01 1.964577
[363,] 47.84307 9.026227e-01 1.967197
[364,] 47.93656 9.054558e-01 1.969804
[365,] 48.03005 9.082818e-01 1.972397
[366,] 48.12354 9.111006e-01 1.974977
[367,] 48.21703 9.139121e-01 1.977545
[368,] 48.31052 9.167161e-01 1.980100
[369,] 48.40401 9.195126e-01 1.982642
[370,] 48.49750 9.223015e-01 1.985173
[371,] 48.59098 9.250827e-01 1.987693
[372,] 48.68447 9.278560e-01 1.990201
[373,] 48.77796 9.306213e-01 1.992697
[374,] 48.87145 9.333787e-01 1.995183
[375,] 48.96494 9.361278e-01 1.997659
[376,] 49.05843 9.388686e-01 2.000124
[377,] 49.15192 9.416010e-01 2.002580
[378,] 49.24541 9.443249e-01 2.005026
[379,] 49.33890 9.470401e-01 2.007462
[380,] 49.43239 9.497465e-01 2.009889
[381,] 49.52588 9.524440e-01 2.012308
[382,] 49.61937 9.551325e-01 2.014718
[383,] 49.71285 9.578118e-01 2.017120
[384,] 49.80634 9.604817e-01 2.019515
[385,] 49.89983 9.631422e-01 2.021902
[386,] 49.99332 9.657932e-01 2.024281
[387,] 50.08681 9.684344e-01 2.026654
[388,] 50.18030 9.710657e-01 2.029020
[389,] 50.27379 9.736871e-01 2.031380
[390,] 50.36728 9.762983e-01 2.033734
[391,] 50.46077 9.788993e-01 2.036082
[392,] 50.55426 9.814898e-01 2.038425
[393,] 50.64775 9.840698e-01 2.040763
[394,] 50.74124 9.866391e-01 2.043096
[395,] 50.83472 9.891975e-01 2.045425
[396,] 50.92821 9.917450e-01 2.047749
[397,] 51.02170 9.942813e-01 2.050070
[398,] 51.11519 9.968064e-01 2.052388
[399,] 51.20868 9.993201e-01 2.054702
[400,] 51.30217 1.001822e+00 2.057013
[401,] 51.39566 1.004313e+00 2.059322
[402,] 51.48915 1.006791e+00 2.061629
[403,] 51.58264 1.009258e+00 2.063934
[404,] 51.67613 1.011713e+00 2.066237
[405,] 51.76962 1.014155e+00 2.068539
[406,] 51.86311 1.016585e+00 2.070840
[407,] 51.95659 1.019002e+00 2.073140
[408,] 52.05008 1.021407e+00 2.075440
[409,] 52.14357 1.023799e+00 2.077740
[410,] 52.23706 1.026179e+00 2.080041
[411,] 52.33055 1.028545e+00 2.082341
[412,] 52.42404 1.030898e+00 2.084643
[413,] 52.51753 1.033237e+00 2.086946
[414,] 52.61102 1.035564e+00 2.089250
[415,] 52.70451 1.037877e+00 2.091556
[416,] 52.79800 1.040176e+00 2.093864
[417,] 52.89149 1.042461e+00 2.096174
[418,] 52.98497 1.044733e+00 2.098487
[419,] 53.07846 1.046990e+00 2.100803
[420,] 53.17195 1.049234e+00 2.103122
[421,] 53.26544 1.051463e+00 2.105444
[422,] 53.35893 1.053678e+00 2.107770
[423,] 53.45242 1.055878e+00 2.110100
[424,] 53.54591 1.058064e+00 2.112435
[425,] 53.63940 1.060236e+00 2.114774
[426,] 53.73289 1.062392e+00 2.117117
[427,] 53.82638 1.064534e+00 2.119466
[428,] 53.91987 1.066661e+00 2.121820
[429,] 54.01336 1.068774e+00 2.124180
[430,] 54.10684 1.070871e+00 2.126545
[431,] 54.20033 1.072953e+00 2.128917
[432,] 54.29382 1.075020e+00 2.131295
[433,] 54.38731 1.077071e+00 2.133679
[434,] 54.48080 1.079107e+00 2.136071
[435,] 54.57429 1.081128e+00 2.138469
[436,] 54.66778 1.083134e+00 2.140874
[437,] 54.76127 1.085124e+00 2.143287
[438,] 54.85476 1.087099e+00 2.145708
[439,] 54.94825 1.089058e+00 2.148137
[440,] 55.04174 1.091001e+00 2.150573
[441,] 55.13523 1.092929e+00 2.153018
[442,] 55.22871 1.094841e+00 2.155472
[443,] 55.32220 1.096737e+00 2.157934
[444,] 55.41569 1.098618e+00 2.160405
[445,] 55.50918 1.100483e+00 2.162886
[446,] 55.60267 1.102333e+00 2.165375
[447,] 55.69616 1.104166e+00 2.167874
[448,] 55.78965 1.105984e+00 2.170383
[449,] 55.88314 1.107786e+00 2.172901
[450,] 55.97663 1.109573e+00 2.175430
[451,] 56.07012 1.111344e+00 2.177968
[452,] 56.16361 1.113099e+00 2.180517
[453,] 56.25710 1.114838e+00 2.183076
[454,] 56.35058 1.116562e+00 2.185646
[455,] 56.44407 1.118270e+00 2.188227
[456,] 56.53756 1.119963e+00 2.190818
[457,] 56.63105 1.121640e+00 2.193420
[458,] 56.72454 1.123301e+00 2.196034
[459,] 56.81803 1.124947e+00 2.198658
[460,] 56.91152 1.126578e+00 2.201294
[461,] 57.00501 1.128194e+00 2.203942
[462,] 57.09850 1.129794e+00 2.206601
[463,] 57.19199 1.131380e+00 2.209271
[464,] 57.28548 1.132950e+00 2.211954
[465,] 57.37896 1.134505e+00 2.214648
[466,] 57.47245 1.136045e+00 2.217354
[467,] 57.56594 1.137571e+00 2.220072
[468,] 57.65943 1.139082e+00 2.222802
[469,] 57.75292 1.140578e+00 2.225544
[470,] 57.84641 1.142060e+00 2.228298
[471,] 57.93990 1.143527e+00 2.231065
[472,] 58.03339 1.144980e+00 2.233844
[473,] 58.12688 1.146419e+00 2.236635
[474,] 58.22037 1.147844e+00 2.239438
[475,] 58.31386 1.149256e+00 2.242254
[476,] 58.40735 1.150653e+00 2.245082
[477,] 58.50083 1.152037e+00 2.247923
[478,] 58.59432 1.153408e+00 2.250776
[479,] 58.68781 1.154765e+00 2.253642
[480,] 58.78130 1.156109e+00 2.256521
[481,] 58.87479 1.157441e+00 2.259411
[482,] 58.96828 1.158759e+00 2.262315
[483,] 59.06177 1.160066e+00 2.265230
[484,] 59.15526 1.161359e+00 2.268159
[485,] 59.24875 1.162640e+00 2.271099
[486,] 59.34224 1.163909e+00 2.274053
[487,] 59.43573 1.165165e+00 2.277018
[488,] 59.52922 1.166410e+00 2.279996
[489,] 59.62270 1.167642e+00 2.282985
[490,] 59.71619 1.168862e+00 2.285987
[491,] 59.80968 1.170069e+00 2.289001
[492,] 59.90317 1.171265e+00 2.292027
[493,] 59.99666 1.172449e+00 2.295066
[494,] 60.09015 1.173621e+00 2.298115
[495,] 60.18364 1.174781e+00 2.301177
[496,] 60.27713 1.175929e+00 2.304251
[497,] 60.37062 1.177066e+00 2.307336
[498,] 60.46411 1.178191e+00 2.310433
[499,] 60.55760 1.179304e+00 2.313542
[500,] 60.65109 1.180406e+00 2.316662
[501,] 60.74457 1.181496e+00 2.319793
[502,] 60.83806 1.182575e+00 2.322936
[503,] 60.93155 1.183643e+00 2.326091
[504,] 61.02504 1.184699e+00 2.329256
[505,] 61.11853 1.185744e+00 2.332433
[506,] 61.21202 1.186778e+00 2.335621
[507,] 61.30551 1.187800e+00 2.338820
[508,] 61.39900 1.188812e+00 2.342031
[509,] 61.49249 1.189813e+00 2.345252
[510,] 61.58598 1.190802e+00 2.348484
[511,] 61.67947 1.191781e+00 2.351727
[512,] 61.77295 1.192749e+00 2.354980
[513,] 61.86644 1.193707e+00 2.358245
[514,] 61.95993 1.194654e+00 2.361520
[515,] 62.05342 1.195590e+00 2.364806
[516,] 62.14691 1.196515e+00 2.368102
[517,] 62.24040 1.197430e+00 2.371409
[518,] 62.33389 1.198335e+00 2.374726
[519,] 62.42738 1.199230e+00 2.378053
[520,] 62.52087 1.200114e+00 2.381391
[521,] 62.61436 1.200988e+00 2.384739
[522,] 62.70785 1.201851e+00 2.388097
[523,] 62.80134 1.202705e+00 2.391466
[524,] 62.89482 1.203549e+00 2.394844
[525,] 62.98831 1.204382e+00 2.398232
[526,] 63.08180 1.205206e+00 2.401630
[527,] 63.17529 1.206020e+00 2.405038
[528,] 63.26878 1.206824e+00 2.408456
[529,] 63.36227 1.207618e+00 2.411883
[530,] 63.45576 1.208403e+00 2.415320
[531,] 63.54925 1.209178e+00 2.418767
[532,] 63.64274 1.209944e+00 2.422223
[533,] 63.73623 1.210700e+00 2.425689
[534,] 63.82972 1.211447e+00 2.429164
[535,] 63.92321 1.212185e+00 2.432648
[536,] 64.01669 1.212913e+00 2.436142
[537,] 64.11018 1.213632e+00 2.439644
[538,] 64.20367 1.214342e+00 2.443156
[539,] 64.29716 1.215043e+00 2.446677
[540,] 64.39065 1.215735e+00 2.450207
[541,] 64.48414 1.216418e+00 2.453746
[542,] 64.57763 1.217092e+00 2.457294
[543,] 64.67112 1.217757e+00 2.460851
[544,] 64.76461 1.218413e+00 2.464416
[545,] 64.85810 1.219061e+00 2.467990
[546,] 64.95159 1.219700e+00 2.471573
[547,] 65.04508 1.220331e+00 2.475165
[548,] 65.13856 1.220952e+00 2.478764
[549,] 65.23205 1.221566e+00 2.482373
[550,] 65.32554 1.222171e+00 2.485990
[551,] 65.41903 1.222768e+00 2.489615
[552,] 65.51252 1.223356e+00 2.493248
[553,] 65.60601 1.223936e+00 2.496890
[554,] 65.69950 1.224509e+00 2.500540
[555,] 65.79299 1.225073e+00 2.504197
[556,] 65.88648 1.225628e+00 2.507863
[557,] 65.97997 1.226176e+00 2.511537
[558,] 66.07346 1.226716e+00 2.515219
[559,] 66.16694 1.227249e+00 2.518909
[560,] 66.26043 1.227773e+00 2.522606
[561,] 66.35392 1.228289e+00 2.526312
[562,] 66.44741 1.228798e+00 2.530025
[563,] 66.54090 1.229299e+00 2.533745
[564,] 66.63439 1.229793e+00 2.537474
[565,] 66.72788 1.230279e+00 2.541209
[566,] 66.82137 1.230758e+00 2.544953
[567,] 66.91486 1.231229e+00 2.548703
[568,] 67.00835 1.231693e+00 2.552461
[569,] 67.10184 1.232149e+00 2.556227
[570,] 67.19533 1.232599e+00 2.559999
[571,] 67.28881 1.233041e+00 2.563779
[572,] 67.38230 1.233476e+00 2.567566
[573,] 67.47579 1.233903e+00 2.571360
[574,] 67.56928 1.234324e+00 2.575161
[575,] 67.66277 1.234738e+00 2.578969
[576,] 67.75626 1.235145e+00 2.582784
[577,] 67.84975 1.235545e+00 2.586606
[578,] 67.94324 1.235938e+00 2.590435
[579,] 68.03673 1.236325e+00 2.594270
[580,] 68.13022 1.236704e+00 2.598112
[581,] 68.22371 1.237077e+00 2.601961
[582,] 68.31720 1.237444e+00 2.605817
[583,] 68.41068 1.237803e+00 2.609679
[584,] 68.50417 1.238157e+00 2.613547
[585,] 68.59766 1.238503e+00 2.617422
[586,] 68.69115 1.238844e+00 2.621304
[587,] 68.78464 1.239178e+00 2.625192
[588,] 68.87813 1.239505e+00 2.629086
[589,] 68.97162 1.239827e+00 2.632986
[590,] 69.06511 1.240142e+00 2.636893
[591,] 69.15860 1.240451e+00 2.640806
[592,] 69.25209 1.240754e+00 2.644725
[593,] 69.34558 1.241051e+00 2.648650
[594,] 69.43907 1.241341e+00 2.652581
[595,] 69.53255 1.241626e+00 2.656518
[596,] 69.62604 1.241905e+00 2.660461
[597,] 69.71953 1.242178e+00 2.664410
[598,] 69.81302 1.242445e+00 2.668365
[599,] 69.90651 1.242706e+00 2.672326
[600,] 70.00000 1.242962e+00 2.676292
Is the effect of cause still statistically significant? Hint: focus on the overall p-value, based on a 4 df test.
anova(fit2, test = "Chisq")
Analysis of Deviance Table
Model: binomial, link: logit
Response: d.unfav
Terms added sequentially (first to last)
Df Deviance Resid. Df Resid. Dev Pr(>Chi)
NULL 2158 2895.5
cause 4 18.517 2154 2877.0 0.0009777 ***
age 1 82.926 2153 2794.1 < 2.2e-16 ***
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
Dropping cause will significantly decrease the goodness of fit, so yes.
How do the effects of the different causes change?
require(tidyr)
fits <- list(without_age = fit, with_age = fit2)
fits %>%
map_df(tidy, .id = "model") %>%
select(estimate, model, term) %>%
spread(key = c("model"), value = "estimate")
term with_age without_age
1 (Intercept) -1.31867601 -0.21268442
2 age 0.03132549 NA
3 causeAssault -0.42137365 -0.40307610
4 causedomestic/fall -0.13704111 0.05016549
5 causeMotorbike -0.35021397 -0.44848846
6 causeRoad traffic accident -0.12902965 -0.29814120
Assault and domestic become lower risk, motorbike and road traffic higher risk
What is your conclusion on the effect of “cause of injury”? Do you think “cause of injury” should be used for predictive purposes?
Based on these data, yes. However, when adding more covariates, this may change.
We will now develop a simple prediction model with three predictors: motor score, pupillary reactivity and age.
Give some descriptive statistics of motor score, pupillary reactivity and age.
tbi %>%
select(d.motor, d.pupil, age) %>%
summary()
d.motor d.pupil age
Min. :1.000 both reactive :1430 Min. :14.00
1st Qu.:3.000 one reactive : 279 1st Qu.:22.00
Median :4.000 no reactive pupils: 327 Median :30.00
Mean :3.991 NA's : 123 Mean :33.21
3rd Qu.:5.000 3rd Qu.:43.00
Max. :6.000 Max. :79.00
- Assess the univariable effects of motor score, pupillary reactivity and age on the outcome (d.unfav) with a logistic regression model.
terms <- c("d.motor", "d.pupil", "age")
fits_uni <- terms %>%
map(function(term) glm(reformulate(term, "d.unfav"),
data = tbi, family = "binomial"))
names(fits_uni) <- terms
fits_uni %>%
map_df(tidy)
term estimate std.error statistic
1 (Intercept) 2.27399375 0.177962699 12.777923
2 d.motor -0.68865430 0.044140229 -15.601512
3 (Intercept) -0.87071813 0.057980412 -15.017454
4 d.pupilone reactive 0.97834880 0.133192368 7.345382
5 d.pupilno reactive pupils 1.71947266 0.133912687 12.840252
6 (Intercept) -1.49482050 0.121854369 -12.267270
7 age 0.03161422 0.003328418 9.498274
p.value
1 2.178073e-37
2 7.109036e-55
3 5.643444e-51
4 2.051725e-13
5 9.755630e-38
6 1.357530e-34
7 2.133986e-21
What is the univariable effect of age on the risk of unfavorable outcome?
exp(coef(fits_uni[["age"]]))
(Intercept) age
0.2242889 1.0321193
What would be a good way to express the effect, using a linear scale? Hint: think of recoding age by decade.
tbi %<>%
mutate(age_cat = cut(age, breaks = seq(from = 10*floor(min(age / 10)),
to = 10*ceiling(max(age / 10)), by = 10)))
tbi %>%
glm(formula = d.unfav ~ age_cat, family = "binomial") %>%
summary()
Call:
glm(formula = d.unfav ~ age_cat, family = "binomial", data = .)
Deviance Residuals:
Min 1Q Median 3Q Max
-1.5829 -0.9244 -0.8640 1.1334 1.5273
Coefficients:
Estimate Std. Error z value Pr(>|z|)
(Intercept) -0.789162 0.106175 -7.433 1.06e-13 ***
age_cat(20,30] -0.003851 0.133811 -0.029 0.9770
age_cat(30,40] 0.313901 0.145060 2.164 0.0305 *
age_cat(40,50] 0.976200 0.155716 6.269 3.63e-10 ***
age_cat(50,60] 0.893522 0.174017 5.135 2.83e-07 ***
age_cat(60,70] 1.319790 0.253405 5.208 1.91e-07 ***
age_cat(70,80] 1.705453 0.843370 2.022 0.0432 *
---
Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
(Dispersion parameter for binomial family taken to be 1)
Null deviance: 2895.5 on 2158 degrees of freedom
Residual deviance: 2797.1 on 2152 degrees of freedom
AIC: 2811.1
Number of Fisher Scoring iterations: 4
If the effect of age were linear (on the log-odds scale), there would be a constant difference between each consecutive category
Alternatively, we could plot the log-odds per age category
logit <- function(x) log(x / (1-x))
tbi %>%
group_by(age_cat) %>%
mutate(p_unfav = mean(d.unfav),
p_unfav_lo = binom.confint_logical(d.unfav)$lower,
p_unfav_hi = binom.confint_logical(d.unfav)$upper,
logit_unfav = logit(p_unfav),
logit_unfav_lo = logit(p_unfav_lo),
logit_unfav_hi = logit(p_unfav_hi)
) %>%
ggplot(aes(x = age_cat, ymin = logit_unfav_lo, y = logit_unfav, ymax = logit_unfav_hi)) +
geom_errorbar()
Linearity does not look too bad on logit scale.
I’m not sure whether calculating a confidence interval for the proportion and then transforming with logit is the best way to go.
Now fit a multivariable model. Include in your model: motor score, pupillary reactivity and age (continuous). Note that there are missing values in the variable ‘d.pupil’, which have been filled in with a statistical imputation procedure in ‘pupil.i’. Perform the analyses twice: once with ‘pupillary reactivity’ including missing values (d.pupil) and once with missing values imputed (pupil.i). What are the numbers of patients in each analysis? Are there differences in prognostic effects?
fit_mis <- glm(d.unfav ~ d.motor + d.pupil + age, data = tbi, family = "binomial")
fit_imp <- glm(d.unfav ~ d.motor + pupil.i + age, data = tbi, family = "binomial")
fits <- list(with_missings = fit_mis, imputed = fit_imp)
fits %>%
map_df(tidy, .id = "model") %>%
select(model, term, estimate) %>%
spread(model, estimate)
term imputed with_missings
1 (Intercept) 0.35269097 0.43316634
2 age 0.03818688 0.03835085
3 d.motor -0.60344181 -0.62470224
4 d.pupilno reactive pupils NA 1.28566748
5 d.pupilone reactive NA 0.56684531
6 pupil.ino reactive pupils 1.27120434 NA
7 pupil.ione reactive 0.59098094 NA
The estimate for age stays the same, for motor is a little different.
Overall the estimates are pretty much the same
fits %>% map_df("df.null") + 1
with_missings imputed
1 2036 2159
Can we interpret the change in age coefficient from univariable analysis to multivariable analysis if the number of subjects between the two analyses differ? Therefore: which variable for ‘pupillary reactivity’ do you prefer for modeling? Use this as the final multivariable model.
The number of missings is relatively low compare to the total number of observations (around 5%). If the missings are random and / or not associated with age or the outcome, the coefficients would not change. Precision decreases a little bit. Best would be to use the imputed variable, but difference will be small.
fits <- list(univariate = fits_uni[["age"]],
with_missings = fit_mis, imputed = fit_imp)
fits %>%
map_df(tidy, .id = "model") %>%
select(term, model, estimate) %>%
spread(model, estimate)
term imputed univariate with_missings
1 (Intercept) 0.35269097 -1.49482050 0.43316634
2 age 0.03818688 0.03161422 0.03835085
3 d.motor -0.60344181 NA -0.62470224
4 d.pupilno reactive pupils NA NA 1.28566748
5 d.pupilone reactive NA NA 0.56684531
6 pupil.ino reactive pupils 1.27120434 NA NA
7 pupil.ione reactive 0.59098094 NA NA
And for the p-value (precision):
fits %>%
map_df(tidy, .id = "model") %>%
select(term, model, p.value) %>%
spread(model, p.value)
term imputed univariate with_missings
1 (Intercept) 1.214577e-01 1.357530e-34 6.789803e-02
2 age 5.853446e-25 2.133986e-21 2.955927e-23
3 d.motor 1.431263e-35 NA 2.731862e-35
4 d.pupilno reactive pupils NA NA 2.458664e-18
5 d.pupilone reactive NA NA 1.139903e-04
6 pupil.ino reactive pupils 2.407998e-19 NA NA
7 pupil.ione reactive 2.916081e-05 NA NA
How many missing values are imputed in the variable ‘pupil.i’? How many more cases can be analyzed by using ‘pupil.i’ rather than ‘d.pupil’?
See above
The regression coefficients of the logistic model can be used to calculate the individual predicted risk of unfavorable outcome. Fit the model (as in d) again and use the option
. Note that your dataset (not the output screen) shows an extra column.
Predicted probabilities are stored in the R object
fit_imp$fitted.values[1:10]
1 2 3 4 5 6 7
0.1062243 0.1785122 0.1785122 0.1785122 0.1062243 0.2843433 0.1062243
8 9 10
0.1062243 0.1785122 0.1766922
Look at the descriptives of the predicted risks. Is the range very narrow / reasonably wide?
summary(fit_imp$fitted.values)
Min. 1st Qu. Median Mean 3rd Qu. Max.
0.06326 0.20825 0.34340 0.39416 0.54199 0.95926
Looks like it covers a whide range of the 0-1 interval
If the model is well calibrated, groups of patients with low predicted risks will include only few patients with unfavorable outcomes; groups of patients with high predicted risks many. To check this, group the patients by predicted risk (use “recode into different variable”): 1: 0.00 - 0.15 2: 0.15 - 0.30 3: 0.30 - 0.40 4: 0.40 - 0.60 5: 0.60 - 1.00 Give the observed proportions of patients with unfavourable outcome for each group (use crosstabs with option cells, percentage).
Add predicted to data.frame
tbi %<>% mutate(
pred_unfav = fit_imp$fitted.values,
pred_unfav_cat = cut(pred_unfav, breaks = c(0, .15, .3, .4, .6, 1)))
View results
tbi %>%
group_by(pred_unfav_cat) %>%
summarize(observed_prob = mean(d.unfav))
# A tibble: 5 x 2
pred_unfav_cat observed_prob
<fct> <dbl>
1 (0,0.15] 0.135
2 (0.15,0.3] 0.227
3 (0.3,0.4] 0.314
4 (0.4,0.6] 0.487
5 (0.6,1] 0.766
These line up OK
Let’s plot calibration
tbi %>%
mutate(pred_unfav_deciles = quant(pred_unfav, n.tiles = 10)) %>%
group_by(pred_unfav_deciles) %>%
summarize(observed_prob = mean(d.unfav),
observed_prob_lo = binom.confint_logical(d.unfav)$lower,
observed_prob_hi = binom.confint_logical(d.unfav)$upper) %>%
ggplot(aes(x = seq(0.05, .95, length.out = 10),
ymin = observed_prob_lo, ymax = observed_prob_hi,
y = observed_prob)) +
geom_errorbar() + geom_point() +
geom_abline(aes(slope = 1, intercept = 0), lty = 2) +
lims(y = c(0,1))
Use the Hosmer-Lemeshow test to assess the calibration of the model as fitted in step d). By default, this test groups patients by deciles of risk. Does the test give a statistically significant result? Is that to be expected when a model is fitted and tested for fit in the same data?
ResourceSelection::hoslem.test(x = tbi$d.unfav, y = tbi$pred_unfav, g = 10)
Hosmer and Lemeshow goodness of fit (GOF) test
data: tbi$d.unfav, tbi$pred_unfav
X-squared = 5.1937, df = 8, p-value = 0.7367
No rejection of null-hypothesis. Seems to fit OK.
Howerever the fit was based on the data, so this may be overfitted.
What do you think would happen with calibration at external validation, i.e. predictions are made for another data set?
Probably, calibration will be worse.
However, we have 851 cases, and fitted 5 degrees of freedom, so overfitting should be limited
Study the discriminative ability of the model with the ROC curve. Use
. For comparison, make also a ROC curve for age alone as a single predictor.
logit_roc <- function(fit, add = F, ...) {
if (!("glm" %in% class(fit)) & fit$family$family == "binomial" & fit$family$link == "logit") {
stop("only works for glm fits with family = binomial(link = 'logit')")
}
formula0 = formula(fit)
all_vars = all.vars(formula0)
response = all_vars[1]
all_terms = all_vars[-1]
roc <- pROC::roc(fit$data[[response]], fit$fitted.values, ci = T)
pROC:::plot.roc.roc(roc, ci = T, add = add, ...)
}
logit_roc(fit_imp, main = "multivariate model")
logit_roc(fit_imp, main = "comparison with univariate model of only age")
logit_roc(fits_uni[["age"]], add = T)
What would you expect for the area under the ROC-curve) if the model were applied in a new data set (external validation)?
A little worse.
What would you expect if the prognostic model is developed in a selection of patients with very narrow inclusion criteria with respect to important predictors such as age and motor score?
It will do a worse, since contrasts are smaller.
The aim of this practical is to introduce you to some of the more commonly used techniques to perform shrinkage and validation. In this practical you will use R (SPSS doesn’t have sufficient functionality).
Start R This practical can be done in R-studio. Once R-studio is started make sure you open a new R-script (File -> New File -> R script). Start your script by including the following code to load the appropriate R-packages (copy, paste and run (Control+R)):
require(glmnet)
require(rms)
require(ggplot2)
require(logistf)
If these packages haven’t been installed already, you may install them by using: install.packages(“package name”) Data The easiest way to work with data in R is to first assign a so-called working directory. For instance, if you want to use the folder: “H/prognostic_research/practical3”, you can use the following code:
setwd(“H:/prognostic_research/practical3”) We are going to use a slightly modified version of the Traumatic Brain Injury (TBI) data (TBIday3.RDS). If you haven’t done so already, download this data set (epidemiology-education.nl) and store it in your working directory. Then use these code lines to load these data in R and split it up into a development and a validation data set:
tbi2 <- readRDS(here("data", "TBIday3.RDS"))
devdata <- tbi2[tbi2$trial=="Tirilazad US",-1]
valdata <- tbi2[tbi2$trial=="Tirilazad International",-1]
We are going to model mortality at 6 months (d.mort) as a function of ten prognostic factors: 1. age 2. hypoxia 3. hypotens 4. glucose 5. hb 6. d.sysbpt 7. edh 8. tsah 9. shift 10. pupil.dich (dichotomized to “both pupils reactive” (0) and “not both pupils reactive”" (1)). The definitions of these prognostic factor are found in the documents of the practical of day 2.
Q1: What are the numbers of events in both data sets? Are the numbers of events sufficient to warrant development and validation? Code:
table(devdata$d.mort)
0 1
816 225
table(valdata$d.mort)
0 1
840 278
Answer: The number of events in the development data is 225 and in the validation data is 278. There seems enough data for a validation study (>200 events). Given that we have 10 candidate predictors (and we assume no interaction or non-linear relationships with the outcome), EPV = 22.5 for development, which is higher than the often suggested minimum (EPV=10) but lower than EPV=50. One may therefore expect that these models may suffer from some level of overfitting, especially when variable selection strategies are employed.
Study the the distribution of the continuous predictor variables in the development data: age, glucose, hb and d.sysbpt. Are these variable normally distributed? If not: should we make adjustments? Code:
# ggplot(aes(age, colour = as.factor(d.mort)),data=devdata)+geom_density()
# ggplot(aes(glucose, colour = as.factor(d.mort)),data=devdata)+geom_density()
# ggplot(aes(hb, colour = as.factor(d.mort)),data=devdata)+geom_density()
# ggplot(aes(d.sysbpt, colour = as.factor(d.mort)),data=devdata)+geom_density()
Or we can do:
pred_vars <- c("age", "hypoxia", "hypotens", "glucose", "hb", "d.sysbpt", "edh",
"tsah", "shift", "pupil.dich")
num_vars <- c("age", "glucose", "hb", "d.sysbpt")
resp_var <- c("d.mort")
tbi2 %>%
select(num_vars, resp_var) %>%
mutate(d.mort = factor(d.mort)) %>%
gather(-d.mort, key = "variable", value = "value") %>%
ggplot(aes(x = value, col = d.mort)) +
geom_density() +
facet_wrap(~variable, scales = "free")
Answer: not all continuous predictor variables are normally distributed. There are, however, no direct assumptions made about the distribution of predictor variables in a logistic regression. No adjustments have to be made at this point.
Look at the univariable associations between the d.mort and binary predictor variables in the development data. Make cross-tables. Code:
table(devdata$hypoxia,devdata$d.mort)
0 1
0 615 123
1 201 102
table(devdata$hypotens,devdata$d.mort)
0 1
0 664 145
1 152 80
table(devdata$edh,devdata$d.mort)
0 1
0 742 210
1 74 15
table(devdata$tsah,devdata$d.mort)
0 1
0 505 86
1 311 139
table(devdata$shift,devdata$d.mort)
0 1
0 698 166
1 118 59
table(devdata$pupil.dich,devdata$d.mort)
0 1
0 602 99
1 214 126
Develop prediction models by maximum likelihood
Develop a prediction model by maximum likelihood with all 10 predictors incorporated (no selection or shrinkage). What is the apparent area under the ROC curve (C statistic) for this model? Code:
# logistic regression model (Full model)
m1 <- d.mort~age+hypoxia+hypotens+glucose+hb+d.sysbpt+edh+tsah+shift+pupil.dich
full_model <- lrm(as.formula(m1),data=devdata,x=T,y=T)
full_model
Logistic Regression Model
lrm(formula = as.formula(m1), data = devdata, x = T, y = T)
Model Likelihood Discrimination Rank Discrim.
Ratio Test Indexes Indexes
Obs 1041 LR chi2 189.85 R2 0.257 C 0.774
0 816 d.f. 10 g 1.260 Dxy 0.548
1 225 Pr(> chi2) <0.0001 gr 3.524 gamma 0.548
max |deriv| 1e-10 gp 0.188 tau-a 0.186
Brier 0.136
Coef S.E. Wald Z Pr(>|Z|)
Intercept -0.2539 0.9437 -0.27 0.7879
age 0.0183 0.0065 2.80 0.0051
hypoxia 0.5319 0.1809 2.94 0.0033
hypotens 0.2444 0.2089 1.17 0.2421
glucose 0.0853 0.0235 3.63 0.0003
hb -0.0579 0.0416 -1.39 0.1638
d.sysbpt -0.0218 0.0058 -3.79 0.0002
edh -0.3965 0.3277 -1.21 0.2264
tsah 0.7960 0.1724 4.62 <0.0001
shift 0.7009 0.2039 3.44 0.0006
pupil.dich 0.9201 0.1713 5.37 <0.0001
Answer: the apparent C statistics is 0.774.
Now use step-wise backward selection with alpha = .05. What is the apparent area under the ROC curve (C statistic) of this model? Which variables are dropped from the model? How does it compare to the full model? Code:
# logistic regression model (backward selection)
selection <- fastbw(full_model,rule="p",sls=.05)
bw_model <- lrm(as.formula(paste("d.mort~",paste(selection$names.kept,collapse="+"))),data=devdata)
bw_model
Logistic Regression Model
lrm(formula = as.formula(paste("d.mort~", paste(selection$names.kept,
collapse = "+"))), data = devdata)
Model Likelihood Discrimination Rank Discrim.
Ratio Test Indexes Indexes
Obs 1041 LR chi2 183.73 R2 0.250 C 0.772
0 816 d.f. 7 g 1.226 Dxy 0.543
1 225 Pr(> chi2) <0.0001 gr 3.409 gamma 0.543
max |deriv| 3e-11 gp 0.185 tau-a 0.184
Brier 0.138
Coef S.E. Wald Z Pr(>|Z|)
Intercept -0.6634 0.7995 -0.83 0.4067
age 0.0201 0.0065 3.11 0.0019
hypoxia 0.5914 0.1774 3.33 0.0009
glucose 0.0946 0.0231 4.10 <0.0001
d.sysbpt -0.0252 0.0055 -4.60 <0.0001
tsah 0.7593 0.1706 4.45 <0.0001
shift 0.6492 0.2021 3.21 0.0013
pupil.dich 0.9359 0.1706 5.49 <0.0001
Answer: hypotens, edh and hb were deleted from the final model. The apparent C statistic of the final model is 0.771, very close to the apparent C statistic of the full model.
Perform internal validation using bootstrap Background Optimism (due to overfitting) can be investigated using the bootstrap. One bootstrap sample is a random sample with replacement of the original data. A bootstrap sample has the same dimensions, i.e. sample size, as the orginal data set. To study optimism: multiple bootstrap samples are generated (say, 1000 bootstrap samples). The prediction model is fitted on each of those samples. If variable selection is applied: this procedure is executed on each bootstrap sample. The predictive performance of these (final) bootstrap models are evaluated on the original data sample. The average of the bootrap performances provides an estimate of performance of the model in the original data sample that is corrected for optimism.
The above described bootstrap procedure is implemented in the validate function (rms library).
Perform an internal validation of the full model (Q4). For computational time reasons: take 200 bootstrap samples.
# internal validation full model
internalfull_model <- validate(full_model,B=200)
internalfull_model
index.orig training test optimism index.corrected n
Dxy 0.5481 0.5628 0.5378 0.0250 0.5231 200
R2 0.2573 0.2709 0.2463 0.0246 0.2327 200
Intercept 0.0000 0.0000 -0.0693 0.0693 -0.0693 200
Slope 1.0000 1.0000 0.9367 0.0633 0.9367 200
Emax 0.0000 0.0000 0.0269 0.0269 0.0269 200
D 0.1814 0.1922 0.1729 0.0193 0.1621 200
U -0.0019 -0.0019 0.0010 -0.0029 0.0010 200
Q 0.1833 0.1941 0.1719 0.0222 0.1611 200
B 0.1365 0.1344 0.1385 -0.0041 0.1406 200
g 1.2596 1.3124 1.2232 0.0892 1.1704 200
gp 0.1876 0.1920 0.1835 0.0085 0.1791 200
Calculate the optimism corrected C statistic for the full model. Make use of the fact that: C = (Dxy/2)+0.5.
internalfull_model[1,] / 2 + 0.5
index.orig training test optimism
0.7740577 0.7814089 0.7689249 0.5124840
index.corrected n
0.7615738 100.5000000
Answer: the bootstrap corrected estimate for the C statistic is about .762. This may vary slightly between executions because bootstrap sampling is a random process.
Perform an internal validation of the backward selection model. For computational time reasons: take 100 bootstrap samples. Hint: the selection must be executed on each bootstrap sample. Code:
# internal validation backward selection model
internvalbw_model <- validate(full_model,bw=T,rule="p",sls=.05,B=200)
Backwards Step-down - Original Model
Deleted Chi-Sq d.f. P Residual d.f. P AIC
hypotens 1.37 1 0.2421 1.37 1 0.2421 -0.63
edh 1.55 1 0.2135 2.92 2 0.2327 -1.08
hb 3.11 1 0.0776 6.03 3 0.1102 0.03
Approximate Estimates after Deleting Factors
Coef S.E. Wald Z P
Intercept -0.64942 0.801306 -0.8104 4.177e-01
age 0.01996 0.006477 3.0816 2.059e-03
hypoxia 0.58772 0.178144 3.2992 9.697e-04
glucose 0.09406 0.023146 4.0638 4.828e-05
d.sysbpt -0.02512 0.005497 -4.5691 4.898e-06
tsah 0.75275 0.171494 4.3893 1.137e-05
shift 0.64423 0.202169 3.1866 1.439e-03
pupil.dich 0.92939 0.171200 5.4287 5.677e-08
Factors in Final Model
[1] age hypoxia glucose d.sysbpt tsah shift
[7] pupil.dich
internvalbw_model
index.orig training test optimism index.corrected n
Dxy 0.5432 0.5553 0.5277 0.0277 0.5155 200
R2 0.2497 0.2633 0.2369 0.0264 0.2233 200
Intercept 0.0000 0.0000 -0.0789 0.0789 -0.0789 200
Slope 1.0000 1.0000 0.9328 0.0672 0.9328 200
Emax 0.0000 0.0000 0.0297 0.0297 0.0297 200
D 0.1755 0.1868 0.1657 0.0211 0.1544 200
U -0.0019 -0.0019 0.0007 -0.0027 0.0007 200
Q 0.1775 0.1887 0.1650 0.0238 0.1537 200
B 0.1380 0.1364 0.1398 -0.0034 0.1414 200
g 1.2265 1.2781 1.1839 0.0942 1.1323 200
gp 0.1849 0.1905 0.1798 0.0107 0.1742 200
Factors Retained in Backwards Elimination
age hypoxia hypotens glucose hb d.sysbpt edh tsah shift pupil.dich
* * * * * * *
* * * * * * *
* * * * * *
* * * * *
* * * * * *
* * * * * *
* * * * * *
* * * * * * * *
* * * * * * *
* * * * * * *
* * * * * *
* * * * * * *
* * * * * * * *
* * * * * * *
* * * * *
* * * * * *
* * * * * * * *
* * * * * *
* * * * * * *
* * * * * * *
* * * * * *
* * * * * * *
* * * * * *
* * * * * * *
* * * * * * * *
* * * * * * *
* * * * * * *
* * * * * *
* * * * * * * *
* * * * * *
* * * * *
* * * * * *
* * * * * * *
* * * * * *
* * * * * * * *
* * * * * * * *
* * * * * *
* * * * * * *
* * * * * * *
* * * * * * * *
* * * * * * *
* * * * * * * *
* * * * * *
* * * * * * * * *
* * * * * * *
* * * * * * * *
* * * * * * * *
* * * * * * * * *
* * * * * * * *
* * * * * * *
* * * * * * * *
* * * * * * *
* * * * *
* * * * * * *
* * * * * * *
* * * * * * *
* * * * * * * *
* * * * * * * * *
* * * * *
* * * * * *
* * * * * * *
* * * * * * * *
* * * * * *
* * * * * *
* * * * * * *
* * * * * *
* * * * *
* * * * * * *
* * * * * *
* * * * * * *
* * * * * *
* * * * * * * *
* * * * * *
* * * * * * *
* * * * * *
* * * * * *
* * * * * * *
* * * * * * *
* * * * * * *
* * * * * *
* * * * * * *
* * * * * *
* * * * * * *
* * * * * * *
* * * * * * *
* * * * * * * *
* * * * * * * *
* * * * * * *
* * * * * * * *
* * * * * * * *
* * * * * * *
* * * * * * *
* * * * * *
* * * * * * * *
* * * * * *
* * * * * * *
* * * * * * * * *
* * * * * * *
* * * * * *
* * * * *
* * * * * * *
* * * * * * * *
* * * * * * *
* * * * * * * *
* * * * * * *
* * * * * *
* * * * * * *
* * * * * * * *
* * * * * * *
* * * * * * *
* * * * * * *
* * * * * * *
* * * * * * * *
* * * * * * *
* * * * * *
* * * * * *
* * * * * * *
* * * * * * * *
* * * * * *
* * * * * * * *
* * * * * * *
* * * * * * *
* * * * * * * *
* * * * * * * *
* * * * * * *
* * * * * * *
* * * * * * * *
* * * * * * *
* * * * * * *
* * * * * * * *
* * * * *
* * * * * * *
* * * * * * * *
* * * * * * *
* * * * *
* * * * * * *
* * * * * * *
* * * * * * * *
* * * * * * *
* * * * * *
* * * * * * * *
* * * * * * * *
* * * * * * *
* * * * * *
* * * * * * *
* * * * * * *
* * * * * * *
* * * * * * *
* * * * * * * *
* * * * *
* * * * * * *
* * * * * * * *
* * * * * * *
* * * * * * * *
* * * * * * *
* * * * * * *
* * * * * *
* * * * * * *
* * * * * *
* * * * * * *
* * * * * *
* * * * * * * *
* * * * * * * *
* * * * * * * *
* * * * * * * *
* * * * * * * * *
* * * * * * * *
* * * * * * * *
* * * * * *
* * * * * * *
* * * * *
* * * * * * *
* * * * * * * *
* * * * * * * *
* * * * * * *
* * * * * * *
* * * * * *
* * * * * * *
* * * * * * * *
* * * * * * *
* * * * * * *
* * * * * * * * *
* * * * * * *
* * * * * * *
* * * * * * * * *
* * * * * * * *
* * * * * * *
* * * * * * *
* * * * * * *
* * * * * * * *
* * * * * * *
* * * * * * *
* * * * * * *
* * * * * * *
* * * * * * * *
* * * * * * *
* * * * * *
* * * * * * * * *
* * * * *
* * * * * * *
Frequencies of Numbers of Factors Retained
5 6 7 8 9
12 41 89 50 8
Calculate the optimism corrected C statistic for the backward selection model.
internvalbw_model[1, ] / 2 + .5
index.orig training test optimism
0.7715795 0.7776680 0.7638259 0.5138421
index.corrected n
0.7577374 100.5000000
Answer: the bootstrap corrected estimate for the C statistic is about .756.
Does this internal validation exercise provide evidence of over- or underfitting? If so, which of these models are affected? Answer: the corrected calibration slopes (“Slope” in the output of the validation function) are around .942 (full model) and .928 (after backward selection). This indicates that the models suffer from some overfitting (Slope = 1 indicates no overfitting, Slope > 1 indicates underfitting).
maximum likelihood with or without stepwise selection is still the most commonly used approach for developing prediction models. However, it is long known that maximum likelihood estimation is not optimal for prediction purposes. For this, the maximum likelihood estimates need shrinkage.
Maximum likelihood estimation of the logistic model proceeds by maximizing the log-likelihood function: logL=∑iyilogπi+(1−yi)log(1−πi), where i stands for individual i, yi is the observed outcome for individual i and πi is the predicted outcome for individual i. There are several methods available to perform shrinkage. We here discuss three methods that are also known as penalized likelihood models. Each of these methods have the form: logL−p(⋅), where p(⋅) stands for a penalty function. Firth’s correction gives penalty: −1/2log|I(θ)|, where log|I(θ)| denotes the Fisher information matrix; Ridge gives penalty: p(⋅)=λ∑j=1β2j, where λ denotes a so-called tuning parameter and βj denotes regression coefficient j; and Lasso gives penalty p(⋅)=λ∑dfj=1|βj|. Estimating the tuning parameters (λ) for Lasso and Ridge is often done using a cross-validation approach. Details appear in a book called “The Elements of Statistical Learning” by Trevor Hastie et al.
Develop a model using Firth’s correction (all ten variables included). Compare the estimated regression coefficients to the full model (Q4). Code:
require(logistf)
# logistic regression model with Firth's correction (Full model)
firth_model <- logistf(as.formula(m1),data=devdata,firth=T)
firth_model
logistf(formula = as.formula(m1), data = devdata, firth = T)
Model fitted by Penalized ML
Confidence intervals and p-values by Profile Likelihood
coef se(coef) lower 0.95 upper 0.95 Chisq
(Intercept) -0.27013768 0.939190337 -2.101399116 1.57069188 0.08331424
age 0.01814908 0.006520080 0.005389062 0.03088027 7.74391852
hypoxia 0.52743276 0.180101206 0.173874221 0.87796228 8.49154615
hypotens 0.24459416 0.208049998 -0.165820246 0.64744970 1.37668381
glucose 0.08375997 0.023409911 0.038669281 0.13009905 13.35964125
hb -0.05700403 0.041403896 -0.138035973 0.02390951 1.90824278
d.sysbpt -0.02138948 0.005724866 -0.032728457 -0.01033559 14.60276525
edh -0.37102361 0.324303059 -1.035117305 0.23610654 1.39187602
tsah 0.78491191 0.171403312 0.451409218 1.12256368 21.42151200
shift 0.69423880 0.203061630 0.294512223 1.08807414 11.39642590
pupil.dich 0.90856021 0.170509291 0.575653707 1.24227119 28.38276217
p
(Intercept) 7.728553e-01
age 5.389373e-03
hypoxia 3.568005e-03
hypotens 2.406668e-01
glucose 2.570975e-04
hb 1.671586e-01
d.sysbpt 1.327197e-04
edh 2.380885e-01
tsah 3.686122e-06
shift 7.358554e-04
pupil.dich 9.954776e-08
Likelihood ratio test=188.3889 on 10 df, p=0, n=1041
Answer: the coefficients are slightly ‘shrunken’ towards to zero effect as compared to the orignal full model estimated with maximum likelihood.
cbind(max_likelihood = coef(full_model), firth = coef(firth_model))
max_likelihood firth
Intercept -0.25394896 -0.27013768
age 0.01834929 0.01814908
hypoxia 0.53189103 0.52743276
hypotens 0.24443357 0.24459416
glucose 0.08526859 0.08375997
hb -0.05791975 -0.05700403
d.sysbpt -0.02179392 -0.02138948
edh -0.39645485 -0.37102361
tsah 0.79602655 0.78491191
shift 0.70091175 0.69423880
pupil.dich 0.92005534 0.90856021
Develop a model using Ridge (all ten variables included). Compare the estimated regression coefficients to the full model (Q4). Note: estimating this model may take some time. Code:
# logistic Ridge regression model using leave-one-out cross-validation
ridge_tuning_parameter <-cv.glmnet(
x=as.matrix(devdata[,-1]),
y=as.matrix(devdata[,1]),
family="binomial",type.measure="mse",
alpha=0,nfolds=nrow(devdata))$lambda.min
Warning: Option grouped=FALSE enforced in cv.glmnet, since < 3 observations
per fold
ridge_model <-glmnet(
x=as.matrix(devdata[,-1]),
y=as.matrix(devdata[,1]),
family="binomial",
lambda=ridge_tuning_parameter,alpha=0)
Answer: the coeffients can be seen using coef(ridge_model). When compared to the full model estimated with maximum likelihood, these coefficients are shrunken.
cbind(max_likelihood = coef(full_model),
firth = coef(firth_model),
ridge = coef(ridge_model))
11 x 3 sparse Matrix of class "dgCMatrix"
max_likelihood firth s0
(Intercept) -0.25394896 -0.27013768 -0.40058227
age 0.01834929 0.01814908 0.01616024
hypoxia 0.53189103 0.52743276 0.48558752
hypotens 0.24443357 0.24459416 0.25767463
glucose 0.08526859 0.08375997 0.07768179
hb -0.05791975 -0.05700403 -0.05604972
d.sysbpt -0.02179392 -0.02138948 -0.01870553
edh -0.39645485 -0.37102361 -0.31366884
tsah 0.79602655 0.78491191 0.69161047
shift 0.70091175 0.69423880 0.61588306
pupil.dich 0.92005534 0.90856021 0.82148405
Develop a model using Lasso (all ten variables included). Note: estimating this model may take some time. Code:
# logistic Ridge regression model using leave-one-out cross-validation
lasso_tuning_parameter <- cv.glmnet(
x=as.matrix(devdata[,-1]),
y=as.matrix(devdata[,1]),
family="binomial",type.measure="mse",
alpha=1,nfolds=nrow(devdata))$lambda.min
Warning: Option grouped=FALSE enforced in cv.glmnet, since < 3 observations
per fold
lasso_model <-glmnet(
x=as.matrix(devdata[,-1]),
y=as.matrix(devdata[,1]),
family="binomial",
lambda=lasso_tuning_parameter,
alpha=1)
# fits <- list(
# max_likelihood = (full_model),
# firth = (firth_model),
# ridge = (ridge_model),
# lasso = (lasso_model))
coef(lasso_model)
11 x 1 sparse Matrix of class "dgCMatrix"
s0
(Intercept) -0.28116274
age 0.01786273
hypoxia 0.52155721
hypotens 0.23911329
glucose 0.08413384
hb -0.05593184
d.sysbpt -0.02140241
edh -0.36101560
tsah 0.78016590
shift 0.68173252
pupil.dich 0.90901921
Answer: the coefficients can be seen using coef(lasso_model). Notice that, unlike ridge, Lasso may perform automated selection by shrinking some variables to zero.
study possible case-mix differences. Are there differences between the development and validation data? Example code
# compare case-mix between development and validation data
by(tbi2,tbi2$trial,function(x)colMeans(x[,-1]))
tbi2$trial: Tirilazad International
d.mort age hypoxia hypotens glucose hb
0.2486583 33.6127013 0.1556351 0.1440072 8.6079944 11.8764758
d.sysbpt edh tsah shift pupil.dich
130.0886673 0.1681574 0.5277281 0.2262970 0.2754919
--------------------------------------------------------
tbi2$trial: Tirilazad US
d.mort age hypoxia hypotens glucose
0.21613833 32.78001921 0.29106628 0.22286263 9.22702530
hb d.sysbpt edh tsah shift
13.05240154 132.86975024 0.08549472 0.43227666 0.17002882
pupil.dich
0.32660903
by(tbi2,tbi2$trial,function(x)table(x$hypoxia))
tbi2$trial: Tirilazad International
0 1
944 174
--------------------------------------------------------
tbi2$trial: Tirilazad US
0 1
738 303
ggplot(tbi2,aes(hb,colour=trial))+geom_density()
Answer: yes, there are differences in case-mix. For instance, hypoxia is more common in the development data than in the validation data.
Or, using data.table
setDT(tbi2)
tbi2[, id := .I]
tbi2 %>%
data.table::melt(id.vars = c("trial", "id")) %>%
.[, list(mean = mean(value), sd = sd(value)), by = c("trial", "variable")] %>%
data.table::melt(id.vars = c("trial", "variable"),
variable.name = "measure", value.name = "value") %>%
data.table::dcast(variable~trial+measure)
variable Tirilazad International_mean Tirilazad International_sd
1: d.mort 0.2486583 0.4324287
2: age 33.6127013 14.5212252
3: hypoxia 0.1556351 0.3626713
4: hypotens 0.1440072 0.3512541
5: glucose 8.6079944 3.5763887
6: hb 11.8764758 2.7624993
7: d.sysbpt 130.0886673 16.0103923
8: edh 0.1681574 0.3741734
9: tsah 0.5277281 0.4994540
10: shift 0.2262970 0.4186208
11: pupil.dich 0.2754919 0.4469618
Tirilazad US_mean Tirilazad US_sd
1: 0.21613833 0.4118075
2: 32.78001921 12.4776194
3: 0.29106628 0.4544723
4: 0.22286263 0.4163669
5: 9.22702530 3.4573085
6: 13.05240154 2.1583300
7: 132.86975024 15.5752142
8: 0.08549472 0.2797509
9: 0.43227666 0.4956304
10: 0.17002882 0.3758387
11: 0.32660903 0.4691983
evaluate the predictive performance at external validiation Code
# external performance full model
predfull_model <- predict(full_model, newdata=valdata, type = "fitted")
val.prob(predfull_model,valdata$d.mort)
Dxy C (ROC) R2 D D:Chi-sq
5.613866e-01 7.806933e-01 2.626398e-01 1.940168e-01 2.179108e+02
D:p U U:Chi-sq U:p Q
NA 3.776862e-03 6.222531e+00 4.454454e-02 1.902400e-01
Brier Intercept Slope Emax E90
1.498881e-01 2.421173e-01 1.048814e+00 6.937963e-02 6.591928e-02
Eavg S:z S:p
2.628233e-02 1.519696e+00 1.285874e-01
# external performance backward selection model
predbw_model <- 1/(1+exp(-predict(bw_model, newdata=valdata)))
val.prob(predbw_model,valdata$d.mort)
Dxy C (ROC) R2 D D:Chi-sq
0.53039997 0.76519998 0.23259825 0.16969935 190.72387412
D:p U U:Chi-sq U:p Q
NA 0.00486532 7.43942750 0.02424091 0.16483403
Brier Intercept Slope Emax E90
0.15412083 0.20965170 0.99847298 0.06774686 0.06643357
Eavg S:z S:p
0.02993510 2.18963239 0.02855091
# external validation Firth's model
X <- model.matrix(firth_model, data=valdata)
predfirth_model <- 1/(1+exp(-X %*% coef(firth_model)))
val.prob(predfirth_model,valdata$d.mort)
Dxy C (ROC) R2 D D:Chi-sq
5.610997e-01 7.805498e-01 2.631810e-01 1.944604e-01 2.184067e+02
D:p U U:Chi-sq U:p Q
NA 3.261937e-03 5.646845e+00 5.940228e-02 1.911985e-01
Brier Intercept Slope Emax E90
1.498282e-01 2.434940e-01 1.063815e+00 6.752707e-02 6.474306e-02
Eavg S:z S:p
2.445818e-02 1.265256e+00 2.057797e-01
# external validation Ridge model
predridge_model <- 1/(1+exp(-predict(ridge_model, newx=as.matrix(valdata[,-1]))))
val.prob(predridge_model,valdata$d.mort)
Dxy C (ROC) R2 D D:Chi-sq
5.621489e-01 7.810744e-01 2.585266e-01 1.906522e-01 2.141492e+02
D:p U U:Chi-sq U:p Q
NA 7.314438e-03 1.017754e+01 6.165594e-03 1.833378e-01
Brier Intercept Slope Emax E90
1.507464e-01 3.874152e-01 1.176764e+00 8.917287e-02 8.295792e-02
Eavg S:z S:p
2.881835e-02 8.649746e-01 3.870527e-01
# external validation Lasso model
predlasso_model <- 1/(1+exp(-predict(lasso_model, newx=as.matrix(valdata[,-1]))))
val.prob(predlasso_model,valdata$d.mort)
Dxy C (ROC) R2 D D:Chi-sq
5.608342e-01 7.804171e-01 2.615932e-01 1.931597e-01 2.169526e+02
D:p U U:Chi-sq U:p Q
NA 4.103177e-03 6.587352e+00 3.711715e-02 1.890565e-01
Brier Intercept Slope Emax E90
1.500865e-01 2.637481e-01 1.068436e+00 7.152309e-02 6.881476e-02
Eavg S:z S:p
2.645759e-02 1.397576e+00 1.622404e-01
In this practical we are going to develop a prognostic model based on time-to-event data. The dataset is an abstraction of the SMART dataset, and contains the following variables: - TEVENT describing follow-up duration in days - EVENT occurrence of new vascular complications - SEX - AGE age of the patient at baseline, in years - CARDIAC previous cardiac vascular problems - BMI - HDL - CREAT creatinine value - IMT intima media thickness We use the dataset SMARTdev.sav for model development, and SMARTval.sav for validation of the model.
Prepare the dataset before modelling. Check the distribution of the timing of events and the timing of censoring. Do you expect any problems in the modeling and in assessment of performance? Check the distribution of other variables. Create standardized versions of the continuous variables in the model (Age, BMI, HDL, CREAT, IMT) by normalization.
Load data
require(haven)
dev <- haven::read_sav(here("data", "SMARTdev.sav"))
val <- haven::read_sav(here("data", "SMARTval.sav"))
For curation, we can conbine the datasets
smart <- rbindlist(list(dev = dev, val = val), idcol = "set")
str(smart)
Classes 'data.table' and 'data.frame': 3735 obs. of 10 variables:
$ set : chr "dev" "dev" "dev" "dev" ...
$ TEVENT : atomic 3466 3465 3465 2445 3463 ...
..- attr(*, "label")= chr "FU-duur tot vasculaire complicatie (klinisch) EP (in dagen)"
..- attr(*, "format.spss")= chr "F5.0"
$ EVENT :Class 'labelled' atomic [1:3735] 0 0 0 1 0 0 1 1 1 1 ...
.. ..- attr(*, "label")= chr "Niet/wel vasculaire complicatie (klinisch) EP"
.. ..- attr(*, "format.spss")= chr "F1.0"
.. ..- attr(*, "labels")= Named num [1:4] 0 1 2 3
.. .. ..- attr(*, "names")= chr [1:4] "FU afgekapt" "Vasculaire complicatie (klinisch) EP" "Lost to FU" "Overleden (anderszins)"
$ SEX : atomic 1 2 1 1 1 2 1 1 1 1 ...
..- attr(*, "label")= chr "Geslacht"
..- attr(*, "format.spss")= chr "F1.0"
$ AGE : atomic 71 46 59 76 57 52 66 72 75 53 ...
..- attr(*, "label")= chr "Leeftijd (aantal voltooide jaren)"
..- attr(*, "format.spss")= chr "F3.0"
$ CARDIAC: atomic 1 0 1 1 1 0 0 1 1 0 ...
..- attr(*, "label")= chr "Ooit cardiaal vaatlijden (inclusief voorgeschiedenis)"
..- attr(*, "format.spss")= chr "F1.0"
$ BMI : atomic 23.9 24 26.2 29.6 29.6 ...
..- attr(*, "format.spss")= chr "F4.1"
$ HDL : atomic 0.94 1.26 1.28 1 0.81 0.95 0.87 1.32 0.58 0.95 ...
..- attr(*, "format.spss")= chr "F4.2"
$ CREAT : atomic 95 66 93 79 91 180 81 101 157 86 ...
..- attr(*, "format.spss")= chr "F6.2"
$ IMT : atomic 0.82 0.57 0.83 1.45 1.07 0.97 1.15 0.7 1.48 0.87 ...
..- attr(*, "format.spss")= chr "F4.2"
- attr(*, ".internal.selfref")=<externalptr>
Convert ‘labelled’ to factors
smart %<>% mutate_if(is.labelled, as_factor)
str(smart)
'data.frame': 3735 obs. of 10 variables:
$ set : chr "dev" "dev" "dev" "dev" ...
$ TEVENT : atomic 3466 3465 3465 2445 3463 ...
..- attr(*, "label")= chr "FU-duur tot vasculaire complicatie (klinisch) EP (in dagen)"
..- attr(*, "format.spss")= chr "F5.0"
$ EVENT : Factor w/ 4 levels "FU afgekapt",..: 1 1 1 2 1 1 2 2 2 2 ...
..- attr(*, "label")= chr "Niet/wel vasculaire complicatie (klinisch) EP"
$ SEX : atomic 1 2 1 1 1 2 1 1 1 1 ...
..- attr(*, "label")= chr "Geslacht"
..- attr(*, "format.spss")= chr "F1.0"
$ AGE : atomic 71 46 59 76 57 52 66 72 75 53 ...
..- attr(*, "label")= chr "Leeftijd (aantal voltooide jaren)"
..- attr(*, "format.spss")= chr "F3.0"
$ CARDIAC: atomic 1 0 1 1 1 0 0 1 1 0 ...
..- attr(*, "label")= chr "Ooit cardiaal vaatlijden (inclusief voorgeschiedenis)"
..- attr(*, "format.spss")= chr "F1.0"
$ BMI : atomic 23.9 24 26.2 29.6 29.6 ...
..- attr(*, "format.spss")= chr "F4.1"
$ HDL : atomic 0.94 1.26 1.28 1 0.81 0.95 0.87 1.32 0.58 0.95 ...
..- attr(*, "format.spss")= chr "F4.2"
$ CREAT : atomic 95 66 93 79 91 180 81 101 157 86 ...
..- attr(*, "format.spss")= chr "F6.2"
$ IMT : atomic 0.82 0.57 0.83 1.45 1.07 0.97 1.15 0.7 1.48 0.87 ...
..- attr(*, "format.spss")= chr "F4.2"
Let’s look at the event variable
tabl(smart$EVENT)
FU afgekapt Vasculaire complicatie (klinisch) EP
3292 443
Lost to FU Overleden (anderszins)
0 0
<NA>
0
We only have 2 levels for the event variable. Lets recode this to a logical
smart %<>%
mutate(vasc = EVENT == "Vasculaire complicatie (klinisch) EP") %>%
as.data.table()
tabl(smart$EVENT, smart$vasc)
FALSE TRUE <NA>
FU afgekapt 3292 0 0
Vasculaire complicatie (klinisch) EP 0 443 0
Lost to FU 0 0 0
Overleden (anderszins) 0 0 0
<NA> 0 0 0
Survival of uncensored observations
require(ggfortify)
smart %>%
filter(vasc == T) %>%
survfit(Surv(TEVENT, vasc)~1, data = .) %>%
autoplot(ylim = c(0,1))
Censoring
require(ggfortify)
smart %>%
filter(vasc == F) %>%
survfit(Surv(TEVENT, !vasc)~1, data = .) %>%
autoplot(ylim = c(0,1))
Both in a picture
smart %>%
mutate(dummy_event = T,
status = factor(vasc, levels = c(T, F), labels = c("event", "censored"))) %>%
survfit(Surv(TEVENT, dummy_event)~status, data = .) %>%
autoplot()
Looks like most of the events occur before censoring, so should be OK
Check distributions of covariates
all_covs <- c("SEX", "AGE", "CARDIAC", "BMI", "HDL", "CREAT", "IMT")
num_covs <- c("AGE", "BMI", "HDL", "CREAT", "IMT")
cat_covs <- setdiff(all_covs, num_covs)
allvars <- c(all_covs, "vasc")
smart %>%
melt.data.table(id.vars = c("set", "EVENT"), measure.vars = num_covs) %>%
ggplot(aes(x = value, fill = set)) +
geom_histogram(alpha = 0.8) +
facet_wrap(~variable, scales = "free")
Most variables look approximately normally distributed, except for creat, which is right-skewed. We can take the log.
smart %<>%
mutate(log_creat = log(CREAT)) %>%
as.data.table()
smart %>%
ggplot(aes(x = log_creat, fill = set)) +
geom_histogram(alpha = 0.8)
num_covs <- c(num_covs, "log_creat")
std_covs <- paste0(num_covs, "_std")
smart[, (std_covs):=map(.SD, scale), .SDcols = num_covs]
smart[1:10, .SD, .SDcols = c(num_covs, std_covs)]
AGE BMI HDL CREAT IMT log_creat AGE_std BMI_std HDL_std
1: 71 23.88 0.94 95 0.82 4.553877 1.09365050 -0.7306108 -0.78628233
2: 46 24.00 1.26 66 0.57 4.189655 -1.29076193 -0.6992909 0.08384594
3: 59 26.17 1.28 93 0.83 4.532599 -0.05086747 -0.1329235 0.13822896
4: 76 29.59 1.00 79 1.45 4.369448 1.57053298 0.7596924 -0.62313328
5: 57 29.64 0.81 91 1.07 4.510860 -0.24162046 0.7727423 -1.13977194
6: 52 34.41 0.95 180 0.97 5.192957 -0.71850294 2.0177066 -0.75909082
7: 66 24.24 0.87 81 1.15 4.394449 0.61676801 -0.6366512 -0.97662288
8: 72 31.25 1.32 101 0.70 4.615121 1.18902700 1.1929504 0.24699500
9: 75 27.68 0.58 157 1.48 5.056246 1.47515649 0.2611847 -1.76517663
10: 53 22.40 0.95 86 0.87 4.454347 -0.62312645 -1.1168890 -0.75909082
CREAT_std IMT_std log_creat_std
1: -0.05264070 -0.4294153 0.11213674
2: -0.49909020 -1.3655967 -1.09632628
3: -0.08343032 -0.3919680 0.04153985
4: -0.29895766 1.9297620 -0.49978561
5: -0.11421994 0.5067662 -0.03059187
6: 1.25591818 0.1322936 2.23255829
7: -0.26816804 0.8063442 -0.41683308
8: 0.03972817 -0.8787824 0.31533870
9: 0.90183754 2.0421038 1.77896083
10: -0.19119399 -0.2421790 -0.21809533
Tables of categoric variables
cat_covs
[1] "SEX" "CARDIAC"
map(cat_covs, function(var) tabl(smart[[var]], smart$set))
[[1]]
dev val <NA>
1 1989 815 0
2 630 301 0
<NA> 0 0 0
[[2]]
dev val <NA>
0 1149 485 0
1 1470 631 0
<NA> 0 0 0
For now we will skip any variable selection, and include all potential predictors in the model. Create a Cox model using the variables in their original scale, and apply transformations (splines, log, etc) to continuous variables if necessary.
Go back to development and validation set
dev <- smart %>% filter(set == "dev")
val <- smart %>% filter(set == "val")
Let’s assess any non-linear efffets of the continous predictors. We will do this for all variables at the same time.
Also possible would be to do this for each predictor separately.
rcspline.plot(x = dev$AGE_std, y = dev$TEVENT, event = dev$vasc,
model = "cox", adj = dev %>%
select(BMI_std, HDL_std, log_creat_std, IMT_std,
SEX, CARDIAC))
Warning in coxph.fit(cbind(x, adj), cbind(y, event), strata = NULL, offset
= NULL, : Loglik converged before variable 6 ; beta may be infinite.
x BMI_std
-0.12315128 1.45776093 -7.56504433 13.92811919 -0.04051853
HDL_std log_creat_std IMT_std SEX CARDIAC
-0.13550992 0.28844034 0.27571495 -0.05712248 0.12548775
[1] -2262.595 -2179.042
cbind(xe, lower, upper)
xe
[1,] -3.102915369 -1.663278701 2.4275347
[2,] -3.094635607 -1.658840439 2.4210571
[3,] -3.086355844 -1.654402176 2.4145795
[4,] -3.078076081 -1.649963914 2.4081019
[5,] -3.069796318 -1.645525651 2.4016243
[6,] -3.061516556 -1.641087389 2.3951468
[7,] -3.053236793 -1.636649126 2.3886692
[8,] -3.044957030 -1.632210864 2.3821916
[9,] -3.036677268 -1.627772601 2.3757140
[10,] -3.028397505 -1.623334339 2.3692364
[11,] -3.020117742 -1.618896076 2.3627588
[12,] -3.011837980 -1.614457814 2.3562812
[13,] -3.003558217 -1.610019551 2.3498036
[14,] -2.995278454 -1.605581289 2.3433260
[15,] -2.986998692 -1.601143026 2.3368485
[16,] -2.978718929 -1.596704764 2.3303709
[17,] -2.970439166 -1.592266501 2.3238933
[18,] -2.962159404 -1.587828239 2.3174157
[19,] -2.953879641 -1.583389976 2.3109381
[20,] -2.945599878 -1.578951714 2.3044605
[21,] -2.937320116 -1.574513451 2.2979829
[22,] -2.929040353 -1.570075189 2.2915053
[23,] -2.920760590 -1.565636926 2.2850277
[24,] -2.912480828 -1.561198664 2.2785501
[25,] -2.904201065 -1.556760401 2.2720726
[26,] -2.895921302 -1.552322139 2.2655950
[27,] -2.887641540 -1.547883876 2.2591174
[28,] -2.879361777 -1.543445614 2.2526398
[29,] -2.871082014 -1.539007351 2.2461622
[30,] -2.862802252 -1.534569089 2.2396846
[31,] -2.854522489 -1.530130826 2.2332070
[32,] -2.846242726 -1.525692564 2.2267294
[33,] -2.837962963 -1.521254301 2.2202518
[34,] -2.829683201 -1.516816039 2.2137743
[35,] -2.821403438 -1.512377776 2.2072967
[36,] -2.813123675 -1.507939514 2.2008191
[37,] -2.804843913 -1.503501251 2.1943415
[38,] -2.796564150 -1.499062989 2.1878639
[39,] -2.788284387 -1.494624726 2.1813863
[40,] -2.780004625 -1.490186464 2.1749087
[41,] -2.771724862 -1.485748201 2.1684311
[42,] -2.763445099 -1.481309939 2.1619535
[43,] -2.755165337 -1.476871676 2.1554760
[44,] -2.746885574 -1.472433414 2.1489984
[45,] -2.738605811 -1.467995151 2.1425208
[46,] -2.730326049 -1.463556889 2.1360432
[47,] -2.722046286 -1.459118626 2.1295656
[48,] -2.713766523 -1.454680364 2.1230880
[49,] -2.705486761 -1.450242101 2.1166104
[50,] -2.697206998 -1.445803839 2.1101328
[51,] -2.688927235 -1.441365576 2.1036552
[52,] -2.680647473 -1.436927314 2.0971776
[53,] -2.672367710 -1.432489051 2.0907001
[54,] -2.664087947 -1.428050789 2.0842225
[55,] -2.655808185 -1.423612526 2.0777449
[56,] -2.647528422 -1.419174264 2.0712673
[57,] -2.639248659 -1.414736001 2.0647897
[58,] -2.630968897 -1.410297739 2.0583121
[59,] -2.622689134 -1.405859476 2.0518345
[60,] -2.614409371 -1.401421214 2.0453569
[61,] -2.606129609 -1.396982951 2.0388793
[62,] -2.597849846 -1.392544689 2.0324018
[63,] -2.589570083 -1.388106426 2.0259242
[64,] -2.581290320 -1.383668164 2.0194466
[65,] -2.573010558 -1.379229901 2.0129690
[66,] -2.564730795 -1.374791639 2.0064914
[67,] -2.556451032 -1.370353376 2.0000138
[68,] -2.548171270 -1.365915114 1.9935362
[69,] -2.539891507 -1.361476851 1.9870586
[70,] -2.531611744 -1.357038589 1.9805810
[71,] -2.523331982 -1.352600326 1.9741035
[72,] -2.515052219 -1.348162064 1.9676259
[73,] -2.506772456 -1.343723801 1.9611483
[74,] -2.498492694 -1.339285539 1.9546707
[75,] -2.490212931 -1.334847276 1.9481931
[76,] -2.481933168 -1.330409014 1.9417155
[77,] -2.473653406 -1.325970751 1.9352379
[78,] -2.465373643 -1.321532489 1.9287603
[79,] -2.457093880 -1.317094226 1.9222827
[80,] -2.448814118 -1.312655964 1.9158052
[81,] -2.440534355 -1.308217701 1.9093276
[82,] -2.432254592 -1.303779439 1.9028500
[83,] -2.423974830 -1.299341176 1.8963724
[84,] -2.415695067 -1.294902914 1.8898948
[85,] -2.407415304 -1.290464651 1.8834172
[86,] -2.399135542 -1.286026389 1.8769396
[87,] -2.390855779 -1.281588126 1.8704620
[88,] -2.382576016 -1.277149864 1.8639844
[89,] -2.374296254 -1.272711601 1.8575068
[90,] -2.366016491 -1.268273339 1.8510293
[91,] -2.357736728 -1.263835076 1.8445517
[92,] -2.349456965 -1.259396814 1.8380741
[93,] -2.341177203 -1.254958551 1.8315965
[94,] -2.332897440 -1.250520289 1.8251189
[95,] -2.324617677 -1.246082026 1.8186413
[96,] -2.316337915 -1.241643764 1.8121637
[97,] -2.308058152 -1.237205501 1.8056861
[98,] -2.299778389 -1.232767239 1.7992085
[99,] -2.291498627 -1.228328976 1.7927310
[100,] -2.283218864 -1.223890714 1.7862534
[101,] -2.274939101 -1.219452451 1.7797758
[102,] -2.266659339 -1.215014189 1.7732982
[103,] -2.258379576 -1.210575926 1.7668206
[104,] -2.250099813 -1.206137664 1.7603430
[105,] -2.241820051 -1.201699401 1.7538654
[106,] -2.233540288 -1.197261139 1.7473878
[107,] -2.225260525 -1.192822876 1.7409102
[108,] -2.216980763 -1.188384614 1.7344327
[109,] -2.208701000 -1.183946351 1.7279551
[110,] -2.200421237 -1.179508089 1.7214775
[111,] -2.192141475 -1.175069826 1.7149999
[112,] -2.183861712 -1.170631564 1.7085223
[113,] -2.175581949 -1.166193301 1.7020447
[114,] -2.167302187 -1.161755039 1.6955671
[115,] -2.159022424 -1.157316776 1.6890895
[116,] -2.150742661 -1.152878514 1.6826119
[117,] -2.142462899 -1.148440251 1.6761343
[118,] -2.134183136 -1.144001989 1.6696568
[119,] -2.125903373 -1.139563726 1.6631792
[120,] -2.117623610 -1.135125464 1.6567016
[121,] -2.109343848 -1.130687201 1.6502240
[122,] -2.101064085 -1.126248939 1.6437464
[123,] -2.092784322 -1.121810676 1.6372688
[124,] -2.084504560 -1.117372414 1.6307912
[125,] -2.076224797 -1.112934151 1.6243136
[126,] -2.067945034 -1.108495889 1.6178360
[127,] -2.059665272 -1.104057626 1.6113585
[128,] -2.051385509 -1.099619364 1.6048809
[129,] -2.043105746 -1.095181101 1.5984033
[130,] -2.034825984 -1.090742839 1.5919257
[131,] -2.026546221 -1.086304576 1.5854481
[132,] -2.018266458 -1.081866314 1.5789705
[133,] -2.009986696 -1.077428052 1.5724929
[134,] -2.001706933 -1.072989789 1.5660153
[135,] -1.993427170 -1.068551527 1.5595377
[136,] -1.985147408 -1.064113264 1.5530602
[137,] -1.976867645 -1.059675002 1.5465826
[138,] -1.968587882 -1.055236739 1.5401050
[139,] -1.960308120 -1.050798477 1.5336274
[140,] -1.952028357 -1.046360214 1.5271498
[141,] -1.943748594 -1.041921952 1.5206722
[142,] -1.935468832 -1.037483689 1.5141946
[143,] -1.927189069 -1.033045427 1.5077170
[144,] -1.918909306 -1.028607164 1.5012394
[145,] -1.910629544 -1.024168902 1.4947618
[146,] -1.902349781 -1.019730639 1.4882843
[147,] -1.894070018 -1.015292377 1.4818067
[148,] -1.885790256 -1.010854114 1.4753291
[149,] -1.877510493 -1.006415852 1.4688515
[150,] -1.869230730 -1.001977589 1.4623739
[151,] -1.860950967 -0.997539327 1.4558963
[152,] -1.852671205 -0.993101064 1.4494187
[153,] -1.844391442 -0.988662802 1.4429411
[154,] -1.836111679 -0.984224539 1.4364635
[155,] -1.827831917 -0.979786277 1.4299860
[156,] -1.819552154 -0.975348014 1.4235084
[157,] -1.811272391 -0.970909752 1.4170308
[158,] -1.802992629 -0.966471489 1.4105532
[159,] -1.794712866 -0.962033227 1.4040756
[160,] -1.786433103 -0.957594964 1.3975980
[161,] -1.778153341 -0.953156702 1.3911204
[162,] -1.769873578 -0.948718439 1.3846428
[163,] -1.761593815 -0.944280191 1.3781653
[164,] -1.753314053 -0.939842109 1.3716887
[165,] -1.745034290 -0.935404415 1.3652140
[166,] -1.736754527 -0.930967336 1.3587426
[167,] -1.728474765 -0.926531097 1.3522755
[168,] -1.720195002 -0.922095922 1.3458140
[169,] -1.711915239 -0.917662037 1.3393591
[170,] -1.703635477 -0.913229666 1.3329122
[171,] -1.695355714 -0.908799034 1.3264742
[172,] -1.687075951 -0.904370368 1.3200465
[173,] -1.678796189 -0.899943891 1.3136301
[174,] -1.670516426 -0.895519830 1.3072263
[175,] -1.662236663 -0.891098410 1.3008363
[176,] -1.653956901 -0.886679855 1.2944611
[177,] -1.645677138 -0.882264392 1.2881019
[178,] -1.637397375 -0.877852247 1.2817600
[179,] -1.629117612 -0.873443645 1.2754365
[180,] -1.620837850 -0.869038812 1.2691326
[181,] -1.612558087 -0.864637976 1.2628493
[182,] -1.604278324 -0.860241362 1.2565880
[183,] -1.595998562 -0.855849197 1.2503498
[184,] -1.587718799 -0.851461710 1.2441358
[185,] -1.579439036 -0.847079126 1.2379473
[186,] -1.571159274 -0.842701675 1.2317854
[187,] -1.562879511 -0.838329586 1.2256512
[188,] -1.554599748 -0.833963086 1.2195459
[189,] -1.546319986 -0.829602406 1.2134708
[190,] -1.538040223 -0.825247776 1.2074270
[191,] -1.529760460 -0.820899427 1.2014156
[192,] -1.521480698 -0.816557591 1.1954379
[193,] -1.513200935 -0.812222499 1.1894950
[194,] -1.504921172 -0.807894386 1.1835880
[195,] -1.496641410 -0.803573486 1.1777183
[196,] -1.488361647 -0.799260033 1.1718869
[197,] -1.480081884 -0.794954265 1.1660951
[198,] -1.471802122 -0.790656418 1.1603439
[199,] -1.463522359 -0.786366731 1.1546347
[200,] -1.455242596 -0.782085445 1.1489685
[201,] -1.446962834 -0.777812801 1.1433466
[202,] -1.438683071 -0.773549042 1.1377702
[203,] -1.430403308 -0.769294412 1.1322404
[204,] -1.422123546 -0.765049159 1.1267585
[205,] -1.413843783 -0.760813529 1.1213255
[206,] -1.405564020 -0.756587774 1.1159428
[207,] -1.397284257 -0.752372145 1.1106116
[208,] -1.389004495 -0.748166897 1.1053329
[209,] -1.380724732 -0.743972286 1.1001081
[210,] -1.372444969 -0.739788571 1.0949382
[211,] -1.364165207 -0.735616014 1.0898246
[212,] -1.355885444 -0.731454880 1.0847685
[213,] -1.347605681 -0.727305434 1.0797709
[214,] -1.339325919 -0.723167947 1.0748333
[215,] -1.331046156 -0.719042693 1.0699567
[216,] -1.322766393 -0.714929947 1.0651424
[217,] -1.314486631 -0.710829989 1.0603917
[218,] -1.306206868 -0.706743102 1.0557057
[219,] -1.297927105 -0.702669574 1.0510856
[220,] -1.289647343 -0.698609695 1.0465328
[221,] -1.281367580 -0.694563760 1.0420485
[222,] -1.273087817 -0.690532068 1.0376338
[223,] -1.264808055 -0.686514922 1.0332901
[224,] -1.256528292 -0.682512631 1.0290186
[225,] -1.248248529 -0.678525507 1.0248205
[226,] -1.239968767 -0.674553868 1.0206972
[227,] -1.231689004 -0.670598036 1.0166499
[228,] -1.223409241 -0.666658340 1.0126798
[229,] -1.215129479 -0.662735111 1.0087882
[230,] -1.206849716 -0.658828691 1.0049765
[231,] -1.198569953 -0.654939422 1.0012458
[232,] -1.190290191 -0.651067656 0.9975976
[233,] -1.182010428 -0.647213750 0.9940330
[234,] -1.173730665 -0.643378066 0.9905535
[235,] -1.165450903 -0.639560975 0.9871603
[236,] -1.157171140 -0.635762852 0.9838547
[237,] -1.148891377 -0.631984082 0.9806382
[238,] -1.140611614 -0.628225053 0.9775119
[239,] -1.132331852 -0.624486164 0.9744772
[240,] -1.124052089 -0.620767819 0.9715356
[241,] -1.115772326 -0.617070430 0.9686883
[242,] -1.107492564 -0.613394418 0.9659367
[243,] -1.099212801 -0.609740210 0.9632823
[244,] -1.090933038 -0.606108242 0.9607262
[245,] -1.082653276 -0.602498957 0.9582701
[246,] -1.074373513 -0.598912808 0.9559152
[247,] -1.066093750 -0.595350254 0.9536629
[248,] -1.057813988 -0.591811765 0.9515148
[249,] -1.049534225 -0.588297817 0.9494721
[250,] -1.041254462 -0.584808896 0.9475363
[251,] -1.032974700 -0.581345496 0.9457089
[252,] -1.024694937 -0.577908121 0.9439913
[253,] -1.016415174 -0.574497282 0.9423850
[254,] -1.008135412 -0.571113499 0.9408914
[255,] -0.999855649 -0.567757302 0.9395120
[256,] -0.991575886 -0.564429229 0.9382482
[257,] -0.983296124 -0.561129826 0.9371017
[258,] -0.975016361 -0.557859649 0.9360738
[259,] -0.966736598 -0.554619263 0.9351661
[260,] -0.958456836 -0.551409238 0.9343801
[261,] -0.950177073 -0.548230157 0.9337174
[262,] -0.941897310 -0.545082608 0.9331794
[263,] -0.933617548 -0.541967189 0.9327676
[264,] -0.925337785 -0.538884505 0.9324838
[265,] -0.917058022 -0.535835170 0.9323293
[266,] -0.908778259 -0.532819803 0.9323058
[267,] -0.900498497 -0.529839034 0.9324149
[268,] -0.892218734 -0.526893497 0.9326581
[269,] -0.883938971 -0.523983834 0.9330370
[270,] -0.875659209 -0.521110693 0.9335532
[271,] -0.867379446 -0.518274730 0.9342083
[272,] -0.859099683 -0.515476604 0.9350039
[273,] -0.850819921 -0.512716982 0.9359416
[274,] -0.842540158 -0.509996533 0.9370230
[275,] -0.834260395 -0.507315935 0.9382498
[276,] -0.825980633 -0.504675865 0.9396236
[277,] -0.817700870 -0.502077009 0.9411460
[278,] -0.809421107 -0.499520052 0.9428186
[279,] -0.801141345 -0.497005684 0.9446430
[280,] -0.792861582 -0.494534596 0.9466210
[281,] -0.784581819 -0.492107482 0.9487541
[282,] -0.776302057 -0.489725036 0.9510440
[283,] -0.768022294 -0.487387953 0.9534923
[284,] -0.759742531 -0.485096929 0.9561006
[285,] -0.751462769 -0.482852658 0.9588707
[286,] -0.743183006 -0.480655834 0.9618041
[287,] -0.734903243 -0.478507148 0.9649025
[288,] -0.726623481 -0.476407289 0.9681674
[289,] -0.718343718 -0.474356945 0.9716006
[290,] -0.710063955 -0.472356798 0.9752037
[291,] -0.701784193 -0.470407526 0.9789782
[292,] -0.693504430 -0.468509804 0.9829259
[293,] -0.685224667 -0.466664302 0.9870482
[294,] -0.676944905 -0.464871680 0.9913469
[295,] -0.668665142 -0.463132597 0.9958235
[296,] -0.660385379 -0.461447701 1.0004796
[297,] -0.652105616 -0.459817635 1.0053168
[298,] -0.643825854 -0.458243031 1.0103367
[299,] -0.635546091 -0.456724516 1.0155409
[300,] -0.627266328 -0.455262704 1.0209309
[301,] -0.618986566 -0.453858177 1.0265081
[302,] -0.610706803 -0.452510912 1.0322710
[303,] -0.602427040 -0.451220257 1.0382151
[304,] -0.594147278 -0.449985489 1.0443357
[305,] -0.585867515 -0.448805835 1.0506281
[306,] -0.577587752 -0.447680475 1.0570875
[307,] -0.569307990 -0.446608540 1.0637091
[308,] -0.561028227 -0.445589119 1.0704881
[309,] -0.552748464 -0.444621252 1.0774194
[310,] -0.544468702 -0.443703940 1.0844982
[311,] -0.536188939 -0.442836136 1.0917195
[312,] -0.527909176 -0.442016755 1.0990781
[313,] -0.519629414 -0.441244674 1.1065691
[314,] -0.511349651 -0.440518728 1.1141873
[315,] -0.503069888 -0.439837717 1.1219275
[316,] -0.494790126 -0.439200408 1.1297847
[317,] -0.486510363 -0.438605532 1.1377534
[318,] -0.478230600 -0.438051789 1.1458285
[319,] -0.469950838 -0.437537850 1.1540048
[320,] -0.461671075 -0.437062358 1.1622768
[321,] -0.453391312 -0.436623929 1.1706393
[322,] -0.445111550 -0.436221155 1.1790868
[323,] -0.436831787 -0.435852604 1.1876141
[324,] -0.428552024 -0.435516825 1.1962156
[325,] -0.420272261 -0.435212345 1.2048860
[326,] -0.411992499 -0.434937675 1.2136198
[327,] -0.403712736 -0.434691309 1.2224115
[328,] -0.395432973 -0.434471727 1.2312558
[329,] -0.387153211 -0.434277397 1.2401470
[330,] -0.378873448 -0.434106773 1.2490797
[331,] -0.370593685 -0.433958303 1.2580484
[332,] -0.362313923 -0.433830423 1.2670475
[333,] -0.354034160 -0.433721563 1.2760716
[334,] -0.345754397 -0.433630150 1.2851151
[335,] -0.337474635 -0.433554603 1.2941724
[336,] -0.329194872 -0.433493341 1.3032381
[337,] -0.320915109 -0.433444779 1.3123065
[338,] -0.312635347 -0.433407334 1.3213723
[339,] -0.304355584 -0.433379421 1.3304297
[340,] -0.296075821 -0.433359460 1.3394733
[341,] -0.287796059 -0.433345872 1.3484975
[342,] -0.279516296 -0.433337081 1.3574969
[343,] -0.271236533 -0.433331519 1.3664658
[344,] -0.262956771 -0.433327623 1.3753988
[345,] -0.254677008 -0.433323835 1.3842903
[346,] -0.246397245 -0.433318607 1.3931349
[347,] -0.238117483 -0.433310401 1.4019270
[348,] -0.229837720 -0.433297686 1.4106612
[349,] -0.221557957 -0.433278944 1.4193320
[350,] -0.213278195 -0.433252667 1.4279339
[351,] -0.204998432 -0.433217361 1.4364615
[352,] -0.196718669 -0.433171543 1.4449094
[353,] -0.188438906 -0.433113747 1.4532720
[354,] -0.180159144 -0.433042520 1.4615441
[355,] -0.171879381 -0.432956424 1.4697202
[356,] -0.163599618 -0.432854040 1.4777949
[357,] -0.155319856 -0.432733964 1.4857629
[358,] -0.147040093 -0.432594812 1.4936188
[359,] -0.138760330 -0.432435217 1.5013574
[360,] -0.130480568 -0.432253835 1.5089732
[361,] -0.122200805 -0.432049339 1.5164611
[362,] -0.113921042 -0.431820426 1.5238158
[363,] -0.105641280 -0.431565815 1.5310320
[364,] -0.097361517 -0.431284248 1.5381046
[365,] -0.089081754 -0.430974493 1.5450282
[366,] -0.080801992 -0.430635340 1.5517979
[367,] -0.072522229 -0.430265609 1.5584084
[368,] -0.064242466 -0.429864145 1.5648546
[369,] -0.055962704 -0.429429824 1.5711314
[370,] -0.047682941 -0.428961548 1.5772339
[371,] -0.039403178 -0.428458255 1.5831569
[372,] -0.031123416 -0.427918910 1.5888955
[373,] -0.022843653 -0.427342515 1.5944447
[374,] -0.014563890 -0.426728105 1.5997996
[375,] -0.006284128 -0.426074754 1.6049554
[376,] 0.001995635 -0.425381571 1.6099071
[377,] 0.010275398 -0.424647707 1.6146500
[378,] 0.018555160 -0.423872352 1.6191793
[379,] 0.026834923 -0.423054741 1.6234903
[380,] 0.035114686 -0.422194154 1.6275783
[381,] 0.043394448 -0.421289917 1.6314387
[382,] 0.051674211 -0.420341541 1.6350680
[383,] 0.059953974 -0.419349375 1.6384694
[384,] 0.068233737 -0.418314091 1.6416488
[385,] 0.076513499 -0.417236369 1.6446117
[386,] 0.084793262 -0.416116902 1.6473640
[387,] 0.093073025 -0.414956391 1.6499114
[388,] 0.101352787 -0.413755541 1.6522598
[389,] 0.109632550 -0.412515061 1.6544149
[390,] 0.117912313 -0.411235664 1.6563824
[391,] 0.126192075 -0.409918058 1.6581682
[392,] 0.134471838 -0.408562951 1.6597780
[393,] 0.142751601 -0.407171045 1.6612176
[394,] 0.151031363 -0.405743037 1.6624927
[395,] 0.159311126 -0.404279614 1.6636091
[396,] 0.167590889 -0.402781452 1.6645726
[397,] 0.175870651 -0.401249216 1.6653888
[398,] 0.184150414 -0.399683557 1.6660635
[399,] 0.192430177 -0.398085108 1.6666025
[400,] 0.200709939 -0.396454487 1.6670112
[401,] 0.208989702 -0.394792290 1.6672956
[402,] 0.217269465 -0.393099094 1.6674611
[403,] 0.225549227 -0.391375451 1.6675134
[404,] 0.233828990 -0.389621890 1.6674581
[405,] 0.242108753 -0.387838914 1.6673007
[406,] 0.250388515 -0.386026995 1.6670468
[407,] 0.258668278 -0.384186578 1.6667020
[408,] 0.266948041 -0.382318077 1.6662716
[409,] 0.275227803 -0.380421872 1.6657612
[410,] 0.283507566 -0.378498309 1.6651761
[411,] 0.291787329 -0.376547698 1.6645217
[412,] 0.300067092 -0.374570311 1.6638034
[413,] 0.308346854 -0.372566383 1.6630265
[414,] 0.316626617 -0.370536105 1.6621962
[415,] 0.324906380 -0.368479631 1.6613177
[416,] 0.333186142 -0.366397067 1.6603962
[417,] 0.341465905 -0.364288478 1.6594368
[418,] 0.349745668 -0.362153882 1.6584447
[419,] 0.358025430 -0.359993249 1.6574248
[420,] 0.366305193 -0.357806503 1.6563822
[421,] 0.374584956 -0.355593515 1.6553217
[422,] 0.382864718 -0.353354110 1.6542484
[423,] 0.391144481 -0.351088059 1.6531669
[424,] 0.399424244 -0.348795080 1.6520821
[425,] 0.407704006 -0.346474839 1.6509988
[426,] 0.415983769 -0.344126949 1.6499216
[427,] 0.424263532 -0.341750966 1.6488551
[428,] 0.432543294 -0.339346392 1.6478039
[429,] 0.440823057 -0.336912673 1.6467725
[430,] 0.449102820 -0.334449198 1.6457654
[431,] 0.457382582 -0.331955301 1.6447869
[432,] 0.465662345 -0.329430259 1.6438414
[433,] 0.473942108 -0.326873289 1.6429332
[434,] 0.482221870 -0.324283557 1.6420666
[435,] 0.490501633 -0.321660166 1.6412456
[436,] 0.498781396 -0.319002168 1.6404744
[437,] 0.507061158 -0.316308556 1.6397571
[438,] 0.515340921 -0.313578268 1.6390977
[439,] 0.523620684 -0.310810188 1.6385001
[440,] 0.531900447 -0.308003146 1.6379682
[441,] 0.540180209 -0.305155920 1.6375059
[442,] 0.548459972 -0.302267234 1.6371169
[443,] 0.556739735 -0.299335765 1.6368050
[444,] 0.565019497 -0.296360138 1.6365738
[445,] 0.573299260 -0.293338933 1.6364271
[446,] 0.581579023 -0.290270683 1.6363685
[447,] 0.589858785 -0.287153879 1.6364014
[448,] 0.598138548 -0.283986971 1.6365294
[449,] 0.606418311 -0.280768369 1.6367560
[450,] 0.614698073 -0.277496447 1.6370846
[451,] 0.622977836 -0.274169547 1.6375186
[452,] 0.631257599 -0.270785978 1.6380614
[453,] 0.639537361 -0.267344024 1.6387164
[454,] 0.647817124 -0.263841944 1.6394869
[455,] 0.656096887 -0.260277974 1.6403761
[456,] 0.664376649 -0.256650336 1.6413873
[457,] 0.672656412 -0.252957237 1.6425239
[458,] 0.680936175 -0.249196877 1.6437892
[459,] 0.689215937 -0.245367449 1.6451862
[460,] 0.697495700 -0.241467146 1.6467184
[461,] 0.705775463 -0.237494165 1.6483890
[462,] 0.714055225 -0.233446718 1.6502012
[463,] 0.722334988 -0.229323595 1.6521566
[464,] 0.730614751 -0.225124787 1.6542534
[465,] 0.738894513 -0.220850491 1.6564895
[466,] 0.747174276 -0.216500965 1.6588626
[467,] 0.755454039 -0.212076527 1.6613707
[468,] 0.763733801 -0.207577554 1.6640118
[469,] 0.772013564 -0.203004484 1.6667840
[470,] 0.780293327 -0.198357807 1.6696853
[471,] 0.788573090 -0.193638069 1.6727139
[472,] 0.796852852 -0.188845871 1.6758681
[473,] 0.805132615 -0.183981863 1.6791460
[474,] 0.813412378 -0.179046746 1.6825460
[475,] 0.821692140 -0.174041269 1.6860665
[476,] 0.829971903 -0.168966226 1.6897058
[477,] 0.838251666 -0.163822459 1.6934625
[478,] 0.846531428 -0.158610852 1.6973350
[479,] 0.854811191 -0.153332329 1.7013218
[480,] 0.863090954 -0.147987857 1.7054216
[481,] 0.871370716 -0.142578439 1.7096329
[482,] 0.879650479 -0.137105116 1.7139545
[483,] 0.887930242 -0.131568964 1.7183849
[484,] 0.896210004 -0.125971094 1.7229229
[485,] 0.904489767 -0.120312646 1.7275673
[486,] 0.912769530 -0.114594792 1.7323168
[487,] 0.921049292 -0.108818735 1.7371703
[488,] 0.929329055 -0.102985700 1.7421266
[489,] 0.937608818 -0.097096943 1.7471846
[490,] 0.945888580 -0.091153741 1.7523432
[491,] 0.954168343 -0.085157394 1.7576012
[492,] 0.962448106 -0.079109223 1.7629576
[493,] 0.970727868 -0.073010569 1.7684114
[494,] 0.979007631 -0.066862791 1.7739615
[495,] 0.987287394 -0.060667262 1.7796069
[496,] 0.995567156 -0.054425374 1.7853467
[497,] 1.003846919 -0.048138528 1.7911798
[498,] 1.012126682 -0.041808141 1.7971052
[499,] 1.020406445 -0.035435638 1.8031221
[500,] 1.028686207 -0.029022453 1.8092294
[501,] 1.036965970 -0.022570030 1.8154262
[502,] 1.045245733 -0.016079817 1.8217116
[503,] 1.053525495 -0.009553267 1.8280846
[504,] 1.061805258 -0.002991839 1.8345443
[505,] 1.070085021 0.003603010 1.8410898
[506,] 1.078364783 0.010229817 1.8477201
[507,] 1.086644546 0.016887123 1.8544344
[508,] 1.094924309 0.023573470 1.8612316
[509,] 1.103204071 0.030287401 1.8681109
[510,] 1.111483834 0.037027465 1.8750713
[511,] 1.119763597 0.043792215 1.8821119
[512,] 1.128043359 0.050580211 1.8892318
[513,] 1.136323122 0.057390021 1.8964299
[514,] 1.144602885 0.064220220 1.9037053
[515,] 1.152882647 0.071069394 1.9110571
[516,] 1.161162410 0.077936139 1.9184842
[517,] 1.169442173 0.084819064 1.9259857
[518,] 1.177721935 0.091716788 1.9335605
[519,] 1.186001698 0.098627947 1.9412077
[520,] 1.194281461 0.105551188 1.9489262
[521,] 1.202561223 0.112485177 1.9567149
[522,] 1.210840986 0.119428594 1.9645728
[523,] 1.219120749 0.126380137 1.9724988
[524,] 1.227400511 0.133338521 1.9804918
[525,] 1.235680274 0.140302480 1.9885507
[526,] 1.243960037 0.147270768 1.9966743
[527,] 1.252239800 0.154242159 2.0048614
[528,] 1.260519562 0.161215446 2.0131109
[529,] 1.268799325 0.168189445 2.0214216
[530,] 1.277079088 0.175162994 2.0297923
[531,] 1.285358850 0.182134951 2.0382217
[532,] 1.293638613 0.189104199 2.0467085
[533,] 1.301918376 0.196069644 2.0552514
[534,] 1.310198138 0.203030214 2.0638492
[535,] 1.318477901 0.209984864 2.0725005
[536,] 1.326757664 0.216932571 2.0812039
[537,] 1.335037426 0.223872337 2.0899580
[538,] 1.343317189 0.230803189 2.0987614
[539,] 1.351596952 0.237724181 2.1076127
[540,] 1.359876714 0.244634390 2.1165105
[541,] 1.368156477 0.251532919 2.1254531
[542,] 1.376436240 0.258418898 2.1344392
[543,] 1.384716002 0.265291481 2.1434671
[544,] 1.392995765 0.272149849 2.1525353
[545,] 1.401275528 0.278993208 2.1616422
[546,] 1.409555290 0.285820791 2.1707862
[547,] 1.417835053 0.292631855 2.1799657
[548,] 1.426114816 0.299425684 2.1891789
[549,] 1.434394578 0.306201587 2.1984241
[550,] 1.442674341 0.312958898 2.2076997
[551,] 1.450954104 0.319696978 2.2170040
[552,] 1.459233866 0.326415210 2.2263350
[553,] 1.467513629 0.333113004 2.2356911
[554,] 1.475793392 0.339789794 2.2450703
[555,] 1.484073154 0.346445166 2.2544713
[556,] 1.492352917 0.353079131 2.2638937
[557,] 1.500632680 0.359691783 2.2733375
[558,] 1.508912443 0.366283216 2.2828024
[559,] 1.517192205 0.372853524 2.2922885
[560,] 1.525471968 0.379402805 2.3017956
[561,] 1.533751731 0.385931155 2.3113236
[562,] 1.542031493 0.392438672 2.3208724
[563,] 1.550311256 0.398925454 2.3304420
[564,] 1.558591019 0.405391601 2.3400323
[565,] 1.566870781 0.411837211 2.3496430
[566,] 1.575150544 0.418262386 2.3592742
[567,] 1.583430307 0.424667227 2.3689258
[568,] 1.591710069 0.431051836 2.3785975
[569,] 1.599989832 0.437416314 2.3882894
[570,] 1.608269595 0.443760765 2.3980014
[571,] 1.616549357 0.450085291 2.4077332
[572,] 1.624829120 0.456389997 2.4174849
[573,] 1.633108883 0.462674987 2.4272563
[574,] 1.641388645 0.468940367 2.4370473
[575,] 1.649668408 0.475186240 2.4468578
[576,] 1.657948171 0.481412712 2.4566877
[577,] 1.666227933 0.487619891 2.4665369
[578,] 1.674507696 0.493807881 2.4764052
[579,] 1.682787459 0.499976790 2.4862927
[580,] 1.691067221 0.506126724 2.4961991
[581,] 1.699346984 0.512257791 2.5061245
[582,] 1.707626747 0.518370099 2.5160685
[583,] 1.715906509 0.524463755 2.5260312
[584,] 1.724186272 0.530538866 2.5360125
[585,] 1.732466035 0.536595542 2.5460122
[586,] 1.740745798 0.542633890 2.5560302
[587,] 1.749025560 0.548654019 2.5660665
[588,] 1.757305323 0.554656038 2.5761208
[589,] 1.765585086 0.560640055 2.5861932
[590,] 1.773864848 0.566606178 2.5962834
[591,] 1.782144611 0.572554518 2.6063915
[592,] 1.790424374 0.578485183 2.6165172
[593,] 1.798704136 0.584398281 2.6266605
[594,] 1.806983899 0.590293921 2.6368212
[595,] 1.815263662 0.596172213 2.6469993
[596,] 1.823543424 0.602033266 2.6571946
[597,] 1.831823187 0.607877188 2.6674071
[598,] 1.840102950 0.613704087 2.6776365
[599,] 1.848382712 0.619514074 2.6878829
[600,] 1.856662475 0.625307256 2.6981461
Let’s start with splines with 5 degrees of freedom for each predictor
require(rms)
fit_all <- cph(Surv(TEVENT, vasc) ~
rcs(AGE_std, 5) +
rcs(BMI_std, 5) +
rcs(HDL_std, 5) +
rcs(log_creat_std, 5) +
rcs(IMT_std, 5) +
SEX + CARDIAC,
data = dev,
x = T, y = T, surv = T)
anova(fit_all, test = "Chisq")
Wald Statistics Response: Surv(TEVENT, vasc)
Factor Chi-Square d.f. P
AGE_std 38.17 4 <.0001
Nonlinear 5.75 3 0.1243
BMI_std 6.89 4 0.1420
Nonlinear 6.27 3 0.0992
HDL_std 4.41 4 0.3531
Nonlinear 0.61 3 0.8931
log_creat_std 54.55 4 <.0001
Nonlinear 3.10 3 0.3761
IMT_std 28.55 4 <.0001
Nonlinear 0.78 3 0.8538
SEX 0.86 1 0.3534
CARDIAC 1.69 1 0.1934
TOTAL NONLINEAR 17.31 15 0.3009
TOTAL 200.68 22 <.0001
For Age, bmi, hdl, creat and imt, the non-linear effects do not seem significant.
fit <- cph(Surv(TEVENT, vasc) ~
rcs(AGE_std, 4) +
rcs(BMI_std, 4) +
rcs(HDL_std, 4) +
rcs(log_creat_std, 4) +
rcs(IMT_std, 4) +
SEX + CARDIAC,
data = dev,
x = T, y = T, surv = T)
anova(fit, test = "Chisq")
Wald Statistics Response: Surv(TEVENT, vasc)
Factor Chi-Square d.f. P
AGE_std 37.41 3 <.0001
Nonlinear 3.94 2 0.1392
BMI_std 6.56 3 0.0875
Nonlinear 5.99 2 0.0501
HDL_std 3.78 3 0.2858
Nonlinear 0.07 2 0.9633
log_creat_std 53.12 3 <.0001
Nonlinear 2.45 2 0.2938
IMT_std 27.42 3 <.0001
Nonlinear 0.33 2 0.8473
SEX 0.89 1 0.3465
CARDIAC 1.60 1 0.2053
TOTAL NONLINEAR 13.39 10 0.2024
TOTAL 198.46 17 <.0001
fit <- cph(Surv(TEVENT, vasc) ~
rcs(AGE_std, 3) +
rcs(BMI_std, 3) +
rcs(HDL_std, 3) +
rcs(log_creat_std, 3) +
rcs(IMT_std, 3) +
SEX + CARDIAC,
data = dev,
x = T, y = T, surv = T)
anova(fit, test = "Chisq")
Wald Statistics Response: Surv(TEVENT, vasc)
Factor Chi-Square d.f. P
AGE_std 38.78 2 <.0001
Nonlinear 3.88 1 0.0490
BMI_std 2.40 2 0.3016
Nonlinear 2.10 1 0.1471
HDL_std 3.84 2 0.1464
Nonlinear 0.03 1 0.8671
log_creat_std 54.53 2 <.0001
Nonlinear 1.22 1 0.2700
IMT_std 27.00 2 <.0001
Nonlinear 0.15 1 0.7023
SEX 0.88 1 0.3479
CARDIAC 1.37 1 0.2414
TOTAL NONLINEAR 7.91 5 0.1611
TOTAL 194.28 12 <.0001
For age, we should keep the splines
fit <- cph(Surv(TEVENT, vasc) ~
rcs(AGE_std, 3) +
BMI_std +
HDL_std +
log_creat_std +
IMT_std +
SEX + CARDIAC,
data = dev,
x = T, y = T, surv = T)
anova(fit, test = "Chisq")
Wald Statistics Response: Surv(TEVENT, vasc)
Factor Chi-Square d.f. P
AGE_std 37.70 2 <.0001
Nonlinear 4.22 1 0.0401
BMI_std 0.37 1 0.5432
HDL_std 3.68 1 0.0551
log_creat_std 50.02 1 <.0001
IMT_std 26.08 1 <.0001
SEX 0.13 1 0.7200
CARDIAC 1.10 1 0.2939
TOTAL 185.82 8 <.0001
We will keep all variables in the model, and only for age a non-linear term
Check the PH assumption and residuals, and assess overall model fit. Which are the most important predictors in the model? EXTRA: Fit the model again, but now using the standardized version of the continuous variables. Where do you see the differences?
PH assumption is proportionality of linear predictor for all time points
We can test this with cox.zph, but this does not work for cph objects
cox.zph(coxph(Surv(TEVENT, vasc) ~
rcs(AGE_std, 3) +
BMI_std +
HDL_std +
log_creat_std +
IMT_std +
SEX + CARDIAC,
data = dev)) %>%
plot()
The betas seem pretty equal accross time-points, except for cardiac
Check overall fit
anova(fit, test = "Chisq")
Wald Statistics Response: Surv(TEVENT, vasc)
Factor Chi-Square d.f. P
AGE_std 37.70 2 <.0001
Nonlinear 4.22 1 0.0401
BMI_std 0.37 1 0.5432
HDL_std 3.68 1 0.0551
log_creat_std 50.02 1 <.0001
IMT_std 26.08 1 <.0001
SEX 0.13 1 0.7200
CARDIAC 1.10 1 0.2939
TOTAL 185.82 8 <.0001
Removing log_creat_std
will hurt the model fit the most.
Overall fit is significant.
Let us assume that we are particularly interested in 3yr and 5yr risk estimates. Extract these predicted probabilities from your model for the patients in the development set. (Suggestion: look at the equation in slide 30)
Get survival function
surv_fun <- Survival(fit)
pred_survivals <- surv_fun(times = c(3,5) * 365, lp = fit$linear.predictors)
Plot baseline survival function for average patient with all covariates zero (which is not actually average since the categorical variables are not standardized)
base_surv <- surv_fun(times = seq(from = min(dev$TEVENT), to = max(dev$TEVENT),
by = 1), lp = 0)
plot(basehaz(fit)[, 2], 1-basehaz(fit)[, 1], type = "l", col = "blue")
lines(base_surv, type = "l")
Or get baseline hazard
base_haz <- basehaz(fit)
plot(base_haz$time, base_haz$hazard, type = "l")
EXTRA: Use the same approach, but then for the model using the standardized variables. Where do you see the differences?
Assess discrimination using Harrell’s c statistic, both for the overall model, and for the 3 and 5 year predicted probabilities specifically.
rcorr.cens(x = -fit$linear.predictors, S = fit$y)[1]
C Index
0.6861534
rcorr.cens(x = -pred_survivals[,1], S = dev$vasc & dev$TEVENT < 3*365)[1]
C Index
0.7042643
rcorr.cens(x = -pred_survivals[,2], S = dev$vasc & dev$TEVENT < 5*365)[1]
C Index
0.7105941
With survAUC
require(survAUC)
chambless <- AUC.cd(
Surv.rsp = fit$y,
Surv.rsp.new = fit$y,
lp = fit$linear.predictors,
lpnew = fit$linear.predictors,
times = seq(1, 10*365, by = 1))
plot(chambless)
Assess calibration of the model by plotting the calibration plot using the groupkm function, for the 3 and 5 yr predictions. EXTRA: use the survAUC package to calculate the c-statistic using Chambless c at 3 years and 5 years, and to calculate the c statistic over time. Calculate the prediction error over time using the functions in the pec package. 5. Presentation of the model Present the results of the model in a risk score or a nomogram B MODEL VALIDATION To validate our model, we use the dataset SMARTval.sav. 1. Check characteristics Import the dataset SMARTval.sav and check the characteristics of the predictors and the events plus timing. Compare to the development dataset. Which differences do you see? Where do you think this might lead to problems? 2. Calculate predicted probabilities i) Calculate predicted probabilities for those patients using the model based on the original variables, for 3yr and 5 yr predictions. ii) Calculate these predicted probabilities after updating the average patient in your model. Now compare the probabilities obtained by i) and ii). Calculate the c statistics based on i) and ii) and compare. Create calibration plots based on i) and ii) and compare. 3. Recalibration and revision i) Recalibrate the model by updating the baseline hazard and the ‘average patient’. Calculate again the predicted probabilities for 3 yrs and 5 yrs, and the performance measures. ii) Revise the model by updating the baseline hazard, the ‘average patient’ and the slope. Calculate again the predicted probabilities for 3yrs and 5 yr, and the performance measures.
sessionInfo()
R version 3.4.3 (2017-11-30)
Platform: x86_64-apple-darwin15.6.0 (64-bit)
Running under: macOS Sierra 10.12.6
Matrix products: default
BLAS: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRblas.0.dylib
LAPACK: /Library/Frameworks/R.framework/Versions/3.4/Resources/lib/libRlapack.dylib
locale:
[1] en_US.UTF-8/en_US.UTF-8/en_US.UTF-8/C/en_US.UTF-8/en_US.UTF-8
attached base packages:
[1] stats graphics grDevices utils datasets methods base
other attached packages:
[1] survAUC_1.0-5 ggfortify_0.4.2 logistf_1.22
[4] rms_5.1-2 SparseM_1.77 Hmisc_4.1-1
[7] Formula_1.2-2 survival_2.41-3 lattice_0.20-35
[10] glmnet_2.0-13 foreach_1.4.4 Matrix_1.2-12
[13] tidyr_0.8.0 bindrcpp_0.2 haven_1.1.1
[16] broom_0.4.3 epistats_0.1.0 ggplot2_2.2.1
[19] here_0.1 purrr_0.2.4 magrittr_1.5
[22] data.table_1.10.4-3 dplyr_0.7.4
loaded via a namespace (and not attached):
[1] binom_1.1-1 splines_3.4.3
[3] gtools_3.5.0 assertthat_0.2.0
[5] latticeExtra_0.6-28 yaml_2.1.16
[7] pillar_1.1.0 backports_1.1.2
[9] quantreg_5.35 glue_1.2.0
[11] pROC_1.10.0 digest_0.6.15
[13] RColorBrewer_1.1-2 checkmate_1.8.5
[15] ResourceSelection_0.3-2 colorspace_1.3-2
[17] sandwich_2.4-0 htmltools_0.3.6
[19] plyr_1.8.4 psych_1.7.8
[21] pkgconfig_2.0.1 gmodels_2.16.2
[23] mvtnorm_1.0-7 scales_0.5.0
[25] gdata_2.18.0 MatrixModels_0.4-1
[27] git2r_0.21.0 htmlTable_1.11.2
[29] tibble_1.4.2 mgcv_1.8-23
[31] TH.data_1.0-8 nnet_7.3-12
[33] lazyeval_0.2.1 cli_1.0.0
[35] mnormt_1.5-5 crayon_1.3.4
[37] polspline_1.1.12 evaluate_0.10.1
[39] mice_2.46.0 nlme_3.1-131
[41] MASS_7.3-48 forcats_0.2.0
[43] foreign_0.8-69 tools_3.4.3
[45] hms_0.4.1 multcomp_1.4-8
[47] stringr_1.2.0 munsell_0.4.3
[49] cluster_2.0.6 compiler_3.4.3
[51] rlang_0.1.6 grid_3.4.3
[53] iterators_1.0.9 rstudioapi_0.7
[55] htmlwidgets_1.0 labeling_0.3
[57] base64enc_0.1-3 rmarkdown_1.8
[59] codetools_0.2-15 gtable_0.2.0
[61] reshape2_1.4.3 R6_2.2.2
[63] gridExtra_2.3 zoo_1.8-1
[65] knitr_1.19 utf8_1.1.3
[67] bindr_0.1 rprojroot_1.3-2
[69] readr_1.1.1 stringi_1.1.6
[71] parallel_3.4.3 Rcpp_0.12.15
[73] rpart_4.1-12 acepack_1.4.1
[75] tidyselect_0.2.3
This R Markdown site was created with workflowr