Model Distillation: Kidney SOFA Example

author:

• Michael Hagmann
• Stefan Riezler

date: 30. April 2021

download: https://github.com/StatNLP/empirical_methods

---

knitr::opts_chunk$set(echo = TRUE, 
                      cache.lazy = FALSE,
                      dev = c("svglite", "pdf", "png"),
                      dpi = 300,
                      fig.path = 'figures/',
                      fig.keep = "high")

#added from: https://github.com/yihui/knitr-examples/blob/master/077-wrap-output.Rmd
library(knitr)
hook_output = knit_hooks$get('output')
knit_hooks$set(output = function(x, options) {
  # this hook is used only when the linewidth option is not NULL
  if (!is.null(n <- options$linewidth)) {
    x = knitr:::split_lines(x)
    # any lines wider than n should be wrapped
    if (any(nchar(x) > n)) x = strwrap(x, width = n)
    x = paste(x, collapse = '\n')
  }
  hook_output(x, options)
})

library(tidyverse) # A collection of packages for data science. More about it on 
                   # www.tidyverse.com
library(magrittr)  # A package that provides pipe operators like %>%
library(mgcv)      # A package written by Simon Wood that implements his 
                   # (spline based) GAM specification.
library(glue)      # A package that provides interpolated string functions.

#When you want to reproduce the examples, you have to run this code chunk in advance!!
powerSet <- function(set) {
  
    apply(X      = rep(list(c(T,F)), length = length(set)) %>% expand.grid(.),
          MARGIN = 1,
          FUN    = function(...) set[...])
}

accuracy <- 
  function(y, y_hat, filter = TRUE) mean((y == y_hat)[filter])

class_accuracy <-
  function(y, y_hat, classes = sort(unique(y))) 
    sapply(as.character(classes), function(c) accuracy(y, y_hat, filter = (y == as.numeric(c))))

train_data <- readRDS("data/train_set_with_sofa_predictions.rds")

test_data <- readRDS("data/test_set_with_sofa_predictions.rds")

Learner: Feed forward neural network with seven layers (number of neurons per layer: 128, 161, 203 , 256, 203 , 161, 128) and relu activation function. The network was trained with SGD (batch size 64) and learning rate .01 for 5 epochs and a dropoutrate of .2 for inter hidden layer weights. All other optimizer settings are default values of pyTorch's SGD optimizer. The objective function was MSE.

In a second step thresholds to turn the network output (score) into class predictions were learned . For this purpose a ordinal regression model was trained (R package "ordinal").

Train set size: nrow(train_data) measurement points. Test set size: nrow(test_data) measurement points.

Performance Evaluation

Circular Model - Train set accuracy (overall and per class)

with(train_data, accuracy(target, sofa_circ))
with(train_data, class_accuracy(target, sofa_circ))

0.956886742279007

0

0.994546390142152

1

0.966884840717686

2

0.817238679368394

3

0.635302617684655

4

0.932277710109622

Confusion Matrix:

xtabs(~ target + sofa_circ, data = train_data)

      sofa_circ
target      0      1      2      3      4
     0 216467    901    142    142      2
     1   1083  42249    307     57      0
     2      0    669  12525   2119     13
     3    778   1494   1814  11431   2476
     4      2     16     14   1914  26789

Circular Model - Test set accuracy (overall and per class)

with(test_data, accuracy(target, sofa_circ))
with(test_data, class_accuracy(target, sofa_circ))

0.905579452343469

0

0.98656943839791

1

0.90375219855384

2

0.739289558665231

3

0.656794998354722

4

0.843968327899394

Confusion Matrix:

xtabs(~ target + sofa_circ, data = test_data)

      sofa_circ
target     0     1     2     3     4
     0 45323   468    63    83     3
     1   681  9249   244    53     7
     2     0   453  3434   753     5
     3   358   500   874  5988  1397
     4     1     2    15  1657  9060

Non-circular Model - Train set accuracy (overall and per class)

with(train_data, accuracy(target, sofa_nocirc))
with(train_data, class_accuracy(target, sofa_nocirc))

0.963318326303942

0

0.989170885901477

1

0.951277004760161

2

0.897429205272087

3

0.752292558217084

