Choisir un modèle pour les réadmissions à l'hôpital

Étude de cas 5
120 min
Facultative · 1 étude de cas requise pour le certificat

Étude de cas 5 sur 7 · Apprentissage automatique : du problème à un modèle fiableChoisir un modèle pour les réadmissions à l'hôpital

Aperçu du cours

1. Comprendre l'apprentissage automatique

2. Cadrer le problème de prédiction

3. Évaluer les modèles supervisés

4. Créer les premiers modèles

5. Essayer d'autres familles de modèles

6. Régler et finaliser le modèle

7. Comprendre et choisir le modèle

8. Préparer la passation du modèle

9. Études de cas

Ce que vous saurez faire

Traitez un vrai problème de prédiction, du premier examen des données jusqu’à un modèle utilisable par une équipe hospitalière. Vous comparerez quatre familles de modèles selon un même protocole équitable. Vous déciderez, preuves à l’appui, si le meilleur score justifie sa complexité.

Étude de casPièce de portfolio

Jeu de données: Diabetes 130-US hospitals, 1999–2008 (UCI)
Licence: CC BY 4.0
Durée: Environ 2 h
S’exécute dans: Colab, Kaggle ou VS Code
Compétences utilisées: Préparation des données
Séparation par groupe de patients
Quatre familles de modèles
Intervalles de confiance par bootstrap
Seuil de capacité
Sous-groupes pour vérifier l'équité

Ouvrir dans Colab Ouvrir dans Kaggle Voir sur GitHub

Les carnets sont en anglais. Colab et Kaggle ouvrent votre propre copie. Pour utiliser VS Code, téléchargez le carnet depuis GitHub. Vous êtes bloqué ? Consultez le carnet de solution.

Le carnet (notebook) Jupyter exécute tous les blocs de code de cette page, avec cinq tâches et des vérifications. Ses cellules de préparation téléchargent les données (3 Mo) et les vérifient. Il ne nécessite donc aucune installation locale. Les cellules des modèles prennent environ une minute sur un portable doté de deux cœurs de processeur.

Avant de commencerCe qu’il faut savoir et avoir sous la mainEnviron 120 minutes

Ce que vous devez savoir

Cette étude de cas utilise ce que vous avez appris dans les modules 1 à 8. Voici les idées les plus utilisées :

Les données d’entraînement, de validation et de test, et la validation croisée, dans Données d’entraînement, de validation et de test.
La précision moyenne et la précision selon une capacité (precision@k), dans Le compromis entre précision et rappel.
Un pipeline avec encodage à chaud et mise à l’échelle, dans Prétraiter avec un pipeline.
Le journal d’expériences, dans L’apprentissage automatique est itératif : le journal d’expériences.
Les arbres de décision, les forêts aléatoires et le boosting (renforcement), et le coût de la complexité, dans Les arbres de décision et Les méthodes d’ensemble : le coût de la complexité.
La recherche par grille et la règle d’une erreur-type, dans Régler les hyperparamètres : la règle d’une erreur-type.
Les sous-groupes, et les situations où un modèle simple est préférable, dans L’analyse par sous-groupes et Quand un modèle simple est meilleur.
Un intervalle par bootstrap (rééchantillonnage avec remise), dans Estimer l’incertitude : le bootstrap, du cours Mathématiques et statistiques essentielles pour l’IA.

Ce qu’il vous faut sur votre ordinateur (uniquement pour l’exécution locale; le carnet n’en a pas besoin)

Un environnement Python avec pandas 3.0.6 et scikit-learn 1.9.1, comme le .venv de votre dossier de projet décrit dans Définir une référence de base : préparer votre projet. Cette étude de cas n’utilise pas les fichiers de Larkfield.
Un nouveau dossier vide pour cette étude de cas, par exemple ~/projects/readmissions, contenant le carnet. La cellule de préparation du carnet télécharge les données dans readmissions/data/. Pour recommencer, supprimez ce dossier et exécutez à nouveau la cellule de préparation.

Testé sur macOS avec Python 3.14.6, pandas 3.0.6 et scikit-learn 1.9.1, dans Jupyter. Windows, Colab et Kaggle n’ont pas encore été testés.

La situation

Nadia Haddad dirige l’équipe de transition des soins d’un groupe hospitalier. Son personnel infirmier appelle les patients dans les premiers jours suivant leur sortie de l’hôpital. Pendant l’appel, une infirmière vérifie les médicaments du patient, réserve une visite de suivi et demande qui l’aide à la maison. Ces appels peuvent prévenir un problème avant que le patient doive revenir.

L’équipe ne peut pas appeler tout le monde. Avec quatre infirmières, Nadia peut appeler environ un patient sorti de l’hôpital sur dix. Aujourd’hui, les infirmières utilisent une règle simple : appeler les patients ayant le plus de séjours hospitaliers dans la dernière année. Nadia vous pose une question simple : « Un modèle peut-il mieux choisir qui nous appelons? »

Felix Brandt, le directeur de la qualité, ajoute sa question : « Et si oui, quel modèle? J’entends dire que les méthodes les plus récentes gagnent toutes les compétitions. J’entends aussi qu’il faut pouvoir expliquer ce qu’on utilise. Montrez-moi les preuves. »

Chaque type d’erreur a un coût :

Un appel à un patient qui ne reviendrait pas prend environ 20 minutes du temps d’une infirmière. Il ne nuit pas au patient.
Un patient qui revient sans avoir reçu d’appel représente l’erreur coûteuse. Une réadmission est mauvaise pour le patient et coûte un lit et de l’argent à l’hôpital.

Comme le nombre d’appels est fixe, la question n’est pas « Ce patient est-il à risque? ». C’est plutôt « Quels sont les 10 % les plus à risque? ». C’est un problème de classement, comme la capacité de 20 % de Larkfield.

Ce que vous allez créer

Vous créerez un modèle qui classe une liste de patients sortants selon leur risque de revenir à l’hôpital dans les 30 jours. L’équipe appelle les 10 % en tête. Vous comparerez quatre familles de modèles connues : régression logistique, arbre de décision, forêt aléatoire et boosting par gradient. Vous les comparerez selon un même protocole : même séparation, mêmes plis, même métrique. Vous consignerez chaque essai. Ensuite, vous utiliserez un intervalle de confiance et la règle d’une erreur-type pour décider si le meilleur modèle justifie sa complexité. Vous terminerez par des vérifications par sous-groupes, un test final et des recommandations pour Nadia et Felix. Vous aurez aussi un dépôt à montrer à un employeur.

Les données sont réelles : 101 766 séjours de patients diabétiques, dans 130 hôpitaux américains entre 1999 et 2008. Le groupe hospitalier de Nadia et son équipe sont fictifs, créés pour cette étude de cas. Une vraie équipe utiliserait ses propres données récentes, mais suivrait les mêmes étapes.

Les données

Les données proviennent du jeu de données Diabetes 130-US hospitals for years 1999–2008 du UCI Machine Learning Repository (doi:10.24432/C5230J). Ses auteurs sont John Clore, Krzysztof Cios, Jon DeShazo et Beata Strack. Il est publié sous la licence Creative Commons Attribution 4.0, qui autorise la réutilisation avec attribution et mention des changements. L’étude qui le décrit est celle de Strack et ses collègues (2014), « Impact of HbA1c measurement on hospital readmission rates », dans BioMed Research International (doi:10.1155/2014/781670).

Nextia Learning héberge une copie dans le dépôt des laboratoires, avec une fiche du jeu de données qui répertorie chaque changement. Il y a un seul changement : dans deux colonnes, le texte None (« le test n’a pas été effectué ») devient une cellule vide. La source a dépersonnalisé les données : elles ne contiennent ni noms ni dates. Ce sont tout de même des données de santé sur de vraies personnes. Nous ne tentons donc jamais de découvrir leur identité.

Les colonnes que vous utiliserez

Une ligne représente une rencontre : un séjour hospitalier d’un patient, d’une durée de 1 à 14 jours. Le fichier compte 50 colonnes. Vous utiliserez celles-ci :

Colonne	Signification
`encounter_id`, `patient_nbr`	Le numéro du séjour et celui du patient. Un patient peut avoir plusieurs séjours.
`readmitted`	`<30` : retour à l’hôpital dans les 30 jours; `>30` : retour plus tard; `NO` (non) : aucun retour enregistré.
`age`	Tranche d’âge de 10 ans, comme `[70-80)`.
`admission_type_id`, `admission_source_id`, `discharge_disposition_id`	Codes indiquant comment le patient est arrivé et où il est allé après le séjour. Un fichier distinct les traduit.
`time_in_hospital`	Nombre de jours à l’hôpital.
`num_lab_procedures`, `num_procedures`, `num_medications`, `number_diagnoses`	Nombres relevés pendant le séjour.
`number_outpatient`, `number_emergency`, `number_inpatient`	Visites dans l’année avant le séjour : consultations externes, visites à l’urgence, séjours hospitaliers.
`medical_specialty`	La spécialité du médecin qui a admis le patient.
`diag_1`	Le diagnostic principal, sous forme de code ICD-9 (la liste internationale des maladies utilisée alors par les hôpitaux américains).
`A1Cresult`, `max_glu_serum`	Résultats de deux tests de glycémie, ou cellule vide si aucune mesure n’a été prise.
`insulin`, `change`, `diabetesMed`	Changements à l’insuline, tout changement aux médicaments contre le diabète et administration ou non d’un médicament contre le diabète.
`race`, `gender`, `payer_code`	Colonnes sensibles : utilisées uniquement pour vérifier le modèle, jamais comme entrées (voir Cadrer le problème).

Premières vérifications

Chargez le fichier. La source représente les valeurs manquantes par ?. Indiquez donc à pandas de lire ? comme une valeur manquante :

raw = pd.read_csv("data/diabetes_130_hospitals.csv.gz", na_values="?", low_memory=False)
print(raw.shape)
print("encounters:", raw["encounter_id"].nunique(), "| patients:", raw["patient_nbr"].nunique())
print("duplicate rows:", raw.duplicated().sum())

Le résultat :

(101766, 50)
encounters: 101766 | patients: 71518
duplicate rows: 0

Ce que cela nous apprend. Chaque ligne représente un séjour différent : 101 766 lignes, 101 766 numéros de rencontre et aucune ligne en double. Mais il n’y a que 71 518 patients. Environ 30 000 lignes représentent donc un deuxième, troisième séjour ou davantage des mêmes personnes. Gardez cela en tête. Ce fait détermine comment nous séparons les données.

Ensuite, voici la proportion de valeurs manquantes dans chaque colonne qui en contient :

missing = raw.isna().mean()
missing[missing > 0].sort_values(ascending=False).round(3)

Le résultat :

weight               0.969
max_glu_serum        0.947
A1Cresult            0.833
medical_specialty    0.491
payer_code           0.396
race                 0.022
diag_3               0.014
diag_2               0.004
diag_1               0.000
dtype: float64

Ce que cela nous apprend. Cette liste cache trois types de valeurs « manquantes ». Chacun demande une décision différente :

weight est absent dans 97 % des séjours. Une colonne presque toujours vide ne peut pas aider. Nous la supprimons donc.
max_glu_serum et A1Cresult ne sont pas vraiment manquants. Une cellule vide signifie « le test n’a pas été effectué », ce qui constitue une information. Strack et ses collègues ont justement étudié si le test HbA1c change les réadmissions. Nous créons une valeur distincte pour « non mesuré ».
medical_specialty est absent dans la moitié des séjours, probablement parce que personne ne l’a renseigné. Nous conservons « manquant » comme valeur pour que le modèle puisse apprendre si cela compte. diag_1 est absent dans seulement 21 lignes.

Passons à la cible :

raw["readmitted"].value_counts()

Le résultat :

readmitted
NO     54864
>30    35545
<30    11357
Name: count, dtype: int64

Ce que cela nous apprend. Seuls 11 357 séjours (11 %) sont suivis d’un retour dans les 30 jours. Nadia appelle dans les premiers jours après le séjour. Un retour après deux mois n’est donc pas ce qu’elle peut prévenir. Nous prédisons <30 et traitons >30 et NO comme « aucun retour dans les 30 jours ». La classe positive est rare. L’exactitude serait donc trompeuse. Un modèle qui dit « personne ne revient » a raison dans 89 % des cas, mais il est inutile. La leçon Les métriques de classification l’a montré.

