Mucosal microbiota alterations in primary sclerosis cholangitis persist after liver transplantation and are associated with clinical features independently of geography

ML models - overfitting check

---

1 Introduction

To prevent potential training errors and overfitting, a validity check was performed for all models by randomly shuffling the sample labels. The performance of these models should not reach significantly high values.

Importing libraries and custom functions built for this analysis

source("custom_functions.R")

2 Data Import

Importing ASV, taxa and metadata tables for both Czech and Norway samples.

Czech

path = "../../data/analysis_ready_data/ikem/"
asv_tab_ikem <- as.data.frame(fread(file.path(path,"asv_table_ikem.csv"),
                                    check.names = FALSE))
taxa_tab_ikem <- as.data.frame(fread(file.path(path,"taxa_table_ikem.csv"),
                                     check.names = FALSE))
metadata_ikem <- as.data.frame(fread(file.path(path,"metadata_ikem.csv"),
                                     check.names = FALSE))

Norway

path = "../../data/analysis_ready_data/norway/"
asv_tab_norway <- as.data.frame(fread(file.path(path,"asv_table_norway.csv"),
                                    check.names = FALSE))
taxa_tab_norway <- as.data.frame(fread(file.path(path,"taxa_table_norway.csv"),
                                    check.names = FALSE))
metadata_norway <- as.data.frame(fread(file.path(path,"metadata_norway.csv"),
                                    check.names = FALSE))

3 Merging

Terminal ileum

ileum_data <- merging_data(asv_tab_1=asv_tab_ikem,
                           asv_tab_2=asv_tab_norway,
                           taxa_tab_1=taxa_tab_ikem,
                           taxa_tab_2=taxa_tab_norway,
                           metadata_1=metadata_ikem,
                           metadata_2=metadata_norway,
                           segment="TI",Q="Q1")

Removing 1498 ASV(s)

Removing 1834 ASV(s)

Merging at ASV level

Finding inconsistencies in taxonomy, trying to keep the ones that have better taxonomy assignment

ileum_asv_tab <- ileum_data[[1]]
ileum_taxa_tab <- ileum_data[[2]]
ileum_metadata <- ileum_data[[3]]

Colon

colon_data <- merging_data(asv_tab_1=asv_tab_ikem,
                           asv_tab_2=asv_tab_norway,
                           taxa_tab_1=taxa_tab_ikem,
                           taxa_tab_2=taxa_tab_norway,
                           metadata_1=metadata_ikem,
                           metadata_2=metadata_norway,
                           segment="colon",Q="Q1")

Removing 739 ASV(s)

Removing 266 ASV(s)

Merging at ASV level

Finding inconsistencies in taxonomy, trying to keep the ones that have better taxonomy assignment

colon_asv_tab <- colon_data[[1]]
colon_taxa_tab <- colon_data[[2]]
colon_metadata <- colon_data[[3]]

4 Data Analysis - Terminal ileum

segment="terminal_ileum"

4.1 Machine learning

path = "../results/Q1/models_overfitting_check"

4.1.1 ENET

model="enet"

4.1.1.1 Genus level

level="genus"

Aggregate taxa

genus_data <- aggregate_taxa(ileum_asv_tab,
                             ileum_taxa_tab,
                             taxonomic_level = level)

ileum_genus_tab <- genus_data[[1]]
ileum_genus_taxa_tab <- genus_data[[2]]

4.1.1.1.1 pre_ltx vs healthy

group <- c("pre_ltx","healthy")
comparison_name <- paste0(group[1], " vs ",group[2])

model_name <- paste(comparison_name,level,segment)

# prepare the data
filt_ileum_uni_data <- binomial_prep(ileum_genus_tab,
                                     ileum_genus_taxa_tab,
                                     ileum_metadata,
                                     group, 
                                     usage="ml_clr")

Removing 84 ASV(s)

Removing 10 ASV(s)

# fit the model
enet_model <- glmnet_binomial(filt_ileum_uni_data,
                              sample_method = "atypboot",
                              outcome="Group", 
                              N=500,
                              reuse=TRUE,
                              file=model_name,
                              Q="Q1",
                              overfitting_check = TRUE)

# ROC
roc_c <- roc_curve(enet_model, group)

models_summ <- list()
models_cm <- list()
roc_cs <- list()
betas <- list()
# save the results
models_summ[[model_name]] <- enet_model$model_summary
models_cm[[model_name]] <- enet_model$conf_matrices$original
roc_cs[[model_name]] <- enet_model$kfold_rocobjs
betas[[model_name]] <- as.matrix(enet_model$betas)

# see the results
enet_model$model_summary %>% t()

                                       [,1]
alpha                             0.0000000
lambda                          122.5503879
auc                               0.7645793
auc_czech                         0.6958874
auc_no                            0.8234234
auc_optimism_corrected            0.4484626
auc_optimism_corrected_CIL        0.3242591
auc_optimism_corrected_CIU        0.5949898
accuracy                          0.5104895
accuracy_czech                          NaN
accuracy_no                       0.5487805
accuracy_optimism_corrected       0.4543236
accuracy_optimism_corrected_CIL   0.3396226
accuracy_optimism_corrected_CIU   0.5742809

enet_model$conf_matrices

$original
   0  
0 73 0
1 70 0

$czech
   0  
0 28 0
1 33 0

$no
   0  
0 45 0
1 37 0

enet_model$plot

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

roc_c

4.1.1.1.2 pre_ltx vs post_ltx

group <- c("pre_ltx","post_ltx")
comparison_name <- paste0(group[1], " vs ",group[2])

model_name <- paste(comparison_name,level,segment)

# prepare the data
filt_ileum_uni_data <- binomial_prep(ileum_genus_tab,
                                     ileum_genus_taxa_tab,
                                     ileum_metadata,
                                     group, 
                                     usage="ml_clr")

Removing 46 ASV(s)

Removing 6 ASV(s)

# fit the model
enet_model <- glmnet_binomial(filt_ileum_uni_data,
                              sample_method = "atypboot",
                              outcome="Group", N=500,
                              reuse=TRUE,
                              file=model_name,
                              Q="Q1",
                              overfitting_check = TRUE)

# ROC
roc_c <- roc_curve(enet_model, group)