4

0.953088567948495

Confusion Matrix:

xtabs(~ target + sofa_nocirc, data = train_data)

      sofa_nocirc
target      0      1      2      3      4
     0 215297   1619    402    336      0
     1    923  41567    807    399      0
     2      0    743  13754    819     10
     3     81    431    920  13536   3025
     4      0      0     16   1332  27387

Non-circular Model - Test set accuracy (overall and per class)

with(test_data, accuracy(target, sofa_nocirc))
with(test_data, class_accuracy(target, sofa_nocirc))

0.92558664203989

0

0.967501088376143

1

0.920949775258941

2

0.815931108719053

3

0.762641219699463

4

0.936469492314858

Confusion Matrix:

xtabs(~ target + sofa_nocirc, data = test_data)

      sofa_nocirc
target     0     1     2     3     4
     0 44447  1005   258   228     2
     1   217  9425   397   195     0
     2     0   522  3790   327     6
     3    61   284   543  6953  1276
     4     0     7     4   671 10053

Correlation pattern

library(circlize)
#to install package "ComplexHeatmap":
# install.packages("BiocManager")
# BiocManager::install("ComplexHeatmap")
library(ComplexHeatmap)
library(gridBase)

Loading required package: grid

========================================
ComplexHeatmap version 2.6.2
Bioconductor page: http://bioconductor.org/packages/ComplexHeatmap/
Github page: https://github.com/jokergoo/ComplexHeatmap
Documentation: http://jokergoo.github.io/ComplexHeatmap-reference

If you use it in published research, please cite:
Gu, Z. Complex heatmaps reveal patterns and correlations in multidimensional 
  genomic data. Bioinformatics 2016.

This message can be suppressed by:
  suppressPackageStartupMessages(library(ComplexHeatmap))
========================================

method <- "pearson"

corr_matrix <-
  cbind(
    "label"        = sapply(select(test_data, amv:horovitz), 
                            cor,
                            method = method,
                            y = test_data$target),
    "circular"     = sapply(select(test_data, amv:horovitz), 
                            cor,
                            method = method,
                            y = test_data$sofa_circ),
    "non-circular" = sapply(select(test_data, amv:horovitz), 
                            cor,
                            method = method,
                            y = test_data$sofa_nocirc))

plot.new()
circle_size = unit(1, "snpc") # snpc unit gives you a square region

pushViewport(viewport(x = 0, 
                      y = 0.5, 
                      width = circle_size, 
                      height = circle_size, 
                      just = c("left", "center")))


col_scale <- 
  colorRamp2(c(-1, 0, 1), c("#67001F", "#FFFFFF", "#053061"), space ="LAB")

par(omi = gridOMI(), new = TRUE)
par(mar = c(0.1, 0.1, 0.1, 0.1))
circos.par(gap.after = 35, cell.padding = c(0, 0, 0, 0))


  circos.heatmap(corr_matrix[, 3, drop = FALSE], 
                 col = col_scale, 
                 track.height = .1,
                 rownames.side = "outside",
                 rownames.font = 2)
                 #rownames.cex = 0.5 + .25 * rownames(corr_matrix) %in% vars)
  circos.track(track.index = 2, 
               panel.fun = function(x, y) {
                  cn = colnames(corr_matrix)[3]
                  n  = length(cn)
                  circos.text(rep(CELL_META$cell.xlim[2], n) + convert_x(1, "mm"), 
                              1:n - 0.5, cn, 
                              cex = 0.75, 
                              adj = c(0, 0.5), 
                              facing = "bending.inside",
                              niceFacing = TRUE,
                              font = 2)},
              bg.border = NA)
  
                
  circos.heatmap(corr_matrix[,2 , drop = FALSE], 
               col = col_scale, 
               track.height = .1)
  circos.track(track.index = 3,
               panel.fun = function(x, y) {
                  cn = colnames(corr_matrix)[2]
                  n = length(cn)
                  circos.text(rep(CELL_META$cell.xlim[2], n) + convert_x(1, "mm"), 
                              1:n - 0.5, cn, 
                              cex = 0.75, 
                              adj = c(0, 0.5), 
                              facing = "bending.inside",
                              niceFacing = TRUE,
                              font = 2)},
              bg.border = NA)
  
  
  circos.heatmap(corr_matrix[,1 , drop = FALSE], 
               col = col_scale, 
               track.height = .1)
  circos.track(track.index = 4,
               panel.fun = function(x, y) {
                  cn = colnames(corr_matrix)[1]
                  n = length(cn)
                  circos.text(rep(CELL_META$cell.xlim[2], n) + convert_x(1, "mm"), 
                              1:n - 0.5, cn, 
                              cex = 0.75, 
                              adj = c(0, 0.5), 
                              facing = "bending.inside",
                              niceFacing = TRUE,
                              font = 2)},
              bg.border = NA)
  
  circos.clear()