Traduire les codes

Trois colonnes contiennent des nombres, par exemple discharge_disposition_id = 11. Leur signification figure dans IDS_mapping.csv. Ce fichier contient trois petits tableaux successifs, séparés par une ligne vide. Le bloc suivant les lit pour créer un dictionnaire par colonne :

import csv

codes, column = {}, None
with open("data/IDS_mapping.csv", newline="") as file:
    for row in csv.reader(file):
        if not row or not row[0]:
            continue                          # the empty line between two tables
        if not row[0].isdigit():
            column = row[0]                   # a header line starts a new table
            codes[column] = {}
        else:
            codes[column][int(row[0])] = row[1].strip()
print({name: len(table) for name, table in codes.items()})
print(codes["admission_type_id"])

Le résultat :

{'admission_type_id': 8, 'discharge_disposition_id': 30, 'admission_source_id': 25}
{1: 'Emergency', 2: 'Urgent', 3: 'Elective', 4: 'Newborn', 5: 'Not Available', 6: 'NULL', 7: 'Trauma Center', 8: 'Not Mapped'}

Ce que cela nous apprend. Les codes se répartissent dans trois tableaux de 8, 30 et 25 valeurs. Regardez les types d’admission : trois des huit codes signifient tous « nous ne savons pas ». Ce sont Not Available (non disponible), NULL (valeur nulle) et Not Mapped (sans correspondance). Un modèle les traiterait comme trois éléments différents. Nous les regrouperons sous une seule valeur, unknown (inconnu).

Retirer les séjours sans possibilité de retour

Certains séjours se terminent par le décès du patient ou son transfert dans un établissement de soins de fin de vie. Ces patients ne peuvent pas être réadmis, et l’équipe de Nadia ne les appellerait pas. Si nous les gardons, le modèle apprend que « décédé à l’hôpital » prédit « aucun retour ». C’est vrai, mais inutile. Strack et ses collègues les ont aussi retirés.

ENDED = [11, 13, 14, 19, 20, 21]   # the patient died or went to a hospice
for code in ENDED:
    print(code, codes["discharge_disposition_id"][code])
data = raw[~raw["discharge_disposition_id"].isin(ENDED)].copy()
print(len(raw) - len(data), "encounters removed,", len(data), "left")

Le résultat :

11 Expired
13 Hospice / home
14 Hospice / medical facility
19 Expired at home. Medicaid only, hospice.
20 Expired in a medical facility. Medicaid only, hospice.
21 Expired, place unknown. Medicaid only, hospice.
2423 encounters removed, 99343 left

Ce que cela nous apprend. 2 423 séjours se sont terminés par un décès ou des soins de fin de vie. Nous en conservons 99 343. Nous affichons les six codes du tableau plutôt que de nous fier à une liste mémorisée. Tout le monde peut ainsi vérifier cette étape.

Créez maintenant la cible : 1 pour un retour dans les 30 jours :

TARGET = "readmitted_30"
data[TARGET] = (data["readmitted"] == "<30").astype(int)
print(data[TARGET].sum(), "encounters were followed by a readmission within 30 days:",
      f"{data[TARGET].mean():.1%}")

Le résultat :

11314 encounters were followed by a readmission within 30 days: 11.4%

Ce que cela nous apprend. Le taux est de 11,4 %. Si l’équipe appelait 10 % des patients au hasard, environ 11 appels sur 100 joindraient un patient qui revient. Chaque modèle doit faire mieux.

Cadrer le problème

Avant de créer un modèle, notez exactement ce que nous prédisons, et à quel moment. Ces choix protègent le projet des erreurs les plus courantes.

Unité : une rencontre, au moment où le patient quitte l’hôpital. Un patient ayant trois séjours représente trois lignes. Chaque séjour reçoit son propre score.
Cible : readmitted_30, un retour avec hospitalisation dans les 30 jours.
Moment de la prédiction : le jour de la sortie. Une entrée est permise uniquement si l’équipe la connaît ce jour-là. Tout ce qui concerne le séjour est connu : durée, tests, médicaments et destination du patient. Les visites dans l’année avant le séjour sont aussi connues.
Indisponible au moment de la prédiction (fuites) : readmitted lui-même, et tout ce qui se produit après la sortie. Un exemple tentant est « combien de séjours ce patient a dans le fichier ». Le fichier contient aussi des séjours après celui-ci. Ce nombre utilise donc le futur. Le cas d’échec montre son effet.
Jamais des entrées : encounter_id et patient_nbr identifient un séjour ou une personne. Ce ne sont pas des faits sur la santé du patient.

Colonnes sensibles

Trois colonnes décrivent qui est le patient, pas son état de santé : race, gender et payer_code. Cette dernière indique le type d’assurance, qui révèle aussi le revenu. Nous les gardons hors des entrées. Un modèle classant les patients par race ou assurance accorderait moins d’appels à certains groupes pour des raisons indéfendables pour l’équipe. Il existe un exemple réel de ce préjudice. En 2019, Obermeyer et ses collègues ont publié dans Science une étude sur un algorithme de santé largement utilisé aux États-Unis. Il attribuait aux patients noirs des scores de risque inférieurs à ceux de patients blancs aussi malades. La raison : il prédisait les coûts de santé, et moins d’argent avait été dépensé pour les patients noirs.

Exclure ces colonnes ne suffit pas, car d’autres entrées peuvent transmettre la même information. Nous les conservons donc pour les vérifications par sous-groupes. Après l’entraînement, nous vérifions si le modèle fonctionne aussi bien et attribue une part équitable des appels à chaque groupe. L’âge est différent. C’est un facteur de risque médical, et les infirmières en tiennent compte. Il fait donc partie des entrées. Nous vérifions tout de même le modèle pour chaque tranche d’âge.

Préparer les variables

Le diagnostic principal compte des centaines de codes ICD-9. Strack et ses collègues les ont répartis en neuf groupes. Par exemple : maladies circulatoires (codes 390 à 459, et 785) et diabète (250.xx). Nous utilisons les mêmes groupes :

GROUPS = [("circulatory", 390, 460, "785"), ("respiratory", 460, 520, "786"),
          ("digestive", 520, 580, "787"), ("injury", 800, 1000, None),
          ("musculoskeletal", 710, 740, None), ("genitourinary", 580, 630, "788"),
          ("neoplasms", 140, 240, None)]


def diagnosis_group(code):
    """The group of the main diagnosis (an ICD-9 code), as in Strack et al. (2014)."""
    if pd.isna(code):
        return "missing"
    if code.startswith(("V", "E")):
        return "other"                        # visits for other reasons, and external causes
    if code.startswith("250"):
        return "diabetes"
    number = float(code)
    for name, low, high, extra in GROUPS:
        if low <= number < high or (extra and code.startswith(extra)):
            return name
    return "other"


print([diagnosis_group(code) for code in ["428", "250.83", "786", "V57", "996"]])

Le résultat :

['circulatory', 'diabetes', 'respiratory', 'other', 'injury']

Ce que cela nous apprend. La fonction donne le groupe attendu pour cinq codes de test. 428 (insuffisance cardiaque) correspond aux maladies circulatoires, et 250,83 au diabète. V57 (visite de réadaptation) correspond au groupe « autre », et 996 (complication liée à un dispositif) aux blessures. Tester une fonction sur quelques valeurs connues prend une ligne et détecte beaucoup d’erreurs.

Traduisez maintenant les codes, regroupez les codes « inconnu » et créez les autres variables :

UNKNOWN = ["NULL", "Not Available", "Not Mapped", "Unknown/Invalid"]
for column, name in [("admission_type_id", "admission_type"), ("admission_source_id", "admission_source"),
                     ("discharge_disposition_id", "discharge_to")]:
    data[name] = data[column].map(codes[column]).replace(UNKNOWN, "unknown").fillna("unknown")
data["diagnosis_group"] = data["diag_1"].map(diagnosis_group)
data["age_years"] = data["age"].str.extract(r"(\d+)")[0].astype(int) + 5   # "[70-80)" becomes 75
data["specialty"] = data["medical_specialty"].fillna("missing")
tests = ["A1Cresult", "max_glu_serum"]
data[tests] = data[tests].fillna("not measured")
data[["admission_type", "discharge_to", "diagnosis_group", "age_years", "A1Cresult"]].head()

Le résultat :

  admission_type        discharge_to diagnosis_group  age_years     A1Cresult
0        unknown             unknown        diabetes          5  not measured
1      Emergency  Discharged to home           other         15  not measured
2      Emergency  Discharged to home           other         25  not measured
3      Emergency  Discharged to home           other         35  not measured
4      Emergency  Discharged to home       neoplasms         45  not measured

Ce que cela nous apprend. Les codes sont maintenant des mots qu’une infirmière peut lire. La tranche d’âge est devenue un nombre, son point milieu. « Plus âgé » signifie donc « plus grand ». .fillna("unknown") prend en charge les codes totalement absents du tableau.

Enfin, notez les listes de variables. Les colonnes sensibles ont leur propre liste. Le code vérifie qu’aucune ne fait partie des entrées :

CATEGORIES = ["admission_type", "admission_source", "discharge_to", "specialty", "diagnosis_group",
              "A1Cresult", "max_glu_serum", "insulin", "change", "diabetesMed"]
NUMBERS = ["age_years", "time_in_hospital", "num_lab_procedures", "num_procedures", "num_medications",
           "number_outpatient", "number_emergency", "number_inpatient", "number_diagnoses"]
FEATURES = CATEGORIES + NUMBERS
SENSITIVE = ["race", "gender", "payer_code"]   # for slice checks only, never inputs
print(len(FEATURES), "features; sensitive columns among them:", set(SENSITIVE) & set(FEATURES))

Le résultat :

19 features; sensitive columns among them: set()

Ce que cela nous apprend. Le modèle verra 19 variables : 10 catégorielles et 9 numériques. set() est un ensemble vide : aucune colonne sensible ne fait partie des entrées. Pour limiter la taille de l’étude de cas, nous excluons les deuxième et troisième diagnostics, ainsi que 22 colonnes de médicaments individuels. La liste Aller plus loin dans Faites-en un projet de portfolio suggère de les ajouter.

Explorer les besoins de la décision

La règle simple de Nadia utilise un seul fait : les séjours hospitaliers dans la dernière année. Avant de créer quoi que ce soit, vérifiez la force de ce fait. Une règle déjà efficace impose une exigence élevée au modèle.

prior_stays = data["number_inpatient"].clip(upper=5)   # 5 means "5 or more"
data.groupby(prior_stays)[TARGET].agg(["size", "mean"]).round(3)

Le résultat :

                   size   mean
number_inpatient
0                 66245  0.086
1                 18984  0.133
2                  7300  0.179
3                  3271  0.210
4                  1574  0.241
5                  1969  0.371

La proportion de séjours suivis d'une réadmission dans les 30 jours, selon le nombre de séjours hospitaliers dans l'année précédente. Le taux augmente avec chaque séjour antérieur.

Description textuelle du graphique

Un graphique à colonnes. L’axe horizontal indique le nombre de séjours hospitaliers dans l’année précédente : 0, 1, 2, 3, 4 et 5 ou plus. L’axe vertical indique la proportion de réadmissions dans les 30 jours : 8,6 % (66 245 séjours), 13,3 % (18 984), 17,9 % (7 300), 21,0 % (3 271), 24,1 % (1 574) et 37,1 % (1 969).

Ce que cela nous apprend. La règle simple est bonne : le taux passe de 8,6 % sans séjour antérieur à 37,1 % avec cinq ou plus. Mais regardez les effectifs. Deux séjours sur trois (66 245) n’avaient aucun séjour antérieur. Leur taux est de 8,6 %, pas zéro. Ces patients représentent une grande part des retours, et la règle ne peut pas les distinguer. Un modèle peut seulement battre la règle si d’autres faits aident pour ces patients.

La destination du patient après le séjour est le deuxième fait connu d’une infirmière au moment de la sortie :