# save the results
models_summ[[model_name]] <- enet_model$model_summary
models_cm[[model_name]] <- enet_model$conf_matrices$original
roc_cs[[model_name]] <- enet_model$kfold_rocobjs
betas[[model_name]] <- as.matrix(enet_model$betas)

# see the results
enet_model$model_summary %>% t()

                                      [,1]
alpha                           1.00000000
lambda                          0.09461542
auc                             0.50000000
auc_czech                       0.50000000
auc_no                          0.50000000
auc_optimism_corrected          0.44402219
auc_optimism_corrected_CIL      0.32593459
auc_optimism_corrected_CIU      0.56059550
accuracy                        0.66346154
accuracy_czech                         NaN
accuracy_no                     0.60975610
accuracy_optimism_corrected     0.58967912
accuracy_optimism_corrected_CIL 0.47400472
accuracy_optimism_corrected_CIU 0.69620253

enet_model$conf_matrices

$original
    0  
0 138 0
1  70 0

$czech
   0  
0 88 0
1 38 0

$no
   0  
0 50 0
1 32 0

enet_model$plot

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

roc_c

4.1.1.1.3 post_ltx vs healthy

group <- c("post_ltx","healthy")
comparison_name <- paste0(group[1], " vs ",group[2])

model_name <- paste(comparison_name,level,segment)

# prepare the data
filt_ileum_uni_data <- binomial_prep(ileum_genus_tab,
                                     ileum_genus_taxa_tab,
                                     ileum_metadata,group, 
                                     usage="ml_clr")

Removing 45 ASV(s)

Removing 2 ASV(s)

# fit the model
enet_model <- glmnet_binomial(filt_ileum_uni_data,
                              sample_method = "atypboot",
                              outcome="Group", 
                              N=500,
                              reuse=TRUE,
                              file=model_name,
                              Q="Q1",
                              overfitting_check = TRUE)

# ROC
roc_c <- roc_curve(enet_model, group)

# save the results
models_summ[[model_name]] <- enet_model$model_summary
models_cm[[model_name]] <- enet_model$conf_matrices$original
roc_cs[[model_name]] <- enet_model$kfold_rocobjs
betas[[model_name]] <- as.matrix(enet_model$betas)

# see the results
enet_model$model_summary %>% t()

                                      [,1]
alpha                            0.0000000
lambda                          81.9195172
auc                              0.7156045
auc_czech                        0.6625641
auc_no                           0.7407121
auc_optimism_corrected           0.5175379
auc_optimism_corrected_CIL       0.3928180
auc_optimism_corrected_CIU       0.6223579
accuracy                         0.6540284
accuracy_czech                         NaN
accuracy_no                      0.5277778
accuracy_optimism_corrected      0.6017558
accuracy_optimism_corrected_CIL  0.5136918
accuracy_optimism_corrected_CIU  0.7006818

enet_model$conf_matrices

$original
    1  
0  73 0
1 138 0

$czech
    1  
0  39 0
1 100 0

$no
   1  
0 34 0
1 38 0

enet_model$plot

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

roc_c

4.1.1.2 Saving results

models_summ_df_ileum <- do.call(rbind, 
  models_summ[grep(segment,names(models_summ),value = TRUE)])

write.csv(models_summ_df_ileum,file.path(path,paste0("elastic_net_",segment,".csv")))

4.1.2 Supplementary models

supplements_models <- list()

4.1.2.1 CLR-transformed data

4.1.2.1.1 kNN

model="knn"

4.1.2.1.1.1 Genus level

level="genus"

Aggregate taxa

genus_data <- aggregate_taxa(ileum_asv_tab,
                             ileum_taxa_tab,
                             taxonomic_level = level)

ileum_genus_tab <- genus_data[[1]]
ileum_genus_taxa_tab <- genus_data[[2]]

pre_ltx vs healthy

group <- c("pre_ltx","healthy")
comparison_name <- paste0(group[1], " vs ",group[2])

model_name <- paste(comparison_name,level,segment)

# prepare the data
filt_ileum_uni_data <- binomial_prep(ileum_genus_tab,
                                     ileum_genus_taxa_tab,
                                     ileum_metadata,
                                     group, 
                                     usage="ml_clr")

Removing 84 ASV(s)

Removing 10 ASV(s)

# fit the model
knn_model <- knn_binomial(filt_ileum_uni_data,
                              sample_method = "atypboot",
                              outcome="Group", 
                              N=500,
                              reuse=TRUE,
                              file=model_name,
                              Q="Q1",
                              overfitting_check = TRUE)

# ROC curve
roc_c <- roc_curve(knn_model, group)

# save the results
supplements_models[["models_summ"]][["knn_model"]][[model_name]] <- knn_model$model_summary
supplements_models[["roc_cs"]][["knn_model"]][[model_name]] <- knn_model$kfold_rocobjs


# see the results
knn_model$model_summary %>% t()

                                      [,1]
k                               13.0000000
auc                              0.6862035
auc_optimism_corrected           0.5110905
auc_optimism_corrected_CIL       0.3610817
auc_optimism_corrected_CIU       0.6615794
accuracy                         0.5874126
accuracy_optimism_corrected      0.4977140
accuracy_optimism_corrected_CIL  0.3796379
accuracy_optimism_corrected_CIU  0.6226415

roc_c

pre_ltx vs post_ltx

group <- c("pre_ltx","post_ltx")
comparison_name <- paste0(group[1], " vs ",group[2])

model_name <- paste(comparison_name,level,segment)

# prepare the data
filt_ileum_uni_data <- binomial_prep(ileum_genus_tab,
                                     ileum_genus_taxa_tab,
                                     ileum_metadata,
                                     group, 
                                     usage="ml_clr")

Removing 46 ASV(s)

Removing 6 ASV(s)

# fit the model
knn_model <- knn_binomial(filt_ileum_uni_data,
                              sample_method = "atypboot",
                              outcome="Group", 
                              N=500,
                              reuse=TRUE,
                              file=model_name,
                              Q="Q1",
                              overfitting_check = TRUE)

# ROC curve
roc_c <- roc_curve(knn_model, group)

# save the results
supplements_models[["models_summ"]][["knn_model"]][[model_name]] <- knn_model$model_summary
supplements_models[["roc_cs"]][["knn_model"]][[model_name]] <- knn_model$kfold_rocobjs


# see the results
knn_model$model_summary %>% t()

                                      [,1]