upViewport()

lgd_color <- 
  Legend(at = seq(-1, 1, .2), 
         col_fun = col_scale, 
         title_position = "topcenter", 
         title = "r",
         border = "black")
draw(lgd_color,  x = circle_size, just = "left")

Note: 1 point is out of plotting region in sector 'group', track '2'.

Note: 1 point is out of plotting region in sector 'group', track '2'.

Note: 12 points are out of plotting region in sector 'group', track
'2'.

Note: 1 point is out of plotting region in sector 'group', track '3'.

Note: 1 point is out of plotting region in sector 'group', track '3'.

Note: 8 points are out of plotting region in sector 'group', track '3'.

Note: 1 point is out of plotting region in sector 'group', track '4'.

Note: 1 point is out of plotting region in sector 'group', track '4'.

Note: 5 points are out of plotting region in sector 'group', track '4'.

Nullification

list_wo_circ <- 
  paste0(c(glue("s({names(test_data)[-c(10, 34, 42, 45:47)]})"), 
           names(test_data)[42]),
        collapse = " + ")

student_nocirc_all <- 
  gam(as.formula(paste0("sofa_nocirc ~ s(urine24, k = 150) + s(crea, k = 125) + ", list_wo_circ)), 
      data = test_data)

student_circ_all <- 
  gam(as.formula(paste0("sofa_circ ~ s(urine24, k = 150) + s(crea, k = 125) + ", list_wo_circ)), 
      data = test_data)

list_wo_circ <- 
  paste0(c(glue("s({names(test_data)[-c(10, 34, 42, 45:47)]})"), 
           names(test_data)[42]),
        collapse = " + ")

student_nocirc_noBili <- 
  gam(as.formula(paste0("sofa_nocirc ~ ", list_wo_circ)), 
      data = test_data)

student_circ_noBili <- 
  gam(as.formula(paste0("sofa_circ ~ ",list_wo_circ)), 
      data = test_data)

TEACHER: circular, STUDENT: all features

summary(student_circ_all)

Family: gaussian 
Link function: identity 

Formula:
sofa_circ ~ s(urine24, k = 150) + s(crea, k = 125) + s(amv) + 
    s(heartrate) + s(leuko) + s(temp) + s(pco2) + s(rr) + s(artph) + 
    s(bili) + s(bun) + s(dbp) + s(fio2) + s(mbp) + s(pao2) + 
    s(sbp) + s(thrombo) + s(lactate) + s(bicarb) + s(crp) + s(haemo) + 
    s(hzv) + s(lympho) + s(natrium) + s(panklip) + s(pct) + s(quinr) + 
    s(saO2) + s(glucose) + s(baseex) + s(chlorid) + s(calcium) + 
    s(kalium) + s(szvo2) + s(netto) + s(alat) + s(asat) + s(svol) + 
    s(svri) + s(adrenaline) + s(dobutamine) + s(noradrenaline) + 
    s(horovitz) + dopamine

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 1.083611   0.001389 780.369  < 2e-16 ***
dopamine    0.120183   0.040510   2.967  0.00301 ** 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
                     edf  Ref.df        F p-value    