by_discharge = data.groupby("discharge_to")[TARGET].agg(["size", "mean"])
by_discharge[by_discharge["size"] >= 1000].sort_values("mean", ascending=False).round(3)

Le résultat :

                                                     size   mean
discharge_to
Discharged/transferred to another rehab fac inc...   1993  0.277
Discharged/transferred to another type of inpat...   1184  0.209
Discharged/transferred to another short term ho...   2128  0.161
Discharged/transferred to SNF                       13954  0.147
Discharged/transferred to home with home health...  12902  0.127
unknown                                              4680  0.118
Discharged to home                                  60234  0.093

Ce que cela nous apprend. Les patients transférés dans une unité de réadaptation reviennent dans 27,7 % des cas. C’est trois fois plus que ceux qui rentrent chez eux (9,3 %). Les patients transférés dans un établissement de soins infirmiers spécialisés (SNF) ou recevant des soins à domicile se situent entre les deux. La règle n’utilise pas ce fait. Un modèle combinant les deux devrait donc mieux classer les patients. La suite de l’étude de cas mesure s’il y parvient, et de combien.

Répartir par patient

Pourquoi ne pas répartir par ligne

Dans ces données, 16 341 des 69 990 patients ont plus d’un séjour :

encounters = data["patient_nbr"].value_counts()
print(encounters.size, "patients;", (encounters > 1).sum(), "of them have more than one encounter")
print("their encounters:", encounters[encounters > 1].sum(), "of", len(data), "rows")
print("most encounters for one patient:", encounters.max())

La sortie :

69990 patients; 16341 of them have more than one encounter
their encounters: 45694 of 99343 rows
most encounters for one patient: 40

Ce que cela nous apprend. Près de la moitié des lignes (45 694) appartiennent à des patients qui apparaissent plusieurs fois. Un patient compte 40 séjours. Les lignes ne sont pas indépendantes. Deux séjours d’une même personne partagent sa maladie, son milieu de vie et ses habitudes.

Reprenons l’exemple d’un examen. Si les questions reprennent les exercices avec de nouveaux nombres, l’élève qui a mémorisé les réponses semble meilleur qu’il l’est. Une répartition par ligne fait la même chose. Elle peut placer le séjour de mars d’un patient dans les données d’entraînement, et celui de mai dans les données de test. Le modèle a alors « vu » ce patient. Le score de test répond à une question plus facile que celle de Nadia. Son modèle rencontrera beaucoup de patients pour la première fois. Nous répartissons donc par patient : chaque patient va, avec tous ses séjours, dans une seule partie.

La répartition

La répartition attribue une partie à chaque patient, plutôt qu’à chaque ligne : 60 % pour l’entraînement, 20 % pour la validation, 20 % pour le test. Un générateur aléatoire avec une graine fixe rend cette répartition reproductible.

patients = data["patient_nbr"].unique()
rng = np.random.default_rng(0)
part_of = dict(zip(patients, rng.choice(["train", "valid", "test"], size=len(patients), p=[0.6, 0.2, 0.2])))
data["part"] = data["patient_nbr"].map(part_of)
train, valid, test = (data[data["part"] == name].copy() for name in ["train", "valid", "test"])
for name, part in [("train", train), ("valid", valid), ("test", test)]:
    print(f"{name}: {len(part)} encounters, {part['patient_nbr'].nunique()} patients,"
          f" readmitted {part[TARGET].mean():.3f}")
print("patients in more than one part:", (data.groupby("patient_nbr")["part"].nunique() > 1).sum())

La sortie :

train: 59556 encounters, 41999 patients, readmitted 0.114
valid: 20179 encounters, 14139 patients, readmitted 0.113
test: 19608 encounters, 13852 patients, readmitted 0.116
patients in more than one part: 0

Ce que cela nous apprend. Les trois parties contiennent environ 60 %, 20 % et 20 % des lignes. Leur taux de réadmission est presque identique (11,3 % à 11,6 %). Aucun patient ne figure dans deux parties. Désormais, les données de test restent scellées : nous ne les consultons pas avant le test final.

Pourquoi ne pas répartir selon le temps, comme chez Larkfield? Un déploiement réel prédit toujours l’avenir. Une répartition temporelle serait donc préférable. Mais ce fichier ne contient aucune date. Les numéros des rencontres augmentent probablement avec le temps, mais la source ne le confirme pas. Nous ne nous basons pas sur une supposition. Nous signalons cette limite.

Ce qu’une répartition par ligne donnerait

Pour voir le problème, répartissez les données d’entraînement par ligne, comme le fait train_test_split par défaut. Les données de test restent scellées :

from sklearn.model_selection import train_test_split

fit_rows, check_rows = train_test_split(train, test_size=0.2, random_state=0)
seen = check_rows["patient_nbr"].isin(fit_rows["patient_nbr"])
print("patients in both parts:", check_rows.loc[seen, "patient_nbr"].nunique())
print(f"check rows of a patient that the model has seen: {seen.sum()} of {len(check_rows)} ({seen.mean():.0%})")

La sortie :

patients in both parts: 4046
check rows of a patient that the model has seen: 4915 of 11912 (41%)

Ce que cela nous apprend. Avec une répartition par ligne, 41 % des lignes que nous « vérifions » appartiennent à 4 046 patients ayant servi à l’entraînement. Tout score sur ces données de vérification mélange deux questions. « À quel point le modèle connaît-il ces patients? » et « À quel point prédit-il correctement pour de nouveaux patients? ». Seule la deuxième compte pour Nadia. Plus loin, vous mesurerez à quel point cela change le score sur ces données.

Références de base et protocole

Un protocole pour tous les modèles

« Comparer équitablement » signifie utiliser les mêmes données, les mêmes plis, la même métrique et les mêmes caractéristiques. Nous fixons ces quatre éléments maintenant, avant tout modèle, pour n’en avantager aucun :

Données : uniquement les données d’entraînement. La validation confirme le choix; les données de test restent scellées.
Plis : 5 plis, par patient. StratifiedGroupKFold de scikit-learn s’en charge. « Group » signifie que toutes les lignes d’un patient vont dans le même pli. « Stratified » signifie que chaque pli reçoit environ le même taux de réadmission.
Métrique : la précision moyenne (AP), comme dans Le compromis entre précision et rappel. Elle résume la capacité du modèle à placer les rares cas positifs en tête, sur l’ensemble des seuils. C’est la bonne métrique pour choisir un modèle. Pour la décision, nous ajouterons plus tard la précision parmi les 10 % les mieux classés.
Caractéristiques et préparation : les 19 caractéristiques, dans un même pipeline pour tous les modèles.

Le protocole. Les 99 343 rencontres de 69 990 patients, après exclusion des décès et des soins palliatifs, passent à Répartir par patient. Cette étape place chaque patient dans une seule partie. Entraînement, 59 556 rencontres, passe à 5 plis par patient : régler chaque famille, comparer, règle d'une erreur-type. Le résultat passe à Confirmer sur les données de validation. Cette étape reçoit aussi Validation, 20 179 rencontres : estimer l'écart par bootstrap et fixer le seuil de 10 %. Confirmer sur les données de validation et Test, 19 608 rencontres, passent à Tester une seule fois, avec le modèle choisi. — Le protocole de cette étude de cas. Chaque modèle est choisi sur les mêmes 5 plis d'entraînement, répartis par patient, puis confirmé sur les données de validation. Seul le modèle choisi rencontre les données de test, une seule fois.

Le bloc suivant crée les plis et le pipeline. OneHotEncoder transforme chaque catégorie en colonnes de 0 et de 1, comme dans Prétraiter avec un pipeline. Deux réglages sont nouveaux. min_frequency=100 regroupe les catégories de moins de 100 lignes dans une colonne « peu fréquente ». Une spécialité rare avec 12 séjours n’enseigne que du bruit. handle_unknown="infrequent_if_exist" envoie dans cette même colonne toute catégorie qui apparaît seulement dans les nouvelles données, plutôt que de provoquer une erreur.

from sklearn.compose import ColumnTransformer
from sklearn.model_selection import StratifiedGroupKFold
from sklearn.pipeline import Pipeline
from sklearn.preprocessing import OneHotEncoder, StandardScaler

folds = list(StratifiedGroupKFold(5, shuffle=True, random_state=0)
             .split(train, train[TARGET], groups=train["patient_nbr"]))


def build_pipeline(model, numbers=NUMBERS):
    """The same preparation for every model: one-hot categories (rare ones grouped) and scaled numbers."""
    prepare = ColumnTransformer([
        ("categories", OneHotEncoder(handle_unknown="infrequent_if_exist", min_frequency=100,
                                     sparse_output=False), CATEGORIES),
        ("numbers", StandardScaler(), numbers)])
    return Pipeline([("prepare", prepare), ("model", model)])

Vérifiez que les plis séparent bien les patients :

for number, (fit, check) in enumerate(folds, 1):
    shared = set(train["patient_nbr"].iloc[fit]) & set(train["patient_nbr"].iloc[check])
    print(f"fold {number}: {len(check)} check rows, readmitted {train[TARGET].iloc[check].mean():.3f},"
          f" patients also in the fit rows: {len(shared)}")

La sortie :

fold 1: 11912 check rows, readmitted 0.114, patients also in the fit rows: 0
fold 2: 11910 check rows, readmitted 0.114, patients also in the fit rows: 0
fold 3: 11912 check rows, readmitted 0.113, patients also in the fit rows: 0
fold 4: 11911 check rows, readmitted 0.114, patients also in the fit rows: 0
fold 5: 11911 check rows, readmitted 0.114, patients also in the fit rows: 0

Ce que cela nous apprend. Chaque pli vérifie environ 11 900 lignes avec le même taux (11,3 % à 11,4 %). Aucun patient ne figure des deux côtés. Tous les modèles ci-dessous utilisent exactement ces folds (plis). La liste est fixe, donc tous les modèles sont évalués sur les mêmes lignes.

Références de base

Le journal des expériences reprend le tableau de L’apprentissage automatique est itératif. Il contient les colonnes de ce cours : identifiant, date, changement, AP en validation croisée et son écart-type, AP de validation, décision et raison. La fonction utilitaire log_experiment() du carnet (notebook) Jupyter écrit une ligne dans experiment_log.csv. Elle conserve aussi les scores des cinq plis dans FOLD_SCORES, pour comparer les modèles pli par pli.

Deux références de base viennent d’abord. Le plancher donne le même score à chaque patient. Il ne classe personne, donc son AP égale le taux de réadmission. La règle des séjours antérieurs traduit la règle pratique de Nadia en score. Le score est le nombre de séjours dans la dernière année; les visites aux urgences départagent les égalités. Une règle n’a pas besoin d’entraînement. Nous la mesurons donc seulement sur chaque pli de vérification.

from sklearn.metrics import average_precision_score


def stays_rule(df):
    """The rule of thumb: most hospital stays in the past year first; emergency visits break ties."""
    return df["number_inpatient"] + df["number_emergency"].clip(upper=99) / 100


def rule_folds(score):
    """The AP of a fixed score (nothing to train) on each of the 5 check folds."""
    return np.array([average_precision_score(train[TARGET].iloc[check], score.iloc[check])
                     for _, check in folds])


log_experiment("E01", "Baseline: the same score for every patient", rule_folds(train["number_inpatient"] * 0.0),
               decision="baseline", reason="The floor: its AP is the readmission rate.")
log_experiment("E02", "Baseline: the prior-stays rule", rule_folds(stays_rule(train)),
               decision="baseline", reason="The rule of thumb that a model must beat.")
show_log()

La sortie :

 id                                     change  cv_ap  cv_std valid_ap decision
E01 Baseline: the same score for every patient  0.114   0.000          baseline
E02             Baseline: the prior-stays rule  0.182   0.009          baseline

Ce que cela nous apprend. Le plancher est de 0,114, soit le taux de réadmission. La règle atteint 0,182, soit 60 % de plus que le plancher. C’est le niveau à dépasser : un modèle qui ne dépasse pas 0,182 ne vaut pas la peine d’être construit.

Quatre familles de modèles sous un même protocole

