Code du mémoire — pipeline R et Python
1 Objectif du document
Ce document reprend le code utilisé dans le cadre du mémoire : “Développement d’un score de discrimination de la dégradation clinique chez les patients hospitalisés au CHU de Liège : vers une amélioration de la sécurité des patients”. Il présente le pipeline de préparation des données, de modélisation et d’évaluation sous R et Python.
Les blocs de code sont affichés, mais ne sont pas exécutés lors du rendu. Aucun résultat, tableau ou graphique n’est donc généré dans cette annexe.
2 Pipeline R — préparation des données et régression logistique
Cette partie correspond au traitement principal réalisé dans RStudio : importation de la base pseudonymisée, création des snapshots, définition de la variable cible y_event, retrait des admissions directes en USI, construction du jeu de données au niveau séjour, séparation entraînement/test au niveau patient, prétraitement via recipes, modélisation par régression logistique pénalisée Elastic Net, transformation des probabilités en score et évaluation du modèle.
################################################################################
# MÉMOIRE MASTER 2
################################################################################
# ==============================================================================
# 00. CONFIGURATION GÉNÉRALE
# ==============================================================================
CONFIG <- list(
seed = 1234,
input_rds = "/Users/Sarah/Desktop/code_memoire/Database_pseudonymise.rds",
out_dir_database = "/Users/Sarah/Desktop/Memoire/database",
out_dir_python = "/Users/Sarah/Desktop/Memoire/python_file",
out_dir_results = "/Users/Sarah/Desktop/result_sarah",
split_prop = 0.80,
cv_folds = 5,
# Les seuils du score ne sont pas fixés à l'avance.
# Ils sont calculés automatiquement à partir des prédictions out-of-fold
# du jeu d'entraînement, en divisant les probabilités en 5 groupes égaux.
score_cuts = NULL,
th_score = 5,
tune_grid_size = 100
)
dir.create(CONFIG$out_dir_database, recursive = TRUE, showWarnings = FALSE)
dir.create(CONFIG$out_dir_python, recursive = TRUE, showWarnings = FALSE)
dir.create(CONFIG$out_dir_results, recursive = TRUE, showWarnings = FALSE)
# ==============================================================================
# 01. LIBRAIRIES
# ==============================================================================
library(data.table)
library(glmnet)
library(recipes)
library(readxl)
library(dplyr)
library(readr)
library(tidymodels)
library(tidyr)
library(janitor)
library(lubridate)
library(skimr)
library(ggplot2)
library(stringr)
library(tibble)
library(purrr)
# ==============================================================================
# 02. FONCTIONS UTILITAIRES GÉNÉRALES
# ==============================================================================
load_and_clean_data <- function(path) {
readRDS(path) %>%
clean_names() %>%
remove_empty("rows")
}
basic_controls <- function(df) {
list(
n_patients = n_distinct(df$pat_id_num),
n_sejours = n_distinct(df$sejour_id_num),
dim = dim(df),
classes = tibble(
variable = names(df),
type = vapply(df, function(x) paste(class(x), collapse = ","), character(1))
),
skim = skimr::skim(df)
)
}
find_cols_with_less_than <- function(df) {
tibble(
colonne = names(df),
inf = vapply(
df,
function(x) any(grepl("<", as.character(x), fixed = TRUE)),
logical(1)
)
) %>%
filter(inf)
}
to_num <- function(df, cols) {
df %>%
mutate(
across(
{{ cols }},
~ as.numeric(gsub(",", ".", as.character(.x)))
)
)
}
log_to_char <- function(df, cols) {
df %>%
mutate(across({{ cols }}, as.character))
}
convert_initial_types <- function(df) {
df %>%
to_num(starts_with(c("ews", "fc", "pad", "pas", "spo2", "temp")) & !any_of("spo2_j1")) %>%
log_to_char(starts_with(c("eoda")))
}
safe_max_datetime <- function(x) {
if (all(is.na(x))) as.POSIXct(NA) else max(x, na.rm = TRUE)
}
safe_min_datetime <- function(x) {
if (all(is.na(x))) as.POSIXct(NA) else min(x, na.rm = TRUE)
}
set_y_event_factor <- function(df) {
if (!"y_event" %in% names(df)) {
stop("La colonne y_event n'existe pas dans le dataframe.")
}
df %>%
mutate(
y_event = as.character(y_event),
y_event = factor(y_event, levels = c("1", "0"))
)
}
# ==============================================================================
# 03. SNAPSHOTS : PASSAGE WIDE -> LONG -> SNAPSHOT
# ==============================================================================
get_static_variables <- function() {
c(
"pat_id_num", "sejour_id_num", "datetime_adm", "campus_urgence",
"flg_sortie", "flg_deces", "sexe_id", "discipline_post_urg",
"campus_hospi", "dcd", "date_de_sortie", "age", "flg_adm_usi",
"first_adm_usi", "sortie_usi", "last_adm_sortie_usi", "type_last_adm_usi"
)
}
get_dynamic_pattern <- function() {
var_dyn_prefixe <- c(
"adm_o2", "eoda", "ews", "fc", "glasgow", "pad",
"pas", "spo2", "temp", "crp", "creatinine", "glucose", "potassium",
"troponine_i", "uree"
)
var_dyn_suffix <- c(
"adm_usi", "amd_usi", "j1", "j2", "j3", "j4", "j5", "j6", "j7",
"adm_urg", "amd_urg", "j_1", "j_2", "j_3", "j_4", "j_5", "j_6",
"j_7", "meth", "icu_j_1", "icu_j_2", "icu_j_3", "icu_j1", "icu_j2", "icu_j3"
)
paste0(
"^(",
paste(var_dyn_prefixe, collapse = "|"),
").*_(",
paste(var_dyn_suffix, collapse = "|"),
")$"
)
}
get_dynamic_cols <- function(df) {
pattern_dyn <- get_dynamic_pattern()
df %>%
select(matches(pattern_dyn)) %>%
names()
}
pivot_dynamic_long <- function(df) {
dyn_cols <- get_dynamic_cols(df)
df %>%
mutate(across(all_of(dyn_cols), as.character)) %>%
pivot_longer(
cols = all_of(dyn_cols),
names_to = c("var", "suffix"),
names_pattern = paste0(
"^(.*)_(",
"adm_usi|amd_usi|j1|j2|j3|j4|j5|j6|j7|",
"adm_urg|amd_urg|j_1|j_2|j_3|j_4|j_5|j_6|j_7|",
"meth|icu_j_1|icu_j_2|icu_j_3|icu_j1|icu_j2|icu_j3",
")$"
),
values_to = "valeur"
)
}
add_numeric_value_and_flag <- function(df_long) {
df_long %>%
mutate(
flag_signe = if_else(str_starts(valeur, "<"), 1L, 0L),
val_num = valeur %>%
str_remove("^<") %>%
str_replace(",", ".") %>%
as.numeric()
)
}
add_jour_index <- function(df_long) {
df_long %>%
mutate(
suffix_clean = tolower(trimws(suffix)),
var_clean = tolower(var)
) %>%
mutate(
jour_index = case_when(
suffix_clean %in% c("adm_urg", "amd_urg") ~ 0L,
suffix_clean %in% c("j1", "j_1") & !str_detect(var_clean, "(icu|usi)") ~ 1L,
suffix_clean %in% c("j2", "j_2") & !str_detect(var_clean, "(icu|usi)") ~ 2L,
suffix_clean %in% c("j3", "j_3") & !str_detect(var_clean, "(icu|usi)") ~ 3L,
(
suffix_clean %in% c("j_3", "j3") & str_detect(var_clean, "(icu|usi)")
) | (
suffix_clean %in% c("j_4", "j4") & !str_detect(var_clean, "(icu|usi)")
) ~ 4L,
(
suffix_clean %in% c("j_2", "j2") & str_detect(var_clean, "(icu|usi)")
) | (
suffix_clean %in% c("j_5", "j5") & !str_detect(var_clean, "(icu|usi)")
) ~ 5L,
(
suffix_clean %in% c("j_1", "j1") & str_detect(var_clean, "(icu|usi)")
) | (
suffix_clean %in% c("j6", "j_6") & !str_detect(var_clean, "(icu|usi)")
) ~ 6L,
suffix_clean %in% c("adm_usi", "amd_usi", "meth", "j7", "j_7") ~ 7L,
TRUE ~ NA_integer_
)
)
}
deduplicate_pre_icu_priority <- function(df_long) {
pre_icu_suffix <- c("j_3", "j3", "j_2", "j2", "j_1", "j1")
DT <- as.data.table(df_long)
DT[, var_l := tolower(var)]
DT[, suffix_l := tolower(suffix)]
DT[, is_pre_icu := suffix_l %in% pre_icu_suffix & grepl("(icu|usi|uci)", var_l)]
setorder(DT, pat_id_num, sejour_id_num, jour_index, var, -is_pre_icu)
DT <- DT[, .SD[1], by = .(pat_id_num, sejour_id_num, jour_index, var)]
DT[, c("var_l", "suffix_l", "is_pre_icu") := NULL]
as_tibble(DT)
}
pivot_snapshot_wide <- function(df_long) {
df_snapshot <- df_long %>%
filter(!is.na(jour_index)) %>%
select(
pat_id_num, sejour_id_num, jour_index,
var, val_num, flag_signe, valeur
) %>%
pivot_wider(
id_cols = c(pat_id_num, sejour_id_num, jour_index),
names_from = var,
values_from = c(flag_signe, val_num, valeur),
names_glue = "{var}_{.value}"
)
id_cols <- c("pat_id_num", "sejour_id_num", "jour_index")
dyn_cols <- setdiff(names(df_snapshot), id_cols)
ordre_dyn <- tibble(nom = dyn_cols) %>%
separate(
nom,
into = c("var", "type"),
sep = "_(?=(flag_signe|val_num|valeur)$)",
remove = FALSE,
extra = "merge",
fill = "right"
) %>%
arrange(
var,
match(type, c("flag_signe", "val_num", "valeur"))
) %>%
pull(nom)
df_snapshot %>%
select(all_of(id_cols), all_of(ordre_dyn))
}
add_static_and_followup <- function(df_snapshot, df_clean) {
var_stats <- get_static_variables()
df_stats <- df_clean %>%
select(all_of(var_stats))
df_snapshot %>%
left_join(df_stats, by = c("pat_id_num", "sejour_id_num")) %>%
mutate(
snapshot_time = as_datetime(datetime_adm) + days(jour_index),
fin_suivi = pmin(date_de_sortie, dcd, na.rm = TRUE)
) %>%
filter(is.na(fin_suivi) | snapshot_time <= fin_suivi)
}
create_event_table <- function(df_snapshot) {
df_snapshot %>%
group_by(pat_id_num, sejour_id_num) %>%
summarise(
fin_suivi = safe_max_datetime(fin_suivi),
dcd = safe_max_datetime(dcd),
first_adm_usi = safe_min_datetime(first_adm_usi),
.groups = "drop"
) %>%
mutate(
date_deces_valid = case_when(
is.na(dcd) ~ as.POSIXct(NA),
dcd <= fin_suivi ~ dcd,
TRUE ~ as.POSIXct(NA)
),
date_usi_valid = case_when(
is.na(first_adm_usi) ~ as.POSIXct(NA),
first_adm_usi <= fin_suivi ~ first_adm_usi,
TRUE ~ as.POSIXct(NA)
),
event_time = case_when(
!is.na(date_deces_valid) & !is.na(date_usi_valid) ~ pmin(date_deces_valid, date_usi_valid),
!is.na(date_deces_valid) ~ date_deces_valid,
!is.na(date_usi_valid) ~ date_usi_valid,
TRUE ~ as.POSIXct(NA)
)
)
}
create_snapshot_dataset <- function(df_clean) {
df_long <- df_clean %>%
pivot_dynamic_long() %>%
add_numeric_value_and_flag() %>%
add_jour_index() %>%
deduplicate_pre_icu_priority()
df_snapshot <- df_long %>%
pivot_snapshot_wide() %>%
add_static_and_followup(df_clean)
evt_sejour <- create_event_table(df_snapshot)
list(
df_snapshot = df_snapshot,
evt_sejour = evt_sejour,
df_long = df_long
)
}
# ==============================================================================
# 04. PRÉPARATION ML : LEAKAGE, ICU DIRECT, FUSION ICU/NON-ICU
# ==============================================================================
remove_direct_icu_stays <- function(df_snapshot) {
sejours_direct <- df_snapshot %>%
select(pat_id_num, sejour_id_num, datetime_adm, first_adm_usi) %>%
distinct() %>%
mutate(
date_urg = as.Date(datetime_adm),
date_icu = as.Date(first_adm_usi),
direct_icu = !is.na(date_icu) & (date_icu == date_urg)
) %>%
filter(direct_icu) %>%
select(pat_id_num, sejour_id_num)
df_snapshot %>%
anti_join(sejours_direct, by = c("pat_id_num", "sejour_id_num"))
}
remove_leaky_columns <- function(df_snapshot) {
df_snapshot %>%
select(
-any_of(c(
"flg_sortie",
"flg_deces",
"flg_adm_usi",
"dcd",
"date_de_sortie",
"first_adm_usi",
"sortie_usi",
"fin_suivi",
"date_deces_valid",
"date_usi_valid"
))
)
}
merge_icu_and_non_icu_columns <- function(df_ml) {
icu_cols <- names(df_ml) %>%
str_subset("_icu_")
match_to_nonicu <- tibble(icu = icu_cols) %>%
mutate(non_icu = str_replace(icu, "_icu_", "_")) %>%
filter(non_icu %in% names(df_ml))
reduce(
.x = seq_len(nrow(match_to_nonicu)),
.init = df_ml,
.f = function(d, i) {
nonc <- match_to_nonicu$non_icu[i]
ic <- match_to_nonicu$icu[i]
d %>%
mutate(!!nonc := coalesce(.data[[ic]], .data[[nonc]])) %>%
select(-all_of(ic))
}
)
}
prepare_ml_snapshot <- function(df_snapshot) {
df_snapshot %>%
remove_direct_icu_stays() %>%
remove_leaky_columns() %>%
filter(jour_index != 7) %>%
merge_icu_and_non_icu_columns()
}
# ==============================================================================
# 05. PASSAGE SNAPSHOT-LEVEL -> SÉJOUR-LEVEL
# ==============================================================================
mode_value <- function(x) {
x <- x[!is.na(x) & x != ""]
if (length(x) == 0) NA_character_ else names(sort(table(x), decreasing = TRUE))[1]
}
create_sejour_level_dataset <- function(df_ml, evt_sejour) {
sejour_event <- evt_sejour %>%
select(pat_id_num, sejour_id_num, date_usi_valid, date_deces_valid, event_time) %>%
mutate(
date_event_valid = event_time,
is_event = !is.na(date_event_valid),
event_type = case_when(
!is.na(date_usi_valid) & !is.na(date_deces_valid) & date_usi_valid <= date_deces_valid ~ "USI",
!is.na(date_usi_valid) & !is.na(date_deces_valid) & date_deces_valid < date_usi_valid ~ "deces",
!is.na(date_usi_valid) ~ "USI",
!is.na(date_deces_valid) ~ "deces",
TRUE ~ "aucun"
)
)
df_ml2 <- df_ml %>%
left_join(sejour_event, by = c("pat_id_num", "sejour_id_num")) %>%
mutate(
is_event = if_else(is.na(is_event), FALSE, is_event),
pre_event = is_event & !is.na(snapshot_time) & snapshot_time < date_event_valid
)
cat_cols <- c("adm_o2_valeur", "eoda_valeur")
num_cols <- names(df_ml2) %>% str_subset("_val_num$")
flag_cols <- names(df_ml2) %>% str_subset("_flag_signe$")
stat_var <- c(
"campus_urgence", "sexe_id", "discipline_post_urg", "campus_hospi",
"age", "datetime_adm", "last_adm_sortie_usi", "type_last_adm_usi",
"snapshot_time", "date_usi_valid", "date_deces_valid",
"date_event_valid", "event_time", "event_type", "is_event", "pre_event"
)
event_last_data <- df_ml2 %>%
filter(pre_event) %>%
arrange(pat_id_num, sejour_id_num, snapshot_time) %>%
group_by(pat_id_num, sejour_id_num) %>%
slice_tail(n = 1) %>%
ungroup() %>%
mutate(y_event = 1L)
non_event_mean <- df_ml2 %>%
filter(!is_event) %>%
group_by(pat_id_num, sejour_id_num) %>%
summarise(
across(all_of(num_cols), ~ mean(.x, na.rm = TRUE)),
across(all_of(flag_cols), ~ max(.x, na.rm = TRUE)),
across(all_of(intersect(cat_cols, names(df_ml2))), mode_value),
.groups = "drop"
) %>%
mutate(
across(all_of(num_cols), ~ ifelse(is.nan(.x), NA_real_, .x)),
across(all_of(flag_cols), ~ ifelse(is.infinite(.x), NA_real_, .x)),
y_event = 0L
)
df_stat_snapshot <- df_ml2 %>%
arrange(pat_id_num, sejour_id_num, snapshot_time) %>%
group_by(pat_id_num, sejour_id_num) %>%
slice_head(n = 1) %>%
ungroup() %>%
select(pat_id_num, sejour_id_num, all_of(intersect(stat_var, names(df_ml2))))
bind_rows(non_event_mean, event_last_data) %>%
distinct(pat_id_num, sejour_id_num, .keep_all = TRUE) %>%
select(-any_of(stat_var)) %>%
left_join(df_stat_snapshot, by = c("pat_id_num", "sejour_id_num")) %>%
select(-any_of(c("is_event", "pre_event", "jour_index")))
}
make_patient_event_table <- function(df_final) {
df_final %>%
group_by(pat_id_num) %>%
summarise(
any_event = as.integer(any(y_event == 1, na.rm = TRUE)),
.groups = "drop"
)
}
split_train_test_by_patient <- function(df_final, prop = 0.8, seed = 1234) {
patient_event <- make_patient_event_table(df_final)
set.seed(seed)
split_pat <- initial_split(
patient_event,
prop = prop,
strata = any_event
)
train_pat <- training(split_pat)
test_pat <- testing(split_pat)
df_train <- df_final %>%
semi_join(train_pat, by = "pat_id_num")
df_test <- df_final %>%
semi_join(test_pat, by = "pat_id_num")
checks <- list(
no_overlap_df = length(intersect(df_train$pat_id_num, df_test$pat_id_num)) == 0,
no_overlap_pat = length(intersect(train_pat$pat_id_num, test_pat$pat_id_num)) == 0,
all_patients = length(unique(patient_event$pat_id_num)) ==
length(unique(train_pat$pat_id_num)) + length(unique(test_pat$pat_id_num)),
all_stays = nrow(df_final) == nrow(df_train) + nrow(df_test),
prevalence_final = prop.table(table(df_final$y_event)) * 100,
prevalence_train = prop.table(table(df_train$y_event)) * 100,
prevalence_test = prop.table(table(df_test$y_event)) * 100,
prevalence_patient_final = prop.table(table(patient_event$any_event)) * 100,
prevalence_patient_train = prop.table(table(train_pat$any_event)) * 100,
prevalence_patient_test = prop.table(table(test_pat$any_event)) * 100
)
list(
df_train = df_train,
df_test = df_test,
patient_event = patient_event,
train_pat = train_pat,
test_pat = test_pat,
checks = checks
)
}
# ==============================================================================
# 06. NETTOYAGE PRÉ-RECIPE : FLAGS, VALEUR/VAL_NUM, MISSINGNESS
# ==============================================================================
drop_useless_flags <- function(df_final, df_train, df_test) {
flag_count <- df_train %>%
select(ends_with("_flag_signe")) %>%
summarise(across(everything(), ~ n_distinct(., na.rm = TRUE)))
flag_vals <- as.numeric(flag_count[1, ])
flag_names <- colnames(flag_count)
flag_to_drop <- flag_names[flag_vals <= 1]
list(
df_final = df_final %>% select(-all_of(flag_to_drop)),
df_train = df_train %>% select(-all_of(flag_to_drop)),
df_test = df_test %>% select(-all_of(flag_to_drop)),
flag_to_drop = flag_to_drop
)
}
drop_redundant_value_columns <- function(df_final, df_train, df_test) {
dyn_cols <- names(df_final) %>%
str_subset("_(valeur|val_num|flag_signe)$")
dyn_prefixes <- dyn_cols %>%
str_replace("_(valeur|val_num|flag_signe)$", "") %>%
unique()
num_prefix <- dyn_prefixes[
sapply(dyn_prefixes, function(v) {
col <- paste0(v, "_val_num")
col %in% names(df_final) && any(!is.na(df_final[[col]]))
})
]
text_prefix <- setdiff(dyn_prefixes, num_prefix)
drop_valeur_num <- paste0(num_prefix, "_valeur")
drop_val_num_text <- paste0(text_prefix, "_val_num")
cols_to_drop <- c(drop_valeur_num, drop_val_num_text)
list(
df_final = df_final %>% select(-all_of(intersect(cols_to_drop, names(df_final)))),
df_train = df_train %>% select(-all_of(intersect(cols_to_drop, names(df_train)))),
df_test = df_test %>% select(-all_of(intersect(cols_to_drop, names(df_test)))),
cols_to_drop = cols_to_drop,
num_prefix = num_prefix,
text_prefix = text_prefix
)
}
make_missingness_table <- function(df_train, var_sur_indication) {
df_train %>%
summarise(across(everything(), ~ mean(is.na(.x)) * 100)) %>%
pivot_longer(everything(), names_to = "variable", values_to = "pct_na") %>%
mutate(
pct_na = round(pct_na, 0),
var_type = if_else(variable %in% var_sur_indication, "sur_indication", "routine")
) %>%
arrange(desc(pct_na))
}
select_missingness_drops <- function(na_table) {
na_table %>%
filter(
(var_type == "routine" & pct_na > 5) |
(var_type == "sur_indication" & pct_na > 10)
) %>%
pull(variable) %>%
unique()
}
prepare_recipe_inputs <- function(df_final, df_train, df_test) {
cleaned_flags <- drop_useless_flags(df_final, df_train, df_test)
cleaned_values <- drop_redundant_value_columns(
cleaned_flags$df_final,
cleaned_flags$df_train,
cleaned_flags$df_test
)
var_sur_indication <- c(
"glasgow_total_val_num", "uree_val_num", "ews_val_num", "eoda_valeur",
"troponine_i_us_flag_signe", "troponine_i_us_val_num", "glucose_val_num",
"potassium_val_num", "creatinine_enzymatique_val_num",
"crp_flag_signe", "crp_val_num"
)
na_table <- make_missingness_table(cleaned_values$df_train, var_sur_indication)
vars_to_drop <- select_missingness_drops(na_table)
na_table_keep <- na_table %>%
filter(!variable %in% vars_to_drop)
list(
df_final = cleaned_values$df_final,
df_train = cleaned_values$df_train,
df_test = cleaned_values$df_test,
flag_to_drop = cleaned_flags$flag_to_drop,
cols_to_drop = cleaned_values$cols_to_drop,
na_table = na_table,
vars_to_drop = vars_to_drop,
na_table_keep = na_table_keep,
num_sur_indic = character(0),
cat_sur_indic = character(0)
)
}
# ==============================================================================
# 07. RECIPE TIDYMODELS
# ==============================================================================
add_clinical_bounds <- function(rec, df_train) {
if ("age" %in% names(df_train)) {
rec <- rec %>%
step_mutate(age = ifelse(!is.na(age) & (age < 18 | age > 120), NA, age))
}
if ("creatinine_enzymatique_val_num" %in% names(df_train)) {
rec <- rec %>%
step_mutate(
creatinine_enzymatique_val_num = ifelse(
!is.na(creatinine_enzymatique_val_num) &
(creatinine_enzymatique_val_num < 0.1 | creatinine_enzymatique_val_num > 25),
NA,
creatinine_enzymatique_val_num
)
)
}
if ("crp_val_num" %in% names(df_train)) {
rec <- rec %>%
step_mutate(crp_val_num = ifelse(!is.na(crp_val_num) & (crp_val_num < 0 | crp_val_num > 1000), NA, crp_val_num))
}
if ("ews_val_num" %in% names(df_train)) {
rec <- rec %>%
step_mutate(ews_val_num = ifelse(!is.na(ews_val_num) & (ews_val_num < 0 | ews_val_num > 20), NA, ews_val_num))
}
if ("fc_val_num" %in% names(df_train)) {
rec <- rec %>%
step_mutate(fc_val_num = ifelse(!is.na(fc_val_num) & (fc_val_num < 20 | fc_val_num > 250), NA, fc_val_num))
}
if ("glasgow_total_val_num" %in% names(df_train)) {
rec <- rec %>%
step_mutate(
glasgow_total_val_num = ifelse(
!is.na(glasgow_total_val_num) &
(glasgow_total_val_num < 3 | glasgow_total_val_num > 15),
NA,
glasgow_total_val_num
)
)
}
if ("glucose_val_num" %in% names(df_train)) {
rec <- rec %>%
step_mutate(glucose_val_num = ifelse(!is.na(glucose_val_num) & (glucose_val_num < 20 | glucose_val_num > 1000), NA, glucose_val_num))
}
if ("pad_val_num" %in% names(df_train)) {
rec <- rec %>%
step_mutate(pad_val_num = ifelse(!is.na(pad_val_num) & (pad_val_num < 20 | pad_val_num > 200), NA, pad_val_num))
}
if ("pas_val_num" %in% names(df_train)) {
rec <- rec %>%
step_mutate(pas_val_num = ifelse(!is.na(pas_val_num) & (pas_val_num < 40 | pas_val_num > 300), NA, pas_val_num))
}
if ("potassium_val_num" %in% names(df_train)) {
rec <- rec %>%
step_mutate(potassium_val_num = ifelse(!is.na(potassium_val_num) & (potassium_val_num < 1.0 | potassium_val_num > 10), NA, potassium_val_num))
}
if ("spo2_val_num" %in% names(df_train)) {
rec <- rec %>%
step_mutate(spo2_val_num = ifelse(!is.na(spo2_val_num) & (spo2_val_num < 30 | spo2_val_num > 100), NA, spo2_val_num))
}
if ("temp_val_num" %in% names(df_train)) {
rec <- rec %>%
step_mutate(temp_val_num = ifelse(!is.na(temp_val_num) & (temp_val_num < 25 | temp_val_num > 44), NA, temp_val_num))
}
if ("troponine_i_us_val_num" %in% names(df_train)) {
rec <- rec %>%
step_mutate(
troponine_i_us_val_num = ifelse(
!is.na(troponine_i_us_val_num) &
(troponine_i_us_val_num < 0 | troponine_i_us_val_num > 100000),
NA,
troponine_i_us_val_num
)
)
}
if ("uree_val_num" %in% names(df_train)) {
rec <- rec %>%
step_mutate(uree_val_num = ifelse(!is.na(uree_val_num) & (uree_val_num < 0 | uree_val_num > 500), NA, uree_val_num))
}
rec
}
build_recipe <- function(df_train, vars_to_drop, num_sur_indic, cat_sur_indic) {
df_train <- set_y_event_factor(df_train)
rec <- recipe(y_event ~ ., data = df_train) %>%
update_role(pat_id_num, sejour_id_num, new_role = "id")
rec <- add_clinical_bounds(rec, df_train)
rec %>%
step_rm(
snapshot_time,
date_usi_valid,
date_deces_valid,
date_event_valid,
event_time,
event_type,
any_of(vars_to_drop)
) %>%
step_mutate(
saison = case_when(
month(datetime_adm) %in% c(3, 4, 5) ~ "printemps",
month(datetime_adm) %in% c(6, 7, 8) ~ "été",
month(datetime_adm) %in% c(9, 10, 11) ~ "automne",
month(datetime_adm) %in% c(1, 2, 12) ~ "hiver",
TRUE ~ NA_character_
),
sexe = case_when(
sexe_id == "M" ~ 1,
sexe_id == "F" ~ 0,
TRUE ~ NA_real_
),
campus_urgence = case_when(
campus_urgence == "ST" ~ 1,
campus_urgence == "BY" ~ 0,
TRUE ~ NA_real_
)
) %>%
step_rm(datetime_adm, sexe_id) %>%
step_mutate_at(
ends_with("_flag_signe"),
fn = ~ as.integer(replace(., is.na(.), 0))
) %>%
step_impute_median(all_numeric_predictors()) %>%
step_string2factor(all_nominal_predictors()) %>%
step_novel(all_nominal_predictors()) %>%
step_impute_mode(all_nominal_predictors()) %>%
step_other(all_nominal_predictors(), threshold = 0.01, other = "autres") %>%
step_dummy(all_nominal_predictors(), one_hot = TRUE) %>%
step_zv(all_predictors()) %>%
step_normalize(all_numeric_predictors())
}
bake_and_export <- function(rec, df_train, df_test, out_dir_python) {
df_train <- set_y_event_factor(df_train)
df_test <- set_y_event_factor(df_test)
rec_prep <- prep(rec, training = df_train, retain = TRUE)
train_baked <- bake(rec_prep, new_data = df_train)
test_baked <- bake(rec_prep, new_data = df_test)
stopifnot(identical(names(train_baked), names(test_baked)))
stopifnot("y_event" %in% names(train_baked))
stopifnot("y_event" %in% names(test_baked))
stopifnot(max(colSums(is.na(train_baked))) == 0)
train_baked_py <- train_baked %>%
mutate(y_event = as.integer(as.character(y_event)))
test_baked_py <- test_baked %>%
mutate(y_event = as.integer(as.character(y_event)))
write.csv(train_baked_py, file.path(out_dir_python, "train_baked.csv"), row.names = FALSE)
write.csv(test_baked_py, file.path(out_dir_python, "test_baked.csv"), row.names = FALSE)
list(
rec_prep = rec_prep,
train_baked = train_baked,
test_baked = test_baked,
train_baked_py = train_baked_py,
test_baked_py = test_baked_py
)
}
# ==============================================================================
# 08. MODÉLISATION RÉGRESSION LOGISTIQUE ELASTIC NET
# ==============================================================================
#3 Pipeline Python — sélection de variables, PyCaret et évaluation
Cette partie correspond au pipeline Python réalisé à partir des données préparées sous RStudio. Les données importées sont déjà imputées, normalisées, encodées et séparées en jeux d’entraînement et de test. Les vérifications portent sur la cohérence des colonnes, l’absence de fuite patient, l’absence de variables leaky, le format numérique des prédicteurs et la présence de la variable cible y_event.
# ==============================================================================
# 00. LIBRAIRIES ET IMPORT DES DONNÉES
# ==============================================================================
import pandas as pd
import numpy as np
from sklearn.ensemble import RandomForestClassifier
from sklearn.model_selection import StratifiedGroupKFold, cross_val_score
from sklearn.metrics import (
average_precision_score,
make_scorer,
roc_auc_score,
confusion_matrix
)
from sklearn.calibration import calibration_curve
import matplotlib.pyplot as plt
from pycaret.classification import (
setup,
compare_models,
create_model,
tune_model,
finalize_model,
predict_model,
pull,
add_metric,
get_metrics
)
train = pd.read_csv("/Users/do.un.ia_/Desktop/python_file/train_baked.csv")
test = pd.read_csv("/Users/do.un.ia_/Desktop/python_file/test_baked.csv")
TARGET = "y_event"
# ==============================================================================
# 01. VÉRIFICATIONS DE BASE
# ==============================================================================
print("=" * 80)
print("VÉRIFICATIONS DE BASE")
print("=" * 80)
print("\nDimensions :")
print("Train :", train.shape)
print("Test :", test.shape)
if TARGET not in train.columns:
raise ValueError(f"La variable cible {TARGET} est absente du jeu d'entraînement.")
if TARGET not in test.columns:
raise ValueError(f"La variable cible {TARGET} est absente du jeu de test.")
print("\nVariable cible utilisée :", TARGET)
for name, df in [("train", train), ("test", test)]:
df[TARGET] = pd.to_numeric(df[TARGET], errors="coerce").astype(int)
print(f"\nDistribution cible — {name}")
print(df[TARGET].value_counts(dropna=False).sort_index())
print("Prévalence événement :", round(df[TARGET].mean() * 100, 2), "%")
unique_values = sorted(df[TARGET].dropna().unique())
if unique_values != [0, 1]:
raise ValueError(f"La cible dans {name} n'est pas binaire 0/1 : {unique_values}")
for df in [train, test]:
df["pat_id_num"] = df["pat_id_num"].astype(str)
df["sejour_id_num"] = df["sejour_id_num"].astype(str)
overlap_patients = set(train["pat_id_num"]).intersection(set(test["pat_id_num"]))
print("\nPatients communs train/test :", len(overlap_patients))
if len(overlap_patients) > 0:
print("ATTENTION : fuite patient train/test")
print(list(overlap_patients)[:20])
else:
print("OK : aucun patient commun entre train et test.")
train_cols = set(train.columns)
test_cols = set(test.columns)
only_train = sorted(train_cols - test_cols)
only_test = sorted(test_cols - train_cols)
print("\nColonnes uniquement dans train :", only_train)
print("Colonnes uniquement dans test :", only_test)
if len(only_train) > 0 or len(only_test) > 0:
raise ValueError("Train et test n'ont pas exactement les mêmes colonnes.")
print("OK : train et test ont les mêmes colonnes.")
vars_leaky = [
"snapshot_time",
"date_usi_valid",
"date_deces_valid",
"date_event_valid",
"event_time",
"event_type",
"datetime_adm",
"sexe_id",
"first_adm_usi",
"dcd",
"flg_deces",
"flg_adm_usi",
"date_de_sortie"
]
leaky_train = [v for v in vars_leaky if v in train.columns]
leaky_test = [v for v in vars_leaky if v in test.columns]
print("\nVariables leaky encore dans train :", leaky_train)
print("Variables leaky encore dans test :", leaky_test)
if len(leaky_train) > 0 or len(leaky_test) > 0:
raise ValueError("Certaines variables potentiellement leaky sont encore présentes.")
print("OK : aucune variable leaky détectée.")
ignore_cols = ["pat_id_num", "sejour_id_num", TARGET]
features_col = [c for c in train.columns if c not in ignore_cols]
X_train = train[features_col].copy()
X_test = test[features_col].copy()
y = train[TARGET].astype(int).values
groups = train["pat_id_num"].astype(str).values
print("\nVariables exclues des prédicteurs :", ignore_cols)
print("Présence de y_event dans features_col :", "y_event" in features_col)
if TARGET in features_col:
raise ValueError("ERREUR : la cible est encore présente dans les prédicteurs.")
else:
print("OK : la cible n'est pas présente dans les prédicteurs.")
non_numeric_train = X_train.select_dtypes(exclude=[np.number]).columns.tolist()
non_numeric_test = X_test.select_dtypes(exclude=[np.number]).columns.tolist()
print("\nColonnes non numériques train :", non_numeric_train)
print("Colonnes non numériques test :", non_numeric_test)
if len(non_numeric_train) > 0 or len(non_numeric_test) > 0:
raise ValueError("Il reste des colonnes non numériques dans les prédicteurs.")
print("OK : tous les prédicteurs sont numériques.")
na_train = X_train.isna().sum().sum()
na_test = X_test.isna().sum().sum()
print("\nNombre total de NA dans X_train :", na_train)
print("Nombre total de NA dans X_test :", na_test)
if na_train > 0 or na_test > 0:
print("ATTENTION : il reste des NA.")
print("\nTop NA train :")
print(X_train.isna().sum().sort_values(ascending=False).head(20))
print("\nTop NA test :")
print(X_test.isna().sum().sort_values(ascending=False).head(20))
else:
print("OK : aucun NA dans les prédicteurs.")
print("\n" + "=" * 80)
print("RÉSUMÉ")
print("=" * 80)
print("Variable cible :", TARGET)
print("Nombre de variables prédictives :", len(features_col))
print("Prévalence train :", round(train[TARGET].mean() * 100, 2), "%")
print("Prévalence test :", round(test[TARGET].mean() * 100, 2), "%")
print("Patients communs train/test :", len(overlap_patients))
print("NA train :", na_train)
print("NA test :", na_test)
print("\nVérifications terminées.")
# ==============================================================================
# 02. MAPPING DES COLONNES TRANSFORMÉES VERS LES VARIABLES ORIGINALES
# ==============================================================================
known_base_vars = [
"age", "sexe", "campus_urgence", "campus_hospi", "discipline_post_urg",
"saison", "last_adm_sortie_usi", "type_last_adm_usi", "adm_o2_valeur",
"eoda_valeur", "crp_flag_signe", "troponine_i_us_flag_signe", "crp_val_num",
"creatinine_enzymatique_val_num", "ews_val_num", "fc_val_num",
"glasgow_total_val_num", "glucose_val_num", "pad_val_num", "pas_val_num",
"potassium_val_num", "spo2_val_num", "temp_val_num",
"troponine_i_us_val_num", "uree_val_num"
]
known_base_vars = [
v for v in known_base_vars
if any(
c == v or c.startswith(v + "_") or c == "na_ind_" + v
for c in features_col
)
]
print("\nVariables originales reconnues :")
print(known_base_vars)
print("\nNombre de variables originales reconnues :", len(known_base_vars))
def map_column_to_base_variable(col, known_bases):
if col.startswith("na_ind_"):
stripped = col.replace("na_ind_", "", 1)
if stripped in known_bases:
return stripped
matches = [
b for b in known_bases
if stripped == b or stripped.startswith(b + "_")
]
if len(matches) > 0:
return max(matches, key=len)
return stripped
if col in known_bases:
return col
matches = [
b for b in known_bases
if col.startswith(b + "_")
]
if len(matches) > 0:
return max(matches, key=len)
return col
col_to_base = {
col: map_column_to_base_variable(col, known_base_vars)
for col in features_col
}
mapping_df = pd.DataFrame({
"column": list(col_to_base.keys()),
"base_variable": list(col_to_base.values())
}).sort_values(["base_variable", "column"])
print("\nMapping colonnes transformées vers variables originales :")
print(mapping_df.head(60).to_string(index=False))
print("\nNombre de colonnes transformées :", mapping_df["column"].nunique())
print("Nombre de variables originales après mapping :", mapping_df["base_variable"].nunique())
def aggregate_importance_by_base(feature_importance_df, col_to_base):
out = feature_importance_df.copy()
out["base_variable"] = out["feature"].map(col_to_base)
base_importance = (
out
.groupby("base_variable", as_index=False)
.agg(
importance=("importance", "sum"),
n_columns=("feature", "count"),
columns=("feature", lambda x: list(x))
)
.sort_values("importance", ascending=False)
.reset_index(drop=True)
)
return base_importance
def columns_for_base_variables(selected_base_vars, col_to_base):
return [
col for col, base in col_to_base.items()
if base in selected_base_vars
]
def make_pr_auc_scorer():
try:
return make_scorer(
average_precision_score,
response_method="predict_proba"
)
except TypeError:
return make_scorer(
average_precision_score,
needs_proba=True
)
# ==============================================================================
# 03. RANDOM FOREST POUR SÉLECTION DE VARIABLES
# ==============================================================================
outer_cv = StratifiedGroupKFold(
n_splits=5,
shuffle=True,
random_state=1451
)
inner_cv = StratifiedGroupKFold(
n_splits=5,
shuffle=True,
random_state=1453
)
candidate_k = [5, 7, 10]
pr_auc_scorer = make_pr_auc_scorer()
outer_results = []
for outer_fold, (train_idx, valid_idx) in enumerate(
outer_cv.split(X_train, y, groups),
start=1
):
print("\n" + "=" * 50)
print(f"OUTER FOLD {outer_fold}")
print("=" * 50)
X_outer_train = X_train.iloc[train_idx].copy()
y_outer_train = y[train_idx]
groups_outer_train = groups[train_idx]
X_outer_valid = X_train.iloc[valid_idx].copy()
y_outer_valid = y[valid_idx]
rf_selector = RandomForestClassifier(
n_estimators=500,
max_depth=None,
min_samples_leaf=5,
class_weight="balanced",
random_state=1516 + outer_fold,
n_jobs=-1
)
rf_selector.fit(X_outer_train, y_outer_train)
importance_columns_df = (
pd.DataFrame({
"feature": features_col,
"importance": rf_selector.feature_importances_
})
.sort_values("importance", ascending=False)
.reset_index(drop=True)
)
importance_base_df = aggregate_importance_by_base(
importance_columns_df,
col_to_base
)
print("\nTop 15 variables originales dans ce fold externe :")
print(
importance_base_df[
["base_variable", "importance", "n_columns"]
].head(15).to_string(index=False)
)
inner_results = {}
for k in candidate_k:
selected_base_vars = importance_base_df["base_variable"].iloc[:k].tolist()
selected_cols = columns_for_base_variables(selected_base_vars, col_to_base)
rf_inner = RandomForestClassifier(
n_estimators=500,
max_depth=None,
min_samples_leaf=5,
class_weight="balanced",
random_state=1526,
n_jobs=1
)
inner_scores = cross_val_score(
rf_inner,
X_outer_train[selected_cols],
y_outer_train,
cv=inner_cv,
groups=groups_outer_train,
scoring=pr_auc_scorer,
n_jobs=-1
)
inner_results[k] = {
"selected_base_vars": selected_base_vars,
"selected_cols": selected_cols,
"inner_pr_auc_mean": inner_scores.mean(),
"inner_pr_auc_std": inner_scores.std(),
"n_columns": len(selected_cols)
}
print(
f"\nTop {k:2d} variables originales "
f"({len(selected_cols)} colonnes transformées) "
f"Inner PR-AUC = {inner_scores.mean():.4f} ± {inner_scores.std():.4f}"
)
best_k_raw = max(
inner_results,
key=lambda kk: inner_results[kk]["inner_pr_auc_mean"]
)
best_score_inner = inner_results[best_k_raw]["inner_pr_auc_mean"]
best_se_inner = (
inner_results[best_k_raw]["inner_pr_auc_std"] /
np.sqrt(inner_cv.get_n_splits())
)
threshold_1se = best_score_inner - best_se_inner
eligible_k = [
kk for kk in sorted(inner_results)
if inner_results[kk]["inner_pr_auc_mean"] >= threshold_1se
]
best_k = min(eligible_k)
best_base_vars = inner_results[best_k]["selected_base_vars"]
best_cols = inner_results[best_k]["selected_cols"]
print(f"\nMeilleur k brut : Top {best_k_raw}")
print(f"Meilleur PR-AUC interne : {best_score_inner:.4f}")
print(f"Erreur standard du meilleur score : {best_se_inner:.4f}")
print(f"Seuil one standard error : {threshold_1se:.4f}")
print(f"k éligibles : {eligible_k}")
print(f"k choisi avec règle parcimonieuse : Top {best_k}")
rf_outer_final = RandomForestClassifier(
n_estimators=500,
max_depth=None,
min_samples_leaf=5,
class_weight="balanced",
random_state=1536,
n_jobs=-1
)
rf_outer_final.fit(X_outer_train[best_cols], y_outer_train)
p_valid = rf_outer_final.predict_proba(X_outer_valid[best_cols])[:, 1]
outer_pr_auc = average_precision_score(y_outer_valid, p_valid)
print(f"OUTER PR-AUC fold {outer_fold} = {outer_pr_auc:.4f}")
outer_results.append({
"outer_fold": outer_fold,
"best_k_raw": best_k_raw,
"best_se_inner": best_se_inner,
"threshold_1se": threshold_1se,
"eligible_k": eligible_k,
"best_k": best_k,
"best_score_inner": best_score_inner,
"score_top5": inner_results.get(5, {}).get("inner_pr_auc_mean", np.nan),
"score_top7": inner_results.get(7, {}).get("inner_pr_auc_mean", np.nan),
"score_top10": inner_results.get(10, {}).get("inner_pr_auc_mean", np.nan),
"best_base_vars": best_base_vars,
"best_cols": best_cols,
"n_columns": len(best_cols),
"inner_pr_auc_mean": inner_results[best_k]["inner_pr_auc_mean"],
"inner_pr_auc_std": inner_results[best_k]["inner_pr_auc_std"],
"outer_pr_auc": outer_pr_auc,
})
nested_results_df = pd.DataFrame(outer_results)
print("\n" + "=" * 50)
print("Résultats nested CV")
print("=" * 50)
print(
nested_results_df[
[
"outer_fold",
"best_k",
"score_top5",
"score_top7",
"score_top10",
"best_score_inner",
"n_columns",
"inner_pr_auc_mean",
"inner_pr_auc_std",
"outer_pr_auc"
]
].to_string(index=False)
)
print("\nPR-AUC externe moyenne :", round(nested_results_df["outer_pr_auc"].mean(), 4))
print("PR-AUC externe écart-type :", round(nested_results_df["outer_pr_auc"].std(), 4))
print("\nNombre de fois où chaque top k est sélectionné :")
print(nested_results_df["best_k"].value_counts().sort_index())
# ==============================================================================
# 04. SÉLECTION FINALE DES VARIABLES
# ==============================================================================
final_k = nested_results_df["best_k"].mode()[0]
print("\nStratégie choisie :")
print(f"Top {final_k} variables originales")
rf_selector_final = RandomForestClassifier(
n_estimators=500,
max_depth=None,
min_samples_leaf=5,
class_weight="balanced",
random_state=1740,
n_jobs=-1
)
rf_selector_final.fit(X_train, y)
importance_columns_final_df = (
pd.DataFrame({
"feature": features_col,
"importance": rf_selector_final.feature_importances_,
})
.sort_values("importance", ascending=False)
.reset_index(drop=True)
)
importance_base_final_df = aggregate_importance_by_base(
importance_columns_final_df,
col_to_base
)
print("\nTop variables originales finales :")
print(
importance_base_final_df[
["base_variable", "importance", "n_columns"]
].head(final_k).to_string(index=False)
)
final_base_vars = importance_base_final_df["base_variable"].iloc[:final_k].tolist()
final_cols = columns_for_base_variables(final_base_vars, col_to_base)
final_cols = [c for c in final_cols if c not in ["y_event"]]
print("\nVariables originales finales :")
print(final_base_vars)
print("\nColonnes transformées envoyées à PyCaret :")
print(final_cols)
print("\nNombre de variables originales :", len(final_base_vars))
print("Nombre de colonnes transformées :", len(final_cols))
print("Présence de y_event dans final_cols :", "y_event" in final_cols)
if "y_event" in final_cols:
raise ValueError("ERREUR : la cible est présente dans les colonnes envoyées à PyCaret.")
else:
print("OK : aucune cible envoyée à PyCaret comme prédicteur.")
def make_subset_for_pycaret(df, selected_cols):
return df[
["pat_id_num", "sejour_id_num"] + selected_cols + ["y_event"]
].copy()
train_selected = make_subset_for_pycaret(train, final_cols)
test_selected = make_subset_for_pycaret(test, final_cols)
print("\nShape train_selected :", train_selected.shape)
print("Shape test_selected :", test_selected.shape)
# ==============================================================================
# 05. COMPARAISON DES MODÈLES AVEC PYCARET
# ==============================================================================
from pycaret.classification import (
setup, compare_models, create_model, tune_model,
finalize_model, predict_model, pull, add_metric, get_metrics
)
from sklearn.model_selection import StratifiedGroupKFold
from sklearn.metrics import average_precision_score, roc_auc_score, confusion_matrix
from sklearn.calibration import calibration_curve
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
# Dataset provenant de la sélection de variables par Random Forest
train_pycaret = train_selected.copy()
test_pycaret = test_selected.copy()
# Vérifications
for df in [train_pycaret, test_pycaret]:
df["pat_id_num"] = df["pat_id_num"].astype(str)
df["sejour_id_num"] = df["sejour_id_num"].astype(str)
df["y_event"] = pd.to_numeric(df["y_event"], errors="coerce").astype(int)
print(train_pycaret.shape)
print(test_pycaret.shape)
print(
"Patients communs train/test :",
len(set(train_pycaret["pat_id_num"]).intersection(set(test_pycaret["pat_id_num"])))
)
print("Prévalence train :", train_pycaret["y_event"].mean())
print("Prévalence test :", test_pycaret["y_event"].mean())
# ------------------------------------------------------------------------------
# Setup PyCaret
# ------------------------------------------------------------------------------
models_setting = setup(
data=train_pycaret,
target="y_event",
session_id=1504,
fold_strategy=StratifiedGroupKFold(n_splits=5, shuffle=True, random_state=1507),
fold_groups="pat_id_num",
ignore_features=["pat_id_num", "sejour_id_num"],
fix_imbalance=False,
verbose=True,
use_gpu=False
)
# ------------------------------------------------------------------------------
# Ajouter la métrique PR-AUC
# ------------------------------------------------------------------------------
add_metric(
id="pr_auc",
name="PR-AUC",
score_func=average_precision_score,
greater_is_better=True
)
print(get_metrics()[["Name", "Display Name"]].to_string(index=False))
# ------------------------------------------------------------------------------
# Comparaison des modèles
# ------------------------------------------------------------------------------
best_models = compare_models(
sort="pr_auc",
n_select=5,
fold=5
)
results_compare = pull()
print(results_compare.head(20))
model_n1 = best_models[0]
# ------------------------------------------------------------------------------
# Vérification de la PR-AUC out-of-fold du modèle sélectionné
# ------------------------------------------------------------------------------
m = create_model(model_n1)
oof = predict_model(m, raw_score=True)
y_oof = oof["y_event"].astype(int).values
p_oof = oof["prediction_score_1"].astype(float).values
print("PR-AUC train OOF manual :", average_precision_score(y_oof, p_oof))
results_compare = pull()
print(results_compare.head(20))
# ------------------------------------------------------------------------------
# Tuning du modèle sélectionné
# ------------------------------------------------------------------------------
tuned_lgbm = tune_model(
model_n1,
optimize="pr_auc",
fold=5,
choose_better=True
)
tuning_results = pull()
print(tuning_results)
# ------------------------------------------------------------------------------
# Réentraînement final sur tout le jeu d'entraînement
# ------------------------------------------------------------------------------
final_lgbm = finalize_model(tuned_lgbm)
# ------------------------------------------------------------------------------
# Prédictions sur le jeu de test
# ------------------------------------------------------------------------------
test_pred = predict_model(final_lgbm, data=test_pycaret, raw_score=True)
print(test_pred["prediction_score_1"].min(), test_pred["prediction_score_1"].max())
y = test_pred["y_event"].astype(int).values
p = test_pred["prediction_score_1"].astype(float).values
print("PR-AUC test manual :", average_precision_score(y, p))
print(test_pred.columns.tolist())
cand = [
c for c in test_pred.columns
if c.lower() in ["score_1", "probability_1", "prediction_score_1"]
]
if not cand:
cand = [
c for c in test_pred.columns
if c.lower().endswith("_1") and ("score" in c.lower() or "prob" in c.lower())
]
if not cand:
raise ValueError(
f"Aucune colonne de proba/score classe 1 trouvée. "
f"Colonnes disponibles : {test_pred.columns.tolist()}"
)
p1_col = cand[0]
p1 = test_pred[p1_col].astype(float).values
print("Colonne proba utilisée =", p1_col)
y_true = test_pred["y_event"].astype(int).values
test_pred["y_true"] = y_true
# ==============================================================================
# 06. SCORE, MÉTRIQUES ET ÉVALUATION TEST
# ==============================================================================
print(models_setting.get_config("X_train").columns)
# Transformation des probabilités en score 1 à 5 à partir des quintiles OOF
score_cuts_quintiles_python = np.quantile(
p_oof,
[0.20, 0.40, 0.60, 0.80]
)
print("Seuils score Python par quintiles OOF :", score_cuts_quintiles_python)
test_pred["score_5"] = pd.cut(
pd.Series(p1, index=test_pred.index),
bins=[-np.inf, *score_cuts_quintiles_python, np.inf],
labels=[1, 2, 3, 4, 5],
right=True,
include_lowest=True
).astype(int)
print(test_pred["score_5"].value_counts().sort_index())
# Règle binaire : score 5 = profil à haut risque
test_pred["y_pred_score"] = (test_pred["score_5"] == 5).astype(int)
pr_auc_test = average_precision_score(y_true, p1)
roc_auc_test = roc_auc_score(y_true, p1)
print("PR-AUC :", round(pr_auc_test, 4))
print("ROC-AUC :", round(roc_auc_test, 4))
test_pred["y_true"] = test_pred["y_event"].astype(int)
test_pred["p"] = p1
brier = np.mean((test_pred["p"] - test_pred["y_true"]) ** 2)
print("Brier score :", round(brier, 4))
# ------------------------------------------------------------------------------
# Courbe de calibration
# ------------------------------------------------------------------------------
test_pred["bin"] = pd.qcut(test_pred["p"], 10, duplicates="drop")
calib_tbl = (
test_pred
.groupby("bin")
.agg(
p_mean=("p", "mean"),
obs_rate=("y_true", "mean"),
n=("y_true", "size")
)
.reset_index()
)
print(calib_tbl)
plt.figure()
plt.plot(calib_tbl["p_mean"], calib_tbl["obs_rate"], marker="o")
plt.plot([0, 1], [0, 1], linestyle="--")
plt.xlabel("Probabilité prédite")
plt.ylabel("Proportion observée")
plt.title("Courbe de calibration — test")
plt.show()
prob_true, prob_pred = calibration_curve(
y_true,
p1,
n_bins=10,
strategy="quantile"
)
plt.figure()
plt.plot(prob_pred, prob_true, marker="o")
plt.plot([0, 1], [0, 1], linestyle="--")
plt.xlabel("Probabilité prédite")
plt.ylabel("Proportion observée")
plt.title("Calibration curve — sklearn")
plt.show()
# ------------------------------------------------------------------------------
# Matrice de confusion
# ------------------------------------------------------------------------------
tn, fp, fn, tp = confusion_matrix(
test_pred["y_true"],
test_pred["y_pred_score"]
).ravel()
sens = tp / (tp + fn)
spec = tn / (tn + fp)
ppv = tp / (tp + fp) if (tp + fp) > 0 else float("nan")
npv = tn / (tn + fn)
print("TP:", tp, "FP:", fp, "TN:", tn, "FN:", fn)
print("Sensibilité:", round(sens, 3))
print("Spécificité:", round(spec, 3))
print("VPP:", round(ppv, 3))
print("VPN:", round(npv, 3))
# ==============================================================================
# 07. ANALYSE D'ERREURS, TENDANCE DU SCORE ET COURBE DÉCISIONNELLE
# ==============================================================================
# ------------------------------------------------------------------------------
# Analyse exploratoire des erreurs
# ------------------------------------------------------------------------------
variable_erreur = [
"crp_val_num",
"fc_val_num",
"creatinine_enzymatique_val_num",
"pad_val_num"
]
variable_erreur = [v for v in variable_erreur if v in test_pred.columns]
df_err = test_pred.copy()
print(test_pred.columns.tolist())
df_err["type_erreur"] = np.select(
[
(df_err["y_true"] == 1) & (df_err["y_pred_score"] == 1),
(df_err["y_true"] == 0) & (df_err["y_pred_score"] == 0),
(df_err["y_true"] == 0) & (df_err["y_pred_score"] == 1),
(df_err["y_true"] == 1) & (df_err["y_pred_score"] == 0),
],
["TP", "TN", "FP", "FN"],
default="?"
)
if len(variable_erreur) > 0:
analyse_explo_erreur = (
df_err[df_err["type_erreur"].isin(["FP", "TN", "TP", "FN"])]
.groupby("type_erreur")[variable_erreur]
.median(numeric_only=True)
.reset_index()
)
print(analyse_explo_erreur)
else:
print("Aucune des variables choisies n'est présente dans le dataset test.")
# ------------------------------------------------------------------------------
# Tendance du score
# ------------------------------------------------------------------------------
score_tbl_test = (
test_pred
.groupby("score_5", as_index=False)
.agg(
n=("score_5", "size"),
event_rate=("y_true", "mean"),
mean_p=("p", "mean")
)
.sort_values("score_5")
)
print(score_tbl_test)
# ------------------------------------------------------------------------------
# Courbe décisionnelle
# ------------------------------------------------------------------------------
y = test_pred["y_true"].astype(int).values
p = test_pred["p"].astype(float).values
N = len(y)
prev = y.mean()
thresholds = np.arange(0.001, 0.051, 0.001)
nb_model = []
nb_all = []
nb_none = []
for t in thresholds:
pred = (p >= t).astype(int)
tn, fp, fn, tp = confusion_matrix(y, pred).ravel()
nb = (tp / N) - (fp / N) * (t / (1 - t))
nb_model.append(nb)
nb_none.append(0.0)
nb_all.append(prev - (1 - prev) * (t / (1 - t)))
plt.figure()
plt.plot(thresholds, nb_model, label="Modèle")
plt.plot(thresholds, nb_all, linestyle="--", label="Traiter tous")
plt.plot(thresholds, nb_none, linestyle="--", label="Traiter aucun")
plt.xlabel("Seuil de décision")
plt.ylabel("Bénéfice net")
plt.title("Courbe décisionnelle — modèle Python")
plt.legend(loc="lower left")
plt.show()
# ==============================================================================
# 07B. VÉRIFICATIONS SUPPLÉMENTAIRES — LEAKAGE / COHÉRENCE PYCARET
# ==============================================================================
# Vérifier que la colonne de probabilité de la classe positive est bien présente
proba_col = [
c for c in test_pred.columns
if c.startswith("prediction_score_") and c.endswith("_1")
]
if proba_col:
p1 = test_pred[proba_col[0]].astype(float).values
elif "prediction_score_1" in test_pred.columns:
p1 = test_pred["prediction_score_1"].astype(float).values
else:
raise ValueError("Impossible de trouver la probabilité de la classe 1 dans test_pred.")
print("Colonne de probabilité utilisée :", proba_col[0] if proba_col else "prediction_score_1")
# ------------------------------------------------------------------------------
# 1. Vérifier les variables réellement utilisées par PyCaret
# ------------------------------------------------------------------------------
X_train_cols = list(models_setting.get_config("X_train").columns)
X_test_cols = list(models_setting.get_config("X_test").columns)
leak_cols = {"pat_id_num", "sejour_id_num", "y_event"}
print("Leak in X_train:", leak_cols.intersection(X_train_cols))
print("Leak in X_test :", leak_cols.intersection(X_test_cols))
if len(leak_cols.intersection(X_train_cols)) > 0 or len(leak_cols.intersection(X_test_cols)) > 0:
raise ValueError("ERREUR : une variable d'identification ou la cible est présente dans les prédicteurs PyCaret.")
else:
print("OK : aucune variable d'identification ni cible dans les prédicteurs PyCaret.")
# ------------------------------------------------------------------------------
# 2. Vérifier que train/test ont les mêmes colonnes hors cible
# ------------------------------------------------------------------------------
train_cols_check = set(train_pycaret.columns) - {"y_event"}
test_cols_check = set(test_pycaret.columns) - {"y_event"}
print("Columns only in train:", sorted(train_cols_check - test_cols_check)[:20])
print("Columns only in test :", sorted(test_cols_check - train_cols_check)[:20])
if len(train_cols_check - test_cols_check) > 0 or len(test_cols_check - train_cols_check) > 0:
raise ValueError("ERREUR : train_pycaret et test_pycaret n'ont pas les mêmes colonnes hors cible.")
else:
print("OK : mêmes colonnes dans train_pycaret et test_pycaret hors cible.")
# ------------------------------------------------------------------------------
# 3. Comparer PR-AUC OOF et PR-AUC test
# ------------------------------------------------------------------------------
pr_auc_oof = average_precision_score(y_oof, p_oof)
pr_auc_test = average_precision_score(y_true, p1)
print("PR-AUC train OOF :", round(pr_auc_oof, 4))
print("PR-AUC test :", round(pr_auc_test, 4))
print("Différence OOF - test :", round(pr_auc_oof - pr_auc_test, 4))
# ==============================================================================
# 08. ANALYSE SHAP — MODÈLE LIGHTGBM FINAL
# ==============================================================================
import shap
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
# ------------------------------------------------------------------------------
# Pipeline PyCaret + estimateur LightGBM
# ------------------------------------------------------------------------------
pipe = final_lgbm
lgbm = pipe.steps[-1][1]
# Colonnes réellement utilisées par PyCaret
feature_cols_model = list(models_setting.get_config("X_train").columns)
X_test_external = test_pycaret[feature_cols_model].copy()
print("Nombre de variables utilisées par le modèle :", len(feature_cols_model))
print("Dimensions X_test_external :", X_test_external.shape)
# ------------------------------------------------------------------------------
# Importance LightGBM
# ------------------------------------------------------------------------------
imp = pd.Series(
lgbm.feature_importances_,
index=feature_cols_model
).sort_values(ascending=False)
top20 = imp.head(20)
print(f"Top 20 variables selon l'importance LightGBM :\n{top20.to_string()}")
# ------------------------------------------------------------------------------
# SHAP global
# ------------------------------------------------------------------------------
explainer = shap.TreeExplainer(lgbm)
shap_values = explainer.shap_values(X_test_external)
if isinstance(shap_values, list):
shap_values_test = shap_values[1]
else:
shap_values_test = shap_values
# ------------------------------------------------------------------------------
# Barplot SHAP global : importance moyenne des variables
# ------------------------------------------------------------------------------
plt.figure()
shap.summary_plot(
shap_values_test,
X_test_external,
plot_type="bar",
show=False
)
plt.title("Importance globale des variables par SHAP — test externe")
plt.tight_layout()
plt.show()
# ------------------------------------------------------------------------------
# Beeswarm SHAP : sens et dispersion des effets
# ------------------------------------------------------------------------------
plt.figure()
shap.summary_plot(
shap_values_test,
X_test_external,
show=False
)
plt.tight_layout()
plt.show()
# ------------------------------------------------------------------------------
# Analyse des erreurs à partir des variables SHAP les plus importantes
# ------------------------------------------------------------------------------
vars_check = list(top20.index[:8])
summary_fpfn = (
df_err[df_err["type_erreur"].isin(["FP", "FN", "TP", "TN"])]
.groupby("type_erreur")[vars_check]
.median(numeric_only=True)
.reset_index()
)
print(summary_fpfn.to_string())
def compare_groups(df_err, g1, g2, vars_):
a = df_err[df_err.type_erreur == g1][vars_].median(numeric_only=True)
b = df_err[df_err.type_erreur == g2][vars_].median(numeric_only=True)
out = (
pd.DataFrame({
"median_" + g1: a,
"median_" + g2: b,
"diff": a - b
})
.sort_values("diff", key=np.abs, ascending=False)
)
return out
print(compare_groups(df_err, "FP", "TN", vars_check).head(15))
print(compare_groups(df_err, "FN", "TP", vars_check).head(15))
# ------------------------------------------------------------------------------
# SHAP individuel sur quelques faux positifs et faux négatifs
# ------------------------------------------------------------------------------
shap_test_external = shap_values_test
base_value = explainer.expected_value
if isinstance(base_value, list):
base_value = base_value[1]
idx_fp = np.where(df_err["type_erreur"] == "FP")[0][:3]
idx_fn = np.where(df_err["type_erreur"] == "FN")[0][:3]
for tag, idx in [("FP", idx_fp), ("FN", idx_fn)]:
for k, pos in enumerate(idx):
print(
f"\n--- {tag} position {pos} "
f"proba {float(df_err.iloc[pos]['p']):.4f} "
f"y_true {int(df_err.iloc[pos]['y_true'])}"
)
plt.close("all")
plt.figure(figsize=(14, 6))
shap.plots.bar(
shap.Explanation(
values=shap_test_external[pos],
base_values=base_value,
data=X_test_external.iloc[pos],
feature_names=X_test_external.columns
),
max_display=15,
show=False
)
plt.gcf().subplots_adjust(left=0.45)
plt.tight_layout()
plt.show()
# ==============================================================================
# 09. BEESWARM SHAP — NOMS DE VARIABLES PROPRES
# ==============================================================================
nice_labels = {
"ews_val_num": "EWS",
"adm_o2_valeur_oui": "Oxygène à l’admission (oui)",
"adm_o2_valeur_non": "Oxygène à l’admission (non)",
"crp_val_num": "CRP",
"pad_val_num": "PAD",
"pas_val_num": "PAS",
"glucose_val_num": "Glucose",
"spo2_val_num": "SpO₂",
"creatinine_enzymatique_val_num": "Créatinine",
"fc_val_num": "Fréquence cardiaque",
"temp_val_num": "Température",
"potassium_val_num": "Potassium",
"uree_val_num": "Urée",
"troponine_i_us_val_num": "Troponine I",
"glasgow_val_num": "Score de Glasgow"
}
X_test_external_clean = X_test_external.copy()
X_test_external_clean = X_test_external_clean.rename(columns=nice_labels)
plt.figure(figsize=(9, 6.5))
shap.summary_plot(
shap_values_test,
X_test_external_clean,
show=False,
max_display=15
)
plt.title(
"Distribution des valeurs SHAP",
fontsize=16,
fontweight="bold"
)
plt.xlabel(
"Valeur SHAP (impact sur la prédiction du modèle)",
fontsize=13,
fontweight="bold"
)
plt.tight_layout()
plt.show()