Un cadre général d’apprentissage automatique pour l’analyse de survie

L’analyse de survie est une branche des statistiques qui fournit un cadre pour l’analyse des données de temps avant événement, c’est-à-dire que le résultat est défini par le temps qu’il faut jusqu’à ce qu’un événement se produise. L’analyse de ces données nécessite des techniques spécialisées car, contrairement aux tâches de régression ou de classification standard,

• le résultat ne peut souvent pas être pleinement observé (censure, troncature),

• les caractéristiques peuvent changer de valeur au cours de la période d’observation (caractéristiques variant dans le temps (TVF),

• l’association de la ou des caractéristiques avec le résultat évolue dans le temps (effets variant dans le temps (TVE)),

• un ou plusieurs autres événements se produisent qui rendent impossible l’observation de l’événement d’intérêt (risques concurrents (CR)),

• plus généralement, dans un contexte multi-états, les unités d’observation peuvent passer de et vers différents états (modèles multi-états (MSM)).

Le fait de ne pas prendre ces questions en considération entraîne généralement des estimations biaisées, une interprétation incorrecte des effets des caractéristiques sur le résultat, une perte de précision prédictive ou une combinaison de ces éléments.
Dans ce travail, nous utilisons une reformulation de la tâche de survie en une tâche de régression standard qui fournit une approche holistique de l’analyse de survie. Dans ce cadre, la censure, la troncature et les fonctionnalités variables dans le temps (TVF) peuvent être incorporées par des transformations de données spécifiques et des extensions d’effets variables dans le temps (TVE) ainsi que des risques concurrents et des modèles multi-états peuvent être ré-exprimés en termes de effets d’interaction. Cette abstraction de la tâche de survie loin des algorithmes spécialisés est illustrée schématiquement dans la figure

1. Le prétraitement approprié à la tâche (sous-graphe le plus à gauche) produit un format de données normalisé qui permet d’estimer les taux de risque conditionnel en utilisant n’importe quel algorithme d’apprentissage qui peut minimiser le log-vraisemblance négatif de Poisson, tel que GBT, réseaux de neurones profonds (DNN) , méthodes basées sur la régularisation et autres (sous-graphe du milieu).

Dans le contexte de l’analyse de survie, une observation consiste généralement en un tuple $(t_{\mathrm{je}},{\delta }_{\mathrm{je}},X_{\mathrm{je}})$, où $t_{\mathrm{je}}$ est l’heure de l’événement observé pour l’unité d’observation $\mathrm{je}=1,\dots ,n$, ${\delta }_{\mathrm{je}}\in \{0,1\}$ est l’indicateur d’événement ou d’état (c’est-à-dire 1 si l’événement s’est produit, 0 si l’observation est censurée) et $X_{\mathrm{je}}$ est le $p$

-vecteur de caractéristiques dimensionnelles. La présence de la censure nécessite des techniques d’estimation spéciales, car le temps de l’événement ne peut pas être observé lorsque la censure se produit avant l’événement d’intérêt. Ainsi

$t_{\mathrm{je}}=min(T_{\mathrm{je}},C_{\mathrm{je}})$, où $T_{\mathrm{je}}$ et $C_{\mathrm{je}}$variables aléatoires de l’heure de l’événement et de l’heure de censure, respectivement. Un exemple classique est le temps jusqu’à la mort lorsque la censure se produit lorsque les patients abandonnent l’étude (sans rapport avec l’événement d’intérêt, $T_{\mathrm{je}}\perp C_{\mathrm{je}}$

). La troncature à gauche se produit lorsque l’événement d’intérêt s’est déjà produit avant que le sujet puisse être inclus dans l’échantillon et présente ainsi une forme de biais d’échantillonnage.
Dans certains contextes, un autre événement pourrait empêcher l’observation de l’événement d’intérêt ou modifier la probabilité de son occurrence. Dans ce cas on parle de risques concurrents (CR), donc l’observation consiste en

$(t_{\mathrm{je}},{\delta }_{\mathrm{je}},k,X_{\mathrm{je}})$, où $k=1,\dots ,K$ indique le type d’événement survenu à, $t_{\mathrm{je}}$ si ${\delta }_{\mathrm{je}}=1$. Plus généralement, il peut y avoir plusieurs états depuis et vers lesquels les unités d’observation peuvent faire la transition. On parle alors de modèles multi-états (MSM) et $k$ est un indicateur de différentes transitions (cf. Eq. (8)).

En général, le but de l’analyse de survie est d’estimer la distribution conditionnelle des temps d’événement définie par la probabilité de survie $S\left(t|X\right)=P(T>t|X)$. Alors que certaines méthodes se concentrent sur l’estimation de $S\left(t|X\right)$ directement, il est souvent plus pratique d’estimer le danger (log-)

à partir duquel $S\left(t|X\right)$ suit comme

Ici, nous représentons (1) passant par

où $g$ est une fonction générale de potentiellement TVF

Dans ce travail, nous approchons (3) en utilisant le modèle exponentiel par morceaux [13]. Laisser $t_{\mathrm{je}}$ l’événement observé ou l’heure de censure et ${\delta }_{\mathrm{je}}\in \{0,1\}$ l’indicateur de censure ou d’événement respectif pour les unités d’observation $\mathrm{je}=1,\dots ,n$. La distribution des temps de censure peut dépendre des caractéristiques mais est supposée être indépendante du processus de temps de l’événement $T$. En répartissant le suivi, c’est-à-dire le laps de temps sous enquête, en $j=1,\dots ,J$ intervalles avec points de coupure ${\kappa }_0=0<\cdots <{\kappa }_J$ et cloisons
$({\kappa }_0,{\kappa }_1],\dots ,({\kappa }_{j–1},{\kappa }_j],\dots ({\kappa }_{J–1},{\kappa }_J]$,
nous pouvons réécrire (3) en utilisant des taux de risque constants par pièce

avec $t_j$ une représentation du temps dans l’intervalle $j$, par exemple., $t_j:={\kappa }_j$ et $X_{\mathrm{je}j}$ la valeur du TVF en intervalle $j$. En fonction de la résolution souhaitée, des points de coupure supplémentaires peuvent être introduits à chaque instant auquel les valeurs de caractéristiques sont mises à jour, sinon plusieurs valeurs de caractéristiques doivent être agrégées dans un intervalle. Ce modèle suppose que seule la valeur actuelle de $X_{\mathrm{je}j}$ affecte le danger dans l’intervalle $j$, mais des approches plus sophistiquées ont été proposées dans ce cadre qui prennent en compte toute l’histoire de TVF [3].
Les dangers constants par morceau impliquent des contributions de vraisemblance log exponentielle par morceau

où $J_{\mathrm{je}}$ est le dernier intervalle dans lequel l’unité d’observation $\mathrm{je}$ a été observée, de telle sorte que $t_{\mathrm{je}}\in ({\kappa }_{J_{\mathrm{je}}–1},{\kappa }_{J_{\mathrm{je}}}]$ et

Exemples concrets du type de transformations de données nécessaires pour obtenir (6) pour les données censurées à droite (y compris TVF) sont fournies dans [2] (cf. Tableaux 1 et 2.

Utilisation de l’hypothèse de travail ${\delta }_{\mathrm{je}j}\stackrel{\mathrm{je}\mathrm{je}\mathrm{ré}}{\sim }Po\mathrm{je}sson({\mu }_{\mathrm{je}j}={\lambda }_{\mathrm{je}j}{t}_{\mathrm{je}j})$ et avec $F\left({\delta }_{\mathrm{je}j}\right)$ la fonction de densité de Poisson, [13] a montré que la log-vraisemblance de Poisson

est proportionnel à (5) et donc le premier peut être minimisé à l’aide de la régression de Poisson. Noter que (7) peut être directement étendu au paramètre avec des temps d’événement tronqués à gauche [16] en remplaçant $j=1$ avec $j_{\mathrm{je}}$, le premier intervalle dans lequel l’unité d’observation $\mathrm{je}$ fait partie de l’ensemble des risques. L’espérance est définie par ${\mu }_{\mathrm{je}j}={\lambda }_{\mathrm{je}j}{t}_{\mathrm{je}j}=\exp (g(Xjej,tj)+Journal(tjej))$. Pour l’estimation, $\mathrm{Journal}(tjej)$ est inclus comme compensation, donc le taux de risque $\frac{{\mu }_{\mathrm{je}j}}{t_{\mathrm{je}j}}$ est définie comme l’espérance conditionnelle d’avoir un événement dans l’intervalle $j$ divisé par le temps à risque. Notez que l’hypothèse de Poisson est simplement un véhicule de calcul pour l’estimation de l’aléa (4) plutôt qu’une hypothèse sur la distribution des heures de l’événement. Malgré la partition du suivi en intervalles, il s’agit d’une méthode pour les temps d’événements continus car les informations sur le temps sous risque dans chaque intervalle sont contenues dans le décalage et donc utilisées lors de l’estimation. Le nombre et l’emplacement des points de coupure contrôlent l’approximation du danger et pourraient donc être considérés comme un paramètre de réglage potentiel. D’après notre expérience, cependant, fixer des seuils aux moments uniques de l’événement $\{t_{\mathrm{je}}:{\delta }_{\mathrm{je}}=1,\mathrm{je}=1,\dots ,n\}$ dans les données de formation conduit toujours à une bonne approximation (au moins avec une régularisation suffisante) car le nombre de points de coupure augmentera dans les zones avec de nombreux événements. Pour des ensembles de données plus volumineux, cependant, nous recommandons de définir ces seuils sur un sous-échantillon représentatif plus petit de l’ensemble de données (cf. 4).

Pour l’extension de (3) aux HSH, nous définissons

comme danger spécifique à la transition pour la transition indexée par $k$, c’est à dire., $k$ est un index des transitions $m_k\to m_k^{\prime }$ où $m_k$ est l’état initial et $m_k^{\prime }$ un état transitoire ou absorbant. L’ensemble des transitions possibles est donné par $\{m_k\to m_k^{\prime }:k=1,\dots ,K\}\subseteq \{m\to m^{\prime }:m,m^{\prime }\in \{0,\dots ,M\},m\ne m^{\prime }\}$, où $M+1$ désigne le nombre total d’états possibles.

Tableau 1 montre comment les données doivent être transformées pour estimer (8) via les PEM pour l’établissement des risques concurrents, c’est-à-dire $k=1,\dots ,K$ est un index des transitions $m_k=0\to m_k^{\prime },m_k^{\prime }\in \{1,\dots ,M\}$; un exemple concret est donné dans le tableau 2. Pour chaque $\mathrm{je}=1,\dots ,n$, il y a une ligne pour chaque intervalle où l’unité d’observation était exposée à un risque pour une transition spécifique. Ainsi, un ensemble de données est créé pour chaque transition de sorte que les transitions vers l’état $m_k^{\prime }$ sont codés comme $1$ et tout le reste, c’est-à-dire la censure et la transition vers d’autres états, est codé comme $0$. Ces ensembles de données spécifiques à la transition, chacun contenant un vecteur d’entités avec l’index de transition $k$, sont ensuite concaténées. Notez que nous avons utilisé les mêmes points de partage d’intervalle ${\kappa }_j$ pour toutes les transitions dans le tableau 1. Cependant, il serait également possible de choisir des seuils spécifiques à la transition ${\kappa }_{j,k}$, ou, plus généralement, utiliser plusieurs échelles de temps [20]. Dans le contexte général multi-états, le nombre d’unités d’observation à risque peut dépendre de la transition et des intervalles visités par $\mathrm{je}$ sont définis par $(t_{\mathrm{je},m},{\kappa }_{j_{\mathrm{je},k}}],\dots ,({\kappa }_{J_{\mathrm{je},k}–1},{\kappa }_{J_{\mathrm{je},K}}]$, où $t_{\mathrm{je},m}$ est le moment auquel $\mathrm{je}$ entre en état $m$.

Dans la section 3 nous évaluons l’approche suggérée en utilisant une implémentation basée sur GBT que nous appelons GBT (PEM). En tant que moteur de calcul concret, nous avons utilisé la bibliothèque de boosting de gradient extrême (XGBoost) [8] sans aucune modification de l’algorithme. Par conséquent, toutes les fonctionnalités de la bibliothèque peuvent être utilisées directement lors de l’estimation du danger sur l’ensemble de données transformé. Notez cependant qu’en fonction de l’algorithme utilisé, il faut être capable de spécifier un décalage lors de l’estimation et éventuellement d’autres paramètres spécifiques à un algorithme ou à une implémentation. Par exemple, lorsque vous utilisez XGBoost pour estimer le GBT (PEM), la fonction objectif doit être définie sur l’objectif de Poisson et le score de base doit être défini sur 1, car la valeur par défaut de 0,5 impliquerait un décalage erroné, tandis que $\mathrm{Journal}(1)=0$. Le décalage ($\mathrm{Journal}(tjej)$) doit être attaché aux données via l’argument de marge de base lors de l’estimation. En revanche, pour la prédiction de l’aléa conditionnel $\lambda \left(t|X_{\mathrm{je}}\right)={\lambda }_{\mathrm{je}j}$ basé sur de nouveaux points de données, le décalage doit être omis, sinon l’algorithme prédira ${\hat{\mu }}_{\mathrm{je}j}={\hat{\lambda }}_{\mathrm{je}j}\cdot t_{\mathrm{je}j}$ (l’attente) au lieu de ${\hat{\lambda }}_j\left(X\right)$ (le danger). Cependant, lors de la prédiction du danger cumulatif ou de la probabilité de survie, le temps sous risque dans chaque intervalle doit être pris en compte, de sorte que

Nous réalisons un ensemble d’expériences de référence avec des ensembles de données du monde réel et synthétiques, en utilisant exclusivement des données disponibles ouvertement et directement, y compris un sous-ensemble d’ensembles de données de publications récentes sur les forêts de survie aléatoire oblique (ORSF, [23]) et DeepHit [28]. Il a été démontré que DeepHit et ORSF surpassent d’autres approches telles que RSF [22], forêts conditionnelles [18], régression de Cox régularisée [12] et DeepSurv [31]. Nous comparons notre approche dans les benchmarks avec les deux algorithmes, qui sont évalués séparément sur la base des mesures d’évaluation utilisées dans les publications respectives pour assurer la comparabilité.
Tout le code pour effectuer des analyses respectives ainsi que des fichiers supplémentaires supplémentaires sont fournis dans un GitHub
dépôt: https://github.com/adibender/machine-learning-for-survival-ecml2020.

Les ensembles de données utilisés pour les comparaisons d’événements uniques sont répertoriés dans le tableau 3.
L’ensemble de données «synthétique (TVE)» est créé sur la base d’un prédicteur additif $g(X,t)=F_0(X_0,t)\cdot 6–\mathrm{0,1}\cdot X_1+F_2(X_2,t)+F_3(X_3,t)$, où $F_0$, $F_2$ et $F_3$ sont des fonctions bivariées et non linéaires des entrées (voir le référentiel de code pour plus de détails) et $X_0,\dots ,X_3$ colonnes de caractéristiques comprises dans $X$

. De plus, 20 variables de bruit sont tirées de la distribution uniforme

$U(0,1)$.

Pour la comparaison avec ORSF, nous utilisons le score de Brier intégré $\text{SCI}\left(\tau \right)=\frac{1}{\tau }$, où $\widehat{\text{BS}}(t,\hat{S})$ est le score Brier estimé au moment $t$ pondéré par la probabilité inverse de censurer les poids [15] et $\hat{S}$ la fonction de probabilité de survie estimée de l’algorithme respectif. En plus, [23] rapporter l’indice C dépendant du temps [14]. Nous ne considérons ici que l’IBS car il mesure l’étalonnage ainsi que la discrimination, tandis que l’indice C ne mesure que cette dernière. Notez que l’IBS dépend du temps d’évaluation spécifique $\tau$

et différentes méthodes pourraient mieux fonctionner à différents moments d’évaluation. Par conséquent, nous calculons l’IBS pour trois moments différents, les quantiles de 25%, 50% et 75% des temps d’événement dans les données de test, ci-après désignés par Q25, Q50 et Q75, respectivement.

Pour la comparaison avec DeepHit, nous utilisons l’ensemble de données métabriques (cf. [28]) pour la comparaison d’un seul événement ainsi que deux ensembles de données CR. Les données «MGUS 2» sont décrites dans [26]. L’ensemble de données «synthétiques (TVE CR)» est simulé à l’aide d’un prédicteur additif identique à celui utilisé pour la simulation de données «synthétiques (TVE)» pour la première cause. Le prédicteur de la deuxième cause a une structure plus simple $F_0(X_0,t)+2\cdot X_4–.1\cdot X_5$, cependant, avec un risque de référence non proportionnel par rapport à $X_0\in \{–1,1\}$. Le nombre de variables de bruit est limité à 10 pour ce paramètre. Ici, nous rapportons l’indice C pondéré à côté du score de Brier pondéré, car il s’agissait de la principale mesure rapportée dans [28]. L’approche GBT (PEM) proposée pour la CR (cf. 2) est un modèle de risque spécifique à une cause, cependant, les paramètres des deux causes sont estimés conjointement et les risques des deux causes peuvent avoir des effets partagés (voir la figure 2). Le paramètre de simulation constitue donc une configuration difficile car il n’y a pas d’effets et d’optimisation partagés. la première cause favorisera les paramètres qui permettent des modèles flexibles tandis que l’optimisation w.r.t. la seconde cause favorise les modèles clairsemés et donc les paramètres qui restreindraient la flexibilité.

We now briefly describe algorithmic details and discuss the complexity of the resulting algorithms when using the proposed framework.

We have presented a general machine learning framework for time-to-event analysis based on a data augmentation strategy that reduces a large variety of survival analysis tasks to the optimization of a Poisson likelihood. We demonstrated its versatility and state-of-the-art performance. The availability of Poisson regression for most machine learning frameworks provides additional practical advantages. For example, photon-ML [37] is a scalable machine learning library for Apache Spark [36]

that has no native support for survival analysis, but implements generalized linear mixed models. Therefore, survival modeling with high cardinality random effects (frailty) is directly available using our framework. Similarly, lightGBM

[24], a high-performance implementation of GBT, currently has no implementation of survival methods, but could be also used for high-dimensional survival tasks based on PEMs, including reliability analysis or churn analysis with intermediate states.

This work has been funded by the German Federal Ministry of Education and Research (BMBF) under Grant No. 01IS18036A. The authors of this work take full responsibilities for its content.