Personne ne peut dire à l’avance quelle famille convient le mieux à ces données. L’apprentissage automatique est itératif l’expliquait. Chaque famille suppose une forme différente de relation. Nous en essayons donc quatre, de la plus simple à la plus complexe :

La régression logistique suppose que chaque fait ajoute ou retire du risque de façon indépendante. Elle est rapide et compacte, et chaque score se décompose en contributions exactes.
Un arbre de décision pose quelques questions dans un ordre précis. Il peut trouver des combinaisons (« nombreux séjours et orientation vers la réadaptation »). Un petit arbre peut être imprimé pour le personnel infirmier. Mais un seul arbre est peu précis et instable.
Une forêt aléatoire fait la moyenne de nombreux arbres profonds. Chacun est entraîné sur un échantillon aléatoire de lignes et de caractéristiques. La moyenne réduit la variance des arbres individuels.
Le gradient boosting (renforcement par gradient) ajoute de petits arbres l’un après l’autre. Chacun corrige les erreurs des précédents. Il offre souvent le meilleur classement sur les données tabulaires.

Régler chaque famille avec modération

Chaque famille a un ou deux réglages particulièrement importants. Nous réglons un paramètre par famille avec une petite grille, sur les mêmes plis, comme dans Régler les hyperparamètres. Nous gardons volontairement les grilles petites. Chaque réglage supplémentaire essayé sur les mêmes plis augmente un peu le risque de s’ajuster à leur bruit. La fonction utilitaire tune() applique la règle d’une erreur-type au sein de chaque famille. Elle retient le réglage le plus simple parmi ceux dont l’AP moyenne est à moins d’une erreur-type du meilleur. Chaque grille présente son réglage le plus simple en premier.

from sklearn.model_selection import GridSearchCV


def tune(pipeline, grid):
    """Score each setting of `grid` on the 5 patient folds. Keep the simplest setting whose mean is
    within one standard error of the best (each grid lists the simplest setting first)."""
    search = GridSearchCV(pipeline, grid, cv=folds, scoring="average_precision", refit=False)
    results = pd.DataFrame(search.fit(train[FEATURES], train[TARGET]).cv_results_)
    results["se"] = results["std_test_score"] / np.sqrt(len(folds))
    best = results["mean_test_score"].idxmax()
    bar = results.loc[best, "mean_test_score"] - results.loc[best, "se"]
    chosen = results.index[results["mean_test_score"] >= bar][0]
    print(results[["params", "mean_test_score", "std_test_score"]].round(4).to_string())
    print("one-SE bar:", round(bar, 4), "| chosen:", results.loc[chosen, "params"])
    pipeline.set_params(**results.loc[chosen, "params"])
    return np.array([results.loc[chosen, f"split{i}_test_score"] for i in range(len(folds))])

Régression logistique

Le réglage est C, qui détermine la force de la pénalité sur les grands coefficients. Un petit C signifie une forte pénalité et un modèle plus simple.

from sklearn.linear_model import LogisticRegression

logistic = build_pipeline(LogisticRegression(max_iter=1000))
scores = tune(logistic, {"model__C": [0.01, 0.1, 1.0]})
log_experiment("E03", "Logistic regression; C tuned over 0.01, 0.1, 1", scores, decision="keep",
               reason="Well above the rule; the strongest penalty is within one SE of the best.")

La sortie :

               params  mean_test_score  std_test_score
0  {'model__C': 0.01}           0.2100          0.0101
1   {'model__C': 0.1}           0.2116          0.0101
2   {'model__C': 1.0}           0.2114          0.0098
one-SE bar: 0.2071 | chosen: {'model__C': 0.01}

Ce que cela nous apprend. Les trois réglages obtiennent entre 0,2100 et 0,2116, soit un écart de 0,0016. C’est bien moins qu’une erreur-type (environ 0,0045). La règle conserve donc le plus simple, C = 0.01. Avec 60 000 lignes et seulement 19 caractéristiques, la pénalité compte très peu. La régression logistique dépasse déjà nettement la règle : 0,210 contre 0,182.

Arbre de décision

Le réglage est la profondeur maximale, avec au moins 100 séjours dans chaque feuille. Ainsi, le taux de chaque feuille repose sur assez de patients.

from sklearn.tree import DecisionTreeClassifier

tree = build_pipeline(DecisionTreeClassifier(min_samples_leaf=100, random_state=0))
scores = tune(tree, {"model__max_depth": [4, 6, 8, 10]})
log_experiment("E04", "Decision tree, min leaf 100; depth tuned over 4, 6, 8, 10", scores, decision="keep",
               reason="Above the rule, below logistic regression.")

La sortie :

                     params  mean_test_score  std_test_score
0   {'model__max_depth': 4}           0.1946          0.0071
1   {'model__max_depth': 6}           0.2020          0.0069
2   {'model__max_depth': 8}           0.2020          0.0061
3  {'model__max_depth': 10}           0.2002          0.0066
one-SE bar: 0.1989 | chosen: {'model__max_depth': 6}

Ce que cela nous apprend. Les profondeurs 6 et 8 sont à égalité à 0,2020. La profondeur 10 fait déjà moins bien : les arbres plus profonds commencent à s’ajuster au bruit. La règle conserve la profondeur 6. Même le meilleur arbre (0,202) reste sous la régression logistique (0,210). Un arbre de profondeur 6 compte jusqu’à 64 feuilles, déjà trop pour qu’une infirmière puisse le lire. Il découpe aussi chaque nombre en quelques paliers.

Forêt aléatoire

Nous limitons la forêt à 100 arbres pour que le carnet s’exécute en quelques minutes. Ajouter des arbres change très peu le score. Le réglage est le nombre minimal de séjours par feuille : de plus grandes feuilles donnent des arbres plus simples.

from sklearn.ensemble import RandomForestClassifier

forest = build_pipeline(RandomForestClassifier(n_estimators=100, n_jobs=-1, random_state=0))
scores = tune(forest, {"model__min_samples_leaf": [50, 20, 10]})
log_experiment("E05", "Random forest, 100 trees; min leaf tuned over 50, 20, 10", scores, decision="keep",
               reason="Above logistic regression on the folds.")

La sortie :

                            params  mean_test_score  std_test_score
0  {'model__min_samples_leaf': 50}           0.2163          0.0080
1  {'model__min_samples_leaf': 20}           0.2206          0.0099
2  {'model__min_samples_leaf': 10}           0.2237          0.0108
one-SE bar: 0.2188 | chosen: {'model__min_samples_leaf': 20}

Ce que cela nous apprend. Les petites feuilles obtiennent de meilleurs scores : 0,2163, 0,2206, 0,2237. Le meilleur réglage est une feuille de 10. Mais une feuille de 20 reste à moins d’une erreur-type, donc la règle la conserve. À 0,2206, la forêt est le premier modèle à dépasser la régression logistique.

Gradient boosting

Nous utilisons HistGradientBoostingClassifier de scikit-learn avec 100 arbres. early_stopping=False force la construction des 100 arbres à chaque fois. Chaque entraînement effectue donc la même quantité de travail. Le réglage est la profondeur de chaque petit arbre.

from sklearn.ensemble import HistGradientBoostingClassifier

boosting = build_pipeline(HistGradientBoostingClassifier(early_stopping=False, random_state=0))
scores = tune(boosting, {"model__max_depth": [2, 3, None]})
log_experiment("E06", "Gradient boosting, 100 trees; depth tuned over 2, 3, no limit", scores, decision="keep",
               reason="The best family on the folds, with the forest.")

La sortie :

                       params  mean_test_score  std_test_score
0     {'model__max_depth': 2}           0.2188          0.0085
1     {'model__max_depth': 3}           0.2226          0.0087
2  {'model__max_depth': None}           0.2212          0.0119
one-SE bar: 0.2188 | chosen: {'model__max_depth': 2}

Ce que cela nous apprend. La profondeur 3 est la meilleure (0,2226). La profondeur 2 atteint tout juste le niveau requis (0,2188), donc la règle conserve les arbres de profondeur 2. Chacun de ces arbres pose seulement deux questions. Le gradient boosting dépasse aussi la régression logistique.

Comparer les familles

Placez maintenant côte à côte le réglage choisi pour chaque famille, avec l’erreur-type de chaque moyenne :

compare = pd.DataFrame({id: FOLD_SCORES[id] for id in ["E02", "E03", "E04", "E05", "E06"]},
                       index=[f"fold {i}" for i in range(1, 6)])
summary = pd.DataFrame({"mean": compare.mean(), "se": compare.std(ddof=0) / np.sqrt(len(folds))})
best = summary["mean"].idxmax()
bar = summary.loc[best, "mean"] - summary.loc[best, "se"]
print(summary.round(4).T)
print("best:", best, "| one-SE bar:", round(bar, 4), "| within one SE:", summary.index[summary["mean"] >= bar].tolist())
print("boosting minus logistic, per fold:", (compare["E06"] - compare["E03"]).round(4).tolist())

La sortie :

         E02     E03     E04     E05     E06
mean  0.1824  0.2100  0.2020  0.2206  0.2188
se    0.0040  0.0045  0.0031  0.0044  0.0038
best: E05 | one-SE bar: 0.2162 | within one SE: ['E05', 'E06']
boosting minus logistic, per fold: [0.0127, 0.007, 0.0136, 0.0065, 0.004]

Précision moyenne en validation croisée de chaque expérience sur les mêmes 5 plis d'entraînement répartis par patient, avec une erreur-type de chaque côté. Les deux méthodes d'ensemble (orange) atteignent le niveau requis par la règle d'une erreur-type; la régression logistique et l'arbre restent en dessous.

Description textuelle du graphique

Un graphique d’intervalles. Pour chaque expérience, un point montre la précision moyenne en validation croisée. Une barre montre une erreur-type de chaque côté. E01, le même score pour tous : 0,114. E02, la règle des séjours antérieurs : 0,182. E03, la régression logistique : 0,210. E04, l’arbre de décision : 0,202. E05, la forêt aléatoire : 0,221. E06, le gradient boosting : 0,219. Une ligne verticale marque le niveau requis par la règle d’une erreur-type (« bar », niveau requis) à 0,216. Seuls la forêt et le gradient boosting l’atteignent.

Ce que cela nous apprend. Lisez ce résultat avec la règle d’une erreur-type. Conservez le modèle le plus simple dont la moyenne est à moins d’une erreur-type du meilleur.

La forêt a la meilleure moyenne (0,2206). Son erreur-type est de 0,0044, donc le niveau requis est de 0,2162. Seuls la forêt et le gradient boosting l’atteignent. Ils sont à égalité.
Entre les deux méthodes d’ensemble, conservez la plus simple : le gradient boosting. Il compte 100 arbres de profondeur 2. Son fichier est environ 190 fois plus petit, comme nous le mesurerons plus bas. C’est le modèle concurrent.
La régression logistique reste sous le niveau requis. Sa moyenne est de 0,2100. L’écart avec le gradient boosting est de 0,0088, soit environ deux erreurs-types. Cet écart ne vient pas de la chance sur un seul pli. Le gradient boosting fait mieux sur les cinq plis, de 0,004 à 0,014. Ici, contrairement à Larkfield dans Quand un modèle simple est meilleur, le modèle complexe gagne au-delà du bruit.

La règle ne nous permet donc pas de parler d’égalité. Le gradient boosting classe mieux. Une question demeure, tirée de Les méthodes d’ensemble : le coût de la complexité. Pour la décision qui compte, le gain est-il assez grand pour justifier la complexité?

Deux expériences que vous devez rejeter

Avant les données de validation, deux autres lignes entrent dans le journal. Toutes deux semblent montrer une amélioration, et toutes deux sont erronées.

E07 : plis par ligne. À quel point une répartition par ligne changerait-elle le score? Évaluez le même pipeline de gradient boosting avec une validation croisée ordinaire à 5 plis, par ligne :

from sklearn.model_selection import StratifiedKFold, cross_val_score