s(urine24)       121.389 136.408 1497.337  <2e-16 ***
s(crea)          123.897 123.999  508.270  <2e-16 ***
s(amv)             6.969   8.037    9.911  <2e-16 ***
s(heartrate)       8.915   8.997   25.192  <2e-16 ***
s(leuko)           8.994   9.000  102.721  <2e-16 ***
s(temp)            8.814   8.987   99.944  <2e-16 ***
s(pco2)            8.552   8.933   41.482  <2e-16 ***
s(rr)              8.618   8.932   17.297  <2e-16 ***
s(artph)           8.679   8.961   36.090  <2e-16 ***
s(bili)            8.991   9.000  216.846  <2e-16 ***
s(bun)             8.746   8.979  319.919  <2e-16 ***
s(dbp)             9.000   9.000   92.228  <2e-16 ***
s(fio2)            8.803   8.987   23.920  <2e-16 ***
s(mbp)             7.816   8.614   30.081  <2e-16 ***
s(pao2)            8.631   8.958   11.942  <2e-16 ***
s(sbp)             8.488   8.913    6.589  <2e-16 ***
s(thrombo)         8.992   9.000   95.454  <2e-16 ***
s(lactate)         7.225   8.209  144.790  <2e-16 ***
s(bicarb)          8.647   8.938   16.788  <2e-16 ***
s(crp)             8.935   8.998   62.565  <2e-16 ***
s(haemo)           8.818   8.988   36.431  <2e-16 ***
s(hzv)             6.701   7.635   47.006  <2e-16 ***
s(lympho)          8.818   8.987   36.340  <2e-16 ***
s(natrium)         8.361   8.826   69.806  <2e-16 ***
s(panklip)         8.904   8.996   66.658  <2e-16 ***
s(pct)             8.995   9.000  309.859  <2e-16 ***
s(quinr)           8.871   8.992  172.790  <2e-16 ***
s(saO2)            8.741   8.974   14.744  <2e-16 ***
s(glucose)         8.589   8.940   13.183  <2e-16 ***
s(baseex)          9.000   9.000   28.168  <2e-16 ***
s(chlorid)         8.983   9.000  102.063  <2e-16 ***
s(calcium)         8.801   8.985   64.826  <2e-16 ***
s(kalium)          8.843   8.989   30.444  <2e-16 ***
s(szvo2)           8.628   8.950   60.681  <2e-16 ***
s(netto)           4.620   5.777    8.274  <2e-16 ***
s(alat)            8.991   9.000   67.694  <2e-16 ***
s(asat)            8.867   8.993   26.943  <2e-16 ***
s(svol)            8.683   8.954   36.059  <2e-16 ***
s(svri)            7.990   8.675   17.733  <2e-16 ***
s(adrenaline)      4.264   4.462    0.156   0.946    
s(dobutamine)      8.990   9.000   53.987  <2e-16 ***
s(noradrenaline)   4.264   4.462    0.156   0.946    
s(horovitz)        8.955   8.998   20.935  <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Rank: 643/644
R-sq.(adj) =   0.93   Deviance explained =   93%
GCV = 0.15664  Scale est. = 0.1555    n = 80671

TEACHER: circular, STUDENT: no bili

summary(student_circ_noBili)

Family: gaussian 
Link function: identity 

Formula:
sofa_circ ~ s(amv) + s(heartrate) + s(leuko) + s(temp) + s(pco2) + 
    s(rr) + s(artph) + s(bili) + s(bun) + s(dbp) + s(fio2) + 
    s(mbp) + s(pao2) + s(sbp) + s(thrombo) + s(lactate) + s(bicarb) + 
    s(crp) + s(haemo) + s(hzv) + s(lympho) + s(natrium) + s(panklip) + 
    s(pct) + s(quinr) + s(saO2) + s(glucose) + s(baseex) + s(chlorid) + 
    s(calcium) + s(kalium) + s(szvo2) + s(netto) + s(alat) + 
    s(asat) + s(svol) + s(svri) + s(adrenaline) + s(dobutamine) + 
    s(noradrenaline) + s(horovitz) + dopamine

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 1.083600   0.003262 332.151   <2e-16 ***
dopamine    0.137979   0.093929   1.469    0.142    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
                   edf Ref.df       F p-value    