k                               27.0000000
auc                              0.5889752
auc_optimism_corrected           0.4685397
auc_optimism_corrected_CIL       0.3548958
auc_optimism_corrected_CIU       0.5926814
accuracy                         0.6634615
accuracy_optimism_corrected      0.6286877
accuracy_optimism_corrected_CIL  0.5218434
accuracy_optimism_corrected_CIU  0.7080288

roc_c

post_ltx vs healthy

group <- c("post_ltx","healthy")
comparison_name <- paste0(group[1], " vs ",group[2])

model_name <- paste(comparison_name,level,segment)

# prepare the data
filt_ileum_uni_data <- binomial_prep(ileum_genus_tab,
                                     ileum_genus_taxa_tab,
                                     ileum_metadata,
                                     group, 
                                     usage="ml_clr")

Removing 45 ASV(s)

Removing 2 ASV(s)

# fit the model
knn_model <- knn_binomial(filt_ileum_uni_data,
                              sample_method = "atypboot",
                              outcome="Group", 
                              N=500,
                              reuse=TRUE,
                              file=model_name,
                              Q="Q1",
                              overfitting_check = TRUE)

# ROC curve
roc_c <- roc_curve(knn_model, group)

# save the results
supplements_models[["models_summ"]][["knn_model"]][[model_name]] <- knn_model$model_summary
supplements_models[["roc_cs"]][["knn_model"]][[model_name]] <- knn_model$kfold_rocobjs


# see the results
knn_model$model_summary %>% t()

                                      [,1]
k                               21.0000000
auc                              0.5996625
auc_optimism_corrected           0.4792983
auc_optimism_corrected_CIL       0.3537347
auc_optimism_corrected_CIU       0.5968061
accuracy                         0.6872038
accuracy_optimism_corrected      0.5839843
accuracy_optimism_corrected_CIL  0.4670779
accuracy_optimism_corrected_CIU  0.6857902

roc_c

4.1.2.1.2 RF

model="rf"

4.1.2.1.2.1 Genus level

level="genus"

Aggregate taxa

genus_data <- aggregate_taxa(ileum_asv_tab,
                             ileum_taxa_tab,
                             taxonomic_level = level)

ileum_genus_tab <- genus_data[[1]]
ileum_genus_taxa_tab <- genus_data[[2]]

pre_ltx vs healthy

group <- c("pre_ltx","healthy")
comparison_name <- paste0(group[1], " vs ",group[2])

model_name <- paste(comparison_name,level,segment)

# prepare the data
filt_ileum_uni_data <- binomial_prep(ileum_genus_tab,
                                     ileum_genus_taxa_tab,
                                     ileum_metadata,
                                     group, 
                                     usage="ml_clr")

Removing 84 ASV(s)

Removing 10 ASV(s)

# fit the model
rf_model <- rf_binomial(filt_ileum_uni_data,
                              sample_method = "atypboot",
                              outcome="Group", 
                              N=500,
                              reuse=TRUE,
                              file=model_name,
                              Q="Q1",
                              overfitting_check = TRUE)

# ROC curve
roc_c <- roc_curve(rf_model, group)

# save the results
supplements_models[["models_summ"]][["rf_model"]][[model_name]] <- rf_model$model_summary
supplements_models[["roc_cs"]][["rf_model"]][[model_name]] <- rf_model$kfold_rocobjs

# see the results
rf_model$model_summary %>% t()

                                [,1]       
mtry                            "177"      
splitrule                       "gini"     
min.node.size                   "5"        
auc                             "1"        
auc_optimism_corrected          "0.4369043"
auc_optimism_corrected_CIL      "0.2979512"
auc_optimism_corrected_CIU      "0.5754973"
accuracy                        "1"        
accuracy_optimism_corrected     "0.4466668"
accuracy_optimism_corrected_CIL "0.3333333"
accuracy_optimism_corrected_CIU "0.5643573"

roc_c

pre_ltx vs post_ltx

group <- c("pre_ltx","post_ltx")
comparison_name <- paste0(group[1], " vs ",group[2])

model_name <- paste(comparison_name,level,segment)

# prepare the data
filt_ileum_uni_data <- binomial_prep(ileum_genus_tab,
                                     ileum_genus_taxa_tab,
                                     ileum_metadata,
                                     group, 
                                     usage="ml_clr")

Removing 46 ASV(s)

Removing 6 ASV(s)

# fit the model
rf_model <- rf_binomial(filt_ileum_uni_data,
                              sample_method = "atypboot",
                              outcome="Group", 
                              N=500,
                              reuse=TRUE,
                              file=model_name,
                              Q="Q1",
                              overfitting_check = TRUE)

# ROC curve
roc_c <- roc_curve(rf_model, group)

# save the results
supplements_models[["models_summ"]][["rf_model"]][[model_name]] <- rf_model$model_summary
supplements_models[["roc_cs"]][["rf_model"]][[model_name]] <- rf_model$kfold_rocobjs

# see the results
rf_model$model_summary %>% t()

                                [,1]       
mtry                            "25"       
splitrule                       "gini"     
min.node.size                   "5"        
auc                             "1"        
auc_optimism_corrected          "0.4166664"
auc_optimism_corrected_CIL      "0.3120122"
auc_optimism_corrected_CIU      "0.5166617"
accuracy                        "1"        
accuracy_optimism_corrected     "0.5963562"
accuracy_optimism_corrected_CIL "0.5152902"
accuracy_optimism_corrected_CIU "0.6762738"

roc_c

post_ltx vs healthy

group <- c("post_ltx","healthy")
comparison_name <- paste0(group[1], " vs ",group[2])

model_name <- paste(comparison_name,level,segment)

# prepare the data
filt_ileum_uni_data <- binomial_prep(ileum_genus_tab,
                                     ileum_genus_taxa_tab,
                                     ileum_metadata,
                                     group, 
                                     usage="ml_clr")

Removing 45 ASV(s)

Removing 2 ASV(s)

# fit the model
rf_model <- rf_binomial(filt_ileum_uni_data,
                              sample_method = "atypboot",
                              outcome="Group", 
                              N=500,
                              reuse=TRUE,
                              file=model_name,
                              Q="Q1",
                              overfitting_check = TRUE)

# ROC curve
roc_c <- roc_curve(rf_model, group)