row_folds = StratifiedKFold(5, shuffle=True, random_state=0)
row_scores = cross_val_score(boosting, train[FEATURES], train[TARGET], cv=row_folds, scoring="average_precision")
print("boosting, folds by row:    ", row_scores.round(4), round(row_scores.mean(), 4))
print("boosting, folds by patient:", FOLD_SCORES["E06"].round(4), round(FOLD_SCORES["E06"].mean(), 4))
log_experiment("E07", "E06 scored with folds by row (patients in both parts)", row_scores, decision="reject",
               reason="Not a valid protocol: it scores patients that the model has seen.")

La sortie :

boosting, folds by row:     [0.2087 0.2179 0.2315 0.2171 0.2211] 0.2193
boosting, folds by patient: [0.2191 0.2111 0.2186 0.2109 0.2343] 0.2188

Ce que cela nous apprend. Sur ces données, la répartition par ligne change la moyenne de seulement 0,0005, bien moins que le bruit. C’est un résultat honnête, qui s’explique. Nos entrées décrivent le séjour : sa durée, ses examens, les visites qui le précèdent. Elles changent d’un séjour à l’autre, donc le modèle ne peut pas reconnaître un patient. Supposons qu’une entrée reste fixe pour une personne, ou que les données contiennent un historique détaillé de chaque patient. La même erreur produit alors un grand gain trompeur. La répartition par patient est correcte grâce à ce qu’elle mesure, pas à la taille du changement observé ici. Un résultat obtenu par répartition par ligne répond à une question que Nadia n’a pas posée.

E08 ajoute une caractéristique qui semble très utile. C’est le cas d’échec à la fin de cette page. Sa ligne dans le journal indique « reject : a leak » (rejeter : fuite de données).

Métrique et point de fonctionnement

De la précision moyenne aux appels

L’AP a servi à choisir les familles, car elle résume tout le classement. La décision de Nadia utilise seulement le haut du classement : les 10 % que l’équipe appelle. La métrique de décision est donc la précision@10 % : parmi les patients appelés, la proportion qui reviendrait. Vous l’avez vue sous le nom de précision@k dans Le compromis entre précision et rappel. À capacité fixe, le nombre d’appels reste fixe. La précision et le rappel augmentent et diminuent donc ensemble. La précision est plus facile à utiliser pour Nadia : « 1 appel sur 4 rejoint un patient qui reviendrait ».

Maintenant, et seulement maintenant, entraînez chaque pipeline réglé sur toutes les données d’entraînement, puis calculez les scores sur les données de validation :

from sklearn.metrics import roc_auc_score


def precision_at(y, score, share=0.10):
    """The precision among the `share` of patients with the highest scores: the patients the team calls."""
    y, score = np.asarray(y), np.asarray(score)
    top = np.argsort(-score, kind="stable")[: int(share * len(y))]
    return y[top].mean()


models = {"E03 logistic": logistic, "E04 tree": tree, "E05 forest": forest, "E06 boosting": boosting}
valid_scores = {"E02 rule": stays_rule(valid).to_numpy()}
for name, pipeline in models.items():
    valid_scores[name] = pipeline.fit(train[FEATURES], train[TARGET]).predict_proba(valid[FEATURES])[:, 1]
table = pd.DataFrame({name: {"AP": average_precision_score(valid[TARGET], score),
                             "ROC AUC": roc_auc_score(valid[TARGET], score),
                             "precision@10%": precision_at(valid[TARGET], score)}
                      for name, score in valid_scores.items()}).T
table.round(3)

Le résultat :

                 AP  ROC AUC  precision@10%
E02 rule      0.175    0.614          0.240
E03 logistic  0.206    0.671          0.261
E04 tree      0.195    0.655          0.251
E05 forest    0.219    0.678          0.258
E06 boosting  0.214    0.673          0.268

Précision@10 % sur les données de validation : la proportion des patients appelés qui reviennent dans les 30 jours. Les appels au hasard atteignent 0,113, la règle 0,240, et les quatre modèles de 0,251 à 0,268.

Description textuelle du graphique

Un diagramme à colonnes de la précision@10 % sur les 20 179 séjours de validation. Appels au hasard : 0,113. Règle des séjours antérieurs : 0,240. Régression logistique : 0,261. Arbre de décision : 0,251. Forêt aléatoire : 0,258. Amplification du gradient : 0,268.

Ce que cela nous apprend. L’ordre des familles est le même qu’en validation croisée. Les méthodes d’ensemble classent le mieux (AP de 0,219 et 0,214). Viennent ensuite la régression logistique (0,206), l’arbre (0,195) et la règle (0,175). L’aire sous la courbe ROC donne le même résultat (de 0,614 à 0,678). Mais pour les 10 % en tête, les différences diminuent. Les quatre modèles atteignent de 0,251 à 0,268, et même la règle atteint 0,240. La forêt a la meilleure AP, mais seulement la troisième meilleure précision@10 %. Une différence d’AP ne se traduit pas toujours dans la décision.

Convertissez la précision en nombre de patients. Les données de validation représentent une période de 20 179 sorties. L’équipe fait donc 2 017 appels :

calls = int(0.10 * len(valid))
for name in ["E02 rule", "E03 logistic", "E06 boosting"]:
    reached = round(precision_at(valid[TARGET], valid_scores[name]) * calls)
    print(f"{name:13s} {calls} calls reach {reached} of the {valid[TARGET].sum()} patients who return")

Le résultat :

E02 rule      2017 calls reach 484 of the 2273 patients who return
E03 logistic  2017 calls reach 526 of the 2273 patients who return
E06 boosting  2017 calls reach 540 of the 2273 patients who return

Ce que cela nous apprend. Pour les mêmes 2 017 appels, la régression logistique rejoint 42 patients de plus qui reviendraient que la règle (526 contre 484). L’amplification en rejoint 14 de plus que la régression logistique (540 contre 526). Ces différences sont-elles réelles, ou dues au hasard de ces données de validation? C’est la prochaine question.

Le gain est-il réel? Un intervalle par bootstrap

Les données de validation forment un échantillon de patients. Un autre échantillon donnerait des nombres légèrement différents. Le bootstrap (rééchantillonnage) estime l’ampleur de cette variation. Tirez des données de validation un nouvel échantillon de même taille, avec remise, puis mesurez à nouveau l’écart. Répétez cette opération 1 000 fois. Les 95 % des écarts au centre forment un intervalle à 95 %. C’est la méthode de Comparer des options : un intervalle pour la différence dans le cours de mathématiques. Un détail compte ici : nous tirons des patients, pas des lignes. Les séjours d’un patient ne sont pas indépendants. Ils restent donc ensemble, pour la même raison que nous séparons les données par patient.

def bootstrap_gap(metric, a, b, repeats=1000):
    """metric(b) - metric(a) on 1,000 resampled validation sets. It resamples patients, not rows,
    because the encounters of one patient are not independent."""
    rows_of = list(valid.groupby("patient_nbr").indices.values())
    y = valid[TARGET].to_numpy()
    rng = np.random.default_rng(0)
    gaps = []
    for _ in range(repeats):
        rows = np.concatenate([rows_of[i] for i in rng.integers(0, len(rows_of), len(rows_of))])
        gaps.append(metric(y[rows], b[rows]) - metric(y[rows], a[rows]))
    return np.array(gaps)


def show_gap(name, metric, a, b):
    gaps = bootstrap_gap(metric, a, b)
    low, high = np.percentile(gaps, [2.5, 97.5])
    observed = metric(valid[TARGET], b) - metric(valid[TARGET], a)
    print(f"{name}: gap {observed:+.4f}, 95% interval {low:+.4f} to {high:+.4f}")

Exécutez ce calcul pour trois écarts. Il prend environ 10 secondes :

logistic_score, boosting_score = valid_scores["E03 logistic"], valid_scores["E06 boosting"]
show_gap("boosting - logistic, AP           ", average_precision_score, logistic_score, boosting_score)
show_gap("boosting - logistic, precision@10%", precision_at, logistic_score, boosting_score)
show_gap("logistic - rule, precision@10%    ", precision_at, valid_scores["E02 rule"], logistic_score)

Le résultat :

boosting - logistic, AP           : gap +0.0085, 95% interval +0.0039 to +0.0132
boosting - logistic, precision@10%: gap +0.0069, 95% interval -0.0055 to +0.0168
logistic - rule, precision@10%    : gap +0.0208, 95% interval +0.0064 to +0.0360

Intervalles par bootstrap (1 000 rééchantillonnages de patients) pour les écarts sur les données de validation. L'amplification classe mieux que la régression logistique dans l'ensemble, mais à une capacité de 10 %, l'intervalle inclut zéro. La régression logistique dépasse la règle à cette capacité.

Description textuelle du graphique

Deux panneaux. Le premier, classement selon la précision moyenne : amplification moins régression logistique donne +0,0085, avec un intervalle à 95 % de +0,0039 à +0,0132, entièrement au-dessus de zéro. Le second, précision des 10 % en tête : amplification moins régression logistique donne +0,0069, avec un intervalle de −0,0055 à +0,0168, qui inclut zéro. Régression logistique moins règle donne +0,0208, avec un intervalle de +0,0064 à +0,0360, entièrement au-dessus de zéro.

Ce que cela nous apprend. Trois réponses, chacune avec ses preuves :

L’amplification classe mieux dans l’ensemble. Son écart d’AP de +0,0085 a un intervalle de +0,004 à +0,013, entièrement au-dessus de zéro. Cela confirme les résultats des plis.
À la capacité de Nadia, nous ne pouvons pas les départager. L’écart de précision@10 % est de +0,007 (14 patients sur 2 017 appels). Mais l’intervalle va de −0,006 à +0,017. Sur un autre échantillon de patients, la régression logistique pourrait tout aussi bien être en tête.
La régression logistique dépasse clairement la règle. Son gain de +0,021 (42 patients) a un intervalle de +0,006 à +0,036. Utiliser un modèle en vaut la peine.

Le gain justifie-t-il la complexité?

Utilisez la liste de vérification de Quand un modèle simple est meilleur, question par question.

Une personne doit-elle recevoir des raisons, ou un réviseur doit-il vérifier le modèle? Nous ne connaissons aucune loi qui exige des raisons pour un appel de suivi. L’appel est une offre de soins, pas un refus (ce n’est pas un avis juridique). Mais les infirmières ont besoin d’une raison pour commencer chaque appel (« vous avez eu quatre séjours à l’hôpital cette année »). Et la gouvernance clinique de l’hôpital doit approuver tout nouvel outil. Pour cette raison, les hôpitaux utilisent souvent des scores simples à points pour le risque de réadmission. Deux exemples connus sont l’indice LACE (van Walraven et al., 2010) et le score HOSPITAL (Donzé et al., 2013).
Les entrées sont-elles des images, du texte ou du son? Non, un tableau de 19 colonnes.
Le modèle complexe gagne-t-il par plus d’une erreur-type? Pour le classement, oui : environ deux erreurs-types en validation croisée, et un intervalle par bootstrap au-dessus de zéro en validation.
Le gain vaut-il le coût? Pour la décision qui compte, le gain de précision est de +0,007, et son intervalle inclut zéro. Regardez maintenant le coût :

import io
import time

import joblib


def size_and_speed(pipeline):
    """The size of the saved pipeline (KB) and the time to score the validation set (seconds)."""
    buffer = io.BytesIO()
    joblib.dump(pipeline, buffer)
    start = time.perf_counter()
    pipeline.predict_proba(valid[FEATURES])
    return round(buffer.tell() / 1024), round(time.perf_counter() - start, 3)


for name in ["E03 logistic", "E05 forest", "E06 boosting"]:
    print(name, "size KB, seconds to score valid:", size_and_speed(models[name]))

Le résultat (les temps dépendent de votre ordinateur) :

E03 logistic size KB, seconds to score valid: (12, 0.02)
E05 forest size KB, seconds to score valid: (14936, 0.048)
E06 boosting size KB, seconds to score valid: (79, 0.035)

Ce que cela nous apprend. Les trois sont rapides : moins d’un dixième de seconde pour 20 000 patients. La taille et la vitesse ne sont pas le vrai coût ici. La forêt occupe 15 Mo, contre 79 Ko pour l’amplification et 12 Ko pour la régression logistique. Cela confirme que l’amplification était le bon modèle concurrent. Son vrai coût est ailleurs. Son score vient de 100 arbres. Elle ne peut donc pas donner une liste exacte de raisons pour chaque appel. Après un changement à l’hôpital, il est plus difficile de voir ce qui a changé dans 100 arbres que dans une somme pondérée. Et la gouvernance clinique doit l’examiner.