s(amv)           8.729  8.971  23.599  <2e-16 ***
s(heartrate)     7.841  8.564  53.538  <2e-16 ***
s(leuko)         8.972  9.000 119.455  <2e-16 ***
s(temp)          8.952  8.999 205.966  <2e-16 ***
s(pco2)          7.943  8.671  19.966  <2e-16 ***
s(rr)            7.945  8.637 100.212  <2e-16 ***
s(artph)         8.779  8.977  57.518  <2e-16 ***
s(bili)          8.989  9.000 198.342  <2e-16 ***
s(bun)           8.161  8.804 816.376  <2e-16 ***
s(dbp)           8.943  8.998  54.852  <2e-16 ***
s(fio2)          8.945  8.999  26.779  <2e-16 ***
s(mbp)           7.913  8.664  17.054  <2e-16 ***
s(pao2)          8.434  8.901  10.882  <2e-16 ***
s(sbp)           8.381  8.877  21.086  <2e-16 ***
s(thrombo)       8.974  9.000  82.331  <2e-16 ***
s(lactate)       8.779  8.983  40.398  <2e-16 ***
s(bicarb)        7.564  8.335  15.875  <2e-16 ***
s(crp)           8.931  8.998  54.837  <2e-16 ***
s(haemo)         8.873  8.993  79.931  <2e-16 ***
s(hzv)           8.314  8.825   9.737  <2e-16 ***
s(lympho)        8.839  8.990  94.743  <2e-16 ***
s(natrium)       9.000  9.000 121.101  <2e-16 ***
s(panklip)       8.894  8.996  61.948  <2e-16 ***
s(pct)           8.847  8.992  73.306  <2e-16 ***
s(quinr)         8.845  8.989 192.689  <2e-16 ***
s(saO2)          8.770  8.979  14.367  <2e-16 ***
s(glucose)       7.460  8.338  78.953  <2e-16 ***
s(baseex)        8.926  8.981   9.636  <2e-16 ***
s(chlorid)       8.752  8.971 311.549  <2e-16 ***
s(calcium)       8.958  8.999 268.796  <2e-16 ***
s(kalium)        7.576  8.378  78.095  <2e-16 ***
s(szvo2)         8.958  8.999  63.294  <2e-16 ***
s(netto)         8.971  9.000  67.850  <2e-16 ***
s(alat)          8.530  8.931  66.171  <2e-16 ***
s(asat)          8.910  8.997 114.951  <2e-16 ***
s(svol)          8.916  8.996  45.572  <2e-16 ***
s(svri)          8.894  8.995  22.108  <2e-16 ***
s(adrenaline)    4.208  4.442   0.427   0.851    
s(dobutamine)    8.576  8.939  34.478  <2e-16 ***
s(noradrenaline) 4.208  4.442   0.427   0.851    
s(horovitz)      8.801  8.980   9.793  <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Rank: 370/371
R-sq.(adj) =  0.613   Deviance explained = 61.4%
GCV = 0.86201  Scale est. = 0.85831   n = 80671

TEACHER: non-circular, STUDENT: all features

summary(student_nocirc_all)

Family: gaussian 
Link function: identity 

Formula:
sofa_nocirc ~ s(urine24, k = 150) + s(crea, k = 125) + s(amv) + 
    s(heartrate) + s(leuko) + s(temp) + s(pco2) + s(rr) + s(artph) + 
    s(bili) + s(bun) + s(dbp) + s(fio2) + s(mbp) + s(pao2) + 
    s(sbp) + s(thrombo) + s(lactate) + s(bicarb) + s(crp) + s(haemo) + 
    s(hzv) + s(lympho) + s(natrium) + s(panklip) + s(pct) + s(quinr) + 
    s(saO2) + s(glucose) + s(baseex) + s(chlorid) + s(calcium) + 
    s(kalium) + s(szvo2) + s(netto) + s(alat) + s(asat) + s(svol) + 
    s(svri) + s(adrenaline) + s(dobutamine) + s(noradrenaline) + 
    s(horovitz) + dopamine

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 1.136660   0.001222 930.121   <2e-16 ***
dopamine    0.030698   0.035666   0.861    0.389    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
                     edf  Ref.df        F p-value    