# save the results
supplements_models[["models_summ"]][["rf_model"]][[model_name]] <- rf_model$model_summary
supplements_models[["roc_cs"]][["rf_model"]][[model_name]] <- rf_model$kfold_rocobjs

# see the results
rf_model$model_summary %>% t()

                                [,1]       
mtry                            "1"        
splitrule                       "gini"     
min.node.size                   "5"        
auc                             "1"        
auc_optimism_corrected          "0.4759634"
auc_optimism_corrected_CIL      "0.3675048"
auc_optimism_corrected_CIU      "0.5945952"
accuracy                        "1"        
accuracy_optimism_corrected     "0.595825" 
accuracy_optimism_corrected_CIL "0.4938978"
accuracy_optimism_corrected_CIU "0.6783764"

roc_c

4.1.2.1.3 GBoost

model="gb"

4.1.2.1.3.1 Genus level

level="genus"

Aggregate taxa

genus_data <- aggregate_taxa(ileum_asv_tab,
                             ileum_taxa_tab,
                             taxonomic_level = level)

ileum_genus_tab <- genus_data[[1]]
ileum_genus_taxa_tab <- genus_data[[2]]

pre_ltx vs healthy

group <- c("pre_ltx","healthy")
comparison_name <- paste0(group[1], " vs ",group[2])

model_name <- paste(comparison_name,level,segment)

# prepare the data
filt_ileum_uni_data <- binomial_prep(ileum_genus_tab,
                                     ileum_genus_taxa_tab,
                                     ileum_metadata,
                                     group, 
                                     usage="ml_clr")

Removing 84 ASV(s)

Removing 10 ASV(s)

# fit the model
gbm_model <- gbm_binomial(filt_ileum_uni_data,
                              sample_method = "atypboot",
                              outcome="Group", 
                              N=500,
                              reuse=TRUE,
                              file=model_name,
                              Q="Q1",
                              overfitting_check = TRUE)

# ROC curve
roc_c <- roc_curve(gbm_model, group)

# save the results
supplements_models[["models_summ"]][["gbm_model"]][[model_name]] <- gbm_model$model_summary
supplements_models[["roc_cs"]][["gbm_model"]][[model_name]] <- gbm_model$kfold_rocobjs

# see the results
gbm_model$model_summary %>% t()

                                       [,1]
n.trees                         500.0000000
interaction.depth                 1.0000000
shrinkage                         0.1000000
n.minobsinnode                   20.0000000
auc                               1.0000000
auc_optimism_corrected            0.4316582
auc_optimism_corrected_CIL        0.2936529
auc_optimism_corrected_CIU        0.5382335
accuracy                          1.0000000
accuracy_optimism_corrected       0.4507677
accuracy_optimism_corrected_CIL   0.3274216
accuracy_optimism_corrected_CIU   0.5660377

roc_c

pre_ltx vs post_ltx

group <- c("pre_ltx","post_ltx")
comparison_name <- paste0(group[1], " vs ",group[2])

model_name <- paste(comparison_name,level,segment)

# prepare the data
filt_ileum_uni_data <- binomial_prep(ileum_genus_tab,
                                     ileum_genus_taxa_tab,
                                     ileum_metadata,
                                     group, 
                                     usage="ml_clr")

Removing 46 ASV(s)

Removing 6 ASV(s)

# fit the model
gbm_model <- gbm_binomial(filt_ileum_uni_data,
                              sample_method = "atypboot",
                              outcome="Group", 
                              N=500,
                              reuse=TRUE,
                              file=model_name,
                              Q="Q1",
                              overfitting_check = TRUE)

# ROC curve
roc_c <- roc_curve(gbm_model, group)

# save the results
supplements_models[["models_summ"]][["gbm_model"]][[model_name]] <- gbm_model$model_summary
supplements_models[["roc_cs"]][["gbm_model"]][[model_name]] <- gbm_model$kfold_rocobjs

# see the results
gbm_model$model_summary %>% t()

                                       [,1]
n.trees                         500.0000000
interaction.depth                 1.0000000
shrinkage                         0.1000000
n.minobsinnode                   20.0000000
auc                               1.0000000
auc_optimism_corrected            0.5006084
auc_optimism_corrected_CIL        0.3838311
auc_optimism_corrected_CIU        0.6216013
accuracy                          1.0000000
accuracy_optimism_corrected       0.6020980
accuracy_optimism_corrected_CIL   0.5068451
accuracy_optimism_corrected_CIU   0.6918566

roc_c

post_ltx vs healthy

group <- c("post_ltx","healthy")
comparison_name <- paste0(group[1], " vs ",group[2])

model_name <- paste(comparison_name,level,segment)

# prepare the data
filt_ileum_uni_data <- binomial_prep(ileum_genus_tab,
                                     ileum_genus_taxa_tab,
                                     ileum_metadata,
                                     group, 
                                     usage="ml_clr")

Removing 45 ASV(s)

Removing 2 ASV(s)

# fit the model
gbm_model <- gbm_binomial(filt_ileum_uni_data,
                              sample_method = "atypboot",
                              outcome="Group", 
                              N=500,
                              reuse=TRUE,
                              file=model_name,
                              Q="Q1",
                              overfitting_check = TRUE)

# ROC curve
roc_c <- roc_curve(gbm_model, group)

# save the results
supplements_models[["models_summ"]][["gbm_model"]][[model_name]] <- gbm_model$model_summary
supplements_models[["roc_cs"]][["gbm_model"]][[model_name]] <- gbm_model$kfold_rocobjs

# see the results
gbm_model$model_summary %>% t()

                                       [,1]
n.trees                         200.0000000
interaction.depth                 1.0000000
shrinkage                         0.1000000
n.minobsinnode                   10.0000000
auc                               0.9998015
auc_optimism_corrected            0.5728982
auc_optimism_corrected_CIL        0.4579271
auc_optimism_corrected_CIU        0.6886252
accuracy                          0.9715640
accuracy_optimism_corrected       0.5907432
accuracy_optimism_corrected_CIL   0.4967593
accuracy_optimism_corrected_CIU   0.6845890

roc_c

4.1.2.2 Saving results

models_list <- list()

for (model_name in names(supplements_models$models_summ)){
  df <- do.call(rbind, supplements_models$models_summ[[model_name]])
  models_list[[model_name]] <- df
}