La décision. Felix voulait des preuves. Les voici en trois phrases. L’amplification classe mieux les patients que la régression logistique, avec un écart supérieur au bruit. À une capacité de 10 %, le gain est de 14 patients sur 2 017 appels, et l’intervalle inclut zéro. Pour ce gain incertain, l’équipe renoncerait aux raisons exactes de chaque appel. Nous choisissons donc la régression logistique (E03) pour un projet pilote. Nous gardons l’amplification (E06) comme modèle concurrent documenté. Ce modèle aura une deuxième chance quand le projet pilote aura de nouvelles données (voir Recommandations).

La règle d’une erreur-type et la liste de vérification de Quand un modèle simple est meilleur fonctionnent ici ensemble. La règle nous a dit que les méthodes d’ensemble classent vraiment mieux. Nous ne prétendons donc pas qu’il y a égalité. La liste nous a dit que « mieux classer » ne signifie pas encore « mieux pour Nadia ». Mettez le journal à jour avec les décisions et l’AP de validation des deux finalistes :

log_experiment("E03", "Logistic regression; C tuned over 0.01, 0.1, 1", FOLD_SCORES["E03"],
               valid_ap=table.loc["E03 logistic", "AP"], decision="keep",
               reason="Chosen: at 10% capacity it ties boosting, and each call gets reasons.")
log_experiment("E06", "Gradient boosting, 100 trees; depth tuned over 2, 3, no limit", FOLD_SCORES["E06"],
               valid_ap=table.loc["E06 boosting", "AP"], decision="reject",
               reason="Better AP beyond one SE, but no clear gain at 10%; kept as the challenger.")
log_experiment("E04", "Decision tree, min leaf 100; depth tuned over 4, 6, 8, 10", FOLD_SCORES["E04"],
               decision="reject", reason="Below logistic regression on every fold.")
log_experiment("E05", "Random forest, 100 trees; min leaf tuned over 50, 20, 10", FOLD_SCORES["E05"],
               decision="reject", reason="Ties boosting on the folds; bigger and slower than boosting.")
show_log()

Le résultat :

 id                                                        change  cv_ap  cv_std valid_ap decision
E01                    Baseline: the same score for every patient  0.114   0.000          baseline
E02                                Baseline: the prior-stays rule  0.182   0.009          baseline
E03                Logistic regression; C tuned over 0.01, 0.1, 1  0.210   0.010    0.206     keep
E04     Decision tree, min leaf 100; depth tuned over 4, 6, 8, 10  0.202   0.007            reject
E05      Random forest, 100 trees; min leaf tuned over 50, 20, 10  0.221   0.010            reject
E06 Gradient boosting, 100 trees; depth tuned over 2, 3, no limit  0.219   0.009    0.214   reject
E07         E06 scored with folds by row (patients in both parts)  0.219   0.007            reject
E08     E03 plus encounters_in_file (counts later encounters too)  0.244   0.008            reject

Ce que cela nous apprend. Le journal raconte toute l’histoire en huit lignes, y compris les deux pièges. Chaque lecteur peut voir ce qui a été essayé, les scores obtenus et les raisons de conserver ou de rejeter chaque essai. Dans le fichier, la colonne des raisons indique que E06 est le modèle concurrent.

Le seuil pour 10 %

L’équipe a besoin d’un seuil : appeler chaque patient dont le score atteint ou dépasse ce seuil. Fixez-le sur les données de validation, comme dans Seuils et calibration : le score du 2 017e patient.

score = valid_scores["E03 logistic"]
threshold = np.sort(score)[::-1][calls - 1]
valid["score"] = score
valid["called"] = valid["score"] >= threshold
print("threshold:", round(threshold, 3), "| called:", valid["called"].sum(), "of", len(valid))

Le résultat :

threshold: 0.18 | called: 2017 of 20179

Ce que cela nous apprend. Un patient ayant un score de 0,18 ou plus reçoit un appel. Exactement 2 017 patients des données de validation atteignent ce seuil. Sur de nouvelles données, la proportion variera un peu. Le test final montre de combien.

Une raison pour chaque appel

Avec la régression logistique, chaque score se décompose en contributions, comme dans Expliquer les prédictions. Prenez un patient appelé typique, celui du milieu de la liste d’appels :

def contributions(pipeline, patient):
    """Coefficient times prepared value, for each column of one patient (as in Explain predictions)."""
    prepare, model = pipeline.named_steps["prepare"], pipeline.named_steps["model"]
    values = np.asarray(prepare.transform(patient[FEATURES]))[0]
    names = [name.split("__", 1)[1][:40] for name in prepare.get_feature_names_out()]
    table = pd.DataFrame({"value": values, "coefficient": model.coef_[0]}, index=names)
    table["contribution"] = table["value"] * table["coefficient"]
    return table