s(urine24)       126.248 139.812 2089.634  <2e-16 ***
s(crea)          123.877 123.999  601.264  <2e-16 ***
s(amv)             7.737   8.569   11.115  <2e-16 ***
s(heartrate)       8.741   8.971   22.112  <2e-16 ***
s(leuko)           8.992   9.000   75.686  <2e-16 ***
s(temp)            8.749   8.976   13.143  <2e-16 ***
s(pco2)            8.594   8.941   40.944  <2e-16 ***
s(rr)              8.120   8.720   29.506  <2e-16 ***
s(artph)           8.687   8.959   28.334  <2e-16 ***
s(bili)            9.000   9.000   79.759  <2e-16 ***
s(bun)             8.693   8.971  289.187  <2e-16 ***
s(dbp)             8.987   9.000   65.109  <2e-16 ***
s(fio2)            8.818   8.989   31.968  <2e-16 ***
s(mbp)             8.553   8.933   20.778  <2e-16 ***
s(pao2)            8.706   8.973   17.220  <2e-16 ***
s(sbp)             8.616   8.951    9.641  <2e-16 ***
s(thrombo)         8.890   8.996   60.536  <2e-16 ***
s(lactate)         7.923   8.672   67.496  <2e-16 ***
s(bicarb)          6.653   7.697   41.026  <2e-16 ***
s(crp)             8.906   8.997   34.006  <2e-16 ***
s(haemo)           8.459   8.900   72.280  <2e-16 ***
s(hzv)             8.790   8.981   18.996  <2e-16 ***
s(lympho)          8.819   8.987   61.013  <2e-16 ***
s(natrium)         8.983   9.000   44.996  <2e-16 ***
s(panklip)         8.854   8.992  440.856  <2e-16 ***
s(pct)             8.994   9.000  212.210  <2e-16 ***
s(quinr)           8.884   8.993   70.667  <2e-16 ***
s(saO2)            8.854   8.991   34.379  <2e-16 ***
s(glucose)         8.765   8.980   21.289  <2e-16 ***
s(baseex)          9.000   9.000   45.697  <2e-16 ***
s(chlorid)         8.892   8.994   40.391  <2e-16 ***
s(calcium)         8.229   8.804   91.472  <2e-16 ***
s(kalium)          8.988   9.000   28.672  <2e-16 ***
s(szvo2)           8.733   8.974   96.264  <2e-16 ***
s(netto)           8.295   8.863   49.580  <2e-16 ***
s(alat)            8.545   8.935   92.962  <2e-16 ***
s(asat)            8.859   8.992   44.692  <2e-16 ***
s(svol)            8.976   9.000   65.472  <2e-16 ***
s(svri)            8.983   9.000   25.457  <2e-16 ***
s(adrenaline)      4.197   4.439    0.123   0.933    
s(dobutamine)      8.825   8.989  104.676  <2e-16 ***
s(noradrenaline)   4.197   4.439    0.123   0.933    
s(horovitz)        8.987   9.000   34.222  <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Rank: 643/644
R-sq.(adj) =  0.947   Deviance explained = 94.8%
GCV = 0.12134  Scale est. = 0.12044   n = 80671

TEACHER: non-circular, STUDENT: no bili

summary(student_nocirc_noBili)

Family: gaussian 
Link function: identity 

Formula:
sofa_nocirc ~ s(amv) + s(heartrate) + s(leuko) + s(temp) + s(pco2) + 
    s(rr) + s(artph) + s(bili) + s(bun) + s(dbp) + s(fio2) + 
    s(mbp) + s(pao2) + s(sbp) + s(thrombo) + s(lactate) + s(bicarb) + 
    s(crp) + s(haemo) + s(hzv) + s(lympho) + s(natrium) + s(panklip) + 
    s(pct) + s(quinr) + s(saO2) + s(glucose) + s(baseex) + s(chlorid) + 
    s(calcium) + s(kalium) + s(szvo2) + s(netto) + s(alat) + 
    s(asat) + s(svol) + s(svri) + s(adrenaline) + s(dobutamine) + 
    s(noradrenaline) + s(horovitz) + dopamine

Parametric coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) 1.136642   0.003268 347.811   <2e-16 ***
dopamine    0.059909   0.094072   0.637    0.524    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Approximate significance of smooth terms:
                   edf Ref.df       F p-value    