write.xlsx(models_list,
           file=file.path(path,paste0("supplements_models_",segment,".xlsx")),
           rowNames=TRUE)

5 Data Analysis - Colon

segment="colon"

5.1 Machine learning

path = "../results/Q1/models_overfitting_check"

5.1.1 ENET

model="enet"

5.1.1.1 Genus level

level="genus"

Aggregate taxa

genus_data <- aggregate_taxa(colon_asv_tab,
                             colon_taxa_tab,
                             taxonomic_level = level)

colon_genus_tab <- genus_data[[1]]
colon_genus_taxa_tab <- genus_data[[2]]

5.1.1.1.1 pre_ltx vs healthy

group <- c("pre_ltx","healthy")
comparison_name <- paste0(group[1], " vs ",group[2])

model_name <- paste(comparison_name,level,segment)

# prepare the data
filt_colon_uni_data <- binomial_prep(colon_genus_tab,
                                     colon_genus_taxa_tab,
                                     colon_metadata,
                                     group, 
                                     usage="ml_clr",
                                     patient = TRUE)

Removing 135 ASV(s)

Removing 10 ASV(s)

# fit the model
enet_model <- glmnet_binomial(filt_colon_uni_data,
                              sample_method = "atypboot",
                              outcome="Group", 
                              N=500,
                              clust_var="Patient",
                              reuse=TRUE,
                              file=model_name,
                              Q="Q1",overfitting_check = TRUE)


# ROC
roc_c <- roc_curve(enet_model, group)

# save the results
models_summ[[model_name]] <- enet_model$model_summary
models_cm[[model_name]] <- enet_model$conf_matrices$original
roc_cs[[model_name]] <- enet_model$kfold_rocobjs
betas[[model_name]] <- as.matrix(enet_model$betas)

# see the results
enet_model$model_summary %>% t()

                                     [,1]
alpha                           0.6000000
lambda                          0.1142179
auc                             0.5653292
auc_czech                       0.5312098
auc_no                          0.5910540
auc_optimism_corrected          0.5004172
auc_optimism_corrected_CIL      0.4204680
auc_optimism_corrected_CIU      0.5823541
accuracy                        0.6082725
accuracy_czech                        NaN
accuracy_no                     0.6423358
accuracy_optimism_corrected     0.5655580
accuracy_optimism_corrected_CIL 0.4813607
accuracy_optimism_corrected_CIU 0.6341662

enet_model$conf_matrices

$original
    1  
0 161 0
1 250 0

$czech
   1  
0 63 0
1 74 0

$no
    1  
0  98 0
1 176 0

enet_model$plot

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

roc_c

5.1.1.1.2 pre_ltx vs post_ltx

group <- c("pre_ltx","post_ltx")
comparison_name <- paste0(group[1], " vs ",group[2])

model_name <- paste(comparison_name,level,segment)

# prepare the data
filt_colon_uni_data <- binomial_prep(colon_genus_tab,
                                     colon_genus_taxa_tab,
                                     colon_metadata,
                                     group, 
                                     usage="ml_clr",
                                     patient = TRUE)

Removing 17 ASV(s)

Removing 10 ASV(s)

# fit the model
enet_model <- glmnet_binomial(filt_colon_uni_data,
                              sample_method = "atypboot",
                              outcome="Group", 
                              N=500,
                              clust_var="Patient",
                              reuse=TRUE,
                              file=model_name,
                              Q="Q1",
                              overfitting_check = TRUE)

# ROC
roc_c <- roc_curve(enet_model, group)

# save the results
models_summ[[model_name]] <- enet_model$model_summary
models_cm[[model_name]] <- enet_model$conf_matrices$original
roc_cs[[model_name]] <- enet_model$kfold_rocobjs
betas[[model_name]] <- as.matrix(enet_model$betas)

# see the results
enet_model$model_summary %>% t()

                                      [,1]
alpha                            0.0000000
lambda                          55.8015092
auc                              0.6231771
auc_czech                        0.6506435
auc_no                           0.6103118
auc_optimism_corrected           0.4986184
auc_optimism_corrected_CIL       0.4393670
auc_optimism_corrected_CIU       0.5579272
accuracy                         0.5833333
accuracy_czech                         NaN
accuracy_no                      0.6017192
accuracy_optimism_corrected      0.5620984
accuracy_optimism_corrected_CIL  0.5000000
accuracy_optimism_corrected_CIU  0.6269051

enet_model$conf_matrices

$original
    0  
0 350 0
1 250 0

$czech
    0  
0 140 0
1 111 0

$no
    0  
0 210 0
1 139 0

enet_model$plot

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

roc_c

5.1.1.1.3 post_ltx vs healthy

group <- c("post_ltx","healthy")
comparison_name <- paste0(group[1], " vs ",group[2])

model_name <- paste(comparison_name,level,segment)

# prepare the data
filt_colon_uni_data <- binomial_prep(colon_genus_tab,
                                     colon_genus_taxa_tab,
                                     colon_metadata,
                                     group, 
                                     usage="ml_clr",
                                     patient = TRUE)

Removing 71 ASV(s)

Removing 5 ASV(s)

# fit the model
enet_model <- glmnet_binomial(filt_colon_uni_data,
                              sample_method = "atypboot",
                              outcome="Group",
                              N=500,
                              clust_var="Patient",
                              reuse=TRUE,
                              file=model_name,
                              Q="Q1",
                              overfitting_check = TRUE)

# ROC
roc_c <- roc_curve(enet_model, group)

# save the results
models_summ[[model_name]] <- enet_model$model_summary
models_cm[[model_name]] <- enet_model$conf_matrices$original
roc_cs[[model_name]] <- enet_model$kfold_rocobjs
betas[[model_name]] <- as.matrix(enet_model$betas)

# see the results
enet_model$model_summary %>% t()

                                      [,1]
alpha                            0.0000000
lambda                          59.1733127
auc                              0.6228039
auc_czech                        0.5983574
auc_no                           0.6619718
auc_optimism_corrected           0.4917990
auc_optimism_corrected_CIL       0.4149382
auc_optimism_corrected_CIU       0.5685243
accuracy                         0.6849315
accuracy_czech                         NaN
accuracy_no                      0.6859903
accuracy_optimism_corrected      0.6668773
accuracy_optimism_corrected_CIL  0.6040824
accuracy_optimism_corrected_CIU  0.7232358

