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

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


1 introduction

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).

Figure 1: Une abstraction de l’analyse de survie pour différentes tâches. La structure
des données exponentielles par morceaux (PED) dépend des exigences de la tâche, par exemple,
troncature à gauche ou risques concurrents. Compte tenu du prétraitement approprié, le
l’étape d’estimation est indépendante du point de vue du calcul de la tâche de survie, sauf pour une utilisation appropriée des termes d’interaction.

Nos contributions

Nous définissons un cadre général d’apprentissage automatique pour l’analyse de survie basé sur des modèles exponentiels par morceaux (cf. Section 2). Dans ce cadre, différents concepts spécifiques à l’analyse des données temporelles peuvent être compris en termes d’augmentation des données et d’inclusion de termes d’interaction. En ré-exprimant la tâche de survie comme une tâche de régression de Poisson, une grande variété d’algorithmes devient disponible pour l’analyse de survie. Sur la base de l’approche proposée, nous implémentons un algorithme d’arbres boostés par gradient avec un effort de développement relativement faible et montrons qu’il atteint des performances de pointe (cf. section 3).

Travaux connexes

La communauté d’apprentissage automatique a développé de nombreuses méthodes très efficaces pour les paramètres de haute dimension dans différents domaines, y compris l’analyse de survie. Cependant, les méthodes et les implémentations individuelles ne prennent souvent en charge qu’un sous-ensemble des cas pertinents pour l’analyse temporelle des événements mentionnés ci-dessus. Par exemple, la forêt de survie aléatoire (RSF) proposée dans [22] a ensuite été étendu à la configuration des risques concurrents [21]

, mais ne prend pas en charge la troncature à gauche, les modèles TVF et TVE, ni les modèles multi-états. Une autre mise en œuvre populaire de forêts aléatoires

[35] prend uniquement en charge les données censurées à droite et les modèles de risques proportionnels. Une extension de RSF, le RSF oblique (ORSF, [23]) s’est avéré surpasser les autres algorithmes basés sur RSF, mais a les mêmes limites. En ce qui concerne TVF et TVE, un examen des méthodes d’analyse de survie basées sur les arbres et les forêts a indiqué que «la modélisation des caractéristiques variant dans le temps et des effets variant dans le temps mérite beaucoup plus d’attention» [6]. De même, un examen plus récent des méthodes d’apprentissage automatique pour l’analyse de survie [34] répertorie uniquement le modèle de Cox dépendant du temps [25, Ch. 9.2], et leurs extensions régularisées L1 et L2, comme possibilité d’inclusion de TVF.

Les méthodes basées sur l’apprentissage en profondeur pour les données temporelles des événements ont également suscité beaucoup d’attention ces derniers temps. Une utilisation précoce des réseaux de neurones pour les modèles de type Cox a été proposée par [10]. Plus récemment, [31] a présenté un cadre d’analyse de la survie d’un événement unique profond basé sur un processus latent conjoint pour les caractéristiques et les temps de survie à l’aide de familles exponentielles profondes. Pour les données de risques concurrentes, un cadre d’apprentissage en profondeur basé sur des processus gaussiens a été décrit dans [1]

. Un autre framework récent est DeepHit, qui peut gérer des risques concurrents à l’aide d’une fonction de perte personnalisée.

[28] et a été étendu pour gérer TVF [27], mais n’a pas
discuter de la troncature à gauche, des modèles multi-états et du TVE.

La stimulation a également été une technique populaire pour l’analyse de survie en haute dimension. Par example, [5] proposent une approche de boosting de type Cox pour l’estimation des aléas de sous-distribution proportionnelle. Un modèle multi-états flexible basé sur la vraisemblance partielle de Cox stratifiée dans le contexte du boosting basé sur un modèle [17] est présenté dans [32]

. De plus, une implémentation d’arbres à gradient boosté (GBT) pour le modèle Cox PH est également disponible pour l’implémentation populaire XGBoost.

[8], qui s’est également révélé performant par rapport à l’ORSF [23]. Récemment, [29] ont dérivé un algorithme personnalisé pour les arbres boostés par gradient qui prennent en charge TVF et démontrent que leur inclusion améliore les performances prédictives par rapport aux algorithmes boosters qui ne prennent pas en compte TVF.

Par rapport aux méthodes basées sur la régression de Cox, peu de publications ont développé des méthodes basées sur le modèle exponentiel par morceaux, sur lequel repose le cadre proposé ici. Parmi eux, une application précoce des réseaux de neurones à l’analyse de survie suggérée dans [30] et prolongé par [4]