s(amv)           8.655  8.958  26.275  <2e-16 ***
s(heartrate)     6.125  7.256  54.553  <2e-16 ***
s(leuko)         8.970  9.000  81.739  <2e-16 ***
s(temp)          8.949  8.999 147.489  <2e-16 ***
s(pco2)          7.676  8.516  12.607  <2e-16 ***
s(rr)            8.226  8.785 110.209  <2e-16 ***
s(artph)         8.894  8.993  53.273  <2e-16 ***
s(bili)          8.998  9.000 125.727  <2e-16 ***
s(bun)           8.310  8.863 817.307  <2e-16 ***
s(dbp)           8.937  8.997  46.533  <2e-16 ***
s(fio2)          8.908  8.997  19.330  <2e-16 ***
s(mbp)           7.968  8.691  18.577  <2e-16 ***
s(pao2)          8.391  8.884   7.398  <2e-16 ***
s(sbp)           8.395  8.882  26.801  <2e-16 ***
s(thrombo)       8.934  8.998  71.730  <2e-16 ***
s(lactate)       8.787  8.984  36.897  <2e-16 ***
s(bicarb)        8.311  8.760  23.120  <2e-16 ***
s(crp)           8.907  8.997  44.540  <2e-16 ***
s(haemo)         8.824  8.987  74.802  <2e-16 ***
s(hzv)           5.388  6.425  26.743  <2e-16 ***
s(lympho)        8.787  8.983 118.889  <2e-16 ***
s(natrium)       8.998  9.000 130.535  <2e-16 ***
s(panklip)       8.827  8.989  80.814  <2e-16 ***
s(pct)           8.858  8.993  76.317  <2e-16 ***
s(quinr)         8.932  8.998 144.125  <2e-16 ***
s(saO2)          8.695  8.964   7.049  <2e-16 ***
s(glucose)       7.299  8.217  82.998  <2e-16 ***
s(baseex)        8.593  8.873  12.448  <2e-16 ***
s(chlorid)       8.720  8.963 316.544  <2e-16 ***
s(calcium)       8.919  8.997 290.277  <2e-16 ***
s(kalium)        7.150  8.040 117.336  <2e-16 ***
s(szvo2)         8.965  8.999  77.284  <2e-16 ***
s(netto)         8.973  9.000  58.513  <2e-16 ***
s(alat)          8.930  8.998  71.220  <2e-16 ***
s(asat)          8.844  8.991  91.756  <2e-16 ***
s(svol)          8.951  8.998  49.970  <2e-16 ***
s(svri)          8.725  8.970  21.731  <2e-16 ***
s(adrenaline)    4.146  4.417   0.575   0.773    
s(dobutamine)    7.972  8.673  50.203  <2e-16 ***
s(noradrenaline) 4.146  4.417   0.575   0.773    
s(horovitz)      8.644  8.944  11.190  <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Rank: 370/371
R-sq.(adj) =  0.622   Deviance explained = 62.3%
GCV = 0.86492  Scale est. = 0.86127   n = 80671

corr_matrix[order(abs(corr_matrix[,"circular"]), decreasing=TRUE),][1:5,]

A matrix: 5 × 3 of type dbl

	label	circular	non-circular
crea	0.6273276	0.6365420	0.6287862
urine24	-0.5038256	-0.4884345	-0.5146991
artph	-0.3339940	-0.3614628	-0.3553907
bun	0.3378248	0.3377452	0.3371343
quinr	0.2957017	0.3222822	0.2899301

z <- plot(student_circ_all, select = 2)
shift_crea   <- abs(mean(z[[2]]$fit[between(z[[2]]$x, 0, 1)]))
shift_urine24 <- coef(student_circ_all)[1] - shift_crea

dev_nocirc   <- round(summary(student_circ_noBili)$dev.expl * 100)
dev_withcirc <- round(summary(student_circ_all)$dev.expl * 100)

