a)

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 

b)

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

c)

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 

d)

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 | 
----------------------|-----------|-----------|-----------|

 

e)

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"

f)

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).

g)

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

h)

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

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

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

i)

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.

a)

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  

2. 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.

c)

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

d)

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

e)

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

f)

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 

g)

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

h)

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))

i)

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

j)

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)

k)

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.

dependencies

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:

load development and validation data

tbi2 <- readRDS(here("data", "TBIday3.RDS"))
devdata <- tbi2[tbi2$trial=="Tirilazad US",-1]
valdata <- tbi2[tbi2$trial=="Tirilazad International",-1]

Descriptive statistics

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.

Q2

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.

Q3:

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

Q4

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.

Q5

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).

Q6

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

Q7

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.

Q8

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 

Q9

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.

Q10

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).

Background

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.

Q11

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

Q13

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

Q14

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.

Q15

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

Q16

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 

1.

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

Check survival and censoring distribution

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)

Create standardized variables

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

Try to plot RCsplines

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.