called = valid[valid["called"]].sort_values("score")
patient = called.iloc[[len(called) // 2]]            # a typical called patient: the middle of the list
print("score:", round(patient["score"].iloc[0], 3), "| readmitted:", patient[TARGET].iloc[0])
print(patient[["age", "number_inpatient", "number_emergency", "num_medications"]].to_string(index=False))
reasons = contributions(logistic, patient)
print(reasons[reasons["value"] != 0].sort_values("contribution", ascending=False).head(5).round(2))

Le résultat :

score: 0.226 | readmitted: 1
    age  number_inpatient  number_emergency  num_medications
[80-90)                 4                 0               20
                                     value  coefficient  contribution
number_inpatient                      2.69         0.33          0.89
diabetesMed_Yes                       1.00         0.07          0.07
age_years                             1.21         0.06          0.07
admission_source_Physician Referral   1.00         0.07          0.07
number_diagnoses                      0.82         0.07          0.06

Ce que cela nous apprend. Ce patient, dans la quatre-vingtaine, a eu quatre séjours à l’hôpital dans la dernière année. Ce seul fait ajoute +0,89 à la somme du modèle, plus de dix fois toute autre contribution. L’infirmière peut commencer l’appel avec ce fait. Rappelez-vous l’avertissement d’Expliquer les prédictions : une contribution décrit le modèle, par rapport à un patient moyen. Ce n’est pas la cause d’une réadmission. Ce patient est effectivement revenu dans les 30 jours.

À vous de jouer : une capacité de 20 %

Changez une condition. Le budget de Nadia pourrait augmenter et permettre à l’équipe d’appeler 20 % des patients. Dans le carnet (notebook) Jupyter, calculez precision_at(..., share=0.20) sur les données de validation pour la règle, la régression logistique et l’amplification. Avant de l’exécuter, faites une prédiction : l’écart entre les modèles augmente-t-il ou diminue-t-il?

Afficher la réponse

La règle atteint 0,201, la régression logistique 0,212 et l’amplification 0,219. Les trois diminuent, car les 10 % suivants de la liste comptent moins de patients à risque. L’écart entre l’amplification et la régression logistique est de 0,007, le même qu’à 10 %. Avec deux fois plus d’appels, cela représente environ 28 patients de plus rejoints sur 4 035 appels. Avant de décider, refaites le bootstrap à 20 %. Un gain de même taille reste faible par rapport au bruit. L’idée : quand une décision dépend de la capacité, il faut la vérifier à nouveau si la capacité change.

Erreurs et sous-groupes

Un score pour toutes les données de validation peut cacher des groupes pour lesquels le modèle fonctionne mal. L’analyse par sous-groupes a montré comment les examiner. Chaque tableau par sous-groupe comporte ici six colonnes. Ce sont le nombre de séjours, le taux de réadmission, la proportion appelée, la précision et le rappel des appels, et l’aire sous la courbe ROC du groupe.

Selon l’âge

def slice_table(df, column):
    """Per group: encounters, readmission rate, share called, precision and recall of the calls, AUC."""
    groups = df.groupby(column, observed=True)
    return pd.DataFrame({
        "encounters": groups.size(),
        "readmitted": groups[TARGET].mean(),
        "called": groups["called"].mean(),
        "precision": groups.apply(lambda g: g.loc[g["called"], TARGET].mean()),
        "recall": groups.apply(lambda g: g.loc[g[TARGET] == 1, "called"].mean()),
        "AUC": groups.apply(lambda g: roc_auc_score(g[TARGET], g["score"]) if g[TARGET].nunique() == 2 else np.nan),
    }).round(3)


valid["age_band"] = valid["age"].replace({"[0-10)": "[0-40)", "[10-20)": "[0-40)", "[20-30)": "[0-40)",
                                          "[30-40)": "[0-40)"})
slice_table(valid, "age_band")

Le résultat :

          encounters  readmitted  called  precision  recall    AUC
age_band
[0-40)          1201       0.097   0.081      0.381   0.319  0.762
[40-50)         1985       0.102   0.078      0.277   0.213  0.713
[50-60)         3427       0.092   0.074      0.249   0.199  0.680
[60-70)         4354       0.110   0.092      0.252   0.211  0.665
[70-80)         5339       0.129   0.113      0.255   0.222  0.647
[80-90)         3343       0.125   0.126      0.276   0.278  0.643
[90-100)         530       0.106   0.172      0.143   0.232  0.559

Qualité du classement des patients dans chaque tranche d'âge des données de validation. Le classement se dégrade avec l'âge. Pour les patients de 90 ans ou plus, il est proche du hasard (0,5).

Description textuelle du graphique

Un diagramme à colonnes de l’aire sous la courbe ROC de la régression logistique dans chaque tranche d’âge des données de validation : 0 à 40 ans, 0,76 (1 201 séjours); 40 à 50 ans, 0,71 (1 985); 50 à 60 ans, 0,68 (3 427); 60 à 70 ans, 0,67 (4 354); 70 à 80 ans, 0,65 (5 339); 80 à 90 ans, 0,64 (3 343); 90 à 100 ans, 0,56 (530).

Ce que cela nous apprend. Le modèle classe mieux les patients jeunes (aire sous la courbe de 0,76 avant 40 ans). Il se dégrade à chaque décennie, jusqu’à 0,56 pour les patients de 90 ans ou plus, près du hasard. Deux choses se produisent en même temps pour le groupe le plus âgé. Le modèle les appelle le plus souvent (17,2 %), car l’âge et les séjours antérieurs augmentent le score. Mais ses appels sont les moins précis dans ce groupe. Seulement 14,3 % des appels rejoignent un patient qui revient, près du taux du groupe lui-même, soit 10,6 %. Ces 530 séjours forment un petit sous-groupe. Les nombres varient donc sous l’effet du hasard. Mais la tendance dans les sept tranches est claire. Pour les patients les plus âgés, les infirmières devraient peu se fier au score. Leur propre jugement compte davantage.

Selon les colonnes sensibles

La race, le genre et le code du payeur ne sont pas des entrées. Changent-ils quand même qui reçoit un appel, ou la qualité des appels?

valid["race_group"] = valid["race"].fillna("missing")
valid["payer_group"] = valid["payer_code"].fillna("missing").where(
    valid["payer_code"].isin(["MC", "HM", "SP", "BC", "MD"]) | valid["payer_code"].isna(), "other")
for column in ["race_group", "gender", "payer_group"]:
    print(slice_table(valid, column).to_string(), end="\n\n")

Le résultat :

                 encounters  readmitted  called  precision  recall    AUC
race_group
AfricanAmerican        3849       0.113   0.101      0.237   0.211  0.665
Asian                   121       0.083   0.083      0.100   0.100  0.724
Caucasian             15058       0.115   0.103      0.269   0.241  0.673
Hispanic                419       0.088   0.062      0.385   0.270  0.674
Other                   293       0.092   0.051      0.200   0.111  0.551
missing                 439       0.084   0.068      0.133   0.108  0.659

                 encounters  readmitted  called  precision  recall    AUC
gender
Female                10914       0.113   0.106      0.247   0.233  0.669
Male                   9264       0.112   0.092      0.280   0.229  0.673
Unknown/Invalid           1       0.000   0.000        NaN     NaN    NaN

             encounters  readmitted  called  precision  recall    AUC
payer_group
BC                  900       0.080   0.050      0.200   0.125  0.685
HM                 1323       0.121   0.106      0.264   0.231  0.650
MC                 6365       0.120   0.134      0.249   0.278  0.649
MD                  659       0.096   0.117      0.273   0.333  0.675
SP                  989       0.100   0.080      0.253   0.202  0.706
missing            8038       0.116   0.089      0.282   0.217  0.688
other              1905       0.096   0.056      0.226   0.131  0.636

Ce que cela nous apprend. Lisez chaque tableau en tenant compte des effectifs, comme l’enseignait L’analyse par sous-groupes.

Race. Les deux grands groupes reviennent au même taux (11,3 % et 11,5 %). Le modèle les appelle aussi au même taux (10,1 % et 10,3 %). Mais les appels sont moins précis pour les patients afro-américains (0,237 contre 0,269). Ils rejoignent aussi moins de ceux qui reviennent (rappel de 0,211 contre 0,241). Avec 3 849 séjours, cet écart mérite une surveillance. Les entrées du modèle semblent décrire un peu mieux le risque des patients blancs. Les petits groupes (121 à 439 séjours) sont trop petits pour porter un jugement.
Genre. Les femmes et les hommes reviennent au même taux (11,3 % et 11,2 %). Le rappel est presque le même (0,233 et 0,229). Le modèle appelle les femmes un peu plus souvent (10,6 % contre 9,2 %). Leurs appels sont donc un peu moins précis. Un séjour a pour genre Unknown/Invalid (inconnu ou invalide).
Payeur. Le fichier ne précise pas le sens des codes. MC désigne probablement Medicare, l’assurance américaine pour les personnes de 65 ans ou plus. MD désigne probablement Medicaid, pour les personnes à faible revenu. Le modèle appelle le plus le groupe MC (13,4 %), car il est plus âgé. Il appelle le moins BC et « autre » (5,0 % et 5,6 %), avec un faible rappel. Cela suit en partie leurs taux de réadmission plus faibles (8,0 % et 9,6 %).

Aucun de ces écarts ne vient d’une entrée sensible, car il n’y en a aucune. Ils passent par d’autres entrées, comme l’âge et les séjours antérieurs, et par ce que les données ne consignent pas. La bonne réponse n’est pas de cacher le tableau. Présentez-le, surveillez l’écart selon la race pendant le projet pilote et laissez les infirmières ajouter les patients absents de la liste.

Quels patients le modèle manque-t-il?

Les appels rejoignent 23 % des patients qui reviennent. Qui sont les autres?

returned = valid[valid[TARGET] == 1]
first = returned["number_inpatient"] == 0
print(f"patients who return: {len(returned)}; with no stay in the past year: {first.sum()} ({first.mean():.0%})")
print(f"recall among them: {returned.loc[first, 'called'].mean():.3f};"
      f" among the others: {returned.loc[~first, 'called'].mean():.3f}")

Le résultat :

patients who return: 2273; with no stay in the past year: 1116 (49%)
recall among them: 0.058; among the others: 0.398

Ce que cela nous apprend. C’est l’erreur la plus importante du modèle. La moitié des patients qui reviennent (1 116 sur 2 273) n’ont eu aucun séjour à l’hôpital dans la dernière année. Le modèle appelle seulement 5,8 % d’entre eux, contre 39,8 % des patients ayant des séjours antérieurs. Le modèle, comme la règle, trouve surtout les patients qui reviennent souvent. Pour un patient à son premier séjour, le fichier contient peu d’éléments qui prédisent un retour. C’est une limite des données, pas de la famille de modèles : les quatre familles dépendent des mêmes colonnes.

Le test final, une seule fois

Tout est décidé : les caractéristiques, le modèle (régression logistique, C = 0.01, entraîné sur les données d’entraînement) et le seuil (0,18). Ouvrez maintenant les données de test réservées, une seule fois, comme dans Seuils et calibration : tester une seule fois. Quel que soit le résultat, nous ne revenons pas modifier le modèle.

test_score = logistic.predict_proba(test[FEATURES])[:, 1]
called = test_score >= threshold
returned = test[TARGET].to_numpy() == 1
print(f"test: {len(test)} encounters, {returned.sum()} returned within 30 days")
print(f"called with the valid threshold: {called.sum()} ({called.mean():.1%})")
print(f"precision {returned[called].mean():.3f}, recall {called[returned].mean():.3f}")
print(f"AP {average_precision_score(returned, test_score):.3f}, ROC AUC {roc_auc_score(returned, test_score):.3f},"
      f" precision@10% {precision_at(returned, test_score):.3f}")
print(f"the rule on test: precision@10% {precision_at(returned, stays_rule(test)):.3f}")

Le résultat :

test: 19608 encounters, 2278 returned within 30 days
called with the valid threshold: 1870 (9.5%)
precision 0.270, recall 0.221
AP 0.217, ROC AUC 0.663, precision@10% 0.264
the rule on test: precision@10% 0.247

Ce que cela nous apprend. Sur 19 608 séjours de patients qu’aucune étape n’a vus, le modèle se comporte comme en validation :

Le seuil de validation fait appeler 9,5 % des patients de test, près des 10 % prévus. En pratique, la proportion varie d’une semaine à l’autre. L’équipe devrait donc compter les appels et s’ajuster.
Environ 27 appels sur 100 rejoignent un patient qui revient (précision de 0,270). Au hasard, ce serait environ 11 sur 100; avec la règle, 25. Les appels rejoignent 22 % de tous les patients qui reviennent.
L’AP de 0,217 et l’aire sous la courbe ROC de 0,663 sont proches des valeurs de validation (0,206 et 0,671). Le modèle n’a pas fait de surapprentissage des choix faits sur les données de validation.

Voici le nombre à présenter : une précision de 0,27 à une capacité de 10 %, contre 0,25 pour la règle simple et 0,11 au hasard.

Recommandations

Pour Nadia et Felix, dans leurs mots :

Lancez un projet pilote avec le modèle de régression logistique. Appelez les patients dont le score atteint 0,18, soit environ 1 sur 10. Prévoyez qu’environ 1 appel sur 4 joindra un patient qui reviendrait, contre 1 sur 9 au hasard. Par rapport à la règle actuelle, cela représente environ 40 patients de plus joints par 2 000 appels.
Expliquez les raisons de chaque appel. Montrez au personnel infirmier les deux ou trois contributions principales avec les valeurs du patient, par exemple « 4 hospitalisations cette année ». Ce sont les raisons du modèle, pas un diagnostic.
Laissez le personnel infirmier ajouter des patients. Le modèle manque la plupart des patients qui reviennent sans hospitalisation antérieure. Un membre du personnel infirmier qui connaît une raison (aucune aide à domicile, nouveau traitement à l’insuline) doit pouvoir ajouter un patient.
Mesurez correctement les résultats du projet pilote. Notez le résultat pour chaque patient appelé. Une comparaison fournit les meilleures preuves. Pendant quelques semaines, choisissez certains appels avec le modèle et d’autres avec la règle, puis comparez les taux de réadmission. Décidez ensuite si vous gardez le modèle.
Donnez une deuxième chance au modèle concurrent. Le boosting (amélioration successive de modèles) a mieux classé les patients dans chaque comparaison du classement complet. Mais le gain à 10 % restait dans la marge du bruit. Quand le projet pilote aura plusieurs milliers de nouveaux résultats, refaites la comparaison sur les nouvelles données. Changez seulement si le gain à la capacité disponible dépasse clairement zéro et justifie la perte de raisons exactes.
Surveillez les sous-groupes. Chaque mois, présentez la précision et le rappel par âge, race, genre et payeur. Si l’écart pour les patients afro-américains augmente, ou si les patients les plus âgés sont souvent appelés avec une faible précision, réexaminez le modèle.
Réentraînez le modèle sur vos propres données récentes. Ce modèle a appris à partir d’hôpitaux américains de 1999 à 2008. Avant toute utilisation réelle, un hôpital doit l’entraîner et le tester sur ses propres patients.

À quoi ne pas utiliser le modèle. Ne l’utilisez pas pour décider qui peut quitter l’hôpital, qui reçoit moins de soins ou combien un patient paie. Il ne pose aucun diagnostic. Ne l’utilisez pas pour les patients sans diabète, pour les enfants, ni dans un hôpital dont les données n’ont pas servi à le tester. Un faible score ne signifie pas « sans danger » : la moitié des patients qui reviennent n’ont aucune hospitalisation antérieure, et la plupart obtiennent de faibles scores.

Limites et prochaines étapes

Soyez honnêtes sur le plafond. Ces données prédisent faiblement la réadmission. Le meilleur modèle de cette étude de cas atteint une aire sous la courbe ROC d’environ 0,67. Une valeur de 0,5 correspond au hasard, et 1,0 est parfaite. Trois faits montrent que la limite vient des données, pas du choix du modèle :

Quatre familles très différentes, toutes réglées, terminent à moins de 0,02 d’AP les unes des autres (0,202 à 0,221). Si ces colonnes contenaient une forte tendance cachée, les modèles flexibles en auraient trouvé beaucoup plus.
La règle simple, avec une seule colonne, obtient déjà une bonne partie du résultat : 0,240 contre 0,261 à la capacité disponible.
La moitié des réadmissions concernent des patients sans hospitalisation antérieure, et le modèle n’en trouve presque aucun.

La réadmission dépend de facteurs que ce fichier ne consigne pas. Par exemple, qui aide à domicile, si le patient peut obtenir ses médicaments et si une visite de suivi a lieu. Le fichier renseigne aussi peu sur la gravité réelle de la maladie du patient : il n’a aucun résultat de laboratoire, sauf deux tests de diabète. De meilleures données aideraient davantage qu’un meilleur algorithme.

Autres limites :

Aucune date. Nous n’avons pas pu séparer les données dans le temps, donc nous ignorons comment le modèle vieillit. Un vrai projet doit tester le modèle sur les mois les plus récents.
Données anciennes. Hôpitaux américains, 1999–2008, codes ICD-9. Les soins, la codification et les assurances ont changé.
Une seule étiquette. <30 (moins de 30 jours) indique un retour dans un des hôpitaux des données. Un retour dans un autre hôpital peut apparaître comme « aucun retour ».
Hospitalisations corrélées. Nous avons conservé toutes les hospitalisations d’un patient. Un modèle limité aux premières hospitalisations, comme chez Strack et ses collègues, répondrait à une question légèrement différente.

Prochaines étapes, par ordre de valeur attendue :

Ajoutez les deuxième et troisième diagnostics et le nombre de maladies chroniques (un score de comorbidité, comme dans LACE).
Ajoutez les colonnes de médicaments.
Calibrez le score pour que 0,25 signifie « 1 sur 4 », comme dans Seuils et calibration.
Planifiez le projet pilote comme une expérience.

Cas d’échec : un compte qui regarde dans le futur

Erreur courante

Une variable qui utilise des hospitalisations ultérieures

Symptôme : Dev, un ingénieur d’une autre équipe, regarde votre carnet (notebook) Jupyter. « Les patients qui ont beaucoup d’hospitalisations reviennent plus souvent. Pourquoi ne pas compter toutes leurs hospitalisations? » Il ajoute une variable, encounters_in_file : le nombre de lignes du patient dans le fichier complet. L’AP en validation croisée passe de 0,210 à 0,244, un gain supérieur à celui de toute famille de modèles :

in_file = data["patient_nbr"].map(data["patient_nbr"].value_counts())   # all encounters, also later ones
train["encounters_in_file"] = in_file.loc[train.index]
leaky = build_pipeline(LogisticRegression(C=0.01, max_iter=1000), numbers=NUMBERS + ["encounters_in_file"])
leaky_scores = cross_val_score(leaky, train[FEATURES + ["encounters_in_file"]], train[TARGET], cv=folds,
                               scoring="average_precision")
print("with encounters_in_file:", leaky_scores.round(4), round(leaky_scores.mean(), 4))
print("without (E03):          ", FOLD_SCORES["E03"].round(4), round(FOLD_SCORES["E03"].mean(), 4))
log_experiment("E08", "E03 plus encounters_in_file (counts later encounters too)", leaky_scores,
               decision="reject", reason="A leak: at discharge, nobody knows the later encounters.")

with encounters_in_file: [0.2543 0.2307 0.2389 0.2463 0.2495] 0.2439
without (E03):           [0.2064 0.204  0.205  0.2044 0.2302] 0.21

Cause : le compte inclut les hospitalisations qui arrivent après celle pour laquelle nous faisons une prédiction. Un patient qui revient dans les 30 jours a, par définition, au moins une hospitalisation de plus dans le fichier. La variable contient donc une partie de la réponse. Les plis par patient ne protègent pas contre cette fuite, car elle se produit dans le temps, pas entre les patients. Au congé, le personnel infirmier ne peut pas connaître le nombre d’hospitalisations futures d’un patient.

Solution : demandez, pour chaque variable : « Connaîtrions-nous cette valeur le jour du congé? » Utilisez seulement des comptes du passé, comme number_inpatient, que la source a calculé pour l’année avant l’hospitalisation. Un gain beaucoup plus grand que toute modification honnête du journal est un signal d’alerte : cherchez une fuite avant de célébrer. Notez dans le journal que l’expérience est rejetée, avec la raison, pour que personne ne la tente à nouveau.

Conséquence dans le travail en IA : pendant le projet pilote, les hospitalisations futures d’un nouveau patient n’existent pas encore. La variable vaudrait donc 1 pour la plupart des patients. Les scores baisseraient, et les appels iraient aux mauvais patients. L’équipe ferait confiance à un nombre (0,244) qui n’a jamais été réel. Les fuites temporelles expliquent souvent pourquoi un modèle prometteur dans un carnet échoue en utilisation réelle.

Faites-en un projet de portfolio

Cette étude de cas est un projet complet qu’un employeur peut lire en dix minutes. Mettez-le dans un dépôt public.

Quoi mettre dans le dépôt :

readmission-model/
  README.md               the summary below
  MODEL_CARD.md           intended use, not for, data, evaluation, slices, limits
  DATA.md                 the dataset card: source, DOI, licence, attribution, your changes
  notebook.ipynb          the notebook, run from top to bottom, with its outputs
  experiment_log.csv      every experiment, E01 to E08
  requirements.txt        pandas==3.0.6, scikit-learn==1.9.1, matplotlib

Vous n’avez pas besoin d’inclure le fichier de données lui-même dans un commit (enregistrement des modifications). CC BY 4.0 le permet, avec attribution et mention de vos modifications, mais un lien et une somme de contrôle suffisent.

Un modèle de README :

Afficher le modèle de README

# Who should the care team call? A 30-day readmission model

**Problem.** A care-transitions team can call 10% of discharged patients with diabetes.
Which patients are most likely to come back to hospital within 30 days?

**Data.** Diabetes 130-US hospitals 1999–2008 (UCI, doi:10.24432/C5230J, CC BY 4.0;
Strack et al. 2014). 99,343 stays of 69,990 patients after removing deaths and hospice
discharges. Changes: see DATA.md. Race, gender and payer are used for checks only.

**Approach.** Split by patient (60/20/20). Baselines: the same score for all, and the
rule of thumb "most hospital stays in the past year". Four model families in one
pipeline, tuned on the same 5 patient folds (average precision), one-standard-error rule.
Bootstrap intervals (by patient) for the gaps on validation. Test set used once.

**Results (test set, 19,608 stays).** At a 10% capacity, precision 0.27 (random 0.11,
rule of thumb 0.25); ROC AUC 0.66. Gradient boosting ranked better than logistic
regression (AP +0.009), but at the capacity the gain was within the noise, so the
readable model was chosen.

**Limits.** Weak signal (half of the readmissions have no earlier stay and are rarely
found); no dates, so no time split; US data from 1999–2008. Not for decisions about
discharge or care.

La fiche du modèle : rédigez MODEL_CARD.md avec les sections de Documenter les limites d’utilisation :

Utilisation prévue : classer les patients diabétiques ayant reçu leur congé pour un appel de suivi, avec une capacité de 10 % et une révision de la liste par le personnel infirmier.
Utilisations exclues : décisions de congé, réduction des soins, prix, autres patients.
Données : la fiche du jeu de données et vos modifications.
Pourquoi ce modèle : la décision de cette étude de cas, avec les chiffres.
Évaluation : les chiffres du test et les références de base.
Sous-groupes : le tableau par âge et les trois tableaux de variables sensibles.
Limites : le plafond, l’absence de dates, les données anciennes.

Comment le décrire en entrevue, en deux phrases : « J’ai créé un modèle de réadmission à 30 jours sur 100 000 hospitalisations réelles, pour une équipe pouvant appeler 10 % des patients. J’ai comparé quatre familles de modèles avec un protocole par patient et choisi une régression logistique lisible, car le petit gain du boosting à la capacité de l’équipe restait dans la marge du bruit. »

Attendez-vous à la question suivante : « Pourquoi pas le meilleur modèle? » Votre réponse se trouve dans les preuves de Le gain justifie-t-il la complexité?.

Allez plus loin :

Ajoutez diag_2, diag_3 et les 22 colonnes de médicaments. Une famille gagne-t-elle plus d’une erreur-type?
Créez un score de comorbidité à partir des trois diagnostics et testez un score à points de type LACE comme troisième référence de base.
Calibrez la régression logistique et tracez un graphique de fiabilité.
Entraînez le modèle uniquement sur les premières hospitalisations, comme Strack et ses collègues, puis comparez.
Utilisez encounter_id comme ordre temporel approximatif, indiquez clairement cette hypothèse et testez sur les dernières 20 % des hospitalisations.
Rédigez un plan de projet pilote d’une page : durée en semaines, comparaison du modèle et de la règle, résultat qui entraînerait l’arrêt.

Vérification des connaissances

Vérifiez votre compréhension

Le boosting a une AP en validation croisée supérieure à celle de la régression logistique d'environ deux erreurs-types. Il est aussi meilleur sur les cinq plis. Felix dit : alors, la règle d'une erreur-type choisit le boosting, donc nous devons l'utiliser. Qu'est-ce qui est juste et qu'est-ce qui manque dans son argument?

Afficher la réponse

Il a raison : la règle ne permet pas d’égalité ici. La régression logistique (0,2100) reste sous la barre de 0,2162. L’intervalle bootstrap pour l’écart d’AP sur les données de validation (+0,004 à +0,013) est supérieur à zéro. Le boosting classe vraiment mieux les patients dans l’ensemble. Ce qui manque, c’est la décision. La règle compare la métrique globale; l’équipe utilise seulement les premiers 10 %. À ce niveau, le gain de précision est de +0,007 (14 patients sur 2 017 appels). Son intervalle (−0,006 à +0,017) inclut zéro. La liste de vérification de Quand un modèle simple est meilleur demande alors si le gain vaut le coût. Pour un gain incertain, le boosting abandonnerait les raisons exactes de chaque appel et compliquerait la surveillance et la révision. Les preuves appuient donc la régression logistique pour le projet pilote, avec le boosting comme concurrent.

Vérifiez votre compréhension

Un collègue sépare les hospitalisations par ligne et rapporte une AP de test. Ici, les plis par ligne ont changé l'AP du boosting en validation croisée de seulement 0,0005. Une séparation par ligne est-elle acceptable pour ce projet? Expliquez.

Afficher la réponse

Non. Avec une séparation par ligne, 41 % des lignes vérifiées appartiennent à des patients utilisés pour entraîner le modèle. Le score répond alors à « à quel point le modèle connaît-il ces patients? » plutôt qu’à « prédit-il bien de nouveaux patients? ». Sur ces données, le changement est faible parce que les entrées décrivent l’hospitalisation, pas la personne. Mais c’est une propriété de ces colonnes, et personne ne peut s’y fier à l’avance. Une entrée personnelle fixe ou un historique plus riche rendrait la fuite importante. Le protocole doit mesurer la question de l’équipe : séparez donc par patient (et par période, quand il y a des dates).

Vérifiez votre compréhension

Nadia demande : le modèle ne produit qu'environ 27 appels pertinents sur 100. Devrait-on attendre un meilleur algorithme? Utilisez les preuves de cette étude de cas.

Afficher la réponse

Un meilleur algorithme aiderait probablement peu. Quatre familles très différentes terminent à moins de 0,02 d’AP les unes des autres. Les familles flexibles gagnent peu à la capacité disponible. La règle simple avec une colonne atteint 0,247 sur les données de test, près du 0,270 du modèle. La moitié des patients qui reviennent n’avaient aucune hospitalisation antérieure, et le modèle n’en trouve que 5,8 %. Le fichier ne consigne pas ce qui prédit leur retour : soutien à domicile, médicaments, soins de suivi. De meilleures données aideraient davantage. En attendant, le modèle est clairement meilleur que le hasard (0,11) et que la règle (intervalle de +0,006 à +0,036 à la capacité disponible). Un projet pilote vaut donc la peine, en permettant au personnel infirmier d’ajouter des patients.

Sources

Clore, J., Cios, K., DeShazo, J., & Strack, B. (2014). Diabetes 130-US Hospitals for Years 1999-2008 [Jeu de données]. UCI Machine Learning Repository. doi:10.24432/C5230J. CC BY 4.0.
Strack, B., DeShazo, J. P., Gennings, C., Olmo, J. L., Ventura, S., Cios, K. J., & Clore, J. N. (2014). Effet de la mesure de l’HbA1c sur les taux de réadmission à l’hôpital : analyse de 70 000 dossiers de patients d’une base de données cliniques. BioMed Research International, 781670. doi:10.1155/2014/781670.
Obermeyer, Z., Powers, B., Vogeli, C., & Mullainathan, S. (2019). Analyse du biais racial dans un algorithme utilisé pour gérer la santé des populations. Science, 366(6464), 447–453. doi:10.1126/science.aax2342.
van Walraven, C., et al. (2010). Élaboration et validation d’un indice pour prédire le décès précoce ou la réadmission non planifiée après le congé de l’hôpital vers la communauté (l’indice LACE). CMAJ, 182(6), 551–557. doi:10.1503/cmaj.091117.
Donzé, J., Aujesky, D., Williams, D., & Schnipper, J. L. (2013). Réadmissions à l’hôpital potentiellement évitables à 30 jours chez les patients en médecine : élaboration et validation d’un modèle de prédiction (le score HOSPITAL). JAMA Internal Medicine, 173(8), 632–638. doi:10.1001/jamainternmed.2013.3023.

Récapitulatif

Préparez des données réelles avec une raison pour chaque étape. ? (valeur manquante) indique une absence, un test vide signifie « non mesuré », les codes deviennent des mots, et les hospitalisations se terminant par un décès ou des soins palliatifs sont retirées.
Cadrez la prédiction au congé : entrées connues ce jour-là, aucun compte futur, colonnes sensibles réservées aux vérifications par sous-groupes.
Séparez par patient, car 46 % des lignes appartiennent à des patients ayant plusieurs hospitalisations. Ici, une séparation par ligne changeait peu le score, mais elle mesure la mauvaise question.
Un seul protocole pour quatre familles : les mêmes 5 plis par patient, la même métrique (AP), le même pipeline, chaque essai dans le journal. Le boosting et la forêt dépassent la régression logistique d’environ deux erreurs-types.
Fondez le choix sur la décision. À la capacité de 10 %, l’intervalle bootstrap du gain du boosting (+0,007) inclut zéro. Le modèle lisible a donc gagné. Sur les données de test : précision de 0,27, contre 0,25 pour la règle et 0,11 au hasard.
Soyez honnêtes sur le plafond : ces données prédisent faiblement la réadmission, et la moitié des retours concernent des patients que le modèle trouve rarement.

Ensuite, une étude de cas sans cible. Vous trouvez des groupes de clients d’une boutique en ligne selon la récence, la fréquence et le montant de leurs achats. Puis vous comparez les groupes à ce que les clients font ensuite.

Leçon suivante : Segmenter les clients d'une boutique en ligne