enet_model$conf_matrices

$original
    1  
0 161 0
1 350 0

$czech
    1  
0  96 0
1 208 0

$no
    1  
0  65 0
1 142 0

enet_model$plot

`geom_smooth()` using formula = 'y ~ x'

`geom_smooth()` using formula = 'y ~ x'

roc_c

5.1.1.2 Saving results

models_summ_df_colon <- do.call(rbind, 
  models_summ[grep(segment,names(models_summ),value = TRUE)])

write.csv(models_summ_df_colon,file.path(path,paste0("elastic_net_",segment,".csv")))

5.1.2 Supplementary models

5.1.2.1 CLR-transformed data

5.1.2.1.1 kNN

model="knn"

5.1.2.1.1.1 Genus level

level="genus"

Aggregate taxa

genus_data <- aggregate_taxa(colon_asv_tab,
                             colon_taxa_tab,
                             taxonomic_level = level)

colon_genus_tab <- genus_data[[1]]
colon_genus_taxa_tab <- genus_data[[2]]

pre_ltx vs healthy

group <- c("pre_ltx","healthy")
comparison_name <- paste0(group[1], " vs ",group[2])

model_name <- paste(comparison_name,level,segment)

# prepare the data
filt_colon_uni_data <- binomial_prep(colon_genus_tab,
                                     colon_genus_taxa_tab,
                                     colon_metadata,
                                     group, 
                                     usage="ml_clr",
                                     patient = TRUE)

Removing 135 ASV(s)

Removing 10 ASV(s)

# fit the model
knn_model <- knn_binomial(filt_colon_uni_data,
                              sample_method = "atypboot",
                              outcome="Group", 
                              N=500,
                          clust_var="Patient",
                              reuse=TRUE,
                              file=model_name,
                              Q="Q1",
                          overfitting_check = TRUE)

# ROC curve
roc_c <- roc_curve(knn_model, group)

# save the results
supplements_models[["models_summ"]][["knn_model"]][[model_name]] <- knn_model$model_summary
supplements_models[["roc_cs"]][["knn_model"]][[model_name]] <- knn_model$kfold_rocobjs


# see the results
knn_model$model_summary %>% t()

                                      [,1]
k                               27.0000000
auc                              0.6238261
auc_optimism_corrected           0.4942465
auc_optimism_corrected_CIL       0.3970386
auc_optimism_corrected_CIU       0.5861549
accuracy                         0.6277372
accuracy_optimism_corrected      0.5315829
accuracy_optimism_corrected_CIL  0.4477612
accuracy_optimism_corrected_CIU  0.6114949

roc_c

pre_ltx vs post_ltx

group <- c("pre_ltx","post_ltx")
comparison_name <- paste0(group[1], " vs ",group[2])

model_name <- paste(comparison_name,level,segment)

# prepare the data
filt_colon_uni_data <- binomial_prep(colon_genus_tab,
                                     colon_genus_taxa_tab,
                                     colon_metadata,
                                     group, 
                                     usage="ml_clr",
                                     patient = TRUE)

Removing 17 ASV(s)

Removing 10 ASV(s)

# fit the model
knn_model <- knn_binomial(filt_colon_uni_data,
                              sample_method = "atypboot",
                              outcome="Group", 
                              N=500,
                          clust_var="Patient",
                              reuse=TRUE,
                              file=model_name,
                              Q="Q1",
                          overfitting_check = TRUE)

# ROC curve
roc_c <- roc_curve(knn_model, group)

# save the results
supplements_models[["models_summ"]][["knn_model"]][[model_name]] <- knn_model$model_summary
supplements_models[["roc_cs"]][["knn_model"]][[model_name]] <- knn_model$kfold_rocobjs


# see the results
knn_model$model_summary %>% t()

                                      [,1]
k                               28.0000000
auc                              0.5659029
auc_optimism_corrected           0.4891374
auc_optimism_corrected_CIL       0.4230055
auc_optimism_corrected_CIU       0.5574708
accuracy                         0.6050000
accuracy_optimism_corrected      0.5260455
accuracy_optimism_corrected_CIL  0.4639652
accuracy_optimism_corrected_CIU  0.5877226

roc_c

post_ltx vs healthy

group <- c("post_ltx","healthy")
comparison_name <- paste0(group[1], " vs ",group[2])

model_name <- paste(comparison_name,level,segment)

# prepare the data
filt_colon_uni_data <- binomial_prep(colon_genus_tab,
                                     colon_genus_taxa_tab,
                                     colon_metadata,
                                     group, 
                                     usage="ml_clr",
                                     patient = TRUE)

Removing 71 ASV(s)

Removing 5 ASV(s)

# fit the model
knn_model <- knn_binomial(filt_colon_uni_data,
                              sample_method = "atypboot",
                              outcome="Group", 
                              N=500,
                          clust_var="Patient",
                              reuse=TRUE,
                              file=model_name,
                              Q="Q1",
                          overfitting_check = TRUE)

# ROC curve
roc_c <- roc_curve(knn_model, group)

# save the results
supplements_models[["models_summ"]][["knn_model"]][[model_name]] <- knn_model$model_summary
supplements_models[["roc_cs"]][["knn_model"]][[model_name]] <- knn_model$kfold_rocobjs


# see the results
knn_model$model_summary %>% t()

                                      [,1]
k                               30.0000000
auc                              0.5929636
auc_optimism_corrected           0.4840235
auc_optimism_corrected_CIL       0.3986432
auc_optimism_corrected_CIU       0.5643297
accuracy                         0.6829746
accuracy_optimism_corrected      0.6282546
accuracy_optimism_corrected_CIL  0.5403005
accuracy_optimism_corrected_CIU  0.6962623

roc_c

5.1.2.1.2 RF

model="rf"

5.1.2.1.2.1 Genus level

level="genus"

Aggregate taxa

genus_data <- aggregate_taxa(colon_asv_tab,
                             colon_taxa_tab,
                             taxonomic_level = level)

colon_genus_tab <- genus_data[[1]]
colon_genus_taxa_tab <- genus_data[[2]]

pre_ltx vs healthy

group <- c("pre_ltx","healthy")
comparison_name <- paste0(group[1], " vs ",group[2])