. Ce dernier offre un cadre général basé sur la représentation de modèles linéaires généralisés via des réseaux de neurones à réaction, mais ne traite pas du MSM. Des arbres exponentiels par morceau avec TVF et des divisions basées sur la fonction de survie exponentielle par morceau ont été suggérés par

[19]. Une estimation basée sur les splines de la fonction de risque a été discutée dans [7], qui pourrait également être représenté via des réseaux de neurones (cf. [11]). Une estimation flexible des modèles multi-états basés sur un modèle exponentiel par morceaux avec des effets partagés à l’aide de la fusion structurée Lasso a été développée en [33]. Toutes ces méthodes peuvent être considérées comme des cas particuliers dans le cadre proposé. Par example, [4] pourrait être étendu à différentes architectures de réseaux neuronaux et MSM, [19] pourrait être étendu aux forêts.

2 Analyse de survie comme régression de Poisson

Dans le contexte de l’analyse de survie, une observation consiste généralement en un tuple , où est l’heure de l’événement observé pour l’unité d’observation , 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 est le

-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

, où et 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,

). 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

, où indique le type d’événement survenu à, si . 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 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(t|X)=P(T>t|X). Alors que certaines méthodes se concentrent sur l’estimation de directement, il est souvent plus pratique d’estimer le danger (log-)

(1)

à partir duquel suit comme

(2)

Ici, nous représentons (1) passant par

(3)

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 l’événement observé ou l’heure de censure et l’indicateur de censure ou d’événement respectif pour les unités d’observation . 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 . En répartissant le suivi, c’est-à-dire le laps de temps sous enquête, en intervalles avec points de coupure et cloisons
,
nous pouvons réécrire (3) en utilisant des taux de risque constants par pièce

(4)

avec une représentation du temps dans l’intervalle , par exemple., et la valeur du TVF en intervalle . 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 affecte le danger dans l’intervalle , 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

(5)

est le dernier intervalle dans lequel l’unité d’observation a été observée, de telle sorte que et

(6)

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 et avec la fonction de densité de Poisson, [13] a montré que la log-vraisemblance de Poisson

(7)

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 avec , le premier intervalle dans lequel l’unité d’observation fait partie de l’ensemble des risques. L’espérance est définie par . Pour l’estimation, est inclus comme compensation, donc le taux de risque est définie comme l’espérance conditionnelle d’avoir un événement dans l’intervalle 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 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

(8)

comme danger spécifique à la transition pour la transition indexée par , c’est à dire., est un index des transitions est l’état initial et un état transitoire ou absorbant. L’ensemble des transitions possibles est donné par , où désigne le nombre total d’états possibles. 6) nous définissons

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 est un index des transitions ; un exemple concret est donné dans le tableau 2. Pour chaque , 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 sont codés comme et tout le reste, c’est-à-dire la censure et la transition vers d’autres états, est codé comme . Ces ensembles de données spécifiques à la transition, chacun contenant un vecteur d’entités avec l’index de transition , sont ensuite concaténées. Notez que nous avons utilisé les mêmes points de partage d’intervalle pour toutes les transitions dans le tableau 1. Cependant, il serait également possible de choisir des seuils spécifiques à la transition , 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 sont définis par , où est le moment auquel entre en état .