par(mfrow = c(5,2),  mar=c(2,3,2,3), oma = c(1,1,1,1))
 #first row: urine24
   plot(x = 0:1000, y = seq(0,4, length.out=1001),
       xlim   = c(0,1000), 
       ylim   = c(0, 4), 
       main   = glue("Student without Creatinine and Urine Output (D\U00B2={dev_nocirc}%)"),
       ylab   = "", 
       xlab   = "",
       type   = "n")
  text(labels = "not included in student model",
       x      = mean(c(0,1000)),
       y      = mean(c(0,4)),
       adj    = .5)
  
  plot(student_circ_all, 
     se     = FALSE, 
     xlim   = c(0,1000), 
     ylim   = c(0,4), 
     main   = glue("Student with all Features (D\U00B2={dev_withcirc}%)"),
     ylab   = "", 
     xlab   = "",
     select = 1,
     shift  = shift_urine24)
   rug(unique(test_data$urine24))
  mtext("urine ouput", side = 4, line = 1)
  
  #second row: crea
   plot(x = seq(.15, 5.2, by = .01), y = seq(0, 4, length.out=506),
       xlim   = c(.15, 5.2), 
       ylim   = c(0, 4), 
       ylab   = "", 
       xlab   = "",
       type   = "n")
  text(labels = "not included in student model",
       x      = mean(c(.15, 5.2)),
       y      = mean(c(0,4)),
       adj    = .5)
  
  plot(student_circ_all, 
       ylab   = "",
       xlab   = "",
       xlim   = c(.15, 5.2), 
       ylim   = c(0, 4), 
       select = 2,
       shift = shift_crea )
    abline(h=0, lty = 3)
      rug(unique(test_data$crea))
  mtext("creatinine", side = 4, line = 1)
  
  #third row: calcium
  plot(student_circ_noBili, 
       se     = FALSE, 
       #rug    = TRUE, 
       ylab   = "",
       xlab   = "",
       ylim   = c(-1,4), 
       xlim   = c(.5, 1.6), 
       select = 30)
    abline(h=0, lty = 3)
    rug(unique(test_data$calcium))
  
  plot(student_circ_all, 
       se     = FALSE, 
       #rug    = TRUE, 
       ylab   = "",
       xlab   = "",
       ylim   = c(-1,4), 
       xlim   = c(.5, 1.6), 
       select = 32)
    abline(h=0, lty = 3)
    rug(unique(test_data$calcium))
  mtext("calcium", side = 4, line = 1)
  
  #forth row: bun
  plot(student_circ_noBili, 
       se     = FALSE, 
       #rug    = TRUE, 
       #ylab   = "", 
       xlab   = "",
       ylim   = c(-1,4),
       xlim   = c(10, 180), 
       select = 9)
    rug(unique(test_data$bun))
    abline(h=0, lty = 3)
  
  plot(student_circ_all, 
       se     = FALSE, 
       #rug    = TRUE, 
       ylab   = "", 
       #xlab   = "",
       ylim   = c(-1,4),
       xlim   = c(10, 180),
       select = 11)
    abline(h=0, lty = 3)
    rug(unique(test_data$bun))
  mtext("bun", side = 4, line = 1)
  
  #fifth row: chlorid
  plot(student_circ_noBili, 
       se     = FALSE, 
       #rug    = TRUE, 
       ylab   = "",
       #xlab   = "",
       ylim   = c(-1,4),
       xlim   = c(95,140),
       select = 29)
    abline(h=0, lty = 3)
    rug(unique(test_data$chlorid))
  
  plot(student_circ_all, 
       se     = FALSE, 
       #rug    = TRUE, 
       ylab   = "",
       xlab   = "",
       ylim   = c(-1,4),
       xlim   = c(95,140),
       select = 31)
    abline(h=0, lty = 3)
    rug(unique(test_data$chlorid))
  mtext("chlorid", side = 4, line = 1)
par(mfrow = c(1,1), mar = c(5,4,4,2))

Warning message in rug(unique(test_data$urine24)):
“some values will be clipped”
Warning message in rug(unique(test_data$crea)):
“some values will be clipped”
Warning message in rug(unique(test_data$calcium)):
“some values will be clipped”
Warning message in rug(unique(test_data$calcium)):
“some values will be clipped”
Warning message in rug(unique(test_data$bun)):
“some values will be clipped”
Warning message in rug(unique(test_data$bun)):
“some values will be clipped”
Warning message in rug(unique(test_data$chlorid)):
“some values will be clipped”
Warning message in rug(unique(test_data$chlorid)):
“some values will be clipped”