model_name <- paste(comparison_name,level,segment)

# prepare the data
filt_colon_uni_data <- binomial_prep(colon_genus_tab,
                                     colon_genus_taxa_tab,
                                     colon_metadata,
                                     group, 
                                     usage="ml_clr",
                                     patient = TRUE)

Removing 135 ASV(s)

Removing 10 ASV(s)

# fit the model
rf_model <- rf_binomial(filt_colon_uni_data,
                              sample_method = "atypboot",
                              outcome="Group", 
                              N=500,
                        clust_var="Patient",
                              reuse=TRUE,
                              file=model_name,
                              Q="Q1",
                        overfitting_check = TRUE)

# ROC curve
roc_c <- roc_curve(rf_model, group)

# save the results
supplements_models[["models_summ"]][["rf_model"]][[model_name]] <- rf_model$model_summary
supplements_models[["roc_cs"]][["rf_model"]][[model_name]] <- rf_model$kfold_rocobjs

# see the results
rf_model$model_summary %>% t()

                                [,1]       
mtry                            "171"      
splitrule                       "gini"     
min.node.size                   "5"        
auc                             "1"        
auc_optimism_corrected          "0.5431729"
auc_optimism_corrected_CIL      "0.4415343"
auc_optimism_corrected_CIU      "0.5874833"
accuracy                        "1"        
accuracy_optimism_corrected     "0.5147448"
accuracy_optimism_corrected_CIL "0.4917254"
accuracy_optimism_corrected_CIU "0.5364672"

roc_c

pre_ltx vs post_ltx

group <- c("pre_ltx","post_ltx")
comparison_name <- paste0(group[1], " vs ",group[2])

model_name <- paste(comparison_name,level,segment)

# prepare the data
filt_colon_uni_data <- binomial_prep(colon_genus_tab,
                                     colon_genus_taxa_tab,
                                     colon_metadata,
                                     group, 
                                     usage="ml_clr",
                                     patient = TRUE)

Removing 17 ASV(s)

Removing 10 ASV(s)

# fit the model
rf_model <- rf_binomial(filt_colon_uni_data,
                              sample_method = "atypboot",
                              outcome="Group", 
                              N=500,
                        clust_var="Patient",
                              reuse=TRUE,
                              file=model_name,
                              Q="Q1",
                        overfitting_check = TRUE)

# ROC curve
roc_c <- roc_curve(rf_model, group)

# save the results
supplements_models[["models_summ"]][["rf_model"]][[model_name]] <- rf_model$model_summary
supplements_models[["roc_cs"]][["rf_model"]][[model_name]] <- rf_model$kfold_rocobjs

# see the results
rf_model$model_summary %>% t()

                                [,1]       
mtry                            "43"       
splitrule                       "gini"     
min.node.size                   "2"        
auc                             "1"        
auc_optimism_corrected          "0.5455831"
auc_optimism_corrected_CIL      "0.4784397"
auc_optimism_corrected_CIU      "0.6023098"
accuracy                        "1"        
accuracy_optimism_corrected     "0.5610166"
accuracy_optimism_corrected_CIL "0.502336" 
accuracy_optimism_corrected_CIU "0.6167833"

roc_c

post_ltx vs healthy

group <- c("post_ltx","healthy")
comparison_name <- paste0(group[1], " vs ",group[2])

model_name <- paste(comparison_name,level,segment)

# prepare the data
filt_colon_uni_data <- binomial_prep(colon_genus_tab,
                                     colon_genus_taxa_tab,
                                     colon_metadata,
                                     group, 
                                     usage="ml_clr",
                                     patient = TRUE)

Removing 71 ASV(s)

Removing 5 ASV(s)

# fit the model
rf_model <- rf_binomial(filt_colon_uni_data,
                              sample_method = "atypboot",
                              outcome="Group", 
                              N=500,
                        clust_var="Patient",
                              reuse=TRUE,
                              file=model_name,
                              Q="Q1",
                        overfitting_check = TRUE)

# ROC curve
roc_c <- roc_curve(rf_model, group)

# save the results
supplements_models[["models_summ"]][["rf_model"]][[model_name]] <- rf_model$model_summary
supplements_models[["roc_cs"]][["rf_model"]][[model_name]] <- rf_model$kfold_rocobjs

# see the results
rf_model$model_summary %>% t()

                                [,1]       
mtry                            "1"        
splitrule                       "gini"     
min.node.size                   "2"        
auc                             "1"        
auc_optimism_corrected          "0.4941543"
auc_optimism_corrected_CIL      "0.4149008"
auc_optimism_corrected_CIU      "0.5638735"
accuracy                        "1"        
accuracy_optimism_corrected     "0.6589997"
accuracy_optimism_corrected_CIL "0.6017515"
accuracy_optimism_corrected_CIU "0.7087798"

roc_c

5.1.2.1.3 GBoost

model="gb"

5.1.2.1.3.1 Genus level

level="genus"

Aggregate taxa

genus_data <- aggregate_taxa(colon_asv_tab,
                             colon_taxa_tab,
                             taxonomic_level = level)

colon_genus_tab <- genus_data[[1]]
colon_genus_taxa_tab <- genus_data[[2]]

pre_ltx vs healthy

group <- c("pre_ltx","healthy")
comparison_name <- paste0(group[1], " vs ",group[2])

model_name <- paste(comparison_name,level,segment)

# prepare the data
filt_colon_uni_data <- binomial_prep(colon_genus_tab,
                                     colon_genus_taxa_tab,
                                     colon_metadata,
                                     group, 
                                     usage="ml_clr",
                                     patient = TRUE)

Removing 135 ASV(s)

Removing 10 ASV(s)

# fit the model
gbm_model <- gbm_binomial(filt_colon_uni_data,
                              sample_method = "atypboot",
                              outcome="Group", 
                              N=500,
                          clust_var="Patient",
                              reuse=TRUE,
                              file=model_name,
                              Q="Q1",
                          overfitting_check = TRUE)

# ROC curve
roc_c <- roc_curve(gbm_model, group)

# save the results
supplements_models[["models_summ"]][["gbm_model"]][[model_name]] <- gbm_model$model_summary
supplements_models[["roc_cs"]][["gbm_model"]][[model_name]] <- gbm_model$kfold_rocobjs