Tableau 1: Structure des données après transformation au format de données exponentiel par morceaux dans le cadre des risques concurrents. Les lignes horizontales indiquent une nouvelle unité d’observation . Les doubles lignes horizontales indiquent une nouvelle transition indexée . Comme avant, est l’intervalle dans lequel a été observé en dernier, c.-à-d. . Caractéristiques dépend du temps via 1 1 0 1 1 1 2 0 0,3 1 2 1 0 0,5 1 3 1 0 1 1 3 2 0 0,5 1 3 3 1 1.2 1
1 1 0 1 2
1 2 1 0,3 2
2 1 0 0,5 2
3 1 0 1 2
3 2 0 0,5 2
3 3 0 1.2 2
Tableau 2: Transformation des données pour un exemple de risques concurrents hypotéthiques () pour 3 matières, . Matière
vécu un événement de type à , matière vécu un événement de type au moment , matière a été censuré à . Les tableaux présentent les données transformées avec des intervalles pour des causes (à gauche) et (droite). Ceux-ci peuvent être utilisés pour estimer la cause des dangers spécifiques, en appliquant l’algorithme à chacun des tableaux séparément ou provoquer des dangers spécifiques avec des effets potentiellement partagés, en empilant les tableaux et en utilisant comme une fonctionnalité.

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 . Le décalage () 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 basé sur de nouveaux points de données, le décalage doit être omis, sinon l’algorithme prédira (l’attente) au lieu de (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 https://github.com/adibender/pem.xgb.

3 Expériences

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 , où , et 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 colonnes de caractéristiques comprises dans

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

.

Pour la comparaison avec ORSF, nous utilisons le score de Brier intégré , où est le score Brier estimé au moment pondéré par la probabilité inverse de censurer les poids [15] et 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

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.

Tableau 3: Ensembles de données utilisés dans des expériences de référence pour comparaison avec ORSF.

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 , cependant, avec un risque de référence non proportionnel par rapport à . 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é.

Tableau 4: Ensembles de données utilisés pour la comparaison avec DeepHit. MGUS2 et Synthetic (TVE CR) sont des ensembles de données avec deux risques concurrents et une censure à droite supplémentaire.

3,1 Évaluation

Nous comparons quatre algorithmes, l’estimation non paramétrique de Kaplan Meier (référence) comme base de référence minimale, le modèle de risques proportionnels de Cox [9] (référence pour les effets linéaires, constants de temps), la forêt de survie aléatoire oblique (ORSF) [23] et DeepHit [28]. Pour chaque réplication expérimentale pour un ensemble de données spécifique, 70% des données sont attribuées au hasard en tant que données d’apprentissage et les 30% restants sont utilisés pour calculer les mesures d’évaluation à trois moments Q25, Q50 et Q75. Les algorithmes sont ajustés sur les données d’entraînement en utilisant une recherche aléatoire avec un budget fixe et une validation croisée quadruple. Chaque algorithme est ensuite recyclé sur l’ensemble de données d’apprentissage en utilisant le meilleur ensemble de paramètres avant de faire des prédictions finales sur l’ensemble de test.
La recherche aléatoire consiste en 20 itérations pour chaque algorithme. Pour le GBT (PEM), nous définissons l’espace de recherche comme suit (plage possible entre parenthèses): profondeur maximale de l’arbre , réduction minimale des pertes [0, 5], poids minimum de l’enfant , pourcentage de sous-échantillon (lignes) [0.5, 1], pourcentage de sous-échantillon d’entités dans chaque arborescence [.5, 1], Régularisation L2 [1, 3]. Le taux d’apprentissage est fixé à 0,05 et le nombre de tours à 5000, avec un arrêt prématuré après 50 tours sans amélioration. Pour l’ORSF, nous ajustons le paramètre de mélange de filet élastique (0, 1),
le paramètre qui pénalise la complexité du prédicteur linéaire dans chaque nœud (0,25, 0,75), nombre minimum d’événements à séparer le nœud et observations minimales pour diviser le nœud

. Pour DeepHit, nous utilisons 50 itérations de recherche aléatoires, où nous recherchons dans {1,2,3,5} couches partagées avec {50, 100, 200, 300} dimensions, {1,2,3,5} réseau spécifique à la cause couches aux dimensions {50,100,200,300}, ReLU, eLU ou Tanh comme fonction d’activation dans ces couches, une taille de lot en {32,64,128}, un maximum de 50000 itérations, un taux d’abandon de 0,6 (tiré du papier d’origine) et un apprentissage taux de 0,0001. Les paramètres spécifiques au réseau

et sont également choisis en fonction du papier d’origine et réglés sur et , respectivement, tandis que le paramètre spécifique au réseau varie dans la recherche aléatoire avec des valeurs possibles dans {0,1, 0,5, 1, 3, 5}.

3.2 Résultats

Les résultats des expériences basées sur la comparaison de scénarios d’événement unique avec ORSF
sont résumés dans le tableau 5. La méthode proposée effectue
bien dans de nombreux contextes par rapport à ORSF. Notamment, les deux algorithmes ne sont souvent pas beaucoup mieux
que les modèles Cox PH indiquant que l’hypothèse PH n’est pas violée fortement dans ces ensembles de données et que la taille de l’échantillon peut être trop petite pour détecter de petits écarts
w.r.t. à la non-linéarité des effets de caractéristiques, des effets d’interaction et des variations dans le temps
effets. Le paramètre «synthétique (TVE)» illustre qu’en présence de
TVE forte, non linéaire et non linéaire, notre approche surpasse clairement les autres méthodes.
Pour les données PBC, nous avons également effectué une analyse incluant TVF avec GBT (PEM). Dans ce cas, l’inclusion de TVF s’est traduite par une moins bonne performance (IBS de 4,3 (Q25), 6,4 (Q50) et 9,2 (Q75)), ce qui indique que l’inclusion de TVF conduit à un surajustement ou que la simple inclusion du dernier observé la valeur et le report de la dernière valeur ne sont pas appropriés dans ce cadre.

Tableau 5: Résultats d’expériences de référence pour des données d’événement unique comparant GBT (PEM) avec ORSF. Les chiffres en gras indiquent les meilleures performances pour chaque paramètre.

Tableau 6 résume les résultats des comparaisons avec DeepHit. Le GBT (PEM) affiche à nouveau de bonnes performances globales. Pour l’ensemble de données synthétiques, notre méthode surpasse clairement les autres approches car elle est capable d’estimer la non-linéarité ainsi que la variation dans le temps. On the MGUS 2 data set, DeepHit shows the best performance for cause 1, while GBT (PEM) outperforms the other approaches for cause 2. On the synthetic data, the cause-specific Cox-PH model shows good discrimination (C-Index) for the second cause, but is worse than GBT (PEM) and DeepHit w.r.t. to the Brier Score.

Table 6: Results of benchmark experiments comparing GBT (PEM) with DeepHit for single event and CR data.
Bold numbers indicate the best performance for each setting.

4 Algorithmic Details and Complexity Analysis

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

Algorithmic Details

The proposed framework is general in the sense that it transforms a survival task into a regression task. Nevertheless, different methods (and algorithms) have different strengths and weaknesses and different strategies can be applied to specify various alternative models within this framework. For example, in tree based methods, time-variation of feature effects could be controlled by allowing interactions of the time variable only with a subset of features, e.g. based on prior information, and similarly in order to control shared vs. transition specific effects in the multi-state setting. Tree-based methods are particularly intuitive when it comes to understanding the integration of TVE and extension to multi-state models via interaction terms into the model. This is illustrated in Figure 2. For example, in panel (A) of Figure 2, features and split points before the split w.r.t. time indicate feature effects common to all time-points. Once the data in panel (A) is split w.r.t. temps , the predicted hazard will be different for intervals with et for observations with . Similarly, in a multi-state setting (panel (B) in Figure 2), splits above the split w.r.t. indicate shared effects for all transitions, while splits below indicate different effects for transitions vs. . Forcing a split w.r.t. à at the root node would be equivalent to an estimation of cause specific hazards on each subset and no shared effects.

Figure 2: Illustration of how TVE and shared vs. transitions specific effects can be understood in terms of feature interactions in tree based models.

Neural networks are particularly flexible when it comes to the specification of different PEMs. For example, the network could be split in two subnetworks, one for the temporal component, one for features, which is equivalent to the specification of a proportional hazards model, while allowing for non-linearity and high-dimensional interactions in feature effects. Similarly, defining subnetworks of the time variable for each category of a categorical feature would imply a stratified proportional hazards model.

Complexity Analysis

As described in the literature review, various approaches exist that account for special survival characteristics like TVF, CR or continuous time-scale prediction by altering the underlying method.
While adapting the structure of the algorithm itself potentially increases the complexity of the method,
our approach leaves the algorithm of choice unchanged as different time points and transitions are simply included as features. This allows to employ commonly used prediction methods without introducing further algorithmic complexity.
We note, however, that our approach might be improved upon in terms of scaling with respect to the number of intervals relative to the number of observations . In the worst case, the number of total data points is quadratic in (or more precisely ) when one interval cut-point is introduced for each observed event or censoring time.
We therefore propose a refinement of the presented method that improves run-times without forfeiting performance. Instead of setting cut-points at all unique event times, we suggest to define cut points more sparsely, for example, based on a sub-sample of the original data.

To investigate this strategy we conduct a scaling experiment where the sample size was consecutively doubled starting from up to . For each sample size, ten replications of one experiment as described for the “synthetic TVE” setting in Section 3 were performed and the elapsed time (hours) as well as performance (IBS) for two different strategies of cut-point selection was measured. The first strategy (full) uses all event times ( where ) as cut-points. The second strategy (sub-sample) is equivalent to the first strategy, but event times were chosen based on a sub-sample of , selected randomly from the training data in each iteration.
Results in Table 7 show that the “sub-sample” strategy leads to an approximately linear increase in computation time while the performance remains virtually unchanged. Potentially, a sparser choice of cut-points could also lead to a more robust and thus improved hazard estimation, as more events are available in each interval, but we did not conduct a formal investigation in that regard.

Table 7: Results of the scaling experiments with the number of observation in the simulated data. “strategy” refers to the way interval cut-points were selected with “full” splits at all unique event times ( where ) and “sub-sample” refers to the selection of cut-points based on unique event times based on a random sub-sample of size (but all observations were used for estimation). Mean time (in hours) and IBS over 10 replications are reported for each setting.

5 Conclusion

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.

Acknowledgements:

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.



Source de l’article

A découvrir