# see the results
gbm_model$model_summary %>% t()

                                       [,1]
n.trees                         100.0000000
interaction.depth                 1.0000000
shrinkage                         0.1000000
n.minobsinnode                   10.0000000
auc                               0.8847329
auc_optimism_corrected            0.5070577
auc_optimism_corrected_CIL        0.4250892
auc_optimism_corrected_CIU        0.5821484
accuracy                          0.7639903
accuracy_optimism_corrected       0.5066920
accuracy_optimism_corrected_CIL   0.4313082
accuracy_optimism_corrected_CIU   0.5769537

roc_c

pre_ltx vs post_ltx

group <- c("pre_ltx","post_ltx")
comparison_name <- paste0(group[1], " vs ",group[2])

model_name <- paste(comparison_name,level,segment)

# prepare the data
filt_colon_uni_data <- binomial_prep(colon_genus_tab,
                                     colon_genus_taxa_tab,
                                     colon_metadata,
                                     group, 
                                     usage="ml_clr",
                                     patient = TRUE)

Removing 17 ASV(s)

Removing 10 ASV(s)

# fit the model
gbm_model <- gbm_binomial(filt_colon_uni_data,
                              sample_method = "atypboot",
                              outcome="Group", 
                              N=500,
                          clust_var="Patient",
                              reuse=TRUE,
                              file=model_name,
                              Q="Q1",
                          overfitting_check = TRUE)

# ROC curve
roc_c <- roc_curve(gbm_model, group)

# save the results
supplements_models[["models_summ"]][["gbm_model"]][[model_name]] <- gbm_model$model_summary
supplements_models[["roc_cs"]][["gbm_model"]][[model_name]] <- gbm_model$kfold_rocobjs

# see the results
gbm_model$model_summary %>% t()

                                       [,1]
n.trees                         100.0000000
interaction.depth                 3.0000000
shrinkage                         0.1000000
n.minobsinnode                   10.0000000
auc                               0.9571657
auc_optimism_corrected            0.4865905
auc_optimism_corrected_CIL        0.4188913
auc_optimism_corrected_CIU        0.5515051
accuracy                          0.8816667
accuracy_optimism_corrected       0.5557342
accuracy_optimism_corrected_CIL   0.5023909
accuracy_optimism_corrected_CIU   0.6059283

roc_c

post_ltx vs healthy

group <- c("post_ltx","healthy")
comparison_name <- paste0(group[1], " vs ",group[2])

model_name <- paste(comparison_name,level,segment)

# prepare the data
filt_colon_uni_data <- binomial_prep(colon_genus_tab,
                                     colon_genus_taxa_tab,
                                     colon_metadata,
                                     group, 
                                     usage="ml_clr",
                                     patient = TRUE)

Removing 71 ASV(s)

Removing 5 ASV(s)

# fit the model
gbm_model <- gbm_binomial(filt_colon_uni_data,
                              sample_method = "atypboot",
                              outcome="Group", 
                              N=500,
                          clust_var="Patient",
                              reuse=TRUE,
                              file=model_name,
                              Q="Q1",
                          overfitting_check = TRUE)

# ROC curve
roc_c <- roc_curve(gbm_model, group)

# save the results
supplements_models[["models_summ"]][["gbm_model"]][[model_name]] <- gbm_model$model_summary
supplements_models[["roc_cs"]][["gbm_model"]][[model_name]] <- gbm_model$kfold_rocobjs

# see the results
gbm_model$model_summary %>% t()

                                       [,1]
n.trees                         100.0000000
interaction.depth                 1.0000000
shrinkage                         0.1000000
n.minobsinnode                   20.0000000
auc                               0.8581012
auc_optimism_corrected            0.4927075
auc_optimism_corrected_CIL        0.4162465
auc_optimism_corrected_CIU        0.5710246
accuracy                          0.7592955
accuracy_optimism_corrected       0.5687969
accuracy_optimism_corrected_CIL   0.4943182
accuracy_optimism_corrected_CIU   0.6405416

roc_c

5.1.2.2 Saving results

models_list <- list()

for (model_name in names(supplements_models$models_summ)){
  df <- do.call(rbind, supplements_models$models_summ[[model_name]])
  models_list[[model_name]] <- df
}

write.xlsx(models_list,
           file=file.path(path,paste0("supplements_models_",segment,".xlsx")),
           rowNames=TRUE)

6 Overview

models_summ_df <- do.call(rbind, 
  models_summ)

# Build final dataframe
models_list[["enet_model"]] <- models_summ_df
final_df <- tibble(row_names = rownames(models_list[[1]]))

# Loop through models and extract required values
for (model_name in names(models_list)) {
  model_df <- models_list[[model_name]]
  
  # Combine AUC_optimism_corrected with its CI values
  final_df[[model_name]] <- paste0(
    round(model_df$auc_optimism_corrected, 3), 
    " (", round(model_df$auc_optimism_corrected_CIL, 3), "; ", 
    round(model_df$auc_optimism_corrected_CIU, 3), ")"
  )
}

# see the results
knitr::kable(final_df, caption="All models")

All models

row_names	knn_model	rf_model	gbm_model	enet_model
pre_ltx vs healthy genus terminal_ileum	0.511 (0.361; 0.662)	0.437 (0.298; 0.575)	0.432 (0.294; 0.538)	0.448 (0.324; 0.595)
pre_ltx vs post_ltx genus terminal_ileum	0.469 (0.355; 0.593)	0.417 (0.312; 0.517)	0.501 (0.384; 0.622)	0.444 (0.326; 0.561)
post_ltx vs healthy genus terminal_ileum	0.479 (0.354; 0.597)	0.476 (0.368; 0.595)	0.573 (0.458; 0.689)	0.518 (0.393; 0.622)
pre_ltx vs healthy genus colon	0.494 (0.397; 0.586)	0.543 (0.442; 0.587)	0.507 (0.425; 0.582)	0.5 (0.42; 0.582)
pre_ltx vs post_ltx genus colon	0.489 (0.423; 0.557)	0.546 (0.478; 0.602)	0.487 (0.419; 0.552)	0.499 (0.439; 0.558)
post_ltx vs healthy genus colon	0.484 (0.399; 0.564)	0.494 (0.415; 0.564)	0.493 (0.416; 0.571)	0.492 (0.415; 0.569)